feat: Add ND Hybercube array

This functionality is required for ND cubic interpolation

Signed-off-by: Sietze van Buuren <s.van.buuren@gmail.com>

Author: swvanbuuren

include/arrayn.hpp

@@ -0,0 +1,53 @@
+#pragma once
+
+#include <mdspan/mdspan.hpp>
+#include <array>
+#include <cstddef>
+
+
+namespace cip {
+
+
+constexpr std::size_t power(std::size_t base, std::size_t exp) {
+    return exp == 0 ? 1 : base * power(base, exp - 1);
+}
+
+
+template <std::size_t Value, std::size_t... I>
+constexpr std::extents<std::size_t, ((void)I, Value)...>
+make_extents(std::index_sequence<I...>) {
+    return {};
+}
+
+
+template<typename T, std::size_t N, std::size_t Size>
+class NDHyperArray {
+    static constexpr std::size_t total_size = power(Size, N);
+    using Array = std::array<T, total_size>;
+    using ExtentsType = decltype(make_extents<Size>(std::make_index_sequence<N>{}));
+    using Mdspan  = std::mdspan<T, ExtentsType>;
+public:
+    NDHyperArray()
+      : mdspan(data.data())
+    {}
+
+    ~NDHyperArray() = default;
+
+    template<typename... Index>
+    T& operator()(Index... idx) {
+        return mdspan(idx...);
+    }
+
+    template<typename... Index>
+    const T& operator()(Index... idx) const {
+        return mdspan(idx...);
+    }
+
+private:
+    Array data;
+    Mdspan mdspan;
+
+};
+
+
+} // namespace cip

include/cubinterpp.hpp

@@ -3,6 +3,7 @@
 #include "utils.hpp"
 #include "slopes.hpp"
 #include "vectorn.hpp"
+#include "arrayn.hpp"
 #include "linear_interp.hpp"
 #include "cubic_spline.hpp"
 #include "cubic_splines_1d.hpp"

tests/CMakeLists.txt

@@ -22,6 +22,7 @@ FetchContent_Declare(
 FetchContent_MakeAvailable(googletest)
 
 cubinterpp_add_test(test_vectorn)
+cubinterpp_add_test(test_arrayn)
 cubinterpp_add_test(test_linear_interp)
 cubinterpp_add_test(test_cubic_splines_1d)
 cubinterpp_add_test(test_cubic_splines_2d)

tests/test_arrayn.cpp

@@ -0,0 +1,49 @@
+#include <gtest/gtest.h>
+#include <cubinterpp.hpp>
+
+
+// Test a 1-dimensional NDHyperArray.
+TEST(NDHyperArrayTest, OneDArray) {
+    cip::NDHyperArray<int, 1, 4> arr;
+    arr(0) = 10;
+    arr(1) = 20;
+    arr(2) = 30;
+    arr(3) = 40;
+    EXPECT_EQ(arr(0), 10);
+    EXPECT_EQ(arr(1), 20);
+    EXPECT_EQ(arr(2), 30);
+    EXPECT_EQ(arr(3), 40);
+}
+
+
+// Test a 2-dimensional NDHyperArray.
+TEST(NDHyperArrayTest, TwoDArray) {
+    cip::NDHyperArray<int, 2, 4> arr;
+    arr(0, 0) = 1;
+    arr(0, 1) = 2;
+    arr(1, 0) = 3;
+    arr(3, 3) = 4;
+    EXPECT_EQ(arr(0, 0), 1);
+    EXPECT_EQ(arr(0, 1), 2);
+    EXPECT_EQ(arr(1, 0), 3);
+    EXPECT_EQ(arr(3, 3), 4);
+}
+
+
+// Test a 3-dimensional NDHyperArray.
+TEST(NDHyperArrayTest, ThreeDArray) {
+    cip::NDHyperArray<double, 3, 4> arr;
+    arr(1, 2, 3) = 3.14159;
+    arr(0, 0, 0) = 2.71828;
+    EXPECT_NEAR(arr(1, 2, 3), 3.14159, 1e-6);
+    EXPECT_NEAR(arr(0, 0, 0), 2.71828, 1e-6);
+}
+
+
+// Test const access to ensure the const overload of operator() works.
+TEST(NDHyperArrayTest, ConstAccess) {
+    cip::NDHyperArray<int, 2, 4> arr;
+    arr(2, 2) = 99;
+    const cip::NDHyperArray<int, 2, 4>& constArr = arr;
+    EXPECT_EQ(constArr(2, 2), 99);
+}