flyspeck_needs "../formal_lp/formal_interval/interval_m/recurse.hl";;
flyspeck_needs "../formal_lp/formal_interval/interval_m/recurse0.hl";;

module Verifier = struct

open Interval;;
open Univariate;;
open Line_interval;;
open Taylor;;
open Recurse;;



(****************************)
(* Interval functions for the native OCaml arithmetic *)

type int_fun =
  | F_int_var of int
  | F_int_const of interval
  | F_int_pow of int * int_fun
  | F_int_neg of int_fun
  | F_int_add of int_fun * int_fun
  | F_int_sub of int_fun * int_fun
  | F_int_mul of int_fun * int_fun;;


let ipow = Arith_misc.gen_pow imul Interval.one;;    


let eval_int_fun i_fun =
  fun x ->
    let rec eval_rec f =
      match f with
	| F_int_var i -> List.nth x (i - 1)
	| F_int_const int -> int
	| F_int_neg f1 -> ineg (eval_rec f1)
	| F_int_pow (n,f1) -> ipow n (eval_rec f1)
	| F_int_add (f1,f2) -> iadd (eval_rec f1) (eval_rec f2)
	| F_int_sub (f1,f2) -> isub (eval_rec f1) (eval_rec f2)
	| F_int_mul (f1,f2) -> imul (eval_rec f1) (eval_rec f2) in
      eval_rec i_fun;;

(**********************************)
let run_test f x z min_flag min_max allow_d convex_flag mono_pass_flag raw_int_flag eps =
  let pad = replicate 0.0 (8 - length x) in
  let xx = x @ pad and zz = z @ pad in
  let mone = mk_interval(-1.0,-1.0) in
  let neg_f = Scale(f, mone) in
  let ff = if min_flag then
    Plus(neg_f, Scale(unit,mk_interval(min_max, min_max)))
  else 
    Plus(f, Scale(unit, ineg (mk_interval(min_max, min_max)))) in
  let opt =  {
    only_check_deriv1_negative = false;
    is_using_dihmax =false;
    is_using_bigface126 =false;
    width_cutoff =0.05;
    allow_sharp =false;
    allow_derivatives =allow_d;
    iteration_count =0;
    iteration_limit =0;
    recursion_depth =200;
    mono_pass = mono_pass_flag;
    convex_flag = convex_flag;
    raw_int_flag = raw_int_flag;
    eps = eps;
  } in
    recursive_verifier(0,xx,zz,xx,zz,ff,opt);;



open Recurse0;;


let run_test0 f x z min_flag min_max allow_d convex_flag mono_pass_flag eps =
  let pad = replicate 0.0 (8 - length x) in
  let xx = x @ pad and zz = z @ pad in
  let mone = mk_interval(-1.0,-1.0) in
  let neg_f = Scale(f, mone) in
  let ff = if min_flag then
    Plus(neg_f, Scale(unit,mk_interval(min_max, min_max)))
  else 
    Plus(f, Scale(unit, ineg (mk_interval(min_max, min_max)))) in
  let opt =  {
    only_check_deriv1_negative = false;
    is_using_dihmax =false;
    is_using_bigface126 =false;
    width_cutoff =0.05;
    allow_sharp =false;
    allow_derivatives =allow_d;
    iteration_count =0;
    iteration_limit =0;
    recursion_depth =200;
    mono_pass = mono_pass_flag;
    convex_flag = convex_flag;
    raw_int_flag = true;
    eps = eps;
  } in
    recursive_verifier0(0,xx,zz,xx,zz,ff,opt);;



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




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

let domain_str x z =
  let s1 = map string_of_float x and
      s2 = map string_of_float z in
    sprintf "[%s], [%s]" (String.concat "; " s1) (String.concat "; " s2);;

let path_str p =
  String.concat "," (map (fun s, j -> sprintf "%s(%d)" s j) p);;


(* get_results0 *)
(* This function finds all subtrees of the given solution tree which can be
   veified immediately (no Result_pass_mono). These subtrees are added to
   the accumulator. Paths to the roots of all subtrees are also saved in
   the accumulator. The third returned value is a solution tree where all
   found subtrees are replaced with Result_pass_ref j, with j = #of the corresponding
   subtree in the accumulator (1-based) *)


