1 | export class MonotonicInterpolant {
|
2 | constructor(xs, ys) {
|
3 | const { length } = xs;
|
4 |
|
5 | const indexes = [];
|
6 | for (let i = 0; i < length; i++) {
|
7 | indexes.push(i);
|
8 | }
|
9 | indexes.sort((a, b) => (xs[a] < xs[b] ? -1 : 1));
|
10 |
|
11 | const dys = [];
|
12 | const dxs = [];
|
13 | const ms = [];
|
14 | let dx;
|
15 | let dy;
|
16 | for (let i = 0; i < length - 1; i++) {
|
17 | dx = xs[i + 1] - xs[i];
|
18 | dy = ys[i + 1] - ys[i];
|
19 | dxs.push(dx);
|
20 | dys.push(dy);
|
21 | ms.push(dy / dx);
|
22 | }
|
23 |
|
24 | const c1s = [ms[0]];
|
25 | for (let i = 0; i < dxs.length - 1; i++) {
|
26 | const m2 = ms[i];
|
27 | const mNext = ms[i + 1];
|
28 | if (m2 * mNext <= 0) {
|
29 | c1s.push(0);
|
30 | }
|
31 | else {
|
32 | dx = dxs[i];
|
33 | const dxNext = dxs[i + 1];
|
34 | const common = dx + dxNext;
|
35 | c1s.push((3 * common) / ((common + dxNext) / m2 + (common + dx) / mNext));
|
36 | }
|
37 | }
|
38 | c1s.push(ms[ms.length - 1]);
|
39 |
|
40 | const c2s = [];
|
41 | const c3s = [];
|
42 | let m;
|
43 | for (let i = 0; i < c1s.length - 1; i++) {
|
44 | m = ms[i];
|
45 | const c1 = c1s[i];
|
46 | const invDx = 1 / dxs[i];
|
47 | const common = c1 + c1s[i + 1] - m - m;
|
48 | c2s.push((m - c1 - common) * invDx);
|
49 | c3s.push(common * invDx * invDx);
|
50 | }
|
51 | this.xs = xs;
|
52 | this.ys = ys;
|
53 | this.c1s = c1s;
|
54 | this.c2s = c2s;
|
55 | this.c3s = c3s;
|
56 | }
|
57 | interpolate(x) {
|
58 | const { xs, ys, c1s, c2s, c3s } = this;
|
59 |
|
60 | let i = xs.length - 1;
|
61 | if (x === xs[i]) {
|
62 | return ys[i];
|
63 | }
|
64 |
|
65 | let low = 0;
|
66 | let high = c3s.length - 1;
|
67 | let mid;
|
68 | while (low <= high) {
|
69 | mid = Math.floor(0.5 * (low + high));
|
70 | const xHere = xs[mid];
|
71 | if (xHere < x) {
|
72 | low = mid + 1;
|
73 | }
|
74 | else if (xHere > x) {
|
75 | high = mid - 1;
|
76 | }
|
77 | else {
|
78 | return ys[mid];
|
79 | }
|
80 | }
|
81 | i = Math.max(0, high);
|
82 |
|
83 | const diff = x - xs[i];
|
84 | const diffSq = diff * diff;
|
85 | return ys[i] + c1s[i] * diff + c2s[i] * diffSq + c3s[i] * diff * diffSq;
|
86 | }
|
87 | }
|