Report a bug
If you spot a problem with this page, click here to create a Github issue.
Improve this page
Quickly fork, edit online, and submit a pull request for this page. Requires a signed-in GitHub account. This works well for small changes. If you'd like to make larger changes you may want to consider using a local clone.

mir.interpolate.spline

Cubic Spline Interpolation

See Also:
Authors:
Ilya Yaroshenko
Examples:
import std.math: approxEqual;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: vmap;
import mir.ndslice.algorithm: all;

auto x = [-1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22].idup.sliced;
auto y = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6].idup.sliced;

auto interpolant = spline!double(x, y); // constructs Spline
auto xs = x + 0.5;                     // input X values for cubic spline

/// not-a-knot (default)
assert(xs.vmap(interpolant).approxEqual([
    -0.68361541,   7.28568719,  10.490694  ,   0.36192032,
    11.91572713,  16.44546433,  17.66699525,   4.52730869,
    19.22825394,  -2.3242592 ]));

/// natural cubic spline
interpolant = spline!double(x, y, SplineBoundaryType.secondDerivative);
assert(xs.vmap(interpolant).approxEqual([
    10.85298372,   5.26255911,  10.71443229,   0.1824536 ,
    11.94324989,  16.45633939,  17.59185094,   4.86340188,
    17.8565408 ,   2.81856494]));

/// set both boundary derivatives to 3
interpolant = spline!double(x, y, SplineBoundaryType.firstDerivative, 3);
assert(xs.vmap(interpolant).approxEqual([
    16.45728263,   4.27981687,  10.82295092,   0.09610695,
    11.95252862,  16.47583126,  17.49964521,   5.26561539,
    16.21803478,   8.96104974]));

/// different boundary conditions
interpolant = spline!double(x, y,
        SplineBoundaryCondition!double(SplineBoundaryType.firstDerivative, 3),
        SplineBoundaryCondition!double(SplineBoundaryType.secondDerivative, -5));
assert(xs.vmap(interpolant).approxEqual([
        16.45730084,   4.27966112,  10.82337171,   0.09403945,
        11.96265209,  16.44067375,  17.6374694 ,   4.67438921,
        18.6234452 ,  -0.05582876]));
// ditto
interpolant = spline!double(x, y,
        SplineBoundaryCondition!double(SplineBoundaryType.secondDerivative, -5),
        SplineBoundaryCondition!double(SplineBoundaryType.firstDerivative, 3));
assert(xs.vmap(interpolant).approxEqual([
        12.37135558,   4.99638066,  10.74362441,   0.16008641,
        11.94073593,  16.47908148,  17.49841853,   5.26600921,
        16.21796051,   8.96102894]));
Examples:
import std.math: approxEqual;
import mir.ndslice.allocation: uninitSlice;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: vmap, map;
import mir.ndslice.algorithm: all;
import mir.functional: aliasCall;

auto x = [-1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22].idup.sliced;
auto y = [
   8.77842512,
   7.96429686,
   7.77074363,
   1.10838032,
   2.69925191,
   1.84922654,
   1.48167283,
   2.8267636 ,
   0.40200172,
   7.78980608].sliced;

auto interpolant = x.spline!double(y); // default boundary condition is 'not-a-knot'

auto xs = x + 0.5;
auto ys = xs.vmap(interpolant);

auto r =
  [ 5.56971848,
    9.30342403,
    4.44139761,
   -0.74740285,
    3.00994108,
    1.50750417,
    1.73144979,
    2.64860361,
    0.64413911,
   10.81768928];

assert(all!approxEqual(ys, r));

// first derivative
auto d1 = xs.vmap(interpolant.aliasCall!"withDerivative").map!"a[1]";
auto r1 =
   [-4.51501279,
     2.15715986,
    -7.28363308,
    -2.14050449,
     0.03693092,
    -0.49618999,
     0.58109933,
    -0.52926703,
     0.7819035 ,
     6.70632693];
assert(all!approxEqual(d1, r1));

// second derivative
auto d2 = xs.vmap(interpolant.aliasCall!"withTwoDerivatives").map!"a[2]";
auto r2 =
   [ 7.07104751,
    -2.62293241,
    -0.01468508,
     5.70609505,
    -2.02358911,
     0.72142061,
     0.25275483,
    -0.6133589 ,
     1.26894416,
     2.68067146];
