CubicSpline

src/interpolator/n/strategies.rs

@@ -225,3 +225,23 @@ where
         false
     }
 }
+
+impl<D> StrategyND<D> for CubicSpline<D::Elem>
+where
+    D: Data + RawDataClone + Clone,
+    D::Elem: Float + Debug,
+{
+    fn interpolate(
+        &self,
+        data: &InterpDataND<D>,
+        point: &[D::Elem],
+    ) -> Result<D::Elem, InterpolateError> {
+        let grids: Vec<ArrayView1<D::Elem>> = data.grid.iter().map(|g| g.view()).collect();
+        spline_eval_nd_recursive(&grids, data.values.view(), point, &self.boundary_conditions)
+    }
+
+    /// Returns `true`: the boundary cubic polynomials extend naturally.
+    fn allow_extrapolate(&self) -> bool {
+        true
+    }
+}

src/interpolator/n/tests.rs

@@ -1,5 +1,47 @@
 use super::*;
 
+#[test]
+fn test_cubic_spline() {
+    // f(x, y) = 2x + y: linear, reproduced exactly by any spline
+    let interp = InterpND::new(
+        vec![array![0., 1., 2.], array![0., 1., 2.]],
+        array![[0., 1., 2.], [2., 3., 4.], [4., 5., 6.]].into_dyn(),
+        strategy::CubicSpline::natural(),
+        Extrapolate::Enable,
+    )
+    .unwrap();
+    assert_approx_eq!(interp.interpolate(&[0.5, 0.5]).unwrap(), 1.5);
+    assert_approx_eq!(interp.interpolate(&[1., 1.]).unwrap(), 3.);
+    assert_approx_eq!(interp.interpolate(&[3., 1.]).unwrap(), 7.); // extrapolation
+}
+
+#[test]
+fn test_cubic_spline_knot_exactness() {
+    let interp = InterpND::new(
+        vec![array![0., 1., 2., 3.], array![0., 1., 2., 3.]],
+        array![
+            [0., 1., 4., 9.],
+            [1., 2., 5., 10.],
+            [4., 5., 8., 13.],
+            [9., 10., 13., 18.],
+        ]
+        .into_dyn(),
+        strategy::CubicSpline::natural(),
+        Extrapolate::Error,
+    )
+    .unwrap();
+    for i in 0..4usize {
+        for j in 0..4usize {
+            let xi = interp.data.grid[0][i];
+            let yj = interp.data.grid[1][j];
+            assert_approx_eq!(
+                interp.interpolate(&[xi, yj]).unwrap(),
+                interp.data.values[[i, j]]
+            );
+        }
+    }
+}
+
 #[test]
 fn test_linear_0d() {
     let interp = InterpND::new(

src/interpolator/one/strategies.rs

@@ -114,3 +114,48 @@ where
         false
     }
 }
+
+impl<D> Strategy1D<D> for CubicSpline<D::Elem>
+where
+    D: Data + RawDataClone + Clone,
+    D::Elem: Float + Debug,
+{
+    /// Computes and caches `M[0..=n]` using [`compute_m`] for the configured BC.
+    fn init(&mut self, data: &InterpData1D<D>) -> Result<(), ValidateError> {
+        let new_m = compute_m(data.grid[0].view(), data.values.view(), self.bc_for_dim(0))?;
+        self.m = new_m;
+        Ok(())
+    }
+
+    fn interpolate(
+        &self,
+        data: &InterpData1D<D>,
+        point: &[D::Elem; 1],
+    ) -> Result<D::Elem, InterpolateError> {
+        let grid = data.grid[0].view();
+        // For extrapolation, clamp to the boundary interval.
+        let i = if &point[0] < grid.first().unwrap() {
+            0
+        } else if &point[0] > grid.last().unwrap() {
+            grid.len() - 2
+        } else {
+            find_nearest_index(grid, &point[0])
+        };
+        let two = D::Elem::one() + D::Elem::one();
+        let six = two + two + two;
+        let h = grid[i + 1] - grid[i];
+        let dx = point[0] - grid[i]; // t - x[i]
+        let dx_r = h - dx; // x[i+1] - t
+        let six_h = six * h;
+        let h2_over_six = h * h / six;
+        Ok(self.m[i] * dx_r * dx_r * dx_r / six_h
+            + self.m[i + 1] * dx * dx * dx / six_h
+            + (data.values[i] - self.m[i] * h2_over_six) * dx_r / h
+            + (data.values[i + 1] - self.m[i + 1] * h2_over_six) * dx / h)
+    }
+
+    /// Returns `true`: the boundary cubic polynomials extend naturally.
+    fn allow_extrapolate(&self) -> bool {
+        true
+    }
+}

src/interpolator/one/tests.rs

@@ -1,5 +1,60 @@
 use super::*;
 
