4 (* port of recurse.cc *)
7 This is the code that verifies a disjunct of nonlinear inequalities.
8 The are given as a list (tf:tfunction list). If tf = [f1;....;fk], then
9 the list represents the inequality (f1 < 0 \/ f2 < 0 .... fk < 0).
11 The end user should only need to define a cell option,
12 and then call recursive_verifier, which recursively bisects the domain
13 until a partition of the domain is found on which verifier_cell gives
14 a pass on each piece of the partition.
18 flyspeck_needs "../formal_lp/formal_interval/interval_m/types.hl";;
19 flyspeck_needs "../formal_lp/formal_interval/interval_m/report.hl";;
20 flyspeck_needs "../formal_lp/formal_interval/interval_m/interval.hl";;
21 flyspeck_needs "../formal_lp/formal_interval/interval_m/univariate.hl";;
22 flyspeck_needs "../formal_lp/formal_interval/interval_m/line_interval.hl";;
23 flyspeck_needs "../formal_lp/formal_interval/interval_m/taylor.hl";;
25 module Recurse = struct
33 only_check_deriv1_negative : bool;
34 is_using_dihmax : bool;
35 is_using_bigface126 : bool;
38 allow_derivatives : bool;
39 mutable iteration_count : int;
40 iteration_limit : int;
41 recursion_depth : int;
48 (* cell verification is complex, and we use exceptions to
49 exit as soon as the status has been determined. *)
60 | Cell_pass of mono_status list list * bool
61 | Cell_pass_mono of mono_status list list * mono_status
63 | Cell_inconclusive_ti of (mono_status list list * taylor_interval * float list * float list * float list * float list)
64 | Cell_inconclusive of (mono_status list list * float list * float list * float list * float list);;
66 exception Return of cell_status;;
69 | Result_false of (float list * float list)
70 | Result_pass of (bool * float list * float list)
71 | Result_pass_mono of mono_status
72 | Result_pass_ref of int
73 | Result_mono of mono_status list * result_tree
74 (* variable, convex_flag, r1, r2 *)
75 | Result_glue of (int * bool * result_tree * result_tree);;
82 | P_result_pass of p_status * bool
83 | P_result_mono of p_status * mono_status list * p_result_tree
84 | P_result_glue of p_status * int * bool * p_result_tree * p_result_tree
85 | P_result_ref of int;;
87 let rec result_size r =
89 | Result_false _ -> failwith "False result detected"
90 | Result_mono (_,r1) -> result_size r1
91 | Result_glue (_, _, r1, r2) -> result_size r1 + result_size r2
92 | Result_pass_mono _ -> 1
96 let rec p_result_size r =
98 | P_result_pass _ -> 1
99 | P_result_mono (_, _, r1) -> p_result_size r1
100 | P_result_glue (_, _, _, r1, r2) -> p_result_size r1 + p_result_size r2
104 let return c = raise (Return c);;
107 (* error checking and reporting functions *)
109 let string_of_domain x =
111 Printf.sprintf "{%f, %f, %f, %f, %f, %f, %f, %f}" (n x 0) (n x 1) (n x 2) (n x 3) (n x 4) (n x 5) (n x 6) (n x 7);;
113 let string3 (x,z,s) = (string_of_domain x ^"\n"^ string_of_domain z ^ "\n" ^ s);;
115 let boolify _ = true;;
117 let report_current = boolify o Report.report_timed o string3;;
119 let report_error = boolify o Report.report_error o string3;;
121 let report_fatal = boolify o Report.report_fatal o string3;;
123 (* let t = [0.1;0.2;0.3;0.4;0.5;0.6] in report_error (t,t,"ok");; *)
126 let end_count = ref 0 in
128 let _ = end_count := !end_count + 1 in
129 (0 = ( !end_count mod 1000));;
131 let check_limit opt depth =
132 let _ = opt.iteration_count <- opt.iteration_count + 1 in
133 ( opt.iteration_count < opt.iteration_limit or opt.iteration_limit = 0 ) &&
134 (depth < opt.recursion_depth);;
136 let sgn x = if (x.lo > 0.0) then 1 else if (x.hi < 0.0) then -1 else 0;;
138 let rec same_sgn x y = (x = []) or (sgn (hd x) = sgn (hd y) && same_sgn (tl x) (tl y));;
143 let rec has_monotone opt tf ti domain0 x z x0 z0 is found = match is with
144 | [] -> (x,z,x0,z0,List.rev found)
145 | j::js when (mth x j >= mth z j) ->
146 has_monotone opt tf ti domain0 x z x0 z0 js found
149 if opt.raw_int_flag then
150 try evalf0 tf (j + 1) (fst domain0) (snd domain0)
151 with Unstable -> mk_interval (-1.0,1.0)
153 mk_interval (-1.0, 1.0) in
154 let allpos_df0, allpos_ti = df_int.lo >= opt.eps, lower_partial ti j >= opt.eps in
155 let allneg_df0, allneg_ti = df_int.hi < ~-.(opt.eps), upper_partial ti j < ~-.(opt.eps) in
156 if (allpos_df0 or allpos_ti) then
158 {variable = j + 1; decr_flag = false; df0_flag = allpos_df0; ti_flag = allpos_ti} in
159 if opt.mono_pass && mth z j < mth z0 j then return (Cell_pass_mono ([], status))
161 let setj u = table (fun i -> (if i=j then mth z j else mth u i)) in
162 has_monotone opt tf ti domain0 (setj x) (setj z)
163 (setj x0) (setj z0) js (status :: found)
164 else if (allneg_df0 or allneg_ti) then
166 {variable = j + 1; decr_flag = true; df0_flag = allneg_df0; ti_flag = allneg_ti} in
167 if opt.mono_pass && mth x j > mth x0 j then return (Cell_pass_mono ([], status))
169 let setj u = table (fun i -> (if i=j then mth x j else mth u i)) in
170 has_monotone opt tf ti domain0 (setj x) (setj z)
171 (setj x0) (setj z0) js (status :: found)
172 else has_monotone opt tf ti domain0 x z x0 z0 js found;;
174 (* loop as long as monotonicity keeps making progress. *)
176 let rec going_strong(x,z,x0,z0,tf,opt,mono) =
177 let (y,w) = center_form (x,z) in
178 let maxwidth = maxl w in
180 if opt.raw_int_flag then
181 try evalf0 tf 0 x z with Unstable -> one
184 let _ = target0.hi >= ~-.(opt.eps) or return (Cell_pass (mono, true)) in
186 try evalf tf x z with Unstable -> return (Cell_inconclusive (mono,x,z,x0,z0)) in
187 let _ = upper_bound target >= ~-.(opt.eps) or return (Cell_pass (mono, false)) in
188 let _ = lower_bound target < 0.0 or return Cell_counterexample in
189 let epsilon_width = 1.0e-8 in
190 let _ = (maxwidth >= epsilon_width) or return Cell_counterexample in
191 let (x,z,x0,z0,strong) =
192 if (opt.allow_derivatives) then
194 has_monotone opt tf target (x,z) x z x0 z0 iter8 []
195 with Return (Cell_pass_mono (_, status)) -> return (Cell_pass_mono (mono, status))
196 else (x,z,x0,z0,[]) in
197 if (strong <> []) then
198 going_strong(x,z,x0,z0,tf,opt,mono @ [strong])
200 (target,x,z,x0,z0,maxwidth,mono);;
204 This procedure is mostly guided by heuristics that don't require formal
205 verification. In particular, no justification is required for tossing out inequalities
206 (since they appear as disjuncts, we can choose which one to prove).
208 Formal verification is required whenever a Cell_passes is issued,
209 and whenever the domain (x,z) is restricted.
211 The record (x0,z0) of the current outer boundary must be restricted to (x,z)
212 whenever an inequality is tossed out.
215 let rec verify_cell (x,z,x0,z0,tf,opt) =
217 let _ = not(periodic_count ()) or report_current (x,z,"periodic report") in
218 let (ti,x,z,x0,z0,maxwidth,mono) = going_strong(x,z,x0,z0,tf,opt,[]) in
219 if opt.convex_flag then
220 Cell_inconclusive_ti (mono,ti,x,z,x0,z0)
222 Cell_inconclusive (mono,x,z,x0,z0)
226 let rec recursive_verifier (depth,x,z,x0,z0,tf,opt) =
227 let _ = check_limit opt depth or report_fatal(x,z,Printf.sprintf "depth %d" depth) in
228 let split_and_verify j x z x0 z0 convex_flag =
229 let ( ++ ), ( / ) = up(); upadd, updiv in
230 let yj = (mth x j ++ mth z j) / 2.0 in
231 let delta b v = table (fun i-> if (i = j && b) then yj else mth v i) in
234 x, table (fun i -> if i = j then mth x i else mth z i)
236 delta false x, delta true z in
239 table (fun i -> if i = j then mth z i else mth x i), z
241 delta true x, delta false z in
242 let r1 = recursive_verifier(depth+1,x1,z1,x0,z0,tf,opt) in
244 | Result_false t -> Result_false t
246 (let r2 = recursive_verifier(depth+1,x2,z2,x0,z0,tf,opt) in
248 | Result_false t -> Result_false t
249 | _ -> Result_glue (j, convex_flag, r1, r2)) in
251 let add_mono mono r1 =
252 itlist (fun m r -> Result_mono (m, r)) mono r1 in
254 match verify_cell(x,z,x0,z0,tf,opt) with
255 | Cell_counterexample -> Result_false (x,z)
256 | Cell_pass (mono, f0_flag) -> add_mono mono (Result_pass (f0_flag,x,z))
257 | Cell_pass_mono (mono, status) -> add_mono mono (Result_pass_mono status)
258 | Cell_inconclusive_ti(mono,ti,x,z,x0,z0) ->
259 let dds = map (fun i -> mth (mth ti.dd i) i, i) iter8 in
260 let convex_dds = filter (fun dd, i -> dd.lo >= opt.eps && mth x i < mth z i) dds in
261 let convex_i = map snd convex_dds in
262 let w2 = List.map2 upsub z x in
263 let convex_flag, ws, ws_i =
264 if convex_dds = [] then
267 true, map (mth w2) convex_i, convex_i in
268 let maxwidth2 = maxl ws in
269 let j_wide = try( find (fun i -> mth w2 i = maxwidth2) ws_i) with
270 | _ -> failwith "recursive_verifier find" in
271 add_mono mono (split_and_verify j_wide x z x0 z0 convex_flag)
273 | Cell_inconclusive(mono,x,z,x0,z0) ->
274 let w2 = List.map2 upsub z x in
275 let maxwidth2 = maxl w2 in
276 let j_wide = try( find (fun i -> mth w2 i = maxwidth2) iter8) with
277 | _ -> failwith "recursive_verifier find" in
278 add_mono mono (split_and_verify j_wide x z x0 z0 false);;