Skip to content

Commit a7004c3

Browse files
9il9il
committed
rework findInterval
1 parent 8965698 commit a7004c3

File tree

4 files changed

+138
-111
lines changed

4 files changed

+138
-111
lines changed

‎source/mir/interpolate/constant.d‎

Lines changed: 24 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -116,6 +116,13 @@ struct Constant(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators = Rep
116116
import mir.utility: min, max;
117117
package enum alignment = min(64u, F.sizeof).max(size_t.sizeof);
118118

119+
package ref shared(sizediff_t) counter() @trusted @property
120+
{
121+
assert(_ownsData);
122+
auto p = cast(shared sizediff_t*) _data.ptr;
123+
return *(p - 1);
124+
}
125+
119126
///
120127
this(this) @safe nothrow @nogc
121128
{
@@ -124,13 +131,6 @@ struct Constant(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators = Rep
124131
counter.atomicOp!"+="(1);
125132
}
126133

127-
///
128-
GridVectors[dimension] grid(size_t dimension = 0)() @property
129-
if (dimension < N)
130-
{
131-
return _grid[dimension].sliced(_data._lengths[dimension]);
132-
}
133-
134134
/++
135135
Frees _data if $(LREF Spline._freeData) is true.
136136
+/
@@ -144,13 +144,6 @@ struct Constant(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators = Rep
144144
alignedFree(cast(void*)(_data.ptr) - alignment);
145145
}
146146

147-
package ref shared(sizediff_t) counter() @trusted @property
148-
{
149-
assert(_ownsData);
150-
auto p = cast(shared sizediff_t*) _data.ptr;
151-
return *(p - 1);
152-
}
153-
154147
/++
155148
+/
156149
this()(GridVectors grid) @trusted nothrow @nogc
@@ -191,53 +184,53 @@ struct Constant(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators = Rep
191184
}
192185

193186
///
194-
size_t[N] dataShape()() @property
187+
GridVectors[dimension] grid(size_t dimension = 0)() const @property
188+
if (dimension < N)
195189
{
196-
return _data.shape;
190+
return _grid[dimension].sliced(_data._lengths[dimension]);
197191
}
198192

199193
/++
200-
Interval index for x.
194+
Returns: intervals count.
201195
+/
202-
size_t interval(size_t dimension = 0, X)(in X x)
196+
size_t intervalCount(size_t dimension = 0)() const @property
203197
{
204-
import std.range: assumeSorted;
205-
return _data._lengths[dimension] - 1 - _grid[dimension].sliced(_data._lengths[dimension])[1 .. $]
206-
.assumeSorted
207-
.upperBound(x)
208-
.length;
198+
assert(_data._lengths[dimension] > 0);
199+
return _data._lengths[dimension] - 0;
209200
}
210201

211202
///
212-
enum uint maxDerivative = 0;
203+
size_t[N] gridShape()() const @property
204+
{
205+
return _data.shape;
206+
}
207+
208+
///
209+
enum uint derivativeOrder = 0;
213210

214211
///
215212
template opCall(uint derivative = 0)
216-
if (derivative <= maxDerivative)
213+
if (derivative <= derivativeOrder)
217214
{
218215
@trusted:
219216
/++
220217
`(x)` and `[x]` operators.
221218
Complexity:
222219
`O(log(grid.length))`
223220
+/
224-
auto opCall(X...)(in X xs)
221+
auto opCall(X...)(in X xs) const
225222
if (X.length == N)
226223
// @FUTURE@
227224
// X.length == N || derivative == 0 && X.length && X.length <= N
228225
{
229-
import mir.ndslice.topology: iota;
230-
231226
size_t[N] indexes = void;
232227

233-
enum rp2d = derivative;
234-
235228
foreach(i; Iota!N)
236229
{
237230
static if (isInterval!(typeof(xs[i])))
238231
indexes[i] = xs[i][1];
239232
else
240-
indexes[i] = interval!i(xs[i]);
233+
indexes[i] = this.findInterval!i(xs[i]);
241234
}
242235
return _data[indexes];
243236
}

‎source/mir/interpolate/linear.d‎

Lines changed: 24 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -90,14 +90,6 @@ version(mir_test)
9090
assert(approxEqual(xs.sliced.map!interpolation, data, 1e-4, 1e-4));
9191
}
9292

93-
unittest
94-
{
95-
auto x = [0.0, 1, 2].sliced;
96-
auto y = [10.0, 2, 4].sliced;
97-
auto interpolation = linear!double(x, y);
98-
assert(interpolation.interval(1.0) == 1);
99-
}
100-
10193
/// R^2 -> R: Bilinear interpolaiton
10294
unittest
10395
{
@@ -207,6 +199,13 @@ struct Linear(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators...)
207199
import mir.utility: min, max;
208200
package enum alignment = min(64u, F.sizeof).max(size_t.sizeof);
209201

202+
package ref shared(sizediff_t) counter() @trusted @property
203+
{
204+
assert(_ownsData);
205+
auto p = cast(shared sizediff_t*) _data.ptr;
206+
return *(p - 1);
207+
}
208+
210209
///
211210
this(this) @safe nothrow @nogc
212211
{
@@ -215,13 +214,6 @@ struct Linear(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators...)
215214
counter.atomicOp!"+="(1);
216215
}
217216

218-
///
219-
GridVectors[dimension] grid(size_t dimension = 0)() @property
220-
if (dimension < N)
221-
{
222-
return _grid[dimension].sliced(_data._lengths[dimension]);
223-
}
224-
225217
/++
226218
Frees _data if $(LREF Spline._freeData) is true.
227219
+/
@@ -235,13 +227,6 @@ struct Linear(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators...)
235227
alignedFree(cast(void*)(_data.ptr) - alignment);
236228
}
237229

238-
package ref shared(sizediff_t) counter() @trusted @property
239-
{
240-
assert(_ownsData);
241-
auto p = cast(shared sizediff_t*) _data.ptr;
242-
return *(p - 1);
243-
}
244-
245230
/++
246231
+/
247232
this()(GridVectors grid) @trusted nothrow @nogc
@@ -285,37 +270,40 @@ struct Linear(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators...)
285270
}
286271

287272
///
288-
size_t[N] dataShape()() @property
273+
GridVectors[dimension] grid(size_t dimension = 0)() const @property
274+
if (dimension < N)
289275
{
290-
return _data.shape;
276+
return _grid[dimension].sliced(_data._lengths[dimension]);
291277
}
292278

293279
/++
294-
Interval index for x.
280+
Returns: intervals count.
295281
+/
296-
size_t interval(size_t dimension = 0, X)(in X x)
282+
size_t intervalCount(size_t dimension = 0)() const @property
283+
{
284+
assert(_data._lengths[dimension] > 1);
285+
return _data._lengths[dimension] - 1;
286+
}
287+
288+
///
289+
size_t[N] gridShape()() const @property
297290
{
298-
import std.range: assumeSorted;
299-
return _data._lengths[dimension] - 2 - _grid[dimension].sliced(_data._lengths[dimension])[1 .. $ - 1]
300-
.assumeSorted
301-
.upperBound(x)
302-
.length;
291+
return _data.shape;
303292
}
304293

305294
///
306-
enum uint maxDerivative = 1;
295+
enum uint derivativeOrder = 1;
307296

308297
///
309298
template opCall(uint derivative = 0)
310-
if (derivative <= maxDerivative)
299+
if (derivative <= derivativeOrder)
311300
{
312-
@trusted:
313301
/++
314302
`(x)` and `[x]` operators.
315303
Complexity:
316304
`O(log(grid.length))`
317305
+/
318-
auto opCall(X...)(in X xs)
306+
auto opCall(X...)(in X xs) const @trusted
319307
if (X.length == N)
320308
// @FUTURE@
321309
// X.length == N || derivative == 0 && X.length && X.length <= N
@@ -338,7 +326,7 @@ struct Linear(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators...)
338326
else
339327
{
340328
alias x = xs[i];
341-
indexes[i] = interval!i(x);
329+
indexes[i] = this.findInterval!i(x);
342330
}
343331
kernels[i] = Kernel(_grid[i][indexes[i]], _grid[i][indexes[i] + 1], x);
344332
}

‎source/mir/interpolate/package.d‎

Lines changed: 58 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,51 @@ package template Repeat(ulong n, L...)
4141
alias Repeat = L[0 .. 0];
4242
}
4343

