(* ========================================================================== *)
(* FLYSPECK - BOOK FORMALIZATION                                              *)
(*                                                                            *)
(* Chapter: nonlinear inequalities                                            *)
(* Author:  Thomas Hales     *)
(* Date: 2010-08-30                                                           *)
(* ========================================================================== *)


(*
Input is a nonlinear inequality.
preprocess_split_idq generates cases and does general preprocessing.
Output is executable interval arithmetic code in C++.
*)

(*
flyspeck_needs "general/sphere.hl";;
flyspeck_needs "nonlinear/ineq.hl";;
flyspeck_needs "nonlinear/lemma.hl";;
*)

module Optimize = struct

open Lib;;
open Nonlinear_lemma;;


let length = List.length;;


let svn_version = Flyspeck_lib.svn_version;;

let ineq = Sphere.ineq;;
let all_forall = Sphere.all_forall;;
let add = Ineq.add;;


let nub = Flyspeck_lib.nub;;

let join_comma = Flyspeck_lib.join_comma;;
let join_space = Flyspeck_lib.join_space;;
let join_lines = Flyspeck_lib.join_lines;;

let string_of_num' = Flyspeck_lib.string_of_num';;
let dest_decimal = Flyspeck_lib.dest_decimal;;

let strip_let_tm = Nonlinear_lemma.strip_let_tm;;

let output_filestring = Flyspeck_lib.output_filestring;;

let BY = Hales_tactic.BY;;
let unlist = Hales_tactic.unlist;;

(* from glpk_link.ml *)

let load_and_close_channel_true = Flyspeck_lib.load_and_close_channel_true;;

let load_file = Flyspeck_lib.load_file;;

(* ========================================================================= *)
(*    GENERAL PROCEDURES                                                  *)
(* ========================================================================= *)

(* The marker NONLIN is used in mk_all_ineq.hl, 
   where it is replaced with a constant equivalent to the conjunction of all the
   nonlinear ineqs.   *)

let NONLIN = new_definition `NONLIN = T`;;
let FROZEN_REWRITE_TAC ths = 
   if (ths=[]) then REWRITE_TAC [] else
     let th = end_itlist CONJ ths in
      FREEZE_THEN (fun t -> REWRITE_TAC[t]) th;;

let quoted s = let q = "\"" in (q^s^q);;

let paren s = "("^ s ^")";;

let real_ty = `:real`;;
let mk_x i = mk_var("x"^string_of_int i,real_ty);;

let x9list =  map mk_x (1--9);;
let x6list = map mk_x (1--6);;
let x_backlist = map mk_x [7;2;3;4;8;9];;

let xspec = SPECL x6list;;


let dest_ineq ineq = 
  let t = snd(strip_forall ineq) in
  let (vs,i) = dest_comb t in
  let (_,vs) = dest_comb vs in
  let vs = dest_list vs in
  let vs = map (fun t -> let (a,b) = dest_pair t in (a,dest_pair b)) vs in
  let vs = map (fun (a,(b,c)) -> (a, b, c)) vs in
    (t,vs,disjuncts i);;

(* renames variables as x1,x2,x3,.... *)

let (X_RENAME_TAC:tactic) = 
   fun  gl  ->
     let (_,loc,_) = dest_ineq (goal_concl gl) in
     let loc' = map (fun (_,b,_) -> b) loc in
     let xvar = map mk_x (1 -- (length loc')) in
       EVERY (map SPEC_TAC (zip loc' xvar)) gl ;;


(* ========================================================================= *)
(*    SPLITTING INEQUALITIES AT h0                                           *)
(* ========================================================================= *)

let get_split_tags idq = 
  let ts = idq.tags in
  let rec gs = (function
    | [] -> []
    | (Split a::_) -> a 
    | _ :: rs -> gs rs) in
    gs ts;;

let split_2h0 = 
  
let split_interval = 
prove(
    `! a b (y:real). (a <= y /\ y <= b) ==>
      ((a <= y /\ y <= m) \/ (m <= y /\ y <= b) )
      `,
  REAL_ARITH_TAC) in
    INST [(`&2 * h0`,`m:real`)] split_interval;;
let split_2h0_strict = 
  
let split_interval = 
prove(
    `! a b (y:real). (a <= y /\ y <= b) ==>
      ((a <= y /\ y <= m) \/ (m < y /\ y <= b) )
      `,
  REAL_ARITH_TAC) in
    INST [(`&2 * h0`,`m:real`)] split_interval;;

(*  changed 2012-06-20
let SPLIT_H0_TAC pos =  
 (REPEAT GEN_TAC) THEN
 (REWRITE_TAC[ineq_expand6]) THEN
 (REPLICATE_TAC (pos) DISCH_TAC) THEN
 (DISCH_THEN (fun t -> MP_TAC (MATCH_MP split_2h0 t))) THEN
 DISCH_THEN DISJ_CASES_TAC    THEN 
 (REPEAT (POP_ASSUM MP_TAC)) THEN 
 (REWRITE_TAC[GSYM ineq_expand6]);;
*)


let h0cutA = 
prove_by_refinement(
  `!y. (y <= &2 * h0) ==> ((h0cut y) = &1)`,

  (* {{{ proof *)
  [
  MESON_TAC[Sphere.h0cut]
  ]);;
  (* }}} *)

let h0cutB = 
prove_by_refinement(
  `!y. (&2 * h0 < y) ==> ((h0cut y) = &0)`,

  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.h0cut];
  REPEAT STRIP_TAC;
  BY(COND_CASES_TAC THEN REPEAT (FIRST_X_ASSUM MP_TAC) THEN REAL_ARITH_TAC)
  ]);;
  (* }}} *)

let h0cutC = 
prove_by_refinement(
  `!y. (y <= &2 * hminus) ==> (h0cut y = &1)`,

  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  MATCH_MP_TAC h0cutA;
  MP_TAC Nonlinear_lemma.hminus_lt_h0;
  POP_ASSUM MP_TAC;
  BY(REAL_ARITH_TAC)
  ]);;
  (* }}} *)


let SPLIT_H0_TAC pos =  
  let destrict = REAL_ARITH `a < y /\ y <= b ==> a <= y /\ y <= b` in
 (REPEAT GEN_TAC) THEN
 (REWRITE_TAC[ineq_expand6]) THEN
 (REPLICATE_TAC (pos) DISCH_TAC) THEN
 (DISCH_THEN (fun t -> MP_TAC (MATCH_MP split_2h0_strict t))) THEN
 DISCH_THEN DISJ_CASES_TAC    THEN 
 ASM_SIMP_TAC[h0cutA;h0cutB;h0cutC] THENL
   [ALL_TAC;POP_ASSUM ( ASSUME_TAC o (MATCH_MP destrict))] THEN
 (REPEAT (POP_ASSUM MP_TAC)) THEN 
 (REWRITE_TAC[GSYM ineq_expand6]);;


(* WARNING: This disturbs the goal *)

