formal_lp/old/formal_interval/interval_m/verifier.hl

   1 flyspeck_needs "../formal_lp/formal_interval/interval_m/recurse.hl";;
   2 flyspeck_needs "../formal_lp/formal_interval/interval_m/recurse0.hl";;
   3
   4 module Verifier = struct
   5
   6 open Interval;;
   7 open Univariate;;
   8 open Line_interval;;
   9 open Taylor;;
  10 open Recurse;;
  11
  12
  13
  14 (****************************)
  15 (* Interval functions for the native OCaml arithmetic *)
  16
  17 type int_fun =
  18   | F_int_var of int
  19   | F_int_const of interval
  20   | F_int_pow of int * int_fun
  21   | F_int_neg of int_fun
  22   | F_int_add of int_fun * int_fun
  23   | F_int_sub of int_fun * int_fun
  24   | F_int_mul of int_fun * int_fun;;
  25
  26
  27 let ipow = Arith_misc.gen_pow imul Interval.one;;
  28
  29
  30 let eval_int_fun i_fun =
  31   fun x ->
  32     let rec eval_rec f =
  33       match f with
  34         | F_int_var i -> List.nth x (i - 1)
  35         | F_int_const int -> int
  36         | F_int_neg f1 -> ineg (eval_rec f1)
  37         | F_int_pow (n,f1) -> ipow n (eval_rec f1)
  38         | F_int_add (f1,f2) -> iadd (eval_rec f1) (eval_rec f2)
  39         | F_int_sub (f1,f2) -> isub (eval_rec f1) (eval_rec f2)
  40         | F_int_mul (f1,f2) -> imul (eval_rec f1) (eval_rec f2) in
  41       eval_rec i_fun;;
  42
  43 (**********************************)
  44 let run_test f x z min_flag min_max allow_d convex_flag mono_pass_flag raw_int_flag eps =
  45   let pad = replicate 0.0 (8 - length x) in
  46   let xx = x @ pad and zz = z @ pad in
  47   let mone = mk_interval(-1.0,-1.0) in
  48   let neg_f = Scale(f, mone) in
  49   let ff = if min_flag then
  50     Plus(neg_f, Scale(unit,mk_interval(min_max, min_max)))
  51   else
  52     Plus(f, Scale(unit, ineg (mk_interval(min_max, min_max)))) in
  53   let opt =  {
  54     only_check_deriv1_negative = false;
  55     is_using_dihmax =false;
  56     is_using_bigface126 =false;
  57     width_cutoff =0.05;
  58     allow_sharp =false;
  59     allow_derivatives =allow_d;
  60     iteration_count =0;
  61     iteration_limit =0;
  62     recursion_depth =200;
  63     mono_pass = mono_pass_flag;
  64     convex_flag = convex_flag;
  65     raw_int_flag = raw_int_flag;
  66     eps = eps;
  67   } in
  68     recursive_verifier(0,xx,zz,xx,zz,ff,opt);;
  69
  70
  71
  72 open Recurse0;;
  73
  74
  75 let run_test0 f x z min_flag min_max allow_d convex_flag mono_pass_flag eps =
  76   let pad = replicate 0.0 (8 - length x) in
  77   let xx = x @ pad and zz = z @ pad in
  78   let mone = mk_interval(-1.0,-1.0) in
  79   let neg_f = Scale(f, mone) in
  80   let ff = if min_flag then
  81     Plus(neg_f, Scale(unit,mk_interval(min_max, min_max)))
  82   else
  83     Plus(f, Scale(unit, ineg (mk_interval(min_max, min_max)))) in
  84   let opt =  {
  85     only_check_deriv1_negative = false;
  86     is_using_dihmax =false;
  87     is_using_bigface126 =false;
  88     width_cutoff =0.05;
  89     allow_sharp =false;
  90     allow_derivatives =allow_d;
  91     iteration_count =0;
  92     iteration_limit =0;
  93     recursion_depth =200;
  94     mono_pass = mono_pass_flag;
  95     convex_flag = convex_flag;
  96     raw_int_flag = true;
  97     eps = eps;
  98   } in
  99     recursive_verifier0(0,xx,zz,xx,zz,ff,opt);;
 100
 101
 102
 103 (****************************************)
 104
 105
 106
 107
 108 (*************)
 109
 110 let domain_str x z =
 111   let s1 = map string_of_float x and
 112       s2 = map string_of_float z in
 113     sprintf "[%s], [%s]" (String.concat "; " s1) (String.concat "; " s2);;
 114
 115 let path_str p =
 116   String.concat "," (map (fun s, j -> sprintf "%s(%d)" s j) p);;
 117
 118
 119 (* get_results0 *)
 120 (* This function finds all subtrees of the given solution tree which can be
 121    veified immediately (no Result_pass_mono). These subtrees are added to
 122    the accumulator. Paths to the roots of all subtrees are also saved in
 123    the accumulator. The third returned value is a solution tree where all
 124    found subtrees are replaced with Result_pass_ref j, with j = #of the corresponding
 125    subtree in the accumulator (1-based) *)
 126
 127
 128 let get_results0 path r acc =
 129   let dummy_tree = Result_false ([], []) in
 130   let is_ref r = match r with Result_pass_ref _ -> true | _ -> false in
 131
 132   let rec get_rec path r acc =
 133     match r with
 134       | Result_mono (mono, r1) ->
 135           let get_m m = (if m.decr_flag then "ml" else "mr"), m.variable in
 136           let path' = rev_itlist (fun m l -> get_m m :: l) mono path in
 137           let flag, acc', tree = get_rec path' r1 acc in
 138             if flag then true, acc', dummy_tree
 139             else false, acc', Result_mono (mono, tree)
 140       | Result_glue (j, convex_flag, r1, r2) ->
 141           let s1, s2 = if convex_flag then "ml", "mr" else "l", "r" in
 142           let p1, p2 = ((s1, j + 1) :: path), ((s2, j + 1) :: path) in
 143           let flag1, acc1, tree1 = get_rec p1 r1 acc in
 144           let flag2, acc', tree2 = get_rec p2 r2 acc1 in
 145           let n = (length acc' + 1) in
 146             if flag1 then
 147               if flag2 then
 148                 true, acc', dummy_tree
 149               else if is_ref r1 then
 150                 false, acc', Result_glue (j, convex_flag, r1, tree2)
 151               else
 152                 false, acc' @ [rev p1, r1], Result_glue (j, convex_flag, Result_pass_ref n, tree2)
 153             else
 154               if flag2 then
 155                 if is_ref r2 then
 156                   false, acc', Result_glue (j, convex_flag, tree1, r2)
 157                 else
 158                   false, acc' @ [rev p2, r2], Result_glue (j, convex_flag, tree1, Result_pass_ref n)
 159               else
 160                 false, acc', Result_glue (j, convex_flag, tree1, tree2)
 161
 162       | Result_pass_mono _ -> false, acc, r
 163       | _ -> true, acc, dummy_tree in
 164
 165     get_rec path r acc;;
 166
 167
 168
 169
 170 (* transform_result *)
 171
 172
 173 let transform_result x z r =
 174   (* get_domain *)
 175   (* Subdivides the given domain (x,z) according to the given path *)
 176   let domain_hash = Hashtbl.create 1000 in
 177   let find_hash, mem_hash, add_hash =
 178     Hashtbl.find domain_hash, Hashtbl.mem domain_hash, Hashtbl.add domain_hash in
 179
 180   let get_domain path =
 181     let n = length x in
 182     let table f = map f (0--(n - 1)) in
 183     let rec rec_domain (x, z) path hash =
 184       match path with
 185         | [] -> x, z
 186         | (s, j) :: ps ->
 187             let hash' = hash^s^(string_of_int j) in
 188               if mem_hash hash' then
 189                 rec_domain (find_hash hash') ps hash'
 190               else
 191                 let j = j - 1 in
 192                 let domain' =
 193                   if s = "l" or s = "r" then
 194                     let ( ++ ), ( / ) = up(); upadd, updiv in
 195                     let yj = (mth x j ++ mth z j) / 2.0 in
 196                     let delta b v = table (fun i -> if i = j && b then yj else mth v i) in
 197                       if s = "l" then
 198                         delta false x, delta true z
 199                       else
 200                         delta true x, delta false z
 201                   else
 202                     if s = "ml" then
 203                       x, table (fun i -> if i = j then mth x i else mth z i)
 204                     else
 205                       table (fun i -> if i = j then mth z i else mth x i), z in
 206                 let _ = add_hash hash' domain' in
 207                   rec_domain domain' ps hash' in
 208       rec_domain (x,z) path "" in
 209
 210   (* sub_domain *)
 211   (* Verifies if interval [x',z'] SUBSET interval [x,z] *)
 212   let sub_domain (x',z') (x,z) =
 213     let le a b = itlist2 (fun a b c -> c & (a <= b)) a b true in
 214       le x x' & le z' z in
 215
 216   (* transform_pass_mono *)
 217   (* Replaces all (Result_pass_mono m) with (Result_mono [m] (Result_ref j)) where
 218      j is the reference to the corresponding domain *)
 219   let transform_pass_mono x z domains r =
 220     let domains_i = zip domains (1--length domains) in
 221
 222     let find_domain x' z' =
 223       try find (fun d, _ -> sub_domain (x', z') d) domains_i with Failure _ -> (x,z), -1 in
 224
 225     let get_m m = (if m.decr_flag then "ml" else "mr"), m.variable in
 226
 227     let rec rec_transform path r =
 228       match r with
 229         | Result_mono (mono, r1) ->
 230             let path' = rev_itlist (fun m l -> get_m m :: l) mono path in
 231               Result_mono (mono, rec_transform path' r1)
 232         | Result_glue (j, convex_flag, r1, r2) ->
 233             let s1, s2 = if convex_flag then "ml", "mr" else "l", "r" in
 234             let p1, p2 = ((s1, j + 1) :: path), ((s2, j + 1) :: path) in
 235             let t1 = rec_transform p1 r1 in
 236             let t2 = rec_transform p2 r2 in
 237               Result_glue (j, convex_flag, t1, t2)
 238         | Result_pass_mono m ->
 239             let path' = rev (get_m m :: path) in
 240             let x', z' = get_domain path' in
 241             let _, i = find_domain x' z' in
 242               (*          let _ = report (sprintf "p = %s, d = %s, found: %d"
 243                           (domain_str x' z') (path_str path') i) in *)
 244               if i >= 0 then Result_mono ([m], Result_pass_ref (-i)) else r
 245         | _ -> r in
 246
 247       rec_transform [] r in
 248
 249   let rec transform acc r =
 250     let flag, rs, r' = get_results0 [] r acc in
 251       if flag then (rs @ [[], r])
 252       else
 253         let domains = map (fun p, _ -> get_domain p) rs in
 254         let r_next = transform_pass_mono x z domains r' in
 255         let _ = r_next <> r' or failwith "transform_result: deadlock" in
 256           transform rs r_next in
 257     transform [] r;;
 258
 259
 260
 261
 262
 263
 264 (*
 265 let pp = 8;;
 266 let n = 7;;
 267 let xx = `[-- &1; -- &1; -- &1; -- &1; -- &1; -- &1; -- &1]` and
 268     zz = `[&1; &1; &1; &1; &1; &1; &1]`;;
 269
 270 let xx1 = convert_to_float_list pp true xx and
 271     zz1 = convert_to_float_list pp false zz;;
 272 let xx_float = map float_of_float_tm (dest_list xx1) and
 273     zz_float = map float_of_float_tm (dest_list zz1);;
 274
 275 let eval0_magnetism, eval_magnetism, tf_magnetism =
 276   mk_verification_functions pp magnetism_poly true magnetism_min;;
 277 let c1 = run_test tf_magnetism xx_float zz_float false 0.0 true false true;;
 278 let c0 = run_test0 tf_magnetism xx_float zz_float false 0.0 true false true;;
 279 (* 77 *)
 280 result_size c1;;
 281 (* 121 *)
 282 result_size c0;;
 283
 284 let r = transform_result xx_float zz_float c1;;
 285 r;;
 286 length r;;
 287 let x = map (fun _, r ->
 288        match r with | Result_pass_ref j -> 1 | _ -> 0) r;;
 289 itlist (+) x 0;;
 290
 291 map (fun _, r ->
 292        match r with
 293          | Result_pass _ -> "pass"
 294          | Result_mono _ -> "mono"
 295          | Result_glue _ -> "glue"
 296          | Result_false _ -> "_|_"
 297          | Result_pass_mono _ -> "pass_mono"
 298          | Result_pass_ref _ -> "ref"
 299     ) r;;
 300
 301
 302 length it;;
 303 transform_result xx_float zz_float c0;;
 304
 305 let pp = 5;;
 306 let n = 2;;
 307
 308 let poly_tm = expr_to_vector_fun `x1 pow 2 + x2 pow 2`;;
 309
 310 let xx = `[-- &6; -- &6]` and
 311     zz = `[&1; &2]`;;
 312
 313 let xx1 = convert_to_float_list pp true xx and
 314     zz1 = convert_to_float_list pp false zz;;
 315 let xx_float = map float_of_float_tm (dest_list xx1) and
 316     zz_float = map float_of_float_tm (dest_list zz1);;
 317
 318 let eval0_poly, eval_poly, tf_poly = mk_verification_functions pp poly_tm true `-- &1`;;
 319 let c1 = run_test tf_poly xx_float zz_float false 0.0 true false true;;
 320
 321 transform_result xx_float zz_float c1;;
 322
 323 let c0 = run_test0 tf_poly xx_float zz_float false 0.0 true false true;;
 324 (* 22 *)
 325 result_size c1;;
 326 (* 38 *)
 327 result_size c0;;
 328 m_verify_raw n pp eval0_poly eval_poly c1 xx1 zz1;;
 329 m_verify_raw n pp eval0_poly eval_poly c1 xx1 zz1;;
 330
 331
 332
 333 transform_result xx_float zz_float c1;;
 334
 335 let flag, rs1, r1 = get_results0 [] c1 [];;
 336 rs1;;
 337 r1;;
 338 let domains = map (fun p, _ -> get_domain xx_float zz_float p) rs1;;
 339 let r2 = transform_pass_mono xx_float zz_float domains r1;;
 340
 341 let flag, rs3, r3 = get_results0 [] r2 rs1;;
 342 rs3;;
 343 r3;;
 344 let domains = map (fun p, _ -> get_domain xx_float zz_float p) rs3;;
 345 let r4 = transform_pass_mono xx_float zz_float domains r3;;
 346
 347 let flag, rs5, r5 = get_results0 [] r4 rs3;;
 348 rs5;;
 349 r5;;
 350
 351
 352
 353 test 1000 (get_results0 [] c1) [];;
 354 test 1000 (map (fun p, _ -> get_domain xx_float zz_float p)) rs1;;
 355 test 1000 (transform_pass_mono xx_float zz_float domains) r1;;
 356
 357 let flag, rs3, r3 = get_results0 [] r2 rs1;;
 358 let domains = map (fun p, _ -> get_domain xx_float zz_float p) rs3;;
 359 let r4 = transform_pass_mono xx_float zz_float domains r3;;
 360
 361 let flag, rs5, r5 = get_results0 [] r4 rs3;;
 362 *)
 363
 364
 365
 366 let result_stat result =
 367   let pass = ref 0 and
 368       mono = ref 0 and
 369       glue = ref 0 and
 370       pass_mono = ref 0 and
 371       pass_raw = ref 0 and
 372       glue_convex = ref 0 in
 373
 374   let rec count r =
 375     match r with
 376       | Result_false _ -> failwith "False result"
 377       | Result_pass (flag, _, _) ->
 378           pass := !pass + 1;
 379           if flag then pass_raw := !pass_raw + 1 else ()
 380       | Result_pass_mono _ -> pass_mono := !pass_mono + 1
 381       | Result_mono (_, r1) -> mono := !mono + 1; count r1
 382       | Result_glue (_, flag, r1, r2) ->
 383           glue := !glue + 1;
 384           if flag then glue_convex := !glue_convex + 1 else ();
 385           count r1; count r2 in
 386
 387   let _ = count result in
 388   let s = sprintf "pass = %d (pass_raw = %d)\nmono = %d\nglue = %d (glue_convex = %d)\npass_mono = %d"
 389     !pass !pass_raw !mono !glue !glue_convex !pass_mono in
 390     report s;;
 391
 392
 393 let result_p_stat glue_flag p_result =
 394   let p_table = Hashtbl.create 10 in
 395   let add1 p =
 396     let c = if Hashtbl.mem p_table p then Hashtbl.find p_table p else 0 in
 397       Hashtbl.replace p_table p (succ c) in
 398
 399   let rec count r =
 400     match r with
 401       | P_result_ref _ -> ()
 402       | P_result_pass (pp, _) -> add1 pp.pp
 403       | P_result_mono (pp, _, r1) -> add1 pp.pp; count r1
 404       | P_result_glue (pp, _, _, r1, r2) ->
 405           if glue_flag then add1 pp.pp else ();
 406           count r1; count r2 in
 407
 408   let _ = count p_result in
 409   let s = Hashtbl.fold
 410     (fun p c s -> (sprintf "p = %d: %d\n" p c) ^ s) p_table "" in
 411     report s;;
 412
 413
 414 end;;