text_formalization/nonlinear/parse_ineq.hl

   1 (* ========================================================================== *)
   2 (* FLYSPECK - BOOK FORMALIZATION                                              *)
   3 (*                                                                            *)
   4 (* Chapter: nonlinear inequalities                                            *)
   5 (* Author:  Thomas Hales starting from Roland Zumkeller's ccform function     *)
   6 (* Date: 2010-05-09                                                           *)
   7 (* ========================================================================== *)
   8
   9 (*
  10
  11    Use HOL Light specifications of functions and inequalities to
  12    automatically generate code in other languages.
  13
  14    ** cfsqp C++ code,
  15    ** glpk for linear programming
  16    ** Objective Caml export of functions
  17    ** LaTeX documentation of inequalities.
  18
  19    C++ code generation for interval arithmetic is done
  20    separately in optimize.hl.
  21    (This file is independent from that one.)
  22
  23    CFSQP NOTES
  24    Executing the cfsqp produces three output files:
  25    /tmp/cfsqp_err.txt,
  26    tmp/t.cc tmp/t.o (in the cfsqp directory)
  27    Use macro_add to add new function (that is expanded with rewrites)
  28    Use autogen_add to add new function (automatic C code generation)
  29 *)
  30
  31 (*
  32 flyspeck_needs "general/flyspeck_lib.hl";;
  33 flyspeck_needs "general/sphere.hl";;
  34 flyspeck_needs "nonlinear/ineq.hl";;
  35 flyspeck_needs "nonlinear/lemma.hl";;
  36 *)
  37
  38 module Parse_ineq = struct
  39
  40   open Sphere;;
  41   open Nonlinear_lemma;;
  42   open Functional_equation;;
  43
  44 let cfsqp_trialcount = ref 500;;
  45
  46 let dest_decimal = Flyspeck_lib.dest_decimal;;
  47
  48 let string_of_num' = Flyspeck_lib.string_of_num';;
  49
  50 let join_comma  = Flyspeck_lib.join_comma;;
  51
  52 let join_lines  = Flyspeck_lib.join_lines;;
  53
  54 let join_space  = Flyspeck_lib.join_space;;
  55
  56 let paren s =  "("^ s ^")";;
  57
  58 let nub = Flyspeck_lib.nub;;
  59
  60 (* moved to lemmas.hl.
  61 let strip_let_tm t = snd(dest_eq(concl(REDEPTH_CONV let_CONV t)));;
  62
  63 let strip_let t = REWRITE_RULE[REDEPTH_CONV let_CONV (concl t )] t;;
  64 *)
  65
  66 (* from HOL Light lib.ml *)
  67
  68 let rec (--) = fun m n -> if m > n then [] else m::((m + 1) -- n);;
  69
  70
  71 (* renamed from output_string to avoid Pervasives module name clash *)
  72
  73 let output_filestring = Flyspeck_lib.output_filestring;;
  74
  75
  76 (* ========================================
  77 cfsqp  generation
  78 ======================================== *)
  79
  80 let c_string_of_term t =
  81   let rec soh t =
  82     if is_var t then fst (dest_var t)
  83     else
  84       let (f,xs) = strip_comb t in
  85       let _ =  (is_const f) or
  86         failwith ("Non-const function name: " ^ string_of_term f) in
  87       let ifix i  =
  88         let [a;b] = xs in paren(soh a ^ " " ^ i ^ " " ^ soh b) in
  89       let ifix_swapped i  =
  90         let [b;a] = xs in paren(soh a ^ " " ^ i ^ " " ^ soh b) in
  91         (* soh: *)
  92         match fst (dest_const f) with
  93           | "real_gt" | "real_ge" | "real_sub" -> ifix "-"
  94           | "real_lt" | "real_le" -> ifix_swapped "-"
  95           | "real_add" -> ifix "+"
  96           | "real_mul" -> ifix "*"
  97           | "real_div" -> ifix "/"
  98           | "\\/" -> ifix "\\/"
  99           | "real_neg" -> let [a] = xs in paren("-" ^ soh a)
 100           | "acs" -> let [a] = xs in paren ("acos" ^ (paren (soh a)))
 101           | "real_of_num" -> let [a] = xs in soh a
 102           | "NUMERAL" -> let [_] = xs in string_of_num' (dest_numeral t)
 103           | "<" -> let [a;b] = xs in paren(soh a ^ " < " ^ soh b)
 104           | ">" -> let [a;b] = xs in paren(soh a ^ " > " ^ soh b)
 105           | "+" -> let [a;b] = xs in paren(soh a ^ " + " ^ soh b)
 106           | "*" -> let [a;b] = xs in paren(soh a ^ " * " ^ soh b)
 107           | "DECIMAL" ->  string_of_num' (dest_decimal t)
 108           | "COND" ->
 109               let [a;b;c] = xs in
 110                 paren(soh a ^ " ? " ^ soh b ^ " : " ^ soh c)
 111           | "atn2"      ->
 112               let [ab] = xs in
 113               let (a,b) = dest_pair ab in
 114                 paren("atn2( " ^ soh a ^ "," ^ soh b ^ ")")
 115           | s -> paren(s ^  "(" ^ join_comma(map soh xs) ^ ")") in
 116    (* c_string_of_term: *)
 117     try
 118       soh t
 119     with Failure s -> failwith (s^" .......   "^string_of_term t);;
 120
 121 let constant_names i =
 122   let rec cn b =function
 123     | Const (s,_) -> s::b
 124     | Comb (s,t) -> (cn [] s) @ (cn[] t) @ b
 125     | Abs(x,t) -> (cn b t)
 126     | _ -> b in
 127   nub (sort (<) (cn [] i ));;
 128
 129
 130 (* Rewrite unusual terms to prepare for C++ conversion *)
 131
 132 (* function calls are dealt with three different ways:
 133       - native_c: use the native C code definition of the function.
 134       - autogen: automatically generate a C style function
 135       - macro_expand: use rewrites to eliminate the function.
 136 *)
 137
 138 (* Native is the default case.  There is no need to list them, except
 139    as documentation. *)
 140
 141 let native_c =  [
 142   ",";"BIT0";"BIT1";"CONS";"DECIMAL"; "NIL"; "NUMERAL"; "_0"; "acs";
 143   "ineq";  "pi"; "adodec"; "bdodec";"real_add"; "real_div";"real_pow";"cos";
 144   "real_ge"; "real_mul"; "real_of_num"; "real_sub"; "machine_eps";
 145   (* -- *)
 146   "delta_x";"sol_y";"dih_y";"rhazim";
 147   "lmfun";"lnazim";"hminus";
 148   "wtcount3_y";"wtcount6_y";"beta_bumpA_y";"matan";"sqp";"sqn";
 149   (* added May 2011 *)
 150   (* "taum_m_diff_quotient";"taum_m_diff_quotient2"; *)
 151   (* added 2012 *)
 152   "ell_uvx";
 153   "root";
 154   "h0cut";
 155    ];;
 156
 157 let autogen = ref[];;
 158
 159 let prep_autogen = (function b -> snd(strip_forall (concl (strip_let b))));;
 160
 161 let autogen_remove thm =
 162   let thm' = prep_autogen thm in
 163     autogen := filter (fun t -> not (t = thm')) !autogen;;
 164
 165 let autogen_add thm =
 166   let thm' = prep_autogen thm in
 167   let _ = autogen_remove thm in (* duplicates create C compile error *)
 168     autogen := !autogen @ [thm'];;
 169
 170 (* cfsqp autogeneration *)
 171
 172 autogen := map prep_autogen
 173   [sol0;tau0;hplus;mm1;mm2;vol_x;sqrt8;sqrt2;sqrt3;rho_x;
 174    rad2_x;ups_x;eta_x;eta_y;volR;
 175    norm2hh;arclength;regular_spherical_polygon_area;
 176    beta_bump_force_y;
 177    a_spine5;b_spine5;beta_bump_lb;marchal_quartic;vol2r;
 178    tame_table_d;delta_x4;dih_x_alt;sol_x;delta4_squared_x;x1_delta_x;
 179    quad_root_plus_curry;
 180    edge_flat_rewrite;const1;taum;flat_term;
 181    taum_y1;taum_y2;taum_y1_y2;
 182    arclength_y1;arclength_y2;arc_hhn;asn797k;asnFnhk;
 183    lfun_y1;arclength_x_123;
 184
 185    acs_sqrt_x1_d4;acs_sqrt_x2_d4;
 186    tauq;enclosed_rewrite;
 187    Nonlin_def.sol_euler_x_div_sqrtdelta;
 188    Nonlin_def.dih_x_div_sqrtdelta_posbranch;
 189    dih5_y;dih6_y;
 190    num1;Nonlin_def.dnum1;
 191    flat_term_x;
 192
 193    eulerA_x;
 194
 195 Functional_equation.nonf_gamma3_x;
 196 Functional_equation.nonf_gamma23_full8_x;
 197 Functional_equation.nonf_gamma23_keep135_x;
 198 Nonlin_def.gamma2_x_div_azim_v2;
 199 Functional_equation.nonf_gamma2_x1_div_a_v2;
 200 (* Added 2013-05-20 *)
 201 Nonlin_def.mu_y;Nonlin_def.delta_x1;Nonlin_def.taud;
 202 Functional_equation.nonf_ups_126;
 203    Functional_equation.nonf_gamma3f_x_div_sqrtdelta;
 204    ];;
 205
 206 (* cfsqp definition macros expansion *)
 207
 208 let macros = ref  [Sphere.gamma4f];;
 209
 210 let macro_remove thm =
 211   let _ = macros := filter (fun t -> not (t = thm)) !macros in
 212     ();;
 213
 214 let macro_add thm =
 215   let _ = macro_remove thm in
 216     macros := !macros @ [thm];;
 217
 218 let get_macro_expand() = (!macros);;
 219
 220 macros := (
 221   [gamma4f;gamma4fgcy;vol4f;y_of_x_e;vol_y_e;rad2_y_e;
 222    vol3f;vol3r;vol2f;delta4_y;delta4_squared_y;x1_delta_y;
 223    gamma3f;gamma23f;  (* gamma23f_126_w1;gamma23f_red; *)
 224    gamma23f_red_03;gamma23f_126_03;
 225    GSYM quadratic_root_plus_curry;REAL_MUL_LZERO;
 226    REAL_MUL_RZERO;FST;SND;pathL;pathR;node2_y;node3_y;
 227    rhazim2;rhazim3;rhazim_x;rhazim2_x;rhazim3_x;
 228    (* functional code *)
 229    rotate2;rotate3;rotate4;rotate5;rotate6;
 230    uni;constant6;promote3_to_6;promote1_to_6;
 231    dummy6;
 232    mul6;add6;sub6;scalar6;compose6;I_DEF;
 233
 234    functional_proj_x1;functional_proj_y1;
 235    functional_proj_x2;functional_proj_y2;
 236    functional_proj_x3;functional_proj_y3;
 237    functional_proj_x4;functional_proj_y4;
 238    functional_proj_x5;functional_proj_y5;
 239    functional_proj_x6;functional_proj_y6;
 240    Nonlin_def.sol_euler345_x_div_sqrtdelta;
 241    Nonlin_def.sol_euler156_x_div_sqrtdelta;
 242    Nonlin_def.sol_euler246_x_div_sqrtdelta;
 243
 244
 245    Nonlin_def.dih4_x_div_sqrtdelta_posbranch;
 246    Nonlin_def.ldih_x_div_sqrtdelta_posbranch;
 247    Nonlin_def.ldih2_x_div_sqrtdelta_posbranch;
 248    Nonlin_def.ldih2_x_div_sqrtdelta_posbranch;
 249    Nonlin_def.ldih3_x_div_sqrtdelta_posbranch;
 250    Nonlin_def.ldih5_x_div_sqrtdelta_posbranch;
 251    Nonlin_def.ldih6_x_div_sqrtdelta_posbranch;
 252    taum_x;
 253    edge_flat2_x;
 254
 255    delta_126_x;
 256    delta_234_x;
 257    delta_135_x;
 258
 259    (* may 2011 additions *)
 260    Nonlin_def.rhazim_x_div_sqrtdelta_posbranch;
 261    Nonlin_def.rhazim2_x_div_sqrtdelta_posbranch;
 262    Nonlin_def.rhazim3_x_div_sqrtdelta_posbranch;
 263    Nonlin_def.tau_residual_x;
 264    edge_flat_x;
 265    delta_sub1_x;
 266
 267
 268    (* 2012 *)
 269    (* may 2013 additions *)
 270    Nonlin_def.mud_135_x;
 271    Nonlin_def.mud_126_x;Nonlin_def.mud_234_x;
 272    Nonlin_def.mudLs_234_x;Nonlin_def.mudLs_135_x;Nonlin_def.mudLs_126_x;
 273    Nonlin_def.taud_x;Functional_equation.nonfunctional_mu6_x;
 274    Functional_equation.nonfunctional_taud_D1;
 275    Functional_equation.nonfunctional_taud_D2;
 276    Functional_equation.nonfunctional_edge2_126_x;
 277    Functional_equation.nonfunctional_edge2_135_x;
 278    Functional_equation.nonfunctional_edge2_234_x;
 279    Nonlin_def.flat_term2_126_x;
 280    Nonlin_def.flat_term2_135_x;Nonlin_def.flat_term2_234_x
 281    ] @ (!Ineq.dart_classes));;
 282
 283 let prep_term t =
 284   let t' = REWRITE_CONV (get_macro_expand()) (strip_let_tm t) in
 285   let (a,b)=  dest_eq (concl t') in
 286     b;;
 287
 288 let cc_function t =
 289   let args xs =
 290     let ls = map (fun s -> "double "^s) xs in join_comma ls in
 291   let (lhs,rhs) = dest_eq (prep_term t) in
 292   let (f,xs) = strip_comb lhs in
 293   let ss = map c_string_of_term xs in
 294   let p = Printf.sprintf in
 295   let s = join_lines [
 296      p"double %s(" (fst (dest_const f)); args ss;
 297      p") { \nreturn ( %s ); \n}\n\n"  (c_string_of_term rhs);
 298      ]
 299   in s;;
 300
 301 let dest_ineq ineq =
 302   let t = snd(strip_forall (prep_term (ineq))) in
 303   let (vs,i) = dest_comb t in
 304   let (_,vs) = dest_comb vs in
 305   let vs = dest_list vs in
 306   let vs = map (fun t -> let (a,b) = dest_pair t in (a,dest_pair b)) vs in
 307   let vs = map (fun (a,(b,c)) -> (a, b, c)) vs in
 308     (t,vs,disjuncts i);;
 309
 310 let c_dest_ineq ineq =
 311   let cs = c_string_of_term in
 312   let (b,vs,i) = dest_ineq ineq in
 313     (cs b, map (fun (a,b,c) -> (cs a, cs b,cs c)) vs,map cs i);;
 314
 315 (* generate c++ code of ineq *)
 316
 317 let case p i j =
 318   Printf.sprintf "case %d: *ret = (%s) - (%s); break;" j (List.nth i j) p;;
 319
 320 let vardecl y vs =
 321   let varname = map (fun (a,b,c) -> b) vs in
 322   let nvs = List.length vs in
 323   let  v j = Printf.sprintf "double %s = %s[%d];"   (List.nth varname j) y j in
 324     join_lines (map v (0-- (nvs-1)));;
 325
 326 let bounds f vs =
 327   let lbs = map f vs in
 328   join_comma lbs;;
 329
 330 let rec geteps =
 331   let getepsf = function
 332     Eps t -> t
 333     | _ -> 0.0
 334   in function
 335       [] -> 0.0
 336   | b::bs -> max(getepsf b) (geteps bs);;
 337
 338 let (has_penalty,penalty) =
 339   let penalties iqd =
 340     filter (function Penalty _ -> true | _ -> false) iqd.tags in
 341   let hasp iqd = (List.length (penalties iqd) >0) in
 342   let onep iqd = if (hasp iqd) then hd(penalties iqd) else Penalty (0.0,0.0) in
 343   (hasp,onep);;
 344
 345 let penalty_var iqd =
 346    let penalty_ub = function Penalty(a,_) -> string_of_float a in
 347    ["0.0","penalty",penalty_ub (penalty iqd)];;
 348
 349 let penalty_wt iqd = if has_penalty iqd then
 350   match (penalty iqd) with
 351       Penalty(_,b) -> (string_of_float b)^" * penalty"
 352 else "0.0";;
 353
 354 let rec cfsqp_branch = function
 355   | [] -> 0
 356   | Cfsqp_branch i ::_ -> i
 357   | _::a -> cfsqp_branch a;;
 358
 359 let move_first i ls =
 360   let (a,b::xs) = chop_list i ls in
 361   b::(a @ xs);;
 362
 363 let cc_main =
 364 "\n\nint main(){
 365   //Mathematica generated test data
 366   assert(near (pi(),4.0*atan(1.0)));
 367   assert(near (sqrt2(),1.41421356237309));
 368   assert(near (sqrt8(),2.828427124746190));
 369   assert(near (sol0(),0.5512855984325308));
 370   assert(near (tau0(),1.54065864570856));
 371   assert(near (acos(0.3),1.26610367277949));
 372   assert(near(hminus(),1.2317544220903216));
 373   assert(near(hplus(),1.3254));
 374   assert(near(mm1(),1.012080868420655));
 375   assert(near(mm2(),0.0254145072695089));
 376   assert(near(real_pow(1.18,2.),1.3924));
 377   assert(near(marchal_quartic(1.18),0.228828103048681825));
 378   assert(near(lmfun(1.18),0.30769230769230793));
 379   assert(near(lmfun(1.27),0.0));
 380   assert(near(rad2_x(4.1,4.2,4.3,4.4,4.5,4.6),1.6333363881302794));
 381   assert(near(dih_y(2.1,2.2,2.3,2.4,2.5,2.6),1.1884801338917963));
 382   assert(near(sol_y(2.1,2.2,2.3,2.4,2.5,2.6), 0.7703577405137815));
 383   assert(near(sol_y(2, 2, 2, 2.52, 3.91404, 3.464),4.560740765722419));
 384   assert(near(taum(2.1,2.2,2.3,2.4,2.5,2.6),tau_m(2.1,2.2,2.3,2.4,2.5,2.6) ));
 385   assert(near(taum(2.1,2.2,2.3,2.4,2.5,2.6),tau_m_alt(2.1,2.2,2.3,2.4,2.5,2.6) ));
 386   assert(near(taum(2.1,2.2,2.3,2.4,2.5,2.6),0.4913685097602183));
 387   assert(near(taum(2, 2, 2, 2.52, 3.91404, 3.464),4.009455167289888));
 388   assert(near(ups_x(4.1,4.2,4.3), 52.88));
 389   assert(near(eta_y(2.1,2.2,2.3), 1.272816758217772));
 390   assert(near(beta_bump_force_y(2.1,2.2,2.3,2.4,2.5,2.6),
 391     -0.04734449962124398));
 392   assert(near(beta_bump_force_y(2.5,2.05,2.1,2.6,2.15,2.2),
 393     beta_bumpA_y(2.5,2.05,2.1,2.6,2.15,2.2)));
 394   assert(near(atn2(1.2,1.3),atan(1.3/1.2)));
 395   assert(near(edge_flat(2.1,2.2,2.3,2.5,2.6),4.273045018670291));
 396   assert(near(flat_term(2.1),-0.4452691371955056));
 397   assert(near(enclosed(2.02,2.04,2.06,2.08,2.1,2.12,2.14,2.16,2.18),
 398     3.426676872737882));
 399 }\n\n";;
 400
 401 let cc_code outs iqd =
 402   let (b,vs,i) = c_dest_ineq iqd.ineq in
 403   let vs = vs @ if (has_penalty iqd) then penalty_var iqd else [] in
 404   let branch = cfsqp_branch iqd.tags in
 405   let i = move_first branch i in
 406   let eps = geteps (iqd.tags) in
 407   let casep = if has_penalty iqd then "max(0.0,penalty)" else "0.0" in
 408   let nvs = List.length vs in
 409   let ni = List.length i in
 410   let y = "y_mangle__" in
 411   let p = Printf.sprintf in
 412   let s = join_lines ([
 413     p"// This code is machine generated ";
 414     p"#include <iomanip.h>\n#include <iostream.h>\n#include <math.h>";
 415     p"#include \"../Minimizer.h\"\n#include \"../numerical.h\"";
 416     p"#include <assert.h>";
 417     p"class trialdata { public: trialdata(Minimizer M,char* s)";
 418     p"{ M.coutReport(s); };};";
 419     p"int trialcount = %d;\n"  (!cfsqp_trialcount);
 420     join_lines(map cc_function (!autogen));
 421     p"void c0(int numargs,int whichFn,double* %s, double* ret,void*) {" y;
 422     vardecl y vs ;
 423     p"switch(whichFn) {";
 424     ] @ map (case casep i) (1-- (-1 + List.length i)) @ [
 425     p"default: *ret = 0; break; }}\n\n";
 426     p"void t0(int numargs,int whichFn,double* %s, double* ret,void*) { " y;
 427     vardecl y vs ;
 428     p"*ret = (%e) + (%s) + (%s);" eps (List.nth i 0) (penalty_wt iqd);
 429     p"}";
 430     p"Minimizer m0() {";
 431     p"  double xmin[%d]= {" nvs;(bounds (function (a,b,c) -> a) vs);
 432     p "};\n  double xmax[%d]= {" nvs; (bounds (function (a,b,c) -> c) vs);
 433     p "};\n     Minimizer M(trialcount,%d,%d,xmin,xmax);" nvs (ni-1);
 434     p"  M.func = t0;";
 435     p"  M.cFunc = c0;";
 436     p"  return M;";
 437     p"}";
 438     p "trialdata d0(m0(),\"%s\");\n\n"  iqd.idv;
 439     p"int near(double x,double y)";
 440     p"  { double eps = 1.0e-8; return (mabs(x-y)<eps); }";
 441     ]) in
 442   Printf.fprintf outs "%s %s" s cc_main;;
 443
 444 let mk_cc tmpfile iqd =
 445   let outs = open_out tmpfile in
 446   let _ = cc_code outs iqd in
 447    close_out outs ;;
 448
 449 let compile =
 450   let err = Filename.temp_file "cfsqp_err" ".txt" in
 451     fun () ->
 452       let e = Sys.command("cd "^flyspeck_dir^"/../cfsqp; make tmp/t.o >& "^err) in
 453       let _ =   (e=0) or (Sys.command ("cat "^err); failwith "compiler error ") in
 454         ();;
 455
 456  let execute_cfsqp =
 457   let cfsqp_out = Filename.temp_file "tmp_cfsqp" ".out" in
 458   let cfsqp_dir = flyspeck_dir^"/../cfsqp" in
 459     fun idq ->
 460       let _ =  mk_cc (cfsqp_dir ^ "/tmp/t.cc") idq in
 461       let _ = try (compile()) with Failure s -> failwith (s^idq.idv) in
 462       let _ = (0=  Sys.command(cfsqp_dir^"/tmp/t.o | tee "^cfsqp_out)) or
 463         failwith ("execution error "^idq.idv) in
 464       let s = process_to_string ("grep 'NEGATIVE' "^cfsqp_out) in
 465       let _ =  (s="") or failwith ("NEGATIVE "^idq.idv) in
 466         ();;
 467
 468  let execute_cfsqp_list xs =
 469    let fails = ref [] in
 470    let _ = for i=0 to (List.length xs - 1)
 471       do
 472         try  (execute_cfsqp (List.nth xs i); Sys.command("date"); ())
 473         with Failure s -> (fails := s :: (!fails)) done in
 474      fails;;
 475
 476
 477 (* ========================================
 478
 479 glpk  generation
 480 ======================================== *)
 481
 482
 483  let yy6 =  [`y1:real`;`y2:real`;`y3:real`;`y4:real`;`y5:real`;`y6:real`];;
 484  let yy9 = yy6 @ [`y7:real`;`y8:real`;`y9:real`];;
 485
 486  let  glpk_translate =ref
 487      [("dih_y","azim[i,j]");("dih2_y","azim2[i,j]");("dih3_y","azim3[i,j]");
 488       ("sol_y","sol[j]");
 489       ("taum","tau[j]");("tauq","tau[j]");("rhazim","rhazim[i,j]");
 490       ("rhazim2","rhazim2[i,j]");("rhazim3","rhazim3[i,j]");
 491       ("sqrt8","sqrt8")];;
 492
 493  let reverse =
 494      let xs =     [("y2","y3");("y5","y6");("rhazim2","rhazim3");
 495       ("dih2_y","dih3_y");] in
 496      let f (a,b) = (b,a) in
 497      let ys =     xs @ map f xs in
 498        fun s -> try (List.assoc s ys) with Not_found -> s;;
 499
 500   let glpk_lookup rev s xs =
 501     let s' = if rev then reverse s else s in
 502     if (xs = yy6) or (xs = yy9)
 503     then
 504       try
 505         (List.assoc s' (!glpk_translate))
 506       with Not_found -> failwith ("glpk_lookup translate" ^ s)
 507     else if xs = []
 508     then
 509       try
 510         (List.assoc s' (!glpk_translate))
 511       with Not_found -> (s'^"[i,j]")
 512     else failwith ("glpk_lookup unknown arg list:" ^ s);;
 513
 514  let rec glpk_form rev t =
 515   let soh = glpk_form rev in
 516   if is_var t then glpk_lookup rev (fst (dest_var t)) [] else
 517   let (f,xs) = strip_comb t in
 518   let ifix i = let [a;b] = xs in paren(soh a ^ " " ^ i ^ " " ^ soh b) in
 519   let ifix_swapped i =
 520     let [b;a] = xs in paren(soh a ^ " " ^ i ^ " " ^ soh b) in
 521   let _ = (is_const f) or
 522       failwith ("glpk constant expected: " ^ string_of_term f) in
 523   match fst (dest_const f) with
 524   | "real_gt" | "real_ge" | "real_sub" -> ifix "-"
 525   | "real_lt" | "real_le" -> ifix_swapped "-"
 526   | "real_add" -> ifix "+"
 527   | "real_mul" -> ifix "*"
 528   | "real_div" -> ifix "/"
 529   | "real_neg" -> let [a] = xs in paren("-" ^ soh a)
 530   | "real_of_num" -> let [a] = xs in soh a
 531   | "NUMERAL" -> let [_] = xs in string_of_num' (dest_numeral t)
 532   | "DECIMAL" ->  string_of_num' (dest_decimal t)
 533   | s -> paren(glpk_lookup rev s xs) ;;
 534
 535 let (counter,counter_reset) =
 536       let ineqcounter = ref 0 in
 537       let counter t =
 538         let u = !ineqcounter in
 539         let _ =  ineqcounter := u+1 in
 540           u in
 541       let counter_reset _ =  (ineqcounter := 0) in
 542         (counter,counter_reset);;
 543
 544 let mk_glpk_ineq rev iqd =
 545     let ineq = iqd.ineq in
 546   let t = snd(strip_forall  (ineq)) in
 547   let (vs,i) = dest_comb t  in
 548   let (_,vs) = dest_comb vs in
 549   let (f,xs) = strip_comb vs in
 550   let (dart,_) = dest_const f in
 551   let i' = hd(disjuncts i) in
 552   let _ = (xs = yy6) or (xs = yy9) or
 553     (print_qterm vs; failwith "dart_class y1 ... y6 expected") in
 554   let p = Printf.sprintf in
 555   let s =   p"ineq%d 'ID[%s]' \n  { (i,j) in %s } : \n  %s >= 0.0;\n\n"
 556     (counter()) iqd.idv dart (glpk_form rev i') in
 557     s;;
 558
 559 let mk_glpk_ineq_id rev s =
 560   let iqd = Ineq.getexact s in mk_glpk_ineq rev (hd iqd);;
 561
 562 let test_domain_symmetry idq =
 563   let ineq = idq.ineq in
 564   let dart = List.nth (snd(strip_comb (snd(strip_forall ineq)))) 0 in
 565   let dom = map(fun (a,(b,c)) -> (a,c)) (map (fun (a,b) -> (a,dest_pair b))
 566      (map dest_pair
 567      (dest_list(snd(dest_eq(
 568                       concl(REWRITE_CONV(!Ineq.dart_classes) dart))))))) in
 569   let nth = List.nth in
 570     ((nth dom 1 = nth dom 2) & (nth dom 4 = nth dom 5)) or
 571       failwith "domain asym";;
 572
 573 let lpstring() =
 574   let _ = counter_reset() in
 575   let _ = map test_domain_symmetry (Ineq.getfield Lpsymmetry) in
 576   join_lines
 577     (("# File automatically generated from nonlinear inequality list "^
 578        "via lpstring().\n\n") ::
 579     (map (mk_glpk_ineq false) (Ineq.getfield Lp)) @
 580     ["# Symmetry section\n\n"] @
 581     (map (mk_glpk_ineq true) (Ineq.getfield Lpsymmetry)));;
 582
 583
 584
 585 (* ========================================
 586
 587 Objective Caml  generation
 588 ======================================== *)
 589
 590
 591 let ocaml_string_of_term t =
 592   let rec soh t =
 593     if is_var t then fst (dest_var t) else
 594       let (f,xs) = strip_comb t in
 595       let ifix i = let [a;b] = xs in paren(soh a ^ " " ^ i ^ " " ^ soh b) in
 596       let fv =
 597         if is_var f
 598         then
 599           let fv = fst(dest_var f) in
 600           let _ = warn true ("variable function name: "^fv) in
 601             fv
 602         else if is_const f then fst (dest_const f)
 603         else
 604           failwith ("var or const expected: " ^ string_of_term f) in
 605         match fv with
 606           | "real_sub" -> ifix "-. "
 607           | "real_add" -> ifix "+."
 608           | "real_mul" -> ifix "*."
 609           | "real_div" -> ifix "/."
 610           | "real_lt" -> ifix "<"
 611           | "real_le" -> ifix "<="
 612           | "real_gt" -> ifix ">"
 613           | "real_ge" -> ifix ">="
 614           | "real_pow" -> ifix " ** "
 615           | "/\\" -> ifix "&&"
 616           | "\\/" -> ifix "or"
 617           | "real_neg" -> let [a] = xs in paren("-. " ^ soh a)
 618           | "acs" ->
 619               (try
 620                 let [a] = xs in "(acos("^soh a^ "))"
 621               with Match_failure _ ->
 622                 let _ = (xs=[]) or failwith "ocaml:acs" in
 623                   paren "acos")
 624           | "real_of_num" -> let [a] = xs in soh a
 625           | "NUMERAL" -> let [_] = xs in string_of_num' (dest_numeral t)
 626           | "=" -> let [a;b] = xs in soh a ^ " = " ^ soh b
 627           | "," -> let [a;b] = xs in paren( soh a ^ "," ^ soh b)
 628           | "DECIMAL" ->  string_of_num' (dest_decimal t)
 629           | "COND" ->
 630               let [a;b;c] = xs in
 631                 paren("if "^( paren(soh a)) ^
 632                         " then " ^ (paren (soh b)) ^
 633                         " else " ^ (paren(soh c)) )
 634           | "atn2"      -> let [ab] = xs in let (a,b) = dest_pair ab in
 635               "atan2 " ^ paren(soh b) ^ " " ^ paren(soh a) ^ "  "
 636                 (* note reverse order atan2 in ocaml *)
 637           | "domain6" -> let [h;f;g] = xs in domain6 soh h f g
 638           | s -> " " ^ s ^ " " ^ join_space(map (paren o soh) xs) ^ " "
 639   and domain6 soh h f g =
 640     let (hv,hbody) = strip_abs h in
 641     let (f1,xs) = strip_comb f in
 642     let vs a = join_space (map (fst o dest_var) a) in
 643     let hbodys = soh hbody in
 644     let fs = fst(dest_const f1) in
 645     let gs = soh g in
 646       (Printf.sprintf
 647          " %s %s = \n  let _ = ( %s ) or failwith \"domain6:%s\" in \n
 648       ( %s ) %s" fs (vs (xs @ hv)) hbodys fs gs (vs hv)) in
 649     try
 650       (soh t)
 651     with Failure s -> failwith (s^" .......   "^string_of_term t);;
 652
 653 let ocaml_function t =
 654   "let " ^ ocaml_string_of_term t ^ ";;\n\n";;
 655
 656 let ocaml_autogen = map
 657   (function b -> snd(strip_forall (concl (strip_let b))))
 658   [sqrt2;sqrt8;delta_x;delta_y;delta_x4;delta_x6;ups_x;eta_x;
 659   eta_y;dih_x;dih_y;dih2_y;dih3_y;sol_x;sol_y;interp;ly;const1;
 660   rho;rhazim;lnazim;taum;arclength];;
 661
 662 let ocaml_code =
 663   "(* code automatically generated from Parse_ineq.ocaml_code *)\n\n"^
 664    "module Sphere_math = struct\n\n"^
 665    "let sqrt = Pervasives.sqrt\n\n" ^
 666    "let pi = 4.0 *. atan(1.0);;\n" ^
 667    (join_lines (map ocaml_function ocaml_autogen)) ^
 668    "end;;\n";;
 669
 670 let sphere_ml = Filename.temp_file "sphere" ".ml";;
 671
 672 output_filestring sphere_ml ocaml_code;;
 673
 674 (* ========================================
 675
 676 LaTeX  generation
 677 ======================================== *)
 678
 679 let mk_texstring (tm,idv,s) = match tm with
 680   | Section -> "\\section{"^idv^" "^s^"}\n"
 681   | Comment -> "\\begin{remark}["^idv^"]\n"^s^"\n\\end{remark}\n"
 682   | Ineqdoc -> "\\begin{calculation}["^idv^"]\n"^s^"\n\\end{calculation}\n";;
 683
 684 let texstring() = join_lines (map mk_texstring (List.rev (!Ineq.ineqdoc)));;
 685
 686 (*
 687 output_filestring "/tmp/ineqdoc.tex" (texstring());;
 688 *)
 689
 690
 691 end;;