Skip to content

Commit feb42e9

Browse files
9il9il
committed
continue rework interpolation
1 parent 046b7b3 commit feb42e9

File tree

4 files changed

+101
-74
lines changed

4 files changed

+101
-74
lines changed

‎source/mir/functional.d‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -718,10 +718,16 @@ template aliasCall(string methodName, TemplateArgs...)
718718
Returns:
719719
wrapped value with implemented opCall and opIndex methods
720720
+/
721-
AliasCall!(T, methodName, TemplateArgs) aliasCall(T)(auto ref T value)
721+
AliasCall!(T, methodName, TemplateArgs) aliasCall(T)(T value) @property
722722
{
723723
return typeof(return)(value);
724724
}
725+
726+
/// ditto
727+
ref AliasCall!(T, methodName, TemplateArgs) aliasCall(T)(return ref T value) @property @trusted
728+
{
729+
return *cast(typeof(return)*) &value;
730+
}
725731
}
726732

727733
///

‎source/mir/interpolate/linear.d‎

Lines changed: 29 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -310,8 +310,9 @@ struct Linear(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators...)
310310
// @FUTURE@
311311
// X.length == N || derivative == 0 && X.length && X.length <= N
312312
{
313+
import mir.functional: AliasCall;
313314
import mir.ndslice.topology: iota;
314-
alias Kernel = LinearKernel!(derivative, F);
315+
alias Kernel = AliasCall!(LinearKernel!F, "opCall", derivative);
315316

316317
size_t[N] indexes = void;
317318
Kernel[N] kernels = void;
@@ -330,7 +331,7 @@ struct Linear(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators...)
330331
alias x = xs[i];
331332
indexes[i] = this.findInterval!i(x);
332333
}
333-
kernels[i] = Kernel(_grid[i][indexes[i]], _grid[i][indexes[i] + 1], x);
334+
kernels[i] = LinearKernel!F(_grid[i][indexes[i]], _grid[i][indexes[i] + 1], x);
334335
}
335336

336337
align(64) F[2 ^^ N][derivative + 1] local = void;
@@ -373,13 +374,14 @@ struct Linear(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators...)
373374
alias withDerivative = opCall!1;
374375
}
375376

376-
struct LinearKernel(uint derivative, X)
377-
if (derivative <= 3)
377+
///
378+
struct LinearKernel(X)
378379
{
379380
X step = 0;
380381
X w0 = 0;
381382
X w1 = 0;
382383

384+
///
383385
this()(X x0, X x1, X x)
384386
{
385387
step = x1 - x0;
@@ -389,22 +391,31 @@ struct LinearKernel(uint derivative, X)
389391
w1 = c1 / step;
390392
}
391393

392-
auto opCall(Y)(in Y y0, in Y y1)
394+
///
395+
template opCall(uint derivative = 0)
396+
if (derivative <= 1)
393397
{
394-
auto r0 = y0 * w1;
395-
auto r1 = y1 * w0;
396-
auto y = r0 + r1;
397-
static if (derivative)
398+
///
399+
auto opCall(Y)(in Y y0, in Y y1)
398400
{
399-
auto diff = y1 - y0;
400-
Y[derivative + 1] ret = 0;
401-
ret[0] = y;
402-
ret[1] = diff / step;
403-
return ret;
404-
}
405-
else
406-
{
407-
return y;
401+
auto r0 = y0 * w1;
402+
auto r1 = y1 * w0;
403+
auto y = r0 + r1;
404+
static if (derivative)
405+
{
406+
auto diff = y1 - y0;
407+
Y[derivative + 1] ret = 0;
408+
ret[0] = y;
409+
ret[1] = diff / step;
410+
return ret;
411+
}
412+
else
413+
{
414+
return y;
415+
}
408416
}
409417
}
418+
419+
///
420+
alias withDerivative = opCall!1;
410421
}

‎source/mir/interpolate/spline.d‎

