let mem_stat () =
  let stat = Gc.stat() in
  let word = float_of_int (Sys.word_size / 8) in
  let free = float_of_int stat.Gc.free_words *. word /. 1024.0 in
  let total = float_of_int stat.Gc.heap_words *. word /. 1024.0 in
  let allocated = total -. free in
  let str = sprintf "allocated = %f (free = %f; total_size = %f; %f)\n" 
    allocated free total (free /. total) in
    print_string str;;


(* allocated = 85862 *)
mem_stat();;

Gc.set { (Gc.get()) with Gc.verbose = 0x05 };;
Gc.compact();;

(* allocated = 60243 *)
mem_stat();;


needs "../formal_lp/formal_interval/m_taylor_arith.hl";;
needs "../formal_lp/formal_interval/m_verifier.hl";;

let reset_all () =
  Arith_cache.reset_cache();
  Arith_cache.reset_stat();
  Arith_float.reset_cache();
  Arith_float.reset_stat();;

Arith_cache.print_stat();;
Arith_float.print_stat();;
reset_all();;


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

(* dummy functions *)

let eval_d0 i t1 t2 = failwith "eval_d0";;
let eval_dd0 i j t1 t2 = failwith "eval_dd0";;
let eval_0 t1 t2 = failwith "eval_0";;


(******************************)
(* real tests *)

let n = 6;;
let pp = 15;;
let pp = 8;;

(* delta_y *)
let delta_y_poly =
  let tm = (rand o concl o SPEC_ALL o REWRITE_RULE[Sphere.delta_x]) Sphere.delta_y in
    expr_to_vector_fun tm;;

(* delta_y4 *)
let delta_y4_poly =
  let tm = (rand o concl o SPECL[`y1*y1`; `y2*y2`; `y3*y3`; `y4*y4`; `y5*y5`; `y6*y6`]) Sphere.delta_x4 in
    expr_to_vector_fun tm;;

(* 4 * y1 * delta_y *)
let four_y1_delta_y_poly = 
  let x_var, tm = dest_abs delta_y_poly in
    mk_abs (x_var, mk_binop mul_op_real `&4` (mk_binop mul_op_real `(x:real^6)$1 * x$1` tm));;

(* - delta_y4 *)
let neg_delta_y4_poly =
  let x_var, tm = dest_abs delta_y4_poly in
    mk_abs (x_var, mk_comb (neg_op_real, tm));;


let (_, _, _, eval_neg_delta_y4), tf_neg_delta_y4 = 
  mk_verification_functions pp neg_delta_y4_poly false `&0`;;

let (_, _, _, eval_4y1_delta_y), tf_4y1_delta_y = 
  mk_verification_functions pp four_y1_delta_y_poly false `&0`;;

let (_, _, _, eval_pi2), tf_pi2 = 
  mk_verification_functions pp `\x:real^6. pi / &2` false `&0`;;


let (_, _, _, eval_pi2_plus), tf_pi2_plus = 
  mk_verification_functions pp `\x:real^6. pi / &2 - #1.893` false `&0`;;


(* dih_y *)
let tf_dih_y_hi =
  let denom = Uni_compose (uinv, Uni_compose (usqrt, tf_4y1_delta_y)) in
    Plus (tf_pi2_plus, Uni_compose (uatan, Product (tf_neg_delta_y4, denom)));;

let eval_dih_y_hi th =
  let inv, atn, sqrt, ( * ), ( + ) = 
    eval_m_taylor_inv n pp, eval_m_taylor_atn n pp, eval_m_taylor_sqrt n pp, 
    eval_m_taylor_mul n pp, eval_m_taylor_add n pp in
  let poly1 = eval_4y1_delta_y th and
      poly2 = eval_neg_delta_y4 th and
      pi2 = eval_pi2_plus th in
    pi2 + atn (poly2 * inv (sqrt (poly1)));;
    


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

let xx = `[&2; &2; &2; &2; &2; &2]` and
    zz = `[#2.52; #2.52; #2.52; #2.52; #2.52; #2.52]`;;

let xx1 = convert_to_float_list pp true xx and
    zz1 = convert_to_float_list pp false zz;;


(* Small domain *)
let xx_s = `[&2; &2; &2; &2; &2; &2]` and
    zz_s = `[#2.1; #2.1; #2.1; #2.1; #2.1; #2.1]`;;


