formal_ineqs/verifier/interval_m/recurse.ml

   1 (* =========================================================== *)
   2 (* OCaml verification procedure                                *)
   3 (* Authors: Thomas C. Hales, Alexey Solovyev                   *)
   4 (* Date: 2012-10-27                                            *)
   5 (* =========================================================== *)
   6
   7 (* port of recurse.cc *)
   8
   9 (*
  10 This is the code that verifies a disjunct of nonlinear inequalities.
  11 The are given as a list   (tf:tfunction list).  If tf = [f1;....;fk], then
  12 the list represents the inequality (f1 < 0 \/ f2 < 0 .... fk < 0).
  13
  14 The end user should only need to define a cell option,
  15 and then call recursive_verifier, which recursively bisects the domain
  16 until a partition of the domain is found on which verifier_cell gives
  17 a pass on each piece of the partition.
  18
  19 *)
  20
  21 needs "verifier/interval_m/taylor.ml";;
  22 needs "verifier/interval_m/report.ml";;
  23 needs "verifier_options.hl";;
  24
  25 module Recurse = struct
  26
  27 open Interval_types;;
  28 open Interval;;
  29 open Univariate;;
  30 open Line_interval;;
  31 open Taylor;;
  32 open Verifier_options;;
  33 open List;;
  34
  35 type cellOption = {
  36   only_check_deriv1_negative : bool;
  37   is_using_dihmax : bool;
  38   is_using_bigface126 : bool;
  39   width_cutoff : float;
  40   allow_sharp : bool;
  41   allow_derivatives : bool;
  42   mutable iteration_count : int;
  43   iteration_limit : int;
  44   recursion_depth : int;
  45   mono_pass : bool;
  46   convex_flag : bool;
  47   raw_int_flag : bool;
  48   eps : float;
  49 };;
  50
  51 (* cell verification is complex, and we use exceptions to
  52     exit as soon as the status has been determined.   *)
  53
  54 type mono_status = {
  55   variable : int;
  56   decr_flag : bool;
  57   df0_flag : bool;
  58   ti_flag : bool;
  59 };;
  60
  61
  62 type  cell_status =
  63   | Cell_pass of mono_status list list * bool
  64   | Cell_pass_mono of mono_status list list * mono_status
  65   | Cell_counterexample
  66   | Cell_inconclusive_ti of (mono_status list list * taylor_interval * float list * float list * float list * float list)
  67   | Cell_inconclusive of (mono_status list list * float list * float list * float list * float list);;
  68
  69 exception Return of cell_status;;
  70
  71 type result_tree =
  72   | Result_false of (float list * float list)
  73   | Result_pass of (bool * float list * float list)
  74   | Result_pass_mono of mono_status
  75   | Result_pass_ref of int
  76   | Result_mono of mono_status list * result_tree
  77       (* variable, convex_flag, r1, r2 *)
  78   | Result_glue of (int * bool * result_tree * result_tree);;
  79
  80 type p_status = {
  81         pp : int;
  82 };;
  83
  84 type p_result_tree =
  85         | P_result_pass of p_status * bool
  86         | P_result_mono of p_status * mono_status list * p_result_tree
  87         | P_result_glue of p_status * int * bool * p_result_tree * p_result_tree
  88         | P_result_ref of int;;
  89
  90 let rec result_size r =
  91   match r with
  92     | Result_false _ -> failwith "False result detected"
  93     | Result_mono (_,r1) -> result_size r1
  94     | Result_glue (_, _, r1, r2) -> result_size r1 + result_size r2
  95     | Result_pass_mono _ -> 1
  96     | Result_pass _ -> 1
  97     | _ -> 0;;
  98
  99 let rec p_result_size r =
 100         match r with
 101                 | P_result_pass _ -> 1
 102                 | P_result_mono (_, _, r1) -> p_result_size r1
 103                 | P_result_glue (_, _, _, r1, r2) -> p_result_size r1 + p_result_size r2
 104                 | _ -> 0;;
 105
 106
 107 let return c = raise (Return c);;
 108
 109
 110 (* error checking and reporting functions *)
 111
 112 let string_of_domain x =
 113   let n = mth in
 114   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);;
 115
 116 let string3 (x,z,s) =  (string_of_domain x ^"\n"^ string_of_domain z ^ "\n" ^ s);;
 117
 118 let boolify _ = true;;
 119
 120 let report_current = boolify o Report.report_timed o string3;;
 121
 122 let report_error = boolify o Report.report_error o string3;;
 123
 124 let report_fatal = boolify o Report.report_fatal o string3;;
 125
 126 (* let t = [0.1;0.2;0.3;0.4;0.5;0.6] in report_error (t,t,"ok");; *)
 127
 128 let periodic_count =
 129   let end_count = ref 0 in
 130     fun () ->
 131       let _ = end_count := !end_count + 1 in
 132         (0 = ( !end_count mod 1000));;
 133
 134 let check_limit opt depth =
 135   let _ = opt.iteration_count <- opt.iteration_count + 1 in
 136    ( opt.iteration_count < opt.iteration_limit or opt.iteration_limit = 0 ) &&
 137      (depth < opt.recursion_depth);;
 138
 139 let sgn x = if (x.lo > 0.0) then 1 else if (x.hi < 0.0) then -1 else 0;;
 140
 141 let rec same_sgn x y = (x = []) or (sgn (hd x) = sgn (hd y) && same_sgn (tl x) (tl y));;
 142
 143
 144 (* has_monotone *)
 145
 146 let rec has_monotone opt tf ti domain0 x z x0 z0 is found = match is with
 147   | [] -> (x,z,x0,z0,List.rev found)
 148   | j::js when (mth x j >= mth z j) ->
 149       has_monotone opt tf ti domain0 x z x0 z0 js found
 150   | j::js ->
 151       let df_int =
 152         if opt.raw_int_flag then
 153           try evalf0 tf (j + 1) (fst domain0) (snd domain0)
 154           with Unstable -> mk_interval (-1.0,1.0)
 155         else
 156           mk_interval (-1.0, 1.0) in
 157       let allpos_df0, allpos_ti = df_int.lo >= opt.eps, lower_partial ti j >= opt.eps in
 158       let allneg_df0, allneg_ti = df_int.hi < ~-.(opt.eps), upper_partial ti j < ~-.(opt.eps) in
 159         if (allpos_df0 or allpos_ti) then
 160           let status =
 161             {variable = j + 1; decr_flag = false; df0_flag = allpos_df0; ti_flag = allpos_ti} in
 162             if opt.mono_pass && mth z j < mth z0 j then return (Cell_pass_mono ([], status))
 163             else
 164               let setj u = table (fun i -> (if i=j then mth z j else mth u i))  in
 165                 has_monotone opt tf ti domain0 (setj x) (setj z)
 166                   (setj x0) (setj z0) js (status :: found)
 167         else if (allneg_df0 or allneg_ti) then
 168           let status =
 169             {variable = j + 1; decr_flag = true; df0_flag = allneg_df0; ti_flag = allneg_ti} in
 170             if opt.mono_pass && mth x j > mth x0 j then return (Cell_pass_mono ([], status))
 171             else
 172               let setj u = table (fun i -> (if i=j then mth x j else mth u i)) in
 173                 has_monotone opt tf ti domain0 (setj x) (setj z)
 174                   (setj x0) (setj z0) js (status :: found)
 175         else has_monotone opt tf ti domain0 x z x0 z0 js found;;
 176
 177 (* loop as long as monotonicity keeps making progress.  *)
 178
 179 let rec going_strong(x,z,x0,z0,tf,opt,mono) =
 180   let (y,w) = center_form (x,z) in
 181   let maxwidth = maxl w in
 182   let target0 =
 183     if opt.raw_int_flag then
 184       try evalf0 tf 0 x z with Unstable -> one
 185     else
 186       one in
 187   let _ = target0.hi >= ~-.(opt.eps) or return (Cell_pass (mono, true)) in
 188   let target =
 189         try evalf tf x z with Unstable -> return (Cell_inconclusive (mono,x,z,x0,z0)) in
 190   let _ = upper_bound target >= ~-.(opt.eps) or return (Cell_pass (mono, false)) in
 191   let _ = lower_bound target < 0.0 or return Cell_counterexample in
 192   let epsilon_width = 1.0e-8 in
 193   let _ = (maxwidth >= epsilon_width) or return Cell_counterexample in
 194   let (x,z,x0,z0,strong) =
 195     if (opt.allow_derivatives) then
 196       try
 197         has_monotone opt tf target (x,z) x z x0 z0 iter8 []
 198       with Return (Cell_pass_mono (_, status)) -> return (Cell_pass_mono (mono, status))
 199     else (x,z,x0,z0,[]) in
 200     if (strong <> []) then
 201       going_strong(x,z,x0,z0,tf,opt,mono @ [strong])
 202     else
 203       (target,x,z,x0,z0,maxwidth,mono);;
 204
 205
 206 (*
 207 This procedure is mostly guided by heuristics that don't require formal
 208 verification. In particular, no justification is required for tossing out inequalities
 209 (since they appear as disjuncts, we can choose which one to prove).
 210
 211 Formal verification is required whenever a Cell_passes is issued,
 212 and whenever the domain (x,z) is restricted.
 213
 214 The record (x0,z0) of the current outer boundary must be restricted to (x,z)
 215 whenever an inequality is tossed out.
 216 *)
 217
 218 let rec verify_cell (x,z,x0,z0,tf,opt) =
 219   try (
 220   let _ = not(periodic_count () && !info_print_level >= 2) or report_current (x,z,"periodic report") in
 221   let (ti,x,z,x0,z0,maxwidth,mono) = going_strong(x,z,x0,z0,tf,opt,[]) in
 222     if opt.convex_flag then
 223       Cell_inconclusive_ti (mono,ti,x,z,x0,z0)
 224     else
 225       Cell_inconclusive (mono,x,z,x0,z0)
 226   )
 227   with Return c -> c;;
 228
 229 let recursive_verifier (x,z,x0,z0,tf,opt) =
 230   let w_init, indices = unzip (filter (fun p -> fst p > 1e-8) (zip (map2 (-.) z x) (1--length x))) in
 231   let ws = map2 (-.) z x in
 232   let total_vol = itlist ( *. ) w_init 1.0 in
 233   let verified_vol = ref 0.0 in
 234   let last_report = ref 0 in
 235   let compute_vol x z w =
 236     let rec compute i indices x z w =
 237       match indices with
 238         | [] -> 1.0
 239         | (r :: t) when r = i ->
 240             let l = hd z -. hd x in
 241               (if l > 1e-8 then l else hd w) *. compute (i + 1) t (tl x) (tl z) (tl w)
 242         | _ -> compute (i + 1) indices (tl x) (tl z) (tl w) in
 243       compute 1 indices x z w in
 244   let update_verified_vol x z w =
 245     if !info_print_level > 0 then
 246       let _ = verified_vol := !verified_vol +. compute_vol x z w in
 247       let verified = int_of_float (!verified_vol /. total_vol *. 100.5) in
 248       if verified > !last_report then
 249         let _ = last_report := verified in report0 (sprintf "%d " !last_report) else ()
 250     else () in
 251
 252   let rec rec_verifier (depth,x,z,x0,z0,w0,tf) =
 253     let _ = check_limit opt depth or report_fatal(x,z,Printf.sprintf "depth %d" depth) in
 254     let split_and_verify j x z x0 z0 convex_flag =
 255       let ( ++ ), ( / ) = up(); upadd, updiv in
 256       let yj = (mth x j ++  mth z j) / 2.0 in
 257       let delta b v = table (fun i-> if (i = j && b) then yj else mth v i) in
 258       let x1, z1 =
 259         if convex_flag then
 260           x, table (fun i -> if i = j then mth x i else mth z i)
 261         else
 262           delta false x, delta true z in
 263       let x2, z2 =
 264         if convex_flag then
 265           table (fun i -> if i = j then mth z i else mth x i), z
 266         else
 267           delta true x, delta false z in
 268       let w1 = table (fun i -> if i = j then mth w0 i / 2.0 else mth w0 i) in
 269       let r1 = rec_verifier(depth+1,x1,z1,x0,z0,w1,tf) in
 270         match r1 with
 271           | Result_false t -> Result_false t
 272           | _ ->
 273               (let r2 = rec_verifier(depth+1,x2,z2,x0,z0,w1,tf) in
 274                  match r2 with
 275                    | Result_false t -> Result_false t
 276                    | _ -> Result_glue (j, convex_flag, r1, r2)) in
 277
 278     let add_mono mono r1 =
 279       itlist (fun m r -> Result_mono (m, r)) mono r1 in
 280
 281       match verify_cell(x,z,x0,z0,tf,opt)  with
 282         | Cell_counterexample -> Result_false (x,z)
 283         | Cell_pass (mono, f0_flag) ->
 284             let _ = update_verified_vol x z w0 in
 285               add_mono mono (Result_pass (f0_flag,x,z))
 286         | Cell_pass_mono (mono, status) ->
 287             let _ = update_verified_vol x z w0 in
 288               add_mono mono (Result_pass_mono status)
 289         | Cell_inconclusive_ti(mono,ti,x,z,x0,z0) ->
 290             let dds = map (fun i -> mth (mth ti.dd i) i, i) iter8 in
 291             let convex_dds = filter (fun dd, i -> dd.lo >= opt.eps && mth x i < mth z i) dds in
 292             let convex_i = map snd convex_dds in
 293             let w2 = List.map2 upsub z x in
 294             let convex_flag, ws, ws_i =
 295               if convex_dds = [] then
 296                 false, w2, iter8
 297               else
 298                 true, map (mth w2) convex_i, convex_i in
 299             let maxwidth2 = maxl ws in
 300             let j_wide =  try( find (fun i -> mth w2 i = maxwidth2) ws_i) with
 301               | _ -> failwith "recursive_verifier find" in
 302               add_mono mono (split_and_verify j_wide x z x0 z0 convex_flag)
 303
 304         | Cell_inconclusive(mono,x,z,x0,z0) ->
 305             let w2 = List.map2 upsub z x in
 306             let maxwidth2 = maxl w2 in
 307             let j_wide =  try( find (fun i -> mth w2 i = maxwidth2) iter8) with
 308               | _ -> failwith "recursive_verifier find" in
 309               add_mono mono (split_and_verify j_wide x z x0 z0 false) in
 310
 311     rec_verifier (0,x,z,x0,z0,ws,tf);;
 312
 313
 314
 315  end;;