(* ========================================================================== *)
(* FLYSPECK - CODE FORMALIZATION                                              *)
(*                                                                            *)
(* Chapter: Linear Programming Inequalities                                          *)
(* Author: Thomas C. Hales                                                    *)
(* Date: 2010-08-01                                                           *)
(* ========================================================================== *)


(* Generate new linear inequalities that are needed for linear
   programs.  This file is part of the scaffolding of the project.  It
   was used extensively to generate new inequalities that have been used
   in the linear programming part of the proof.  These inequalities have
   all been added to ineq.hl.

   This file has served its purpose and is no longer needed.

   Mathematica and Cfsqp both need to be installed before this code can
   be executed.  The mathematica runs are generated with the script found
   in the file whose pathname is given by the string math_initf.

  There is also tight coordinate between this file and the linear programming
   directory glpk/tame_archive/  

   WARNING: For the functions that write external files and reread them into a
   reference, it seems that it is necessary to open this module for it to
   function properly.
*)


(*

  There are several steps required to create an inequality and rerun
  linear programs.  This file automates these steps.

  0/- Analyze output to find a triangle that is bad with
  sorted_azim_weighted_diff darts_of_std_tri amx;; 
  1/- Mathematica.  holineq[Dihedral,{2,2,2,2,2,2},{2.52,2.25,2.52,2.52,2.52,2.52}].
  2/- save, strip "\ " at the end of lines.
  3/- run Parse_ineq.execute_cfsqp new_ineq to check if it seems to be true.
  If not, edit until it is correct.
  4/- change the domain info to something that is recognized by glpk model.
  5/- add new_ineq to put in archive.
  6/- Parse_ineq.lp_string()  to generate model body.mod.
  7/- Save body.mod in directory: glpk/tame_archive.
  8/- Lpproc.make_model() to make the new model.
  9/- Rerun the linear programs.

*)

(* Interaction with inequality generation. Mathematica is needed for
    the next bit.  *)




flyspeck_needs "nonlinear/ineq.hl";;

module Temp_ineq = struct


  let all_forall = Sphere.all_forall;;

let junkval =  { 
 idv = "junk";
 doc="junk";
 tags=[];
 ineq = 
 all_forall `ineq [(&1,y,&2)] (y > &0)`;
 };;

(* ========================================================================== *)
(* Mathematica section. *)
(* ========================================================================== *)

let mathfile = "/Applications/Mathematica\ 5.2.app/Contents/MacOS/MathKernel";;

let math_initf= "/Users/thomashales/Desktop/googlecode/flyspeck/mathematica/auto_mkineq.m";;

let mathtemp = "/tmp/mathtemp4.hl";; (* hardwired into mathematica code *)

let command_string fn lb ub = 
(mathfile^" -run \"<< "^math_initf^ ";\n runQuit["^fn^","^lb^","^ub^"];\"");;

let command_string_p fn p lb ub = 
(mathfile^" -run \"<< "^math_initf^ ";\n runQuitp["^fn^","^p^","^lb^","^ub^"];\"");;

(* ========================================================================== *)
(* cfsqp and ineq creation section.  The mathematica code produces a
   guess of an inequality (that is often inaccurate).  It is run to cfsqp
   to correct the guess, then packaged into an ineq.  *)
(* ========================================================================== *)


 let cfsqp_margin idq = 
  let cfsqp_dir = flyspeck_dir^"/../cfsqp" in
  let _ =  Parse_ineq.mk_cc (cfsqp_dir ^ "/tmp/t.cc") idq in
  let _ = Parse_ineq.compile() in 
  let _ = (0=  Sys.command(cfsqp_dir^"/tmp/t.o > /tmp/cfsqp_out.txt")) or failwith "execution error" in
  let s2 =  process_to_string("grep 'constrained min:' /tmp/cfsqp_out.txt | sed 's/^constrained min://g' | tr -d '\n'") in  float_of_string s2;;

let padded_decimal_of_cfsqp idq = (* add padding to make nonlinear verifications easier *)
   let r =  cfsqp_margin idq in
     let a = if r < 0.001 then 0.001 -. r else 0.0 in
     let b = int_of_float (Pervasives.ceil (10000.0 *. a)) in
     mk_binop `DECIMAL` (mk_numeral (Int b)) (mk_numeral (Int 10000));;