44+
/++
45+
Interval index for x value given.
46+
+/
47+
template findInterval(size_t dimension = 0)
48+
{
49+
/++
50+
Interval index for x value given.
51+
Params:
52+
interpolant = interpolant
53+
x = X value
54+
+/
55+
size_t findInterval(Interpolant, X)(auto ref const Interpolant interpolant, in X x) @trusted
56+
{
57+
import mir.ndslice.slice: sliced;
58+
import std.range: assumeSorted;
59+
static if (dimension)
60+
{
61+
immutable sizediff_t len = interpolant.intervalCount!dimension - 1;
62+
auto grid = interpolant.grid!dimension[1 .. $][0 .. len];
63+
}
64+
else
65+
{
66+
immutable sizediff_t len = interpolant.intervalCount - 1;
67+
auto grid = interpolant.grid[1 .. $][0 .. len];
68+
}
69+
assert(len >= 0);
70+
return len - grid
71+
.assumeSorted
72+
.upperBound(x)
73+
.length;
74+
}
75+
}
76+
77+
///
78+
unittest
79+
{
80+
import mir.ndslice.slice: sliced;
81+
import mir.interpolate.linear;
82+
83+
auto x = [0.0, 1, 2].sliced;
84+
auto y = [10.0, 2, 4].sliced;
85+
auto interpolation = linear!double(x, y);
86+
assert(interpolation.findInterval(1.0) == 1);
87+
}
88+
4489
/++
4590
Lazy interpolation shell with linear complexity.
4691
@@ -61,14 +106,14 @@ auto interp1(Range, Interpolant)(Range range, Interpolant interpolant, size_t in
61106
}
62107

