formal_lp/old/formal_interval/verifier.hl

   1 needs "../formal_lp/formal_interval/second_approx.hl";;
   2 needs "../formal_lp/formal_interval/interval_1d/recurse.hl";;
   3
   4
   5 (***********************************)
   6
   7 let cell_pass = new_definition `cell_pass f int <=> !x. interval_arith x int ==> f x < &0`;;
   8
   9 let dest_cell_pass tm =
  10   let lhs, int_tm = dest_comb tm in
  11     rand lhs, int_tm;;
  12
  13 let CELL_PASS_GLUE = prove(`!f x t z. cell_pass f (x, t) /\ cell_pass f (t, z) ==> cell_pass f (x, z)`,
  14                            REWRITE_TAC[cell_domain; cell_pass; interval_arith] THEN
  15                              REPEAT STRIP_TAC THEN
  16                              DISJ_CASES_TAC (REAL_ARITH `x' <= t \/ t <= x'`) THENL
  17                              [
  18                                FIRST_X_ASSUM ((fun th -> ALL_TAC) o SPEC `x':real`) THEN
  19                                  FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
  20                                FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[]
  21                              ]);;
  22
  23 let CELL_PASS_GLUE' = RULE CELL_PASS_GLUE;;
  24
  25 let cell_pass_glue pass1 pass2 =
  26   let f_tm, int1 = dest_cell_pass (concl pass1) and
  27       int2 = rand (concl pass2) in
  28   let x_tm, t_tm = dest_pair int1 and
  29       z_tm = rand int2 in
  30   let th0 = INST[f_tm, f_var_fun; x_tm, x_var_real; t_tm, t_var_real; z_tm, z_var_real] CELL_PASS_GLUE' in
  31     (MY_PROVE_HYP pass1 o MY_PROVE_HYP pass2) th0;;
  32
  33
  34
  35
  36 let CELL_PASS_LEMMA' = (UNDISCH_ALL o PURE_REWRITE_RULE[GSYM IMP_IMP] o prove)
  37   (`bounded_on_int f int (f_lo, f_hi) /\ (f_hi < &0 <=> T) ==> cell_pass f int`,
  38    REWRITE_TAC[bounded_on_int; cell_pass; interval_arith] THEN
  39      REPEAT STRIP_TAC THEN
  40      MATCH_MP_TAC REAL_LET_TRANS THEN
  41      EXISTS_TAC `f_hi:real` THEN ASM_REWRITE_TAC[] THEN
  42      FIRST_X_ASSUM (MP_TAC o SPEC `x:real`) THEN ASM_SIMP_TAC[]);;
  43
  44
  45 let taylor_cell_pass pp taylor_int_th =
  46   let bound_th = eval_taylor_f_bounds pp taylor_int_th in
  47   let lhs, f_bounds = dest_comb (concl bound_th) in
  48   let lhs2, int_tm = dest_comb lhs in
  49   let f_tm = rand lhs2 in
  50   let f_lo, f_hi = dest_pair f_bounds in
  51   let f_hi_lt0 = float_lt0 f_hi in
  52     if (fst o dest_const o rand o concl) f_hi_lt0 <> "T" then
  53       failwith "taylor_cell_pass: f_hi >= 0"
  54     else
  55       let th0 = INST[f_tm, f_var_fun; int_tm, int_var; f_lo, f_lo_var; f_hi, f_hi_var] CELL_PASS_LEMMA' in
  56         (MY_PROVE_HYP bound_th o MY_PROVE_HYP f_hi_lt0) th0;;
  57
  58
  59 (*************************************)
  60
  61
  62 open Taylor;;
  63 open Univariate;;
  64 open Interval;;
  65 open Recurse;;
  66
  67 (******************************************)
  68
  69
  70 let DECIMAL' = SPEC_ALL DECIMAL;;
  71
  72 (* Unary interval operations *)
  73 let unary_interval_operations =
  74   let table = Hashtbl.create 10 in
  75   let add = Hashtbl.add in
  76     add table `--` (fun pp -> float_interval_neg);
  77     add table `inv` float_interval_inv;
  78     add table `sqrt` float_interval_sqrt;
  79     add table `atn` float_interval_atn;
  80     add table `acs` float_interval_acs;
  81     table;;
  82
  83
  84 (* Binary interval operations *)
  85 let binary_interval_operations =
  86   let table = Hashtbl.create 10 in
  87   let add = Hashtbl.add in
  88     add table `+` float_interval_add;
  89     add table `-` float_interval_sub;
  90     add table `*` float_interval_mul;
  91     add table `/` float_interval_div;
  92     table;;
  93
  94 (* Interval approximations of constants *)
  95 let interval_constants =
  96   let table = Hashtbl.create 10 in
  97   let add = Hashtbl.add in
  98     add table `pi` (fun pp -> pi_approx_array.(pp));
  99     table;;
 100
 101
 102 (* Builds a function for evaluating an interval approximation of the given expression *)
 103 let rec build_eval_interval expr_tm =
 104   (* Constants *)
 105   if is_const expr_tm then
 106     Hashtbl.find interval_constants expr_tm
 107   else
 108     let ltm, r_tm = dest_comb expr_tm in
 109       (* Unary operations *)
 110       if is_const ltm then
 111         (* & *)
 112         if ltm = amp_op_real then
 113           let n = dest_numeral r_tm in
 114             (fun pp -> mk_float_interval_num n)
 115         else
 116           let r_fun = build_eval_interval r_tm in
 117           let f = Hashtbl.find unary_interval_operations ltm in
 118             (fun pp -> f pp (r_fun pp))
 119       else
 120         (* Binary operations *)
 121         let op, l_tm = dest_comb ltm in
 122           (* DECIMAL *)
 123           if (fst o dest_const) op = "DECIMAL" then
 124             let n, d = dest_numeral l_tm, dest_numeral r_tm in
 125               (fun pp ->
 126                  let (/), (!) = float_interval_div pp, mk_float_interval_num in
 127                  let int_th = !n / !d in
 128                  let eq_th = INST[l_tm, x_var_num; r_tm, y_var_num] DECIMAL' in
 129                    norm_interval int_th eq_th)
 130           else
 131             let lhs = build_eval_interval l_tm and
 132                 rhs = build_eval_interval r_tm in
 133             let f = Hashtbl.find binary_interval_operations op in
 134               (fun pp -> f pp (lhs pp) (rhs pp));;
 135
 136
 137 (* Converts a float term to the corresponding rational number *)
 138 let num_of_float_tm tm =
 139   let s, n_tm, e_tm = dest_float tm in
 140   let b = Num.num_of_int base in
 141   let m = Num.num_of_int min_exp in
 142   let ( * ), (^), (-), (!) = ( */ ), ( **/ ), (-/), raw_dest_hash in
 143   let r = !n_tm * (b ^ (!e_tm - m)) in
 144     if s = "T" then minus_num r else r;;
 145
 146
 147 let float_of_float_tm tm =
 148   float_of_num (num_of_float_tm tm);;
 149
 150
 151 (**************************************)
 152
 153 let taylor_interval_div_alt pp th1 th2 =
 154   let domain_th, _, _ = dest_taylor_interval_th th1 in
 155   let rhs = taylor_interval_compose pp eval_taylor_inv (fun _ _ -> th2) in
 156     taylor_interval_mul pp th1 (rhs domain_th);;
 157
 158
 159 let mk_formal_binop taylor_binop arg1 arg2 =
 160   fun pp domain_th ->
 161     let lhs = arg1 pp domain_th and
 162         rhs = arg2 pp domain_th in
 163       taylor_binop pp lhs rhs;;
 164
 165
 166 let rec mk_taylor_functions pp expr_tm =
 167   (* Constant *)
 168   try
 169     let eval_f = build_eval_interval expr_tm in
 170     let int_th = eval_f pp in
 171     let lo_tm, hi_tm = (dest_pair o rand o concl) int_th in
 172     let float_int = mk_interval (float_of_float_tm lo_tm, float_of_float_tm hi_tm) in
 173       Scale (unit, float_int), (fun pp -> eval_taylor_const_interval int_th)
 174   with _ ->
 175     (* Variable *)
 176     if is_var expr_tm then
 177       if (expr_tm = x_var_real) then x1, (fun pp -> eval_taylor_x)
 178       else failwith "mk_taylor_functions: only x variable is allowed"
 179     else
 180       (* Constant (operation) *)
 181       if is_const expr_tm then
 182         let name,_ = dest_const expr_tm in
 183           if name = "real_inv" then
 184             Uni (mk_primitiveU uinv), eval_taylor_inv
 185           else if name = "sqrt" then
 186             Uni (mk_primitiveU usqrt), eval_taylor_sqrt
 187           else if name = "atn" then
 188             Uni (mk_primitiveU uatan), eval_taylor_atn
 189           else if name = "acs" then
 190             Uni (mk_primitiveU uacos), eval_taylor_acs
 191           else
 192             failwith ("mk_taylor_functions: unknown constant: "^name)
 193       else
 194         (* Combination *)
 195         let lhs, r_tm = dest_comb expr_tm in
 196         let r, r2 = mk_taylor_functions pp r_tm in
 197           (* Unary operations *)
 198           if lhs = neg_op_real then
 199             Scale (r, ineg one), (fun pp domain_th -> taylor_interval_neg (r2 pp domain_th))
 200           else
 201             (* Composite *)
 202             try
 203               let l, l2 = mk_taylor_functions pp lhs in
 204                 if r_tm = x_var_real then
 205                   l, l2
 206                 else
 207                   Composite (l,r), (fun pp domain_th -> taylor_interval_compose pp l2 r2 domain_th)
 208             with Failure _ ->
 209               (* Binary operations *)
 210               if is_comb lhs then
 211                 let op, l_tm = dest_comb lhs in
 212                 let l, l2 = mk_taylor_functions pp l_tm in
 213                   if op = add_op_real then
 214                     Plus (l,r), mk_formal_binop taylor_interval_add l2 r2
 215                   else if op = sub_op_real then
 216                     Minus (l,r), mk_formal_binop taylor_interval_sub l2 r2
 217                   else if op = mul_op_real then
 218                     Product (l,r), mk_formal_binop taylor_interval_mul l2 r2
 219                   else if op = div_op_real then
 220                     Product (l, Uni_compose (uinv,r)), mk_formal_binop taylor_interval_div_alt l2 r2
 221                   else
 222                     failwith ("mk_taylor_functions: unknown binary operation "^string_of_term op)
 223               else
 224                 failwith ("mk_taylor_functions: error: "^string_of_term lhs);;
 225
 226 (**********************************)
 227
 228
 229 let rec result_size r =
 230   match r with
 231     | Result_false -> failwith "False result detected"
 232     | Result_glue (r1, r2) -> result_size r1 + result_size r2
 233     | Result_pass _ -> 1;;
 234
 235
 236 let prove_result pp formal_f result lo_tm hi_tm =
 237   let pass_counter = ref 0 in
 238   let r_size = result_size result in
 239   let rec rec_prove result lo_tm hi_tm =
 240     let domain_th = mk_center_domain pp lo_tm hi_tm in
 241       match result with
 242         | Result_false -> failwith "False result"
 243         | Result_pass _ ->
 244             pass_counter := !pass_counter + 1;
 245             let th0 = taylor_cell_pass pp (formal_f pp domain_th) in
 246             let _ = report (sprintf "Result_pass: %d/%d" !pass_counter r_size) in
 247               th0
 248       | Result_glue (r1, r2) ->
 249           let _, y_tm, _, _ = dest_cell_domain (concl domain_th) in
 250           let th1 = rec_prove r1 lo_tm y_tm and
 251               th2 = rec_prove r2 y_tm hi_tm in
 252             cell_pass_glue th1 th2 in
 253     rec_prove result lo_tm hi_tm;;
 254
 255
 256 let CELL_PASS_IMP = prove
 257   (`cell_pass f (lo, hi) ==>
 258      interval_arith x (lo, w1) ==>
 259      interval_arith z (w2, hi) ==>
 260      !t. t IN real_interval [x, z] ==> f t < &0`,
 261   REWRITE_TAC[cell_pass; interval_arith; IN_REAL_INTERVAL] THEN REPEAT STRIP_TAC THEN
 262     FIRST_X_ASSUM MATCH_MP_TAC THEN
 263     ASM_REAL_ARITH_TAC);;
 264
 265
 266 let prove_ineq1d pp expr_tm lo_tm hi_tm =
 267   let f, formal_f = mk_taylor_functions pp expr_tm in
 268   let lo_int, hi_int = build_eval_interval lo_tm pp, build_eval_interval hi_tm pp in
 269   let lo1_tm = (fst o dest_pair o rand o concl) lo_int and
 270       hi1_tm = (snd o dest_pair o rand o concl) hi_int in
 271   let lo, hi = float_of_float_tm lo1_tm, float_of_float_tm hi1_tm in
 272   let _ = report (sprintf "Starting verification from %f to %f" lo hi) in
 273
 274   let result =
 275     let opt = {
 276       allow_sharp = false;
 277       iteration_count = 0;
 278       iteration_limit = 0;
 279       recursion_depth = 50
 280     } in
 281       recursive_verifier (0, lo, hi, f, opt) in
 282
 283   let th0 = prove_result pp formal_f result lo1_tm hi1_tm in
 284     REWRITE_RULE[] (MATCH_MP (MATCH_MP (MATCH_MP CELL_PASS_IMP th0) lo_int) hi_int);;