import mir.interpolate.spline;
import mir.math.common: log;
import mir.ndslice.allocation: rcslice;
import mir.test;

alias g = (double x, uint d = 0) =>
    d == 0 ? 3 * x ^^ 3 + 5 * x ^^ 2 + 0.23 * x + 2 :
    d == 1 ? 9 * x ^^ 2 + 10 * x + 0.23 :
    double.nan;

alias f = (double x, ref scope y)
{
    y[0] = log(x);
    static if (y.length >= 2)
        y[1] = 1 / x;
    static if (y.length >= 3)
        y[2] = -y[1] * y[1];
    static if (y.length >= 4)
        y[3] = -2 * y[1] * y[2];
    static if (y.length >= 5)
        static assert(0, "Not implemented");
};

auto s = spline!double(
    [0.1, 0.4, 0.5, 1.0].rcslice!(immutable double),
    [g(0.1), g(0.4), g(0.5), g(1.0)].rcslice!(const double)
);

auto m = s.interpolationMap!f;

m(0.7).shouldApprox == log(g(0.7));
auto d = m.opCall!3(0.7);
d[0].shouldApprox == log(g(0.7));
d[1].shouldApprox == 1 / g(0.7) * g(0.7, 1);
d[2].shouldApprox == -0.252301;
d[3].shouldApprox == -4.03705;