formal_ineqs/verifier/interval_m/recurse0.ml

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