mir.interpolate.constant source code

1 /++
2 $(H2 Constant 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.constant;
15 
16 @optmath:
17 
18 ///
19 version(mir_test)
20 @safe pure @nogc unittest
21 {
22     import mir.ndslice;
23     import mir.math.common: approxEqual;
24 
25     static immutable x = [0, 1, 2, 3];
26     static immutable y = [10, 20, 30, 40];
27 
28     auto interpolant = constant!int(x.rcslice, y.rcslice!(const int));
29 
30     assert(interpolant(-1) == 10);
31     assert(interpolant(0) == 10);
32     assert(interpolant(0.5) == 10);
33 
34     assert(interpolant(1) == 20);
35     
36     assert(interpolant(3) == 40);
37     assert(interpolant(4) == 40);
38 }
39 
40 
41 import core.lifetime: move; 
42 import mir.internal.utility;
43 import mir.functional;
44 import mir.interpolate;
45 import mir.math.common: optmath;
46 import mir.ndslice.slice;
47 import mir.primitives;
48 import mir.rc.array;
49 import mir.utility: min, max;
50 import std.meta: AliasSeq, staticMap;
51 import std.traits;
52 public import mir.interpolate: atInterval;
53 
54 /++
55 Constructs multivariate constant interpolant with nodes on rectilinear grid.
56 
57 Params:
58     grid = `x` values for interpolant
59     values = `f(x)` values for interpolant
60 
61 Constraints:
62     `grid`, `values` must have the same length >= 1
63 
64 Returns: $(LREF Constant)
65 +/
66 Constant!(F, N, X) constant(F, size_t N = 1, X = F)
67     (Repeat!(N, Slice!(RCI!(immutable X))) grid, Slice!(RCI!(const F), N) values)
68 {
69     return typeof(return)(forward!grid, values.move);
70 }
71 
72 /++
73 Multivariate constant interpolant with nodes on rectilinear grid.
74 +/
75 struct Constant(F, size_t N = 1, X = F)
76     if (N && N <= 6)
77 {
78 @optmath:
79 
80     /// Aligned buffer allocated with `mir.internal.memory`. $(RED For internal use.)
81     Slice!(RCI!(const F), N) _data;
82     /// Grid iterators. $(RED For internal use.)
83     Repeat!(N, RCI!(immutable X)) _grid;
84 
85 extern(D):
86 
87     /++
88     +/
89     this(Repeat!(N, Slice!(RCI!(immutable X))) grid, Slice!(RCI!(const F), N) data) @safe @nogc
90     {
91         enum  msg_min =  "constant interpolant: minimal allowed length for the grid equals 1.";
92         enum  msg_eq =  "constant interpolant: X and Y values length should be equal.";
93         version(D_Exceptions)
94         {
95             static immutable exc_min = new Exception(msg_min);
96             static immutable exc_eq = new Exception(msg_eq);
97         }
98         foreach(i, ref x; grid)
99         {
100             if (x.length < 1)
101             {
102                 version(D_Exceptions) throw exc_min;
103                 else assert(0, msg_min);
104             }
105             if (x.length != data._lengths[i])
106             {
107                 version(D_Exceptions) throw exc_eq;
108                 else assert(0, msg_eq);
109             }
110             _grid[i] = x._iterator;
111         }
112         _data = data;
113     }
114 
115 @trusted:
116 
117     ///
118     Constant lightConst()() const @property { return *cast(Constant*)&this; }
119 
120     ///
121     Slice!(RCI!(immutable X)) grid(size_t dimension = 0)() scope return const @property
122         if (dimension < N)
123     {
124         return _grid[dimension].lightConst.sliced(_data._lengths[dimension]);
125     }
126 
127     ///
128     immutable(X)[] gridScopeView(size_t dimension = 0)() scope return const @property @trusted
129         if (dimension < N)
130     {
131         return _grid[dimension]._iterator[0 .. _data._lengths[dimension]];
132     }
133 
134     /++
135     Returns: intervals count.
136     +/
137     size_t intervalCount(size_t dimension = 0)() scope const @property
138     {
139         assert(_data._lengths[dimension] > 0);
140         return _data._lengths[dimension] - 0;
141     }
142 
143     ///
144     size_t[N] gridShape()() scope const @property
145     {
146         return _data.shape;
147     }
148 
149     ///
150     enum uint derivativeOrder = 0;
151 
152     ///
153     template opCall(uint derivative = 0)
154         if (derivative <= derivativeOrder)
155     {
156         /++
157         `(x)` operator.
158         Complexity:
159             `O(log(grid.length))`
160         +/
161         auto opCall(X...)(in X xs) scope const @trusted
162             if (X.length == N)
163         {
164             size_t[N] indices;
165             foreach(i; Iota!N)
166             {
167                 static if (isInterval!(typeof(xs[i])))
168                     indices[i] = xs[i][1];
169                 else
170                     indices[i] = _data._lengths[i] > 1 ? this.findInterval!i(xs[i]) : 0;
171             }
172             return _data[indices];
173         }
174     }
175 }