formal_ineqs/verifier/interval_m/verifier.ml

   1 (* =========================================================== *)
   2 (* OCaml verification and result transformation functions      *)
   3 (* Author: Alexey Solovyev                                     *)
   4 (* Date: 2012-10-27                                            *)
   5 (* =========================================================== *)
   6
   7 needs "verifier/interval_m/recurse.ml";;
   8 needs "verifier/interval_m/recurse0.ml";;
   9
  10 module Verifier = struct
  11
  12 open Interval_types;;
  13 open Interval;;
  14 open Univariate;;
  15 open Line_interval;;
  16 open Taylor;;
  17 open Recurse;;
  18
  19
  20 type certificate_stats =
  21 {
  22   pass : int;
  23   pass_raw : int;
  24   pass_mono : int;
  25   mono : int;
  26   glue : int;
  27   glue_convex : int;
  28 };;
  29
  30
  31 let dummy_stats =
  32 {
  33   pass = 0; pass_raw = 0; pass_mono = 0;
  34   mono = 0; glue = 0; glue_convex = 0;
  35 };;
  36
  37
  38 (**********************************)
  39 let run_test f x z min_flag min_max allow_d convex_flag mono_pass_flag raw_int_flag eps =
  40   let pad = replicate 0.0 (8 - length x) in
  41   let xx = x @ pad and zz = z @ pad in
  42   let mone = mk_interval(-1.0,-1.0) in
  43   let neg_f = Scale(f, mone) in
  44   let ff = if min_flag then
  45     Plus(neg_f, Scale(unit,mk_interval(min_max, min_max)))
  46   else
  47     Plus(f, Scale(unit, ineg (mk_interval(min_max, min_max)))) in
  48   let opt =  {
  49     only_check_deriv1_negative = false;
  50     is_using_dihmax =false;
  51     is_using_bigface126 =false;
  52     width_cutoff =0.05;
  53     allow_sharp =false;
  54     allow_derivatives =allow_d;
  55     iteration_count =0;
  56     iteration_limit =0;
  57     recursion_depth =200;
  58     mono_pass = mono_pass_flag;
  59     convex_flag = convex_flag;
  60     raw_int_flag = raw_int_flag;
  61     eps = eps;
  62   } in
  63     recursive_verifier(xx,zz,xx,zz,ff,opt);;
  64
  65
  66 (* A verification procedure which uses raw interval arithmetic only *)
  67 (*
  68 open Recurse0;;
  69
  70 let run_test0 f x z min_flag min_max allow_d convex_flag mono_pass_flag eps =
  71   let pad = replicate 0.0 (8 - length x) in
  72   let xx = x @ pad and zz = z @ pad in
  73   let mone = mk_interval(-1.0,-1.0) in
  74   let neg_f = Scale(f, mone) in
  75   let ff = if min_flag then
  76     Plus(neg_f, Scale(unit,mk_interval(min_max, min_max)))
  77   else
  78     Plus(f, Scale(unit, ineg (mk_interval(min_max, min_max)))) in
  79   let opt =  {
  80     only_check_deriv1_negative = false;
  81     is_using_dihmax =false;
  82     is_using_bigface126 =false;
  83     width_cutoff =0.05;
  84     allow_sharp =false;
  85     allow_derivatives =allow_d;
  86     iteration_count =0;
  87     iteration_limit =0;
  88     recursion_depth =200;
  89     mono_pass = mono_pass_flag;
  90     convex_flag = convex_flag;
  91     raw_int_flag = true;
  92     eps = eps;
  93   } in
  94     recursive_verifier0(0,xx,zz,xx,zz,ff,opt);;
  95 *)
  96
  97
  98 (****************************************)
  99
 100 let domain_str x z =
 101   let s1 = map string_of_float x and
 102       s2 = map string_of_float z in
 103     sprintf "[%s], [%s]" (String.concat "; " s1) (String.concat "; " s2);;
 104
 105 let path_str p =
 106   String.concat "," (map (fun s, j -> sprintf "%s(%d)" s j) p);;
 107
 108
 109 (* get_results0 *)
 110 (* This function finds all subtrees of the given solution tree which can be
 111    veified immediately (no Result_pass_mono). These subtrees are added to
 112    the accumulator. Paths to the roots of all subtrees are also saved in
 113    the accumulator. The third returned value is a solution tree where all
 114    found subtrees are replaced with Result_pass_ref j, with j = #of the corresponding
 115    subtree in the accumulator (1-based) *)
 116
 117
 118 let get_results0 path r acc =
 119   let dummy_tree = Result_false ([], []) in
 120   let is_ref r = match r with Result_pass_ref _ -> true | _ -> false in
 121
 122   let rec get_rec path r acc =
 123     match r with
 124       | Result_mono (mono, r1) ->
 125           let get_m m = (if m.decr_flag then "ml" else "mr"), m.variable in
 126           let path' = rev_itlist (fun m l -> get_m m :: l) mono path in
 127           let flag, acc', tree = get_rec path' r1 acc in
 128             if flag then true, acc', dummy_tree
 129             else false, acc', Result_mono (mono, tree)
 130       | Result_glue (j, convex_flag, r1, r2) ->
 131           let s1, s2 = if convex_flag then "ml", "mr" else "l", "r" in
 132           let p1, p2 = ((s1, j + 1) :: path), ((s2, j + 1) :: path) in
 133           let flag1, acc1, tree1 = get_rec p1 r1 acc in
 134           let flag2, acc', tree2 = get_rec p2 r2 acc1 in
 135           let n = (length acc' + 1) in
 136             if flag1 then
 137               if flag2 then
 138                 true, acc', dummy_tree
 139               else if is_ref r1 then
 140                 false, acc', Result_glue (j, convex_flag, r1, tree2)
 141               else
 142                 false, acc' @ [rev p1, r1], Result_glue (j, convex_flag, Result_pass_ref n, tree2)
 143             else
 144               if flag2 then
 145                 if is_ref r2 then
 146                   false, acc', Result_glue (j, convex_flag, tree1, r2)
 147                 else
 148                   false, acc' @ [rev p2, r2], Result_glue (j, convex_flag, tree1, Result_pass_ref n)
 149               else
 150                 false, acc', Result_glue (j, convex_flag, tree1, tree2)
 151
 152       | Result_pass_mono _ -> false, acc, r
 153       | _ -> true, acc, dummy_tree in
 154
 155     get_rec path r acc;;
 156
 157
 158
 159
 160 (* transform_result *)
 161
 162
 163 let transform_result x z r =
 164   (* get_domain *)
 165   (* Subdivides the given domain (x,z) according to the given path *)
 166   let domain_hash = Hashtbl.create 1000 in
 167   let find_hash, mem_hash, add_hash =
 168     Hashtbl.find domain_hash, Hashtbl.mem domain_hash, Hashtbl.add domain_hash in
 169
 170   let get_domain path =
 171     let n = length x in
 172     let table f = map f (0--(n - 1)) in
 173     let rec rec_domain (x, z) path hash =
 174       match path with
 175         | [] -> x, z
 176         | (s, j) :: ps ->
 177             let hash' = hash^s^(string_of_int j) in
 178               if mem_hash hash' then
 179                 rec_domain (find_hash hash') ps hash'
 180               else
 181                 let j = j - 1 in
 182                 let domain' =
 183                   if s = "l" or s = "r" then
 184                     let ( ++ ), ( / ) = up(); upadd, updiv in
 185                     let yj = (mth x j ++ mth z j) / 2.0 in
 186                     let delta b v = table (fun i -> if i = j && b then yj else mth v i) in
 187                       if s = "l" then
 188                         delta false x, delta true z
 189                       else
 190                         delta true x, delta false z
 191                   else
 192                     if s = "ml" then
 193                       x, table (fun i -> if i = j then mth x i else mth z i)
 194                     else
 195                       table (fun i -> if i = j then mth z i else mth x i), z in
 196                 let _ = add_hash hash' domain' in
 197                   rec_domain domain' ps hash' in
 198       rec_domain (x,z) path "" in
 199
 200   (* sub_domain *)
 201   (* Verifies if interval [x',z'] SUBSET interval [x,z] *)
 202   let sub_domain (x',z') (x,z) =
 203     let le a b = itlist2 (fun a b c -> c & (a <= b)) a b true in
 204       le x x' & le z' z in
 205
 206   (* transform_pass_mono *)
 207   (* Replaces all (Result_pass_mono m) with (Result_mono [m] (Result_ref j)) where
 208      j is the reference to the corresponding domain *)
 209   let transform_pass_mono x z domains r =
 210     let domains_i = zip domains (1--length domains) in
 211
 212     let find_domain x' z' =
 213       try find (fun d, _ -> sub_domain (x', z') d) domains_i with Failure _ -> (x,z), -1 in
 214
 215     let get_m m = (if m.decr_flag then "ml" else "mr"), m.variable in
 216
 217     let rec rec_transform path r =
 218       match r with
 219         | Result_mono (mono, r1) ->
 220             let path' = rev_itlist (fun m l -> get_m m :: l) mono path in
 221               Result_mono (mono, rec_transform path' r1)
 222         | Result_glue (j, convex_flag, r1, r2) ->
 223             let s1, s2 = if convex_flag then "ml", "mr" else "l", "r" in
 224             let p1, p2 = ((s1, j + 1) :: path), ((s2, j + 1) :: path) in
 225             let t1 = rec_transform p1 r1 in
 226             let t2 = rec_transform p2 r2 in
 227               Result_glue (j, convex_flag, t1, t2)
 228         | Result_pass_mono m ->
 229             let path' = rev (get_m m :: path) in
 230             let x', z' = get_domain path' in
 231             let _, i = find_domain x' z' in
 232               (*          let _ = report (sprintf "p = %s, d = %s, found: %d"
 233                           (domain_str x' z') (path_str path') i) in *)
 234               if i >= 0 then Result_mono ([m], Result_pass_ref (-i)) else r
 235         | _ -> r in
 236
 237       rec_transform [] r in
 238
 239   let rec transform acc r =
 240     let flag, rs, r' = get_results0 [] r acc in
 241       if flag then (rs @ [[], r])
 242       else
 243         let domains = map (fun p, _ -> get_domain p) rs in
 244         let r_next = transform_pass_mono x z domains r' in
 245         let _ = r_next <> r' or failwith "transform_result: deadlock" in
 246           transform rs r_next in
 247     transform [] r;;
 248
 249
 250 (* Computes result statistics *)
 251
 252 let result_stats result =
 253   let pass = ref 0 and
 254       mono = ref 0 and
 255       glue = ref 0 and
 256       pass_mono = ref 0 and
 257       pass_raw = ref 0 and
 258       glue_convex = ref 0 in
 259
 260   let rec count r =
 261     match r with
 262       | Result_false _ -> failwith "False result"
 263       | Result_pass (flag, _, _) ->
 264           pass := !pass + 1;
 265           if flag then pass_raw := !pass_raw + 1 else ()
 266       | Result_pass_mono _ -> pass_mono := !pass_mono + 1
 267       | Result_mono (_, r1) -> mono := !mono + 1; count r1
 268       | Result_glue (_, flag, r1, r2) ->
 269           glue := !glue + 1;
 270           if flag then glue_convex := !glue_convex + 1 else ();
 271           count r1; count r2 in
 272
 273   let _ = count result in
 274     {pass = !pass; pass_raw = !pass_raw; pass_mono = !pass_mono;
 275      mono = !mono; glue = !glue; glue_convex = !glue_convex};;
 276
 277
 278 let report_stats stats =
 279   let s = sprintf "pass = %d (pass_raw = %d)\nmono = %d\nglue = %d (glue_convex = %d)\npass_mono = %d"
 280     stats.pass stats.pass_raw stats.mono stats.glue stats.glue_convex stats.pass_mono in
 281     report s;;
 282
 283
 284 let result_p_stats glue_flag p_result =
 285   let p_table = Hashtbl.create 10 in
 286   let add1 p =
 287     let c = if Hashtbl.mem p_table p then Hashtbl.find p_table p else 0 in
 288       Hashtbl.replace p_table p (succ c) in
 289
 290   let rec count r =
 291     match r with
 292       | P_result_ref _ -> ()
 293       | P_result_pass (pp, _) -> add1 pp.pp
 294       | P_result_mono (pp, _, r1) -> add1 pp.pp; count r1
 295       | P_result_glue (pp, _, _, r1, r2) ->
 296           if glue_flag then add1 pp.pp else ();
 297           count r1; count r2 in
 298
 299   let _ = count p_result in
 300   let s = Hashtbl.fold
 301     (fun p c s -> (sprintf "p = %d: %d\n" p c) ^ s) p_table "" in
 302     report s;;
 303
 304
 305 end;;