UNPKG

2.74 kBJavaScriptView Raw
1export default class MonotonicInterpolant {
2 constructor(xs, ys) {
3 const { length } = xs;
4 // Rearrange xs and ys so that xs is sorted
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 // Get consecutive differences and slopes
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 // Get degree-1 coefficients
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 // Get degree-2 and degree-3 coefficients
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 // The rightmost point in the dataset should give an exact result
60 let i = xs.length - 1;
61 if (x === xs[i]) {
62 return ys[i];
63 }
64 // Search for the interval x is in, returning the corresponding y if x is one of the original xs
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 // Interpolate
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}