formal_ineqs/informal/informal_m_taylor.hl

   1 (* =========================================================== *)
   2 (* Informal taylor intervals                                   *)
   3 (* Author: Alexey Solovyev                                     *)
   4 (* Date: 2012-10-27                                            *)
   5 (* =========================================================== *)
   6
   7 (* Dependencies *)
   8 needs "informal/informal_arith.hl";;
   9 needs "informal/informal_eval_interval.hl";;
  10
  11
  12 module Informal_taylor = struct
  13
  14 open Informal_interval;;
  15 open Informal_float;;
  16 open Informal_atn;;
  17 open Informal_eval_interval;;
  18
  19
  20 type m_cell_domain =
  21 {
  22   lo : ifloat list;
  23   hi : ifloat list;
  24   y : ifloat list;
  25   w : ifloat list;
  26 };;
  27
  28
  29 type m_taylor_interval =
  30 {
  31   n : int;
  32   domain : m_cell_domain;
  33   f : interval;
  34   df : interval list;
  35   ddf : interval list list;
  36 };;
  37
  38
  39 let float_0 = mk_small_num_float 0 and
  40     float_1 = mk_small_num_float 1 and
  41     float_2 = mk_small_num_float 2;;
  42
  43 let float_inv2 = div_float_lo 1 float_1 float_2;;
  44
  45 (* convert_to_float_list *)
  46 let convert_to_float_list pp lo_flag list_tm =
  47   let tms = dest_list list_tm in
  48   let i_funs = map build_interval_fun tms in
  49   let ints = map (fun f -> eval_interval_fun pp f [] []) i_funs in
  50   let extract = (if lo_flag then fst else snd) o dest_interval in
  51     map extract ints;;
  52
  53
  54 (* mk_m_center_domain *)
  55 let mk_m_center_domain pp x_list z_list =
  56   let y_list =
  57     let ( * ), (+) = mul_float_eq, add_float_hi pp in
  58       map2 (fun x z -> if eq_float x z then x else float_inv2 * (x + z)) x_list z_list in
  59
  60   (* test: x <= y <= z *)
  61   let flag1 = itlist2 (fun x y a -> le_float x y && a) x_list y_list true and
  62       flag2 = itlist2 (fun y z a -> le_float y z && a) y_list z_list true in
  63     if not flag1 or not flag2 then
  64       failwith "mk_m_center_domain: ~(x <= y <= z)"
  65     else
  66       let w_list =
  67         let (-) = sub_float_hi pp in
  68         let w1 = map2 (-) y_list x_list in
  69         let w2 = map2 (-) z_list y_list in
  70           map2 max_float w1 w2 in
  71         {lo = x_list; hi = z_list; y = y_list; w = w_list};;
  72
  73
  74 (* eval_m_taylor (pp0 for initial evaluation of constants) *)
  75 let eval_m_taylor pp0 f_tm partials partials2 =
  76   let build = eval_constants pp0 o build_interval_fun o snd o dest_abs in
  77   let f = build f_tm in
  78   let n = length partials in
  79   (* Verify that the list of second partial derivatives is correct *)
  80   let _ =  map2 (fun i list -> if length list <> i then
  81                    failwith "eval_m_taylor: incorrect partials2" else ()) (1--n) partials2 in
  82   let dfs = map (build o rand o concl) partials in
  83   let d2fs = map (build o rand o concl) (List.flatten partials2) in
  84   let f_dfs_list = find_and_replace_all (f :: dfs) [] in
  85   let rec shape_list dd i =
  86     if i >= n then [dd] else
  87       let l1, l2 = chop_list i dd in
  88         l1 :: shape_list l2 (i + 1) in
  89   let d2fs_list = find_and_replace_all d2fs [] in
  90     fun p_lin p_second domain ->
  91       let y_ints = map (fun y -> mk_interval (y, y)) domain.y in
  92       let xz_ints = map mk_interval (zip domain.lo domain.hi) in
  93       let f_dfs_vals = eval_interval_fun_list p_lin f_dfs_list y_ints in
  94       let d2fs_vals = eval_interval_fun_list p_second d2fs_list xz_ints in
  95         {n = n; domain = domain;
  96          f = hd f_dfs_vals; df = tl f_dfs_vals;
  97         ddf = shape_list d2fs_vals 1};;
  98
  99
 100 (* mk_eval_functionq *)
 101 let mk_eval_function pp0 f_tm =
 102   let build = eval_constants pp0 o build_interval_fun o snd o dest_abs in
 103   let f = build f_tm in
 104   let f_list = find_and_replace_all [f] [] in
 105     fun pp x_list z_list ->
 106       let xz_ints = map mk_interval (zip x_list z_list) in
 107       let f_val = eval_interval_fun_list pp f_list xz_ints in
 108         hd f_val;;
 109
 110
 111 (* error_mul_f2_hi *)
 112 let error_mul_f2_hi pp a int = mul_float_hi pp a (abs_interval int);;
 113
 114
 115 (* eval_m_taylor_error *)
 116 (* sum_{i = 1}^n (w_i * (f_ii * w_i + 2 * sum_{j = 1}^{i - 1} w_j * f_ij)) *)
 117 let eval_m_taylor_error pp ti =
 118   let w = ti.domain.w in
 119   let ns = 1--ti.n in
 120   let ( * ), ( + ) = mul_float_hi pp, add_float_hi pp in
 121   let mul_wdd = map2 (fun list i -> Arith_misc.my_map2 (error_mul_f2_hi pp) w list) ti.ddf ns in
 122   let sums1 = map (end_itlist ( + ) o butlast) (tl mul_wdd) in
 123   let sums2 = (hd o hd) mul_wdd :: map2 (fun list t1 -> last list + float_2 * t1) (tl mul_wdd) sums1 in
 124   let sums = map2 ( * ) w sums2 in
 125     end_itlist ( + ) sums;;
 126
 127
 128 (* eval_m_taylor_upper_bound *)
 129 let eval_m_taylor_upper_bound pp ti =
 130   let f_hi = (snd o dest_interval) ti.f in
 131   let error = eval_m_taylor_error pp ti in
 132   let ( * ), ( + ) = mul_float_hi pp, add_float_hi pp in
 133   let sum2 =
 134     let mul_wd = map2 (error_mul_f2_hi pp) ti.domain.w ti.df in
 135       end_itlist ( + ) mul_wd in
 136   let a = sum2 + float_inv2 * error in
 137     f_hi + a;;
 138
 139 (* eval_m_taylor_lower_bound *)
 140 let eval_m_taylor_lower_bound pp ti =
 141   let f_lo = (fst o dest_interval) ti.f in
 142   let error = eval_m_taylor_error pp ti in
 143   let ( * ), ( + ), ( - ) = mul_float_hi pp, add_float_hi pp, sub_float_lo pp in
 144   let sum2 =
 145     let mul_wd = map2 (error_mul_f2_hi pp) ti.domain.w ti.df in
 146       end_itlist ( + ) mul_wd in
 147   let a = sum2 + float_inv2 * error in
 148     f_lo - a;;
 149
 150
 151 (* eval_m_taylor_bound *)
 152 let eval_m_taylor_bound pp ti =
 153   let f_lo, f_hi = dest_interval ti.f in
 154   let error = eval_m_taylor_error pp ti in
 155   let ( * ), ( + ), ( - ) = mul_float_hi pp, add_float_hi pp, sub_float_lo pp in
 156   let sum2 =
 157     let mul_wd = map2 (error_mul_f2_hi pp) ti.domain.w ti.df in
 158       end_itlist ( + ) mul_wd in
 159   let a = sum2 + float_inv2 * error in
 160   let hi = f_hi + a in
 161   let lo = f_lo - a in
 162     mk_interval (lo, hi);;
 163
 164
 165 (* eval_m_taylor_partial_upper *)
 166 let eval_m_taylor_partial_upper pp i ti =
 167   let df_hi = (snd o dest_interval o List.nth ti.df) (i - 1) in
 168   let dd_list = map (fun j -> if j <= i then
 169                        List.nth (List.nth ti.ddf (i - 1)) (j - 1)
 170                      else
 171                        List.nth (List.nth ti.ddf (j - 1)) (i - 1)) (1--ti.n) in
 172   let sum2 =
 173     let mul_dd = map2 (error_mul_f2_hi pp) ti.domain.w dd_list in
 174       end_itlist (add_float_hi pp) mul_dd in
 175     add_float_hi pp df_hi sum2;;
 176
 177
 178 (* eval_m_taylor_partial_lower *)
 179 let eval_m_taylor_partial_lower pp i ti =
 180   let df_lo = (fst o dest_interval o List.nth ti.df) (i - 1) in
 181   let dd_list = map (fun j -> if j <= i then
 182                        List.nth (List.nth ti.ddf (i - 1)) (j - 1)
 183                      else
 184                        List.nth (List.nth ti.ddf (j - 1)) (i - 1)) (1--ti.n) in
 185   let sum2 =
 186     let mul_dd = map2 (error_mul_f2_hi pp) ti.domain.w dd_list in
 187       end_itlist (add_float_hi pp) mul_dd in
 188     sub_float_lo pp df_lo sum2;;
 189
 190
 191 (* eval_m_taylor_partial_bound *)
 192 let eval_m_taylor_partial_bound pp i ti =
 193   let df_lo, df_hi = (dest_interval o List.nth ti.df) (i - 1) in
 194   let dd_list = map (fun j -> if j <= i then
 195                        List.nth (List.nth ti.ddf (i - 1)) (j - 1)
 196                      else
 197                        List.nth (List.nth ti.ddf (j - 1)) (i - 1)) (1--ti.n) in
 198   let sum2 =
 199     let mul_dd = map2 (error_mul_f2_hi pp) ti.domain.w dd_list in
 200       end_itlist (add_float_hi pp) mul_dd in
 201   let lo = sub_float_lo pp df_lo sum2 in
 202   let hi = add_float_hi pp df_hi sum2 in
 203     mk_interval (lo, hi);;
 204
 205
 206 (* add *)
 207 let eval_m_taylor_add p_lin p_second taylor1 taylor2 =
 208   let ( + ), ( ++ ) = add_interval p_lin, add_interval p_second in
 209     {
 210       n = taylor1.n;
 211       domain = taylor1.domain;
 212       f = taylor1.f + taylor2.f;
 213       df = map2 (+) taylor1.df taylor2.df;
 214       ddf = map2 (map2 (++)) taylor1.ddf taylor2.ddf
 215     };;
 216
 217
 218 (* sub *)
 219 let eval_m_taylor_sub p_lin p_second taylor1 taylor2 =
 220   let ( - ), ( -- ) = sub_interval p_lin, sub_interval p_second in
 221     {
 222       n = taylor1.n;
 223       domain = taylor1.domain;
 224       f = taylor1.f - taylor2.f;
 225       df = map2 (-) taylor1.df taylor2.df;
 226       ddf = map2 (map2 (--)) taylor1.ddf taylor2.ddf
 227     };;
 228
 229
 230 (* mul *)
 231 let eval_m_taylor_mul p_lin p_second ti1 ti2 =
 232   let n = ti1.n in
 233   let ns = 1--n in
 234   let bounds = mul_interval p_lin ti1.f ti2.f in
 235   let df = map2 (fun d1 d2 ->
 236                    let ( * ), ( + ) = mul_interval p_lin, add_interval p_lin in
 237                      d1 * ti2.f + ti1.f * d2) ti1.df ti2.df in
 238   let d1_bounds = map (fun i -> eval_m_taylor_partial_bound p_second i ti1) ns in
 239   let d2_bounds = map (fun i -> eval_m_taylor_partial_bound p_second i ti2) ns in
 240   let f1_bound = eval_m_taylor_bound p_second ti1 in
 241   let f2_bound = eval_m_taylor_bound p_second ti2 in
 242   let ddf =
 243     let ( * ), ( + ) = mul_interval p_second, add_interval p_second in
 244       map2 (fun (list1, list2) i ->
 245               let di1 = List.nth d1_bounds (i - 1) in
 246               let di2 = List.nth d2_bounds (i - 1) in
 247                 map2 (fun (dd1, dd2) j ->
 248                         let dj1 = List.nth d1_bounds (j - 1) in
 249                         let dj2 = List.nth d2_bounds (j - 1) in
 250                           (dd1 * f2_bound + di1 * dj2) + (dj1 * di2 + f1_bound * dd2))
 251                   (zip list1 list2) (1--i)) (zip ti1.ddf ti2.ddf) ns in
 252     {
 253       n = n;
 254       domain = ti1.domain;
 255       f = bounds;
 256       df = df;
 257       ddf = ddf;
 258     };;
 259
 260 (* neg *)
 261 let eval_m_taylor_neg taylor1 =
 262   let neg = neg_interval in
 263     {
 264       n = taylor1.n;
 265       domain = taylor1.domain;
 266       f = neg taylor1.f;
 267       df = map neg taylor1.df;
 268       ddf = map (map neg) taylor1.ddf;
 269     };;
 270
 271 (* inv *)
 272 let eval_m_taylor_inv p_lin p_second ti =
 273   let n = ti.n in
 274   let ns = 1--n in
 275   let f1_bound = eval_m_taylor_bound p_second ti in
 276   let bounds = inv_interval p_lin ti.f in
 277   let u_bounds =
 278     let neg, inv, ( * ) = neg_interval, inv_interval p_lin, mul_interval p_lin in
 279       neg (inv (ti.f * ti.f)) in
 280   let df =
 281     let ( * ) = mul_interval p_lin in
 282       map (fun d -> u_bounds * d) ti.df in
 283   let d1_bounds = map (fun i -> eval_m_taylor_partial_bound p_second i ti) ns in
 284   let d1, d2 =
 285     let inv, ( * ) = inv_interval p_second, mul_interval p_second in
 286     let ff = f1_bound * f1_bound in
 287       inv ff, two_interval * inv (f1_bound * ff) in
 288   let ddf =
 289     let ( * ), ( - ) = mul_interval p_second, sub_interval p_second in
 290       map2 (fun dd_list di1 ->
 291               Arith_misc.my_map2 (fun dd dj1 ->
 292                                     (d2 * dj1) * di1 - d1 * dd) dd_list d1_bounds) ti.ddf d1_bounds in
 293     {
 294       n = n;
 295       domain = ti.domain;
 296       f = bounds;
 297       df = df;
 298       ddf = ddf;
 299     };;
 300
 301
 302 (* sqrt *)
 303 let eval_m_taylor_sqrt p_lin p_second ti =
 304   let n = ti.n in
 305   let ns = 1--n in
 306   let f1_bound = eval_m_taylor_bound p_second ti in
 307   let bounds = sqrt_interval p_lin ti.f in
 308   let u_bounds =
 309     let inv, ( * ) = inv_interval p_lin, mul_interval p_lin in
 310       inv (two_interval * bounds) in
 311   let df =
 312     let ( * ) = mul_interval p_lin in
 313       map (fun d -> u_bounds * d) ti.df in
 314   let d1_bounds = map (fun i -> eval_m_taylor_partial_bound p_second i ti) ns in
 315   let d1, d2 =
 316     let neg, sqrt, inv, ( * ) = neg_interval, sqrt_interval p_second,
 317       inv_interval p_second, mul_interval p_second in
 318     let two_sqrt_f = two_interval * sqrt f1_bound in
 319       inv two_sqrt_f, neg (inv (two_sqrt_f * (two_interval * f1_bound))) in
 320   let ddf =
 321     let ( * ), ( + ) = mul_interval p_second, add_interval p_second in
 322       map2 (fun dd_list di1 ->
 323               Arith_misc.my_map2 (fun dd dj1 ->
 324                                     (d2 * dj1) * di1 + d1 * dd) dd_list d1_bounds) ti.ddf d1_bounds in
 325     {
 326       n = n;
 327       domain = ti.domain;
 328       f = bounds;
 329       df = df;
 330       ddf = ddf;
 331     };;
 332
 333
 334 (* atn *)
 335 let eval_m_taylor_atn =
 336   let neg_two_interval = neg_interval two_interval in
 337     fun p_lin p_second ti ->
 338       let n = ti.n in
 339       let ns = 1--n in
 340       let f1_bound = eval_m_taylor_bound p_second ti in
 341       let bounds = atn_interval p_lin ti.f in
 342       let u_bounds =
 343         let inv, ( + ), ( * ) = inv_interval p_lin, add_interval p_lin, mul_interval p_lin in
 344           inv (one_interval + ti.f * ti.f) in
 345       let df =
 346         let ( * ) = mul_interval p_lin in
 347           map (fun d -> u_bounds * d) ti.df in
 348       let d1_bounds = map (fun i -> eval_m_taylor_partial_bound p_second i ti) ns in
 349       let d1, d2 =
 350         let neg, inv, ( + ), ( * ) = neg_interval, inv_interval p_second,
 351           add_interval p_second, mul_interval p_second in
 352         let pow2 = pow_interval p_second 2 in
 353         let inv_one_ff = inv (one_interval + f1_bound * f1_bound) in
 354           inv_one_ff, (neg_two_interval * f1_bound) * pow2 inv_one_ff in
 355       let ddf =
 356         let ( * ), ( + ) = mul_interval p_second, add_interval p_second in
 357           map2 (fun dd_list di1 ->
 358                   Arith_misc.my_map2 (fun dd dj1 ->
 359                                         (d2 * dj1) * di1 + d1 * dd) dd_list d1_bounds) ti.ddf d1_bounds in
 360         {
 361           n = n;
 362           domain = ti.domain;
 363           f = bounds;
 364           df = df;
 365           ddf = ddf;
 366         };;
 367
 368
 369 (* acs *)
 370 let eval_m_taylor_acs p_lin p_second ti =
 371   let n = ti.n in
 372   let ns = 1--n in
 373   let f1_bound = eval_m_taylor_bound p_second ti in
 374   let bounds = acs_interval p_lin ti.f in
 375   let u_bounds =
 376     let inv, sqrt, neg = inv_interval p_lin, sqrt_interval p_lin, neg_interval in
 377     let ( * ), ( - ) = mul_interval p_lin, sub_interval p_lin in
 378       neg (inv (sqrt (one_interval - ti.f * ti.f))) in
 379   let df =
 380     let ( * ) = mul_interval p_lin in
 381       map (fun d -> u_bounds * d) ti.df in
 382   let d1_bounds = map (fun i -> eval_m_taylor_partial_bound p_second i ti) ns in
 383   let d1, d2 =
 384     let neg, sqrt, inv = neg_interval, sqrt_interval p_second, inv_interval p_second in
 385     let ( - ), ( * ), ( / ) = sub_interval p_second, mul_interval p_second, div_interval p_second in
 386     let pow3 = pow_interval p_second 3 in
 387     let ff_1 = one_interval - f1_bound * f1_bound in
 388       inv (sqrt ff_1), neg (f1_bound / sqrt (pow3 ff_1)) in
 389   let ddf =
 390     let ( * ), ( - ) = mul_interval p_second, sub_interval p_second in
 391       map2 (fun dd_list di1 ->
 392               Arith_misc.my_map2 (fun dd dj1 ->
 393                                     (d2 * dj1) * di1 - d1 * dd) dd_list d1_bounds) ti.ddf d1_bounds in
 394     {
 395       n = n;
 396       domain = ti.domain;
 397       f = bounds;
 398       df = df;
 399       ddf = ddf;
 400     };;
 401
 402
 403 end;;