Source file `owl_maths_interpolate.ml`

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
# 1 "src/base/maths/owl_maths_interpolate.ml"
(*
 * OWL - OCaml Scientific and Engineering Computing
 * Copyright (c) 2016-2020 Liang Wang <liang.wang@cl.cam.ac.uk>
 *)

(** Interpolation and Extrapolation *)

let polint xs ys x =
  let n = Array.length xs in
  let m = Array.length ys in
  let error () =
    let s =
      Printf.sprintf
        "polint requires that xs and ys have the same length, but xs length is %i \
         whereas ys length is %i"
        n
        m
    in
    Owl_exception.INVALID_ARGUMENT s
  in
  Owl_exception.verify (m = n) error;
  let c = Array.copy ys in
  let d = Array.copy ys in
  let ns = ref 0 in
  let dy = ref 0. in
  let dif = ref (abs_float (x -. xs.(0))) in
  for i = 0 to n - 1 do
    let dift = abs_float (x -. xs.(i)) in
    if dift < !dif
    then (
      ns := i;
      dif := dift)
  done;
  let y = ref ys.(!ns) in
  ns := !ns - 1;
  for m = 0 to n - 2 do
    for i = 0 to n - m - 2 do
      let ho = xs.(i) -. x in
      let hp = xs.(i + m + 1) -. x in
      let w = c.(i + 1) -. d.(i) in
      let den = ho -. hp in
      assert (den <> 0.);
      let den = w /. den in
      c.(i) <- ho *. den;
      d.(i) <- hp *. den
    done;
    if (2 * !ns) + 2 < n - m - 1
    then dy := c.(!ns + 1)
    else (
      dy := d.(!ns);
      ns := !ns - 1);
    y := !y +. !dy
  done;
  !y, !dy


(* TODO: not tested yet *)
let ratint xs ys x =
  let n = Array.length xs in
  let m = Array.length ys in
  let error () =
    let s =
      Printf.sprintf
        "ratint requires that xs and ys have the same length, but xs length is %i \
         whereas ys length is %i"
        n
        m
    in
    Owl_exception.INVALID_ARGUMENT s
  in
  Owl_exception.verify (m = n) error;
  let c = Array.copy ys in
  let d = Array.copy ys in
  let hh = ref (abs_float (x -. xs.(0))) in
  let y = ref 0. in
  let dy = ref 0. in
  let ns = ref 0 in
  let eps = 1e-25 in
  try
    for i = 0 to n do
      let h = abs_float (x -. xs.(i)) in
      if h = 0.
      then (
        y := ys.(i);
        dy := 0.;
        raise Owl_exception.FOUND)
      else if h < !hh
      then (
        ns := i;
        hh := h;
        c.(i) <- ys.(i);
        d.(i) <- ys.(i) +. eps)
    done;
    y := ys.(!ns);
    ns := !ns - 1;
    for m = 1 to n - 1 do
      for i = 1 to n - m do
        let w = c.(i) -. d.(i - 1) in
        let h = xs.(i + m - 1) -. x in
        let t = (xs.(i - 1) -. x) *. d.(i - 1) /. h in
        let dd = t -. c.(i) in
        if dd = 0. then failwith "Has a pole";
        let dd = w /. dd in
        d.(i - 1) <- c.(i) *. dd;
        c.(i - 1) <- t *. dd
      done
    done;
    !y, !dy
  with
  | Owl_exception.FOUND -> !y, !dy
  | e                   -> raise e