assert(all!approxEqual(d2, r2));

// third derivative (6 * a)
auto d3 = xs.vmap(interpolant.aliasCall!("opCall", 3)).map!"a[3]";
auto r3 =
   [-3.23132664,
    -3.23132664,
    14.91047457,
    -3.46891432,
     1.88520325,
    -0.16559031,
    -0.44056064,
     0.47057577,
     0.47057577,
     0.47057577];
assert(all!approxEqual(d3, r3));
Examples:
R -> R: Cubic interpolation
import mir.ndslice;
import std.math: approxEqual;

immutable x = [0, 1, 2, 3, 5.00274, 7.00274, 10.0055, 20.0137, 30.0192];
auto y = [0.0011, 0.0011, 0.0030, 0.0064, 0.0144, 0.0207, 0.0261, 0.0329, 0.0356,];
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];

auto interpolation = spline!double(x.sliced, y.sliced);

auto data = 
  [ 0.0011    ,  0.003     ,  0.0064    ,  0.01042622,  0.0144    ,
    0.01786075,  0.0207    ,  0.02293441,  0.02467983,  0.0261    ,
    0.02732764,  0.02840225,  0.0293308 ,  0.03012914,  0.03081002,
    0.03138766,  0.03187161,  0.03227637,  0.03261468,  0.0329    ,
    0.03314357,  0.03335896,  0.03355892,  0.03375674,  0.03396413,
    0.03419436,  0.03446018,  0.03477529,  0.03515072,  0.0356    ];

assert(approxEqual(xs.sliced.vmap(interpolation), data, 1e-4, 1e-4));
Examples:
R^2 -> R: Bicubic interpolation
import std.math: approxEqual;
import mir.ndslice;
alias appreq = (a, b) => approxEqual(a, b, 10e-10, 10e-10);

///// set test function ////
const double y_x0 = 2;
const double y_x1 = -7;
const double y_x0x1 = 3;

// this function should be approximated very well
alias f = (x0, x1) => y_x0 * x0 + y_x1 * x1 + y_x0x1 * x0 * x1 - 11;

///// set interpolant ////
auto x0 = [-1.0, 2, 8, 15].idup.sliced;
auto x1 = [-4.0, 2, 5, 10, 13].idup.sliced;
auto grid = cartesian(x0, x1);

auto interpolant = spline!(double, 2)(x0, x1, grid.map!f);

///// compute test data ////
auto test_grid = cartesian(x0 + 1.23, x1 + 3.23);
// auto test_grid = cartesian(x0 + 0, x1 + 0);
auto real_data = test_grid.map!f;
auto interp_data = test_grid.vmap(interpolant);

///// verify result ////
assert(all!appreq(interp_data, real_data));

//// check derivatives ////
auto z0 = 1.23;
auto z1 = 3.21;
// writeln("-----------------");
// writeln("-----------------");
auto d = interpolant.withDerivative(z0, z1);
assert(appreq(interpolant(z0, z1), f(z0, z1)));
// writeln("d = ", d);
// writeln("interpolant.withTwoDerivatives(z0, z1) = ", interpolant.withTwoDerivatives(z0, z1));
// writeln("-----------------");
// writeln("-----------------");
// writeln("interpolant(z0, z1) = ", interpolant(z0, z1));
// writeln("y_x0 + y_x0x1 * z1 = ", y_x0 + y_x0x1 * z1);
// writeln("y_x1 + y_x0x1 * z0 = ", y_x1 + y_x0x1 * z0);
// writeln("-----------------");
// writeln("-----------------");
// assert(appreq(d[0][0], f(z0, z1)));
// assert(appreq(d[1][0], y_x0 + y_x0x1 * z1));
// assert(appreq(d[0][1], y_x1 + y_x0x1 * z0));
// assert(appreq(d[1][1], y_x0x1));
Examples:
R^3 -> R: Tricubic interpolation
import std.math: approxEqual;
import mir.ndslice;
alias appreq = (a, b) => approxEqual(a, b, 10e-10, 10e-10);