let get_results0 path r acc =
  let dummy_tree = Result_false ([], []) in
  let is_ref r = match r with Result_pass_ref _ -> true | _ -> false in

  let rec get_rec path r acc =
    match r with
      | Result_mono (mono, r1) ->
	  let get_m m = (if m.decr_flag then "ml" else "mr"), m.variable in
	  let path' = rev_itlist (fun m l -> get_m m :: l) mono path in
	  let flag, acc', tree = get_rec path' r1 acc in
	    if flag then true, acc', dummy_tree
	    else false, acc', Result_mono (mono, tree)
      | Result_glue (j, convex_flag, r1, r2) ->
	  let s1, s2 = if convex_flag then "ml", "mr" else "l", "r" in
	  let p1, p2 = ((s1, j + 1) :: path), ((s2, j + 1) :: path) in
	  let flag1, acc1, tree1 = get_rec p1 r1 acc in
	  let flag2, acc', tree2 = get_rec p2 r2 acc1 in
	  let n = (length acc' + 1) in
	    if flag1 then
	      if flag2 then
		true, acc', dummy_tree
	      else if is_ref r1 then
		false, acc', Result_glue (j, convex_flag, r1, tree2)
	      else
		false, acc' @ [rev p1, r1], Result_glue (j, convex_flag, Result_pass_ref n, tree2)
	    else
	      if flag2 then
		if is_ref r2 then
		  false, acc', Result_glue (j, convex_flag, tree1, r2)
		else
		  false, acc' @ [rev p2, r2], Result_glue (j, convex_flag, tree1, Result_pass_ref n)
	      else
		false, acc', Result_glue (j, convex_flag, tree1, tree2)

      | Result_pass_mono _ -> false, acc, r
      | _ -> true, acc, dummy_tree in

    get_rec path r acc;;


    

(* transform_result *)


let transform_result x z r =
  (* get_domain *)
  (* Subdivides the given domain (x,z) according to the given path *)
  let domain_hash = Hashtbl.create 1000 in
  let find_hash, mem_hash, add_hash = 
    Hashtbl.find domain_hash, Hashtbl.mem domain_hash, Hashtbl.add domain_hash in

  let get_domain path =
    let n = length x in
    let table f = map f (0--(n - 1)) in
    let rec rec_domain (x, z) path hash =
      match path with
	| [] -> x, z
	| (s, j) :: ps ->
	    let hash' = hash^s^(string_of_int j) in
	      if mem_hash hash' then
		rec_domain (find_hash hash') ps hash'
	      else
		let j = j - 1 in
		let domain' =
		  if s = "l" or s = "r" then
		    let ( ++ ), ( / ) = up(); upadd, updiv in
		    let yj = (mth x j ++ mth z j) / 2.0 in
		    let delta b v = table (fun i -> if i = j && b then yj else mth v i) in
		      if s = "l" then 
			delta false x, delta true z
		      else
			delta true x, delta false z
		  else
		    if s = "ml" then
		      x, table (fun i -> if i = j then mth x i else mth z i)
		    else
		      table (fun i -> if i = j then mth z i else mth x i), z in
		let _ = add_hash hash' domain' in
		  rec_domain domain' ps hash' in
      rec_domain (x,z) path "" in

  (* sub_domain *)
  (* Verifies if interval [x',z'] SUBSET interval [x,z] *)
  let sub_domain (x',z') (x,z) =
    let le a b = itlist2 (fun a b c -> c & (a <= b)) a b true in
      le x x' & le z' z in

  (* transform_pass_mono *)
  (* Replaces all (Result_pass_mono m) with (Result_mono [m] (Result_ref j)) where
     j is the reference to the corresponding domain *)
  let transform_pass_mono x z domains r =
    let domains_i = zip domains (1--length domains) in

    let find_domain x' z' =
      try find (fun d, _ -> sub_domain (x', z') d) domains_i with Failure _ -> (x,z), -1 in

    let get_m m = (if m.decr_flag then "ml" else "mr"), m.variable in

    let rec rec_transform path r =
      match r with
	| Result_mono (mono, r1) ->
	    let path' = rev_itlist (fun m l -> get_m m :: l) mono path in
	      Result_mono (mono, rec_transform path' r1)
	| Result_glue (j, convex_flag, r1, r2) ->
	    let s1, s2 = if convex_flag then "ml", "mr" else "l", "r" in
	    let p1, p2 = ((s1, j + 1) :: path), ((s2, j + 1) :: path) in
	    let t1 = rec_transform p1 r1 in
	    let t2 = rec_transform p2 r2 in
	      Result_glue (j, convex_flag, t1, t2)
	| Result_pass_mono m -> 
	    let path' = rev (get_m m :: path) in
	    let x', z' = get_domain path' in
	    let _, i = find_domain x' z' in
	      (*	  let _ = report (sprintf "p = %s, d = %s, found: %d" 
			  (domain_str x' z') (path_str path') i) in *)
	      if i >= 0 then Result_mono ([m], Result_pass_ref (-i)) else r
	| _ -> r in

      rec_transform [] r in

  let rec transform acc r =
    let flag, rs, r' = get_results0 [] r acc in
      if flag then (rs @ [[], r])
      else
	let domains = map (fun p, _ -> get_domain p) rs in
	let r_next = transform_pass_mono x z domains r' in
	let _ = r_next <> r' or failwith "transform_result: deadlock" in
	  transform rs r_next in
    transform [] r;;






(*
let pp = 8;;
let n = 7;;
let xx = `[-- &1; -- &1; -- &1; -- &1; -- &1; -- &1; -- &1]` and
    zz = `[&1; &1; &1; &1; &1; &1; &1]`;;

let xx1 = convert_to_float_list pp true xx and
    zz1 = convert_to_float_list pp false zz;;
let xx_float = map float_of_float_tm (dest_list xx1) and
    zz_float = map float_of_float_tm (dest_list zz1);;

let eval0_magnetism, eval_magnetism, tf_magnetism = 
  mk_verification_functions pp magnetism_poly true magnetism_min;;
let c1 = run_test tf_magnetism xx_float zz_float false 0.0 true false true;;
let c0 = run_test0 tf_magnetism xx_float zz_float false 0.0 true false true;;
(* 77 *)
result_size c1;;
(* 121 *)
result_size c0;;

let r = transform_result xx_float zz_float c1;;
r;;
length r;;
let x = map (fun _, r ->
       match r with | Result_pass_ref j -> 1 | _ -> 0) r;;
itlist (+) x 0;;

map (fun _, r ->
       match r with 
	 | Result_pass _ -> "pass"
	 | Result_mono _ -> "mono"
	 | Result_glue _ -> "glue"
	 | Result_false _ -> "_|_"
	 | Result_pass_mono _ -> "pass_mono"
	 | Result_pass_ref _ -> "ref"
    ) r;;


length it;;
transform_result xx_float zz_float c0;;

let pp = 5;;
let n = 2;;

let poly_tm = expr_to_vector_fun `x1 pow 2 + x2 pow 2`;;

let xx = `[-- &6; -- &6]` and
    zz = `[&1; &2]`;;

let xx1 = convert_to_float_list pp true xx and
    zz1 = convert_to_float_list pp false zz;;
let xx_float = map float_of_float_tm (dest_list xx1) and
    zz_float = map float_of_float_tm (dest_list zz1);;

let eval0_poly, eval_poly, tf_poly = mk_verification_functions pp poly_tm true `-- &1`;;
let c1 = run_test tf_poly xx_float zz_float false 0.0 true false true;;

transform_result xx_float zz_float c1;;

let c0 = run_test0 tf_poly xx_float zz_float false 0.0 true false true;;
(* 22 *)
result_size c1;;
(* 38 *)
result_size c0;;
m_verify_raw n pp eval0_poly eval_poly c1 xx1 zz1;;
m_verify_raw n pp eval0_poly eval_poly c1 xx1 zz1;;



transform_result xx_float zz_float c1;;

let flag, rs1, r1 = get_results0 [] c1 [];;
rs1;;
r1;;
let domains = map (fun p, _ -> get_domain xx_float zz_float p) rs1;;
let r2 = transform_pass_mono xx_float zz_float domains r1;;

let flag, rs3, r3 = get_results0 [] r2 rs1;;
rs3;;
r3;;
let domains = map (fun p, _ -> get_domain xx_float zz_float p) rs3;;	
let r4 = transform_pass_mono xx_float zz_float domains r3;;

let flag, rs5, r5 = get_results0 [] r4 rs3;;	
rs5;;
r5;;



test 1000 (get_results0 [] c1) [];;
test 1000 (map (fun p, _ -> get_domain xx_float zz_float p)) rs1;;
test 1000 (transform_pass_mono xx_float zz_float domains) r1;;

let flag, rs3, r3 = get_results0 [] r2 rs1;;
let domains = map (fun p, _ -> get_domain xx_float zz_float p) rs3;;	
let r4 = transform_pass_mono xx_float zz_float domains r3;;

let flag, rs5, r5 = get_results0 [] r4 rs3;;	
*)
	


let result_stat result =
  let pass = ref 0 and
      mono = ref 0 and
      glue = ref 0 and
      pass_mono = ref 0 and
      pass_raw = ref 0 and
      glue_convex = ref 0 in
    
  let rec count r =
    match r with
      | Result_false _ -> failwith "False result"
      | Result_pass (flag, _, _) ->
	  pass := !pass + 1;
	  if flag then pass_raw := !pass_raw + 1 else ()
      | Result_pass_mono _ -> pass_mono := !pass_mono + 1
      | Result_mono (_, r1) -> mono := !mono + 1; count r1
      | Result_glue (_, flag, r1, r2) ->
	  glue := !glue + 1;
	  if flag then glue_convex := !glue_convex + 1 else ();
	  count r1; count r2 in

  let _ = count result in
  let s = sprintf "pass = %d (pass_raw = %d)\nmono = %d\nglue = %d (glue_convex = %d)\npass_mono = %d"
    !pass !pass_raw !mono !glue !glue_convex !pass_mono in
    report s;;


let result_p_stat glue_flag p_result =
  let p_table = Hashtbl.create 10 in
  let add1 p =
    let c = if Hashtbl.mem p_table p then Hashtbl.find p_table p else 0 in
      Hashtbl.replace p_table p (succ c) in

  let rec count r =
    match r with
      | P_result_ref _ -> ()
      | P_result_pass (pp, _) -> add1 pp.pp
      | P_result_mono (pp, _, r1) -> add1 pp.pp; count r1
      | P_result_glue (pp, _, _, r1, r2) ->
	  if glue_flag then add1 pp.pp else ();
	  count r1; count r2 in

  let _ = count p_result in
  let s = Hashtbl.fold 
    (fun p c s -> (sprintf "p = %d: %d\n" p c) ^ s) p_table "" in
    report s;;
				  

end;;