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));
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 typeOfBondaries = SplineBoundaryType.notAKnot, in T valueOfBondaryConditions = 0);
Parameters:
GridVectors grid immutable x values for interpolant
Slice!(ykind, [N], yIterator) values f(x) values for interpolant
SplineBoundaryType typeOfBondaries SplineBoundaryType for both tails (optional).
T valueOfBondaryConditions 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 bondaries);
Parameters:
GridVectors grid immutable x values for interpolant
Slice!(ykind, [N], yIterator) values f(x) values for interpolant
SplineBoundaryCondition!T bondaries 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 _computeDerivatives()(SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc, ContiguousVector!F temp);
Computes derivatives and stores them in _data. _data is assumeed 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;