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