Menu
Search

Library Reference

overview

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

The module provides common C2 splines, monotone (PCHIP) splines, Akima splines and others.
See Also:
SplineType, interp1
License:
Apache-2.0
Authors:
Ilya Yaroshenko
Examples: 
import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;
import mir.ndslice.allocation: rcslice;
import mir.ndslice.topology: vmap;

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

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

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

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

static immutable test_data2 = [
    9.94191781,  5.4223652 , 10.69666392,  0.1971149 , 11.93868415,
    16.46378847, 17.56521661,  4.97656997, 17.39645585, 4.54316446];
/// set both boundary second derivatives to 3
interpolant = spline!double(x, y, SplineBoundaryType.secondDerivative, 3);
assert(xs.vmap(interpolant).all!approxEqual(test_data2));

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

static immutable test_data4 = [
        16.45730084,   4.27966112,  10.82337171,   0.09403945,
        11.96265209,  16.44067375,  17.6374694 ,   4.67438921,
        18.6234452 ,  -0.05582876];
/// different boundary conditions
interpolant = spline!double(x, y,
        SplineBoundaryCondition!double(SplineBoundaryType.firstDerivative, 3),
        SplineBoundaryCondition!double(SplineBoundaryType.secondDerivative, -5));
assert(xs.vmap(interpolant).all!approxEqual(test_data4));


static immutable test_data5 = [
        12.37135558,   4.99638066,  10.74362441,   0.16008641,
        11.94073593,  16.47908148,  17.49841853,   5.26600921,
        16.21796051,   8.96102894];
// ditto
interpolant = spline!double(x, y,
        SplineBoundaryCondition!double(SplineBoundaryType.secondDerivative, -5),
        SplineBoundaryCondition!double(SplineBoundaryType.firstDerivative, 3));
assert(xs.vmap(interpolant).all!approxEqual(test_data5));

static immutable test_data6 = [
    11.40871379,  2.64278898,  9.55774317,  4.84791141, 11.24842121,
    16.16794267, 18.58060557,  5.2531411 , 17.45509005,  1.86992521];
/// Akima spline
interpolant = spline!double(x, y, SplineType.akima);
assert(xs.vmap(interpolant).all!approxEqual(test_data6));

/// Double Quadratic spline
interpolant = spline!double(x, y, SplineType.doubleQuadratic);
import mir.interpolate.utility: ParabolaKernel;
auto kernel1 = ParabolaKernel!double(x[2], x[3], x[4],        y[2], y[3], y[4]);
auto kernel2 = ParabolaKernel!double(      x[3], x[4], x[5],        y[3], y[4], y[5]);
// weighted sum of quadratic functions
auto c = 0.35; // from [0 .. 1]
auto xp = c * x[3] + (1 - c) * x[4];
auto yp = c * kernel1(xp) + (1 - c) * kernel2(xp);
assert(interpolant(xp).approxEqual(yp));
// check parabolic extrapolation of the boundary intervals
kernel1 = ParabolaKernel!double(x[0], x[1], x[2], y[0], y[1], y[2]);
kernel2 = ParabolaKernel!double(x[$ - 3], x[$ - 2], x[$ - 1], y[$ - 3], y[$ - 2], y[$ - 1]);
assert(interpolant(x[0] - 23.421).approxEqual(kernel1(x[0] - 23.421)));
assert(interpolant(x[$ - 1] + 23.421).approxEqual(kernel2(x[$ - 1] + 23.421)));
Examples: 
import mir.rc.array: rcarray;
import mir.algorithm.iteration: all;
import mir.functional: aliasCall;
import mir.math.common: approxEqual;
import mir.ndslice.allocation: uninitSlice;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: vmap, map;

auto x = rcarray!(immutable double)(-1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22).asSlice;
auto y = rcarray(
   8.77842512,
   7.96429686,
   7.77074363,
   1.10838032,
   2.69925191,
   1.84922654,
   1.48167283,
   2.8267636 ,
   0.40200172,
   7.78980608).asSlice;

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.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice;

static immutable x = [0, 1, 2, 3, 5.00274, 7.00274, 10.0055, 20.0137, 30.0192];
static immutable 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.rcslice, 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(all!approxEqual(xs.sliced.vmap(interpolation), data));
Examples:
R^2 -> R: Bicubic interpolation 
import mir.math.common: 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 ////
static immutable x0 = [-1.0, 2, 8, 15];
static immutable x1 = [-4.0, 2, 5, 10, 13];
auto grid = cartesian(x0, x1);

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

///// compute test data ////
auto test_grid = cartesian(x0.sliced + 1.23, x1.sliced + 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 mir.math.common: approxEqual;
import mir.ndslice;
alias appreq = (a, b) => approxEqual(a, b, 10e-10, 10e-10);

///// set test function ////
enum y_x0 = 2;
enum y_x1 = -7;
enum y_x2 = 5;
enum y_x0x1 = 10;
enum 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 ////
static immutable x0 = [-1.0, 2, 8, 15];
static immutable x1 = [-4.0, 2, 5, 10, 13];
static immutable x2 = [3, 3.7, 5];
auto grid = cartesian(x0, x1, x2);

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

///// compute test data ////
auto test_grid = cartesian(x0.sliced + 1.23, x1.sliced + 3.23, x2.sliced - 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));
Examples:
Monotone PCHIP 
import mir.math.common: approxEqual;
import mir.algorithm.iteration: all;
import mir.ndslice.allocation: rcslice;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: vmap;

