mir.interpolate.polynomial source code

1 /++
2 $(H2 Lagrange Barycentric Interpolation)
3 
4 See_also: $(REF_ALTTEXT $(TT interp1), interp1, mir, interpolate)
5 
6 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
7 Copyright: 2020 Ilya Yaroshenko, Kaleidic Associates Advisory Limited, Symmetry Investments
8 Authors: Ilya Yaroshenko
9 
10 Macros:
11 SUBREF = $(REF_ALTTEXT $(TT $2), $2, mir, interpolate, $1)$(NBSP)
12 T2=$(TR $(TDNW $(LREF $1)) $(TD $+))
13 +/
14 module mir.interpolate.polynomial;
15 
16 public import mir.interpolate: atInterval;
17 import core.lifetime: move;
18 import mir.functional: RefTuple;
19 import mir.internal.utility : isFloatingPoint, Iota;
20 import mir.interpolate: findInterval;
21 import mir.math.common;
22 import mir.math.numeric;
23 import mir.math.sum;
24 import mir.ndslice.slice;
25 import mir.ndslice.topology: diff, map, member, iota, triplets;
26 import mir.rc.array;
27 import mir.utility: min, swap;
28 import std.traits: Unqual;
29 
30 @optmath:
31 
32 ///
33 version(mir_test)
34 @safe pure unittest
35 {
36     import mir.algorithm.iteration: all;
37     import mir.math;
38     import mir.ndslice;
39 
40     auto f(double x) { return (x-2) * (x-5) * (x-9); }
41     auto fd(double x) { return (x-5) * (x-9) + (x-2) * (x-5) + (x-2) * (x-9); }
42     auto fd2(double x) { return (x-5) + (x-9) + (x-2) + (x-5) + (x-2) + (x-9); }
43     double[2] fWithDerivative(double x) { return [f(x), fd(x)]; }
44     double[3] fWithTwoDerivatives(double x) { return [f(x), fd(x), fd2(x)]; }
45 
46     // 
47     auto x = [0.0, 3, 5, 10];
48     auto y = x.map!f.slice.field;
49     // `!2` adds first two derivatives: f, f', and f''
50     // default parameter is 0
51     auto l = x.lagrange!2(y);
52 
53     foreach(test; x ~ [2.0, 5, 9] ~ [double(PI), double(E), 10, 100, 1e3])
54     foreach(scale; [-1, +1, 1 + double.epsilon, 1 + double.epsilon.sqrt])
55     foreach(shift; [0, double.min_normal/2, double.min_normal*2, double.epsilon, double.epsilon.sqrt])
56     {
57         auto testX = test * scale + shift;
58 
59         auto lx = l(testX);
60         auto l_ldx = l.withDerivative(testX);
61         auto l_ld_ld2x = l.withTwoDerivatives(testX);
62 
63         assert(lx.approxEqual(f(testX)));
64         assert(l_ldx[].all!approxEqual(fWithDerivative(testX)[]));
65         assert(l_ld_ld2x[].all!approxEqual(fWithTwoDerivatives(testX)[]));
66     }
67 }
68 
69 /++
70 Constructs barycentric lagrange interpolant.
71 +/
72 Lagrange!(T, maxDerivative) lagrange(uint maxDerivative = 0, T, X)(scope const X[] x, scope const T[] y)
73     if (maxDerivative < 16)
74 {
75     import mir.ndslice.allocation: rcslice;
76     return x.rcslice!(immutable X).lagrange!maxDerivative(y.sliced);
77 }
78 
79 /// ditto
80 Lagrange!(Unqual!(Slice!(Iterator, 1, kind).DeepElement), maxDerivative, X)
81     lagrange(uint maxDerivative = 0, X, Iterator, SliceKind kind)(Slice!(RCI!(immutable X)) x, Slice!(Iterator, 1, kind) y) @trusted
82     if (maxDerivative < 16)
83 {
84     alias T = Unqual!(Slice!(Iterator, 1, kind).DeepElement);
85     return Lagrange!(T, maxDerivative)(x.move, rcarray!T(y));
86 }
87 
88 /++
89 +/
90 struct Lagrange(T, uint maxAdditionalFunctions = 0, X = T)
91     if (isFloatingPoint!T && maxAdditionalFunctions < 16)
92 {
93     /// $(RED for internal use only.)
94     Slice!(RCI!(immutable X)) _grid;
95     /// $(RED for internal use only.)
96     RCArray!(immutable T) _inversedBarycentricWeights;
97     /// $(RED for internal use only.)
98     RCArray!T[maxAdditionalFunctions + 1] _normalizedValues;
99     /// $(RED for internal use only.)
100     T[maxAdditionalFunctions + 1] _asums;
101 
102 @optmath @safe pure @nogc extern(D):
103 
104     ///
105     enum uint derivativeOrder = maxAdditionalFunctions;
106 
107     /++
108     Complexity: `O(N)`
109     +/
110     pragma(inline, false)
111     this(Slice!(RCI!(immutable X)) grid, RCArray!T values, RCArray!(immutable T) inversedBarycentricWeights)
112     {
113         import mir.algorithm.iteration: all;
114         assert(grid.length > 1);
115         assert(grid.length == values.length);
116         assert(grid.length == inversedBarycentricWeights.length);
117         enum msg = "Lagrange: grid points are too close or have not finite values.";
118         version(D_Exceptions)
119         {
120             enum T md = T.min_normal / T.epsilon;
121             static immutable exc = new Exception(msg);
122             if (!grid.lightScope.diff.all!(a => md <= a && a <= T.max))
123                 throw exc;
124         }
125         swap(_grid, grid);
126         swap(_inversedBarycentricWeights, inversedBarycentricWeights);
127         swap(_normalizedValues[0], values);
128         auto w = _inversedBarycentricWeights[].sliced;
129         foreach (_; Iota!maxAdditionalFunctions)
130         {
131             auto y = _normalizedValues[_][].sliced;
132             static if (_ == 0)
133                 _asums[_] = y.asumNorm;
134             _normalizedValues[_ + 1] = RCArray!T(_grid.length, true, false);
135             auto d = _normalizedValues[_ + 1][].sliced;
136             polynomialDerivativeValues(d, _grid.lightScope, y, w);
137             _asums[_ + 1] = d.asumNorm * _asums[_];
138         }
139     }
140 
141     /++
142     Complexity: `O(N^^2)`
143     +/
144     this(Slice!(RCI!(immutable X)) grid, RCArray!T values)
145     {
146         import core.lifetime: move; 
147         auto weights = grid.lightScope.inversedBarycentricWeights;
148         this(grid.move, values.move, weights.move);
149     }
150 
151 scope const:
152 
153     ///
154     Lagrange lightConst()() const @property @trusted { return *cast(Lagrange*)&this; }
155     ///
156     Lagrange lightImmutable()() immutable @property @trusted { return *cast(Lagrange*)&this; }
157 
158     @property
159     {
160         ///
161         ref const(Slice!(RCI!(immutable X))) grid() { return _grid; }
162         ///
163         immutable(X)[] gridScopeView() scope return const @property @trusted { return _grid.lightScope.field; }
164         ///
165         ref const(RCArray!(immutable T)) inversedBarycentricWeights() { return _inversedBarycentricWeights; }
166         ///
167         ref const(RCArray!T)[maxAdditionalFunctions + 1] normalizedValues() { return _normalizedValues; }
168         ///
169         ref const(T)[maxAdditionalFunctions + 1] asums() { return _asums; }
170         ///
171         size_t intervalCount(size_t dimension = 0)()
172         {
173             assert(_grid.length > 1);
174             return _grid.length - 1;
175         }
176     }
177 
178     template opCall(uint derivative = 0)
179         if (derivative <= maxAdditionalFunctions)
180     {
181         ///
182         pragma(inline, false)
183         auto opCall(T x)
184         {
185             return opCall!derivative(RefTuple!(size_t, T)(this.findInterval(x), x));
186         }
187 
188         ///
189         pragma(inline, false)
190         auto opCall(RefTuple!(size_t, T) tuple)
191         {
192 
193             const x = tuple[1];
194             const idx = tuple[0];
195             const inversedBarycentricWeights = _inversedBarycentricWeights[].sliced;
196             Slice!(const(T)*)[derivative + 1] values;
197             foreach (i; Iota!(derivative + 1))
198                 values[i] = _normalizedValues[i][].sliced;
199             const T[2] pd = [
200                 T(x - _grid[idx + 0]).fabs,
201                 T(x - _grid[idx + 1]).fabs,
202             ];
203             ptrdiff_t fastIndex =
204                 pd[0] < T.min_normal ? idx + 0 :
205                 pd[1] < T.min_normal ? idx + 1 :
206                 -1;
207             if (fastIndex >= 0)
208             {
209                 static if (derivative == 0)
210                 {
211                     return values[0][fastIndex] * _asums[0];
212                 }
213                 else
214                 {
215                     T[derivative + 1] ret = 0;
216                     foreach (_; Iota!(derivative + 1))
217                         ret[_] = values[_][fastIndex] * _asums[_];
218                     return ret;
219                 }
220             }
221             T d = 0;
222             T[derivative + 1] n = 0;
223             foreach (i; 0 .. _grid.length)
224             {
225                 T w = 1 / (inversedBarycentricWeights[i] * (x - _grid[i]));
226                 d += w;
227                 foreach (_; Iota!(derivative + 1))
228                     n[_] += w * values[_][i];
229             }
230             foreach (_; Iota!(derivative + 1))
231             {
232                 n[_] /= d;
233                 n[_] *= asums[_];
234             }
235             static if (derivative == 0)
236                 return n[0];
237             else
238                 return n;
239         }
240     }
241 
242     static if (maxAdditionalFunctions)
243     {
244         ///
245         alias withDerivative = opCall!1;
246 
247         static if (maxAdditionalFunctions > 1)
248         {
249             ///
250             alias withTwoDerivatives = opCall!2;
251         }
252     }
253 }
254 
255 
256 /++
257 +/
258 pragma(inline, false)
259 @nogc
260 RCArray!(immutable T) inversedBarycentricWeights(T)(Slice!(const(T)*) x)
261     if (isFloatingPoint!T)
262 {
263     
264     const n = x.length;
265     scope prodsa = RCArray!(ProdAccumulator!T)(n);
266     scope p = prodsa[].sliced;
267     foreach (triplet; n.iota.triplets) with(triplet)
268     {
269         foreach (l; left)
270         {
271             auto e = prod!T(x[center] - x[l]);
272             p[l] *= -e;
273             p[center] *= e;
274         }
275     }
276     import mir.algorithm.iteration: reduce;
277     const minExp = long.max.reduce!min(p.member!"exp");
278     foreach (ref e; p)
279         e = e.ldexp(1 - minExp);
280     auto ret = rcarray!(immutable T)(p.member!"prod");
281     return ret;
282 }
283 
284 /++
285 Computes derivative values in the same points
286 Params:
287     d = derivative values (output)
288     y = values
289     x = grid
290     w = inversed barycentric weights
291 Returns:
292     Derivative values in the same points
293 +/
294 @nogc
295 Slice!(T*) polynomialDerivativeValues(T)(return Slice!(T*) d,  Slice!(const(T)*) x,  Slice!(const(T)*) y,  Slice!(const(T)*) w)
296     if (isFloatingPoint!T)
297 {
298     pragma(inline, false);
299     import mir.math.sum;
300 
301     const n = x.length;
302     assert(n == w.length);
303     assert(n == y.length);
304 
305     d.fillCollapseSums!(Summation.fast,
306         (c, s1, s2) => w[c] * s1 + y[c] * s2,
307         (c, _) => y[_] / (w[_] * (x[c] - x[_])),
308         (c, _) => map!"1 / a"(x[c] - x[_]));
309     return d;
310 }
311 
312 ///
313 @nogc
314 Slice!(T*) polynomialDerivativeValues(T)(return Slice!(T*) d,  Slice!(const(T)*) x,  Slice!(const(T)*) y)
315     if (isFloatingPoint!T)
316 {
317     return .polynomialDerivativeValues(d, x, y, x.inversedBarycentricWeights[].sliced);
318 }
319 
320 private T asumNorm(T)(Slice!(T*) v)
321 {
322     pragma(inline, false);
323     auto n = v.map!fabs.sum!"fast";
324     v[] /= n;
325     return n;
326 }