formal_lp/old/formal_interval/interval_m/recurse.hl

   1
   2
   3
   4 (* port of recurse.cc *)
   5
   6 (*
   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).
  10
  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.
  15
  16 *)
  17
  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";;
  24
  25 module Recurse = struct
  26
  27 open Interval;;
  28 open Univariate;;
  29 open Line_interval;;
  30 open Taylor;;
  31
  32 type cellOption = {
  33   only_check_deriv1_negative : bool;
  34   is_using_dihmax : bool;
  35   is_using_bigface126 : bool;
  36   width_cutoff : float;
  37   allow_sharp : bool;
  38   allow_derivatives : bool;
  39   mutable iteration_count : int;
  40   iteration_limit : int;
  41   recursion_depth : int;
  42   mono_pass : bool;
  43   convex_flag : bool;
  44   raw_int_flag : bool;
  45   eps : float;
  46 };;
  47
  48 (* cell verification is complex, and we use exceptions to
  49     exit as soon as the status has been determined.   *)
  50
  51 type mono_status = {
  52   variable : int;
  53   decr_flag : bool;
  54   df0_flag : bool;
  55   ti_flag : bool;
  56 };;
  57
  58
  59 type  cell_status =
  60   | Cell_pass of mono_status list list * bool
  61   | Cell_pass_mono of mono_status list list * mono_status
  62   | Cell_counterexample
  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);;
  65
  66 exception Return of cell_status;;
  67
  68 type result_tree =
  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);;
  76
  77 type p_status = {
  78         pp : int;
  79 };;
  80
  81 type p_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;;
  86
  87 let rec result_size r =
  88   match r with
  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
  93     | Result_pass _ -> 1
  94     | _ -> 0;;
  95
  96 let rec p_result_size r =
  97         match r with
  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
 101                 | _ -> 0;;
 102
 103
 104 let return c = raise (Return c);;
 105
 106
 107 (* error checking and reporting functions *)
 108
 109 let string_of_domain x =
 110   let n = mth in
 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);;
 112
 113 let string3 (x,z,s) =  (string_of_domain x ^"\n"^ string_of_domain z ^ "\n" ^ s);;
 114
 115 let boolify _ = true;;
 116
 117 let report_current = boolify o Report.report_timed o string3;;
 118
 119 let report_error = boolify o Report.report_error o string3;;
 120
 121 let report_fatal = boolify o Report.report_fatal o string3;;
 122
 123 (* let t = [0.1;0.2;0.3;0.4;0.5;0.6] in report_error (t,t,"ok");; *)
 124
 125 let periodic_count =
 126   let end_count = ref 0 in
 127     fun () ->
 128       let _ = end_count := !end_count + 1 in
 129         (0 = ( !end_count mod 1000));;
 130
 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);;
 135
 136 let sgn x = if (x.lo > 0.0) then 1 else if (x.hi < 0.0) then -1 else 0;;
 137
 138 let rec same_sgn x y = (x = []) or (sgn (hd x) = sgn (hd y) && same_sgn (tl x) (tl y));;
 139
 140
 141 (* has_monotone *)
 142
 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
 147   | j::js ->
 148       let df_int =
 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)
 152         else
 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
 157           let status =
 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))
 160             else
 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
 165           let status =
 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))
 168             else
 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;;
 173
 174 (* loop as long as monotonicity keeps making progress.  *)
 175
 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
 179   let target0 =
 180     if opt.raw_int_flag then
 181       try evalf0 tf 0 x z with Unstable -> one
 182     else
 183       one in
 184   let _ = target0.hi >= ~-.(opt.eps) or return (Cell_pass (mono, true)) in
 185   let target =
 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
 193       try
 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])
 199     else
 200       (target,x,z,x0,z0,maxwidth,mono);;
 201
 202
 203 (*
 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).
 207
 208 Formal verification is required whenever a Cell_passes is issued,
 209 and whenever the domain (x,z) is restricted.
 210
 211 The record (x0,z0) of the current outer boundary must be restricted to (x,z)
 212 whenever an inequality is tossed out.
 213 *)
 214
 215 let rec verify_cell (x,z,x0,z0,tf,opt) =
 216   try (
 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)
 221     else
 222       Cell_inconclusive (mono,x,z,x0,z0)
 223   )
 224   with Return c -> c;;
 225
 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
 232     let x1, z1 =
 233       if convex_flag then
 234         x, table (fun i -> if i = j then mth x i else mth z i)
 235       else
 236         delta false x, delta true z in
 237     let x2, z2 =
 238       if convex_flag then
 239         table (fun i -> if i = j then mth z i else mth x i), z
 240       else
 241         delta true x, delta false z in
 242     let r1 = recursive_verifier(depth+1,x1,z1,x0,z0,tf,opt) in
 243       match r1 with
 244         | Result_false t -> Result_false t
 245         | _ ->
 246             (let r2 = recursive_verifier(depth+1,x2,z2,x0,z0,tf,opt) in
 247                match r2 with
 248                  | Result_false t -> Result_false t
 249                  | _ -> Result_glue (j, convex_flag, r1, r2)) in
 250
 251   let add_mono mono r1 =
 252     itlist (fun m r -> Result_mono (m, r)) mono r1 in
 253
 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
 265               false, w2, iter8
 266             else
 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)
 272
 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);;
 279
 280
 281
 282  end;;