Generalize assertion helper for N dimensions

include/linear_interp.hpp

@@ -280,6 +280,16 @@ class LinearInterp3D : public LinearInterpND<T, 3> {
     ~LinearInterp3D() { }
 };
 
+template <typename T>
+class LinearInterp4D : public LinearInterpND<T, 4> {
+    using Vector = std::vector<T>;
+    using Vector4 = cip::VectorN<T, 4>;
+public:
+    explicit LinearInterp4D(const Vector &x, const Vector &y, const Vector &z, const Vector &w, const Vector4 &f)
+    : LinearInterpND<T, 4>(f, x, y, z, w)
+    {}
 
+    ~LinearInterp4D() { }
+};
 
 } // namespace cip

tests/assertion_helpers.hpp

@@ -1,63 +1,147 @@
 #include <vector>
 #include <gtest/gtest.h>
-
+#include <array>
+#include <tuple>
+#include <type_traits>
+#include <sstream>
+#include <functional>
+#include <utility>
 
 namespace cip {
 
-
 using Vector = std::vector<double>;
 using Vector2 = std::vector<Vector>;
 using Vector3 = std::vector<Vector2>;
+using Vector4 = std::vector<Vector3>;
 
 constexpr double TOLERANCE_DEFAULT = 5.0e-12;
 
-template <typename T>
-testing::AssertionResult Interp1DAssertions(Vector x, Vector f, Vector x_fine, Vector f_fine, double tol=TOLERANCE_DEFAULT) {
-    T interp(x, f);
-    for ( auto i = 0; i < x_fine.size(); i++ ) {
-        auto val = interp.eval(x_fine[i]);
-        if (!testing::internal::DoubleNearPredFormat("expected", "actual", "tolerance", val, f_fine[i], tol)  ) {
-            return testing::AssertionFailure()
-                << "for x = " << x_fine[i] << " expected " << f_fine[i] << " but got " << val;
+// Helper class for N-dimensional interpolation assertions
+template <typename T, size_t N, typename FuncValT>
+class InterpNDAssertions {
+private:
+    // Helper to access nested vector value at specific indices
+    template <typename VecT>
+    static double get_value(const VecT& val) {
+        return val;
+    }
+
+    template <typename VecT, typename... Indices>
+    static double get_value(const VecT& vec, size_t idx, Indices... rest) {
+        return get_value(vec[idx], rest...);
+    }
+
+    // Helper for dimension name (x, y, z, w, etc.)
+    static constexpr char dim_name(size_t i) {
+        return i < 3 ? 'x' + i : 'w' + (i - 3);
+    }
+
+    // Recursive template to build nested loops over N dimensions
+    template <size_t Dim>
+    static testing::AssertionResult check_points(
+            const T& interp,
+            const std::array<Vector, N>& fine_grid_vectors,
+            const FuncValT& f_fine,
+            std::array<size_t, N>& indices,
+            std::array<double, N>& coords,
+            double tol) {
+
+        // Base case: we've set up all N dimension indices
+        if constexpr (Dim == N) {
+            // Call eval with the coordinates using apply
+            double val = std::apply([&interp](auto... args) { 
+                return interp.eval(args...); 
+            }, coords);
+
+            // Get expected value from nested vector using indices
+            double expected = std::apply([&f_fine](auto... args) { 
+                return get_value(f_fine, args...); 
+            }, indices);
+
+            // Check if the values match within tolerance
+            if (!testing::internal::DoubleNearPredFormat("expected", "actual", "tolerance", val, expected, tol)) {
+                testing::AssertionResult failure = testing::AssertionFailure();
+
+                // Format dimensions for error message (x=1.23, y=4.56, etc.)
+                for (size_t i = 0; i < N; ++i) {
+                    failure << (i == 0 ? "for " : ", ") << dim_name(i) << " = " << coords[i];
+                }
+
+                failure << " expected " << expected << " but got " << val;
+                return failure;
+            }
+
+            return testing::AssertionSuccess();
+        }
+        // Recursive case: loop over the current dimension
+        else {
+            for (size_t i = 0; i < fine_grid_vectors[Dim].size(); ++i) {
+                indices[Dim] = i;
+                coords[Dim] = fine_grid_vectors[Dim][i];
+
+                // Recursively process the next dimension
+                auto result = check_points<Dim + 1>(interp, fine_grid_vectors, f_fine, indices, coords, tol);
+                if (!result) {
+                    return result;
+                }
+            }
+            return testing::AssertionSuccess();
         }
     }
-    return testing::AssertionSuccess();
-}
 
+public:
+    // Main assertion method
+    static testing::AssertionResult test(
+            const std::array<Vector, N>& grid_vectors,
+            const FuncValT& f,
+            const std::array<Vector, N>& fine_grid_vectors,
+            const FuncValT& f_fine,
+            double tol = TOLERANCE_DEFAULT) {
+
+        // Create the interpolator with a simple lambda + std::apply
+        T interp = std::apply([&f](const auto&... grid_args) {
+            return T(grid_args..., f);
+        }, grid_vectors);
+
+        // Set up indices and coordinates arrays
+        std::array<size_t, N> indices{};
+        std::array<double, N> coords{};
+
+        // Start the recursive dimension traversal
+        return check_points<0>(interp, fine_grid_vectors, f_fine, indices, coords, tol);
+    }
+};
 
+// Convenience wrapper functions to maintain backward compatibility
+template <typename T>
+testing::AssertionResult Interp1DAssertions(Vector x, Vector f, Vector x_fine, Vector f_fine, double tol=TOLERANCE_DEFAULT) {
+    std::array<Vector, 1> grid_vectors = {x};
+    std::array<Vector, 1> fine_grid_vectors = {x_fine};
+    return InterpNDAssertions<T, 1, Vector>::test(grid_vectors, f, fine_grid_vectors, f_fine, tol);
+}
 
 template <typename T>
 testing::AssertionResult Interp2DAssertions(const Vector &x, const Vector &y, const Vector2 &f, const Vector &x_fine, const Vector &y_fine, const Vector2 &f_fine, double tol=TOLERANCE_DEFAULT) {
-    T interp2(x, y, f);
-    for ( auto i = 0; i < x_fine.size(); ++i ) {
-        for ( auto j = 0; j < y_fine.size(); ++j ) {
-            auto val = interp2.eval(x_fine[i], y_fine[j]);
-            if (!testing::internal::DoubleNearPredFormat("expected", "actual", "tolerance", val, f_fine[i][j], tol)  ) {
-                return testing::AssertionFailure()
-                    << "for x = " << x_fine[i] << ", y = " << y_fine[j] << " expected " << f_fine[i][j] << " but got " << val;
-            }
-        }
-    }
-    return testing::AssertionSuccess();
+    std::array<Vector, 2> grid_vectors = {x, y};
+    std::array<Vector, 2> fine_grid_vectors = {x_fine, y_fine};
+    return InterpNDAssertions<T, 2, Vector2>::test(grid_vectors, f, fine_grid_vectors, f_fine, tol);
 }
 
-
 template <typename T>
 testing::AssertionResult Interp3DAssertions(const Vector &x, const Vector &y, const Vector &z, const Vector3 &f, const Vector &x_fine, const Vector &y_fine, const Vector &z_fine, const Vector3 &f_fine, double tol=TOLERANCE_DEFAULT) {
-    T interp3(x, y, z, f);
-    for ( auto i = 0; i < x_fine.size(); ++i ) {
-        for ( auto j = 0; j < y_fine.size(); ++j ) {
-            for ( auto k = 0; k < y_fine.size(); ++k ) {
-                auto val = interp3.eval(x_fine[i], y_fine[j], z_fine[k]);
-                if (!testing::internal::DoubleNearPredFormat("expected", "actual", "tolerance", val, f_fine[i][j][k], tol)  ) {
-                    return testing::AssertionFailure()
-                        << "for x = " << x_fine[i] << ", y = " << y_fine[j] << ", z = " << z_fine[j] << " expected " << f_fine[i][j][k] << " but got " << val;
-                }
-            }
-        }
-    }
-    return testing::AssertionSuccess();
+    std::array<Vector, 3> grid_vectors = {x, y, z};
+    std::array<Vector, 3> fine_grid_vectors = {x_fine, y_fine, z_fine};
+    return InterpNDAssertions<T, 3, Vector3>::test(grid_vectors, f, fine_grid_vectors, f_fine, tol);
 }
 
+template <typename T>
+testing::AssertionResult Interp4DAssertions(const Vector &x, const Vector &y, const Vector &z, const Vector &w, 
+                                           const Vector4 &f, const Vector &x_fine, const Vector &y_fine, 
+                                           const Vector &z_fine, const Vector &w_fine, const Vector4 &f_fine, 
+                                           double tol=TOLERANCE_DEFAULT) {
+    std::array<Vector, 4> grid_vectors = {x, y, z, w};
+    std::array<Vector, 4> fine_grid_vectors = {x_fine, y_fine, z_fine, w_fine};
+    return InterpNDAssertions<T, 4, Vector4>::test(grid_vectors, f, fine_grid_vectors, f_fine, tol);
+}
 
 } // namespace cip

tests/test_linear_interp.cpp

@@ -7,6 +7,7 @@
 using VectorN1 = cip::VectorN<double, 1>;
 using VectorN2 = cip::VectorN<double, 2>;
 using VectorN3 = cip::VectorN<double, 3>;
+using VectorN4 = cip::VectorN<double, 4>;
 using Span = std::span<const double>;
 using Pr = std::pair<size_t, size_t>;
 
@@ -197,3 +198,164 @@ TEST(TestInterp3D, test_linear_interp_3d) {
 ASSERT_TRUE(cip::Interp3DAssertions<cip::LinearInterp3D<double>>(x, y, z, f, x_fine, y_fine, z_fine, f_fine));    
 
 }
+
+/**
+ * Test for 4D linear interpolation
+ * 
+ * This test demonstrates using the InterpNDAssertions class with N=4 dimensions.
+ * We use a simple 4D linear function f(x,y,z,w) = a*x + b*y + c*z + d*w + e
+ * to generate the test data, which should be perfectly interpolated by a linear interpolator.
+ */
+TEST(TestInterp4D, test_linear_interp_4d) {
+    // Create grid vectors for each dimension
+    cip::Vector x = { 0.0, 1.0, 2.0 };
+    cip::Vector y = { 0.0, 1.0, 2.0 };
+    cip::Vector z = { 0.0, 1.0, 2.0 };
+    cip::Vector w = { 0.0, 1.0, 2.0 };
+
+    // Define a linear function coefficients: f(x,y,z,w) = 2*x + 3*y - z + 0.5*w + 1
+    const double a = 2.0;
+    const double b = 3.0;
+    const double c = -1.0;
+    const double d = 0.5;
+    const double e = 1.0;
+
+    // Create 4D function values array
+    cip::Vector4 f;
+    f.resize(x.size());
+    for (size_t i = 0; i < x.size(); ++i) {
+        f[i].resize(y.size());
+        for (size_t j = 0; j < y.size(); ++j) {
+            f[i][j].resize(z.size());
+            for (size_t k = 0; k < z.size(); ++k) {
+                f[i][j][k].resize(w.size());
+                for (size_t l = 0; l < w.size(); ++l) {
+                    // Function: f(x,y,z,w) = 2*x + 3*y - z + 0.5*w + 1
+                    f[i][j][k][l] = a*x[i] + b*y[j] + c*z[k] + d*w[l] + e;
+                }
+            }
+        }
+    }
+
+    // Fine grid for testing interpolation
+    cip::Vector x_fine = { 0.0, 0.5, 1.0, 1.5, 2.0 };
+    cip::Vector y_fine = { 0.0, 0.5, 1.0, 1.5, 2.0 };
+    cip::Vector z_fine = { 0.0, 0.5, 1.0, 1.5, 2.0 };
+    cip::Vector w_fine = { 0.0, 0.5, 1.0, 1.5, 2.0 };
+
+    // Expected values at fine grid points
+    cip::Vector4 f_fine;
+    f_fine.resize(x_fine.size());
+    for (size_t i = 0; i < x_fine.size(); ++i) {
+        f_fine[i].resize(y_fine.size());
+        for (size_t j = 0; j < y_fine.size(); ++j) {
+            f_fine[i][j].resize(z_fine.size());
+            for (size_t k = 0; k < z_fine.size(); ++k) {
+                f_fine[i][j][k].resize(w_fine.size());
+                for (size_t l = 0; l < w_fine.size(); ++l) {
+                    // Function: f(x,y,z,w) = 2*x + 3*y - z + 0.5*w + 1
+                    f_fine[i][j][k][l] = a*x_fine[i] + b*y_fine[j] + c*z_fine[k] + d*w_fine[l] + e;
+                }
+            }
+        }
+    }
+
+    // Test the 4D interpolation using our InterpNDAssertions class
+    ASSERT_TRUE(cip::Interp4DAssertions<cip::LinearInterp4D<double>>(
+        x, y, z, w, f, x_fine, y_fine, z_fine, w_fine, f_fine));
+}
+
+// Direct test using the InterpNDAssertions template with N=4
+TEST(TestInterp4D, test_direct_nd_assertions) {
+    // Create grid vectors for each dimension
+    cip::Vector x = { 0.0, 1.0, 2.0 };
+    cip::Vector y = { 0.0, 1.0, 2.0 };
+    cip::Vector z = { 0.0, 1.0, 2.0 };
+    cip::Vector w = { 0.0, 1.0, 2.0 };
+
+    // Create 4D function values array
+    // Function is f(x,y,z,w) = x + y + z + w
+    cip::Vector4 f;
+    f.resize(x.size());
+    for (size_t i = 0; i < x.size(); ++i) {
+        f[i].resize(y.size());
+        for (size_t j = 0; j < y.size(); ++j) {
+            f[i][j].resize(z.size());
+            for (size_t k = 0; k < z.size(); ++k) {
+                f[i][j][k].resize(w.size());
+                for (size_t l = 0; l < w.size(); ++l) {
+                    // Function: f(x,y,z,w) = x + y + z + w
+                    f[i][j][k][l] = x[i] + y[j] + z[k] + w[l];
+                }
+            }
+        }
+    }
+
+    // Fine grid for testing interpolation
+    cip::Vector x_fine = { 0.0, 0.5, 1.0, 1.5, 2.0 };
+    cip::Vector y_fine = { 0.0, 0.5, 1.0, 1.5, 2.0 };
+    cip::Vector z_fine = { 0.0, 0.5, 1.0, 1.5, 2.0 };
+    cip::Vector w_fine = { 0.0, 0.5, 1.0, 1.5, 2.0 };
+
+    // Expected values at fine grid points
+    cip::Vector4 f_fine;
+    f_fine.resize(x_fine.size());
+    for (size_t i = 0; i < x_fine.size(); ++i) {
+        f_fine[i].resize(y_fine.size());
+        for (size_t j = 0; j < y_fine.size(); ++j) {
+            f_fine[i][j].resize(z_fine.size());
+            for (size_t k = 0; k < z_fine.size(); ++k) {
+                f_fine[i][j][k].resize(w_fine.size());
+                for (size_t l = 0; l < w_fine.size(); ++l) {
+                    // Function: f(x,y,z,w) = x + y + z + w
+                    f_fine[i][j][k][l] = x_fine[i] + y_fine[j] + z_fine[k] + w_fine[l];
+                }
+            }
+        }
+    }
+
+    // Pack grid vectors and fine grid vectors into arrays
+    std::array<cip::Vector, 4> grid_vectors = {x, y, z, w};
+    std::array<cip::Vector, 4> fine_grid_vectors = {x_fine, y_fine, z_fine, w_fine};
+
+    // Test the 4D interpolation directly using the InterpNDAssertions class
+    auto result = cip::InterpNDAssertions<cip::LinearInterp4D<double>, 4, cip::Vector4>::test(
+        grid_vectors, f, fine_grid_vectors, f_fine);
+    ASSERT_TRUE(result);
+}
+
+// Test a non-linear function with 4D interpolation
+TEST(TestInterp4D, test_nonlinear_4d) {
+    // Create grid vectors for each dimension
+    cip::Vector x = { 0.0, 1.0, 2.0 };
+    cip::Vector y = { 0.0, 1.0, 2.0 };
+    cip::Vector z = { 0.0, 1.0, 2.0 };
+    cip::Vector w = { 0.0, 1.0, 2.0 };
+
+    // Create 4D function values array for a quadratic function
+    // Function is f(x,y,z,w) = x^2 + y^2 + z^2 + w^2
+    cip::Vector4 f;
+    f.resize(x.size());
+    for (size_t i = 0; i < x.size(); ++i) {
+        f[i].resize(y.size());
+        for (size_t j = 0; j < y.size(); ++j) {
+            f[i][j].resize(z.size());
+            for (size_t k = 0; k < z.size(); ++k) {
+                f[i][j][k].resize(w.size());
+                for (size_t l = 0; l < w.size(); ++l) {
+                    // Quadratic function
+                    f[i][j][k][l] = x[i]*x[i] + y[j]*y[j] + z[k]*z[k] + w[l]*w[l];
+                }
+            }
+        }
+    }
+
+    // Create interpolator
+    cip::LinearInterp4D<double> interp(x, y, z, w, f);
+
+    // Test a few specific points without using fine grid
+    // For interpolation at grid points, result should be exact
+    EXPECT_DOUBLE_EQ(interp.eval(0.0, 0.0, 0.0, 0.0), 0.0);
+    EXPECT_DOUBLE_EQ(interp.eval(1.0, 1.0, 1.0, 1.0), 4.0);
+    EXPECT_DOUBLE_EQ(interp.eval(2.0, 2.0, 2.0, 2.0), 16.0);
+}