formal_lp/old/formal_interval/interval_1d/taylor.hl

   1 (* port of taylor functions, taylor interval *)
   2
   3 (*
   4 The first part of the file implements basic operations on type taylor_interval.
   5
   6 Then a type tfunction is defined that represents a twice continuously
   7 differentiable function of six variables.  It can be evaluated, which
   8 is the taylor_interval data associated with it.
   9
  10 Sometimes a tfunction f is used to represent an inequality f < 0.
  11 (See recurse.hl.
  12 *)
  13
  14
  15 module Taylor = struct
  16
  17 open Interval;;
  18 open Univariate;;
  19 open Line_interval;;
  20
  21
  22 (* general utilities *)
  23
  24
  25 let center_form(x,z) =
  26   let ( + ) , ( - ), ( / ) = up(); upadd,upsub,updiv in
  27   let y = if x = z then x else (x + z) / 2.0 in
  28   let w = max (z - y) (y - x) in
  29   let _ = (w >= 0.0) or failwith "centerform" in
  30      (y,w);;
  31
  32 (* start with taylor interval operations *)
  33
  34 let make_taylor_interval (l1,w1,dd1) = {l = l1; w = w1; dd=dd1;};;
  35
  36 let ti_add (ti1,ti2) =
  37   let _ = (ti1.w = ti2.w) or failwith ("width mismatch in ti") in
  38    make_taylor_interval( ladd ti1.l ti2.l,ti1.w, iadd ti1.dd ti2.dd);;
  39
  40 let ti_sub (ti1,ti2) =
  41         let _ = (ti1.w = ti2.w) or failwith ("width mistmatch in ti") in
  42                 make_taylor_interval( lsub ti1.l ti2.l, ti1.w, isub ti1.dd ti2.dd );;
  43
  44 let ti_scale (ti,t) =
  45    make_taylor_interval( smul ti.l t, ti.w, imul ti.dd t );;
  46
  47
  48
  49 let taylor_error ti =
  50   let ( + ), ( * ) , ( / )= up(); upadd, upmul, updiv in
  51         (ti.w * ti.w * iabs ti.dd) / 2.0;;
  52
  53 let upper_bound ti =
  54   let e = taylor_error ti in
  55   let ( + ), ( * ) = up(); upadd, upmul in
  56         ti.l.f.hi + ti.w * iabs ti.l.df + e;;
  57
  58 let lower_bound ti =
  59   let e = taylor_error ti in
  60   let ( + ), ( * ),( - ) = down(); downadd,downmul,downsub in
  61     ti.l.f.lo - e  + ( ~-. (ti.w)) * iabs ti.l.df;;
  62
  63 (*    t + List.fold_left2 (fun a b c -> a + ( ~-. b) * iabs c) 0.0 ti.w ti.l.df;; *)
  64
  65 let upper_partial ti =
  66   let ( + ), ( * ) =   up(); upadd,upmul in
  67         let err = ti.w * (max ti.dd.hi (~-. (ti.dd.lo))) in
  68       err + Interval.sup ( ti.l.df );;
  69
  70 let lower_partial ti =
  71   let ( + ), ( * ), ( - ) = down();downadd,downmul,downsub in
  72     let err = ti.w * min ti.dd.lo (~-. (ti.dd.hi)) in
  73       Interval.inf ( ti.l.df ) + err;;
  74
  75
  76
  77
  78 let ti_mul (ti1, ti2) =
  79   let _ = (ti1.w = ti2.w) or failwith ("width mistmatch in ti") in
  80   let f = mk_interval (lower_bound ti1, upper_bound ti1) and
  81       g = mk_interval (lower_bound ti2, upper_bound ti2) and
  82       df = mk_interval (lower_partial ti1, upper_partial ti1) and
  83       dg = mk_interval (lower_partial ti2, upper_partial ti2) and
  84       ( + ), ( * ) = iadd, imul in
  85     make_taylor_interval( lmul ti1.l ti2.l, ti1.w,
  86                           ti1.dd * g + two * df * dg + f * ti2.dd);;
  87
  88
  89 (* primitive A *)
  90
  91 type primitiveA = {
  92   hfn : float -> line;
  93   second : float -> float -> interval;
  94 };;
  95
  96 let make_primitiveA (h1,s1) = {hfn = h1; second = s1; };;
  97
  98 let unitA =
  99         make_primitiveA (
 100       (fun y -> line_unit),
 101       (fun x z -> zero)
 102 );;
 103
 104 let evalf4A pA w x y z =
 105   make_taylor_interval(
 106     pA.hfn y,
 107     w,
 108     pA.second x z
 109   );;
 110
 111 let line_estimateA pA y = pA.hfn y;;
 112
 113 (* primitive U *)
 114
 115 type primitiveU = {
 116   uv: univariate;
 117 };;
 118
 119 let mk_primitiveU uv1 = { uv = uv1; };;
 120
 121 let line_estimateU p y =
 122   let y0 = y in
 123   let t = mk_interval(y0,y0) in
 124   let d = eval p.uv t 1 in
 125         mk_line (eval p.uv t 0, d);;
 126
 127 let evalf4U p w x y z =
 128         let t = mk_interval(x,z) in
 129         make_taylor_interval(
 130                 line_estimateU p y,
 131                 w,
 132                 eval p.uv t 2
 133         );;
 134
 135 type tfunction =
 136   |  Prim_a of primitiveA
 137   |  Uni of primitiveU
 138   |  Plus of tfunction * tfunction
 139   |  Minus of tfunction * tfunction
 140   |  Product of tfunction * tfunction
 141   |  Scale of tfunction * interval
 142   |  Uni_compose of univariate * tfunction
 143   |  Composite of tfunction * tfunction;;
 144
 145 let unit = Prim_a unitA;;
 146
 147 let x1 = Uni (mk_primitiveU ux);;
 148
 149
 150 (* This is one of the most difficult functions in the interval code.
 151    It uses the chain rule to compute the second partial derivatives with
 152    respect to x(i) x(j), of a function composition
 153
 154    F(x1,...,x6) = f(g1(x1,...x6),g2(x1,...x6),...,g6(x1,...x6)).
 155
 156    (F i j) = sum {k m} (f k m) (gk i) (gm j)     + sum {r} (f r) (gr i j).
 157
 158    Fast performance of this function is very important, especially
 159    when many of the functions g* are constant.
 160    There is a bit of imperative programming here, in computing the sums.
 161
 162    Note that ( + ) and ( * ) have different types in various subsections.
 163 *)
 164
 165 let eval_composite =
 166   fun hdr p w ->
 167         (* wide and narrow ranges of p *)
 168     let (aw,bw) = lower_bound p, upper_bound p in
 169     let (a,b) = p.l.f.lo, p.l.f.hi in
 170         (* wide and narrow widths from a to b *)
 171     let (u,wu,wf) =
 172                 let ( + ),( - ),( / ) = up();upadd,upsub,updiv in
 173                 let u = (a + b) / 2.0 in
 174                 let wu = max (bw - u) (u - aw) in
 175                 let wf = max (b - u) (u - a) in
 176                         (u, wu, wf) in
 177         let (fu:taylor_interval) = hdr wu aw u bw in
 178     let fpy =
 179                 let t = make_taylor_interval(fu.l,wf,fu.dd) in
 180                         mk_line (
 181                                 mk_interval(lower_bound t, upper_bound t),
 182                                 mk_interval(lower_partial t, upper_partial t)) in
 183         (* use chain rule imperatively to compute narrow first derivative *)
 184         let lin =
 185                 let ( * ) = imul in
 186                         mk_line (fpy.f, fpy.df * p.l.df) in
 187         (* second derivative init *)
 188         let fW_partial = mk_interval(lower_partial fu, upper_partial fu) in
 189         let pW_partial = mk_interval(lower_partial p, upper_partial p) in
 190         let dcw_list =
 191                 let ( + ) = iadd in
 192                 let ( * ) = imul in
 193                         fW_partial * p.dd + fu.dd * pW_partial * pW_partial in
 194         make_taylor_interval(lin,w,dcw_list);;
 195
 196 let rec evalf4 tf w x y z = match tf with
 197   | Prim_a p -> evalf4A p w x y z
 198   | Uni p -> evalf4U p w x y z
 199   | Plus (tf1,tf2) -> ti_add(evalf4 tf1 w x y z, evalf4 tf2 w x y z)
 200   | Minus (tf1, tf2) -> ti_sub(evalf4 tf1 w x y z, evalf4 tf2 w x y z)
 201   | Product (tf1, tf2) -> ti_mul(evalf4 tf1 w x y z, evalf4 tf2 w x y z)
 202   | Composite(h,g) -> eval_composite (evalf4 h) (evalf4 g w x y z) w
 203   | Scale (tf,t) -> ti_scale ((evalf4 tf w x y z),t)
 204   | Uni_compose (uf,tf) ->
 205       evalf4 (Composite(Uni (mk_primitiveU uf),tf)) w x y z;;
 206
 207 let evalf tf x z =
 208   let (y,w) = center_form (x,z) in
 209     evalf4 tf w x y z;;
 210
 211 let line_estimate_composite =
 212   let ( * ) = imul in
 213     fun h p ->
 214       let (a,b) = (p.f.lo, p.f.hi) in
 215       let fN = evalf h a b in
 216       let fN_partial = mk_interval(lower_partial fN, upper_partial fN) in
 217                 mk_line (fN.l.f, fN_partial * p.df);;
 218
 219 let rec line_estimate tf y = match tf with
 220   | Prim_a p -> line_estimateA p y
 221   | Uni p -> line_estimateU p y
 222   | Plus (p,q) -> ladd (line_estimate p y) (line_estimate q y)
 223   | Minus (p,q) -> lsub (line_estimate p y) (line_estimate q y)
 224   | Product (p,q) -> lmul (line_estimate p y) (line_estimate q y)
 225   | Scale (p,t) -> smul (line_estimate p y) t
 226   | Uni_compose (uf,tf) ->
 227       line_estimate (Composite(Uni (mk_primitiveU uf),tf)) y
 228   | Composite(h,g) -> line_estimate_composite h (line_estimate g y);;
 229
 230 end;;