interpn 0.11.2

N-dimensional interpolation/extrapolation methods, no-std and no-alloc compatible.
Documentation
import numpy as np
import pytest
import interpn


def linear_field(x, y):
    return 1.25 + 2.0 * x - 3.0 * y


def assert_linear_grad(grad, dtype, tol):
    expected = np.array([[2.0], [-3.0]], dtype=dtype)
    np.testing.assert_allclose(grad, expected, atol=tol, rtol=tol)


@pytest.mark.parametrize("dtype,tol", [(np.float64, 1e-12), (np.float32, 1e-5)])
@pytest.mark.parametrize(
    "cls,grid_kind",
    [
        (interpn.MultilinearRegular, "regular"),
        (interpn.MultilinearRectilinear, "rectilinear"),
        (interpn.MulticubicRegular, "regular"),
        (interpn.MulticubicRectilinear, "rectilinear"),
        (interpn.MultiBsplineRegular, "regular"),
        (interpn.MultiBsplineRectilinear, "rectilinear"),
    ],
)
def test_eval_grad_linear_field(cls, grid_kind, dtype, tol):
    if grid_kind == "regular":
        x = np.linspace(-1.0, 2.0, 6).astype(dtype)
        y = np.linspace(0.5, 3.0, 7).astype(dtype)
    else:
        x = np.array([-1.0, -0.2, 0.4, 1.1, 2.0, 3.1], dtype=dtype)
        y = np.array([0.5, 0.7, 1.4, 2.2, 3.0, 4.0, 5.2], dtype=dtype)

    xgrid, ygrid = np.meshgrid(x, y, indexing="ij")
    vals = linear_field(xgrid, ygrid).astype(dtype)
    obs = [np.array([0.23], dtype=dtype), np.array([1.07], dtype=dtype)]

    if grid_kind == "regular":
        dims = [x.size, y.size]
        starts = np.array([x[0], y[0]], dtype=dtype)
        steps = np.array([x[1] - x[0], y[1] - y[0]], dtype=dtype)
        interpolator = cls.new(dims, starts, steps, vals)
    else:
        interpolator = cls.new([x, y], vals)

    grad = interpolator.eval_grad(obs)
    assert grad.shape == (2, 1)
    assert grad.dtype == dtype
    assert_linear_grad(grad, dtype, tol)

    out = np.zeros((2, obs[0].size), dtype=dtype)
    returned = interpolator.eval_grad_unchecked(obs, out)
    assert returned is out
    assert_linear_grad(out, dtype, tol)