text_formalization/nonlinear/break_case_exec.hl

   1
   2 flyspeck_needs "nonlinear/prep.hl";;
   3 flyspeck_needs "nonlinear/scripts.hl";;
   4
   5 module Break_case_exec = struct
   6
   7 (* open Sphere2;; *)
   8
   9
  10 (* PRELIMINARIES *)
  11
  12 let ( MSEC_INC ) = 1000;; (* global parameter for msec target size *)
  13 let ERROR_TOLERANCE = 0.3;;  (* try for times within MSEC_INC * (1 +/- ERROR_TOLERANCE ). *)
  14 let verbose = true;;
  15
  16 (*
  17
  18 The type iargs is used to record a partition of a nonlinear inequality domain into pieces.
  19 The code in this file is used to make the pieces all take approximately MSEC_INC milliseconds
  20 to run.  The field 'a records the interval arithmetic C++ runtime in milliseconds.
  21
  22 The partition is determined from the data and the original domain determined by (x,z), where
  23 x and z are lists of real numbers giving the lower left and upper right hand corners of the domain.
  24
  25 Iarg_leaf represents a domain with no subpartition.  The domain is all of (x,z).
  26
  27 Iarg_bisect (i,left,right) is a cut exactly in the middle along the ith variable.
  28 The fields left and right give the further partition of the left and right sides.
  29 For example if i=1 and (x,z) = [1.0;1.0;1.0],[3.0;3.0;3.0] then the partition goes
  30 into (xleft,zleft) = [1.0;1.0;1.0],[3.0;2.0;3.0] and (xright,zright)=[1.0;2.0;1.0],[3.0;3.0;3.0].
  31
  32 Iarg_facet ((i,side),frac,msec,residual) breaks the domain in two unequal parts along the ith coordinate
  33 specified by a fraction 0.0 <= frac <= 1.0.
  34 The fraction is measured from the left if side=false, and from the right otherwise.
  35
  36 For example, if i,side=1,false frac=0.25 and (x,z)= [x0;x1;x2],[z0;z1;z2], then
  37 (xleft,rleft)= x,[z0;0.75 * x1+0.25 *z1;z2] and (xright,zright)=[x0;0.75*x1+0.25*z1,z].
  38 The number msec is the C++ runtime for (xleft,rleft), and residual is the recursively defined
  39 partition on (xright,zright).
  40
  41 For example, if we change side=true and keep the rest of the data the same, then
  42 (xleft,rleft)= x,[z0;0.25 * x1+0.75 *z1;z2] and (xright,zright)=[x0;0.25*x1+0.75*z1,z].
  43 The number msec is now the C++ runtime for (xright,right), and residual is the recursively defined
  44 partition on (xleft,zleft). So changing the side exchanges right and left all the way along.
  45
  46
  47 *)
  48
  49
  50 type ('a) iargs = Iarg_facet of ((int * bool) * float * 'a * ('a) iargs)
  51              | Iarg_leaf of 'a
  52              | Iarg_bisect of (int * ('a) iargs * ('a) iargs);;
  53
  54 let rec string_of_iargs = function
  55   | Iarg_leaf msec -> Printf.sprintf "\n Iarg_leaf %d" msec
  56   | Iarg_facet ((i,b),frac,msec,a) -> Printf.sprintf "\n Iarg_facet ((%d,%s),%3.4f,%d,%s)"
  57       i (string_of_bool b) frac msec (string_of_iargs a)
  58   | Iarg_bisect(idx,a,b) -> Printf.sprintf "\n Iarg_bisect (%d,%s,%s)"
  59       idx (string_of_iargs a) (string_of_iargs b);;
  60
  61 let output_case s iargs =
  62   let sfull = Printf.sprintf "\nadd_case (\"%s\",%s);;\n" s (string_of_iargs iargs) in
  63   let oc = open_out_gen [Open_append;Open_text] 436 (fullpath "nonlinear/break_case_log_more.hl") in
  64   (Pervasives.output_string oc (sfull); close_out oc);;
  65
  66 let sprintf = Printf.sprintf;;
  67
  68 Random.init 0;;
  69
  70 let getprep s = hd(filter (fun t -> t.idv = s) (!Prep.prep_ineqs));;
  71
  72 (* omit the quad cases for now: *)
  73
  74 let idvlist =
  75   let nonquad = filter (fun t -> not(Optimize.is_quad_cluster t.tags))  (!Prep.prep_ineqs) in
  76     map (fun t -> t.idv) nonquad;;
  77
  78 let nth = List.nth;;
  79
  80 let rec cart a b =
  81   match a with
  82     | [] -> []
  83     | a::rest -> (map (fun x -> (a,x)) b) @ cart rest b;;
  84
  85 let maxlist xs = List.fold_right max xs (nth xs 0);;
  86
  87 let rec trim s =
  88   let white c = mem c [' '; '\012'; '\n'; '\r'; '\t'] in
  89   let n = String.length s in
  90   let subs k = String.sub s k (n-1) in
  91     if (n > 0 && white (s.[0])) then trim (subs 1)
  92     else if (n > 1 && white (s.[n-1])) then trim (subs 0)
  93     else s;;
  94
  95 let msec_inc = float_of_int MSEC_INC;;
  96
  97 let float_cache = ref (fun() -> 0.0);;
  98
  99 let eval_float s =
 100     let (b,r) = Flyspeck_lib.eval_command ~silent:false
 101       ("float_cache := (fun () -> ("^s^"));;") in
 102     let _ = b or (print_string (r^"\n"^s^"\n"); failwith "bad input string") in
 103     let t= (!float_cache)() in
 104       t;;
 105
 106 (* C++ CODE GENERATION. *)
 107
 108 let cpp_template_arg = sprintf "
 109  const char svn[] = %s;
 110  const char ineq_id[] = %s;
 111
 112  int testRun(double x1[6],double z1[6]) // autogenerated code
 113         {
 114         // Warning: not rigorous. The rounding is off by epsilon. Use this only for experiments.
 115         interval tx[6]={interval(x1[0],x1[0]),interval(x1[1],x1[1]),interval(x1[2],x1[2]),
 116                         interval(x1[3],x1[3]),interval(x1[4],x1[4]),interval(x1[5],x1[5])   };
 117         interval tz[6]={interval(z1[0],z1[0]),interval(z1[1],z1[1]),interval(z1[2],z1[2]),
 118                         interval(z1[3],z1[3]),interval(z1[4],z1[4]),interval(z1[5],z1[5])};
 119         domain x = domain::lowerD(tx);
 120         domain z = domain::upperD(tz);
 121         domain x0=x;
 122         domain z0=z;
 123         %s
 124         const Function* I[%d] = {%s}; // len ...
 125         cellOption opt;
 126         opt.allowSharp = %d; // sharp
 127         opt.onlyCheckDeriv1Negative = %d; // checkderiv
 128         %s // other options.
 129         return  prove::recursiveVerifier(0,x,z,x0,z0,I,%d,opt); // len
 130         }";;
 131
 132 let mk_cpp_arg_proc t s tags =
 133   let sharp = if  mem Sharp tags then 1 else 0 in
 134   let checkderiv = if  mem Onlycheckderiv1negative tags then 1 else 0 in
 135   let ifd b s = if b then s else "" in
 136   let (b,f) = Optimize.widthCut tags in
 137   let sWidth = ifd b (sprintf "\topt.widthCutoff = %8.16f;\n" f) in
 138   let c = map Optimize.cpp_string_of_term in
 139   let f (x,y,z) = (c x,c y,c z) in
 140   let (aas,bbs,iis) = f (Optimize.dest_nonlin t) in
 141   let len = length iis in
 142   let sq = Optimize.quoted s in
 143   let svn = (Optimize.quoted(Optimize.svn_version())) in
 144     cpp_template_arg svn sq (Optimize.cpp_template_t "" iis)
 145       len (Optimize.cpp_template_Fc "" len) sharp  checkderiv sWidth len;;
 146
 147 let mkfile_arg =
 148   let cpp_tail = Optimize.join_lines (Optimize.load_file  (flyspeck_dir^"/../interval_code/arg_tail.txt")) in
 149   let cpp_header = Optimize.cpp_header() in
 150     fun t s tags  ->
 151         Flyspeck_lib.output_filestring Optimize.tmpfile
 152           (Optimize.join_lines [cpp_header;Auto_lib.interval_code;(mk_cpp_arg_proc t s tags);cpp_tail]);;
 153
 154 let execute_args ex tags s testineq xlist zlist =
 155    let x = List.nth xlist in
 156    let z = List.nth zlist in
 157    let args = sprintf " %f %f %f %f %f %f   %f %f %f %f %f %f"
 158       (x 0) (x 1) (x 2) (x 3) (x 4) (x 5)     (z 0) (z 1) (z 2) (z 3) (z 4) (z 5) in
 159   let interval_dir = flyspeck_dir^"/../interval_code" in
 160   let _ = mkfile_arg testineq s tags in
 161   let _ = Optimize.compile_cpp() in
 162   let _ = (not ex) or (0=  Sys.command(interval_dir^"/test_auto"^args)) or failwith "interval execution error" in
 163     ();;
 164
 165 let process_and_prep_args ex (s,tags,case) =
 166   let _ = report ("process and prep args: "^s) in
 167   let (s,tags,testineq) = (* preprocess debug *) (s,tags,case) in
 168   let (x,y,_) = Optimize.dest_nonlin testineq in
 169     (execute_args ex tags s testineq , x, y);;
 170
 171 let rerun_timer =
 172   let run_out = Filename.temp_file "run" ".out" in
 173     fun xlist zlist timer ->
 174       let x = List.nth xlist in
 175       let z = List.nth zlist in
 176       let args = sprintf " %f %f %f %f %f %f   %f %f %f %f %f %f %f"
 177         (x 0) (x 1) (x 2) (x 3) (x 4) (x 5)     (z 0) (z 1) (z 2) (z 3) (z 4) (z 5) (float_of_int timer) in
 178       let interval_dir = flyspeck_dir^"/../interval_code" in
 179       let _ = (0=  Sys.command(interval_dir^"/test_auto"^args^" | tee "^run_out)) or
 180         failwith "interval execution error" in
 181       let outs = trim (process_to_string ("grep msecs "^run_out^" | sed 's/^.*msecs=//' | sed 's/;.*$//' ")) in
 182         (*  let _ = report (":"^outs^":") in *)
 183       let msecs = try (int_of_string (outs)) with _ -> timer in
 184         (msecs);;
 185
 186 (* DOMAIN CONSTANTS *)
 187
 188 let scrub_c t =
 189   let th = REWRITE_CONV [Sphere.h0;
 190         GSYM Nonlinear_lemma.sol0_over_pi_EQ_const1;
 191         ASSUME `hminus = #1.2317544220903216`;
 192         ASSUME `pi = #3.1415926535897932385`;
 193         ASSUME `sol0 = #0.55128559843253080794`] t in
 194     rhs(concl th);;
 195
 196 let get_constants_xy s =
 197   let (s,tags,case) = Optimize.idq_fields (getprep s) in
 198   let (x,y,_) = Optimize.dest_nonlin case in
 199     (map scrub_c x,map scrub_c y);;
 200
 201 let get_constants s =
 202   let (x,y) = get_constants_xy s in x @ y;;
 203
 204 let get_float_domain s =
 205   let m = map (eval_float o Parse_ineq.ocaml_string_of_term) in
 206   let (x,y) = get_constants_xy s in
 207     (m x,m y);;
 208
 209 (* they all check out on 2013/7/31
 210 find (fun s -> not(can get_float_domain s)) (snd (chop_list 650 idvlist));;
 211 *)
 212
 213
 214 (* interpolate along edge *)
 215
 216 let index_max_width (x ,z) avoid =
 217   let w = map (fun (xi,zi) -> zi -. xi) (zip x z) in
 218   let avoid_filter = map (fun i -> (mem i avoid)) (0--(List.length x - 1)) in
 219   let w' = map (fun (t,b) -> if b then 0.0 else t) (zip w avoid_filter) in
 220   let wm = maxlist w' in
 221     index wm w';;
 222
 223 let facet_opp (i,side)= (i,not side);;
 224
 225 let invert_side (_,side) (x,z) =
 226   if side then (z,x) else (x,z);;
 227
 228 let invert_interpolate (i,side) frac (x,z) =
 229   let _ = (0.0 <= frac && frac <= 1.0) or failwith "invert_interpolate: frac out of range " in
 230   let _ = List.length x = List.length z or failwith "invert_interpolate: length mismatch" in
 231   let n = List.length x in
 232   let rg = (0--(n-1)) in
 233   let (x,z) = invert_side (i,side) (x,z) in
 234   let c = nth x i *. (1.0-. frac) +. nth z i *. frac in
 235   let modi (xi,j) = if (i=j) then c else xi in
 236   let xm = map modi (zip x rg) in
 237   let zm = map modi (zip z rg) in
 238     (* partition will be x--zm,  xm--z *)
 239     invert_side(i,side) (xm,zm);;
 240
 241 let interpolate_frac (f1,t1) (f2,t2) =
 242   let _ = (0.0 <= f1 && f1 < f2 && f2 <= 1.0 && t1 <= msec_inc && msec_inc <= t2)
 243     or failwith "interpolate: out of range" in
 244     (msec_inc -. t1) *. (f2 -. f1) /. (t2 -. t1) +. f1;;
 245
 246 let run_frac frac facet (x1,z1) timeout =
 247   let (xm,zm) = invert_interpolate facet frac (x1,z1) in
 248   let m3 = if (not(snd facet)) then rerun_timer x1 zm timeout else rerun_timer xm z1 timeout in
 249     ((xm,zm),m3);;
 250
 251 let rec recursive_find_frac facet x1 z1  (f1,t1) (f2,t2) =
 252   (* last case in a batch may have a small t2, so dont use abs_float in the following line *)
 253   if (t2/. msec_inc -. 1.0) < ERROR_TOLERANCE then (f2,int_of_float t2)
 254   else if abs_float(t1/. msec_inc -. 1.0) < ERROR_TOLERANCE then (f1,int_of_float t1)
 255   else
 256
 257     let minwidth = 0.001 in
 258       if (f2 <= f1 +. minwidth && t2 < 1.99 *. msec_inc) then (f2,int_of_float t2)
 259       else if (f2 <= f1 +. minwidth && t1 > 0.25 *. msec_inc) then (f1,int_of_float t1)
 260       else
 261         let _ = f1+. minwidth <= f2 or failwith "find_frac width underflow" in
 262
 263         let f3 = interpolate_frac (f1,t1) (f2,t2) in
 264         let f3 = if (Random.int 4 = 0 && t2 > 1.9*. msec_inc) then
 265           (f1+.f2)/. 2.0 +. Random.float (f2 -. f1) /. 2.0 else f3 in
 266         let f3 = if (Random.int 3 = 0 && t1 < 0.5*. msec_inc) then
 267           f1 +. Random.float (f2 -. f1)/. 2.0 else f3 in
 268         let _ = (f1 <= f3 && f3 <= f2) or failwith "recursive_find_frac: out of range" in
 269         let ((xm,zm),m3) = run_frac f3 facet (x1,z1) (2 * MSEC_INC) in
 270         let _ = if verbose then report (sprintf "recursing fracs %3.3f,%3.3f,%3.3f msecs %d,%d,%d" f1 f3 f2 (int_of_float t1) m3 (int_of_float t2)) in
 271
 272         let t3 = float_of_int m3 in
 273           if abs_float(t3/. msec_inc -. 1.0) < ERROR_TOLERANCE then (f3,m3)
 274           else
 275             let ((f1',t1'),(f2',t2')) = if (t3 <= msec_inc) then ((f3,t3),(f2,t2)) else ((f1,t1),(f3,t3)) in
 276               recursive_find_frac facet x1 z1 (f1',t1') (f2',t2');;
 277
 278 (* DEBUG *)
 279
 280 let rec pass_time = function
 281     Iarg_bisect (_,a,b) -> pass_time a + pass_time b
 282   | Iarg_leaf msec -> msec
 283   | Iarg_facet (_,_,msec,c) -> msec + pass_time c;;
 284
 285 let rec recheck_pass  timeout domain iarg = match iarg with
 286   | Iarg_leaf msec ->
 287       let (x,z) = domain in
 288       let msec' = rerun_timer x z timeout in
 289         Iarg_leaf (msec,msec')
 290   | Iarg_facet(facet,frac,msec,b) ->
 291       let (x,z) = domain in
 292       let (xm,zm) = invert_interpolate facet frac domain in
 293       let residual = if (not(snd facet)) then (xm,z) else (x,zm) in
 294       let (_,msec') = run_frac frac facet domain timeout in
 295         Iarg_facet(facet,frac,(msec,msec'),  recheck_pass timeout residual b)
 296   | Iarg_bisect(idx,a,b) ->
 297       let (x,z) = domain in
 298       let facet = (idx,false) in
 299       let (xm,zm) = invert_interpolate facet 0.5 domain in
 300       let partA,partB = ((x,zm),(xm,z)) in
 301         Iarg_bisect(idx,recheck_pass timeout partA a,recheck_pass timeout partB b);;
 302
 303
 304 (* MAIN PROCEDURES *)
 305
 306 let initial_msec =
 307   let ft_sec = Scripts.finalize Scripts.finished_times in
 308   let ft_msec = Scripts.finalize Scripts.finished_times_msecs in
 309     fun s ->
 310       try assoc s ft_msec
 311       with _ -> 1000 * (assoc s ft_sec );;
 312
 313 let initialize s =
 314   let msec = initial_msec s in
 315   let (x1,z1)=get_float_domain s in
 316   let stc = Optimize.idq_fields (getprep s) in
 317   let (compile,_,_) = process_and_prep_args false stc in
 318   let _ = compile x1 z1 in
 319     (x1,z1,msec);;
 320
 321 let rec pass_revised avoid ((x,z),msec) =
 322   let left = false in
 323   let idx = index_max_width (x,z) avoid in
 324   let fi = float_of_int in
 325   let timeout = 2 * MSEC_INC in
 326   let (s0,s1,s2) = (0,1,2) in
 327   let facet = (idx,left) in
 328   let sz m =
 329     let t = fi m /. fi MSEC_INC in
 330       if t < 1.0 -. ERROR_TOLERANCE then s0 else if t < 1.999 then s1 else s2 in
 331     if msec < timeout  then
 332       let _ = if verbose then report "FOUND leaf..." in
 333       (Iarg_leaf msec)
 334     else
 335       let (_,m3A) = run_frac 0.5 facet (x,z) timeout in
 336       let (_,m3B) = run_frac 0.5 (facet_opp facet) (x,z) timeout in
 337       let (m3min,facet,swap) = if m3A<=m3B then (m3A,facet,false) else (m3B,facet_opp facet,true) in
 338       let smin = sz m3min in
 339         if (smin=s2) then
 340           let (xm,zm) = invert_interpolate (idx,left) 0.5 (x,z) in
 341           let partA,partB = ((x,zm),(xm,z)) in
 342           let _ = if verbose then report "FOUND bisection..." in
 343             Iarg_bisect (idx,pass_revised [] (partA,m3A),pass_revised [] (partB,m3B))
 344         else if (smin = s1) then
 345           let (xm,zm) = invert_interpolate facet 0.5 (x,z) in
 346           let residual= if (snd facet=left) then (xm,z) else (x,zm) in
 347           let m3max = if swap then m3A else m3B in
 348           let _ = if verbose then report "FOUND 0.5 facet..." in
 349             Iarg_facet(facet,0.5,m3min,pass_revised [] (residual,m3max))
 350         else (* smin=s0 *)
 351           try (
 352             let (fracC,mC) = recursive_find_frac facet x z (0.5,fi m3min) (1.0,fi msec) in
 353             let (xm,zm) = invert_interpolate facet fracC (x,z) in
 354             let residual= if (snd facet=left) then (xm,z) else (x,zm) in
 355             let mD = rerun_timer (fst residual) (snd residual) timeout in
 356             let _ = if verbose then report (sprintf "FOUND %3.3f,%d facet..." fracC mC) in
 357               Iarg_facet (facet,fracC,mC,pass_revised [] (residual,mD)))
 358           with Failure fail ->
 359             let _ = report ("failure "^fail) in
 360               pass_revised (idx::avoid) ((x,z),msec);;
 361
 362 let run_one save s =
 363   let (x1,z1,msec1) = initialize s in
 364   let iarg = (pass_revised [] ((x1,z1),msec1)) in
 365   let _  = if (save) then output_case s iarg in
 366     iarg;;
 367
 368 let run_all_over2sec save =
 369   let preptimes = map (fun s -> try s,initial_msec s with _ -> failwith s) (idvlist) in
 370   let over2sec =  (sort (fun (_,s) (_,t) -> s < t) (filter (fun (s,t) -> t > 2000) preptimes)) in
 371   let over2 = map fst over2sec in
 372   let alldata = map (fun s -> try s,run_one save s with Failure msg -> report s; failwith msg) over2 in
 373     alldata;;
 374
 375
 376
 377 end;;