///// set test function ////
const y_x0 = 2;
const y_x1 = -7;
const y_x2 = 5;
const y_x0x1 = 10;
const y_x0x1x2 = 3;

// this function should be approximated very well
alias f = (x0, x1, x2) => y_x0 * x0 + y_x1 * x1 + y_x2 * x2
     + y_x0x1 * x0 * x1 + y_x0x1x2 * x0 * x1 * x2 - 11;

///// set interpolant ////
auto x0 = [-1.0, 2, 8, 15].idup.sliced;
auto x1 = [-4.0, 2, 5, 10, 13].idup.sliced;
auto x2 = [3, 3.7, 5].idup.sliced;
auto grid = cartesian(x0, x1, x2);

auto interpolant = spline!(double, 3)(x0, x1, x2, grid.map!f);

///// compute test data ////
auto test_grid = cartesian(x0 + 1.23, x1 + 3.23, x2 - 3);
auto real_data = test_grid.map!f;
auto interp_data = test_grid.vmap(interpolant);

///// verify result ////
assert(all!appreq(interp_data, real_data));

//// check derivatives ////
auto z0 = 1.23;
auto z1 = 3.23;
auto z2 = -3;
auto d = interpolant.withDerivative(z0, z1, z2);
assert(appreq(interpolant(z0, z1, z2), f(z0, z1, z2)));
assert(appreq(d[0][0][0], f(z0, z1, z2)));

// writeln("-----------------");
// writeln("-----------------");
// auto d = interpolant.withDerivative(z0, z1);
assert(appreq(interpolant(z0, z1, z2), f(z0, z1, z2)));
// writeln("interpolant(z0, z1, z2) = ", interpolant(z0, z1, z2));
// writeln("d = ", d);
// writeln("interpolant.withTwoDerivatives(z0, z1, z2) = ", interpolant.withTwoDerivatives(z0, z1, z2));
// writeln("-----------------");
// writeln("-----------------");
// writeln("interpolant(z0, z1) = ", interpolant(z0, z1));
// writeln("y_x0 + y_x0x1 * z1 = ", y_x0 + y_x0x1 * z1);
// writeln("y_x1 + y_x0x1 * z0 = ", y_x1 + y_x0x1 * z0);
// writeln("-----------------");
// writeln("-----------------");

// writeln("y_x0 + y_x0x1 * z1 + y_x0x1x2 * z1 * z2 = ", y_x0 + y_x0x1 * z1 + y_x0x1x2 * z1 * z2);
// assert(appreq(d[1][0][0], y_x0 + y_x0x1 * z1 + y_x0x1x2 * z1 * z2));
// writeln("y_x1 + y_x0x1 * z0 + y_x0x1x2 * z0 * z2 = ", y_x1 + y_x0x1 * z0 + y_x0x1x2 * z0 * z2);
// assert(appreq(d[0][1][0], y_x1 + y_x0x1 * z0 + y_x0x1x2 * z0 * z2));
// writeln("y_x0x1 + y_x0x1x2 * z2 = ", y_x0x1 + y_x0x1x2 * z2);
// assert(appreq(d[1][1][0], y_x0x1 + y_x0x1x2 * z2));
// writeln("y_x0x1x2 = ", y_x0x1x2);
// assert(appreq(d[1][1][1], y_x0x1x2));
template spline(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIterators = Repeat!(N - 1, FirstGridIterator)) if (isFloatingPoint!T && is(T == Unqual!T) && (N <= 6))
Constructs multivariate cubic spline in symmetrical form with nodes on rectilinear grid. Result has continues second derivatives throughout the curve / nd-surface.
Spline!(T, N, GridIterators) spline(SliceKind ykind, yIterator)(GridVectors grid, scope Slice!(ykind, [N], yIterator) values, SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot, in T valueOfBoundaryConditions = 0);
Parameters:
GridVectors grid immutable x values for interpolant
Slice!(ykind, [N], yIterator) values f(x) values for interpolant
SplineBoundaryType typeOfBoundaries SplineBoundaryType for both tails (optional).
T valueOfBoundaryConditions value of the boundary type (optional).

Constraints grid and values must have the same length >= 3