+#[test]
+fn test_cubic_spline() {
+    // Linear data: any spline reproduces it exactly (all second derivatives = 0)
+    let interp = Interp1D::new(
+        array![0., 1., 2., 3.],
+        array![1., 3., 5., 7.], // f(x) = 2x + 1
+        strategy::CubicSpline::not_a_knot(),
+        Extrapolate::Enable,
+    )
+    .unwrap();
+    // Knot values
+    let x = interp.data.grid[0].clone();
+    for (i, xi) in x.iter().enumerate() {
+        assert_approx_eq!(interp.interpolate(&[*xi]).unwrap(), interp.data.values[i]);
+    }
+    // Midpoints
+    assert_approx_eq!(interp.interpolate(&[0.5]).unwrap(), 2.0);
+    assert_approx_eq!(interp.interpolate(&[1.5]).unwrap(), 4.0);
+    assert_approx_eq!(interp.interpolate(&[2.5]).unwrap(), 6.0);
+    // Extrapolation via boundary polynomials
+    assert_approx_eq!(interp.interpolate(&[-1.0]).unwrap(), -1.0);
+    assert_approx_eq!(interp.interpolate(&[4.0]).unwrap(), 9.0);
+}
+
+#[test]
+fn test_cubic_spline_knot_exactness() {
+    // Values at all knots must be reproduced exactly regardless of data shape
+    let interp = Interp1D::new(
+        array![0., 1., 2., 3., 4.],
+        array![0., 1., 4., 9., 16.], // f(x) = x^2
+        strategy::CubicSpline::not_a_knot(),
+        Extrapolate::Error,
+    )
+    .unwrap();
+    let x = interp.data.grid[0].clone();
+    for (i, xi) in x.iter().enumerate() {
+        assert_approx_eq!(interp.interpolate(&[*xi]).unwrap(), interp.data.values[i]);
+    }
+}
+
+#[test]
+fn test_cubic_spline_two_points() {
+    // Degenerate case: 2 points → degenerates to linear interpolation.
+    // Uses Natural BC since NotAKnot requires ≥ 3 points.
+    let interp = Interp1D::new(
+        array![0., 1.],
+        array![0., 2.],
+        strategy::CubicSpline::natural(),
+        Extrapolate::Enable,
+    )
+    .unwrap();
+    assert_approx_eq!(interp.interpolate(&[0.5]).unwrap(), 1.0);
+    assert_approx_eq!(interp.interpolate(&[2.0]).unwrap(), 4.0); // extrapolation
+}
+
 #[test]
 fn test_invalid_args() {
     let interp = Interp1D::new(

src/interpolator/three/strategies.rs

@@ -241,3 +241,28 @@ where
         false
     }
 }
+
+impl<D> Strategy3D<D> for CubicSpline<D::Elem>
+where
+    D: Data + RawDataClone + Clone,
+    D::Elem: Float + Debug,
+{
+    fn interpolate(
+        &self,
+        data: &InterpData3D<D>,
+        point: &[D::Elem; 3],
+    ) -> Result<D::Elem, InterpolateError> {
+        let grids: Vec<ArrayView1<D::Elem>> = data.grid.iter().map(|g| g.view()).collect();
+        spline_eval_nd_recursive(
+            &grids,
+            data.values.view().into_dyn(),
+            point,
+            &self.boundary_conditions,
+        )
+    }
+
+    /// Returns `true`: the boundary cubic polynomials extend naturally.
+    fn allow_extrapolate(&self) -> bool {
+        true
+    }
+}

src/interpolator/three/tests.rs

@@ -1,5 +1,60 @@
 use super::*;
 
