#include "test_common.h"

#include "./utils/matrix.h"
#include "./utils/vector.h"


#include "./numerics/interpolation1d.h"


// ------------ helpers ------------
template <typename T>
static void make_uniform_xy(utils::Vector<T>& x, utils::Vector<T>& y,
                            std::uint64_t N, T x0, T dx,
                            T a, T b) {
    x.resize(N, T(0));
    y.resize(N, T(0));
    for (std::uint64_t i=0; i<N; ++i) {
        T xi = x0 + dx * T(i);
        x[i] = xi;
        y[i] = a*xi + b;     // y = a x + b
    }
}
template <typename T>
static void make_uniform_xy_fun(utils::Vector<T>& x, utils::Vector<T>& y,
                                std::uint64_t N, T x0, T dx,
                                T (*f)(T)) {
    x.resize(N, T(0));
    y.resize(N, T(0));
    for (std::uint64_t i=0; i<N; ++i) {
        T xi = x0 + dx * T(i);
        x[i] = xi;
        y[i] = f(xi);
    }
}
template <typename T>
static inline bool almost_eq(T a, T b, double tol=1e-12) {
    return std::fabs(double(a)-double(b)) <= tol;
}

// ------------ tests ------------

// 1) All interpolants reproduce the data points exactly (or to tiny FP error)
TEST_CASE(Interp_Reproduces_Data_Nodes) {
    const std::uint64_t N = 25;
    utils::Vd x, y;
    // linear data: y = 0.1 x + 0.9 on x = 1..25
    make_uniform_xy<double>(x, y, N, 1.0, 1.0, 0.1, 0.9);

    numerics::interp_linear<double>       lin(x,y);
    numerics::interp_polynomial<double>   pol(x,y, 3);
    numerics::interp_cubic_spline<double> spl(x,y);
    numerics::interp_rational<double>     rat(x,y, 2);
    numerics::interp_barycentric<double>  bar(x,y, 2);

    for (std::uint64_t i=0; i<N; ++i) {
        double xi = x[i];
        double yi = y[i];
        CHECK(almost_eq(lin.interp(xi), yi), "linear fails at node");
        CHECK(almost_eq(pol.interp(xi), yi), "polynomial fails at node");
        CHECK(almost_eq(spl.interp(xi), yi), "spline fails at node");
        CHECK(almost_eq(rat.interp(xi), yi), "rational fails at node");
        CHECK(almost_eq(bar.interp(xi), yi), "barycentric fails at node");
    }
}

// 2) Linear dataset: all interpolants match the analytic line at off-grid points
TEST_CASE(Interp_Linear_Dataset_OffGrid) {
    const std::uint64_t N = 100;
    utils::Vd x, y;
    // your quick test dataset: x = 1..100, y = 0.1 x + 0.9
    make_uniform_xy<double>(x, y, N, 1.0, 1.0, 0.1, 0.9);

    numerics::interp_linear<double>       lin(x,y);
    numerics::interp_polynomial<double>   pol(x,y, 3);
    numerics::interp_cubic_spline<double> spl(x,y);
    numerics::interp_rational<double>     rat(x,y, 3);
    numerics::interp_barycentric<double>  bar(x,y, 3);

    std::vector<double> qs = { 1.5, 5.5, 5.51, 50.01, 99.9 };
    for (double q : qs) {
        const double ytrue = 0.1*q + 0.9;
        CHECK(almost_eq(lin.interp(q), ytrue, 1e-9), "linear off-grid mismatch");
        CHECK(almost_eq(pol.interp(q), ytrue, 1e-9), "polynomial off-grid mismatch");
        CHECK(almost_eq(spl.interp(q), ytrue, 1e-9), "spline off-grid mismatch");
        CHECK(almost_eq(rat.interp(q), ytrue, 1e-9), "rational off-grid mismatch");
        CHECK(almost_eq(bar.interp(q), ytrue, 1e-9), "barycentric off-grid mismatch");
    }

    // endpoints should match exactly (no extrapolation)
    CHECK(almost_eq(lin.interp(x[0]),   y[0]), "linear endpoint");
    CHECK(almost_eq(lin.interp(x[N-1]), y[N-1]), "linear endpoint");
    CHECK(almost_eq(spl.interp(x[0]),   y[0]), "spline endpoint");
    CHECK(almost_eq(spl.interp(x[N-1]), y[N-1]), "spline endpoint");
}

// 3) Quadratic dataset: degree-2 polynomial interpolation should be exact
static double f_quad(double t) { return t*t - 3.0*t + 2.0; } // (t-1)(t-2)

TEST_CASE(Interp_Polynomial_Deg2_Exact_On_Quadratic) {
    const std::uint64_t N = 21;
    utils::Vd x, y;
    make_uniform_xy_fun<double>(x, y, N, -2.0, 0.5, &f_quad);

    numerics::interp_polynomial<double> pol(x,y, 3);

    std::vector<double> qs = { -1.75, -0.1, 0.3, 1.4, 3.9, 7.25 };
    for (double q : qs) {
        // only test inside the data domain to avoid extrap behavior
        if (q >= x[0] && q <= x[N-1]) {
            const double ytrue = f_quad(q);
            //std::cout << pol.interp(q) << ", " << ytrue << ", "<< q <<std::endl;
            CHECK(almost_eq(pol.interp(q), ytrue, 1e-10), "degree-2 polynomial should be exact on quadratic");
        }
    }
}

// 4) Cross-method consistency: methods agree closely on smooth data (cosine wave)
static double f_cos(double t) { return std::cos(0.1*t); }

TEST_CASE(Interp_Cross_Method_Consistency) {
    const std::uint64_t N = 60;
    utils::Vd x, y;
    make_uniform_xy_fun<double>(x, y, N, 0.0, 0.5, &f_cos);

    numerics::interp_linear<double>       lin(x,y);
    numerics::interp_polynomial<double>   pol(x,y, 3);  // small local degree
    numerics::interp_cubic_spline<double> spl(x,y);
    numerics::interp_rational<double>     rat(x,y, 3);
    numerics::interp_barycentric<double>  bar(x,y, 3);

    // sample a handful of interior points and require all methods to be mutually close
    std::vector<double> qs = { 1.25, 7.75, 12.1, 18.6, 22.75 };
    for (double q : qs) {
        // skip if q is outside just in case
        if (q < x[0] || q > x[N-1]) continue;

        double yl = lin.interp(q);
        double yp = pol.interp(q);
        double ys = spl.interp(q);
        double yr = rat.interp(q);
        double yb = bar.interp(q);

        //std::cout << "lin: " << yl << std::endl;
        //std::cout << "pol: " << yp << std::endl;
        //std::cout << "spl: " << ys << std::endl;
        //std::cout << "rat: " << yr << std::endl;
        //std::cout << "bar: " << yb << std::endl;

        // spline is usually the smoothest; use it as the anchor
        CHECK(almost_eq(yl, ys, 5e-4), "linear vs spline");
        CHECK(almost_eq(yp, ys, 5e-4), "polynomial vs spline");
        CHECK(almost_eq(yr, ys, 5e-4), "rational vs spline");
        CHECK(almost_eq(yb, ys, 5e-4), "barycentric vs spline");
    }
}