let domain_th =
  let xx1_s = convert_to_float_list pp true xx_s and
      zz1_s = convert_to_float_list pp false zz_s in
    mk_m_center_domain n pp xx1_s zz1_s;;

reset_all();;
	
(* 10: 9.121 (pp = 15) *)
(* 300: 5.5 (pp = 8) *)
(* 100 (cached, float_cached, pp = 8): 2.68 *)
test 1 eval_dih_y_hi domain_th;;

Arith_cache.print_stat();;
Arith_float.print_stat();;


(***)

let xx_s = `[&2; &2; &2; &2; &2; &2]` and
    zz_s = `[#2.52; #2.1; #2.1; #2.1; #2.1; #2.1]`;;

let xx_s2 = `[&2; &2; &2; &2; &2; &2]` and
    zz_s2 = `[#2.52; #2.2; #2.2; #2.2; #2.2; #2.2]`;;


let xx1_s = convert_to_float_list pp true xx_s and
    zz1_s = convert_to_float_list pp false zz_s;;

let xx1_s2 = convert_to_float_list pp true xx_s2 and
    zz1_s2 = convert_to_float_list pp false zz_s2;;



let xx_float, zz_float, xx_s_float, zz_s_float, xx_s2_float, zz_s2_float =
  let pad = replicate 0.0 (8 - length (dest_list xx1)) in
    map float_of_float_tm (dest_list xx1) @ pad,
    map float_of_float_tm (dest_list zz1) @ pad,
    map float_of_float_tm (dest_list xx1_s) @ pad,
    map float_of_float_tm (dest_list zz1_s) @ pad,
    map float_of_float_tm (dest_list xx1_s2) @ pad,
    map float_of_float_tm (dest_list zz1_s2) @ pad;;



(***)

let c_dih_y_s = run_test tf_dih_y_hi xx_s_float zz_s_float false 0.0 true false true false 0.0;;
let c_dih_y_s2 = run_test tf_dih_y_hi xx_s2_float zz_s2_float false 0.0 true false true false 0.0;;
let c_dih_y0 = run_test tf_dih_y_hi xx_float zz_float false 0.0 true false false false 0.0;;

(* pass = 4 *)
result_stat c_dih_y_s;;
(* pass = 63 *)
result_stat c_dih_y_s2;;
(* pass = 4294, mono = 10 *)
result_stat c_dih_y0;;

reset_all();;
Gc.compact();;
(* 80642 *)
mem_stat();;

(* 10 (pp = 15): 38.418 *)
(* 300 (pp = 8): 22.289 *)
(* 100 (cached, float_cached, pp = 8): 12.229 *)
let _ =
  let start = Sys.time() in
  let result = m_verify_raw0 n pp (eval_d0, eval_dd0, eval_0, eval_dih_y_hi) c_dih_y_s xx1_s zz1_s in
  let finish = Sys.time() in
  let _ = report 
    (sprintf "Verification time: %f" (finish -. start)) in
    result;;

Arith_cache.print_stat();;
Arith_float.print_stat();;

(* 104321 *)
mem_stat();;
Gc.compact();;
(* 84679 *)
mem_stat();;

reset_all();;
Gc.compact();;
(* 80643 *)
mem_stat();;

(* 100 (cached, float_cached, pp = 8): 233.80 *)
let _ =
  let start = Sys.time() in
  let result = m_verify_raw0 n pp (eval_d0, eval_dd0, eval_0, eval_dih_y_hi) c_dih_y_s2 xx1_s2 zz1_s2 in
  let finish = Sys.time() in
  let _ = report 
    (sprintf "Verification time: %f" (finish -. start)) in
    result;;

Arith_cache.print_stat();;
Arith_float.print_stat();;

(* 135718 (after two runs)  *)
mem_stat();;
Gc.compact();;
(* 101048  *)
mem_stat();;

reset_all();;
Gc.compact();;
(* 81711  *)
mem_stat();;


(* 15202 *)
let _ =
  let start = Sys.time() in
  let result = m_verify_raw0 n pp (eval_d0, eval_dd0, eval_0, eval_dih_y_hi) c_dih_y0 xx1 zz1 in
  let finish = Sys.time() in
  let _ = report 
    (sprintf "Verification time: %f" (finish -. start)) in
    result;;