1 (* port of recurse.cc *)
4 This is the code that verifies a disjunct of nonlinear inequalities.
5 The are given as a list (tf:tfunction list). If tf = [f1;....;fk], then
6 the list represents the inequality (f1 < 0 \/ f2 < 0 .... fk < 0).
8 The end user should only need to define a cell option,
9 and then call recursive_verifier, which recursively bisects the domain
10 until a partition of the domain is found on which verifier_cell gives
11 a pass on each piece of the partition.
15 flyspeck_needs "../port_interval/types.hl";;
16 flyspeck_needs "../port_interval/report.hl";;
17 flyspeck_needs "../port_interval/interval.hl";;
18 flyspeck_needs "../port_interval/univariate.hl";;
19 flyspeck_needs "../port_interval/line_interval.hl";;
20 flyspeck_needs "../port_interval/taylor.hl";;
22 module Recurse = struct
30 only_check_deriv1_negative : bool;
31 is_using_dihmax : bool;
32 is_using_bigface126 : bool;
35 allow_derivatives : bool;
36 mutable iteration_count : int;
37 iteration_limit : int;
38 recursion_depth : int;
41 (* cell verification is complex, and we use exceptions to
42 exit as soon as the status has been determined. *)
47 | Cell_inconclusive of (float list * float list * float list * float list * tfunction list);;
49 exception Return of cell_status;;
51 let return c = raise (Return c);;
54 (* error checking and reporting functions *)
56 let string_of_domain x =
58 Printf.sprintf "{%f, %f, %f, %f, %f, %f}" (n x 0) (n x 1) (n x 2) (n x 3) (n x 4) (n x 5);;
60 let string3 (x,z,s) = (string_of_domain x ^"\n"^ string_of_domain z ^ "\n" ^ s);;
62 let boolify _ = true;;
64 let report_current = boolify o Report.report_timed o string3;;
66 let report_error = boolify o Report.report_error o string3;;
68 let report_fatal = boolify o Report.report_fatal o string3;;
70 (* let t = [0.1;0.2;0.3;0.4;0.5;0.6] in report_error (t,t,"ok");; *)
73 let end_count = ref 0 in
75 let _ = end_count := !end_count + 1 in
76 (0 = ( !end_count mod 80000));;
78 let check_limit opt depth =
79 let _ = opt.iteration_count <- opt.iteration_count + 1 in
80 ( opt.iteration_count < opt.iteration_limit or opt.iteration_limit = 0 ) &&
81 (depth < opt.recursion_depth);;
83 let sgn x = if (x.lo > 0.0) then 1 else if (x.hi < 0.0) then -1 else 0;;
85 let rec same_sgn x y = (x = []) or (sgn (hd x) = sgn (hd y) && same_sgn (tl x) (tl y));;
88 (* a bit tricky because unstable exceptions are common early on,
89 and evaluations are very expensive.
90 We don't want a single unstable disjunct to ruin everything.
92 When boxes are small, then we throw away unstable disjuncts and continue w/o them.
93 When the boxes are still big, we return inconclusive to force a bisection.
95 Starting with this procedure, we can throw an exception to return early,
96 as soon as we are able to determine the cell_status. All these exceptions
97 get caught at the last line of verify_cell.
100 let rec set_targets (x,z,x0,z0,tf,tis,opt,maxwidth,has_unstable) =
103 let _ = not(has_unstable) or return (Cell_inconclusive (x,z,x0,z0,map snd tis)) in
106 let target = evalf (hd tf) x z in
108 let _ = not( opt.only_check_deriv1_negative) or return
109 (if upper_partial target 0 < 0.0 then Cell_pass
110 else if lower_partial target 0 > 0.0 then Cell_counterexample
111 else Cell_inconclusive(x,z,x0,z0,tf)) in
113 let _ = upper_bound target >= 0.0 or return Cell_pass in
115 set_targets(x,z,x0,z0,tl tf,(target,hd tf)::tis,opt,maxwidth,has_unstable));
118 if (2.0 *. maxwidth > opt.width_cutoff)
119 then set_targets(x,z,x0,z0,tl tf,tis,opt,maxwidth,true) (* proclaim unstable *)
120 else set_targets(x,z,x0,z0,tl tf,tis,opt,maxwidth,has_unstable) (* drop silently *);
124 let rec delete_false acc tis =
125 if (tis=[]) then List.rev acc
126 else if (lower_bound (fst (hd tis)) > 0.0) then delete_false acc (tl tis)
127 else delete_false (hd tis::acc) (tl tis);;
129 (* If the function is monotone, then we can push the calculation to the boundary.
130 There is a theorem that justifies the return Cell_pass steps. Here's a sketch.
131 We are doing a verification on a large box x0 z0 and
132 x z represents on of many cells in the box.
133 Suppose that there is a counterexample in the large box x0 z0.
134 The set of c/e is a compact set.
135 Write the inequalities as f1 < 0 \/ ... \/ fk < 0.
136 So at a c/e all f1 >= 0 && ... && fk >=0.
137 Pick out a particular counterexample by maximizing f1+...+fk over the
138 set of c/e's, then among this set pick a c/e with largest y1 coordinate,
139 then lexicographically continuing until the largest y6 coordinate.
140 This fixes a particular point y in the domain x0 z0.
142 I claim that has_monotone will detect y, so that if the inequality is false,
143 we find out about it.
145 In the partition of x0 z0 (of which x z is one rectangle), pick the rectangle
147 whose lower endpoint x1 is as large as possible, then continue lexicographically
150 I claim that has_monotone will detect y, when it comes to this box x z.
151 This is clear by construction. If we have monotonic increasing, then y is at the
152 right edge of the box. By construction the box x z is the rightmost containing y
153 so the box x z meets the right edge of the large containing box x0 z0.
154 Thus, it is enough to recursively search over the right edge of the box. Etc.
158 let rec has_monotone tis x z x0 z0 is found = match is with
159 | [] -> (x,z,x0,z0,found)
160 | j::js when (mth x j >= mth z j) ->
161 has_monotone tis x z x0 z0 js found
163 let allpos = List.fold_left (fun a ti -> a && lower_partial (fst ti) j >= 0.0) true tis in
164 let allneg = List.fold_left (fun a ti -> a && upper_partial (fst ti) j < 0.0) true tis in
166 let _ = (mth z j >= mth z0 j) or return Cell_pass in
167 let setj u = table (fun i -> (if i=j then mth z j else mth u i)) in
168 has_monotone tis (setj x) (setj z) (setj x0) (setj z0) js true
169 else if (allneg) then
170 let _ = (mth x j <= mth x0 j) or return Cell_pass in
171 let setj u = table (fun i -> (if i=j then mth x j else mth u i)) in
172 has_monotone tis (setj x) (setj z) (setj x0) (setj z0) js true
173 else has_monotone tis x z x0 z0 js found;;
175 (* loop as long as monotonicity keeps making progress. *)
177 let rec going_strong(x,z,x0,z0,tf,opt) =
178 let (y,w) = center_form (x,z) in
179 let maxwidth = maxl w in
180 let tis = set_targets(x,z,x0,z0,tf,[],opt,maxwidth,false) in
181 let epsilon_width = 1.0e-8 in
182 let _ = (maxwidth >= epsilon_width) or return
183 (if (opt.allow_sharp) then (Report.inc_corner_count(); Cell_pass)
184 else Cell_counterexample) in
185 let tis = delete_false [] tis in
186 let _ = (List.length tis > 0) or return Cell_counterexample in
187 let (x0,z0) = if (List.length tis < List.length tf) then (x,z) else (x0,z0) in
188 let (x,z,x0,z0,strong) =
189 if (opt.allow_derivatives) then has_monotone tis x z x0 z0 iter6 false
190 else (x,z,x0,z0,false) in
191 if (strong) then going_strong(x,z,x0,z0,map snd tis,opt) else (x,z,x0,z0,maxwidth,tis);;
193 let guess_optimal_corners (x,z,ti,tf) =
194 let mixedsign = List.fold_left (fun a i -> a or (sgn(mth ti.l.df i)=0)) false iter6 in
195 let yyn = table (fun i-> mth (if sgn (mth ti.l.df i) >0 then x else z) i) in
196 let yyu = table (fun i-> mth (if sgn (mth ti.l.df i) >0 then z else x) i) in
197 let cn = line_estimate tf yyn in
198 let cu = line_estimate tf yyu in
199 (mixedsign,yyn,yyu,cn,cu);;
201 let rec drop_numerically_false acc (x,z,x0,z0,tis) =
203 | [] -> (x0,z0,List.rev tis)
205 if (ti.l.f.lo <= 0.0) then drop_numerically_false ((ti,tf)::acc) (x,z,x0,z0,tiss)
207 let (mixedsign,yyn,yyu,cn,cu)= guess_optimal_corners(x,z,ti,tf) in
208 if (mixedsign) then drop_numerically_false ((ti,tf)::acc) (x,z,x0,z0,tiss)
209 else if (min cn.f.lo cu.f.lo > 0.0 &&
210 same_sgn ti.l.df cn.df && same_sgn ti.l.df cu.df)
211 then drop_numerically_false acc (x,z,x,z,tiss)
212 else drop_numerically_false ((ti,tf)::acc) (x,z,x0,z0,tiss);;
214 let rec keep_numerically_true best (x,z,x0,z0,tis,evalue) = match tis with
216 if (evalue < 0.0 && List.length best >0) then (x,z,best)
219 if (ti.l.f.hi > 0.0) then keep_numerically_true best (x,z,x0,z0,tiss,evalue)
221 let (mixedsign,yyn,yyu,cn,cu) = guess_optimal_corners(x,z,ti,tf) in
222 let evalue' = max cn.f.hi cu.f.hi in
223 if ((evalue' < evalue) &&
224 same_sgn ti.l.df cn.df && same_sgn ti.l.df cu.df && not(mixedsign)) then
225 keep_numerically_true [(ti,tf)] (x,z,x0,z0,tiss,evalue')
226 else keep_numerically_true best (x,z,x0,z0,tiss,evalue);;
229 This procedure is mostly guided by heuristics that don't require formal
230 verification. In particular, no justification is required for tossing out inequalities
231 (since they appear as disjuncts, we can choose which one to prove).
233 Formal verification is required whenever a Cell_passes is issued,
234 and whenever the domain (x,z) is restricted.
236 The record (x0,z0) of the current outer boundary must be restricted to (x,z)
237 whenever an inequality is tossed out.
240 let rec verify_cell (x,z,x0,z0,tf,opt) =
242 let _ = not(periodic_count ()) or report_current (x,z,"periodic report") in
243 let _ = not(opt.only_check_deriv1_negative) or (List.length tf <= 1) or
244 failwith "verify_cell: incompatible options" in
245 (* XX skip the implementation of rad2, delta126, delta135, for now. *)
246 let (x,z,x0,z0,maxwidth,tis) = going_strong(x,z,x0,z0,tf,opt) in
247 let (x0,z0,tis) = if (maxwidth < opt.width_cutoff && opt.allow_derivatives)
248 then drop_numerically_false [] (x,z,x0,z0,tis) else (x0,z0,tis) in
249 let _ = (List.length tis >0) or return Cell_counterexample in
251 if ((maxwidth<opt.width_cutoff)&&(List.length tis >1) && opt.allow_derivatives)
252 then keep_numerically_true [] (x,z,x0,z0,tis,0.0) else (x0,z0,tis) in
253 Cell_inconclusive (x,z,x0,z0,map snd tis);
257 let rec recursive_verifier (depth,x,z,x0,z0,tf,opt) =
258 let _ = check_limit opt depth or report_fatal(x,z,Printf.sprintf "depth %d" depth) in
259 match verify_cell(x,z,x0,z0,tf,opt) with
260 | Cell_counterexample -> false
262 | Cell_inconclusive(x,z,x0,z0,tf) ->
263 (let ( ++ ), ( / ) = up(); upadd, updiv in
264 let w2 = List.map2 upsub z x in
265 let maxwidth2 = maxl w2 in
266 let j_wide = try( find (fun i -> mth w2 i = maxwidth2) iter6) with
267 | _ -> failwith "recursive_verifier find" in
268 let yj = (mth x j_wide ++ mth z j_wide) / 2.0 in
269 let delta b v =table (fun i-> if (i=j_wide && b) then yj else mth v i) in
270 recursive_verifier(depth+1, delta false x ,delta true z ,x0,z0,tf,opt) &&
271 recursive_verifier(depth+1, delta true x ,delta false z ,x0,z0,tf,opt));;