CubicC1

CHANGELOG.md

@@ -8,6 +8,17 @@ bump the minor version (`0.x` -> `0.(x+1)`), other changes bump the patch versio
 
 ## [Unreleased]
 
+### Added
+- `strategy::CubicC1`: a C¹ local cubic Hermite spline strategy (finite-difference
+  derivative estimate, no global solve). Cheaper to build than `CubicC2`, matching the
+  local/uncached recipe LHAPDF-style consumers (e.g. neopdf) use by default at 3D+, but
+  not aiming for bit-for-bit LHAPDF parity. `derivative_mode` carries the
+  derivative-estimate method (`FiniteDifference` for now, `#[non_exhaustive]` for future
+  monotonicity-preserving variants); `cache_mode` chooses between precomputing the full
+  corner-derivative tensor at `init()` (`Full`, the default, same mechanism as
+  `CubicC2`) or deriving it fresh from a bounded local neighborhood on every query
+  (`None`) at 2-D and above. Closes #55.
+
 ## [0.11.1] - 2026-08-22
 
 ### Added

src/interpolator/n/strategies.rs

@@ -253,6 +253,51 @@ where
     }
 }
 
+impl<D> StrategyND<D> for CubicC1<D::Elem>
+where
+    D: Data + RawDataClone + Clone,
+    D::Elem: Float + Debug,
+{
+    /// Precomputes the full corner-derivative tensor under [`CubicC1CacheMode::Full`]
+    /// (the default); under [`CubicC1CacheMode::None`], only validates.
+    fn init(&mut self, data: &InterpDataNDBase<D>) -> Result<(), ValidateError> {
+        if data.ndim() == 0 {
+            return Ok(());
+        }
+        if self.cache_mode == CubicC1CacheMode::Full {
+            let data_view = data.view();
+            self.cache = compute_corner_cache_fd(&data_view.grid, data_view.values);
+        }
+        Ok(())
+    }
+
+    fn interpolate(
+        &self,
+        data: &InterpDataNDBase<D>,
+        point: &[D::Elem],
+    ) -> Result<D::Elem, InterpolateError> {
+        if data.ndim() == 0 {
+            return data.values.first().copied().ok_or_else(|| {
+                InterpolateError::Other("internal: 0-D interpolation data has no value".into())
+            });
+        }
+        let grids: Vec<ArrayView1<D::Elem>> = data.grid.iter().map(|g| g.view()).collect();
+        Ok(match self.cache_mode {
+            CubicC1CacheMode::Full => {
+                evaluate_spline_corner_cached(&grids, self.cache.view(), point)
+            }
+            CubicC1CacheMode::None => {
+                evaluate_spline_corner_local(&grids, data.values.view(), point)
+            }
+        })
+    }
+
+    /// Returns `true`: the boundary Hermite patch extends naturally.
+    fn allow_extrapolate(&self) -> bool {
+        true
+    }
+}
+
 impl<D, S> StrategyND<D> for GridTransform<D::Elem, S>
 where
     D: Data + RawDataClone + Clone,

src/interpolator/n/tests.rs

@@ -66,6 +66,18 @@ fn test_cubic_spline_0d() {
     assert_eq!(interp.interpolate(&[]).unwrap(), 0.5);
 }
 
+#[test]
+fn test_cubic_c1_0d() {
+    let interp = InterpND::new(
+        vec![array![]],
+        array![0.5].into_dyn(),
+        strategy::CubicC1::default(),
+        Extrapolate::Error,
+    )
+    .unwrap();
+    assert_eq!(interp.interpolate(&[]).unwrap(), 0.5);
+}
+
 #[test]
 fn test_cubic_c2_periodic_outer_axis() {
     // Smoke test: `Periodic` on a non-innermost axis still interpolates successfully.

src/interpolator/one/strategies.rs

@@ -150,5 +150,36 @@ where
     }
 }
 
+impl<D> Strategy1D<D> for CubicC1<D::Elem>
+where
+    D: Data + RawDataClone + Clone,
+    D::Elem: Float + Debug,
+{
+    /// Caches the finite-difference derivative vector. `cache_mode` is ignored here:
+    /// the cache is already O(1) regardless.
+    fn init(&mut self, data: &InterpData1DBase<D>) -> Result<(), ValidateError> {
+        self.cache = compute_fd_cache(data.grid[0].view(), data.values.view());
+        Ok(())
+    }
+
+    fn interpolate(
+        &self,
+        data: &InterpData1DBase<D>,
+        point: &[D::Elem; 1],
+    ) -> Result<D::Elem, InterpolateError> {
+        evaluate_hermite_1d_cached(
+            data.grid[0].view(),
+            data.values.view(),
+            self.cache.view(),
+            point[0],
+        )
+    }
+
+    /// Returns `true`: the boundary Hermite segment extends naturally.
+    fn allow_extrapolate(&self) -> bool {
+        true
+    }
+}
+
 grid_transform_strategy_impl!(Strategy1D, InterpData1DBase, InterpData1DView, 1);
 values_transform_strategy_impl!(Strategy1D, InterpData1DBase, InterpData1DView, Ix1, 1);