+#[test]
+fn test_cubic_spline() {
+    // f(x, y, z) = x + 2y + 3z: linear, reproduced exactly by any spline
+    let interp = Interp3D::new(
+        array![0., 1., 2.],
+        array![0., 1., 2.],
+        array![0., 1., 2.],
+        array![
+            [[0., 3., 6.], [2., 5., 8.], [4., 7., 10.]],
+            [[1., 4., 7.], [3., 6., 9.], [5., 8., 11.]],
+            [[2., 5., 8.], [4., 7., 10.], [6., 9., 12.]],
+        ],
+        strategy::CubicSpline::natural(),
+        Extrapolate::Enable,
+    )
+    .unwrap();
+    // Knots
+    assert_approx_eq!(interp.interpolate(&[1., 1., 1.]).unwrap(), 6.);
+    // Midpoints
+    assert_approx_eq!(interp.interpolate(&[0.5, 0.5, 0.5]).unwrap(), 3.);
+    assert_approx_eq!(interp.interpolate(&[1.0, 0.5, 1.0]).unwrap(), 5.);
+    // Extrapolation
+    assert_approx_eq!(interp.interpolate(&[3., 1., 1.]).unwrap(), 8.);
+}
+
+#[test]
+fn test_cubic_spline_knot_exactness() {
+    let interp = Interp3D::new(
+        array![0., 1., 2.],
+        array![0., 1., 2.],
+        array![0., 1., 2.],
+        array![
+            [[0., 1., 2.], [3., 4., 5.], [6., 7., 8.]],
+            [[9., 10., 11.], [12., 13., 14.], [15., 16., 17.]],
+            [[18., 19., 20.], [21., 22., 23.], [24., 25., 26.]],
+        ],
+        strategy::CubicSpline::natural(),
+        Extrapolate::Error,
+    )
+    .unwrap();
+    let x = interp.data.grid[0].clone();
+    let y = interp.data.grid[1].clone();
+    let z = interp.data.grid[2].clone();
+    for (i, xi) in x.iter().enumerate() {
+        for (j, yj) in y.iter().enumerate() {
+            for (k, zk) in z.iter().enumerate() {
+                assert_approx_eq!(
+                    interp.interpolate(&[*xi, *yj, *zk]).unwrap(),
+                    interp.data.values[[i, j, k]]
+                );
+            }
+        }
+    }
+}
+
 #[test]
 fn test_linear() {
     let interp = Interp3D::new(

src/interpolator/two/strategies.rs

@@ -172,3 +172,28 @@ where
         false
     }
 }
+
+impl<D> Strategy2D<D> for CubicSpline<D::Elem>
+where
+    D: Data + RawDataClone + Clone,
+    D::Elem: Float + Debug,
+{
+    fn interpolate(
+        &self,
+        data: &InterpData2D<D>,
+        point: &[D::Elem; 2],
+    ) -> Result<D::Elem, InterpolateError> {
+        let grids: Vec<ArrayView1<D::Elem>> = data.grid.iter().map(|g| g.view()).collect();
+        spline_eval_nd_recursive(
+            &grids,
+            data.values.view().into_dyn(),
+            point,
+            &self.boundary_conditions,
+        )
+    }
+
+    /// Returns `true`: the boundary cubic polynomials extend naturally.
+    fn allow_extrapolate(&self) -> bool {
+        true
+    }
+}

src/interpolator/two/tests.rs

@@ -1,5 +1,54 @@
 use super::*;
 
+#[test]
+fn test_cubic_spline() {
+    // f(x, y) = 2x + y: linear in both dims, reproduced exactly by any spline
+    let interp = Interp2D::new(
+        array![0., 1., 2.],
+        array![0., 1., 2.],
+        array![[0., 1., 2.], [2., 3., 4.], [4., 5., 6.]],
+        strategy::CubicSpline::natural(),
+        Extrapolate::Enable,
+    )
+    .unwrap();
+    // Knots
+    assert_approx_eq!(interp.interpolate(&[1., 1.]).unwrap(), 3.);
+    assert_approx_eq!(interp.interpolate(&[2., 0.]).unwrap(), 4.);
+    // Midpoints
+    assert_approx_eq!(interp.interpolate(&[0.5, 0.5]).unwrap(), 1.5);
+    assert_approx_eq!(interp.interpolate(&[1.5, 1.0]).unwrap(), 4.);
+    // Extrapolation
+    assert_approx_eq!(interp.interpolate(&[3., 1.]).unwrap(), 7.);
+    assert_approx_eq!(interp.interpolate(&[1., 3.]).unwrap(), 5.);
+}
+
+#[test]
+fn test_cubic_spline_knot_exactness() {
+    let interp = Interp2D::new(
+        array![0., 1., 2., 3.],
+        array![0., 1., 2., 3.],
+        array![
+            [0., 1., 4., 9.],
+            [1., 2., 5., 10.],
+            [4., 5., 8., 13.],
+            [9., 10., 13., 18.],
+        ], // f(x, y) = x^2 + y
+        strategy::CubicSpline::not_a_knot(),
+        Extrapolate::Error,
+    )
+    .unwrap();
+    let x = interp.data.grid[0].clone();
+    let y = interp.data.grid[1].clone();
+    for (i, xi) in x.iter().enumerate() {
+        for (j, yj) in y.iter().enumerate() {
+            assert_approx_eq!(
+                interp.interpolate(&[*xi, *yj]).unwrap(),
+                interp.data.values[[i, j]]
+            );
+        }
+    }
+}
+
 #[test]
 fn test_linear() {
     let interp = Interp2D::new(

src/lib.rs

@@ -42,7 +42,9 @@ pub use interpolator::data;
 pub(crate) use interpolator::data::*;
 
 pub(crate) use error::*;
+pub(crate) use strategy::spline::*;
 pub(crate) use strategy::traits::*;
+pub(crate) use strategy::utils::*;
 
 pub(crate) use std::fmt::Debug;