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.algorithm.iteration: 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.algorithm.iteration: 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;

()@trusted{

    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    ];

()@trusted{
    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(yIterator, SliceKind ykind)(GridVectors grid, scope Slice!(yIterator, N, ykind) values, SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot, in T valueOfBoundaryConditions = 0);
Parameters:
GridVectors grid immutable x values for interpolant
Slice!(yIterator, N, ykind) 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(yIterator, SliceKind ykind)(GridVectors grid, scope Slice!(yIterator, N, ykind) values, SplineBoundaryCondition!T boundaries);
Parameters:
GridVectors grid immutable x values for interpolant
Slice!(yIterator, N, ykind) 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(yIterator, SliceKind ykind)(GridVectors grid, scope Slice!(yIterator, N, ykind) values, SplineBoundaryCondition!T rBoundary, SplineBoundaryCondition!T lBoundary);
Parameters:
GridVectors grid immutable x values for interpolant
Slice!(yIterator, N, ykind) 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.
mir_slice!(mir_rci!(F[2 ^^ N]), N) _data;
Aligned buffer allocated with mir.internal.memory. For internal use.
GridIterators _grid;
Grid iterators. For internal use.
@nogc @safe this(GridVectors grid);
@property scope @trusted void _values(SliceKind kind, Iterator)(scope Slice!(Iterator, N, kind) values);
Assigns function values to the internal memory. For internal use.
nothrow @nogc scope @trusted void _computeDerivatives()(SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc);

nothrow @nogc scope @system void _computeDerivativesTemp()(SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc, Slice!(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 scope GridVectors[dimension] grid(size_t dimension = 0)() return
if (dimension < N);
const @property scope size_t intervalCount(size_t dimension = 0)();
Returns:
intervals count.
const @property scope 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 scope @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!(IP, 1, gkind) points, scope Slice!(IV, 1, vkind) values, scope Slice!(IS, 1, skind) slopes, scope Slice!(T*) temp, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc);
Piecewise cubic hermite interpolating polynomial.
Parameters:
Slice!(IP, 1, gkind) points x values for interpolant
Slice!(IV, 1, vkind) values f(x) values for interpolant
Slice!(IS, 1, skind) slopes uninitialized ndslice to write slopes into
Slice!(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;