src/interpolator/one/tests.rs

@@ -299,6 +299,66 @@ fn test_cubic_c2_clamped_uses_given_derivative() {
     );
 }
 
+#[test]
+fn test_cubic_c1_linear_exact() {
+    // Linear data: finite differences recover the exact constant slope with no error
+    // term, so the Hermite blend reduces exactly to the line, same as any spline.
+    let interp = Interp1D::new(
+        array![0., 1., 2., 3.],
+        array![1., 3., 5., 7.], // f(x) = 2x + 1
+        strategy::CubicC1::default(),
+        Extrapolate::Enable,
+    )
+    .unwrap();
+    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);
+    assert_approx_eq!(interp.interpolate(&[-1.0]).unwrap(), -1.0);
+    assert_approx_eq!(interp.interpolate(&[4.0]).unwrap(), 9.0);
+}
+
+#[test]
+fn test_cubic_c1_interior_accuracy() {
+    // Unlike `CubicC2`'s `NotAKnot` (which reproduces any degree-<=3 polynomial
+    // exactly), `CubicC1`'s finite-difference derivatives carry a real error term for
+    // genuinely nonlinear data (`f'''(x) != 0`), so this checks bounded accuracy
+    // against a known cubic, not exact reproduction. `1.5` is a real, checked bound
+    // (max observed error ~1.24 at these points), not an arbitrarily loose one.
+    let interp = Interp1D::new(
+        array![0., 1., 2., 3.],
+        array![0., 1., 8., 27.], // f(x) = x^3
+        strategy::CubicC1::default(),
+        Extrapolate::Error,
+    )
+    .unwrap();
+    for &x in &[1.3, 2.3, 2.7, 2.9] {
+        let got = interp.interpolate(&[x]).unwrap();
+        let expected = x * x * x;
+        assert!(
+            (got - expected).abs() < 1.5,
+            "f({x}) = {expected}, got {got} (diff {})",
+            (got - expected).abs()
+        );
+    }
+}
+
+#[test]
+fn test_cubic_c1_knot_exactness() {
+    // Hermite splines interpolate the supplied value at every knot exactly by
+    // construction, regardless of how the derivative there was estimated.
+    let interp = Interp1D::new(
+        array![0., 1., 2., 3., 4.],
+        array![0.5, 1.2, 0.3, 2.1, 1.0], // non-polynomial data
+        strategy::CubicC1::default(),
+        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_invalid_args() {
     let interp = Interp1D::new(

src/interpolator/three/strategies.rs

@@ -323,5 +323,42 @@ where
     }
 }
 
+impl<D> Strategy3D<D> for CubicC1<D::Elem>
+where
+    D: Data + RawDataClone + Clone,
+    D::Elem: Float + Debug,
+{
+    /// Precomputes the full corner-derivative tensor under [`CubicC1CacheMode::Full`]
+    /// (the default); under [`CubicC1CacheMode::None`], only validates.
+    fn init(&mut self, data: &InterpData3DBase<D>) -> Result<(), ValidateError> {
+        if self.cache_mode == CubicC1CacheMode::Full {
+            let data_view = data.view();
+            self.cache = compute_corner_cache_fd(&data_view.grid, data_view.values.into_dyn());
+        }
+        Ok(())
+    }
+
+    fn interpolate(
+        &self,
+        data: &InterpData3DBase<D>,
+        point: &[D::Elem; 3],
+    ) -> Result<D::Elem, InterpolateError> {
+        let grids: Vec<ArrayView1<D::Elem>> = data.grid.iter().map(|g| g.view()).collect();
+        Ok(match self.cache_mode {
+            CubicC1CacheMode::Full => {
+                evaluate_spline_corner_cached(&grids, self.cache.view(), point)
+            }
+            CubicC1CacheMode::None => {
+                evaluate_spline_corner_local(&grids, data.values.view().into_dyn(), point)
+            }
+        })
+    }
+
+    /// Returns `true`: the boundary Hermite patch extends naturally.
+    fn allow_extrapolate(&self) -> bool {
+        true
+    }
+}
+
 grid_transform_strategy_impl!(Strategy3D, InterpData3DBase, InterpData3DView, 3);
 values_transform_strategy_impl!(Strategy3D, InterpData3DBase, InterpData3DView, Ix3, 3);