License

.

Copyright

Copyright © 2017, Kaleidic Associates Advisory Limited and Ilya Yaroshenko

Authors

Ilya Yaroshenko

See Also:

Example

import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: vmap;

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).all!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).all!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).all!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).all!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).all!approxEqual([
        12.37135558,   4.99638066,  10.74362441,   0.16008641,
        11.94073593,  16.47908148,  17.49841853,   5.26600921,
        16.21796051,   8.96102894]));

Example

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 = [-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));
}();

Example

R -> R: Cubic interpolation

import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice;

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(all!approxEqual(xs.sliced.vmap(interpolation), data));
}();

Example

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 ////
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));

Example

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 ////
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));