Lines changed: 47 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -595,7 +595,7 @@ struct Spline(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators...)
595595
// X.length == N || derivative == 0 && X.length && X.length <= N
596596
{
597597
import mir.ndslice.topology: iota;
598-
alias Kernel = SplineKernel!(derivative, F);
598+
alias Kernel = AliasCall!(SplineKernel!F, "opCall", derivative);
599599
enum rp2d = derivative == 3 ? 2 : derivative;
600600

601601
size_t[N] indexes = void;
@@ -613,7 +613,7 @@ struct Spline(F, size_t N = 1, FirstGridIterator = F*, NextGridIterators...)
613613
alias x = xs[i];
614614
indexes[i] = this.findInterval!i(x);
615615
}
616-
kernels[i] = Kernel(_grid[i][indexes[i]], _grid[i][indexes[i] + 1], x);
616+
kernels[i] = SplineKernel!F(_grid[i][indexes[i]], _grid[i][indexes[i] + 1], x);
617617
}
618618

619619
align(64) F[2 ^^ N * 2 ^^ N][2] local = void;
@@ -728,10 +728,10 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
728728
import mir.interpolate.utility;
729729
// static if (packs == [1])
730730
{
731-
auto parabolaDerivative = parabolaKernel!1(points[0], points[1], points[2], values[0], values[1], values[2]);
732-
slopes[0] = parabolaDerivative(points[0]);
733-
slopes[1] = parabolaDerivative(points[1]);
734-
slopes[2] = parabolaDerivative(points[2]);
731+
auto parabola = parabolaKernel(points[0], points[1], points[2], values[0], values[1], values[2]);
732+
slopes[0] = parabola.withDerivative(points[0])[1];
733+
slopes[1] = parabola.withDerivative(points[1])[1];
734+
slopes[2] = parabola.withDerivative(points[2])[1];
735735
}
736736
// else
737737
// {
@@ -868,14 +868,15 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
868868
}
869869
}
870870

871-
struct SplineKernel(uint derivative, X)
872-
if (derivative <= 3)
871+
///
872+
struct SplineKernel(X)
873873
{
874874
X step = 0;
875875
X w0 = 0;
876876
X w1 = 0;
877877
X wq = 0;
878878

879+
///
879880
this()(X x0, X x1, X x)
880881
{
881882
step = x1 - x0;
@@ -886,37 +887,48 @@ struct SplineKernel(uint derivative, X)
886887
wq = w0 * w1;
887888
}
888889

889-
auto opCall(Y)(in Y y0, in Y y1, in Y s0, in Y s1) const
890+
///
891+
template opCall(uint derivative = 0)
892+
if (derivative <= 3)
890893
{
891-
auto diff = y1 - y0;
892-
auto z0 = s0 * step - diff;
893-
auto z1 = s1 * step - diff;
894-
auto a0 = z0 * w1;
895-
auto a1 = z1 * w0;
896-
auto pr = a0 - a1;
897-
auto b0 = y0 * w1;
898-
auto b1 = y1 * w0;
899-
auto pl = b0 + b1;
900-
auto y = pl + wq * pr;
901-
static if (derivative)
894+
///
895+
auto opCall(Y)(in Y y0, in Y y1, in Y s0, in Y s1) const
902896
{
903-
Y[derivative + 1] ret = 0;
904-
ret[0] = y;
905-
auto wd = w1 - w0;
906-
auto zd = z1 + z0;
907-
ret[1] = (diff + (wd * pr - wq * zd)) / step;
908-
static if (derivative > 1)
897+
auto diff = y1 - y0;
898+
auto z0 = s0 * step - diff;
899+
auto z1 = s1 * step - diff;
900+
auto a0 = z0 * w1;
901+
auto a1 = z1 * w0;
902+
auto pr = a0 - a1;
903+
auto b0 = y0 * w1;
904+
auto b1 = y1 * w0;
905+
auto pl = b0 + b1;
906+
auto y = pl + wq * pr;
907+
static if (derivative)
909908
{
910-
auto astep = zd / (step * step);
911-
ret[2] = -3 * wd * astep + (s1 - s0) / step;
912-
static if (derivative > 2)
913-
ret[3] = 6 * astep / step;
909+
Y[derivative + 1] ret = 0;
910+
ret[0] = y;
911+
auto wd = w1 - w0;
912+
auto zd = z1 + z0;
913+
ret[1] = (diff + (wd * pr - wq * zd)) / step;
914+
static if (derivative > 1)
915+
{
916+
auto astep = zd / (step * step);
917+
ret[2] = -3 * wd * astep + (s1 - s0) / step;
918+
static if (derivative > 2)
919+
ret[3] = 6 * astep / step;
920+
}
921+
return ret;
922+
}
923+
else
924+
{
925+
return y;
914926
}
915-
return ret;
916-
}
917-
else
918-
{
919-
return y;
920927
}
921928
}
929+
930+
///
931+
alias withDerivative = opCall!1;
932+
///
933+
alias withTwoDerivatives = opCall!2;
922934
}