let add_decimal ineq a =
  let (quant,raw) = strip_forall ineq in
  let (dom,gt) = dest_binop `ineq` raw in
  let (lhs,rhs) = dest_binop `real_gt` gt in
  let lhs' = mk_binop `real_add` lhs a in
  let gt' = mk_binop `real_gt` lhs' rhs in
  let raw' = mk_binop `ineq` dom gt' in
  let ineq' = list_mk_forall (quant,raw') in
   ineq';;

let sqrt8_sqrt2 = prove_by_refinement(
  `sqrt (&8) = &2 * sqrt2`,
  (* {{{ proof *)
  [
  SIMP_TAC[Sphere.sqrt2;SQRT_MUL;
     REAL_ARITH `&8 = &4 * &2 /\ &0 <= &2 /\ &0 <= &4`;];
  REWRITE_TAC[REAL_ARITH `&4 = &2 pow 2 /\ abs(&2) = &2`;POW_2_SQRT_ABS];
  ]);;
  (* }}} *)

let constant_fix = [sqrt8_sqrt2;Sphere.h0;Sphere.sqrt2;Sphere.sqrt8;
  REAL_ARITH `#2 = &2 /\ #2.0 = &2 /\ &2 * #1.26 = #2.52 /\ #2.00 = &2 `];;

(*
problem: dart_std3 and dart_std3_big have the same "domain".
*)

let reduced_dart_classes() = search_thml (term_match [])  [omit `dart_std3_big`] (map (fun t-> ("",t)) (!Ineq.dart_classes));;

let dart_class_rew() = map (REWRITE_RULE constant_fix)  (map (GSYM o snd ) (reduced_dart_classes()));;

let convert_domain ineq = snd (dest_eq (concl 
  (REWRITE_CONV ( constant_fix @ (dart_class_rew()) ) ineq)));;

let adjust = 
  let zero = `&0` in
  fun s idq ->
    let a = padded_decimal_of_cfsqp idq in
      if (a = zero) then idq else
     {
     idv = idq.idv;
     doc = idq.doc ^ s^ "\n  The inequality has been fitted to cfsqp margins.";
     tags = idq.tags;
     ineq = convert_domain(add_decimal idq.ineq a);
     };;

(* By loading the file mathtemp, the value of the reference tempval is
    changed to the output of the mathematica run *)

  let tempval = ref junkval;;

let generate_ineq_datum = 
  fun  fname slb sub ->
  let example = command_string fname slb sub in 
  let _ = Sys.command(example) in
  let _ = loadt mathtemp in
    adjust 
  ("\n Generated by generate_ineq_datum with input "^
   fname ^" "^slb^" "^sub^".") (!tempval);;

let generate_ineq_datum_p =
  fun fname sp slb sub ->
  let example = command_string_p fname sp slb sub in 
  let _ = Sys.command(example) in
  let _ = loadt mathtemp in
    adjust  ("\n Generated by generate_ineq_datum_p with input "^
   fname ^" "^sp^" "^slb^" "^sub^".") (!tempval);;

(* ========================================================================== *)
(* Ocaml function section. *)
(* ========================================================================== *)


(*

This is a quick hack that converts an Hol-light inequality (>) into a function
that gives the difference between the right and left hand sides of an
ineq  LHS > RHS for given values y1..y6.

It generates the code as a string, writes it to an external file, then
reads it back into the reference tempfn.

It assumes that the inequality has this particular form (six free variables
named y1..y6.
*)

let ocaml_funstring_of_ineq iqd = 
  let t = snd(strip_forall (Parse_ineq.prep_term (iqd.ineq))) in
  let (_,tm) = dest_comb t in
let (lhs,rhs) = dest_binop `real_gt` tm in
let olhs = Parse_ineq.ocaml_string_of_term lhs in
let orhs = Parse_ineq.ocaml_string_of_term rhs in
  "(fun y1 y2 y3 y4 y5 y6 -> " ^ olhs ^ " -. " ^ orhs ^ ");;";;

let tempfn = ref (fun (y1:float) (y2:float) (y3:float) (y4:float) (y5:float) (y6:float) -> 0.0);;

let ocaml_eval =
  let tempfile = Filename.temp_file "ocaml_eval_" ".ml" in
  fun iqd ->
  let _ = Parse_ineq.output_filestring tempfile 
    ("tempfn := " ^ ocaml_funstring_of_ineq iqd) in
  let _ = loadt tempfile in
    !tempfn;;

end;;