formal_lp/old/formal_interval/interval_1d/test1d.hl

   1 needs "../formal_lp/formal_interval/interval_1d/recurse.hl";;
   2 needs "../formal_lp/formal_interval/second_approx.hl";;
   3 needs "../formal_lp/formal_interval/verifier.hl";;
   4
   5
   6 open Taylor;;
   7 open Univariate;;
   8 open Interval;;
   9 open Recurse;;
  10
  11 let x_var_real = `x:real` and
  12     add_op_real = `(+):real->real->real` and
  13     sub_op_real = `(-):real->real->real` and
  14     mul_op_real = `( * ):real->real->real` and
  15     div_op_real = `(/):real->real->real` and
  16     neg_op_real = `(--):real->real`;;
  17
  18 let taylor_interval_div_alt pp th1 th2 =
  19   let domain_th, _, _ = dest_taylor_interval_th th1 in
  20   let rhs = taylor_interval_compose pp eval_taylor_inv (fun _ _ -> th2) in
  21     taylor_interval_mul pp th1 (rhs domain_th);;
  22
  23
  24 let mk_formal_binop taylor_binop arg1 arg2 =
  25   fun pp domain_th ->
  26     let lhs = arg1 pp domain_th and
  27         rhs = arg2 pp domain_th in
  28       taylor_binop pp lhs rhs;;
  29
  30
  31 let rec eval_interval pp expr_tm =
  32   let ltm, r_tm = dest_comb expr_tm in
  33     if is_comb ltm then
  34       let op, l_tm = dest_comb ltm in
  35       let l_th = eval_interval pp l_tm in
  36       let r_th = eval_interval pp r_tm in
  37         if op = plus_op_real then
  38           float_interval_add pp l_th r_th
  39         else if op = mul_op_real then
  40           float_interval_mul pp l_th r_th
  41         else if op = div_op_real then
  42           float_interval_div pp l_th r_th
  43         else if op = minus_op_real then
  44           float_interval_sub pp l_th r_th
  45         else
  46           failwith ("Unknown operation: " ^ (fst o dest_const) op)
  47     else
  48       if ltm = neg_op_real then
  49         let r_th = eval_interval pp r_tm in
  50           float_interval_neg r_th
  51       else
  52         mk_float_interval_num (dest_numeral r_tm);;
  53
  54
  55 let rat_term_of_term tm =
  56   let th0 = REWRITE_CONV[DECIMAL] tm in
  57   let th1 = (REAL_RAT_REDUCE_CONV o rand o concl) th0 in
  58     GEN_REWRITE_RULE RAND_CONV [th1] th0;;
  59
  60
  61 let float_of_term tm =
  62     (float_of_num o rat_of_term o rand o concl o rat_term_of_term) tm;;
  63
  64
  65 let interval_of_term pp tm =
  66   let eq_th = rat_term_of_term tm in
  67   let int_th = eval_interval pp (rand (concl eq_th)) in
  68     REWRITE_RULE[SYM eq_th] int_th;;
  69
  70
  71 let rec mk_taylor_functions pp expr_tm =
  72   (* Variable *)
  73   if is_var expr_tm then
  74     if (expr_tm = x_var_real) then x1, (fun pp -> eval_taylor_x)
  75     else failwith "mk_taylor_functions: only x variable is allowed"
  76   else
  77     (* Constant *)
  78     if is_const expr_tm then
  79       let name,_ = dest_const expr_tm in
  80         if name = "inv" then
  81           Uni (mk_primitiveU uinv), eval_taylor_inv
  82         else if name = "sqrt" then
  83           Uni (mk_primitiveU usqrt), eval_taylor_sqrt
  84         else if name = "atn" then
  85           Uni (mk_primitiveU uatan), eval_taylor_atn
  86         else
  87           failwith ("mk_taylor_functions: unknown constant: "^name)
  88     else
  89       (* Number *)
  90       try
  91         let rat_th = rat_term_of_term expr_tm in
  92         let rat_tm = rand (concl rat_th) in
  93         let n = (float_of_num o rat_of_term) rat_tm in
  94         let int_th =
  95           let th0 = eval_interval pp rat_tm in
  96             REWRITE_RULE[SYM rat_th] th0 in
  97           Scale (unit, mk_interval (n, n)), (fun pp -> eval_taylor_const_interval int_th)
  98       with Failure _ ->
  99         (* Combination *)
 100         let lhs, r_tm = dest_comb expr_tm in
 101         let r, r2 = mk_taylor_functions pp r_tm in
 102           (* Unary operations *)
 103           if lhs = neg_op_real then
 104             Scale (r, ineg one), (fun pp -> eval_taylor_x)
 105           else
 106             (* Composite *)
 107             try
 108               let l, l2 = mk_taylor_functions pp lhs in
 109                 if r_tm = x_var_real then
 110                   l, l2
 111                 else
 112                   Composite (l,r), (fun pp domain_th -> taylor_interval_compose pp l2 r2 domain_th)
 113             with Failure _ ->
 114               (* Binary operations *)
 115               if is_comb lhs then
 116                 let op, l_tm = dest_comb lhs in
 117                 let l, l2 = mk_taylor_functions pp l_tm in
 118                   if op = add_op_real then
 119                     Plus (l,r), mk_formal_binop taylor_interval_add l2 r2
 120                   else if op = sub_op_real then
 121                     Minus (l,r), mk_formal_binop taylor_interval_sub l2 r2
 122                   else if op = mul_op_real then
 123                     Product (l,r), mk_formal_binop taylor_interval_mul l2 r2
 124                   else if op = div_op_real then
 125                     Product (l, Uni_compose (uinv,r)), mk_formal_binop taylor_interval_div_alt l2 r2
 126                   else
 127                     failwith ("mk_taylor_functions: unknown binary operation "^string_of_term op)
 128               else
 129                 failwith ("mk_taylor_functions: error: "^string_of_term lhs);;
 130
 131 (**********************************)
 132
 133 let rec result_size r =
 134   match r with
 135     | Result_false -> failwith "False result detected"
 136     | Result_glue (r1, r2) -> result_size r1 + result_size r2
 137     | Result_pass _ -> 1;;
 138
 139
 140
 141 let atan1 = `let P00, P01, P02 = #4.26665376811246382, #3.291955407318148, #0.281939359003812 in
 142   let Q00, Q01 = #4.26665376811354027, #4.714173328637605 in
 143   let x1 = x * x in
 144   let x2 = x1 * x1 in
 145     x * (x1 * P01 + x2 * P02 + P00) / (x1 * Q01 + x2 + Q00) - atn x - #0.001`;;
 146 let atan1_tm = (rand o concl o REPEATC let_CONV) atan1;;
 147
 148 (* tan(pi/16) = 0.198912367379658 *)
 149
 150 let atn1_f, atn1_formal = mk_taylor_functions 15 atan1_tm;;
 151 let result =
 152   let opt = {
 153     allow_sharp = false;
 154     iteration_count = 0;
 155     iteration_limit = 0;
 156     recursion_depth = 20
 157   } in
 158     recursive_verifier (0, 0.0, 0.2, atn1_f, opt);;
 159
 160 result_size result;;
 161
 162
 163 let prove_result pp formal_f result lo_tm hi_tm =
 164   let pass_counter = ref 0 in
 165   let r_size = result_size result in
 166   let rec rec_prove result lo_tm hi_tm =
 167     let domain_th = mk_center_domain pp lo_tm hi_tm in
 168       match result with
 169         | Result_false -> failwith "False result"
 170         | Result_pass _ ->
 171             pass_counter := !pass_counter + 1;
 172             let th0 = taylor_cell_pass pp (formal_f pp domain_th) in
 173             let _ = report (sprintf "Result_pass: %d/%d" !pass_counter r_size) in
 174               th0
 175       | Result_glue (r1, r2) ->
 176           let _, y_tm, _, _ = dest_cell_domain (concl domain_th) in
 177           let th1 = rec_prove r1 lo_tm y_tm and
 178               th2 = rec_prove r2 y_tm hi_tm in
 179             cell_pass_glue th1 th2 in
 180     rec_prove result lo_tm hi_tm;;
 181
 182
 183
 184 let lo_tm, hi_tm = mk_float 0 min_exp, mk_float 2 (min_exp - 1);;
 185 prove_result 15 atn1_formal result lo_tm hi_tm;;
 186
 187
 188 let CELL_PASS_IMP = prove
 189   (`cell_pass f (lo, hi) ==>
 190      interval_arith x (lo, w1) ==>
 191      interval_arith z (w2, hi) ==>
 192      !t. t IN real_interval [x, z] ==> f t < &0`,
 193   REWRITE_TAC[cell_pass; interval_arith; IN_REAL_INTERVAL] THEN REPEAT STRIP_TAC THEN
 194     FIRST_X_ASSUM MATCH_MP_TAC THEN
 195     ASM_REAL_ARITH_TAC);;
 196
 197
 198 let prove_ineq1d pp expr_tm lo_tm hi_tm =
 199   let f, formal_f = mk_taylor_functions pp expr_tm in
 200   let lo, hi = float_of_term lo_tm, float_of_term hi_tm in
 201   let lo_int, hi_int = interval_of_term pp lo_tm, interval_of_term pp hi_tm in
 202   let lo1_tm = (fst o dest_pair o rand o concl) lo_int and
 203       hi1_tm = (snd o dest_pair o rand o concl) hi_int in
 204   let _ = report (sprintf "Starting verification from %f to %f" lo hi) in
 205
 206   let result =
 207     let opt = {
 208       allow_sharp = false;
 209       iteration_count = 0;
 210       iteration_limit = 0;
 211       recursion_depth = 50
 212     } in
 213       recursive_verifier (0, lo, hi, f, opt) in
 214
 215   let th0 = prove_result pp formal_f result lo1_tm hi1_tm in
 216     REWRITE_RULE[] (MATCH_MP (MATCH_MP (MATCH_MP CELL_PASS_IMP th0) lo_int) hi_int);;
 217
 218
 219 prove_ineq1d 10 `atn x - x` `#0.001` `#30`;;
 220
 221
 222
 223
 224 result_size result;;
 225 let domain_th = mk_center_domain 10 (mk_float 175 (min_exp - 3)) (mk_float 2 (min_exp - 1));;
 226 let th1 = atn1_formal 10 domain_th;;
 227 test 1 (atn1_formal 10) domain_th;;
 228 evalf atn1_f 0.0 0.025;;
 229 eval_taylor_f_bounds 10 th1;;
 230
 231
 232
 233 (**********************************)
 234
 235 let pp = 10;;
 236 let domain_th = mk_center_domain pp (mk_float 1 (min_exp - 1)) (mk_float 8 (min_exp - 1));;
 237
 238 let f, f2 = mk_taylor_functions 10 `atn x / x - &1 / &3`;;
 239
 240
 241
 242
 243
 244 evalf f 0.1 0.8;;
 245 f2 pp domain_th;;
 246
 247
 248 let opt = {
 249   allow_sharp = false;
 250   iteration_count = 0;
 251   iteration_limit = 0;
 252   recursion_depth = 100
 253 } in
 254   recursive_verifier (0, 0.1, 0.8, f, opt);;
 255
 256
 257 let int1 = f2 pp (mk_center_domain pp (mk_float 1 49) (mk_float 14375 45));;
 258 eval_taylor_f_bounds pp int1;;