linear

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

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,];
static immutable 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 interpolant = linear!double(x.rcslice!(immutable double), y.rcslice!(const double));

static immutable data = [0.0011, 0.0030, 0.0064, 0.0104, 0.0144, 0.0176, 0.0207, 0.0225, 0.0243, 0.0261, 0.0268, 0.0274, 0.0281, 0.0288, 0.0295, 0.0302, 0.0309, 0.0316, 0.0322, 0.0329, 0.0332, 0.0335, 0.0337, 0.0340, 0.0342, 0.0345, 0.0348, 0.0350, 0.0353, 0.0356];

assert(xs.sliced.vmap(interpolant).all!((a, b) => approxEqual(a, b, 1e-4, 1e-4))(data));

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_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)
    .map!f
    .rcslice
    .lightConst;

auto interpolant = 
    linear!(double, 2)(
        x0.rcslice!(immutable double),
        x1.rcslice!(immutable double),
        grid
    );

///// compute test data ////
auto test_grid = cartesian(x0.sliced + 1.23, x1.sliced + 3.23);
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;
auto d = interpolant.withDerivative(z0, z1);
assert(appreq(interpolant(z0, z1), f(z0, z1)));
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));

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
static auto f(double x0, double x1, double x2)
{
    return 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)
    .map!f
    .as!(const double)
    .rcslice;

auto interpolant = linear!(double, 3)(
        x0.rcslice!(immutable double),
        x1.rcslice!(immutable double),
        x2.rcslice!(immutable double),
        grid);

///// 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.21;
auto z2 = 4;
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)));
assert(appreq(d[1][0][0], y_x0 + y_x0x1 * z1 + y_x0x1x2 * z1 * z2));
assert(appreq(d[0][1][0], y_x1 + y_x0x1 * z0 + y_x0x1x2 * z0 * z2));
assert(appreq(d[1][1][0], y_x0x1 + y_x0x1x2 * z2));
assert(appreq(d[1][1][1], y_x0x1x2));