needs "../formal_lp/formal_interval/second_approx.hl";;
needs "../formal_lp/formal_interval/interval_1d/recurse.hl";;


(***********************************)

let cell_pass = new_definition `cell_pass f int <=> !x. interval_arith x int ==> f x < &0`;;
let dest_cell_pass tm = let lhs, int_tm = dest_comb tm in rand lhs, int_tm;;
let CELL_PASS_GLUE = 
prove(`!f x t z. cell_pass f (x, t) /\ cell_pass f (t, z) ==> cell_pass f (x, z)`,
REWRITE_TAC[cell_domain; cell_pass; interval_arith] THEN REPEAT STRIP_TAC THEN DISJ_CASES_TAC (REAL_ARITH `x' <= t \/ t <= x'`) THENL [ FIRST_X_ASSUM ((fun th -> ALL_TAC) o SPEC `x':real`) THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[]; FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[] ]);;
let CELL_PASS_GLUE' = RULE CELL_PASS_GLUE;; let cell_pass_glue pass1 pass2 = let f_tm, int1 = dest_cell_pass (concl pass1) and int2 = rand (concl pass2) in let x_tm, t_tm = dest_pair int1 and z_tm = rand int2 in 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 (MY_PROVE_HYP pass1 o MY_PROVE_HYP pass2) th0;; let CELL_PASS_LEMMA' = (UNDISCH_ALL o PURE_REWRITE_RULE[GSYM IMP_IMP] o prove) (`bounded_on_int f int (f_lo, f_hi) /\ (f_hi < &0 <=> T) ==> cell_pass f int`, REWRITE_TAC[bounded_on_int; cell_pass; interval_arith] THEN REPEAT STRIP_TAC THEN MATCH_MP_TAC REAL_LET_TRANS THEN EXISTS_TAC `f_hi:real` THEN ASM_REWRITE_TAC[] THEN FIRST_X_ASSUM (MP_TAC o SPEC `x:real`) THEN ASM_SIMP_TAC[]);; let taylor_cell_pass pp taylor_int_th = let bound_th = eval_taylor_f_bounds pp taylor_int_th in let lhs, f_bounds = dest_comb (concl bound_th) in let lhs2, int_tm = dest_comb lhs in let f_tm = rand lhs2 in let f_lo, f_hi = dest_pair f_bounds in let f_hi_lt0 = float_lt0 f_hi in if (fst o dest_const o rand o concl) f_hi_lt0 <> "T" then failwith "taylor_cell_pass: f_hi >= 0" else 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 (MY_PROVE_HYP bound_th o MY_PROVE_HYP f_hi_lt0) th0;; (*************************************) open Taylor;; open Univariate;; open Interval;; open Recurse;; (******************************************) let DECIMAL' = SPEC_ALL DECIMAL;; (* Unary interval operations *) let unary_interval_operations = let table = Hashtbl.create 10 in let add = Hashtbl.add in add table `--` (fun pp -> float_interval_neg); add table `inv` float_interval_inv; add table `sqrt` float_interval_sqrt; add table `atn` float_interval_atn; add table `acs` float_interval_acs; table;; (* Binary interval operations *) let binary_interval_operations = let table = Hashtbl.create 10 in let add = Hashtbl.add in add table `+` float_interval_add; add table `-` float_interval_sub; add table `*` float_interval_mul; add table `/` float_interval_div; table;; (* Interval approximations of constants *) let interval_constants = let table = Hashtbl.create 10 in let add = Hashtbl.add in add table `pi` (fun pp -> pi_approx_array.(pp)); table;; (* Builds a function for evaluating an interval approximation of the given expression *) let rec build_eval_interval expr_tm = (* Constants *) if is_const expr_tm then Hashtbl.find interval_constants expr_tm else let ltm, r_tm = dest_comb expr_tm in (* Unary operations *) if is_const ltm then (* & *) if ltm = amp_op_real then let n = dest_numeral r_tm in (fun pp -> mk_float_interval_num n) else let r_fun = build_eval_interval r_tm in let f = Hashtbl.find unary_interval_operations ltm in (fun pp -> f pp (r_fun pp)) else (* Binary operations *) let op, l_tm = dest_comb ltm in (* DECIMAL *) if (fst o dest_const) op = "DECIMAL" then let n, d = dest_numeral l_tm, dest_numeral r_tm in (fun pp -> let (/), (!) = float_interval_div pp, mk_float_interval_num in let int_th = !n / !d in let eq_th = INST[l_tm, x_var_num; r_tm, y_var_num] DECIMAL' in norm_interval int_th eq_th) else let lhs = build_eval_interval l_tm and rhs = build_eval_interval r_tm in let f = Hashtbl.find binary_interval_operations op in (fun pp -> f pp (lhs pp) (rhs pp));; (* Converts a float term to the corresponding rational number *) let num_of_float_tm tm = let s, n_tm, e_tm = dest_float tm in let b = Num.num_of_int base in let m = Num.num_of_int min_exp in let ( * ), (^), (-), (!) = ( */ ), ( **/ ), (-/), raw_dest_hash in let r = !n_tm * (b ^ (!e_tm - m)) in if s = "T" then minus_num r else r;; let float_of_float_tm tm = float_of_num (num_of_float_tm tm);; (**************************************) let taylor_interval_div_alt pp th1 th2 = let domain_th, _, _ = dest_taylor_interval_th th1 in let rhs = taylor_interval_compose pp eval_taylor_inv (fun _ _ -> th2) in taylor_interval_mul pp th1 (rhs domain_th);; let mk_formal_binop taylor_binop arg1 arg2 = fun pp domain_th -> let lhs = arg1 pp domain_th and rhs = arg2 pp domain_th in taylor_binop pp lhs rhs;; let rec mk_taylor_functions pp expr_tm = (* Constant *) try let eval_f = build_eval_interval expr_tm in let int_th = eval_f pp in let lo_tm, hi_tm = (dest_pair o rand o concl) int_th in let float_int = mk_interval (float_of_float_tm lo_tm, float_of_float_tm hi_tm) in Scale (unit, float_int), (fun pp -> eval_taylor_const_interval int_th) with _ -> (* Variable *) if is_var expr_tm then if (expr_tm = x_var_real) then x1, (fun pp -> eval_taylor_x) else failwith "mk_taylor_functions: only x variable is allowed" else (* Constant (operation) *) if is_const expr_tm then let name,_ = dest_const expr_tm in if name = "real_inv" then Uni (mk_primitiveU uinv), eval_taylor_inv else if name = "sqrt" then Uni (mk_primitiveU usqrt), eval_taylor_sqrt else if name = "atn" then Uni (mk_primitiveU uatan), eval_taylor_atn else if name = "acs" then Uni (mk_primitiveU uacos), eval_taylor_acs else failwith ("mk_taylor_functions: unknown constant: "^name) else (* Combination *) let lhs, r_tm = dest_comb expr_tm in let r, r2 = mk_taylor_functions pp r_tm in (* Unary operations *) if lhs = neg_op_real then Scale (r, ineg one), (fun pp domain_th -> taylor_interval_neg (r2 pp domain_th)) else (* Composite *) try let l, l2 = mk_taylor_functions pp lhs in if r_tm = x_var_real then l, l2 else Composite (l,r), (fun pp domain_th -> taylor_interval_compose pp l2 r2 domain_th) with Failure _ -> (* Binary operations *) if is_comb lhs then let op, l_tm = dest_comb lhs in let l, l2 = mk_taylor_functions pp l_tm in if op = add_op_real then Plus (l,r), mk_formal_binop taylor_interval_add l2 r2 else if op = sub_op_real then Minus (l,r), mk_formal_binop taylor_interval_sub l2 r2 else if op = mul_op_real then Product (l,r), mk_formal_binop taylor_interval_mul l2 r2 else if op = div_op_real then Product (l, Uni_compose (uinv,r)), mk_formal_binop taylor_interval_div_alt l2 r2 else failwith ("mk_taylor_functions: unknown binary operation "^string_of_term op) else failwith ("mk_taylor_functions: error: "^string_of_term lhs);; (**********************************) let rec result_size r = match r with | Result_false -> failwith "False result detected" | Result_glue (r1, r2) -> result_size r1 + result_size r2 | Result_pass _ -> 1;; let prove_result pp formal_f result lo_tm hi_tm = let pass_counter = ref 0 in let r_size = result_size result in let rec rec_prove result lo_tm hi_tm = let domain_th = mk_center_domain pp lo_tm hi_tm in match result with | Result_false -> failwith "False result" | Result_pass _ -> pass_counter := !pass_counter + 1; let th0 = taylor_cell_pass pp (formal_f pp domain_th) in let _ = report (sprintf "Result_pass: %d/%d" !pass_counter r_size) in th0 | Result_glue (r1, r2) -> let _, y_tm, _, _ = dest_cell_domain (concl domain_th) in let th1 = rec_prove r1 lo_tm y_tm and th2 = rec_prove r2 y_tm hi_tm in cell_pass_glue th1 th2 in rec_prove result lo_tm hi_tm;;
let CELL_PASS_IMP = 
prove (`cell_pass f (lo, hi) ==> interval_arith x (lo, w1) ==> interval_arith z (w2, hi) ==> !t. t IN real_interval [x, z] ==> f t < &0`,
REWRITE_TAC[cell_pass; interval_arith; IN_REAL_INTERVAL] THEN REPEAT STRIP_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REAL_ARITH_TAC);;
let prove_ineq1d pp expr_tm lo_tm hi_tm = let f, formal_f = mk_taylor_functions pp expr_tm in let lo_int, hi_int = build_eval_interval lo_tm pp, build_eval_interval hi_tm pp in let lo1_tm = (fst o dest_pair o rand o concl) lo_int and hi1_tm = (snd o dest_pair o rand o concl) hi_int in let lo, hi = float_of_float_tm lo1_tm, float_of_float_tm hi1_tm in let _ = report (sprintf "Starting verification from %f to %f" lo hi) in let result = let opt = { allow_sharp = false; iteration_count = 0; iteration_limit = 0; recursion_depth = 50 } in recursive_verifier (0, lo, hi, f, opt) in
let th0 = prove_result pp formal_f result lo1_tm hi1_tm in
    REWRITE_RULE[] (MATCH_MP (MATCH_MP (MATCH_MP CELL_PASS_IMP th0) lo_int) hi_int);;