formal_ineqs/verifier/interval_m/verifier.hl

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