‎source/mir/interpolate/utility.d‎

Lines changed: 18 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -1,20 +1,18 @@
11
module mir.interpolate.utility;
22

33
import mir.ndslice.slice;
4+
import std.traits;
45

56
/++
67
ParabolaKernel structure.
78
+/
8-
struct ParabolaKernel(uint derivative, T)
9-
if (derivative <= 2)
9+
struct ParabolaKernel(T)
1010
{
1111
///
1212
T a = 0;
1313
///
14-
static if (derivative <= 1)
1514
T b = 0;
1615
///
17-
static if (derivative == 0)
1816
T c = 0;
1917

2018
/// Builds parabola given three points
@@ -30,35 +28,35 @@ struct ParabolaKernel(uint derivative, T)
3028
auto a3 = bm * a1 + a2;
3129
auto d3 = bm * d1 + d2;
3230
a = d3 / a3;
33-
static if (derivative <= 1)
3431
b = (d1 - a1 * a) / b1;
35-
static if (derivative == 0)
3632
c = y1 - x1 * (a * x1 + b);
3733
}
3834

3935
///
40-
auto opCall(V)(V x)
41-
if (!isSlice!V)
36+
auto opCall(uint derivative = 0)(T x) const
37+
if (derivative <= 2)
4238
{
39+
auto y = (a * x + b) * x + c;
4340
static if (derivative == 0)
44-
return (a * x + b) * x + c;
45-
else
46-
static if (derivative == 1)
47-
return 2 * a * x + b;
41+
return y;
4842
else
49-
return 2 * a;
43+
{
44+
auto y1 = 2 * a * x + b;
45+
static if (derivative == 1)
46+
return cast(T[2])[y, y1];
47+
else
48+
return cast(T[3])[y, y1, 2 * a];
49+
}
5050
}
5151

52-
/// ditto
53-
auto opCall(SliceKind kind, size_t[] packs, Iterator)(Slice!(kind, packs, iterator) x)
54-
{
55-
import mir.ndslice.topology: vmap;
56-
return x.vmap(this);
57-
}
52+
///
53+
alias withDerivative = opCall!1;
54+
///
55+
alias withTwoDerivatives = opCall!2;
5856
}
5957

6058
/// ditto
61-
ParabolaKernel!(derivative, typeof(X.init - Y.init)) parabolaKernel(uint derivative = 0, X, Y)(in X x0, in X x1, in X x2, in Y y0, in Y y1, in Y y2)
59+
ParabolaKernel!(Unqual!(typeof(X.init - Y.init))) parabolaKernel(X, Y)(in X x0, in X x1, in X x2, in Y y0, in Y y1, in Y y2)
6260
{
6361
return typeof(return)(x0, x1, x2, y0, y1, y2);
6462
}

0 commit comments

Comments
 (0)