Returns:
Spline!(T, N, GridIterators) spline(SliceKind ykind, yIterator)(GridVectors grid, scope Slice!(ykind, [N], yIterator) values, SplineBoundaryCondition!T boundaries);
Parameters:
GridVectors grid immutable x values for interpolant
Slice!(ykind, [N], yIterator) values f(x) values for interpolant
SplineBoundaryCondition!T boundaries SplineBoundaryCondition for both tails.

Constraints grid and values must have the same length >= 3

Returns:
Spline!(T, N, GridIterators) spline(SliceKind ykind, yIterator)(GridVectors grid, scope Slice!(ykind, [N], yIterator) values, SplineBoundaryCondition!T rBoundary, SplineBoundaryCondition!T lBoundary);
Parameters:
GridVectors grid immutable x values for interpolant
Slice!(ykind, [N], yIterator) values f(x) values for interpolant
SplineBoundaryCondition!T rBoundary SplineBoundaryCondition for left tail.
SplineBoundaryCondition!T lBoundary SplineBoundaryCondition for right tail.

Constraints grid and values must have the same length >= 3

Returns:
enum SplineBoundaryType: int;
Cubic Spline Boundary Condition Type
periodic
not implemented
notAKnot
(default)
firstDerivative
set the first derivative
secondDerivative
set the second derivative
parabolic
struct SplineBoundaryCondition(T);
Cubic Spline Boundary Condition
SplineBoundaryType type;
type (default is SplineBoundaryType.notAKnot)
T value;
value (default is 0)
struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterators...) if (N && (N <= 6) && (NextGridIterators.length == N - 1));
Multivariate cubic spline with nodes on rectilinear grid.
Slice!(Contiguous, [N], F[2 ^^ N]*) _data;
Aligned buffer allocated with mir.internal.memory. For internal use.
GridIterators _grid;
Grid iterators. For internal use.
nothrow @nogc @trusted this()(GridVectors grid);
@property @trusted void _values(SliceKind kind, Iterator)(scope Slice!(kind, [N], Iterator) values);
Assigns function values to the internal memory. For internal use.
nothrow @nogc @trusted void _computeDerivatives()(SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc);

nothrow @nogc @system void _computeDerivativesTemp()(SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc, ContiguousVector!F temp);
Computes derivatives and stores them in _data. _data is assumed to be preinitialized with function values filled in F[2 ^^ N][0].
Parameters:
SplineBoundaryCondition!F lbc left boundary condition
SplineBoundaryCondition!F rbc right boundary condition
For internal use.
const @property GridVectors[dimension] grid(size_t dimension = 0)()
if (dimension < N);
const @property size_t intervalCount(size_t dimension = 0)();
Returns:
intervals count.
const @property size_t[N] gridShape()();
alias withDerivative = opCall!1;
alias withTwoDerivatives = opCall!2;
enum uint derivativeOrder;
template opCall(uint derivative : 2)
template opCall(uint derivative = 0) if (derivative == 0 || derivative == 1 || derivative == 3)
const @trusted auto opCall(X...)(in X xs)
if (X.length == N);
(x) and [x] operators.

Complexity O(log(points.length))

@trusted void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind skind)(scope Slice!(gkind, [1], IP) points, scope Slice!(vkind, [1], IV) values, scope Slice!(skind, [1], IS) slopes, scope Slice!(Contiguous, [1], T*) temp, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc);
Piecewise cubic hermite interpolating polynomial.
Parameters:
Slice!(gkind, [1], IP) points x values for interpolant
Slice!(vkind, [1], IV) values f(x) values for interpolant
Slice!(skind, [1], IS) slopes uninitialized ndslice to write slopes into
Slice!(Contiguous, [1], T*) temp uninitialized temporary ndslice
SplineBoundaryCondition!F lbc left boundary condition
SplineBoundaryCondition!F rbc right boundary condition

Constraints points, values, and slopes, must have the same length > 3; temp must have length greater or equal to points less minus one.

struct SplineKernel(X);
this()(X x0, X x1, X x);
template opCall(uint derivative = 0) if (derivative <= 3)
const auto opCall(Y)(in Y y0, in Y y1, in Y s0, in Y s1);
alias withDerivative = opCall!1;
alias withTwoDerivatives = opCall!2;