Skip to content

Commit 3ea6130

Browse files
9il9il
committed
simplify spline boundaries implementation
1 parent 4885f9d commit 3ea6130

File tree

3 files changed

+104
-112
lines changed

3 files changed

+104
-112
lines changed

‎source/mir/interpolate/constant.d‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -112,7 +112,6 @@ struct Constant(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIter
112112
GridIterators _grid;
113113

114114
import mir.utility: min, max;
115-
package enum alignment = min(64u, F.sizeof).max(size_t.sizeof);
116115

117116
/++
118117
+/

‎source/mir/interpolate/linear.d‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -200,7 +200,6 @@ struct Linear(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
200200
GridIterators _grid;
201201

202202
import mir.utility: min, max;
203-
package enum alignment = min(64u, F.sizeof).max(size_t.sizeof);
204203

205204
/++
206205
+/

‎source/mir/interpolate/spline.d‎

Lines changed: 104 additions & 110 deletions
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,12 @@ import mir.ndslice.traits;
6060
11.94324989, 16.45633939, 17.59185094, 4.86340188,
6161
17.8565408 , 2.81856494]));
6262

63+
/// set both boundary second derivatives to 3
64+
interpolant = spline!double(x, y, SplineBoundaryType.secondDerivative, 3);
65+
assert(xs.vmap(interpolant).all!approxEqual([
66+
9.94191781, 5.4223652 , 10.69666392, 0.1971149 , 11.93868415,
67+
16.46378847, 17.56521661, 4.97656997, 17.39645585, 4.54316446]));
68+
6369
/// set both boundary derivatives to 3
6470
interpolant = spline!double(x, y, SplineBoundaryType.firstDerivative, 3);
6571
assert(xs.vmap(interpolant).all!approxEqual([
@@ -452,7 +458,6 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
452458
GridIterators _grid;
453459

454460
import mir.utility: min, max;
455-
package enum alignment = min(64u, F[2 ^^ N].sizeof).max(size_t.sizeof);
456461

457462
/++
458463
+/
@@ -504,38 +509,31 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
504509
Params:
505510
lbc = left boundary condition
506511
rbc = right boundary condition
512+
temp = temporal buffer length points count (optional)
507513
508514
$(RED For internal use.)
509515
+/
510516
void _computeDerivatives()(SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc) scope @trusted nothrow @nogc
511517
{
512-
import mir.internal.memory;
513518
import mir.algorithm.iteration: maxLength;
514519
auto ml = this._data.maxLength;
515-
auto temp_ptr = cast(F*) alignedAllocate(F[2 ^^ (N - 1)].sizeof * ml, alignment);
516-
if (temp_ptr is null)
517-
assert(0);
518-
debug
519-
{
520-
temp_ptr.sliced(ml)[] = F.init;
521-
}
522-
_computeDerivativesTemp(lbc, rbc, temp_ptr.sliced(ml));
523-
alignedFree(temp_ptr);
520+
auto temp = RCArray!F(ml);
521+
auto tempSlice = temp[].sliced;
522+
_computeDerivativesTemp(lbc, rbc, tempSlice);
524523
}
525524

526525
/// ditto
527526
pragma(inline, false)
528527
void _computeDerivativesTemp()(SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc, Slice!(F*) temp) scope @system nothrow @nogc
529528
{
530-
import mir.internal.memory;
531529
import mir.algorithm.iteration: maxLength, each;
532530
import mir.ndslice.topology: map, byDim, evertPack;
533531

534532
assert(temp.length >= _data.maxLength);
535533

536534
static if (N == 1)
537535
{
538-
splineSlopes!(F, F)(_grid.sliced(_data._lengths[0]), pickDataSubslice(_data, 0), pickDataSubslice(_data, 1), temp, lbc, rbc);
536+
splineSlopes!(F, F)(_grid.sliced(_data._lengths[0]), pickDataSubslice(_data, 0), pickDataSubslice(_data, 1), temp[0 .. _data._lengths[0]], lbc, rbc);
539537
}
540538
else
541539
foreach_reverse(i, ref x; _grid)
@@ -551,7 +549,7 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
551549
auto y = pickDataSubslice(d, l);
552550
auto s = pickDataSubslice(d, L + l);
553551
// debug printf("ptr = %ld, stride = %ld, stride = %ld, d = %ld i = %ld l = %ld\n", d.iterator, d._stride!0, y._stride!0, d.length, i, l);
554-
splineSlopes!(F, F)(x.sliced(_data._lengths[i]), y, s, temp, lbc, rbc);
552+
splineSlopes!(F, F)(x.sliced(_data._lengths[i]), y, s, temp[0 .. _data._lengths[i]], lbc, rbc);
555553
// debug{
556554
// (cast(void delegate() @nogc)(){
557555
// writeln("y = ", y);
@@ -785,170 +783,166 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
785783
assert (points.length >= 2);
786784
assert (points.length == values.length);
787785
assert (points.length == slopes.length);
788-
assert (temp.length + 1 >= points.length);
786+
assert (temp.length == points.length);
789787

790788
auto n = points.length;
791789

792-
auto pd = points.diff;
793-
auto vd = values.diff;
794790

795-
auto xd = cast() pd.front;
796-
auto yd = cast() vd.front;
797-
auto dd = yd / xd;
791+
typeof(slopes[0]) first, last;
798792

799-
// static if (N == 2)
800-
// {
801-
// if (slopes.length!1 != values.length!1)
802-
// assert(0);
803-
// if (values.empty!1)
804-
// return;
805-
// }
793+
auto xd = points.diff;
794+
auto yd = values.diff;
806795

807796
/// special case
808797
static assert(SplineBoundaryType.notAKnot == 0);
809-
with(SplineBoundaryType)
810-
if (_expect(n == 3 && (rbc.type | lbc.type) == 0, false))
798+
if (n <= 3
799+
&& (rbc.type == SplineBoundaryType.parabolic || rbc.type == SplineBoundaryType.notAKnot)
800+
&& (lbc.type == SplineBoundaryType.parabolic || lbc.type == SplineBoundaryType.notAKnot)
801+
)
811802
{
812803
import mir.interpolate.utility;
813-
// static if (N == 1)
804+
if (n == 3)
814805
{
815806
auto parabola = parabolaKernel(points[0], points[1], points[2], values[0], values[1], values[2]);
816807
slopes[0] = parabola.withDerivative(points[0])[1];
817808
slopes[1] = parabola.withDerivative(points[1])[1];
818809
slopes[2] = parabola.withDerivative(points[2])[1];
819810
}
820-
// else
821-
// {
822-
// foreach (i; 0 .. values.length!1)
823-
// {
824-
// auto parabolaDerivative = parabolaKernel!1(points[0], points[1], points[2], values[0][i], values[1][i], values[2][i]);
825-
// slopes[0][i] = parabolaDerivative(points[0]);
826-
// slopes[1][i] = parabolaDerivative(points[1]);
827-
// slopes[2][i] = parabolaDerivative(points[2]);
828-
// }
829-
// }
811+
else
812+
{
813+
assert(slopes.length == 2);
814+
slopes.back = slopes.front = yd.front / xd.front;
815+
assert(0);
816+
}
830817
return;
831818
}
832819

833-
with(SplineBoundaryType) switch(lbc.type)
820+
with(SplineBoundaryType) final switch(lbc.type)
834821
{
835822
case periodic:
836823

837-
if (n > 2)
838-
{
839-
assert(0);
840-
}
841-
goto lsimple;
824+
assert(0);
842825

843826
case notAKnot:
844827

845828
if (n > 2)
846829
{
847-
auto b = pd[1];
848-
auto c = points.diff!2.front;
849-
auto d = ((xd + 2 * c) * b * dd + xd * xd * (vd[1] / b)) / c;
850-
auto r = b;
851-
temp[0] = c / r;
852-
slopes[0].ndassign = d / r;
830+
auto dx0 = xd[0];
831+
auto dx1 = xd[1];
832+
auto dy0 = yd[0];
833+
auto dy1 = yd[1];
834+
auto dd0 = dy0 / dx0;
835+
auto dd1 = dy1 / dx1;
836+
837+
slopes.front = dx1;
838+
first = dx0 + dx1;
839+
temp.front = ((dx0 + 2 * first) * dx1 * dd0 + dx0 ^^ 2 * dd1) / first;
853840
break;
854841
}
842+
else
843+
{
844+
slopes.front = 1;
845+
first = 0;
846+
temp.front = yd.front / xd.front;
855847

856-
lsimple:
857-
858-
temp.front = 0;
859-
slopes.front.ndassign = dd;
848+
}
860849
break;
861850

862851
case firstDerivative:
863852

864-
temp.front = 0;
865-
slopes.front.ndassign = lbc.value;
853+
slopes.front = 1;
854+
first = 0;
855+
temp.front = lbc.value;
866856
break;
867857

868858
case secondDerivative:
869859

870-
temp[0] = 0.5f;
871-
slopes[0].ndassign = 1.5f * dd - 0.25f * lbc.value * xd;
860+
slopes.front = 2;
861+
first = 1;
862+
temp.front = 3 * (yd.front / xd.front) - 0.5 * lbc.value * xd.front;
872863
break;
873864

874865
case parabolic:
875866

876-
if (n > 2 || rbc.type == SplineBoundaryType.periodic)
877-
{
878-
temp[0] = 1;
879-
slopes[0].ndassign = 2 * dd;
880-
break;
881-
}
882-
883-
slopes[0].ndassign = slopes[1].ndassign = dd;
884-
return;
885-
886-
default: assert(0);
887-
}
888-
889-
foreach (i; 1 .. n - 1)
890-
{
891-
auto xq = pd[i];
892-
auto a = xq;
893-
auto b = 2 * (xd + xq);
894-
auto c = xd;
895-
auto r = b - a * temp[i - 1];
896-
temp[i] = c / r;
897-
auto yq = vd[i];
898-
auto dq = yq / xq;
899-
auto d = 3 * (dq * xd + dd * xq);
900-
slopes[i].ndassign = (d - a * slopes[i - 1]) / r;
901-
xd = xq;
902-
yd = yq;
903-
dd = dq;
867+
slopes.front = 2;
868+
first = 1;
869+
temp.front = 2 * (yd.front / xd.front);
870+
break;
904871
}
905872

906-
with(SplineBoundaryType) switch(rbc.type)
873+
with(SplineBoundaryType) final switch(rbc.type)
907874
{
908875
case periodic:
909-
if (n > 2)
910-
{
911-
assert(0);
912-
}
913-
goto rsimple;
876+
assert(0);
914877

915878
case notAKnot:
916879
if (n > 2)
917880
{
918-
auto a = points.diff!2[n - 3];
919-
auto b = pd[n - 3];
920-
auto r = b - a * temp[n - 2];
921-
auto d = ((xd + 2 * a) * b * dd + xd * xd * (vd[n - 3] / b)) / a;
922-
slopes[n - 1].ndassign = (d - a * slopes[n - 2]) / r;
923-
break;
881+
auto dx0 = xd[$ - 1];
882+
auto dx1 = xd[$ - 2];
883+
auto dy0 = yd[$ - 1];
884+
auto dy1 = yd[$ - 2];
885+
auto dd0 = dy0 / dx0;
886+
auto dd1 = dy1 / dx1;
887+
slopes.back = dx1;
888+
last = dx0 + dx1;
889+
temp.back = ((dx0 + 2 * last) * dx1 * dd0 + dx0 ^^ 2 * dd1) / last;
890+
}
891+
else
892+
{
893+
slopes.back = 1;
894+
last = 0;
895+
temp.back = yd.back / xd.back;
924896
}
925-
926-
rsimple:
927-
928-
slopes[n - 1].ndassign = dd;
929897
break;
930-
898+
931899
case firstDerivative:
932900

933-
slopes[n - 1].ndassign = rbc.value;
901+
slopes.back = 1;
902+
last = 0;
903+
temp.back = rbc.value;
934904
break;
935905

936906
case secondDerivative:
937907

938-
slopes[n - 1].ndassign = (3 * dd + 0.5f * rbc.value * xd - slopes[n - 2]) / (2 - temp[n - 2]);
908+
slopes.back = 2;
909+
last = 1;
910+
temp.back = 3 * (yd.back / xd.back) + 0.5 * rbc.value * xd.back;
939911
break;
940912

941913
case parabolic:
942914

943-
slopes[n - 1].ndassign = (2 * dd - slopes[n - 2]) / (1 - temp[n - 2]);
915+
slopes.back = 2;
916+
last = 1;
917+
temp.back = 2 * (yd.back / xd.back);
944918
break;
919+
}
945920

946-
default: assert(0);
921+
foreach (i; 1 .. n - 1)
922+
{
923+
auto dx0 = xd[i - 1];
924+
auto dx1 = xd[i - 0];
925+
auto dy0 = yd[i - 1];
926+
auto dy1 = yd[i - 0];
927+
slopes[i] = 2 * (dx0 + dx1);
928+
temp[i] = 3 * (dy0 / dx0 * dx1 + dy1 / dx1 * dx0);
947929
}
948930

931+
foreach (i; 0 .. n - 1)
932+
{
933+
auto c = i == 0 ? first : xd[i - 1];
934+
auto a = i == n - 2 ? last : xd[i + 1];
935+
auto w = a / slopes[i];
936+
slopes[i + 1] -= w * c;
937+
temp[i + 1] -= w * temp[i];
938+
}
939+
940+
slopes.back = temp.back / slopes.back;
941+
949942
foreach_reverse (i; 0 .. n - 1)
950943
{
951-
slopes[i].ndassign !"-"= temp[i] * slopes[i + 1];
944+
auto c = i == 0 ? first : xd[i - 1];
945+
slopes[i] = (temp[i] - c * slopes[i + 1]) / slopes[i];
952946
}
953947
}
954948

0 commit comments

Comments
 (0)