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
113
# 1 "src/base/maths/owl_maths_interpolate.ml"
(*
 * OWL - OCaml Scientific and Engineering Computing
 * Copyright (c) 2016-2019 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