2 flyspeck_needs "nonlinear/prep.hl";;
3 flyspeck_needs "nonlinear/scripts.hl";;
5 module Break_case_exec = struct
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 ). *)
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.
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.
25 Iarg_leaf represents a domain with no subpartition. The domain is all of (x,z).
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].
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.
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).
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.
50 type ('a) iargs = Iarg_facet of ((int * bool) * float * 'a * ('a) iargs)
52 | Iarg_bisect of (int * ('a) iargs * ('a) iargs);;
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);;
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);;
66 let sprintf = Printf.sprintf;;
70 let getprep s = hd(filter (fun t -> t.idv = s) (!Prep.prep_ineqs));;
72 (* omit the quad cases for now: *)
75 let nonquad = filter (fun t -> not(Optimize.is_quad_cluster t.tags)) (!Prep.prep_ineqs) in
76 map (fun t -> t.idv) nonquad;;
83 | a::rest -> (map (fun x -> (a,x)) b) @ cart rest b;;
85 let maxlist xs = List.fold_right max xs (nth xs 0);;
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)
95 let msec_inc = float_of_int MSEC_INC;;
97 let float_cache = ref (fun() -> 0.0);;
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
106 (* C++ CODE GENERATION. *)
108 let cpp_template_arg = sprintf "
109 const char svn[] = %s;
110 const char ineq_id[] = %s;
112 int testRun(double x1[6],double z1[6]) // autogenerated code
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);
124 const Function* I[%d] = {%s}; // len ...
126 opt.allowSharp = %d; // sharp
127 opt.onlyCheckDeriv1Negative = %d; // checkderiv
129 return prove::recursiveVerifier(0,x,z,x0,z0,I,%d,opt); // len
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;;
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
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]);;
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
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);;
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
186 (* DOMAIN CONSTANTS *)
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
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);;
201 let get_constants s =
202 let (x,y) = get_constants_xy s in x @ y;;
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
209 (* they all check out on 2013/7/31
210 find (fun s -> not(can get_float_domain s)) (snd (chop_list 650 idvlist));;
214 (* interpolate along edge *)
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
223 let facet_opp (i,side)= (i,not side);;
225 let invert_side (_,side) (x,z) =
226 if side then (z,x) else (x,z);;
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);;
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;;
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
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)
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)
261 let _ = f1+. minwidth <= f2 or failwith "find_frac width underflow" in
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
272 let t3 = float_of_int m3 in
273 if abs_float(t3/. msec_inc -. 1.0) < ERROR_TOLERANCE then (f3,m3)
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');;
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;;
285 let rec recheck_pass timeout domain iarg = match iarg with
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);;
304 (* MAIN PROCEDURES *)
307 let ft_sec = Scripts.finalize Scripts.finished_times in
308 let ft_msec = Scripts.finalize Scripts.finished_times_msecs in
311 with _ -> 1000 * (assoc s ft_sec );;
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
321 let rec pass_revised avoid ((x,z),msec) =
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
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
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
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))
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)))
359 let _ = report ("failure "^fail) in
360 pass_revised (idx::avoid) ((x,z),msec);;
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
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