port_interval/recurse.hl

   1 (* port of recurse.cc *)
   2
   3 (*
   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).
   7
   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.
  12
  13 *)
  14
  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";;
  21
  22 module Recurse = struct
  23
  24 open Interval;;
  25 open Univariate;;
  26 open Line_interval;;
  27 open Taylor;;
  28
  29 type cellOption = {
  30   only_check_deriv1_negative : bool;
  31   is_using_dihmax : bool;
  32   is_using_bigface126 : bool;
  33   width_cutoff : float;
  34   allow_sharp : bool;
  35   allow_derivatives : bool;
  36   mutable iteration_count : int;
  37   iteration_limit : int;
  38   recursion_depth : int;
  39 };;
  40
  41 (* cell verification is complex, and we use exceptions to
  42     exit as soon as the status has been determined.   *)
  43
  44 type  cell_status =
  45   | Cell_pass
  46   | Cell_counterexample
  47   | Cell_inconclusive of (float list * float list * float list * float list * tfunction list);;
  48
  49 exception Return of  cell_status;;
  50
  51 let return c = raise (Return c);;
  52
  53
  54 (* error checking and reporting functions *)
  55
  56 let string_of_domain x =
  57   let n = mth in
  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);;
  59
  60 let string3 (x,z,s) =  (string_of_domain x ^"\n"^ string_of_domain z ^ "\n" ^ s);;
  61
  62 let boolify _ = true;;
  63
  64 let report_current = boolify o Report.report_timed o string3;;
  65
  66 let report_error = boolify o Report.report_error o string3;;
  67
  68 let report_fatal = boolify o Report.report_fatal o string3;;
  69
  70 (* let t = [0.1;0.2;0.3;0.4;0.5;0.6] in report_error (t,t,"ok");; *)
  71
  72 let periodic_count =
  73   let end_count = ref 0 in
  74     fun () ->
  75       let _ = end_count := !end_count + 1 in
  76         (0 = ( !end_count mod 80000));;
  77
  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);;
  82
  83 let sgn x = if (x.lo > 0.0) then 1 else if (x.hi < 0.0) then -1 else 0;;
  84
  85 let rec same_sgn x y = (x = []) or (sgn (hd x) = sgn (hd y) && same_sgn (tl x) (tl y));;
  86
  87
  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.
  91
  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.
  94
  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.
  98 *)
  99
 100 let rec set_targets (x,z,x0,z0,tf,tis,opt,maxwidth,has_unstable) =
 101     try (
 102       if (tf = []) then
 103         let _ = not(has_unstable) or return (Cell_inconclusive (x,z,x0,z0,map snd tis)) in
 104           List.rev tis
 105       else (
 106         let target = evalf (hd tf)  x z in
 107
 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
 112
 113         let _ =  upper_bound target >= 0.0 or return Cell_pass in
 114
 115           set_targets(x,z,x0,z0,tl tf,(target,hd tf)::tis,opt,maxwidth,has_unstable));
 116     )
 117     with 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 *);
 121     );
 122 ;;
 123
 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);;
 128
 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.
 141
 142    I claim that has_monotone will detect y, so that if the inequality is false,
 143    we find out about it.
 144
 145    In the partition of x0 z0 (of which x z is one rectangle), pick the rectangle
 146    containing y
 147    whose lower endpoint x1 is as large as possible, then continue lexicographically
 148    through to x6.
 149
 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.
 155
 156  *)
 157
 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
 162   | j::js ->
 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
 165     if (allpos) then
 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;;
 174
 175 (* loop as long as monotonicity keeps making progress.  *)
 176
 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);;
 192
 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);;
 200
 201 let rec drop_numerically_false acc (x,z,x0,z0,tis) =
 202     match tis with
 203       | [] -> (x0,z0,List.rev tis)
 204       | (ti,tf) :: tiss ->
 205           if (ti.l.f.lo <= 0.0) then drop_numerically_false ((ti,tf)::acc) (x,z,x0,z0,tiss)
 206           else
 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);;
 213
 214 let rec keep_numerically_true best (x,z,x0,z0,tis,evalue) = match tis with
 215   | [] ->
 216        if (evalue < 0.0 && List.length best >0) then (x,z,best)
 217        else (x0,z0,tis)
 218   | (ti,tf)::tiss ->
 219       if (ti.l.f.hi > 0.0) then keep_numerically_true best (x,z,x0,z0,tiss,evalue)
 220       else
 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);;
 227
 228 (*
 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).
 232
 233 Formal verification is required whenever a Cell_passes is issued,
 234 and whenever the domain (x,z) is restricted.
 235
 236 The record (x0,z0) of the current outer boundary must be restricted to (x,z)
 237 whenever an inequality is tossed out.
 238 *)
 239
 240 let rec verify_cell (x,z,x0,z0,tf,opt) =
 241   try (
 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
 250   let (x0,z0,tis) =
 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);
 254   )
 255   with Return c -> c;;
 256
 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
 261       | Cell_pass -> true
 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));;
 272
 273
 274  end;;