155 lines
5.9 KiB
C++
155 lines
5.9 KiB
C++
#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");
|
|
}
|
|
} |