63108
/// ditto
64-
struct Interp1(Range, Interpolation)
109+
struct Interp1(Range, Interpolant)
65110
{
66-
/// Sorted range (descending)
67-
Range _range;
68-
/// Interpolation structure
69-
Interpolation _interpolation;
70-
/// Current interpolation interval
71-
size_t _interval;
111+
/// Sorted range (descending). $(RED For internal use.)
112+
private Range _range;
113+
/// Interpolant structure. $(RED For internal use.)
114+
private Interpolant _interpolant;
115+
/// Current interpolation interval. $(RED For internal use.)
116+
private size_t _interval;
72117

73118
static if (hasLength!Range)
74119
/// Length (optional)
@@ -91,10 +136,12 @@ struct Interp1(Range, Interpolation)
91136
assert(!empty);
92137
auto x = _range.front;
93138
return (x) @trusted {
94-
auto points = _interpolation.grid;
95-
while (x > points[_interval + 1] && points.length > _interval + 2)
139+
auto points = _interpolant.grid;
140+
sizediff_t len = _interpolant.intervalCount - 1;
141+
assert(len >= 0);
142+
while (x > points[_interval + 1] && _interval < len)
96143
_interval++;
97-
return _interpolation(x.atInterval(_interval));
144+
return _interpolant(x.atInterval(_interval));
98145
} (x);
99146
}
100147
}
@@ -185,7 +232,7 @@ version(LDC)
185232
else
186233
enum LDC = false;
187234

188-
auto copyvec(F, size_t N)(ref F[N] from, ref F[N] to)
235+
auto copyvec(F, size_t N)(ref const F[N] from, ref F[N] to)
189236
{
190237
import mir.internal.utility;
191238

0 commit comments

Comments
 (0)