formal_lp/old/formal_interval/interval_m/recurse0.hl

   1 (* Recursive verification of inequalities using the basic interval arithmetic only *)
   2
   3 flyspeck_needs "../formal_lp/formal_interval/interval_m/recurse.hl";;
   4
   5 module Recurse0 = struct
   6
   7 open Interval;;
   8 open Univariate;;
   9 open Line_interval;;
  10 open Taylor;;
  11 open Recurse;;
  12
  13
  14
  15 (* has_monotone *)
  16
  17 let rec has_monotone0 opt tf domain0 x z x0 z0 is found = match is with
  18   | [] -> (x,z,x0,z0,List.rev found)
  19   | j::js when (mth x j >= mth z j) ->
  20       has_monotone0 opt tf domain0 x z x0 z0 js found
  21   | j::js ->
  22       let df_int = try evalf0 tf (j + 1) (fst domain0) (snd domain0) with Unstable -> mk_interval (-1.0,1.0) in
  23       let allpos_df0 = df_int.lo >= opt.eps in
  24       let allneg_df0 = df_int.hi < ~-.(opt.eps) in
  25         if allpos_df0 then
  26           let status =
  27             {variable = j + 1; decr_flag = false; df0_flag = allpos_df0; ti_flag = false} in
  28             if opt.mono_pass && mth z j < mth z0 j then return (Cell_pass_mono ([], status))
  29             else
  30               let setj u = table (fun i -> (if i=j then mth z j else mth u i))  in
  31                 has_monotone0 opt tf domain0 (setj x) (setj z)
  32                   (setj x0) (setj z0) js (status :: found)
  33         else if allneg_df0 then
  34           let status =
  35             {variable = j + 1; decr_flag = true; df0_flag = allneg_df0; ti_flag = false} in
  36             if opt.mono_pass && mth x j > mth x0 j then return (Cell_pass_mono ([], status))
  37             else
  38               let setj u = table (fun i -> (if i=j then mth x j else mth u i)) in
  39                 has_monotone0 opt tf domain0 (setj x) (setj z)
  40                   (setj x0) (setj z0) js (status :: found)
  41         else has_monotone0 opt tf domain0 x z x0 z0 js found;;
  42
  43 (* loop as long as monotonicity keeps making progress.  *)
  44
  45 let rec going_strong0(x,z,x0,z0,tf,opt,mono) =
  46   let (y,w) = center_form (x,z) in
  47   let maxwidth = maxl w in
  48   let target0 = try evalf0 tf 0 x z with Unstable -> return (Cell_inconclusive (mono,x,z,x0,z0)) in
  49   let _ = target0.hi >= ~-.(opt.eps) or return (Cell_pass (mono, true)) in
  50   let epsilon_width = 1.0e-8 in
  51   let _ = (maxwidth >= epsilon_width) or return Cell_counterexample in
  52   let (x,z,x0,z0,strong) =
  53     if (opt.allow_derivatives) then
  54       try
  55         has_monotone0 opt tf (x,z) x z x0 z0 iter8 []
  56       with Return (Cell_pass_mono (_, status)) -> return (Cell_pass_mono (mono, status))
  57     else (x,z,x0,z0,[]) in
  58     if (strong <> []) then
  59       going_strong0(x,z,x0,z0,tf,opt,mono @ [strong])
  60     else
  61       (x,z,x0,z0,maxwidth,mono);;
  62
  63
  64 (*
  65 This procedure is mostly guided by heuristics that don't require formal
  66 verification. In particular, no justification is required for tossing out inequalities
  67 (since they appear as disjuncts, we can choose which one to prove).
  68
  69 Formal verification is required whenever a Cell_passes is issued,
  70 and whenever the domain (x,z) is restricted.
  71
  72 The record (x0,z0) of the current outer boundary must be restricted to (x,z)
  73 whenever an inequality is tossed out.
  74 *)
  75
  76 let rec verify_cell0 (x,z,x0,z0,tf,opt) =
  77   try (
  78     let _ = not(periodic_count ()) or report_current (x,z,"periodic report") in
  79     let (x,z,x0,z0,maxwidth,mono) = going_strong0(x,z,x0,z0,tf,opt,[]) in
  80       if opt.convex_flag then
  81         let ti = try evalf tf x z with Unstable -> return (Cell_inconclusive (mono,x,z,x0,z0)) in
  82           Cell_inconclusive_ti (mono,ti,x,z,x0,z0)
  83       else
  84         Cell_inconclusive (mono,x,z,x0,z0)
  85   )
  86   with Return c -> c;;
  87
  88 let rec recursive_verifier0 (depth,x,z,x0,z0,tf,opt) =
  89   let _ = check_limit opt depth or report_fatal(x,z,Printf.sprintf "depth %d" depth) in
  90   let split_and_verify j x z x0 z0 convex_flag =
  91     let ( ++ ), ( / ) = up(); upadd, updiv in
  92     let yj = (mth x j ++  mth z j) / 2.0 in
  93     let delta b v = table (fun i-> if (i = j && b) then yj else mth v i) in
  94     let x1, z1 =
  95       if convex_flag then
  96         x, table (fun i -> if i = j then mth x i else mth z i)
  97       else
  98         delta false x, delta true z in
  99     let x2, z2 =
 100       if convex_flag then
 101         table (fun i -> if i = j then mth z i else mth x i), z
 102       else
 103         delta true x, delta false z in
 104     let r1 = recursive_verifier0(depth+1,x1,z1,x0,z0,tf,opt) in
 105       match r1 with
 106         | Result_false t -> Result_false t
 107         | _ ->
 108             (let r2 = recursive_verifier0(depth+1,x2,z2,x0,z0,tf,opt) in
 109                match r2 with
 110                  | Result_false t -> Result_false t
 111                  | _ -> Result_glue (j, convex_flag, r1, r2)) in
 112
 113   let add_mono mono r1 =
 114     itlist (fun m r -> Result_mono (m, r)) mono r1 in
 115
 116
 117     match verify_cell0(x,z,x0,z0,tf,opt)  with
 118       | Cell_counterexample -> Result_false (x,z)
 119       | Cell_pass (mono, f0_flag) -> add_mono mono (Result_pass (f0_flag,x,z))
 120       | Cell_pass_mono (mono, status) -> add_mono mono (Result_pass_mono status)
 121       | Cell_inconclusive_ti(mono,ti,x,z,x0,z0) ->
 122           let dds = map (fun i -> mth (mth ti.dd i) i, i) iter8 in
 123           let convex_dds = filter (fun dd, i -> dd.lo >= opt.eps && mth x i < mth z i) dds in
 124           let convex_i = map snd convex_dds in
 125           let w2 = List.map2 upsub z x in
 126           let convex_flag, ws, ws_i =
 127             if convex_dds = [] then
 128               false, w2, iter8
 129             else
 130               true, map (mth w2) convex_i, convex_i in
 131           let maxwidth2 = maxl ws in
 132           let j_wide =  try( find (fun i -> mth w2 i = maxwidth2) ws_i) with
 133             | _ -> failwith "recursive_verifier find" in
 134             add_mono mono (split_and_verify j_wide x z x0 z0 convex_flag)
 135
 136       | Cell_inconclusive(mono,x,z,x0,z0) ->
 137           let w2 = List.map2 upsub z x in
 138           let maxwidth2 = maxl w2 in
 139           let j_wide =  try( find (fun i -> mth w2 i = maxwidth2) iter8) with
 140             | _ -> failwith "recursive_verifier find" in
 141             add_mono mono (split_and_verify j_wide x z x0 z0 false);;
 142
 143
 144
 145  end;;