Skip to content

Commit 71a25e9

Browse files
9il9il
committed
add mir.interpolation
1 parent 58f4d5d commit 71a25e9

File tree

3 files changed

+568
-0
lines changed

3 files changed

+568
-0
lines changed

‎source/mir/interpolation/linear.d

Lines changed: 156 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,156 @@
1+
/++
2+
$(H2 Linear Interpolation)
3+
4+
See_also: $(REF_ALTTEXT $(TT interp1), interp1, mir, interpolation)
5+
6+
License: $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0).
7+
Copyright: Copyright © 2017-, Ilya Yaroshenko
8+
Authors: Ilya Yaroshenko
9+
10+
Macros:
11+
SUBREF = $(REF_ALTTEXT $(TT $2), $2, mir, interpolation, $1)$(NBSP)
12+
T2=$(TR $(TDNW $(LREF $1)) $(TD $+))
13+
+/
14+
module mir.interpolation.linear;
15+
16+
import std.traits;
17+
import mir.array.primitives;
18+
import mir.utility: fastmath;
19+
20+
@fastmath:
21+
22+
/++
23+
Unbounded linear interpolation.
24+
+/
25+
struct LinearInterpolation(RangeG, RangeV)
26+
{
27+
///
28+
RangeG _grid;
29+
///
30+
RangeV _values;
31+
32+
private alias G = Unqual!(ForeachType!RangeG);
33+
private alias V = Unqual!(ForeachType!RangeV);
34+
35+
@fastmath:
36+
37+
this()(RangeG grid, RangeV values)
38+
{
39+
assert(grid.length >= 2);
40+
assert(grid.length == values.length);
41+
this._grid = grid;
42+
this._values = values;
43+
}
44+
45+
/++
46+
Interval index for x.
47+
+/
48+
size_t interval(T)(in T x)
49+
{
50+
import std.range: assumeSorted;
51+
return _grid[1 .. $ - 1]
52+
.assumeSorted
53+
.lowerBound(x)
54+
.length;
55+
}
56+
57+
/++
58+
`(x)` and `[x]` operators.
59+
Complexity:
60+
`O(log(_grid.length))`
61+
+/
62+
auto opCall(T)(in T x)
63+
{
64+
return opCall!T(x, interval(x));
65+
}
66+
67+
/++
68+
`(x, interval)` and `[x, interval]` operators.
69+
Complexity:
70+
`O(1)`
71+
+/
72+
auto opCall(T)(in T x, size_t interval)
73+
{
74+
assert(interval + 1 < _grid.length);
75+
76+
auto x0 = _grid [interval + 0];
77+
auto x1 = _grid [interval + 1];
78+
auto y0 = _values[interval + 0];
79+
auto y1 = _values[interval + 1];
80+
81+
return eval(x0, x1, y0, y1, x);
82+
}
83+
84+
/// ditto
85+
alias opIndex = opCall;
86+
87+
package:
88+
static alias eval = .evalImpl!(G, V);
89+
90+
size_t length ()() @property { return grid.length - 1; }
91+
bool empty ()() @property { return length == 0; }
92+
auto ref front()() @property { return grid.front; }
93+
94+
void popFront()()
95+
{
96+
assert(grid.length > 1);
97+
grid.popFront;
98+
values.popFront;
99+
}
100+
}
101+
102+
private template evalImpl(G, V)
103+
{
104+
@fastmath
105+
auto evalImpl(T)(G x0, G x1, V y0, V y1, in T x)
106+
{
107+
auto step = x1 - x0;
108+
auto w0 = x - x0;
109+
auto w1 = x1 - x;
110+
w0 /= step;
111+
w1 /= step;
112+
y0 *= w1;
113+
y1 *= w0;
114+
return y0 + y1;
115+
}
116+
}
117+
118+
119+
/++
120+
Linear interpolation.
121+
122+
Params:
123+
grid = `x` values for interpolation
124+
values = `f(x)` values for interpolation
125+
126+
Constraints:
127+
`grid`, `values` must have the same length >= 3
128+
129+
Returns: $(LREF LinearInterpolation)
130+
+/
131+
LinearInterpolation!(RangeG, RangeV) linearInterpolation(RangeG, RangeV)(RangeG grid, RangeV values)
132+
{
133+
return typeof(return)(grid, values);
134+
}
135+
136+
///
137+
unittest
138+
{
139+
140+
}
141+
142+
unittest
143+
{
144+
import mir.ndslice;
145+
import std.math: approxEqual;
146+
147+
auto x = [0, 1, 2, 3, 5.00274, 7.00274, 10.0055, 20.0137, 30.0192];
148+
auto y = [0.0011, 0.0011, 0.0030, 0.0064, 0.0144, 0.0207, 0.0261, 0.0329, 0.0356,];
149+
auto xs = [1, 2, 3, 4.00274, 5.00274, 6.00274, 7.00274, 8.00548, 9.00548, 10.0055, 11.0055, 12.0082, 13.0082, 14.0082, 15.0082, 16.011, 17.011, 18.011, 19.011, 20.0137, 21.0137, 22.0137, 23.0137, 24.0164, 25.0164, 26.0164, 27.0164, 28.0192, 29.0192, 30.0192];
150+
151+
auto interpolation = linearInterpolation(x, y);
152+
153+
auto data = [0.0011, 0.0030, 0.0064, 0.0104, 0.0144, 0.0176, 0.0207, 0.0225, 0.0243, 0.0261, 0.0268, 0.0274, 0.0281, 0.0288, 0.0295, 0.0302, 0.0309, 0.0316, 0.0322, 0.0329, 0.0332, 0.0335, 0.0337, 0.0340, 0.0342, 0.0345, 0.0348, 0.0350, 0.0353, 0.0356];
154+
155+
assert(approxEqual(xs.sliced.map!interpolation, data, 1e-4, 1e-4));
156+
}