static immutable x = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22];
static immutable y = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6];
auto interpolant = spline!double(x.rcslice, y.sliced, SplineType.monotone);

auto xs = x[0 .. $ - 1].sliced + 0.5;

auto ys = xs.vmap(interpolant);

assert(ys.all!approxEqual([
    5.333333333333334,
    2.500000000000000,
    10.000000000000000,
    4.288971807628524,
    11.202580845771145,
    16.250000000000000,
    17.962962962962962,
    5.558593750000000,
    17.604662698412699,
    ]));
enum SplineType: int;
Cubic Spline types.
The first derivatives are guaranteed to be continuous for all cubic splines.
c2
Spline with contiguous second derivative.
cardinal
Cardinal and Catmull–Rom splines.
monotone
The interpolant preserves monotonicity in the interpolation data and does not overshoot if the data is not smooth. It is also known as PCHIP in numpy and Matlab.
doubleQuadratic
Weighted sum of two nearbor quadratic functions. It is used in financial analysis.
akima
Akima spline.
template spline(T, size_t N = 1, X = T) 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, X) spline(yIterator, SliceKind ykind)(Repeat!(N, Slice!(RCI!(immutable(X)))) grid, Slice!(yIterator, N, ykind) values, SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot, in T valueOfBoundaryConditions = 0);
Parameters:
Repeat!(N, Slice!(RCI!(immutable(X)))) 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
Spline!(T, N, X) spline(yIterator, SliceKind ykind)(Repeat!(N, Slice!(RCI!(immutable(X)))) grid, Slice!(yIterator, N, ykind) values, SplineBoundaryCondition!T boundaries, SplineType kind = SplineType.c2, in T param = 0);
Parameters:
Repeat!(N, Slice!(RCI!(immutable(X)))) grid immutable x values for interpolant
Slice!(yIterator, N, ykind) values f(x) values for interpolant
SplineBoundaryCondition!T boundaries SplineBoundaryCondition for both tails.
SplineType kind SplineType type of cubic spline.
T param tangent power parameter for cardinal SplineType (ignored by other spline types). Use 1 for zero derivatives at knots and 0 for Catmull–Rom spline.

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

Returns:
Spline
Spline!(T, N, X) spline(yIterator, SliceKind ykind)(Repeat!(N, Slice!(RCI!(immutable(X)))) grid, Slice!(yIterator, N, ykind) values, SplineBoundaryCondition!T rBoundary, SplineBoundaryCondition!T lBoundary, SplineType kind = SplineType.c2, in T param = 0);
Parameters:
Repeat!(N, Slice!(RCI!(immutable(X)))) 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.
SplineType kind SplineType type of cubic spline.
T param tangent power parameter for cardinal SplineType (ignored by other spline types). Use 1 for zero derivatives at knots and 0 for Catmull–Rom spline.

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

Returns:
Spline
enum SplineBoundaryType: int;
Cubic Spline Boundary Condition Type.
See Also:
SplineBoundaryCondition SplineType
periodic
Not implemented.
notAKnot
Not-a-knot (or cubic) boundary condition. It is an aggresive boundary condition that is used only for C2 splines and is default for all API calls. For other then C2 splines, notAKnot is changed internally to a default boundary type for used SplineType.
firstDerivative
Set the first derivative.
secondDerivative
Set the second derivative.
parabolic
Default for Cardinal and Double-Quadratic splines.
monotone
Default for monotone (aka PHCIP ) splines.
akima
Default for Akima splines.
struct SplineBoundaryCondition(T);
Cubic Spline Boundary Condition
See Also:
SplineBoundaryType
SplineBoundaryType type;
type (default is SplineBoundaryType.notAKnot)
T value;
value (default is 0)
struct Spline(F, size_t N = 1, X = F) if (N && (N <= 6));
Multivariate cubic spline with nodes on rectilinear grid.
Slice!(RCI!(F[2 ^^ N]), N) _data;
Aligned buffer allocated with mir.internal.memory. For internal use.
Repeat!(N, RCI!(immutable(X))) _grid;
Grid iterators. For internal use.
@nogc @safe this(Repeat!(N, Slice!(RCI!(immutable(X)))) grid);
@property scope @trusted void _values(SliceKind kind, Iterator)(Slice!(Iterator, N, kind) values);
Assigns function values to the internal memory. For internal use.
nothrow @nogc scope @trusted void _computeDerivatives(SplineType kind, F param, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc);

nothrow @nogc scope @system void _computeDerivativesTemp(SplineType kind, F param, 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
Slice!(F*) temp temporal buffer length points count (optional)
For internal use.
const @property Spline lightConst();
immutable @property Spline lightImmutable();
const @property scope Slice!(RCI!(immutable(X))) grid(size_t dimension = 0)() return
if (dimension < N);
const @property scope @trusted immutable(X)[] gridScopeView(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) operator.

Complexity O(log(points.length))

@trusted void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind skind)(Slice!(IP, 1, gkind) points, Slice!(IV, 1, vkind) values, Slice!(IS, 1, skind) slopes, Slice!(T*) temp, SplineType kind, F param, 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
SplineType kind SplineType type of cubic spline.
F param tangent power parameter for cardinal SplineType (ignored by other spline types). Use 1 for zero derivatives at knots and 0 for Catmull–Rom spline.
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;