Report a bug
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.linear

## Linear Interpolation

Authors:
Ilya Yaroshenko
Linear!(F, N, X) linear(F, size_t N = 1, X = F)(Repeat!(N, Slice!(RCI!(immutable(X)))) grid, Slice!(RCI!(const(F)), N) values);
Constructs multivariate linear interpolant with nodes on rectilinear grid.
Parameters:
 Repeat!(N, Slice!(RCI!(immutable(X)))) grid x values for interpolant Slice!(RCI!(const(F)), N) values f(x) values for interpolant

Constraints grid, values must have the same length >= 2

Returns:
Examples:
R -> R: Linear interpolation
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));
Examples:
R^2 -> R: Bilinear 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_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));
Examples:
R^3 -> R: Trilinear 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
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));
struct Linear(F, size_t N = 1, X = F) if (N && (N <= 6));
Multivariate linear interpolant with nodes on rectilinear grid.
Slice!(RCI!(const(F)), 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, Slice!(RCI!(const(F)), N) data);
const @property Linear lightConst()();
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()();
enum uint derivativeOrder;
template opCall(uint derivative = 0) if (derivative <= derivativeOrder)
const scope @trusted auto opCall(X...)(in X xs)
if (X.length == N);
(x) operator.

Complexity O(log(grid.length))

alias withDerivative = opCall!1;
struct LinearKernel(X);
const @property auto lightConst()();
immutable @property auto lightImmutable()();
this()(X x0, X x1, X x);
template opCall(uint derivative = 0) if (derivative <= 1)
auto opCall(Y)(in Y y0, in Y y1);
alias withDerivative = opCall!1;