let split_h0 (ineq,ss) = 
  let current_goals() = 
    map (fun (_,w) -> w) 
      ((fun (_,b,_) -> b) (hd (!current_goalstack))) in 
  let split_h0_term i inq = 
    let _ = g inq in
    let _ = e(SPLIT_H0_TAC i) in
      map all_forall (current_goals()) in 
  let _ = (length ss > 0)  or failwith "empty split" in
  let ineql = split_h0_term (hd ss) ineq in
   map (fun t-> (t, tl ss)) ineql;;

let rec split_all_h0 = function
  | [] -> []
  | (i,[])::rs -> i :: split_all_h0 rs
  | r::rs -> split_all_h0 (split_h0 r @ rs);;



(* ========================================================================== *)
(*    PARSING INEQUALITIES                                                    *)
(* ========================================================================== *)



let dest_nonlin t = 
  let (_,r,il) = dest_ineq t in 
  let p1 (a,_,_) = a in
  let p2 (_,b,_) = b in
  let p3 (_,_,c) = c in
  let dest x = try dest_binop `(real_lt)` x with Failure _ -> dest_binop `(real_le)` x in
  let (iis,zzs) = unzip( map (dest) il) in
  let zz = nub zzs in
  let _ = (zz = [`&0`]) or failwith "zero expected" in
  let _ = ((map p2 r = x6list) or (map p2 r = x9list)) or failwith "x1..x6 or x1..x9 expected" in
  (map p1 r, map p3 r,iis);;

 (* these names don't change in cpp interval code *)

let idem_assoc = [];; (* make it the default *)

(* these names can change in cpp (hol_name,cpp_name).
    Don't generate function defs for these.
    "NOT_IMPLEMENTED" is a special tag to suppress the
    production of test routines. We would be better
    off using an option. *)

let native_fun = [
  ("unit6","Lib::unit");
  ("two6","Lib::two6");
  ("proj_x1","Lib::x1");
  ("proj_x2","Lib::x2");
  ("proj_x3","Lib::x3");
  ("proj_x4","Lib::x4");
  ("proj_x5","Lib::x5");
  ("proj_x6","Lib::x6");
  ("proj_y1","Lib::y1");
  ("sqrt_x1","Lib::y1");
  ("delta_x","Lib::delta_x");
  ("sol_x","Lib::sol_x");
  ("dih_x","Lib::dih_x");
  ("rad2_x","Lib::rad2_x");
  ("rho_x","NOT_IMPLEMENTED");
  ("delta_y","NOT_IMPLEMENTED");
  ("dih_y","NOT_IMPLEMENTED");
  ("sol_y","NOT_IMPLEMENTED");
  ("delta_x4","Lib::delta_x4");
  ("edge_flat2_x","Lib::edge_flat2_x");
  ("halfbump_x1","Lib::halfbump_x1");
  ("ups_126","Lib::ups_126");
  ("eta2_126","Lib::eta2_126");
  ("compose6","Function::compose");
  ("constant6","Lib::constant6");
  ("uni","Lib::uni");
  ("sqrt","univariate::i_sqrt");
  ("acs","univariate::i_acos");
  ("asn","univariate::i_asin");
  ("atn","univariate::i_atan");
  ("cos","univariate::i_cos");
  ("sin","univariate::i_sin");
  ("pow2","univariate::i_pow2");
  ("matan","univariate::i_matan");
  ("gchi","Lib::i_gchi");
  ("lfun","Lib::i_lfun");
  ("rho","Lib::i_rho");
  ("flat_term_x","Lib::i_flat_term_x");
  ("sqn","univariate::i_sqn");
  ("promote1_to_6","Lib::promote1_to_6");
  ("gamma2_x_div_azim_v2","Lib::i_gamma2_x_div_azim_v2");
  ("num1","Lib::num1");
];;


(* and for these 
gradually eliminate this list.
The default will be to have the same names.
*)

let cpp_assoc = 
  let pref = "" in
  let p (a,b) = (a, (pref^b)) in map p 
  (idem_assoc @ [
  ("sqrt_x2","proj_y2");
  ("sqrt_x3","proj_y3");
  ("sqrt_x4","proj_y4");
  ("sqrt_x5","proj_y5");
  ("sqrt_x6","proj_y6");
  ("promote_pow2","x1square");
  ("promote_pow3","x1cube");

 ]);;

let cpp_fn_name s = 
  try Lib.assoc s native_fun 
  with Failure _ ->
    try Lib.assoc s cpp_assoc
    with Failure _ -> s;;

let cpp0_assoc = 
  [("sqrt8","sqrt8");("pi","pi");("const1","const1");
   ("hminus","hminus");("sqrt3","sqrt3");
   ("arc_hhn","arc_hhn");
   ("adodec","aStrongDodec");
   ("bdodec","bStrongDodec");
   ("ydodec","yStrongDodec");
   ];;

let cpp_of_constant s  = 
  try (assoc s cpp0_assoc) with 
      Failure _ -> failwith (s^" find: cpp_of_constant") ;;

let cpp_of_fun to_string xlist s xs = 
  let (arg,xss) = chop_list (length xs - length xlist) xs in
  let _ = (xss = xlist) or failwith ("standard x list expected in "^s) in
  let hd =  (cpp_fn_name s) in
  let arg_s = map (paren o to_string) arg in
  let p = if (length arg > 0) then paren else I in
    (hd ^ " " ^ p(join_comma arg_s));;

(* 2013-08-14, changed for AWS compiler. was interval::interval( ... ) *)
let i_mk s = "interval("^quoted s ^")";;

let i_mk2 t = (* treat small integers exactly *)
  let n = dest_numeral t in
  let s = if (Num.is_integer_num n) && (Num.sign_num n >= 0) && (n </ Int 100) 
     then string_of_num n else string_of_num' n in
    i_mk s;;

let cpp_string_of_term = 
  let rec soh t = 
    if is_var t then fst (dest_var t) else
      let (f,xs) = strip_comb t in
      let _  = (is_const f) or 
	failwith ("constant expected:" ^ string_of_term f) in
      let ifix i = let [a;b] = xs in paren(soh a ^ " " ^ i ^ " " ^ soh b) in
	match fst (dest_const f) with
	  | "real_add" -> ifix "+"
	  | "real_mul" -> ifix "*"
	  | "real_div" -> ifix "/"
	  | "\\/" -> ifix "\\/"
	  | "real_neg" -> let [a] = xs in "(-" ^ soh a ^ ")"
	  | "real_of_num" -> let [a] = xs in i_mk2( a)  
	  | "NUMERAL" -> let [a] = xs in string_of_num' (dest_numeral t)
	  | "<" -> let [a;b] = xs in paren(soh a ^ " < " ^ soh b)
	  | ">" -> let [a;b] = xs in paren(soh a ^ " > " ^ soh b)
	  | "+" -> let [a;b] = xs in paren(soh a ^ " + " ^ soh b)
	  | "*" -> let [a;b] = xs in paren(soh a ^ " * " ^ soh b)
	  | "DECIMAL" ->  i_mk(string_of_num' (dest_decimal t))
	  | s -> if (xs = []) 
	    then "("^cpp_of_constant s^")"  
	    else paren(cpp_of_fun soh x6list s xs) in
    fun t -> 
      try (soh t) 
      with Failure s -> failwith (s^" .......   "^string_of_term t);;

(* generation of cpp code *)

let cpp_template_taylor c (i,s) = Printf.sprintf 
"       Function F%s%d = %s;" c i s;;

let cpp_template_t c iis = 
  join_lines 
    (map (cpp_template_taylor c) (zip (1--(length iis)) iis));;

let cpp_template_F c i = Printf.sprintf "&F%s%d" c i;;

let cpp_template_Fc c len = join_comma 
  (map (cpp_template_F c) (1-- len));;

let rec delta126min = function 
    | [] -> -1.0
    | Delta126min t :: _ -> t
    | _ :: a -> delta126min a;;

let rec widthCut = function
  | [] -> (false,0.0)
  | Widthcutoff t::_ -> (true,t)
  | _:: a -> widthCut a;;

let rec delta126max = function 
    | [] -> -1.0
    | Delta126max t :: _ -> t
    | _ :: a -> delta126max a;;

let rec delta135min = function 
    | [] -> -1.0
    | Delta135min t :: _ -> t
    | _ :: a -> delta135min a;;

let rec delta135max = function 
    | [] -> -1.0
    | Delta135max t :: _ -> t
    | _ :: a -> delta135max a;;

let cpp_template_gen = Printf.sprintf "
 const char svn[] = %s;
 const char ineq_id[] = %s;

 int testRun() // autogenerated code
	{
	interval tx[6]={%s};
	interval tz[6]={%s};
	domain x = domain::lowerD(tx);
	domain z = domain::upperD(tz);
        domain x0=x;
        domain z0=z;
        %s
        const Function* I[%d] = {%s}; // len ...
        cellOption opt;
        opt.allowSharp = %d; // sharp
        opt.onlyCheckDeriv1Negative = %d; // checkderiv
        %s // other options.
	return  prove::recursiveVerifier(0,x,z,x0,z0,I,%d,opt); // len
	}";;

let mk_cpp_proc t s tags = 
  let sharp = if  mem Sharp tags then 1 else 0 in
  let checkderiv = if  mem Onlycheckderiv1negative tags then 1 else 0 in
  let ifd b s = if b then s else "" in
  let setrad2 = ifd (mem Set_rad2 tags) "\topt.setRad2 = 1;\n" in
  let (b,f) = widthCut tags in
  let sWidth = ifd b (Printf.sprintf "\topt.widthCutoff = %8.16f;\n" f) in 
  let d126min = delta126min tags in
  let s126min = ifd (d126min > -1.0) 
    (Printf.sprintf  "\topt.delta126Min = %8.16f;\n" d126min) in
  let d126max = delta126max tags in
  let s126max = ifd (d126max > -1.0) 
    (Printf.sprintf "\topt.delta126Max = %8.16f;\n" d126max) in
  let d135min = delta135min tags in
  let s135min = ifd (d135min > -1.0) 
    (Printf.sprintf "\topt.delta135Min = %8.16f;\n" d135min) in
  let d135max = delta135max tags in
  let s135max = ifd (d135max > -1.0)  
    (Printf.sprintf "\topt.delta135Max = %8.16f;\n" d135max) in
  let s206A = ifd (s="2065952723 A1" ) 
    (Printf.sprintf "\topt.strategy206A=1;\n") in 
  let disallowderiv = ifd  (mem Disallow_derivatives tags)
    "\topt.allowDerivatives=0;\n" in    
  let others = setrad2 ^ s126min ^ s126max ^ s135min ^ s135max ^ 
    sWidth ^ s206A ^ disallowderiv in
  let c = map cpp_string_of_term in
  let f (x,y,z) = (c x,c y,c z) in
  let (aas,bbs,iis) = f (dest_nonlin t) in
  let len = length iis in
  let sq = quoted s in
  let svn = (quoted(svn_version())) in
  let jaas = join_comma aas in
  let jbbs = join_comma bbs in
    cpp_template_gen svn sq jaas jbbs (cpp_template_t "" iis) 
      len (cpp_template_Fc "" len) sharp  checkderiv others len;;

(* quad clusters *)

let has_subterm sub tm = 
  can (find_term ((=) sub)) tm;; 

let has_cross_diag = has_subterm `quad_cross_diag2_x`;;

let has_delta_top = has_subterm `delta_top_x`;;

let is_quad_cluster tags = (can (find (function Quad_cluster _ -> true | _ -> false)) tags);;

let quad_margin tags = (function Quad_cluster t -> t | _ -> 0.0) 
  (find ((function Quad_cluster _ -> true | _ -> false))  tags);;

let check_quad_partition_term (r,s,t) tm = 
  let tm2 = list_mk_binop `(+)` (r @ s @ t) in
  let t = mk_eq (tm ,tm2) in 
      (REAL_ARITH) t;;

let partition_quad_term tm = 
  let rec split tm = 
    if (has_cross_diag tm) then ([],[],[tm]) else
    if (has_delta_top tm) then ([],[],[tm]) else
      if (subset (frees tm) x6list) then ([tm],[],[]) else
	if (subset (frees tm) x_backlist) then ([],[tm],[]) else
	  let (a1,a2) = dest_binop `(+)` tm in
	  let  (r1,s1,t1) = split a1 in
	  let (r2,s2,t2) = split a2  in
	    (r1 @ r2, s1 @ s2, t1 @ t2) in
  let (r,s,t) = split tm in
  let _ = check_quad_partition_term (r,s,t) tm in
  let v =  vsubst [`x1:real`,`x7:real`;  `x5:real`,`x8:real`; `x6:real`,`x9:real`] in
  let w = function | [] -> 
   `unit6 (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) * #0.0` | xl -> list_mk_binop `(+)` xl in
   (w r,w(map v s), t);;


let cppq_template_gen = Printf.sprintf "
 char* svn = %s;
 char* ineq_id = %s;

 int testRun()  // quad cluster case, autogenerated code
        {
	interval txA[6]={%s};
	interval tzA[6]={%s};
	domain xA = domain::lowerD(txA);
	domain zA = domain::upperD(tzA);
	interval txB[6]={%s};
	interval tzB[6]={%s};
	domain xB = domain::lowerD(txB);
	domain zB = domain::upperD(tzB);
        // Declare FA, FB...
        %s
        %s
	const Function* IA[%d] = {%s};
	const Function* IB[%d] = {%s};
	cellOption opt;
        %s // options.
	return prove::recursiveVerifierQ(0,xA,xB,zA,zB,IA,IB,%d,opt);
        }";;

let get_cross_diag_squared is9 =
  let id9' = filter (fun (_,_,c)-> length c >0 && has_cross_diag (hd c)) is9 in
  let is9'' = (fun (_,_,[c])::_ -> c ) id9' in
  let tm = term_match x9list `quad_cross_diag2_x x1 x2 x3 x4 x5 x6 x7 x8 x9 +
    unit6 x1 x2 x3 x4 x5 x6 * t` is9'' in
  let t =   (fun ([],[(x,_)],[]) -> x) tm in
    mk_binop `( * )` t t;;  (* t is negative, but it gets squared *)

let get_delta_top_squared is9 =
  let id9' = filter (fun (_,_,c)-> length c >0  && has_delta_top (hd c)) is9 in
  let is9'' = (fun (_,_,[c])::_ -> c ) id9' in
  let tm = term_match x9list `delta_top_x t x1 x2 x3 x4 x5 x6 x7 x8 x9` is9'' in
  let t =   (fun ([],[(x,_)],[]) -> x) tm in
    mk_binop `( * )` t t;;  (* gets squared *)

let cde_template = 
  Printf.sprintf "opt.crossDiagMinEnclosed = interMath::inf(%s);\n";;

let cdt_template = 
  Printf.sprintf "opt.crossDiagMinDelta = interMath::inf(%s);\n";;

let mk_cppq_proc t s tags =   
  let svn = (quoted(svn_version())) in
  let ineq_id = quoted s in
  let cpp = cpp_string_of_term in
  let (x,z,is) = dest_nonlin t  in
  let is' = map partition_quad_term is in
  let (is9,is6) = partition (fun (_,_,c) -> length c > 0)   is' in
  let (is6A,is6B) = unzip (map (fun (a,b,_) -> (cpp a,cpp b)) is6) in
  let len = length  is6 in
  let xs = map cpp x in
  let zs = map cpp z in
  let nth = List.nth in
  let a_part u = [nth u 0;nth u 1;nth u 2;nth u 3; nth u 4;nth u 5] in
  let b_part u = [nth u 6;nth u 1;nth u 2;nth u 3;nth u 7; nth u 8] in
  let xA = join_comma(a_part xs) in
  let xB = join_comma(b_part xs) in
  let zA = join_comma(a_part zs) in
  let zB = join_comma(b_part zs) in
  let opt_cross = 
    try 
      (  cde_template (cpp ( get_cross_diag_squared is9) ))
    with _ -> "" in
  let opt_delta_top = 
    try 
      (  cdt_template (cpp ( get_delta_top_squared is9) )) 
    with _ -> "" in
  let backsym = 
    if (mem Dim_red_backsym tags) 
    then "\topt.dimRedBackSym = 1;\n" else "" in
  let sWidth = 
    let (b,f) = widthCut tags in
      if b then Printf.sprintf "\topt.widthCutoff = %8.16f;\n" f else ""  in
  let margin = Printf.sprintf "\topt.margin = %8.16f;\n\n" (quad_margin tags) in
  let options = opt_cross ^ opt_delta_top ^ backsym^margin^sWidth in
  cppq_template_gen svn ineq_id xA zA xB zB 
    (cpp_template_t "A" is6A) (cpp_template_t "B" is6B) 
    len (cpp_template_Fc "A" len) len (cpp_template_Fc "B" len) options len ;;

(*
let t = snd(top_goal());;
let s = "test";;
let tags = [Quad_cluster 3.0];;
*)

(*     next: put together header, proc, tail and run *)

let tmpfile = flyspeck_dir^"/../interval_code/test_auto.cc";;

let cpp_header() = join_lines (load_file  (flyspeck_dir^"/../interval_code/generic_head.txt"));;

let cpp_tail() = join_lines (load_file  (flyspeck_dir^"/../interval_code/generic_tail.txt"));;

let mkfile_cpp  t s tags  = 
  output_filestring tmpfile 
   (join_lines [cpp_header(); (mk_cpp_proc t s tags);cpp_tail()]);;

let mkfile_cppq  t s tags =
  output_filestring tmpfile 
   (join_lines [cpp_header(); (mk_cppq_proc t s tags);cpp_tail()]);;

let compile_cpp = 
  let err = Filename.temp_file "cpp_err" ".txt" in
    fun () -> 
      let e = Sys.command("cd "^flyspeck_dir^"/../interval_code; make test_auto >& "^err) in
      let _ =   (e=0) or (let _ = Sys.command ("cat "^err) in failwith "compiler error") in
	();;

 let execute_interval ex tags s testineq = 
  let interval_dir = flyspeck_dir^"/../interval_code" in
  let quad_cluster = is_quad_cluster tags  in
  let _ = if quad_cluster then
   mkfile_cppq testineq s tags
  else 
      mkfile_cpp testineq s tags in
  let _ = compile_cpp() in 
  let _ = (not ex) or (0=  Sys.command(interval_dir^"/test_auto")) or failwith "interval execution error" in
    ();;


(* ========================================================================== *)
(*    MEGA PREP    TACTICS                                                    *)
(* ========================================================================== *)

(* *************************************************************************** *)
(* PRELIM *)
(* *************************************************************************** *)

let macro_expand = (
   [
     Sphere.x1_delta_y;Sphere.delta4_squared_y;
     vol4f_palt;Nonlinear_lemma.y_of_x_e;
     gamma4fgcy_alt;
     Sphere.vol3r;Sphere.vol2f;Sphere.gamma4f;
     Sphere.gamma3f;  GSYM Nonlinear_lemma.quadratic_root_plus_curry;
     REAL_MUL_LZERO;
     REAL_MUL_RZERO;FST;SND;
     Sphere.node2_y;Sphere.node3_y] @ 
     (!Ineq.dart_classes));;

let PRELIM_REWRITE_TAC = EVERY[
  (REWRITE_TAC[GSYM Sphere.rad2_y]) ;
  (REWRITE_TAC(macro_expand)) ;
];;

(* *************************************************************************** *)
(* BRANCH *)
(* *************************************************************************** *)

(* take care of branching *)

let BRANCH_TAC = EVERY[
  REWRITE_TAC[REAL_ARITH `x / &1 = x /\ &0 * x = &0 /\ &0 +x = x`];
  REWRITE_TAC[Sphere.gamma4f;vol4f_lmfun;Sphere.gamma3f;];
  REWRITE_TAC[ineq_expand6];
  DISCH_TAC;
  REPEAT GEN_TAC;
  REPEAT DISCH_TAC;
  ASSUM_LIST (let rec r = function | [] -> ALL_TAC | th::ths -> (MP_TAC (REPEAT_RULE (OR_RULE (MATCH_MP pathL_bound) (MATCH_MP pathR_bound)) th)) THEN r ths in r);
  REWRITE_TAC[];
  SIMP_TAC[c2001;
	   gcy_low;gcy_low_const;gcy_low_hminus;gcy_high;gcy_high_hplus;
           h0_lt_gt;
           lmdih3_ldih3;lmdih5_ldih5;lmdih_ldih;lmdih5_0;lmdih3_0;lmdih0;
           (* Oct 28, 2010:*) vol3f_lm0;vol3f_lmln;
           (* nov23 *) lmfun0;lmfun_lfun;hm0;
	 (* may 26, 2011. macro_expand was moved to prelim tac, 
	    so additional simps are needed. *)
	 lmdih1_0';lmdih3_0';lmdih5_0';
	 lmdih_ldih';lmdih3_ldih3';lmdih5_ldih5';
   ];
  SIMP_TAC[lmfun_lfun;lmfun0];
  REWRITE_TAC[REAL_ARITH `&0 * x = &0 /\ &0 + x = x`];
  REPEAT (DISCH_THEN (fun t-> ALL_TAC));
  REPEAT (POP_ASSUM MP_TAC);
  REWRITE_TAC[GSYM ineq_expand6];
  DISCH_TAC;
  EVERY (map SPEC_TAC [(`y6:real`,`y6:real`);(`y5:real`,`y5:real`);(`y4:real`,`y4:real`);(`y3:real`,`y3:real`);(`y2:real`,`y2:real`);(`y1:real`,`y1:real`)]);
  POP_ASSUM MP_TAC;
  REWRITE_TAC[Sphere.pathL;Sphere.pathR]; (*  moved May 24, 2011 *)
];;

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

(*
  let SQRT_SQRT_TAC = 
EVERY[
  REPEAT DISCH_TAC;
  SUBGOAL_THEN `sqrt x1 * sqrt x1 = x1 /\ sqrt x2 * sqrt x2 = x2 /\ sqrt x3 * sqrt x3 = x3 /\ sqrt x4 * sqrt x4 = x4 /\ sqrt x5 * sqrt x5 = x5 /\ sqrt x6 * sqrt x6 = x6` (unlist REWRITE_TAC) THENL[ ASM_MESON_TAC[sq_pow2];ALL_TAC] THEN
  REPEAT (POP_ASSUM MP_TAC);
  DISCH_TAC;
];;
*)

  let SQRT_SQRT_TAC = EVERY[
  REPEAT DISCH_TAC;
  SUBGOAL_THEN `sqrt x1 * sqrt x1 = x1 /\ sqrt x2 * sqrt x2 = x2 /\ sqrt x3 * sqrt x3 = x3 /\ 
      sqrt x4 * sqrt x4 = x4 /\ sqrt x5 * sqrt x5 = x5 /\ sqrt x6 * sqrt x6 = x6` 
      (unlist REWRITE_TAC) 
  THENL [ REPEAT (FIRST_X_ASSUM (MP_TAC o (MATCH_MP sq_pow2))) THEN 
	    REPEAT DISCH_TAC THEN BY(ASM_REWRITE_TAC[]);ALL_TAC] THEN
  REPEAT (POP_ASSUM MP_TAC);
  DISCH_TAC;];;



let X_SQRT_COMPOUND_ORDER = 
[`norm2hh (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
     norm2hh_x x1 x2 x3 x4 x5 x6 `;  
`rad2_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
     rad2_x x1 x2 x3 x4 x5 x6 `;  
`delta4_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
     delta_x4 x1 x2 x3 x4 x5 x6 `;  
`delta4_y (sqrt x7) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x8) (sqrt x9) = 
   delta_x4 x7 x2 x3 x4 x8 x9 `; 
`dih2_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
     dih2_x x1 x2 x3 x4 x5 x6 `;  
`dih3_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
    dih3_x x1 x2 x3 x4 x5 x6 `;  
`dih_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
     dih_x x1 x2 x3 x4 x5 x6 `;  
`dih4_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
    dih4_x x1 x2 x3 x4 x5 x6 `;  
`dih5_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
   dih5_x x1 x2 x3 x4 x5 x6 `;  
`dih6_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
   dih6_x x1 x2 x3 x4 x5 x6 `;  
`delta_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
   delta_x x1 x2 x3 x4 x5 x6 `; 
`delta_y (sqrt x7) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x8) (sqrt x9) =
   delta_x x7 x2 x3 x4 x8 x9`;
`vol_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
   vol_x x1 x2 x3 x4 x5 x6 `;  
`eta_y (sqrt x1) (sqrt x2) (sqrt x6) pow 2 = eta2_126 x1 x2 x3 x4 x5 x6 `;  
`eta_y (sqrt x1) (sqrt x3) (sqrt x5) pow 2 = eta2_135 x1 x2 x3 x4 x5 x6 `;  
`eta_y (sqrt x4) (sqrt x5) (sqrt x6) pow 2 = eta2_456 x1 x2 x3 x4 x5 x6 `;  
`vol3f (sqrt x1) (sqrt x2) (sqrt x6)  sqrt2 lfun = vol3f_x_lfun x1 x2 x3 x4 x5 x6 `;  
`vol_y sqrt2 sqrt2 sqrt2 (sqrt x1) (sqrt x2) (sqrt x6)  = vol3_x_sqrt x1 x2 x3 x4 x5 x6 `;  
`vol3f_sqrt2_lmplus (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5)      (sqrt x6) =
   vol3f_x_sqrt2_lmplus x1 x2 x3 x4 x5 x6`;
 ];;


  let rec get_consts tm = 
    match tm with       
        Var(s,t) -> []
      | Const(s,t) -> [tm ]
      | Comb (t1,t2) -> get_consts t1 @ get_consts t2 
      | Abs (t1,t2) -> get_consts t1 @ get_consts t2 ;;

  let get_goal_const gl = 
    let w = goal_concl gl in
      get_consts w;;

(*
  let X_SQRT_COMPOUND_ORDER_TAC:tactic =  
    fun gl ->
      let (asl,a) = dest_goal gl in
      let fr = frees a in
	(* edited next line Nov 2012, svn 2931 *)
      let vos = filter (fun t -> subset (frees (fst (dest_eq t))) fr ) X_SQRT_COMPOUND_ORDER  in
    let X_SQRT_COMPOUND_ORDER_FIL = list_mk_conj vos in
  (SUBGOAL_THEN  X_SQRT_COMPOUND_ORDER_FIL  (fun t-> REWRITE_TAC[t;(GSYM Sphere.rhazim_x); (GSYM Sphere.rhazim2_x);xspec (GSYM Sphere.rhazim3_x)]) THENL[ (GEN_MESON_TAC 0 200 1[Sphere.norm2hh_x;rad2_x_y;delta_x4_delta4_y;dih_x_y;dih2_x_y;delta_x_y;dih3_x_y;tauq_x_y;GSYM Sphere.dih4_x;GSYM Sphere.dih5_x;GSYM Sphere.dih6_x;vol_x_y;Sphere.eta2_126; Sphere.eta2_135;Sphere.eta2_456;GSYM vol3f_x_lfun;GSYM vol3_x_sqrt;GSYM vol3f_x_sqrt2_lmplus]);ALL_TAC]) gl
;;
*)

  let X_SQRT_COMPOUND_ORDER_TAC:tactic =  
    fun gl ->
      let (asl,a) = dest_goal gl in
      let fr = frees a in
      let c = get_goal_const gl in
      let vos = filter (fun t -> subset (frees (fst (dest_eq t))) fr ) X_SQRT_COMPOUND_ORDER  in
      let vos' = filter (fun t -> mem (fst (strip_comb (fst(dest_eq t)))) c) vos in
    let vconj = list_mk_conj vos' in
  (SUBGOAL_THEN  vconj  (fun t-> REWRITE_TAC[t;(GSYM Sphere.rhazim_x); (GSYM Sphere.rhazim2_x);xspec (GSYM Sphere.rhazim3_x)]) THENL[ (GEN_MESON_TAC 0 200 1[Sphere.norm2hh_x;rad2_x_y;delta_x4_delta4_y;dih_x_y;dih2_x_y;delta_x_y;dih3_x_y;tauq_x_y;GSYM Sphere.dih4_x;GSYM Sphere.dih5_x;GSYM Sphere.dih6_x;vol_x_y;Sphere.eta2_126; Sphere.eta2_135;Sphere.eta2_456;GSYM vol3f_x_lfun;GSYM vol3_x_sqrt;GSYM vol3f_x_sqrt2_lmplus]);ALL_TAC]) gl
;;


let X_FROZEN_COMPOUND =(map SPEC_ALL [
 vol3r_126_x;REAL_ARITH `&2 = #2.0`;
 GSYM arclength_x_234;
 GSYM arclength_x_126;
 GSYM Functional_equation.nonf_gamma2_x1_div_a_v2;
 ] ) ;;


let decimal = REWRITE_RULE [REAL_ARITH `&2 = #2.0`];;

let X_COMPOUND_DEF =  [REAL_ARITH `&2 = #2.0`;
 beta_bump_force_x;
 GSYM Sphere.gchi1_x;GSYM Sphere.gchi2_x;
 GSYM Sphere.gchi3_x;GSYM Sphere.gchi4_x;
 GSYM Sphere.gchi5_x;GSYM Sphere.gchi6_x;
 GSYM Sphere.arclength_x1;GSYM Sphere.arclength_x2;
 GSYM Enclosed.quad_cross_diag2_x;
 tauq_x_y;] @ 
  map decimal [
 GSYM Sphere.ldih_x;GSYM Sphere.ldih2_x;
 GSYM Sphere.ldih3_x;GSYM Sphere.ldih6_x;
  ];;

let X_OF_Y_TAC = 
  (
  (DISCH_TAC) THEN
  TRY (MATCH_MP_TAC ineq_square2) THEN TRY (MATCH_MP_TAC ineq_square2_9) THEN
  (REWRITE_TAC basic_constants_nn) THEN
  REWRITE_TAC[GSYM CONJ_ASSOC] THEN
  (REPEAT (CONJ_TAC THENL[MP_TAC hminus_gt THEN MP_TAC sqrt3_nn THEN REWRITE_TAC[Sphere.h0;Sphere.hplus] THEN REAL_ARITH_TAC;ALL_TAC])) THEN
  (REPEAT GEN_TAC) THEN
  (REWRITE_TAC [sol_y_123;taum_123]) THEN
  (REWRITE_TAC[ineq]) THEN
  (SQRT_SQRT_TAC) THEN
    (REPEAT DISCH_TAC) THEN
    (* remove x variable compound *)
    X_SQRT_COMPOUND_ORDER_TAC THEN
    FROZEN_REWRITE_TAC X_FROZEN_COMPOUND THEN 
    REWRITE_TAC X_COMPOUND_DEF THEN
    (* repackage *)
    (REPEAT (FIRST_X_ASSUM MP_TAC)) THEN
    (REWRITE_TAC[GSYM ineq_expand9]) THEN
    (REWRITE_TAC[GSYM ineq_expand6]) THEN
    (REWRITE_TAC[REAL_ARITH `(&2 = #2.0) /\ (x pow 2 = x * x) /\ (#2.0 * #2.0 = #4.0) /\ (#2.18 * #2.18 = #4.7524 ) /\ (#2.52 * #2.52 = #6.3504)`;sqrt8_2;sqrt2_sqrt2;Sphere.h0;Sphere.hplus]) THEN
  ALL_TAC
  );;

let EXPAND_lfun = 
  (SUBGOAL_THEN `lfun x1 = (#1.26 - proj_x1  (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real))/(#0.26)` (fun t-> REWRITE_TAC[t])) THENL [
  REWRITE_TAC[Sphere.lfun;proj_x1;Sphere.h0;REAL_ARITH `#1.26 - &1 = #0.26`];ALL_TAC] ;;

let REMOVE_dummy = SUBGOAL_THEN `!(f:bool). (!(dummy:real). ineq [&1,dummy,&1] f) = (!(x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real). ineq[(&1,x1,&1);(&1,x2,&1);(&1,x3,&1);(&1,x4,&1);(&1,x5,&1);(&1,x6,&1)] f)` (fun t-> REWRITE_TAC[t]) THENL[  REWRITE_TAC[ineq] THEN  MESON_TAC[REAL_ARITH `~(&2 <= &1)`]; ALL_TAC];;

let EXPAND_1var = SUBGOAL_THEN `!(f:real->bool) a b. (!(y:real). ineq [a,y,b] (f y)) = (!(x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real). ineq[(a,x1,b);(&1,x2,&1);(&1,x3,&1);(&1,x4,&1);(&1,x5,&1);(&1,x6,&1)] (f x1))` (fun t-> REWRITE_TAC[t]) THENL [REWRITE_TAC[ineq] THEN  MESON_TAC[REAL_ARITH `(&1 <= &1)`];ALL_TAC] ;;


let REAL_SIMPLIFY_EXPRESSION =  let
  arith = REAL_ARITH `!x y z. (&8 = #8) /\ (x - y = x + (-- #1.0) * y) /\ 
  (x * (y+z) = x * y + x * z) /\ (y+z) * x = y * x + z * x /\ 
  (x + y) + z = x + y + z /\ (-- x * -- y = x * y) /\ (x * -- y = -- x * y) /\ 
  (-- x * y = -- (x * y)) /\ (x * y) * z = x * y * z /\ -- #1.0 * x = -- x /\ 
  -- (x + y) = -- x + (--y) /\ -- (-- x) = x /\ (-- (-- x * y) = x * y) /\
   #0.0 = &0 /\ #0 = &0 /\ &0 * x = &0 /\ x * &0 = &0 /\ (&0 + x = x) /\ 
  -- &0 = &0 /\ (x + &0 = x) /\ (&0 + x = x) /\
  &1 * x = x /\ -- &1 * x = -- x /\ x * sqrt8 = sqrt8 * x    ` in
  (REWRITE_TAC[REAL_POW_MUL;real_div;REAL_MUL_LZERO; REAL_MUL_RZERO;arith]) ;;

let EXPAND_SQRT = 
  (SUBGOAL_THEN `sqrt x1 = sqrt_x1 x1 x2 x3 x4 x5 x6 /\ sqrt x2 = sqrt_x2 x1 x2 x3 x4 x5 x6 /\ sqrt x3 = sqrt_x3 x1 x2 x3 x4 x5 x6 /\ sqrt x4 = sqrt_x4 x1 x2 x3 x4 x5 x6 /\ sqrt x5 = sqrt_x5 x1 x2 x3 x4 x5 x6 /\ sqrt x6  = sqrt_x6 x1 x2 x3 x4 x5 x6 /\ sqrt x8 = sqrt_x5 (x7:real) x2 x3 x4 x8 x9 /\ sqrt x9 = sqrt_x6 x7 x2 x3 x4 x8 x9 ` (fun t->REWRITE_TAC[t])  THENL [REWRITE_TAC[sqrt_x1;sqrt_x2; sqrt_x3;sqrt_x4;sqrt_x5;sqrt_x6];ALL_TAC]) ;;

(* for 1d inequality involving vol2f marchal *)

let EXPAND_vol2 = 
REWRITE_TAC[vol2f_marchal_pow_y;vol2r_y] THEN
  SUBGOAL_THEN `x1 pow 1 = promote pow1 x1 x2 x3 x4 x5 x6 /\ 
  x1 pow 2 = promote_pow2 x1 x2 x3 x4 x5 x6 /\ 
  x1 pow 3 = promote_pow3 (x1:real) (x2:real) (x3:real) 
    (x4:real) (x5:real) (x6:real) /\ 
  x1 pow 4 = promote pow4 x1 x2 x3 x4 x5 x6` 
  (fun t->REWRITE_TAC[t]) 
  THENL[ REWRITE_TAC[promote;pow1;pow2;pow3;pow4;promote_pow2;promote_pow3];
	 REWRITE_TAC[LET_DEF;LET_END_DEF]] ;;

let DEF_expand = 
 [Sphere.a_spine5;Sphere.b_spine5;Sphere.mm1;
  Sphere.flat_term;Sphere.beta_bump_lb;REAL_POW_2;
  Sphere.h0;
  Sphere.mm2;GSYM Sphere.sqrt2;GSYM Sphere.sqrt3;
  GSYM Sphere.sqrt8;sol0_const1;sqrt2_sqrt8;
  Sphere.mm1;Sphere.mm2;Sphere.tau0;Sphere.hplus;
  tame_table_d_values;Sphere.vol2f;
  Sphere.lfun; (* added Oct 17, 2010 *)
 ];;

(* *************************************************************************** *)
(* SERIES3Q1H *)
(* *************************************************************************** *)

let SERIES3Q1H_5D_TAC = 
  let instjx = INST_TYPE [(`:real`,`:A`);(`:real`,`:B`);(`:real`,`:C`);
			  (`:real`,`:D`);(`:real`,`:E`);(`:real`,`:F`)] in
  let PROJX = map instjx [ proj_x1;proj_x2;proj_x3;proj_x4;proj_x5;proj_x6] in
  let projx = list_mk_conj (map (concl o GSYM) PROJX ) in
 (REPEAT STRIP_TAC   THEN MATCH_MP_TAC ineq6_of_ineq5) THEN
 (REWRITE_TAC[Sphere.pathL;Sphere.pathR;Sphere.hplus;Sphere.h0]) THEN
 (REWRITE_TAC[ineq_expand6]) THEN
 (REPEAT GEN_TAC THEN REPEAT DISCH_TAC) THEN
 (SUBGOAL_THEN projx (fun t -> PURE_ONCE_REWRITE_TAC[t])) THENL
 [(REWRITE_TAC PROJX);ALL_TAC] THEN
 (REPEAT (POP_ASSUM MP_TAC)) THEN
 (REWRITE_TAC[GSYM ineq_expand6]);;

(* *************************************************************************** *)
(* STYLIZE and WRAPUP *)
(* *************************************************************************** *)

let STYLIZE_TAC = 
  let real_arith_lt = 
    REAL_ARITH `(x < y <=> (x- y < &0)) /\ (x <= y <=> (x-y <= &0))` in
  let real_arith_gt = 
    REAL_ARITH `(x > y <=> (y - x < &0) ) /\ (x >= y <=> (y-x <= &0))` in
  let real_arith_unit =  
    REAL_ARITH `unit0 * (x + y ) = unit0 * x + unit0 * y /\ 
    unit0 * --x = --(unit0 * x) /\ (unit0 * x  = x  * unit0) /\ 
    (x * y) * z = x * y * z` in
  let subgoal_proj = `x2 * unit0 = 
    proj_x2 (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) /\ 
      x3 * unit0 = proj_x3 (x1) (x2) (x3) (x4) (x5) (x6) /\ 
    x1 * unit0 = proj_x1 (x1) (x2) (x3) (x4) (x5) (x6) /\
    (x4 * unit0 = proj_x4 x1 x2 x3 x4 x5 x6) /\
    (x5 * unit0 = proj_x5 x1 x2 x3 x4 x5 x6) /\ 
    (x6 * unit0 = proj_x6 x1 x2 x3 x4 x5 x6)` in
  let subgoal_unit6 = `unit0 = unit6
    (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real)` in
  REMOVE_dummy  THEN
  EXPAND_1var THEN
  DISCH_TAC THEN REPEAT GEN_TAC THEN
  (* rename added Nov 27, 2010 *)
  X_RENAME_TAC THEN 
  REPEAT GEN_TAC THEN
  EXPAND_vol2 THEN
  (* prev line added Sep 8 , 2010 *)
  REWRITE_TAC  DEF_expand THEN
  (* prev line moved down sep 8 *)
  EXPAND_lfun THEN
  REWRITE_TAC[ineq] THEN
  (REPEAT DISCH_TAC) THEN
  (ONCE_REWRITE_TAC[real_arith_lt]) THEN
  (REWRITE_TAC[real_arith_gt]) THEN
  EXPAND_SQRT THEN
  REAL_SIMPLIFY_EXPRESSION THEN
  SUBGOAL_THEN `!a. x1:real * a = a * x1` (fun t->REWRITE_TAC[t])
  THENL[REAL_ARITH_TAC;ALL_TAC] THEN (* added OCt 17, 2010 *)
  (SUBGOAL_THEN `!x. ((x < &0) <=> (unit0 * x < &0)) /\ 
     ((x <= &0) <=> (unit0 * x <= &0))` (fun t -> ONCE_REWRITE_TAC[t])) 
  THENL [REWRITE_TAC[unit0;REAL_ARITH `&1 * x = x`];ALL_TAC] THEN
  (REWRITE_TAC[REAL_ARITH `f x1 x2 x3 x4 x5 x6 * (y:real) = 
       y * f x1 x2 x3 x4 x5 x6 /\ ((x * y) * z = x * y * z)`]) THEN
  (REWRITE_TAC[real_arith_unit]) THEN
  (REWRITE_TAC[REAL_ARITH `unit0 * x = x * unit0`]) THEN
  (REWRITE_TAC[unit0f]) THEN
    (SUBGOAL_THEN subgoal_proj  (fun t-> REWRITE_TAC[t])
  THENL [REWRITE_TAC[proj_x3;proj_x4;proj_x5;proj_x6;proj_x2;
		    proj_x1;unit0;REAL_ARITH `x * &1 = x`];ALL_TAC]) THEN
  (SUBGOAL_THEN subgoal_unit6  (fun t-> REWRITE_TAC[t]) 
     THENL [REWRITE_TAC[unit0;unit6];ALL_TAC]) THEN
  (REPEAT (FIRST_X_ASSUM MP_TAC)) THEN
  (REWRITE_TAC[GSYM ineq_expand9]) THEN
  (REWRITE_TAC[GSYM ineq_expand6]) THEN DISCH_TAC;;

let WRAPUP_TAC = 
  (REWRITE_TAC[REAL_ARITH `(x * y * (z:real)) = (x * y) * z`]) THEN
  (REWRITE_TAC[REAL_ARITH `(y:real) * f x1 x2 x3 x4 x5 x6 =  
       (f x1 x2 x3 x4 x5 x6) *  y `]) THEN
  (REWRITE_TAC[REAL_ARITH ` -- (x * y) = x * (-- y) `]) THEN
  (REWRITE_TAC[REAL_ARITH ` -- (f x1 x2 x3 x4 x5 x6) = 
       f x1 x2 x3 x4 x5 x6 * -- &1`]) THEN
  (REWRITE_TAC[REAL_ARITH `(x * y) * (z:real) = x * y * z`]) THEN
  (REWRITE_TAC[REAL_ARITH `inv y = (&1/y)`]);;




(* ========================================================================== *)
(*    VERIFYING INEQUALITIES                                                  *)
(* ========================================================================== *)

let idq_fields idq = (idq.idv,idq.tags, idq.ineq);;

(*
let preprocess (s,tags,case) = 
  let is_xconvert = mem Xconvert tags in
  let is_branch = mem Branching tags in
  let strip_let_case = strip_let_tm case in
  let _ = report ("process and exec: "^s) in
  let _ = g (strip_let_case) in
  let _ = e(PRELIM_REWRITE_TAC) in 
  let NONLIN_INTRO = e(MP_TAC (REWRITE_RULE[] NONLIN)) in
  let _ = if (is_branch) then e(BRANCH_TAC) else e(ALL_TAC) in
  let _ = if (is_xconvert) then e (X_OF_Y_TAC) else e(ALL_TAC) in
  let _ =   if (is_branch && not(is_xconvert)) then
    e(SERIES3Q1H_5D_TAC) else e(ALL_TAC) in
  let _ = e (STYLIZE_TAC) in
  let _ = e (WRAPUP_TAC) in
  let testineq = snd(top_goal()) in
    (s,tags,testineq);;
*)

(* rewritten Jan 28, 2013 *)

let preprocess (s,tags,case) = 
  let is_xconvert = mem Xconvert tags in
  let is_branch = mem Branching tags in
  let _ = report ("process and exec: "^s) in
  let LET_ELIM_TAC = CONV_TAC (REDEPTH_CONV let_CONV) in
  let tacl = 
    [LET_ELIM_TAC;
     PRELIM_REWRITE_TAC;
     MP_TAC (REWRITE_RULE[] NONLIN);
     if (is_branch) then BRANCH_TAC else ALL_TAC;
     if (is_xconvert) then X_OF_Y_TAC else ALL_TAC;
     if (is_branch && not(is_xconvert)) then SERIES3Q1H_5D_TAC else ALL_TAC;
     STYLIZE_TAC;
     WRAPUP_TAC] in
  let _ = g (case) in
  let _ = e (EVERY tacl) in
  let testineq = snd(top_goal()) in
    (s,tags,testineq);;

let process_and_exec ex (s,tags,case) = 
  let _ = report ("process and exec: "^s) in
  let (s,tags,testineq) = preprocess (s,tags,case) in
    execute_interval ex tags s testineq;;

(*
let process_and_exec ex (s,tags,case) = 
  let _ = report ("process and exec: "^s) in
  let is_xconvert = mem Xconvert tags in
  let is_branch = mem Branching tags in
  let strip_let_case = strip_let_tm case in
  let _ = report ("process and exec: "^s) in
  let _ = g (strip_let_case) in
  let _ = e(PRELIM_REWRITE_TAC) in 
  let NONLIN_INTRO = e(MP_TAC (REWRITE_RULE[] NONLIN)) in
  let _ = if (is_branch) then e(BRANCH_TAC) else e(ALL_TAC) in
  let _ = if (is_xconvert) then e (X_OF_Y_TAC) else e(ALL_TAC) in
  let _ =   if (is_branch && not(is_xconvert)) then
    e(SERIES3Q1H_5D_TAC) else e(ALL_TAC) in
  let _ = e (STYLIZE_TAC) in
  let _ = e (WRAPUP_TAC) in
  let testineq = snd(top_goal()) in
    execute_interval ex tags s testineq;;
*)


(*

let get_firstexact s = idq_fields (hd(Ineq.getexact s));;

let process_and_exec_by_id ex s = 
  let testcase = get_firstexact  s in 
    process_and_exec ex s testcase;;

testsplit_idq
  ** strips let
  ** splits cases at h0 according to the split tags
  ** ships the cases off to process_and_exec.
*)

(*
let testsplit_idq ex idq = 
  let (s,tags,ineq) = idq_fields idq in
  let ineq_strip_let = strip_let_tm ineq in
  let suffix i n = Printf.sprintf "%s split(%d/%d)" s i n in 
  let ls = get_split_tags idq in
    if (ls = []) then process_and_exec ex (s,tags,ineq_strip_let) else
      let cases = split_all_h0  [(ineq_strip_let, ls)] in
      let n = length cases in
	for i=0 to (n - 1) do 
	  try (process_and_exec ex (suffix i n,tags,(List.nth cases i)))
	  with Failure s -> failwith (s ^ " case fail: " ^(string_of_int i))
	done;;
*)

let preprocess_split_idq idq = 
  let (s,tags,ineq) = idq_fields idq in
  let ineq_strip_let = strip_let_tm ineq in
  let suffix i n = Printf.sprintf "%s split(%d/%d)" s i n in 
  let ls = get_split_tags idq in
    if (ls = []) then [preprocess (s,tags,ineq_strip_let)] else
      let cases = split_all_h0  [(ineq_strip_let, ls)] in
      let n = length cases in
	map (fun i ->
	    try preprocess (suffix i n,tags,List.nth cases i)
	    with Failure s -> failwith (s^ " case fail: " ^ (string_of_int i))) 
	  (0--(n-1));;

let testsplit_idq ex idq = 
  let splits = preprocess_split_idq idq in
    map (fun (s,tags,testineq) -> execute_interval ex tags s testineq) splits;;

let testsplit ex s = testsplit_idq ex (hd (Ineq.getexact s));;


end;;