‎source/mir/interpolation/package.d

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
/++
2+
$(H1 Interpolation Algorithms)
3+
4+
$(BOOKTABLE $(H2 Interpolation modules),
5+
$(TR $(TH Module) $(TH Interpolation kind))
6+
$(T2M linear, Linear Interpolation)
7+
$(T2M pchip, Cubic Hermite Interpolation)
8+
)
9+
10+
Copyright: Copyright © 2017-, Ilya Yaroshenko
11+
Authors: Ilya Yaroshenko
12+
13+
Macros:
14+
SUBREF = $(REF_ALTTEXT $(TT $2), $2, mir, interpolation, $1)$(NBSP)
15+
T2M=$(TR $(TDNW $(MREF mir,interpolation,$1)) $(TD $+))
16+
T2=$(TR $(TDNW $(LREF $1)) $(TD $+))
17+
+/
18+
module mir.interpolation;
19+
20+
import mir.primitives;
21+
22+
/++
23+
Lazy interpolation shell with linear complexity.
24+
25+
Params:
26+
range = sorted range
27+
interpolation = interpolation structure with `._grid` and `.opCall(x, interval)` methods.
28+
Complexity:
29+
`O(range.length + interpolation._grid.length)` to evaluate all elements.
30+
Returns:
31+
Lazy input range.
32+
See_also:
33+
$(SUBREF linear, linearInterpolation),
34+
$(SUBREF pchip, pchip).
35+
+/
36+
auto interp1(Range, Interpolation)(Range range, Interpolation interpolation, size_t interval = 0)
37+
{
38+
return Interp1!(Range, Interpolation)(range, interpolation, interval);
39+
}
40+
41+
/// ditto
42+
struct Interp1(Range, Interpolation)
43+
{
44+
/// Sorted range (descending)
45+
Range _range;
46+
/// Interpolation structure
47+
Interpolation _interpolation;
48+
/// Current interpolation interval
49+
size_t _interval;
50+
51+
static if (hasLength!Range)
52+
/// Length (optional)
53+
size_t length() @property { return _range.length; }
54+
/// Input range primitives
55+
bool empty () @property { return _range.empty; }
56+
/// ditto
57+
void popFront() { _range.popFront; }
58+
/// ditto
59+
auto front() @property
60+
61+
{
62+
assert(!empty);
63+
auto x = _range.front;
64+
while (x > _interpolation._grid[_interval + 1] && _interpolation._grid.length > _interval + 2)
65+
_interval++;
66+
return _interpolation(x, _interval);
67+
}
68+
}
69+
70+
/++
71+
PCHIP interpolation.
72+
73+
Complexity:
74+
`O(x.length + xs.length)`
75+
76+
See_also:
77+
$(MREF mir,_interpolation,pchip)
78+
+/
79+
unittest
80+
{
81+
import std.math: approxEqual;
82+
import mir.interpolation: interp1;
83+
import mir.interpolation.pchip;
84+
85+
auto x = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22];
86+
auto y = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6];
87+
auto interpolation = x.pchip(y);
88+
89+
auto xs = x[0 .. $ - 1].dup;
90+
xs[] += 0.5;
91+
92+
auto ys = xs.interp1(interpolation);
93+
94+
import std.stdio;
95+
writeln(ys);
96+
97+
assert(ys.approxEqual([
98+
5.333333333333334,
99+
2.500000000000000,
100+
10.000000000000000,
101+
4.288971807628524,
102+
11.202580845771145,
103+
16.250000000000000,
104+
17.962962962962962,
105+
5.558593750000000,
106+
17.604662698412699,
107+
]));
108+
109+
}
110+
111+
unittest
112+
{
113+
import mir.interpolation.linear;
114+
import mir.ndslice;
115+
import std.math: approxEqual;
116+
117+
auto x = [0, 1, 2, 3, 5.00274, 7.00274, 10.0055, 20.0137, 30.0192];
118+
auto y = [0.0011, 0.0011, 0.0030, 0.0064, 0.0144, 0.0207, 0.0261, 0.0329, 0.0356,];
119+
auto xs = [1, 2, 3, 4.00274, 5.00274, 6.00274, 7.00274, 8.00548, 9.00548, 10.0055, 11.0055, 12.0082, 13.0082, 14.0082, 15.0082, 16.011, 17.011, 18.011, 19.011, 20.0137, 21.0137, 22.0137, 23.0137, 24.0164, 25.0164, 26.0164, 27.0164, 28.0192, 29.0192, 30.0192];
120+
121+
auto interpolation = linearInterpolation(x, y);
122+
123+
auto data = [0.0011, 0.0030, 0.0064, 0.0104, 0.0144, 0.0176, 0.0207, 0.0225, 0.0243, 0.0261, 0.0268, 0.0274, 0.0281, 0.0288, 0.0295, 0.0302, 0.0309, 0.0316, 0.0322, 0.0329, 0.0332, 0.0335, 0.0337, 0.0340, 0.0342, 0.0345, 0.0348, 0.0350, 0.0353, 0.0356];
124+
125+
assert(approxEqual(xs.interp1(interpolation), data, 1e-4, 1e-4));
126+
}

0 commit comments

Comments
 (0)