formal_ineqs/verifier/interval_m/recurse.hl

   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 (* needs "verifier/interval_m/work1.hl";; *)
  36
  37 type cellOption = {
  38   only_check_deriv1_negative : bool;
  39   is_using_dihmax : bool;
  40   is_using_bigface126 : bool;
  41   width_cutoff : float;
  42   allow_sharp : bool;
  43   allow_derivatives : bool;
  44   mutable iteration_count : int;
  45   iteration_limit : int;
  46   recursion_depth : int;
  47   mono_pass : bool;
  48   convex_flag : bool;
  49   raw_int_flag : bool;
  50   eps : float;
  51 };;
  52
  53
  54 (* cell verification is complex, and we use exceptions to
  55     exit as soon as the status has been determined.   *)
  56
  57 type mono_status = {
  58   variable : int;
  59   decr_flag : bool;
  60   df0_flag : bool;
  61   ti_flag : bool;
  62 };;
  63
  64 type fun_type = {
  65   tf : tfunction;
  66   index : int;
  67 };;
  68
  69 (*
  70 let tfs0 = [{tf = tf1; index = 0;};
  71             {tf = tf2; index = 1;}];;
  72 *)
  73
  74 let all_indices = map (fun (_, f) -> f.index);;
  75
  76 type cell_status =
  77   | Cell_pass of mono_status list list * (int * bool)
  78   | Cell_pass_mono of mono_status list list * (int list * mono_status)
  79   | Cell_counterexample
  80   | Cell_inconclusive_ti of (taylor_interval * fun_type) list * (mono_status list list * float list * float list * float list * float list)
  81   | Cell_inconclusive of fun_type list * (mono_status list list * float list * float list * float list * float list);;
  82
  83 exception Return of cell_status;;
  84
  85 type result_tree =
  86   | Result_false of (float list * float list)
  87         (* function number, raw flag *)
  88   | Result_pass of (int * bool)
  89   | Result_pass_mono of mono_status
  90   | Result_pass_ref of int
  91   | Result_mono of mono_status list * result_tree
  92       (* variable, convex_flag, r1, r2 *)
  93   | Result_glue of (int * bool * result_tree * result_tree);;
  94
  95 type p_status = {
  96   pp : int;
  97 };;
  98
  99 type p_result_tree =
 100         | P_result_pass of p_status * int * bool
 101         | P_result_mono of p_status * mono_status list * p_result_tree
 102         | P_result_glue of p_status * int * bool * p_result_tree * p_result_tree
 103         | P_result_ref of int;;
 104
 105 let rec result_size r =
 106   match r with
 107     | Result_false _ -> failwith "False result detected"
 108     | Result_mono (_,r1) -> result_size r1
 109     | Result_glue (_, _, r1, r2) -> result_size r1 + result_size r2
 110     | Result_pass_mono _ -> 1
 111     | Result_pass _ -> 1
 112     | _ -> 0;;
 113
 114 let rec p_result_size r =
 115         match r with
 116                 | P_result_pass _ -> 1
 117                 | P_result_mono (_, _, r1) -> p_result_size r1
 118                 | P_result_glue (_, _, _, r1, r2) -> p_result_size r1 + p_result_size r2
 119                 | _ -> 0;;
 120
 121
 122 let return c = raise (Return c);;
 123
 124
 125 (* error checking and reporting functions *)
 126
 127 let string_of_domain x =
 128   let n = mth in
 129   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);;
 130
 131 let string3 (x,z,s) =  (string_of_domain x ^"\n"^ string_of_domain z ^ "\n" ^ s);;
 132
 133 let boolify _ = true;;
 134
 135 let report_current = boolify o Report.report_timed o string3;;
 136
 137 let report_error = boolify o Report.report_error o string3;;
 138
 139 let report_fatal = boolify o Report.report_fatal o string3;;
 140
 141 (* let t = [0.1;0.2;0.3;0.4;0.5;0.6] in report_error (t,t,"ok");; *)
 142
 143 let periodic_count =
 144   let end_count = ref 0 in
 145     fun () ->
 146       let _ = end_count := !end_count + 1 in
 147         (0 = ( !end_count mod 1000));;
 148
 149 let check_limit opt depth =
 150   let _ = opt.iteration_count <- opt.iteration_count + 1 in
 151    ( opt.iteration_count < opt.iteration_limit or opt.iteration_limit = 0 ) &&
 152      (depth < opt.recursion_depth);;
 153
 154 let sgn x = if (x.lo > 0.0) then 1 else if (x.hi < 0.0) then -1 else 0;;
 155
 156 let rec same_sgn x y = (x = []) or (sgn (hd x) = sgn (hd y) && same_sgn (tl x) (tl y));;
 157
 158
 159 let dummy_taylor = make_taylor_interval ((mk_line (one, [])), [], []);;
 160
 161 (* computes taylor intervals for all functions *)
 162 let rec set_targets (x,z,x0,z0,tfs,tis,opt,mono,maxwidth,has_unstable) =
 163   if (tfs = []) then
 164     let _ = not(has_unstable) or return (Cell_inconclusive (map snd tis, (mono,x,z,x0,z0))) in
 165       List.rev tis
 166   else
 167     let tf = (hd tfs).tf and
 168         index = (hd tfs).index in
 169     let target0 =
 170       if opt.raw_int_flag then
 171         try evalf0 tf 0 x z with Unstable -> one
 172       else
 173         one in
 174     let _ = target0.hi >= ~-.(opt.eps) or return (Cell_pass (mono, (index, true))) in
 175       try
 176         let target = evalf tf x z in
 177         let _ = upper_bound target >= ~-.(opt.eps) or return (Cell_pass (mono, (index, false))) in
 178           set_targets(x,z,x0,z0,tl tfs,(target,hd tfs)::tis,opt,mono,maxwidth,has_unstable)
 179       with Unstable ->
 180         if (2.0 *. maxwidth > opt.width_cutoff) then
 181           set_targets(x,z,x0,z0,tl tfs,(dummy_taylor, hd tfs)::tis,opt,mono,maxwidth,true)  (* proclaim unstable *)
 182         else
 183           set_targets(x,z,x0,z0,tl tfs,tis,opt,mono,maxwidth,has_unstable) (* drop silently *)
 184 ;;
 185
 186 (*
 187 let tis = set_targets (xx, zz, xx, zz, tfs0, [], opt0, [], 1.0, false);;
 188 let ti = fst (nth tis 0);;
 189 lower_bound ti;;
 190 *)
 191
 192 let rec delete_false acc tis =
 193   if (tis=[]) then List.rev acc
 194   else if (lower_bound (fst (hd tis)) > 0.0) then delete_false acc (tl tis)
 195   else delete_false (hd tis::acc) (tl tis);;
 196
 197
 198 (* has_monotone *)
 199 let rec has_monotone opt tis domain0 x z x0 z0 is found = match is with
 200   | [] -> (x,z,x0,z0,List.rev found)
 201   | j::js when (mth x j >= mth z j) ->
 202       has_monotone opt tis domain0 x z x0 z0 js found
 203   | j::js ->
 204       let df_ints =
 205         if opt.raw_int_flag then
 206           map (fun (ti, f) ->
 207                  (try
 208                     evalf0 f.tf (j + 1) (fst domain0) (snd domain0)
 209                   with Unstable -> mk_interval (-1.0,1.0)))
 210             tis
 211         else
 212           [mk_interval (-1.0, 1.0)] in
 213       let allpos_df0 = List.fold_left (fun a i -> a && i.lo >= opt.eps) true df_ints and
 214           allneg_df0 = List.fold_left (fun a i -> a && i.hi < ~-.(opt.eps)) true df_ints in
 215       let allpos_ti = List.fold_left (fun a ti -> a && lower_partial (fst ti) j >= opt.eps) true tis and
 216           allneg_ti = List.fold_left (fun a ti -> a && upper_partial (fst ti) j < ~-.(opt.eps)) true tis in
 217
 218         if (allpos_df0 or allpos_ti) then
 219           let status =
 220             {variable = j + 1; decr_flag = false; df0_flag = allpos_df0; ti_flag = allpos_ti} in
 221             if opt.mono_pass && mth z j < mth z0 j then
 222               return (Cell_pass_mono ([], (all_indices tis, status)))
 223             else
 224               let setj u = table (fun i -> (if i=j then mth z j else mth u i)) in
 225                 has_monotone opt tis domain0 (setj x) (setj z)
 226                   (setj x0) (setj z0) js (status :: found)
 227         else if (allneg_df0 or allneg_ti) then
 228           let status =
 229             {variable = j + 1; decr_flag = true; df0_flag = allneg_df0; ti_flag = allneg_ti} in
 230             if opt.mono_pass && mth x j > mth x0 j then
 231               return (Cell_pass_mono ([], (all_indices tis, status)))
 232             else
 233               let setj u = table (fun i -> (if i=j then mth x j else mth u i)) in
 234                 has_monotone opt tis domain0 (setj x) (setj z)
 235                   (setj x0) (setj z0) js (status :: found)
 236         else has_monotone opt tis domain0 x z x0 z0 js found;;
 237
 238 (* loop as long as monotonicity keeps making progress.  *)
 239
 240 let rec going_strong(x,z,x0,z0,tfs,opt,mono) =
 241   let (y,w) = center_form (x,z) in
 242   let maxwidth = maxl w in
 243   let tis = set_targets (x,z,x0,z0,tfs,[],opt,mono,maxwidth,false) in
 244   let epsilon_width = 1.0e-8 in
 245   let _ = (maxwidth >= epsilon_width) or return Cell_counterexample in
 246   let tis = delete_false [] tis in
 247   let _ = (List.length tis > 0) or return Cell_counterexample in
 248   let (x0,z0) = if (length tis < length tfs) then (x,z) else (x0,z0) in
 249   let (x,z,x0,z0,strong) =
 250     if (opt.allow_derivatives) then
 251       try
 252         has_monotone opt tis (x,z) x z x0 z0 iter8 []
 253       with Return (Cell_pass_mono (_, i_status)) ->
 254         return (Cell_pass_mono (mono, i_status))
 255     else (x,z,x0,z0,[]) in
 256     if (strong <> []) then
 257       going_strong(x,z,x0,z0,map snd tis,opt,mono @ [strong])
 258     else
 259       (tis,x,z,x0,z0,maxwidth,mono);;
 260
 261 (*
 262 let xx = [0.0; -1.0] @ pad;;
 263 let zz = [1.0; 0.0] @ pad;;
 264 let _, x,z,x0,z0,_,mono = going_strong (xx,zz,xx,zz,tfs0,opt0,[]);;
 265
 266 let tis = set_targets (xx,zz,xx,zz,tfs0,[],opt0,[],1.0,false);;
 267
 268 let j = 0 and
 269     opt = opt0 and
 270     domain0 = (xx,zz);;
 271
 272 evalf0 tf1 1 xx zz;;
 273
 274 let df_ints =
 275   if opt.raw_int_flag then
 276     map (fun (ti, f) ->
 277            (try
 278               evalf0 f.tf (j + 1) (fst domain0) (snd domain0)
 279             with Unstable -> mk_interval (-1.0,1.0)))
 280       tis
 281   else
 282     [mk_interval (-1.0, 1.0)];;
 283       let allpos_df0 = List.fold_left (fun a i -> a && i.lo >= opt.eps) true df_ints and
 284           allneg_df0 = List.fold_left (fun a i -> a && i.hi < ~-.(opt.eps)) true df_ints in
 285       let allpos_ti = List.fold_left (fun a ti -> a && lower_partial (fst ti) j >= opt.eps) true tis and
 286           allneg_ti = List.fold_left (fun a ti -> a && upper_partial (fst ti) j < ~-.(opt.eps)) true tis in
 287 *)
 288
 289
 290 (*
 291 This procedure is mostly guided by heuristics that don't require formal
 292 verification. In particular, no justification is required for tossing out inequalities
 293 (since they appear as disjuncts, we can choose which one to prove).
 294
 295 Formal verification is required whenever a Cell_passes is issued,
 296 and whenever the domain (x,z) is restricted.
 297
 298 The record (x0,z0) of the current outer boundary must be restricted to (x,z)
 299 whenever an inequality is tossed out.
 300 *)
 301
 302 let rec verify_cell (x,z,x0,z0,tfs,opt) =
 303   try (
 304   let _ = not(periodic_count () && !info_print_level >= 2) or report_current (x,z,"periodic report") in
 305   let (tis,x,z,x0,z0,maxwidth,mono) = going_strong(x,z,x0,z0,tfs,opt,[]) in
 306     if opt.convex_flag then
 307       Cell_inconclusive_ti (tis, (mono,x,z,x0,z0))
 308     else
 309       Cell_inconclusive (map snd tis, (mono,x,z,x0,z0))
 310   )
 311   with Return c -> c;;
 312
 313 let recursive_verifier (x,z,x0,z0,tfs,opt) =
 314   let w_init, indices = unzip (filter (fun p -> fst p > 1e-8) (zip (map2 (-.) z x) (1--length x))) in
 315   let ws = map2 (-.) z x in
 316   let total_vol = itlist ( *. ) w_init 1.0 in
 317   let verified_vol = ref 0.0 in
 318   let last_report = ref 0 in
 319   let compute_vol x z w =
 320     let rec compute i indices x z w =
 321       match indices with
 322         | [] -> 1.0
 323         | (r :: t) when r = i ->
 324             let l = hd z -. hd x in
 325               (if l > 1e-8 then l else hd w) *. compute (i + 1) t (tl x) (tl z) (tl w)
 326         | _ -> compute (i + 1) indices (tl x) (tl z) (tl w) in
 327       compute 1 indices x z w in
 328   let update_verified_vol x z w =
 329     if !info_print_level > 0 then
 330       let _ = verified_vol := !verified_vol +. compute_vol x z w in
 331       let verified = int_of_float (!verified_vol /. total_vol *. 100.5) in
 332       if verified > !last_report then
 333         let _ = last_report := verified in report0 (sprintf "%d " !last_report) else ()
 334     else () in
 335
 336   let add_mono mono r1 =
 337     itlist (fun m r -> Result_mono (m, r)) mono r1 in
 338
 339   let rec rec_verifier (depth,x,z,x0,z0,w0,tfs) =
 340     let split_and_verify tfs j x z x0 z0 convex_flag =
 341       let ( ++ ), ( / ) = up(); upadd, updiv in
 342       let yj = (mth x j ++  mth z j) / 2.0 in
 343       let delta b v = table (fun i-> if (i = j && b) then yj else mth v i) in
 344       let x1, z1 =
 345         if convex_flag then
 346           x, table (fun i -> if i = j then mth x i else mth z i)
 347         else
 348           delta false x, delta true z in
 349       let x2, z2 =
 350         if convex_flag then
 351           table (fun i -> if i = j then mth z i else mth x i), z
 352         else
 353           delta true x, delta false z in
 354       let w1 = table (fun i -> if i = j then mth w0 i / 2.0 else mth w0 i) in
 355       let r1 = rec_verifier(depth+1,x1,z1,x0,z0,w1,tfs) in
 356         match r1 with
 357           | Result_false t -> Result_false t
 358           | _ ->
 359               (let r2 = rec_verifier(depth+1,x2,z2,x0,z0,w1,tfs) in
 360                  match r2 with
 361                    | Result_false t -> Result_false t
 362                    | _ -> Result_glue (j, convex_flag, r1, r2)) in
 363
 364     let _ = check_limit opt depth or report_fatal(x,z,Printf.sprintf "depth %d" depth) in
 365       match verify_cell(x,z,x0,z0,tfs,opt)  with
 366         | Cell_counterexample -> Result_false (x,z)
 367         | Cell_pass (mono, (i, f0_flag)) ->
 368             let _ = update_verified_vol x z w0 in
 369               add_mono mono (Result_pass (i, f0_flag))
 370         | Cell_pass_mono (mono, (is, status)) ->
 371             let _ = update_verified_vol x z w0 in
 372               add_mono mono (Result_pass_mono status)
 373 (*
 374 (* TODO: convexity *)
 375   | Cell_inconclusive_ti (tis, (mono,x,z,x0,z0)) ->
 376             let dds = map (fun i -> mth (mth ti.dd i) i, i) iter8 in
 377             let convex_dds = filter (fun dd, i -> dd.lo >= opt.eps && mth x i < mth z i) dds in
 378             let convex_i = map snd convex_dds in
 379             let w2 = List.map2 upsub z x in
 380             let convex_flag, ws, ws_i =
 381               if convex_dds = [] then
 382                 false, w2, iter8
 383               else
 384                 true, map (mth w2) convex_i, convex_i in
 385             let maxwidth2 = maxl ws in
 386             let j_wide =  try( find (fun i -> mth w2 i = maxwidth2) ws_i) with
 387               | _ -> failwith "recursive_verifier find" in
 388               add_mono mono (split_and_verify tfs j_wide x z x0 z0 convex_flag)
 389 *)
 390         | Cell_inconclusive_ti _ -> failwith "Convexity is not supported"
 391         | Cell_inconclusive (tfs, (mono,x,z,x0,z0)) ->
 392             let w2 = List.map2 upsub z x in
 393             let maxwidth2 = maxl w2 in
 394             let j_wide =  try( find (fun i -> mth w2 i = maxwidth2) iter8) with
 395               | _ -> failwith "recursive_verifier find" in
 396               add_mono mono (split_and_verify tfs j_wide x z x0 z0 false) in
 397
 398     rec_verifier (0,x,z,x0,z0,ws,tfs);;
 399
 400 (*
 401 let _, tf3, _ = mk_verification_functions_poly 3 `\x:real^2. --(&1 - x$1 * x$1 - x$2 * x$2)`;;
 402 let tfs1 = [{tf = tf3; index = 0;}];;
 403 let xx = [-0.75; -0.7] @ pad and
 404     zz = [0.75; 0.7] @ pad;;
 405
 406
 407 verify_cell (xx,zz,xx,zz,tfs1,opt0);;
 408 recursive_verifier (xx,zz,xx,zz,tfs1,opt0);;
 409 *)
 410
 411 end;;