formal_ineqs/informal/informal_m_verifier.hl

   1 (* =========================================================== *)
   2 (* Informal verification procedures                            *)
   3 (* Author: Alexey Solovyev                                     *)
   4 (* Date: 2012-10-27                                            *)
   5 (* =========================================================== *)
   6
   7 (* Dependencies *)
   8 needs "informal/informal_m_taylor.hl";;
   9 needs "verifier/interval_m/recurse.hl";;
  10 needs "verifier_options.hl";;
  11
  12 module Informal_verifier = struct
  13
  14 open Informal_float;;
  15 open Informal_interval;;
  16 open Informal_taylor;;
  17 open Recurse;;
  18 open Verifier_options;;
  19
  20
  21 type verification_funs =
  22 {
  23   (* p_lin -> p_second -> dom -> ti *)
  24   taylor : int -> int -> m_cell_domain -> m_taylor_interval;
  25   (* pp -> xx -> zz -> interval *)
  26   f : int -> ifloat list -> ifloat list -> interval;
  27   (* j -> pp -> xx -> zz -> interval *)
  28   df : int -> int -> ifloat list -> ifloat list -> interval;
  29   (* i j -> pp -> xx -> zz -> interval *)
  30   ddf : int -> int -> int -> ifloat list -> ifloat list -> interval;
  31 };;
  32
  33
  34 (* m_subset_interval *)
  35 let m_subset_interval a b c d =
  36   let prove_le l1 l2 = itlist2 (fun x y r -> le_float x y && r) l1 l2 true in
  37     prove_le a c && prove_le d b;;
  38
  39 (* m_taylor_cell_pass *)
  40 let m_taylor_cell_pass pp ti =
  41   let upper = eval_m_taylor_upper_bound pp ti in
  42     lt0_float upper;;
  43
  44 (* m_taylor_cell_pass0 *)
  45 let m_taylor_cell_pass0 int =
  46   (lt0_float o snd o dest_interval) int;;
  47
  48 (* m_cell_pass_subdomain *)
  49 let m_cell_pass_subdomain domain2 pass_domain =
  50   let a, b = pass_domain.lo, pass_domain. hi in
  51   let c, d = domain2.lo, domain2.hi in
  52     m_subset_interval a b c d;;
  53
  54 (* m_incr_pass *)
  55 let m_incr_pass pp j ti =
  56   let partial_bound = eval_m_taylor_partial_lower pp j ti in
  57     ge0_float partial_bound;;
  58
  59 (* m_decr_pass *)
  60 let m_decr_pass pp j ti =
  61   let partial_bound = eval_m_taylor_partial_upper pp j ti in
  62     le0_float partial_bound;;
  63
  64 (* m_mono_pass_gen *)
  65 let m_mono_pass_gen decr_flag bound =
  66   (if decr_flag then le0_float else ge0_float) bound;;
  67
  68 (* m_convex_pass *)
  69 let m_convex_pass int =
  70   (ge0_float o fst o dest_interval) int;;
  71
  72
  73 (* mk_verification_functions *)
  74 let mk_verification_functions_poly pp0 f partials partials2 =
  75   let n = length partials in
  76   let taylor = eval_m_taylor pp0 f partials partials2 in
  77   let eval0 = mk_eval_function pp0 f in
  78   let eval1 = map (fun i -> mk_eval_function pp0 ((rand o concl o List.nth partials) (i - 1))) (1--n) in
  79   let eval2 = map (fun i ->
  80                      map (fun j ->
  81                             let d2 = List.nth (List.nth partials2 (i - 1)) (j - 1) in
  82                               mk_eval_function pp0 ((rand o concl) d2)) (1--i)) (1--n) in
  83     {
  84       taylor = taylor;
  85       f = eval0;
  86       df = (fun i -> List.nth eval1 (i - 1));
  87       ddf = (fun i j -> List.nth (List.nth eval2 (j - 1)) (i - 1));
  88     };;
  89
  90
  91 (* split_domain *)
  92 let split_domain pp j domain =
  93   let n = length domain.w in
  94   let t = List.nth domain.y (j - 1) in
  95   let vv = map (fun i -> if i = j then t else List.nth domain.hi (i - 1)) (1--n) in
  96   let uu = map (fun i -> if i = j then t else List.nth domain.lo (i - 1)) (1--n) in
  97     mk_m_center_domain pp domain.lo vv, mk_m_center_domain pp uu domain.hi;;
  98
  99
 100 (* restrict_domain *)
 101 let restrict_domain j left_flag domain =
 102   let replace list j v = map (fun i -> if i = j then v else List.nth list (i - 1)) (1--length list) in
 103   let t = List.nth (if left_flag then domain.lo else domain.hi) (j - 1) in
 104   let lo = if left_flag then domain.lo else replace domain.lo j t in
 105   let hi = if left_flag then replace domain.hi j t else domain.hi in
 106   let w = replace domain.w j float_0 in
 107   let y = replace domain.y j t in
 108     {lo = lo; hi = hi; w = w; y = y};;
 109
 110
 111 (*****************************)
 112 (* aux list functions *)
 113
 114 (* Merges lists of indices. Input lists must be sorted *)
 115 let merge_indices l1 l2 = uniq (merge (<) l1 l2);;
 116
 117 (* Returns all elements at the given indices *)
 118 let take_all l inds = map (List.nth l) inds;;
 119
 120
 121 (*****************************)
 122 (* m_verify_raw *)
 123
 124 (* Constructs a p_result_tree from the given result_tree *)
 125 let m_verify_raw (report_start, total_size) p_split p_min p_max fs_list certificate domain0 ref_list =
 126   let r_size = result_size certificate in
 127   let r_size2 = float_of_int (if total_size > 0 then total_size else (if r_size > 0 then r_size else 1)) in
 128   let k = ref 0 in
 129   let kk = ref report_start in
 130   let last_report = ref (int_of_float (float_of_int !kk /. r_size2 *. 100.0)) in
 131   let ps = p_min -- p_max in
 132
 133   (* finds an optimal precision value *)
 134   let rec find_p p_fun p_list =
 135     match p_list with
 136       | [] -> failwith "find_p: no good p found"
 137       | p :: ps ->
 138           let flag = (try p_fun p with Failure _ -> false | Division_by_zero -> false) in
 139             if flag then
 140               let _ = if !info_print_level >= 2 then report (sprintf "p = %d" p) else () in p
 141             else find_p p_fun ps in
 142
 143   (* pass_test *)
 144   let pass_test domain (j,f0_flag) pp =
 145     let fs = List.nth fs_list j in
 146       if f0_flag then
 147         m_taylor_cell_pass0 (fs.f pp domain.lo domain.hi)
 148       else
 149         m_taylor_cell_pass pp (fs.taylor pp pp domain) in
 150
 151   (* glue_test *)
 152   let glue_test domain i convex_flag inds pp =
 153     let fss = take_all fs_list inds in
 154       if convex_flag then
 155         forall (fun fs -> m_convex_pass (fs.ddf (i + 1) (i + 1) pp domain.lo domain.hi)) fss
 156       else
 157         true in
 158
 159   (* mono_test *)
 160   let mono_test mono domain domains inds pp =
 161     let fss = take_all fs_list inds in
 162     let xx, zz = domain.lo, domain.hi in
 163     let taylors = map (fun fs -> fs.taylor pp pp domain) fss in
 164     let gen_mono m =
 165       if m.df0_flag then
 166         if m.decr_flag then
 167           map (fun fs -> (snd o dest_interval) (fs.df m.variable pp xx zz)) fss
 168         else
 169           map (fun fs -> (fst o dest_interval) (fs.df m.variable pp xx zz)) fss
 170       else
 171         if m.decr_flag then
 172           map (eval_m_taylor_partial_upper pp m.variable) taylors
 173         else
 174           map (eval_m_taylor_partial_lower pp m.variable) taylors in
 175     let monos = map gen_mono mono in
 176       rev_itlist (fun (m, bounds) pass ->
 177                     let flag = m.decr_flag in
 178                       forall (m_mono_pass_gen flag) bounds && pass) (rev (zip mono monos)) true in
 179
 180   (* mk_domains *)
 181   let rec mk_domains mono dom0 acc =
 182     match mono with
 183       | [] -> rev acc
 184       | m :: ms ->
 185           let j, flag = m.variable, m.decr_flag in
 186           let dom = restrict_domain j flag dom0 in
 187             mk_domains ms dom (dom :: acc) in
 188
 189   (* rec_verify *)
 190   let rec rec_verify domain certificate =
 191     match certificate with
 192       | Result_mono (mono, r1) ->
 193           let _ =
 194             if !info_print_level >= 2 then
 195               let mono_strs =
 196                 map (fun m -> sprintf "%s%d (%b)" (if m.decr_flag then "-" else "")
 197                        m.variable m.df0_flag) mono in
 198                 report (sprintf "Mono: [%s]" (String.concat ";" mono_strs))
 199             else () in
 200           let domains = mk_domains mono domain [] in
 201           let tree1, inds = rec_verify (last domains) r1 in
 202             (try
 203                let pp = find_p (mono_test mono domain domains inds) ps in
 204                  P_result_mono ({pp = pp}, mono, tree1), inds
 205              with Failure _ -> failwith "mono: failed")
 206
 207       | Result_pass (j, f0_flag) ->
 208           let _ = k := !k + 1; kk := !kk + 1 in
 209           let _ =
 210             if !info_print_level >= 2 then
 211               report (sprintf "Verifying: %d/%d (f0_flag = %b)" !k r_size f0_flag)
 212             else () in
 213                 let _ = !info_print_level <> 1 or
 214                         (let r = int_of_float (float_of_int !kk /. r_size2 *. 100.0) in
 215                         let _ = if r <> !last_report then (last_report := r; report0 (sprintf "%d " r)) else () in true) in
 216
 217             (try
 218                let pp = find_p (pass_test domain (j, f0_flag)) ps in
 219                  P_result_pass ({pp = pp}, j, f0_flag), [j]
 220              with Failure _ -> failwith "pass: failed")
 221
 222       | Result_glue (i, convex_flag, r1, r2) ->
 223           let domain1, domain2 =
 224             if convex_flag then
 225               let d1 = restrict_domain (i + 1) true domain in
 226               let d2 = restrict_domain (i + 1) false domain in
 227                 d1, d2
 228             else
 229               split_domain p_split (i + 1) domain in
 230           let tree1, inds1 = rec_verify domain1 r1 in
 231           let tree2, inds2 = rec_verify domain2 r2 in
 232           let inds = merge_indices inds1 inds2 in
 233             (try
 234                let pp = find_p (glue_test domain i convex_flag inds) ps in
 235                  P_result_glue ({pp = pp}, i, convex_flag, tree1, tree2), inds
 236              with Failure _ -> failwith "glue: failed")
 237
 238       | Result_pass_ref i ->
 239           let _ = if !info_print_level >= 2 then report (sprintf "Ref: %d" i) else () in
 240           let pass_flag, inds =
 241             if i > 0 then
 242               let _, inds = List.nth ref_list (i - 1) in
 243                 true, inds
 244             else
 245               let pass_domain, inds = List.nth ref_list (-i - 1) in
 246                 m_cell_pass_subdomain domain pass_domain, inds in
 247             if not pass_flag then
 248               failwith "ref: failed"
 249             else
 250               P_result_ref i, inds
 251
 252       | _ -> failwith "False result" in
 253
 254     rec_verify domain0 certificate;;
 255
 256
 257
 258 (*****************)
 259
 260 (* m_verify_raw0 *)
 261 let m_verify_raw0 p_split p_min p_max fs_list certificate xx zz =
 262   m_verify_raw (0, 0) p_split p_min p_max fs_list certificate (mk_m_center_domain p_split xx zz) [];;
 263
 264
 265 (* m_verify_list *)
 266 let m_verify_list p_split p_min p_max fs_list certificate_list xx zz =
 267   let domain_hash = Hashtbl.create (length certificate_list * 10) in
 268   let mem, find, add = Hashtbl.mem domain_hash,
 269     Hashtbl.find domain_hash, Hashtbl.add domain_hash in
 270
 271   let get_m_cell_domain pp domain0 path =
 272     let rec get_rec domain path hash =
 273       match path with
 274         | [] -> domain
 275         | (s, j) :: ps ->
 276             let hash' = hash^s^(string_of_int j) in
 277               if mem hash' then
 278                 get_rec (find hash') ps hash'
 279               else
 280                 if s = "l" or s = "r" then
 281                   let domain1, domain2 = split_domain pp j domain in
 282                   let hash1 = hash^"l"^(string_of_int j) and
 283                       hash2 = hash^"r"^(string_of_int j) in
 284                   let _ = add hash1 domain1; add hash2 domain2 in
 285                     if s = "l" then
 286                       get_rec domain1 ps hash'
 287                     else
 288                       get_rec domain2 ps hash'
 289                 else
 290                   let l_flag = (s = "ml") in
 291                   let domain' = restrict_domain j l_flag domain in
 292                   let _ = add hash' domain' in
 293                     get_rec domain' ps hash' in
 294       get_rec domain0 path "" in
 295
 296   let domain0 = mk_m_center_domain p_split xx zz in
 297   let size = length certificate_list in
 298   let k = ref 0 in
 299   let kk = ref 0 in
 300   let total_size = end_itlist (+) (map (result_size o snd) certificate_list) in
 301
 302   let rec rec_verify certificate_list dom_list tree_list =
 303     match certificate_list with
 304       | [] -> rev tree_list
 305       | (path, certificate) :: cs ->
 306           let _ = k := !k + 1 in
 307           let _ = !info_print_level < 2 or (report (sprintf "List: %d/%d" !k size); true) in
 308           let domain = get_m_cell_domain p_split domain0 path in
 309           let tree, inds = m_verify_raw (!kk, total_size) p_split p_min p_max fs_list certificate domain dom_list in
 310           let _ = kk := !kk + result_size certificate in
 311             rec_verify cs (dom_list @ [domain, inds]) ((path, tree) :: tree_list) in
 312     rec_verify certificate_list [] [];;
 313
 314 end;;
 315
 316
 317