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);;