From 4755614ae976a3613f2e056c048e4bef9fa778f5 Mon Sep 17 00:00:00 2001 From: Kyle Carow Date: Sun, 23 Aug 2026 11:46:02 -0600 Subject: [PATCH 1/2] first draft --- CHANGELOG.md | 11 + src/interpolator/n/strategies.rs | 45 ++ src/interpolator/one/strategies.rs | 31 + src/interpolator/one/tests.rs | 35 ++ src/interpolator/three/strategies.rs | 37 ++ src/interpolator/three/tests.rs | 100 +++ src/interpolator/two/strategies.rs | 37 ++ src/interpolator/two/tests.rs | 102 +++ src/lib.rs | 1 + src/strategy/cubic/c1.rs | 272 ++++++++ src/strategy/cubic/c2.rs | 886 ++++++++++++++++++++++++++ src/strategy/cubic/mod.rs | 897 +-------------------------- src/strategy/cubic/utils.rs | 8 +- src/strategy/enums/n.rs | 1 + src/strategy/enums/one.rs | 1 + src/strategy/enums/three.rs | 1 + src/strategy/enums/two.rs | 1 + src/strategy/mod.rs | 2 +- tests/serde_strategies.rs | 13 + 19 files changed, 1593 insertions(+), 888 deletions(-) create mode 100644 src/strategy/cubic/c1.rs create mode 100644 src/strategy/cubic/c2.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 568990a..80d82b6 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/interpolator/n/strategies.rs b/src/interpolator/n/strategies.rs index 4db09f2..9fa3642 100644 --- a/src/interpolator/n/strategies.rs +++ b/src/interpolator/n/strategies.rs @@ -253,6 +253,51 @@ where } } +impl StrategyND for CubicC1 +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) -> 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, + point: &[D::Elem], + ) -> Result { + 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> = 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 StrategyND for GridTransform where D: Data + RawDataClone + Clone, diff --git a/src/interpolator/one/strategies.rs b/src/interpolator/one/strategies.rs index e9d1466..6e3b827 100644 --- a/src/interpolator/one/strategies.rs +++ b/src/interpolator/one/strategies.rs @@ -150,5 +150,36 @@ where } } +impl Strategy1D for CubicC1 +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) -> Result<(), ValidateError> { + self.cache = compute_fd_cache(data.grid[0].view(), data.values.view()); + Ok(()) + } + + fn interpolate( + &self, + data: &InterpData1DBase, + point: &[D::Elem; 1], + ) -> Result { + 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); diff --git a/src/interpolator/one/tests.rs b/src/interpolator/one/tests.rs index 3c1f580..b26d5e2 100644 --- a/src/interpolator/one/tests.rs +++ b/src/interpolator/one/tests.rs @@ -299,6 +299,41 @@ 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_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( diff --git a/src/interpolator/three/strategies.rs b/src/interpolator/three/strategies.rs index 0f21f87..a37d9f0 100644 --- a/src/interpolator/three/strategies.rs +++ b/src/interpolator/three/strategies.rs @@ -323,5 +323,42 @@ where } } +impl Strategy3D for CubicC1 +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) -> 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, + point: &[D::Elem; 3], + ) -> Result { + let grids: Vec> = 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); diff --git a/src/interpolator/three/tests.rs b/src/interpolator/three/tests.rs index 7386121..f1a65c7 100644 --- a/src/interpolator/three/tests.rs +++ b/src/interpolator/three/tests.rs @@ -415,6 +415,106 @@ fn test_cubic_c2_clamped_uses_given_derivative() { ); } +#[test] +fn test_cubic_c1_linear_exact() { + // Linear data: finite differences recover the exact constant slope on each axis, so + // the Hermite patch reduces exactly to the plane, under both cache modes. + fn f(x: f64, y: f64, z: f64) -> f64 { + 2. * x + 3. * y - z + 1. + } + let grid = [0., 1., 2., 3.]; + let values = Array3::from_shape_fn((4, 4, 4), |(i, j, k)| f(grid[i], grid[j], grid[k])); + for cache_mode in [CubicC1CacheMode::Full, CubicC1CacheMode::None] { + let interp = Interp3D::new( + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + values.clone(), + strategy::CubicC1::new().with_cache_mode(cache_mode), + Extrapolate::Enable, + ) + .unwrap(); + for &(x, y, z) in &[(0.5, 0.5, 0.5), (1.5, 2.5, 0.25), (-0.5, 3.5, 2.25)] { + assert_approx_eq!(interp.interpolate(&[x, y, z]).unwrap(), f(x, y, z)); + } + } +} + +#[test] +fn test_cubic_c1_full_none_agree() { + // `Full` and `None` build the same corner-derivative tensor at different scales + // (whole grid vs. a local window); they must agree exactly on any query. + fn f(x: f64, y: f64, z: f64) -> f64 { + x * x * y + y * y * z + z * z * x + } + let grid = [0., 1., 2., 3.]; + let values = Array3::from_shape_fn((4, 4, 4), |(i, j, k)| f(grid[i], grid[j], grid[k])); + let interp_full = Interp3D::new( + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + values.clone(), + strategy::CubicC1::default(), + Extrapolate::Error, + ) + .unwrap(); + let interp_none = Interp3D::new( + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + values, + strategy::CubicC1::new().with_cache_mode(CubicC1CacheMode::None), + Extrapolate::Error, + ) + .unwrap(); + for &(x, y, z) in &[(0.5, 0.5, 0.5), (1.5, 2.5, 0.25), (2.25, 0.75, 1.5)] { + assert_eq!( + interp_full.interpolate(&[x, y, z]).unwrap(), + interp_none.interpolate(&[x, y, z]).unwrap() + ); + } +} + +#[test] +fn test_cubic_c1_3d_matches_nd() { + // `Interp3D` and `InterpND` both build their corner-derivative tensor via + // `compute_corner_cache_fd`; this confirms the two wrappers agree on the same + // grid/values, under both cache modes. + fn f(x: f64, y: f64, z: f64) -> f64 { + x * x * y + y * y * z + z * z * x + } + let grid = [0., 1., 2., 3.]; + let values = Array3::from_shape_fn((4, 4, 4), |(i, j, k)| f(grid[i], grid[j], grid[k])); + for cache_mode in [CubicC1CacheMode::Full, CubicC1CacheMode::None] { + let interp3d = Interp3D::new( + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + values.clone(), + strategy::CubicC1::new().with_cache_mode(cache_mode), + Extrapolate::Error, + ) + .unwrap(); + let interp_nd = InterpND::new( + vec![ + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + ], + values.clone().into_dyn(), + strategy::CubicC1::new().with_cache_mode(cache_mode), + Extrapolate::Error, + ) + .unwrap(); + for &(x, y, z) in &[(0.5, 0.5, 0.5), (1.5, 2.5, 0.25), (2.25, 0.75, 1.5)] { + assert_approx_eq!( + interp3d.interpolate(&[x, y, z]).unwrap(), + interp_nd.interpolate(&[x, y, z]).unwrap() + ); + } + } +} + #[test] fn test_invalid_args() { let interp = Interp3D::new( diff --git a/src/interpolator/two/strategies.rs b/src/interpolator/two/strategies.rs index b959f96..544dcd8 100644 --- a/src/interpolator/two/strategies.rs +++ b/src/interpolator/two/strategies.rs @@ -217,5 +217,42 @@ where } } +impl Strategy2D for CubicC1 +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: &InterpData2DBase) -> 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: &InterpData2DBase, + point: &[D::Elem; 2], + ) -> Result { + let grids: Vec> = 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!(Strategy2D, InterpData2DBase, InterpData2DView, 2); values_transform_strategy_impl!(Strategy2D, InterpData2DBase, InterpData2DView, Ix2, 2); diff --git a/src/interpolator/two/tests.rs b/src/interpolator/two/tests.rs index a788aba..d34a4e0 100644 --- a/src/interpolator/two/tests.rs +++ b/src/interpolator/two/tests.rs @@ -443,6 +443,108 @@ fn test_cubic_c2_mixed_endpoints_scipy_oracle() { } } +#[test] +fn test_cubic_c1_linear_exact() { + // Linear data: finite differences recover the exact constant slope on each axis, so + // the Hermite patch reduces exactly to the plane, under both cache modes. + fn f(x: f64, y: f64) -> f64 { + 2. * x + 3. * y + 1. + } + let grid = array![0., 1., 2., 3.]; + let values = array![ + [f(0., 0.), f(0., 1.), f(0., 2.), f(0., 3.)], + [f(1., 0.), f(1., 1.), f(1., 2.), f(1., 3.)], + [f(2., 0.), f(2., 1.), f(2., 2.), f(2., 3.)], + [f(3., 0.), f(3., 1.), f(3., 2.), f(3., 3.)], + ]; + for cache_mode in [CubicC1CacheMode::Full, CubicC1CacheMode::None] { + let interp = Interp2D::new( + grid.clone(), + grid.clone(), + values.clone(), + strategy::CubicC1::new().with_cache_mode(cache_mode), + Extrapolate::Enable, + ) + .unwrap(); + for &(x, y) in &[(0.5, 0.5), (1.5, 2.5), (2.5, 1.5), (-0.5, 3.5)] { + assert_approx_eq!(interp.interpolate(&[x, y]).unwrap(), f(x, y)); + } + } +} + +#[test] +fn test_cubic_c1_full_none_agree() { + // `CubicC1CacheMode::Full` and `::None` build the same corner-derivative tensor at + // different scales (whole grid vs. a local window); they must agree exactly on any + // query, not just approximately. + let grid = array![0., 1., 2., 3.]; + let values = array![ + [0., 1., 4., 9.], + [1., 2., 5., 10.], + [4., 5., 8., 13.], + [9., 10., 13., 18.], + ]; // f(x, y) = x^2 + y + let interp_full = Interp2D::new( + grid.clone(), + grid.clone(), + values.clone(), + strategy::CubicC1::default(), + Extrapolate::Error, + ) + .unwrap(); + let interp_none = Interp2D::new( + grid.clone(), + grid, + values, + strategy::CubicC1::new().with_cache_mode(CubicC1CacheMode::None), + Extrapolate::Error, + ) + .unwrap(); + for &(x, y) in &[(0.5, 0.5), (1.5, 2.5), (2.5, 1.5), (0.25, 2.75)] { + assert_eq!( + interp_full.interpolate(&[x, y]).unwrap(), + interp_none.interpolate(&[x, y]).unwrap() + ); + } +} + +#[test] +fn test_cubic_c1_2d_matches_nd() { + // `Interp2D` and `InterpND` both build their corner-derivative tensor via + // `compute_corner_cache_fd` and evaluate via `evaluate_spline_corner_cached`; this + // confirms the two wrappers agree on the same grid/values, under both cache modes. + let grid = array![0., 1., 2., 3.]; + let values = array![ + [0., 1., 4., 9.], + [1., 2., 5., 10.], + [4., 5., 8., 13.], + [9., 10., 13., 18.], + ]; // f(x, y) = x^2 + y + for cache_mode in [CubicC1CacheMode::Full, CubicC1CacheMode::None] { + let interp2d = Interp2D::new( + grid.clone(), + grid.clone(), + values.clone(), + strategy::CubicC1::new().with_cache_mode(cache_mode), + Extrapolate::Error, + ) + .unwrap(); + let interp_nd = InterpND::new( + vec![grid.clone(), grid.clone()], + values.clone().into_dyn(), + strategy::CubicC1::new().with_cache_mode(cache_mode), + Extrapolate::Error, + ) + .unwrap(); + for &(x, y) in &[(0.5, 0.5), (1.5, 2.5), (2.5, 1.5), (0.25, 2.75)] { + assert_approx_eq!( + interp2d.interpolate(&[x, y]).unwrap(), + interp_nd.interpolate(&[x, y]).unwrap() + ); + } + } +} + #[test] fn test_invalid_args() { let interp = Interp2D::new( diff --git a/src/lib.rs b/src/lib.rs index 90e0101..a88ef1a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -20,6 +20,7 @@ extern crate alloc; /// directions) /// - [`strategy::Linear`] /// - [`strategy::LinearUniform`] +/// - [`strategy::CubicC1`] (C¹ local cubic Hermite spline) /// - [`strategy::CubicC2`] (C² cubic spline) /// - `serde`-compatible strategy enums: [`strategy::enums::Strategy1DEnum`]/etc. /// - The extrapolation setting enum: [`Extrapolate`](`interpolator::Extrapolate`) diff --git a/src/strategy/cubic/c1.rs b/src/strategy/cubic/c1.rs new file mode 100644 index 0000000..7c3b1fc --- /dev/null +++ b/src/strategy/cubic/c1.rs @@ -0,0 +1,272 @@ +//! [`CubicC1`]: C¹ local cubic Hermite spline (finite-difference derivative estimate). + +use super::utils::evaluate_hermite_1d; +use super::*; + +/// Local cubic Hermite spline interpolation (). +/// +/// Constructs a C¹ piecewise cubic Hermite polynomial through all data points, with +/// per-knot derivatives estimated locally (see [`CubicC1DerivativeMode`]) rather than +/// [`CubicC2`]'s global tridiagonal solve. Cheaper to build, at the cost of C¹ (not C²) +/// continuity and no configurable boundary condition. See [`CubicC1CacheMode`] for the +/// caching tradeoff at 2-D and above. +/// +/// Supports [`Extrapolate::Enable`](crate::interpolator::Extrapolate::Enable): +/// evaluation beyond the grid extends the boundary Hermite segment. +/// +/// # Example +/// ``` +/// use ndarray::prelude::*; +/// use ninterp::prelude::*; +/// +/// // f(x) = 2x + 1 (linear: reproduced exactly by any Hermite spline) +/// let interp: Interp1D = Interp1D::new( +/// array![0., 1., 2., 3.], +/// array![1., 3., 5., 7.], +/// strategy::CubicC1::default(), +/// Extrapolate::Enable, +/// ) +/// .unwrap(); +/// assert_eq!(interp.interpolate(&[1.5]).unwrap(), 4.0); +/// assert_eq!(interp.interpolate(&[4.0]).unwrap(), 9.0); // extrapolation +/// ``` +#[derive(Debug, Clone, PartialEq)] +#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] +#[cfg_attr( + feature = "serde", + serde(bound(serialize = "T: Serialize", deserialize = "T: Deserialize<'de>")) +)] +pub struct CubicC1 { + /// How per-knot derivatives are estimated. + pub derivative_mode: CubicC1DerivativeMode, + /// Caching tradeoff at 2-D and above; ignored at 1-D, where the cache is already + /// O(1) either way. + pub cache_mode: CubicC1CacheMode, + /// Precomputed derivative data under [`CubicC1CacheMode::Full`]; same shape as + /// [`CubicC2::cache`]. Empty under [`CubicC1CacheMode::None`] above 1-D. + #[cfg_attr(feature = "serde", serde(skip, default = "empty_cache"))] + pub(crate) cache: ArrayD, +} + +/// Derivative-estimation method for [`CubicC1`]. +#[non_exhaustive] +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] +pub enum CubicC1DerivativeMode { + /// Unclipped (non-monotonicity-preserving) finite differences: central in the + /// interior, one-sided at each boundary. + FiniteDifference, +} + +/// Caching strategy for [`CubicC1`] at 2-D and above. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] +pub enum CubicC1CacheMode { + /// Precompute the full corner-derivative tensor at `init()`, same mechanism as + /// [`CubicC2`]: every query afterward is an O(1) cache lookup. + Full, + /// Skip precomputation: `init()` only validates, and every query derives its + /// corner derivatives fresh from a local neighborhood. Cost depends only on grid + /// dimensionality, never grid resolution; agrees with `Full` exactly. + None, +} + +impl Default for CubicC1 { + /// [`CubicC1DerivativeMode::FiniteDifference`] + [`CubicC1CacheMode::Full`]. + fn default() -> Self { + Self { + derivative_mode: CubicC1DerivativeMode::FiniteDifference, + cache_mode: CubicC1CacheMode::Full, + cache: empty_cache(), + } + } +} + +impl CubicC1 { + /// Equivalent to [`default`](Self::default). + pub fn new() -> Self { + Self::default() + } + + /// Sets [`cache_mode`](Self::cache_mode). + pub fn with_cache_mode(mut self, cache_mode: CubicC1CacheMode) -> Self { + self.cache_mode = cache_mode; + self + } +} + +// `Full` and `None` cache modes below call the same functions, just at different +// scales: `Full` on the whole grid, `None` on a small local window. This works because +// `fd_derivatives` is a width-<=3 local stencil (unlike `CubicC2`'s `compute_m` +// global solve), so windowing never changes which branch (central vs. one-sided) fires +// for a given knot: the window is only ever clipped at a true grid boundary, exactly +// where the one-sided branch should fire anyway. `Full`/`None` agreeing exactly on any +// query is a consequence of this, not a separate property to maintain by hand. + +/// Finite-difference derivative at every knot of a 1-D lane: central in the interior, +/// one-sided at each end. `CubicC1`'s per-axis analog of `CubicC2`'s `compute_m`. +fn fd_derivatives(x: ArrayView1, y: ArrayView1) -> Vec { + let n = x.len(); + (0..n) + .map(|i| { + if i == 0 { + (y[1] - y[0]) / (x[1] - x[0]) + } else if i == n - 1 { + (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]) + } else { + (y[i + 1] - y[i - 1]) / (x[i + 1] - x[i - 1]) + } + }) + .collect() +} + +/// Caches [`fd_derivatives`]'s output for [`Strategy1D::init`](crate::strategy::traits::Strategy1D::init). +pub(crate) fn compute_fd_cache(x: ArrayView1, y: ArrayView1) -> ArrayD { + let d = fd_derivatives(x, y); + ArrayD::from_shape_vec(IxDyn(&[d.len()]), d) + .expect("fd_derivatives's output length matches its own shape") +} + +/// Evaluates [`Strategy1D`](crate::strategy::traits::Strategy1D)'s cached derivatives +/// (from [`compute_fd_cache`]) via [`evaluate_hermite_1d`] at `point`. +pub(crate) fn evaluate_hermite_1d_cached( + x: ArrayView1, + y: ArrayView1, + deriv_cache: ArrayViewD, + point: T, +) -> Result { + let m = deriv_cache.into_dimensionality::().map_err(|_| { + InterpolateError::Other( + "internal: non-1-D derivative cache, Strategy1D::cache invariant broken".into(), + ) + })?; + let i = locate_lower_index(x, &point); + let h = x[i + 1] - x[i]; + let t = (point - x[i]) / h; + Ok(evaluate_hermite_1d(y[i], m[i], y[i + 1], m[i + 1], h, t)) +} + +/// Same role as `CubicC2`'s function of the same name (differences one axis of a field +/// for the corner-derivative tensor), but no boundary condition to branch on, so there's +/// no homogenization pass either: every field (raw values or an already-differentiated +/// field) is just differenced the same way. +fn corner_cache_axis_pass( + grid: ArrayView1, + field: ArrayViewD, + axis: usize, +) -> ArrayD { + let axis = Axis(axis); + let mut out = ArrayD::::zeros(IxDyn(field.shape())); + for (y, mut out_lane) in field.lanes(axis).into_iter().zip(out.lanes_mut(axis)) { + let d = fd_derivatives(grid, y); + out_lane.assign(&ArrayView1::from(&d)); + } + out +} + +/// Finite-difference analog of `CubicC2`'s `compute_corner_cache`: builds the same +/// `2^N`-wide corner-derivative tensor shape, over whatever `grids`/`values` it's given. +/// Called on the whole grid for [`CubicC1CacheMode::Full`], or on a small local window +/// (see [`evaluate_spline_corner_local`]) for [`CubicC1CacheMode::None`]. +pub(crate) fn compute_corner_cache_fd( + grids: &[ArrayView1], + values: ArrayViewD, +) -> ArrayD { + let n_axes = grids.len(); + let n_bits = 1usize << n_axes; + let mask_axis = Axis(n_axes); + + let mut out_shape = values.shape().to_vec(); + out_shape.push(n_bits); + let mut cache = ArrayD::::zeros(IxDyn(&out_shape)); + cache.index_axis_mut(mask_axis, 0).assign(&values); + + for (axis, grid) in grids.iter().enumerate() { + let bit = 1usize << axis; + for mask in 0..bit { + let field = cache.index_axis(mask_axis, mask).to_owned(); + let deriv = corner_cache_axis_pass(*grid, field.view(), axis); + cache.index_axis_mut(mask_axis, mask | bit).assign(&deriv); + } + } + cache +} + +/// [`CubicC1CacheMode::None`]: extracts a local window (up to 4 points per axis) around +/// `point`'s bracket, builds a small corner-derivative tensor scoped to just that window +/// via [`compute_corner_cache_fd`], then evaluates it via [`evaluate_spline_corner_cached`] +/// unchanged. Cost is bounded by `grids.len()`, never by grid resolution. +pub(crate) fn evaluate_spline_corner_local( + grids: &[ArrayView1], + values: ArrayViewD, + point: &[T], +) -> T { + let n_axes = grids.len(); + let mut lowers = vec![0usize; n_axes]; + let mut uppers = vec![0usize; n_axes]; + for axis in 0..n_axes { + let l = locate_lower_index(grids[axis], &point[axis]); + lowers[axis] = l.saturating_sub(1); + uppers[axis] = (l + 2).min(grids[axis].len() - 1); + } + + let window_grids: Vec> = grids + .iter() + .enumerate() + .map(|(axis, g)| g.slice_axis(Axis(0), Slice::from(lowers[axis]..=uppers[axis]))) + .collect(); + + let window_values = values.slice_each_axis(|ax| { + let axis = ax.axis.index(); + Slice::from(lowers[axis]..=uppers[axis]) + }); + + let local_cache = compute_corner_cache_fd(&window_grids, window_values); + evaluate_spline_corner_cached(&window_grids, local_cache.view(), point) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn fd_derivatives_matches_hand_computed() { + // Non-uniform grid, non-polynomial data. + let x = [0.0, 0.4, 1.1, 1.8, 2.9, 4.0]; + let y = [0.5, 1.2, 0.3, 2.1, 1.0, 3.3]; + let xv = ArrayView1::from(&x); + let yv = ArrayView1::from(&y); + + let got = fd_derivatives(xv, yv); + let expected = [ + (y[1] - y[0]) / (x[1] - x[0]), // one-sided, lower boundary + (y[2] - y[0]) / (x[2] - x[0]), // central + (y[3] - y[1]) / (x[3] - x[1]), // central + (y[4] - y[2]) / (x[4] - x[2]), // central + (y[5] - y[3]) / (x[5] - x[3]), // central + (y[5] - y[4]) / (x[5] - x[4]), // one-sided, upper boundary + ]; + for (g, e) in got.iter().zip(expected.iter()) { + assert_approx_eq!(*g, *e, 1e-12); + } + } + + #[test] + #[cfg(feature = "serde")] + fn test_serde() { + assert_eq!( + serde_json::to_string(&CubicC1::::default()).unwrap(), + r#"{"derivative_mode":"FiniteDifference","cache_mode":"Full"}"# + ); + let none_mode = CubicC1::::new().with_cache_mode(CubicC1CacheMode::None); + assert_eq!( + serde_json::to_string(&none_mode).unwrap(), + r#"{"derivative_mode":"FiniteDifference","cache_mode":"None"}"# + ); + assert_eq!( + serde_json::from_str::>(&serde_json::to_string(&none_mode).unwrap()) + .unwrap(), + none_mode + ); + } +} diff --git a/src/strategy/cubic/c2.rs b/src/strategy/cubic/c2.rs new file mode 100644 index 0000000..743274d --- /dev/null +++ b/src/strategy/cubic/c2.rs @@ -0,0 +1,886 @@ +//! [`CubicC2`]: C² piecewise cubic spline (global tridiagonal solve). + +use super::*; + +/// Cubic spline interpolation (). +/// +/// Constructs a C² piecewise cubic polynomial through all data points. +/// The boundary condition is set by [`boundary_conditions`](CubicC2::boundary_conditions). +/// Coefficients are precomputed in [`Strategy1D::init`], called automatically +/// by [`Interp1D::new`](crate::interpolator::Interp1D::new) and +/// [`Interp1D::set_strategy`](crate::interpolator::Interp1D::set_strategy). +/// +/// Supports [`Extrapolate::Enable`](crate::interpolator::Extrapolate::Enable): +/// evaluation beyond the grid extends the boundary cubic polynomials. +/// +/// # Example +/// ``` +/// use ndarray::prelude::*; +/// use ninterp::prelude::*; +/// +/// // f(x) = 2x + 1 (linear: reproduced exactly by any spline) +/// let interp: Interp1D = Interp1D::new( +/// array![0., 1., 2., 3.], +/// array![1., 3., 5., 7.], +/// strategy::CubicC2::not_a_knot(), +/// Extrapolate::Enable, +/// ) +/// .unwrap(); +/// assert_eq!(interp.interpolate(&[1.5]).unwrap(), 4.0); +/// assert_eq!(interp.interpolate(&[4.0]).unwrap(), 9.0); // extrapolation +/// ``` +#[derive(Debug, Clone, PartialEq)] +#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] +#[cfg_attr( + feature = "serde", + serde(bound( + serialize = "T: Serialize + Zero", + deserialize = "T: Deserialize<'de> + Zero" + )) +)] +pub struct CubicC2 { + /// Boundary conditions, one per dimension or a single entry broadcast to all. + // Serializes under the key "CubicC2" rather than "boundary_conditions", making the + // strategy type explicit in the output, consistent with how Linear, Nearest, Step, + // etc. serialize to their type name. + #[cfg_attr(feature = "serde", serde(rename = "CubicC2"))] + pub boundary_conditions: Broadcastable>, + /// Precomputed derivative data, populated by `Strategy1D`/`2D`/`3D`/`ND::init`. Its + /// shape depends on which of those populated it: + /// + /// - `Strategy1D`: cached second derivatives (`M[i] = S''(x_i)`), one [`compute_m`] + /// result for the single 1-D pencil, i.e. shape `[n + 1]` for `n` intervals. + /// - `Strategy2D`/`Strategy3D`/`StrategyND`: the full corner-derivative tensor, + /// shaped like the value grid with one extra trailing axis of length `2^N` (`N` = + /// the grid's dimensionality): see [`compute_corner_cache`]. Every query is then + /// an O(1) lookup, no solving. + /// + /// Not included in the serialized form. After deserializing, call the + /// interpolator's `init_strategy` method (e.g. + /// [`Interp1D::init_strategy`](crate::interpolator::Interp1D::init_strategy)) to + /// recompute this before use. + #[cfg_attr(feature = "serde", serde(skip, default = "empty_cache"))] + pub(crate) cache: ArrayD, +} + +/// Boundary conditions for [`CubicC2`]. +/// +/// [`Endpoints::lower`](CubicC2BoundaryConditions::Endpoints)/`upper` are independent, so +/// mixing types (e.g. [`NotAKnot`](CubicC2Endpoint::NotAKnot) on one side, +/// [`FirstDerivative`](CubicC2Endpoint::FirstDerivative) on the other) is allowed. The +/// common symmetric cases have shorthand constructors: [`not_a_knot`](Self::not_a_knot), +/// [`first_derivative`](Self::first_derivative), [`second_derivative`](Self::second_derivative); +/// [`CubicC2::natural`] additionally shorthands `second_derivative(0, 0)`. +/// +/// Serializes as a bare string for [`NotAKnot`](CubicC2Endpoint::NotAKnot) (both +/// endpoints) and `Periodic`, matching their pre-`Endpoints` representation; a bare +/// `"Natural"` string is accepted and produced as a shorthand for symmetric zero +/// [`SecondDerivative`](CubicC2Endpoint::SecondDerivative) too. Anything else (an +/// explicit value, or asymmetric endpoints) uses the general `{"Endpoints": {"lower": +/// ..., "upper": ...}}` form. +#[derive(Debug, Clone, PartialEq)] +pub enum CubicC2BoundaryConditions { + /// Condition applied independently at each end of the axis. + Endpoints { + /// Condition at the lower endpoint. + lower: CubicC2Endpoint, + /// Condition at the upper endpoint. + upper: CubicC2Endpoint, + }, + /// First and second derivatives match at both endpoints. By convention + /// `values[n]` (the last point along this axis) should equal `values[0]`, since + /// that's what makes the axis periodic. This isn't enforced: `values[n]` is read and + /// used like any other data point (both to build the periodic system and to + /// evaluate the last interval), just never compared against `values[0]`. + /// + /// Deliberately not validated, even approximately: unlike grid coordinates (usually + /// synthetic, so float rounding is the only source of nonuniformity), `values` is + /// often real measured data, where the two ends of a period can legitimately differ + /// by far more than rounding error for reasons that have nothing to do with a + /// mistake (sensor noise, distinct samples at each end of the period, etc.). A + /// tolerance tight enough to catch a genuine mismatch would also reject that valid + /// data. + /// + /// If `values[n] != values[0]`, the spline within `[x[0], x[n]]` itself is still + /// perfectly smooth (both derivatives matched at the endpoints, same as if the axis + /// really were periodic) and passes through every supplied value exactly, `values[n]` + /// included. The mismatch only shows up if the axis is then treated as periodic + /// beyond its own bounds, e.g. via [`Extrapolate::Wrap`](crate::interpolator::Extrapolate::Wrap): + /// crossing the seam is a jump in *value* of exactly `values[n] - values[0]`, with + /// slope and curvature matching continuously on both sides of it. + Periodic, +} + +/// A single endpoint's condition, for [`CubicC2BoundaryConditions::Endpoints`]. +#[derive(Debug, Clone, PartialEq)] +#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] +pub enum CubicC2Endpoint { + /// C³ continuity at the second (from this end) knot; no extra input required. + /// Generally gives better accuracy than natural (a zero + /// [`SecondDerivative`](Self::SecondDerivative)) for smooth functions. Requires at + /// least 3 grid points along this axis (4 if both endpoints are `NotAKnot`). + NotAKnot, + /// Specified first derivative at this endpoint ("clamped"). + FirstDerivative(T), + /// Specified second derivative at this endpoint. Zero, the classic "natural" + /// condition, has a shorthand: [`CubicC2::natural`]. + SecondDerivative(T), +} + +impl CubicC2BoundaryConditions { + /// Not-a-knot at both ends. Requires at least 4 data points per dimension. + pub fn not_a_knot() -> Self { + Self::Endpoints { + lower: CubicC2Endpoint::NotAKnot, + upper: CubicC2Endpoint::NotAKnot, + } + } + + /// Specified first derivative at both endpoints. + pub fn first_derivative(lower: T, upper: T) -> Self { + Self::Endpoints { + lower: CubicC2Endpoint::FirstDerivative(lower), + upper: CubicC2Endpoint::FirstDerivative(upper), + } + } + + /// Specified second derivative at both endpoints. + pub fn second_derivative(lower: T, upper: T) -> Self { + Self::Endpoints { + lower: CubicC2Endpoint::SecondDerivative(lower), + upper: CubicC2Endpoint::SecondDerivative(upper), + } + } +} + +impl From> for CubicC2 { + /// Broadcasts `bc` to all dimensions. + /// + /// Use [`CubicC2::not_a_knot`], [`natural`](CubicC2::natural), + /// [`clamped`](CubicC2::clamped), or [`periodic`](CubicC2::periodic) instead when the + /// desired condition is known at the call site; this is for a `CubicC2BoundaryConditions` + /// value obtained generically (e.g. from runtime config), without matching on it first. + fn from(bc: CubicC2BoundaryConditions) -> Self { + Self { + boundary_conditions: Broadcastable::Broadcast(bc), + cache: empty_cache(), + } + } +} + +impl CubicC2 { + /// Create a cubic spline with a distinct boundary condition per grid dimension. + /// + /// Use [`not_a_knot`](Self::not_a_knot), [`natural`](Self::natural), + /// [`clamped`](Self::clamped), or [`periodic`](Self::periodic) instead when every + /// dimension shares the same condition. + pub fn new(boundary_conditions: Vec>) -> Self { + Self { + boundary_conditions: Broadcastable::Each(boundary_conditions), + cache: empty_cache(), + } + } + + /// Create a cubic spline with not-a-knot boundary conditions. + /// Requires at least 4 data points per dimension. + pub fn not_a_knot() -> Self { + Self { + boundary_conditions: Broadcastable::Broadcast(CubicC2BoundaryConditions::not_a_knot()), + cache: empty_cache(), + } + } + + /// Create a cubic spline with natural (zero second derivative at endpoints) BCs. + pub fn natural() -> Self + where + T: Zero, + { + Self { + boundary_conditions: Broadcastable::Broadcast( + CubicC2BoundaryConditions::second_derivative(T::zero(), T::zero()), + ), + cache: empty_cache(), + } + } + + /// Create a cubic spline with specified first derivatives at both endpoints. + pub fn clamped(lower: T, upper: T) -> Self { + Self { + boundary_conditions: Broadcastable::Broadcast( + CubicC2BoundaryConditions::first_derivative(lower, upper), + ), + cache: empty_cache(), + } + } + + /// Create a cubic spline with periodic boundary conditions. By convention + /// `values[n]` (the last point along each periodic axis) should equal + /// `values[0]`, though this isn't enforced. + pub fn periodic() -> Self { + Self { + boundary_conditions: Broadcastable::Broadcast(CubicC2BoundaryConditions::Periodic), + cache: empty_cache(), + } + } +} + +/// Thomas algorithm (tridiagonal matrix algorithm). Solves `A * x = rhs`. +/// `sub.len() == sup.len() == diag.len() - 1`. +pub(crate) fn thomas(sub: &[T], diag: &[T], sup: &[T], rhs: &[T]) -> Vec { + let n = diag.len(); + let mut cp = vec![T::zero(); n]; + let mut dp = vec![T::zero(); n]; + cp[0] = if n > 1 { sup[0] / diag[0] } else { T::zero() }; + dp[0] = rhs[0] / diag[0]; + for k in 1..n { + let w = diag[k] - sub[k - 1] * cp[k - 1]; + cp[k] = if k < n - 1 { sup[k] / w } else { T::zero() }; + dp[k] = (rhs[k] - sub[k - 1] * dp[k - 1]) / w; + } + let mut x = vec![T::zero(); n]; + x[n - 1] = dp[n - 1]; + for k in (0..n - 1).rev() { + x[k] = dp[k] - cp[k] * x[k + 1]; + } + x +} + +/// Sherman-Morrison cyclic tridiagonal solver. +/// Corner elements `corner` appear at `(0, n-1)` and `(n-1, 0)`. +/// `sub.len() == sup.len() == n - 1`. +pub(crate) fn cyclic_thomas( + sub: &[T], + diag: &[T], + sup: &[T], + rhs: &[T], + corner: T, +) -> Vec { + let n = diag.len(); + if n == 1 { + return vec![rhs[0] / (diag[0] + corner + corner)]; + } + let gamma = -diag[0]; + let c_over_g = corner / gamma; + let mut diag_mod = diag.to_vec(); + diag_mod[0] = diag_mod[0] - gamma; + diag_mod[n - 1] = diag_mod[n - 1] - corner * corner / gamma; + let y = thomas(sub, &diag_mod, sup, rhs); + let mut u_vec = vec![T::zero(); n]; + u_vec[0] = gamma; + u_vec[n - 1] = corner; + let z = thomas(sub, &diag_mod, sup, &u_vec); + let vt_y = y[0] + c_over_g * y[n - 1]; + let vt_z = z[0] + c_over_g * z[n - 1]; + let factor = vt_y / (T::one() + vt_z); + y.into_iter() + .zip(z.iter()) + .map(|(yi, zi)| yi - factor * *zi) + .collect() +} + +/// Computes the second-derivative vector `M[0..=n]` for the given cubic spline BC. +/// +/// Used by [`compute_m_cache`] (called from `Strategy1D::init`, stored in `CubicC2::cache`) +/// and [`corner_cache_axis_pass`] (called from [`compute_corner_cache`]). +pub(crate) fn compute_m( + x: ArrayView1, + y: ArrayView1, + bc: &CubicC2BoundaryConditions, +) -> Vec { + let n = x.len() - 1; + let two = T::one() + T::one(); + let six = two + two + two; + let h: Vec = (0..n).map(|i| x[i + 1] - x[i]).collect(); + let slopes: Vec = (0..n).map(|i| (y[i + 1] - y[i]) / h[i]).collect(); + let u: Vec = (0..n.saturating_sub(1)) + .map(|k| six * (slopes[k + 1] - slopes[k])) + .collect(); + + match bc { + CubicC2BoundaryConditions::Endpoints { lower, upper } => { + // Each endpoint contributes either a direct boundary row (`Some`, referencing + // only its own knot and its immediate interior neighbor) or, for `NotAKnot`, + // `None`: that side's M is eliminated via the third-derivative-continuity + // relation at knot 1 (or n-1), folded into the interior equation there instead + // of appearing as its own unknown. `lower_row`/`upper_row` give + // `(diag, off_diagonal, rhs)` for the direct case. + let lower_row = match lower { + CubicC2Endpoint::NotAKnot => None, + CubicC2Endpoint::SecondDerivative(v) => Some((T::one(), T::zero(), *v)), + CubicC2Endpoint::FirstDerivative(v) => { + Some((two * h[0], h[0], six * (slopes[0] - *v))) + } + }; + let upper_row = match upper { + CubicC2Endpoint::NotAKnot => None, + CubicC2Endpoint::SecondDerivative(v) => Some((T::one(), T::zero(), *v)), + CubicC2Endpoint::FirstDerivative(v) => { + Some((two * h[n - 1], h[n - 1], six * (*v - slopes[n - 1]))) + } + }; + + // `validate_bc_min_points` guarantees `count >= 2` (grid_len >= 3 with one + // `NotAKnot` side, >= 4 with both), so the `j == 0` and `j == count - 1` + // branches below are always distinct rows. + let start = if lower_row.is_some() { 0 } else { 1 }; + let end = if upper_row.is_some() { n } else { n - 1 }; + let count = end - start + 1; + + let mut sub = Vec::with_capacity(count - 1); + let mut diag = Vec::with_capacity(count); + let mut sup = Vec::with_capacity(count - 1); + let mut rhs = Vec::with_capacity(count); + for j in 0..count { + let i = start + j; + if j == 0 { + if let Some((diag_v, sup_v, rhs_v)) = lower_row { + diag.push(diag_v); + sup.push(sup_v); + rhs.push(rhs_v); + } else { + // Folds the eliminated M[0] into knot 1's interior equation. + diag.push((h[0] + h[1]) * (h[0] + two * h[1])); + sup.push(h[1] * h[1] - h[0] * h[0]); + rhs.push(h[1] * u[0]); + } + } else if j == count - 1 { + if let Some((diag_v, sub_v, rhs_v)) = upper_row { + sub.push(sub_v); + diag.push(diag_v); + rhs.push(rhs_v); + } else { + // Mirrors the lower fold, at knot n-1 eliminating M[n]. + sub.push(h[n - 2] * h[n - 2] - h[n - 1] * h[n - 1]); + diag.push((h[n - 2] + h[n - 1]) * (two * h[n - 2] + h[n - 1])); + rhs.push(h[n - 2] * u[n - 2]); + } + } else { + sub.push(h[i - 1]); + diag.push(two * (h[i - 1] + h[i])); + sup.push(h[i]); + rhs.push(u[i - 1]); + } + } + let inner = thomas(&sub, &diag, &sup, &rhs); + + let mut m = Vec::with_capacity(n + 1); + if lower_row.is_none() { + m.push(((h[0] + h[1]) * inner[0] - h[0] * inner[1]) / h[1]); + } + m.extend_from_slice(&inner); + if upper_row.is_none() { + let last = inner.len() - 1; + m.push( + ((h[n - 2] + h[n - 1]) * inner[last] - h[n - 1] * inner[last - 1]) / h[n - 2], + ); + } + m + } + CubicC2BoundaryConditions::Periodic => { + // `y[n]` is read below, via `slopes[n - 1]` in the `rhs[0]` line, and + // `evaluate_spline_from_m` also reads `y[n]` directly when evaluating the last + // interval. By convention `y[n]` should equal `y[0]`, but nothing here + // compares them; `y[n]` is just used as ordinary data either way. + if n < 2 { + vec![T::zero(); n + 1] + } else { + let sub_sup = h[..n - 1].to_vec(); + let mut diag = vec![two * (h[n - 1] + h[0])]; + for k in 1..n { + diag.push(two * (h[k - 1] + h[k])); + } + let mut rhs = vec![six * (slopes[0] - slopes[n - 1])]; + rhs.extend_from_slice(&u); + let corner = h[n - 1]; + let mut m_vals = cyclic_thomas(&sub_sup, &diag, &sub_sup, &rhs, corner); + let m0 = m_vals[0]; + m_vals.push(m0); + m_vals + } + } + } +} + +/// Evaluates the M-form cubic spline at `point` using precomputed second derivatives `m`. +pub(crate) fn evaluate_spline_from_m( + x: ArrayView1, + y: ArrayView1, + m: ArrayView1, + point: T, +) -> T { + let two = T::one() + T::one(); + let six = two + two + two; + let i = locate_lower_index(x, &point); + let h = x[i + 1] - x[i]; + let dx = point - x[i]; + let dx_r = h - dx; + let six_h = six * h; + let h2_over_six = h * h / six; + m[i] * dx_r * dx_r * dx_r / six_h + + m[i + 1] * dx * dx * dx / six_h + + (y[i] - m[i] * h2_over_six) * dx_r / h + + (y[i + 1] - m[i + 1] * h2_over_six) * dx / h +} + +/// Checks `grid_len` against boundary condition `bc`'s minimum point requirement (e.g. +/// [`CubicC2Endpoint::NotAKnot`] needs at least 3 grid points on its own side, 4 if both +/// endpoints are `NotAKnot`), ahead of the real work in [`compute_m`]. +/// +/// Pure, no mutation; used by each dimensionality's `Strategy*D::validate`. +pub(crate) fn validate_bc_min_points( + bc: &CubicC2BoundaryConditions, + grid_len: usize, + dim: usize, +) -> Result<(), ValidateError> { + let CubicC2BoundaryConditions::Endpoints { lower, upper } = bc else { + return Ok(()); + }; + let min = match ( + matches!(lower, CubicC2Endpoint::NotAKnot), + matches!(upper, CubicC2Endpoint::NotAKnot), + ) { + (true, true) => 4, + (true, false) | (false, true) => 3, + (false, false) => return Ok(()), + }; + if grid_len < min { + return Err(ValidateError::Other(format!( + "CubicC2: dim {dim} has {grid_len} grid points; NotAKnot requires at least {min}" + ))); + } + Ok(()) +} + +/// Computes and caches `M[0..=n]` (`CubicC2::cache`) for [`Strategy1D::init`], so +/// [`evaluate_spline_1d_cached`] can look them up in O(1) instead of re-solving on every +/// `interpolate` call. +pub(crate) fn compute_m_cache( + x: ArrayView1, + y: ArrayView1, + bc: &CubicC2BoundaryConditions, +) -> ArrayD { + let m = compute_m(x, y, bc); + ArrayD::from_shape_vec(IxDyn(&[m.len()]), m) + .expect("compute_m's output length matches its own shape") +} + +/// Evaluates [`Strategy1D`]'s cached spline (`m_cache`, from [`compute_m_cache`]) at `point`. +pub(crate) fn evaluate_spline_1d_cached( + x: ArrayView1, + y: ArrayView1, + m_cache: ArrayViewD, + point: T, +) -> Result { + let m = m_cache.into_dimensionality::().map_err(|_| { + InterpolateError::Other( + "internal: non-1-D m_cache, Strategy1D::cache invariant broken".into(), + ) + })?; + Ok(evaluate_spline_from_m(x, y, m, point)) +} + +/// Closed-form first derivative `S'(x_i)` at every knot, from an already-solved moment +/// vector `m` (no extra solve). The companion to [`compute_m`], used to build +/// [`compute_corner_cache`]'s derivative fields. +pub(crate) fn knot_derivatives_from_m( + x: ArrayView1, + y: ArrayView1, + m: ArrayView1, +) -> Vec { + let n = x.len() - 1; + let two = T::one() + T::one(); + let six = two + two + two; + let mut d: Vec = (0..n) + .map(|i| { + let h = x[i + 1] - x[i]; + (y[i + 1] - y[i]) / h - h * (two * m[i] + m[i + 1]) / six + }) + .collect(); + let h_last = x[n] - x[n - 1]; + d.push((y[n] - y[n - 1]) / h_last + h_last * (m[n - 1] + two * m[n]) / six); + d +} + +/// Splines every 1-D lane of `field` along `axis` and replaces it with its knot +/// derivatives (via [`compute_m`] + [`knot_derivatives_from_m`]), returning a new array +/// the same shape as `field`, for [`compute_corner_cache`]. +fn corner_cache_axis_pass( + grid: ArrayView1, + field: ArrayViewD, + axis: usize, + bc: &CubicC2BoundaryConditions, +) -> ArrayD { + let axis = Axis(axis); + let mut out = ArrayD::::zeros(IxDyn(field.shape())); + for (y, mut out_lane) in field.lanes(axis).into_iter().zip(out.lanes_mut(axis)) { + let m = compute_m(grid, y, bc); + let d = knot_derivatives_from_m(grid, y, ArrayView1::from(&m)); + out_lane.assign(&ArrayView1::from(&d)); + } + out +} + +/// Precomputes the full corner-derivative tensor for [`CubicC2`]'s `Strategy2D`/ +/// `Strategy3D` full-cache upgrade: for every grid point, all `2^N` partial-derivative +/// combinations (value, first partials, and mixed partials) needed to evaluate a Hermite +/// patch in O(1) via [`evaluate_spline_corner_cached`]. +/// +/// Returns an array shaped like `values`, with one extra trailing axis of length `2^N` +/// (`N = grids.len()`). Index `k` in that axis is a bitmask over the `N` grid axes: bit +/// `i` set means "this entry is `d/dx_i` of the value", so `k == 0` is the raw value and +/// `k == 2^N - 1` is the full mixed partial. +/// +/// Built by processing axes `0..N` in order: for every mask not yet including axis `i`'s +/// bit, splines the corresponding field along axis `i` ([`corner_cache_axis_pass`]) to +/// populate the entry with that bit added. Boundary condition per axis is +/// `bcs[axis]` when splining the raw values (`mask == 0`, this axis's own +/// first-derivative pass). For every other `mask` (a cross-derivative pass over an +/// already-differentiated field, not the original data), each endpoint's +/// [`FirstDerivative`](CubicC2Endpoint::FirstDerivative)/[`SecondDerivative`](CubicC2Endpoint::SecondDerivative) +/// falls back to the same-order condition with a zero value: a `FirstDerivative`/ +/// `SecondDerivative` endpoint fixes the same scalar derivative at every point along that +/// boundary, so differentiating that (constant-along-the-boundary) field with respect to +/// any other axis must itself be zero at that endpoint, at the same order. `NotAKnot` +/// endpoints pass through unchanged (they carry no value to homogenize). Verified +/// empirically: with non-separable data, only these homogeneous fallbacks make this +/// cached path agree (to float precision) with `StrategyND`'s independent recursive +/// solve; leaving a nonzero endpoint value unchanged in a cross-derivative pass does not. +pub(crate) fn compute_corner_cache( + grids: &[ArrayView1], + values: ArrayViewD, + bcs: &Broadcastable>, +) -> ArrayD { + let n_axes = grids.len(); + let n_bits = 1usize << n_axes; + let mask_axis = Axis(n_axes); + + let mut out_shape = values.shape().to_vec(); + out_shape.push(n_bits); + let mut cache = ArrayD::::zeros(IxDyn(&out_shape)); + cache.index_axis_mut(mask_axis, 0).assign(&values); + + for (axis, grid) in grids.iter().enumerate() { + let bit = 1usize << axis; + let bc_axis = &bcs[axis]; + let cross_endpoint = |e: &CubicC2Endpoint| match e { + CubicC2Endpoint::FirstDerivative(_) => CubicC2Endpoint::FirstDerivative(T::zero()), + CubicC2Endpoint::SecondDerivative(_) => CubicC2Endpoint::SecondDerivative(T::zero()), + other @ CubicC2Endpoint::NotAKnot => other.clone(), + }; + let cross_bc = match bc_axis { + CubicC2BoundaryConditions::Endpoints { lower, upper } => { + CubicC2BoundaryConditions::Endpoints { + lower: cross_endpoint(lower), + upper: cross_endpoint(upper), + } + } + CubicC2BoundaryConditions::Periodic => CubicC2BoundaryConditions::Periodic, + }; + for mask in 0..bit { + let bc = if mask == 0 { bc_axis } else { &cross_bc }; + let field = cache.index_axis(mask_axis, mask).to_owned(); + let deriv = corner_cache_axis_pass(*grid, field.view(), axis, bc); + cache.index_axis_mut(mask_axis, mask | bit).assign(&deriv); + } + } + cache +} + +#[cfg(feature = "serde")] +mod cubic_serde { + use super::*; + use serde::Deserializer; + + /// Bare-string wire form for the two `CubicC2BoundaryConditions::Endpoints` cases + /// that carry no distinguishing per-endpoint data (symmetric `NotAKnot`, symmetric + /// zero `SecondDerivative`, i.e. "natural"), plus `Periodic`. Anything else uses the + /// general `Endpoints` shape instead. + #[derive(Deserialize, Serialize)] + enum BcName { + NotAKnot, + Natural, + Periodic, + } + + /// General `Endpoints` shape, owned form for [`Deserialize`]. + #[derive(Deserialize)] + struct EndpointsOwned { + lower: CubicC2Endpoint, + upper: CubicC2Endpoint, + } + + /// General `Endpoints` shape, borrowed form for [`Serialize`]: avoids requiring `T: + /// Clone` just to serialize. + #[derive(Serialize)] + struct EndpointsRef<'a, T> { + lower: &'a CubicC2Endpoint, + upper: &'a CubicC2Endpoint, + } + + #[derive(Deserialize)] + #[serde(untagged)] + enum BcWireDe { + Named(BcName), + Endpoints { + #[serde(rename = "Endpoints")] + endpoints: EndpointsOwned, + }, + } + + #[derive(Serialize)] + #[serde(untagged)] + enum BcWireSer<'a, T> { + Named(BcName), + Endpoints { + #[serde(rename = "Endpoints")] + endpoints: EndpointsRef<'a, T>, + }, + } + + impl Serialize for CubicC2BoundaryConditions { + fn serialize(&self, serializer: S) -> Result { + match self { + Self::Periodic => BcWireSer::::Named(BcName::Periodic).serialize(serializer), + Self::Endpoints { + lower: CubicC2Endpoint::NotAKnot, + upper: CubicC2Endpoint::NotAKnot, + } => BcWireSer::::Named(BcName::NotAKnot).serialize(serializer), + Self::Endpoints { + lower: CubicC2Endpoint::SecondDerivative(lower), + upper: CubicC2Endpoint::SecondDerivative(upper), + } if lower.is_zero() && upper.is_zero() => { + BcWireSer::::Named(BcName::Natural).serialize(serializer) + } + Self::Endpoints { lower, upper } => BcWireSer::Endpoints { + endpoints: EndpointsRef { lower, upper }, + } + .serialize(serializer), + } + } + } + + impl<'de, T: Deserialize<'de> + Zero> Deserialize<'de> for CubicC2BoundaryConditions { + fn deserialize>(deserializer: D) -> Result { + Ok(match BcWireDe::::deserialize(deserializer)? { + BcWireDe::Named(BcName::NotAKnot) => CubicC2BoundaryConditions::not_a_knot(), + BcWireDe::Named(BcName::Natural) => Self::Endpoints { + lower: CubicC2Endpoint::SecondDerivative(T::zero()), + upper: CubicC2Endpoint::SecondDerivative(T::zero()), + }, + BcWireDe::Named(BcName::Periodic) => Self::Periodic, + BcWireDe::Endpoints { endpoints } => Self::Endpoints { + lower: endpoints.lower, + upper: endpoints.upper, + }, + }) + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn from_bc_matches_named_constructors() { + assert_eq!( + CubicC2::from(CubicC2BoundaryConditions::::not_a_knot()), + CubicC2::not_a_knot() + ); + assert_eq!( + CubicC2::from(CubicC2BoundaryConditions::::second_derivative( + 0.0, 0.0 + )), + CubicC2::natural() + ); + assert_eq!( + CubicC2::from(CubicC2BoundaryConditions::first_derivative(1.0, 2.0)), + CubicC2::clamped(1.0, 2.0) + ); + assert_eq!( + CubicC2::from(CubicC2BoundaryConditions::::Periodic), + CubicC2::periodic() + ); + } + + /// Builds the `(n+1) x (n+1)` moment system directly from `bc`'s mathematical + /// definition, independent of `compute_m`'s reduced/eliminated form, to cross-check + /// it on asymmetric endpoint combinations that this crate's other tests (all + /// symmetric until now) wouldn't have caught a mismatched row for. + fn dense_m_system( + x: &[f64], + y: &[f64], + bc: &CubicC2BoundaryConditions, + ) -> (Vec>, Vec) { + let n = x.len() - 1; + let h: Vec = (0..n).map(|i| x[i + 1] - x[i]).collect(); + let slopes: Vec = (0..n).map(|i| (y[i + 1] - y[i]) / h[i]).collect(); + let mut a = vec![vec![0.0; n + 1]; n + 1]; + let mut b = vec![0.0; n + 1]; + + for i in 1..n { + a[i][i - 1] = h[i - 1]; + a[i][i] = 2.0 * (h[i - 1] + h[i]); + a[i][i + 1] = h[i]; + b[i] = 6.0 * (slopes[i] - slopes[i - 1]); + } + + let CubicC2BoundaryConditions::Endpoints { lower, upper } = bc else { + panic!("dense_m_system only supports Endpoints for this test"); + }; + match lower { + CubicC2Endpoint::NotAKnot => { + a[0][0] = h[1]; + a[0][1] = -(h[0] + h[1]); + a[0][2] = h[0]; + } + CubicC2Endpoint::SecondDerivative(v) => { + a[0][0] = 1.0; + b[0] = *v; + } + CubicC2Endpoint::FirstDerivative(v) => { + a[0][0] = 2.0 * h[0]; + a[0][1] = h[0]; + b[0] = 6.0 * (slopes[0] - v); + } + } + match upper { + CubicC2Endpoint::NotAKnot => { + a[n][n] = h[n - 2]; + a[n][n - 1] = -(h[n - 2] + h[n - 1]); + a[n][n - 2] = h[n - 1]; + } + CubicC2Endpoint::SecondDerivative(v) => { + a[n][n] = 1.0; + b[n] = *v; + } + CubicC2Endpoint::FirstDerivative(v) => { + a[n][n] = 2.0 * h[n - 1]; + a[n][n - 1] = h[n - 1]; + b[n] = 6.0 * (v - slopes[n - 1]); + } + } + (a, b) + } + + /// Plain Gaussian elimination with partial pivoting, independent of `thomas`/ + /// `cyclic_thomas`, for cross-checking their output. + fn dense_solve(mut a: Vec>, mut b: Vec) -> Vec { + let n = b.len(); + for col in 0..n { + let pivot = (col..n) + .max_by(|&r1, &r2| a[r1][col].abs().partial_cmp(&a[r2][col].abs()).unwrap()) + .unwrap(); + a.swap(col, pivot); + b.swap(col, pivot); + for row in (col + 1)..n { + let factor = a[row][col] / a[col][col]; + let pivot_row = a[col].clone(); + for (c, pivot_val) in pivot_row.iter().enumerate().skip(col) { + a[row][c] -= factor * pivot_val; + } + b[row] -= factor * b[col]; + } + } + let mut x = vec![0.0; n]; + for row in (0..n).rev() { + let sum: f64 = (row + 1..n).map(|c| a[row][c] * x[c]).sum(); + x[row] = (b[row] - sum) / a[row][row]; + } + x + } + + #[test] + fn compute_m_matches_dense_solve_for_mixed_endpoints() { + // Non-uniform grid, non-polynomial data: exercises the linear algebra for real, + // rather than something any BC would reproduce exactly regardless of correctness. + let x = [0.0, 0.4, 1.1, 1.8, 2.9, 4.0]; + let y = [0.5, 1.2, 0.3, 2.1, 1.0, 3.3]; + let xv = ArrayView1::from(&x); + let yv = ArrayView1::from(&y); + + let endpoints = [ + CubicC2Endpoint::NotAKnot, + CubicC2Endpoint::SecondDerivative(0.0), + CubicC2Endpoint::SecondDerivative(2.5), + CubicC2Endpoint::FirstDerivative(0.0), + CubicC2Endpoint::FirstDerivative(-1.5), + ]; + for lower in &endpoints { + for upper in &endpoints { + let bc = CubicC2BoundaryConditions::Endpoints { + lower: lower.clone(), + upper: upper.clone(), + }; + let got = compute_m(xv, yv, &bc); + let (a, b) = dense_m_system(&x, &y, &bc); + let expected = dense_solve(a, b); + for (g, e) in got.iter().zip(expected.iter()) { + assert_approx_eq!(*g, *e, 1e-9); + } + } + } + } + + #[test] + #[cfg(feature = "serde")] + fn test_serde() { + // The whole-axis cases collapse to a bare string, keyed by strategy name like + // every other built-in strategy. + assert_eq!( + serde_json::to_string(&CubicC2::::not_a_knot()).unwrap(), + r#"{"CubicC2":"NotAKnot"}"# + ); + assert_eq!( + serde_json::to_string(&CubicC2::::natural()).unwrap(), + r#"{"CubicC2":"Natural"}"# + ); + assert_eq!( + serde_json::to_string(&CubicC2::::periodic()).unwrap(), + r#"{"CubicC2":"Periodic"}"# + ); + // Anything carrying an explicit value, or mixing endpoint types, uses the + // general `Endpoints` form instead. + assert_eq!( + serde_json::to_string(&CubicC2::::clamped(1.0, 2.0)).unwrap(), + r#"{"CubicC2":{"Endpoints":{"lower":{"FirstDerivative":1.0},"upper":{"FirstDerivative":2.0}}}}"# + ); + let mixed = CubicC2::::from(CubicC2BoundaryConditions::Endpoints { + lower: CubicC2Endpoint::NotAKnot, + upper: CubicC2Endpoint::FirstDerivative(3.0), + }); + assert_eq!( + serde_json::to_string(&mixed).unwrap(), + r#"{"CubicC2":{"Endpoints":{"lower":"NotAKnot","upper":{"FirstDerivative":3.0}}}}"# + ); + assert_eq!( + serde_json::from_str::>(&serde_json::to_string(&mixed).unwrap()).unwrap(), + mixed + ); + + // Bare "NotAKnot"/"Natural" round-trip; the fully expanded `Endpoints` form is + // also accepted on read for both, even though only the bare form is ever written. + assert_eq!( + serde_json::from_str::>(r#"{"CubicC2":"NotAKnot"}"#).unwrap(), + CubicC2::::not_a_knot() + ); + assert_eq!( + serde_json::from_str::>(r#"{"CubicC2":"Natural"}"#).unwrap(), + CubicC2::::natural() + ); + assert_eq!( + serde_json::from_str::>( + r#"{"CubicC2":{"Endpoints":{"lower":"NotAKnot","upper":"NotAKnot"}}}"# + ) + .unwrap(), + CubicC2::::not_a_knot() + ); + assert_eq!( + serde_json::from_str::>( + r#"{"CubicC2":{"Endpoints":{"lower":{"SecondDerivative":0.0},"upper":{"SecondDerivative":0.0}}}}"# + ) + .unwrap(), + CubicC2::::natural() + ); + } +} diff --git a/src/strategy/cubic/mod.rs b/src/strategy/cubic/mod.rs index f6a2192..9ae7d19 100644 --- a/src/strategy/cubic/mod.rs +++ b/src/strategy/cubic/mod.rs @@ -1,889 +1,20 @@ -//! Cubic interpolation algorithms shared across all dimensionalities, for [`CubicC2`]. +//! Cubic interpolation algorithms shared across all dimensionalities, +//! for [`CubicC1`] and [`CubicC2`]. use super::*; +mod c1; +mod c2; mod utils; -pub(crate) use utils::evaluate_spline_corner_cached; - -/// Cubic spline interpolation (). -/// -/// Constructs a C² piecewise cubic polynomial through all data points. -/// The boundary condition is set by [`boundary_conditions`](CubicC2::boundary_conditions). -/// Coefficients are precomputed in [`Strategy1D::init`], called automatically -/// by [`Interp1D::new`](crate::interpolator::Interp1D::new) and -/// [`Interp1D::set_strategy`](crate::interpolator::Interp1D::set_strategy). -/// -/// Supports [`Extrapolate::Enable`](crate::interpolator::Extrapolate::Enable): -/// evaluation beyond the grid extends the boundary cubic polynomials. -/// -/// # Example -/// ``` -/// use ndarray::prelude::*; -/// use ninterp::prelude::*; -/// -/// // f(x) = 2x + 1 (linear: reproduced exactly by any spline) -/// let interp: Interp1D = Interp1D::new( -/// array![0., 1., 2., 3.], -/// array![1., 3., 5., 7.], -/// strategy::CubicC2::not_a_knot(), -/// Extrapolate::Enable, -/// ) -/// .unwrap(); -/// assert_eq!(interp.interpolate(&[1.5]).unwrap(), 4.0); -/// assert_eq!(interp.interpolate(&[4.0]).unwrap(), 9.0); // extrapolation -/// ``` -#[derive(Debug, Clone, PartialEq)] -#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] -#[cfg_attr( - feature = "serde", - serde(bound( - serialize = "T: Serialize + Zero", - deserialize = "T: Deserialize<'de> + Zero" - )) -)] -pub struct CubicC2 { - /// Boundary conditions, one per dimension or a single entry broadcast to all. - // Serializes under the key "CubicC2" rather than "boundary_conditions", making the - // strategy type explicit in the output, consistent with how Linear, Nearest, Step, - // etc. serialize to their type name. - #[cfg_attr(feature = "serde", serde(rename = "CubicC2"))] - pub boundary_conditions: Broadcastable>, - /// Precomputed derivative data, populated by `Strategy1D`/`2D`/`3D`/`ND::init`. Its - /// shape depends on which of those populated it: - /// - /// - `Strategy1D`: cached second derivatives (`M[i] = S''(x_i)`), one [`compute_m`] - /// result for the single 1-D pencil, i.e. shape `[n + 1]` for `n` intervals. - /// - `Strategy2D`/`Strategy3D`/`StrategyND`: the full corner-derivative tensor, - /// shaped like the value grid with one extra trailing axis of length `2^N` (`N` = - /// the grid's dimensionality): see [`compute_corner_cache`]. Every query is then - /// an O(1) lookup, no solving. - /// - /// Not included in the serialized form. After deserializing, call the - /// interpolator's `init_strategy` method (e.g. - /// [`Interp1D::init_strategy`](crate::interpolator::Interp1D::init_strategy)) to - /// recompute this before use. - #[cfg_attr(feature = "serde", serde(skip, default = "empty_cache"))] - pub(crate) cache: ArrayD, -} - -/// Boundary conditions for [`CubicC2`]. -/// -/// [`Endpoints::lower`](CubicC2BoundaryConditions::Endpoints)/`upper` are independent, so -/// mixing types (e.g. [`NotAKnot`](CubicC2Endpoint::NotAKnot) on one side, -/// [`FirstDerivative`](CubicC2Endpoint::FirstDerivative) on the other) is allowed. The -/// common symmetric cases have shorthand constructors: [`not_a_knot`](Self::not_a_knot), -/// [`first_derivative`](Self::first_derivative), [`second_derivative`](Self::second_derivative); -/// [`CubicC2::natural`] additionally shorthands `second_derivative(0, 0)`. -/// -/// Serializes as a bare string for [`NotAKnot`](CubicC2Endpoint::NotAKnot) (both -/// endpoints) and `Periodic`, matching their pre-`Endpoints` representation; a bare -/// `"Natural"` string is accepted and produced as a shorthand for symmetric zero -/// [`SecondDerivative`](CubicC2Endpoint::SecondDerivative) too. Anything else (an -/// explicit value, or asymmetric endpoints) uses the general `{"Endpoints": {"lower": -/// ..., "upper": ...}}` form. -#[derive(Debug, Clone, PartialEq)] -pub enum CubicC2BoundaryConditions { - /// Condition applied independently at each end of the axis. - Endpoints { - /// Condition at the lower endpoint. - lower: CubicC2Endpoint, - /// Condition at the upper endpoint. - upper: CubicC2Endpoint, - }, - /// First and second derivatives match at both endpoints. By convention - /// `values[n]` (the last point along this axis) should equal `values[0]`, since - /// that's what makes the axis periodic. This isn't enforced: `values[n]` is read and - /// used like any other data point (both to build the periodic system and to - /// evaluate the last interval), just never compared against `values[0]`. - /// - /// Deliberately not validated, even approximately: unlike grid coordinates (usually - /// synthetic, so float rounding is the only source of nonuniformity), `values` is - /// often real measured data, where the two ends of a period can legitimately differ - /// by far more than rounding error for reasons that have nothing to do with a - /// mistake (sensor noise, distinct samples at each end of the period, etc.). A - /// tolerance tight enough to catch a genuine mismatch would also reject that valid - /// data. - /// - /// If `values[n] != values[0]`, the spline within `[x[0], x[n]]` itself is still - /// perfectly smooth (both derivatives matched at the endpoints, same as if the axis - /// really were periodic) and passes through every supplied value exactly, `values[n]` - /// included. The mismatch only shows up if the axis is then treated as periodic - /// beyond its own bounds, e.g. via [`Extrapolate::Wrap`](crate::interpolator::Extrapolate::Wrap): - /// crossing the seam is a jump in *value* of exactly `values[n] - values[0]`, with - /// slope and curvature matching continuously on both sides of it. - Periodic, -} - -/// A single endpoint's condition, for [`CubicC2BoundaryConditions::Endpoints`]. -#[derive(Debug, Clone, PartialEq)] -#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] -pub enum CubicC2Endpoint { - /// C³ continuity at the second (from this end) knot; no extra input required. - /// Generally gives better accuracy than natural (a zero - /// [`SecondDerivative`](Self::SecondDerivative)) for smooth functions. Requires at - /// least 3 grid points along this axis (4 if both endpoints are `NotAKnot`). - NotAKnot, - /// Specified first derivative at this endpoint ("clamped"). - FirstDerivative(T), - /// Specified second derivative at this endpoint. Zero, the classic "natural" - /// condition, has a shorthand: [`CubicC2::natural`]. - SecondDerivative(T), -} - -impl CubicC2BoundaryConditions { - /// Not-a-knot at both ends. Requires at least 4 data points per dimension. - pub fn not_a_knot() -> Self { - Self::Endpoints { - lower: CubicC2Endpoint::NotAKnot, - upper: CubicC2Endpoint::NotAKnot, - } - } - - /// Specified first derivative at both endpoints. - pub fn first_derivative(lower: T, upper: T) -> Self { - Self::Endpoints { - lower: CubicC2Endpoint::FirstDerivative(lower), - upper: CubicC2Endpoint::FirstDerivative(upper), - } - } - - /// Specified second derivative at both endpoints. - pub fn second_derivative(lower: T, upper: T) -> Self { - Self::Endpoints { - lower: CubicC2Endpoint::SecondDerivative(lower), - upper: CubicC2Endpoint::SecondDerivative(upper), - } - } -} - -impl From> for CubicC2 { - /// Broadcasts `bc` to all dimensions. - /// - /// Use [`CubicC2::not_a_knot`], [`natural`](CubicC2::natural), - /// [`clamped`](CubicC2::clamped), or [`periodic`](CubicC2::periodic) instead when the - /// desired condition is known at the call site; this is for a `CubicC2BoundaryConditions` - /// value obtained generically (e.g. from runtime config), without matching on it first. - fn from(bc: CubicC2BoundaryConditions) -> Self { - Self { - boundary_conditions: Broadcastable::Broadcast(bc), - cache: empty_cache(), - } - } -} - -impl CubicC2 { - /// Create a cubic spline with a distinct boundary condition per grid dimension. - /// - /// Use [`not_a_knot`](Self::not_a_knot), [`natural`](Self::natural), - /// [`clamped`](Self::clamped), or [`periodic`](Self::periodic) instead when every - /// dimension shares the same condition. - pub fn new(boundary_conditions: Vec>) -> Self { - Self { - boundary_conditions: Broadcastable::Each(boundary_conditions), - cache: empty_cache(), - } - } - - /// Create a cubic spline with not-a-knot boundary conditions. - /// Requires at least 4 data points per dimension. - pub fn not_a_knot() -> Self { - Self { - boundary_conditions: Broadcastable::Broadcast(CubicC2BoundaryConditions::not_a_knot()), - cache: empty_cache(), - } - } - - /// Create a cubic spline with natural (zero second derivative at endpoints) BCs. - pub fn natural() -> Self - where - T: Zero, - { - Self { - boundary_conditions: Broadcastable::Broadcast( - CubicC2BoundaryConditions::second_derivative(T::zero(), T::zero()), - ), - cache: empty_cache(), - } - } - - /// Create a cubic spline with specified first derivatives at both endpoints. - pub fn clamped(lower: T, upper: T) -> Self { - Self { - boundary_conditions: Broadcastable::Broadcast( - CubicC2BoundaryConditions::first_derivative(lower, upper), - ), - cache: empty_cache(), - } - } - - /// Create a cubic spline with periodic boundary conditions. By convention - /// `values[n]` (the last point along each periodic axis) should equal - /// `values[0]`, though this isn't enforced. - pub fn periodic() -> Self { - Self { - boundary_conditions: Broadcastable::Broadcast(CubicC2BoundaryConditions::Periodic), - cache: empty_cache(), - } - } -} - -/// Thomas algorithm (tridiagonal matrix algorithm). Solves `A * x = rhs`. -/// `sub.len() == sup.len() == diag.len() - 1`. -pub(crate) fn thomas(sub: &[T], diag: &[T], sup: &[T], rhs: &[T]) -> Vec { - let n = diag.len(); - let mut cp = vec![T::zero(); n]; - let mut dp = vec![T::zero(); n]; - cp[0] = if n > 1 { sup[0] / diag[0] } else { T::zero() }; - dp[0] = rhs[0] / diag[0]; - for k in 1..n { - let w = diag[k] - sub[k - 1] * cp[k - 1]; - cp[k] = if k < n - 1 { sup[k] / w } else { T::zero() }; - dp[k] = (rhs[k] - sub[k - 1] * dp[k - 1]) / w; - } - let mut x = vec![T::zero(); n]; - x[n - 1] = dp[n - 1]; - for k in (0..n - 1).rev() { - x[k] = dp[k] - cp[k] * x[k + 1]; - } - x -} - -/// Sherman-Morrison cyclic tridiagonal solver. -/// Corner elements `corner` appear at `(0, n-1)` and `(n-1, 0)`. -/// `sub.len() == sup.len() == n - 1`. -pub(crate) fn cyclic_thomas( - sub: &[T], - diag: &[T], - sup: &[T], - rhs: &[T], - corner: T, -) -> Vec { - let n = diag.len(); - if n == 1 { - return vec![rhs[0] / (diag[0] + corner + corner)]; - } - let gamma = -diag[0]; - let c_over_g = corner / gamma; - let mut diag_mod = diag.to_vec(); - diag_mod[0] = diag_mod[0] - gamma; - diag_mod[n - 1] = diag_mod[n - 1] - corner * corner / gamma; - let y = thomas(sub, &diag_mod, sup, rhs); - let mut u_vec = vec![T::zero(); n]; - u_vec[0] = gamma; - u_vec[n - 1] = corner; - let z = thomas(sub, &diag_mod, sup, &u_vec); - let vt_y = y[0] + c_over_g * y[n - 1]; - let vt_z = z[0] + c_over_g * z[n - 1]; - let factor = vt_y / (T::one() + vt_z); - y.into_iter() - .zip(z.iter()) - .map(|(yi, zi)| yi - factor * *zi) - .collect() -} - -/// Computes the second-derivative vector `M[0..=n]` for the given cubic spline BC. -/// -/// Used by [`compute_m_cache`] (called from `Strategy1D::init`, stored in `CubicC2::cache`) -/// and [`corner_cache_axis_pass`] (called from [`compute_corner_cache`]). -pub(crate) fn compute_m( - x: ArrayView1, - y: ArrayView1, - bc: &CubicC2BoundaryConditions, -) -> Vec { - let n = x.len() - 1; - let two = T::one() + T::one(); - let six = two + two + two; - let h: Vec = (0..n).map(|i| x[i + 1] - x[i]).collect(); - let slopes: Vec = (0..n).map(|i| (y[i + 1] - y[i]) / h[i]).collect(); - let u: Vec = (0..n.saturating_sub(1)) - .map(|k| six * (slopes[k + 1] - slopes[k])) - .collect(); - - match bc { - CubicC2BoundaryConditions::Endpoints { lower, upper } => { - // Each endpoint contributes either a direct boundary row (`Some`, referencing - // only its own knot and its immediate interior neighbor) or, for `NotAKnot`, - // `None`: that side's M is eliminated via the third-derivative-continuity - // relation at knot 1 (or n-1), folded into the interior equation there instead - // of appearing as its own unknown. `lower_row`/`upper_row` give - // `(diag, off_diagonal, rhs)` for the direct case. - let lower_row = match lower { - CubicC2Endpoint::NotAKnot => None, - CubicC2Endpoint::SecondDerivative(v) => Some((T::one(), T::zero(), *v)), - CubicC2Endpoint::FirstDerivative(v) => { - Some((two * h[0], h[0], six * (slopes[0] - *v))) - } - }; - let upper_row = match upper { - CubicC2Endpoint::NotAKnot => None, - CubicC2Endpoint::SecondDerivative(v) => Some((T::one(), T::zero(), *v)), - CubicC2Endpoint::FirstDerivative(v) => { - Some((two * h[n - 1], h[n - 1], six * (*v - slopes[n - 1]))) - } - }; - - // `validate_bc_min_points` guarantees `count >= 2` (grid_len >= 3 with one - // `NotAKnot` side, >= 4 with both), so the `j == 0` and `j == count - 1` - // branches below are always distinct rows. - let start = if lower_row.is_some() { 0 } else { 1 }; - let end = if upper_row.is_some() { n } else { n - 1 }; - let count = end - start + 1; - - let mut sub = Vec::with_capacity(count - 1); - let mut diag = Vec::with_capacity(count); - let mut sup = Vec::with_capacity(count - 1); - let mut rhs = Vec::with_capacity(count); - for j in 0..count { - let i = start + j; - if j == 0 { - if let Some((diag_v, sup_v, rhs_v)) = lower_row { - diag.push(diag_v); - sup.push(sup_v); - rhs.push(rhs_v); - } else { - // Folds the eliminated M[0] into knot 1's interior equation. - diag.push((h[0] + h[1]) * (h[0] + two * h[1])); - sup.push(h[1] * h[1] - h[0] * h[0]); - rhs.push(h[1] * u[0]); - } - } else if j == count - 1 { - if let Some((diag_v, sub_v, rhs_v)) = upper_row { - sub.push(sub_v); - diag.push(diag_v); - rhs.push(rhs_v); - } else { - // Mirrors the lower fold, at knot n-1 eliminating M[n]. - sub.push(h[n - 2] * h[n - 2] - h[n - 1] * h[n - 1]); - diag.push((h[n - 2] + h[n - 1]) * (two * h[n - 2] + h[n - 1])); - rhs.push(h[n - 2] * u[n - 2]); - } - } else { - sub.push(h[i - 1]); - diag.push(two * (h[i - 1] + h[i])); - sup.push(h[i]); - rhs.push(u[i - 1]); - } - } - let inner = thomas(&sub, &diag, &sup, &rhs); - - let mut m = Vec::with_capacity(n + 1); - if lower_row.is_none() { - m.push(((h[0] + h[1]) * inner[0] - h[0] * inner[1]) / h[1]); - } - m.extend_from_slice(&inner); - if upper_row.is_none() { - let last = inner.len() - 1; - m.push( - ((h[n - 2] + h[n - 1]) * inner[last] - h[n - 1] * inner[last - 1]) / h[n - 2], - ); - } - m - } - CubicC2BoundaryConditions::Periodic => { - // `y[n]` is read below, via `slopes[n - 1]` in the `rhs[0]` line, and - // `evaluate_spline_from_m` also reads `y[n]` directly when evaluating the last - // interval. By convention `y[n]` should equal `y[0]`, but nothing here - // compares them; `y[n]` is just used as ordinary data either way. - if n < 2 { - vec![T::zero(); n + 1] - } else { - let sub_sup = h[..n - 1].to_vec(); - let mut diag = vec![two * (h[n - 1] + h[0])]; - for k in 1..n { - diag.push(two * (h[k - 1] + h[k])); - } - let mut rhs = vec![six * (slopes[0] - slopes[n - 1])]; - rhs.extend_from_slice(&u); - let corner = h[n - 1]; - let mut m_vals = cyclic_thomas(&sub_sup, &diag, &sub_sup, &rhs, corner); - let m0 = m_vals[0]; - m_vals.push(m0); - m_vals - } - } - } -} - -/// Evaluates the M-form cubic spline at `point` using precomputed second derivatives `m`. -pub(crate) fn evaluate_spline_from_m( - x: ArrayView1, - y: ArrayView1, - m: ArrayView1, - point: T, -) -> T { - let two = T::one() + T::one(); - let six = two + two + two; - let i = locate_lower_index(x, &point); - let h = x[i + 1] - x[i]; - let dx = point - x[i]; - let dx_r = h - dx; - let six_h = six * h; - let h2_over_six = h * h / six; - m[i] * dx_r * dx_r * dx_r / six_h - + m[i + 1] * dx * dx * dx / six_h - + (y[i] - m[i] * h2_over_six) * dx_r / h - + (y[i + 1] - m[i + 1] * h2_over_six) * dx / h -} - -/// Checks `grid_len` against boundary condition `bc`'s minimum point requirement (e.g. -/// [`CubicC2Endpoint::NotAKnot`] needs at least 3 grid points on its own side, 4 if both -/// endpoints are `NotAKnot`), ahead of the real work in [`compute_m`]. -/// -/// Pure, no mutation; used by each dimensionality's `Strategy*D::validate`. -pub(crate) fn validate_bc_min_points( - bc: &CubicC2BoundaryConditions, - grid_len: usize, - dim: usize, -) -> Result<(), ValidateError> { - let CubicC2BoundaryConditions::Endpoints { lower, upper } = bc else { - return Ok(()); - }; - let min = match ( - matches!(lower, CubicC2Endpoint::NotAKnot), - matches!(upper, CubicC2Endpoint::NotAKnot), - ) { - (true, true) => 4, - (true, false) | (false, true) => 3, - (false, false) => return Ok(()), - }; - if grid_len < min { - return Err(ValidateError::Other(format!( - "CubicC2: dim {dim} has {grid_len} grid points; NotAKnot requires at least {min}" - ))); - } - Ok(()) -} - -/// Computes and caches `M[0..=n]` (`CubicC2::cache`) for [`Strategy1D::init`], so -/// [`evaluate_spline_1d_cached`] can look them up in O(1) instead of re-solving on every -/// `interpolate` call. -pub(crate) fn compute_m_cache( - x: ArrayView1, - y: ArrayView1, - bc: &CubicC2BoundaryConditions, -) -> ArrayD { - let m = compute_m(x, y, bc); - ArrayD::from_shape_vec(IxDyn(&[m.len()]), m) - .expect("compute_m's output length matches its own shape") -} -/// Evaluates [`Strategy1D`]'s cached spline (`m_cache`, from [`compute_m_cache`]) at `point`. -pub(crate) fn evaluate_spline_1d_cached( - x: ArrayView1, - y: ArrayView1, - m_cache: ArrayViewD, - point: T, -) -> Result { - let m = m_cache.into_dimensionality::().map_err(|_| { - InterpolateError::Other( - "internal: non-1-D m_cache, Strategy1D::cache invariant broken".into(), - ) - })?; - Ok(evaluate_spline_from_m(x, y, m, point)) -} +pub use c1::{CubicC1, CubicC1CacheMode, CubicC1DerivativeMode}; +pub use c2::{CubicC2, CubicC2BoundaryConditions, CubicC2Endpoint}; -/// Closed-form first derivative `S'(x_i)` at every knot, from an already-solved moment -/// vector `m` (no extra solve). The companion to [`compute_m`], used to build -/// [`compute_corner_cache`]'s derivative fields. -pub(crate) fn knot_derivatives_from_m( - x: ArrayView1, - y: ArrayView1, - m: ArrayView1, -) -> Vec { - let n = x.len() - 1; - let two = T::one() + T::one(); - let six = two + two + two; - let mut d: Vec = (0..n) - .map(|i| { - let h = x[i + 1] - x[i]; - (y[i + 1] - y[i]) / h - h * (two * m[i] + m[i + 1]) / six - }) - .collect(); - let h_last = x[n] - x[n - 1]; - d.push((y[n] - y[n - 1]) / h_last + h_last * (m[n - 1] + two * m[n]) / six); - d -} - -/// Splines every 1-D lane of `field` along `axis` and replaces it with its knot -/// derivatives (via [`compute_m`] + [`knot_derivatives_from_m`]), returning a new array -/// the same shape as `field`, for [`compute_corner_cache`]. -fn corner_cache_axis_pass( - grid: ArrayView1, - field: ArrayViewD, - axis: usize, - bc: &CubicC2BoundaryConditions, -) -> ArrayD { - let axis = Axis(axis); - let mut out = ArrayD::::zeros(IxDyn(field.shape())); - for (y, mut out_lane) in field.lanes(axis).into_iter().zip(out.lanes_mut(axis)) { - let m = compute_m(grid, y, bc); - let d = knot_derivatives_from_m(grid, y, ArrayView1::from(&m)); - out_lane.assign(&ArrayView1::from(&d)); - } - out -} - -/// Precomputes the full corner-derivative tensor for [`CubicC2`]'s `Strategy2D`/ -/// `Strategy3D` full-cache upgrade: for every grid point, all `2^N` partial-derivative -/// combinations (value, first partials, and mixed partials) needed to evaluate a Hermite -/// patch in O(1) via [`evaluate_spline_corner_cached`]. -/// -/// Returns an array shaped like `values`, with one extra trailing axis of length `2^N` -/// (`N = grids.len()`). Index `k` in that axis is a bitmask over the `N` grid axes: bit -/// `i` set means "this entry is `d/dx_i` of the value", so `k == 0` is the raw value and -/// `k == 2^N - 1` is the full mixed partial. -/// -/// Built by processing axes `0..N` in order: for every mask not yet including axis `i`'s -/// bit, splines the corresponding field along axis `i` ([`corner_cache_axis_pass`]) to -/// populate the entry with that bit added. Boundary condition per axis is -/// `bcs[axis]` when splining the raw values (`mask == 0`, this axis's own -/// first-derivative pass). For every other `mask` (a cross-derivative pass over an -/// already-differentiated field, not the original data), each endpoint's -/// [`FirstDerivative`](CubicC2Endpoint::FirstDerivative)/[`SecondDerivative`](CubicC2Endpoint::SecondDerivative) -/// falls back to the same-order condition with a zero value: a `FirstDerivative`/ -/// `SecondDerivative` endpoint fixes the same scalar derivative at every point along that -/// boundary, so differentiating that (constant-along-the-boundary) field with respect to -/// any other axis must itself be zero at that endpoint, at the same order. `NotAKnot` -/// endpoints pass through unchanged (they carry no value to homogenize). Verified -/// empirically: with non-separable data, only these homogeneous fallbacks make this -/// cached path agree (to float precision) with `StrategyND`'s independent recursive -/// solve; leaving a nonzero endpoint value unchanged in a cross-derivative pass does not. -pub(crate) fn compute_corner_cache( - grids: &[ArrayView1], - values: ArrayViewD, - bcs: &Broadcastable>, -) -> ArrayD { - let n_axes = grids.len(); - let n_bits = 1usize << n_axes; - let mask_axis = Axis(n_axes); - - let mut out_shape = values.shape().to_vec(); - out_shape.push(n_bits); - let mut cache = ArrayD::::zeros(IxDyn(&out_shape)); - cache.index_axis_mut(mask_axis, 0).assign(&values); - - for (axis, grid) in grids.iter().enumerate() { - let bit = 1usize << axis; - let bc_axis = &bcs[axis]; - let cross_endpoint = |e: &CubicC2Endpoint| match e { - CubicC2Endpoint::FirstDerivative(_) => CubicC2Endpoint::FirstDerivative(T::zero()), - CubicC2Endpoint::SecondDerivative(_) => CubicC2Endpoint::SecondDerivative(T::zero()), - other @ CubicC2Endpoint::NotAKnot => other.clone(), - }; - let cross_bc = match bc_axis { - CubicC2BoundaryConditions::Endpoints { lower, upper } => { - CubicC2BoundaryConditions::Endpoints { - lower: cross_endpoint(lower), - upper: cross_endpoint(upper), - } - } - CubicC2BoundaryConditions::Periodic => CubicC2BoundaryConditions::Periodic, - }; - for mask in 0..bit { - let bc = if mask == 0 { bc_axis } else { &cross_bc }; - let field = cache.index_axis(mask_axis, mask).to_owned(); - let deriv = corner_cache_axis_pass(*grid, field.view(), axis, bc); - cache.index_axis_mut(mask_axis, mask | bit).assign(&deriv); - } - } - cache -} - -#[cfg(feature = "serde")] -mod cubic_serde { - use super::*; - use serde::Deserializer; - - /// Bare-string wire form for the two `CubicC2BoundaryConditions::Endpoints` cases - /// that carry no distinguishing per-endpoint data (symmetric `NotAKnot`, symmetric - /// zero `SecondDerivative`, i.e. "natural"), plus `Periodic`. Anything else uses the - /// general `Endpoints` shape instead. - #[derive(Deserialize, Serialize)] - enum BcName { - NotAKnot, - Natural, - Periodic, - } - - /// General `Endpoints` shape, owned form for [`Deserialize`]. - #[derive(Deserialize)] - struct EndpointsOwned { - lower: CubicC2Endpoint, - upper: CubicC2Endpoint, - } - - /// General `Endpoints` shape, borrowed form for [`Serialize`]: avoids requiring `T: - /// Clone` just to serialize. - #[derive(Serialize)] - struct EndpointsRef<'a, T> { - lower: &'a CubicC2Endpoint, - upper: &'a CubicC2Endpoint, - } - - #[derive(Deserialize)] - #[serde(untagged)] - enum BcWireDe { - Named(BcName), - Endpoints { - #[serde(rename = "Endpoints")] - endpoints: EndpointsOwned, - }, - } - - #[derive(Serialize)] - #[serde(untagged)] - enum BcWireSer<'a, T> { - Named(BcName), - Endpoints { - #[serde(rename = "Endpoints")] - endpoints: EndpointsRef<'a, T>, - }, - } - - impl Serialize for CubicC2BoundaryConditions { - fn serialize(&self, serializer: S) -> Result { - match self { - Self::Periodic => BcWireSer::::Named(BcName::Periodic).serialize(serializer), - Self::Endpoints { - lower: CubicC2Endpoint::NotAKnot, - upper: CubicC2Endpoint::NotAKnot, - } => BcWireSer::::Named(BcName::NotAKnot).serialize(serializer), - Self::Endpoints { - lower: CubicC2Endpoint::SecondDerivative(lower), - upper: CubicC2Endpoint::SecondDerivative(upper), - } if lower.is_zero() && upper.is_zero() => { - BcWireSer::::Named(BcName::Natural).serialize(serializer) - } - Self::Endpoints { lower, upper } => BcWireSer::Endpoints { - endpoints: EndpointsRef { lower, upper }, - } - .serialize(serializer), - } - } - } - - impl<'de, T: Deserialize<'de> + Zero> Deserialize<'de> for CubicC2BoundaryConditions { - fn deserialize>(deserializer: D) -> Result { - Ok(match BcWireDe::::deserialize(deserializer)? { - BcWireDe::Named(BcName::NotAKnot) => CubicC2BoundaryConditions::not_a_knot(), - BcWireDe::Named(BcName::Natural) => Self::Endpoints { - lower: CubicC2Endpoint::SecondDerivative(T::zero()), - upper: CubicC2Endpoint::SecondDerivative(T::zero()), - }, - BcWireDe::Named(BcName::Periodic) => Self::Periodic, - BcWireDe::Endpoints { endpoints } => Self::Endpoints { - lower: endpoints.lower, - upper: endpoints.upper, - }, - }) - } - } -} - -#[cfg(test)] -mod tests { - use super::*; - - #[test] - fn from_bc_matches_named_constructors() { - assert_eq!( - CubicC2::from(CubicC2BoundaryConditions::::not_a_knot()), - CubicC2::not_a_knot() - ); - assert_eq!( - CubicC2::from(CubicC2BoundaryConditions::::second_derivative( - 0.0, 0.0 - )), - CubicC2::natural() - ); - assert_eq!( - CubicC2::from(CubicC2BoundaryConditions::first_derivative(1.0, 2.0)), - CubicC2::clamped(1.0, 2.0) - ); - assert_eq!( - CubicC2::from(CubicC2BoundaryConditions::::Periodic), - CubicC2::periodic() - ); - } - - /// Builds the `(n+1) x (n+1)` moment system directly from `bc`'s mathematical - /// definition, independent of `compute_m`'s reduced/eliminated form, to cross-check - /// it on asymmetric endpoint combinations that this crate's other tests (all - /// symmetric until now) wouldn't have caught a mismatched row for. - fn dense_m_system( - x: &[f64], - y: &[f64], - bc: &CubicC2BoundaryConditions, - ) -> (Vec>, Vec) { - let n = x.len() - 1; - let h: Vec = (0..n).map(|i| x[i + 1] - x[i]).collect(); - let slopes: Vec = (0..n).map(|i| (y[i + 1] - y[i]) / h[i]).collect(); - let mut a = vec![vec![0.0; n + 1]; n + 1]; - let mut b = vec![0.0; n + 1]; - - for i in 1..n { - a[i][i - 1] = h[i - 1]; - a[i][i] = 2.0 * (h[i - 1] + h[i]); - a[i][i + 1] = h[i]; - b[i] = 6.0 * (slopes[i] - slopes[i - 1]); - } - - let CubicC2BoundaryConditions::Endpoints { lower, upper } = bc else { - panic!("dense_m_system only supports Endpoints for this test"); - }; - match lower { - CubicC2Endpoint::NotAKnot => { - a[0][0] = h[1]; - a[0][1] = -(h[0] + h[1]); - a[0][2] = h[0]; - } - CubicC2Endpoint::SecondDerivative(v) => { - a[0][0] = 1.0; - b[0] = *v; - } - CubicC2Endpoint::FirstDerivative(v) => { - a[0][0] = 2.0 * h[0]; - a[0][1] = h[0]; - b[0] = 6.0 * (slopes[0] - v); - } - } - match upper { - CubicC2Endpoint::NotAKnot => { - a[n][n] = h[n - 2]; - a[n][n - 1] = -(h[n - 2] + h[n - 1]); - a[n][n - 2] = h[n - 1]; - } - CubicC2Endpoint::SecondDerivative(v) => { - a[n][n] = 1.0; - b[n] = *v; - } - CubicC2Endpoint::FirstDerivative(v) => { - a[n][n] = 2.0 * h[n - 1]; - a[n][n - 1] = h[n - 1]; - b[n] = 6.0 * (v - slopes[n - 1]); - } - } - (a, b) - } - - /// Plain Gaussian elimination with partial pivoting, independent of `thomas`/ - /// `cyclic_thomas`, for cross-checking their output. - fn dense_solve(mut a: Vec>, mut b: Vec) -> Vec { - let n = b.len(); - for col in 0..n { - let pivot = (col..n) - .max_by(|&r1, &r2| a[r1][col].abs().partial_cmp(&a[r2][col].abs()).unwrap()) - .unwrap(); - a.swap(col, pivot); - b.swap(col, pivot); - for row in (col + 1)..n { - let factor = a[row][col] / a[col][col]; - let pivot_row = a[col].clone(); - for (c, pivot_val) in pivot_row.iter().enumerate().skip(col) { - a[row][c] -= factor * pivot_val; - } - b[row] -= factor * b[col]; - } - } - let mut x = vec![0.0; n]; - for row in (0..n).rev() { - let sum: f64 = (row + 1..n).map(|c| a[row][c] * x[c]).sum(); - x[row] = (b[row] - sum) / a[row][row]; - } - x - } - - #[test] - fn compute_m_matches_dense_solve_for_mixed_endpoints() { - // Non-uniform grid, non-polynomial data: exercises the linear algebra for real, - // rather than something any BC would reproduce exactly regardless of correctness. - let x = [0.0, 0.4, 1.1, 1.8, 2.9, 4.0]; - let y = [0.5, 1.2, 0.3, 2.1, 1.0, 3.3]; - let xv = ArrayView1::from(&x); - let yv = ArrayView1::from(&y); - - let endpoints = [ - CubicC2Endpoint::NotAKnot, - CubicC2Endpoint::SecondDerivative(0.0), - CubicC2Endpoint::SecondDerivative(2.5), - CubicC2Endpoint::FirstDerivative(0.0), - CubicC2Endpoint::FirstDerivative(-1.5), - ]; - for lower in &endpoints { - for upper in &endpoints { - let bc = CubicC2BoundaryConditions::Endpoints { - lower: lower.clone(), - upper: upper.clone(), - }; - let got = compute_m(xv, yv, &bc); - let (a, b) = dense_m_system(&x, &y, &bc); - let expected = dense_solve(a, b); - for (g, e) in got.iter().zip(expected.iter()) { - assert_approx_eq!(*g, *e, 1e-9); - } - } - } - } - - #[test] - #[cfg(feature = "serde")] - fn test_serde() { - // The whole-axis cases collapse to a bare string, keyed by strategy name like - // every other built-in strategy. - assert_eq!( - serde_json::to_string(&CubicC2::::not_a_knot()).unwrap(), - r#"{"CubicC2":"NotAKnot"}"# - ); - assert_eq!( - serde_json::to_string(&CubicC2::::natural()).unwrap(), - r#"{"CubicC2":"Natural"}"# - ); - assert_eq!( - serde_json::to_string(&CubicC2::::periodic()).unwrap(), - r#"{"CubicC2":"Periodic"}"# - ); - // Anything carrying an explicit value, or mixing endpoint types, uses the - // general `Endpoints` form instead. - assert_eq!( - serde_json::to_string(&CubicC2::::clamped(1.0, 2.0)).unwrap(), - r#"{"CubicC2":{"Endpoints":{"lower":{"FirstDerivative":1.0},"upper":{"FirstDerivative":2.0}}}}"# - ); - let mixed = CubicC2::::from(CubicC2BoundaryConditions::Endpoints { - lower: CubicC2Endpoint::NotAKnot, - upper: CubicC2Endpoint::FirstDerivative(3.0), - }); - assert_eq!( - serde_json::to_string(&mixed).unwrap(), - r#"{"CubicC2":{"Endpoints":{"lower":"NotAKnot","upper":{"FirstDerivative":3.0}}}}"# - ); - assert_eq!( - serde_json::from_str::>(&serde_json::to_string(&mixed).unwrap()).unwrap(), - mixed - ); - - // Bare "NotAKnot"/"Natural" round-trip; the fully expanded `Endpoints` form is - // also accepted on read for both, even though only the bare form is ever written. - assert_eq!( - serde_json::from_str::>(r#"{"CubicC2":"NotAKnot"}"#).unwrap(), - CubicC2::::not_a_knot() - ); - assert_eq!( - serde_json::from_str::>(r#"{"CubicC2":"Natural"}"#).unwrap(), - CubicC2::::natural() - ); - assert_eq!( - serde_json::from_str::>( - r#"{"CubicC2":{"Endpoints":{"lower":"NotAKnot","upper":"NotAKnot"}}}"# - ) - .unwrap(), - CubicC2::::not_a_knot() - ); - assert_eq!( - serde_json::from_str::>( - r#"{"CubicC2":{"Endpoints":{"lower":{"SecondDerivative":0.0},"upper":{"SecondDerivative":0.0}}}}"# - ) - .unwrap(), - CubicC2::::natural() - ); - } -} +pub(crate) use c1::{ + compute_corner_cache_fd, compute_fd_cache, evaluate_hermite_1d_cached, + evaluate_spline_corner_local, +}; +pub(crate) use c2::{ + compute_corner_cache, compute_m_cache, evaluate_spline_1d_cached, validate_bc_min_points, +}; +pub(crate) use utils::evaluate_spline_corner_cached; diff --git a/src/strategy/cubic/utils.rs b/src/strategy/cubic/utils.rs index 32c41c4..c1e1e26 100644 --- a/src/strategy/cubic/utils.rs +++ b/src/strategy/cubic/utils.rs @@ -1,6 +1,6 @@ -//! Corner-derivative-tensor evaluation, shared by [`CubicC2`] and (once it exists) -//! `CubicC1`: BC-agnostic, works the same regardless of how the tensor's entries were -//! populated (a solve for `CubicC2`, finite differences for `CubicC1`). +//! Corner-derivative-tensor evaluation, shared by [`c1`](super::c1) and +//! [`c2`](super::c2): BC-agnostic, works the same regardless of how the tensor's +//! entries were populated (a solve for `CubicC2`, finite differences for `CubicC1`). use super::*; use ndarray::Zip; @@ -8,7 +8,7 @@ use ndarray::Zip; /// Standard cubic Hermite basis blend of endpoint values `p0`/`p1` and their derivatives /// `m0`/`m1` (actual derivatives, not yet scaled by the interval width `h`) at /// fractional position `t` in `[0, 1]` between them. -fn evaluate_hermite_1d(p0: T, m0: T, p1: T, m1: T, h: T, t: T) -> T { +pub(crate) fn evaluate_hermite_1d(p0: T, m0: T, p1: T, m1: T, h: T, t: T) -> T { let two = T::one() + T::one(); let three = two + T::one(); let t2 = t * t; diff --git a/src/strategy/enums/n.rs b/src/strategy/enums/n.rs index c124a21..3937232 100644 --- a/src/strategy/enums/n.rs +++ b/src/strategy/enums/n.rs @@ -11,6 +11,7 @@ strategy_enum_impl!( (Step, strategy::Step), (Linear, strategy::Linear), (LinearUniform, strategy::LinearUniform), + (CubicC1, strategy::CubicC1), (CubicC2, strategy::CubicC2), (GridTransform, strategy::GridTransform>>), (ValuesTransform, strategy::ValuesTransform>>), diff --git a/src/strategy/enums/one.rs b/src/strategy/enums/one.rs index eff1892..ff38e37 100644 --- a/src/strategy/enums/one.rs +++ b/src/strategy/enums/one.rs @@ -11,6 +11,7 @@ strategy_enum_impl!( (Step, strategy::Step), (Linear, strategy::Linear), (LinearUniform, strategy::LinearUniform), + (CubicC1, strategy::CubicC1), (CubicC2, strategy::CubicC2), (GridTransform, strategy::GridTransform>>), (ValuesTransform, strategy::ValuesTransform>>), diff --git a/src/strategy/enums/three.rs b/src/strategy/enums/three.rs index 893ac5c..b943feb 100644 --- a/src/strategy/enums/three.rs +++ b/src/strategy/enums/three.rs @@ -11,6 +11,7 @@ strategy_enum_impl!( (Step, strategy::Step), (Linear, strategy::Linear), (LinearUniform, strategy::LinearUniform), + (CubicC1, strategy::CubicC1), (CubicC2, strategy::CubicC2), (GridTransform, strategy::GridTransform>>), (ValuesTransform, strategy::ValuesTransform>>), diff --git a/src/strategy/enums/two.rs b/src/strategy/enums/two.rs index ca3c4fd..8dae8d1 100644 --- a/src/strategy/enums/two.rs +++ b/src/strategy/enums/two.rs @@ -11,6 +11,7 @@ strategy_enum_impl!( (Step, strategy::Step), (Linear, strategy::Linear), (LinearUniform, strategy::LinearUniform), + (CubicC1, strategy::CubicC1), (CubicC2, strategy::CubicC2), (GridTransform, strategy::GridTransform>>), (ValuesTransform, strategy::ValuesTransform>>), diff --git a/src/strategy/mod.rs b/src/strategy/mod.rs index 8c1e72f..99f5e6e 100644 --- a/src/strategy/mod.rs +++ b/src/strategy/mod.rs @@ -10,7 +10,7 @@ pub mod broadcast; use broadcast::Broadcastable; pub mod cubic; -pub use cubic::CubicC2; +pub use cubic::{CubicC1, CubicC1CacheMode, CubicC1DerivativeMode, CubicC2}; pub mod step; pub use step::Step; diff --git a/tests/serde_strategies.rs b/tests/serde_strategies.rs index e67b6a7..603085e 100644 --- a/tests/serde_strategies.rs +++ b/tests/serde_strategies.rs @@ -27,6 +27,13 @@ fn cubic_c2_variants() -> Vec> { ] } +fn cubic_c1_variants() -> Vec> { + vec![ + CubicC1::default(), + CubicC1::new().with_cache_mode(CubicC1CacheMode::None), + ] +} + #[test] fn bare_strategies_round_trip() { round_trip(&Nearest); @@ -38,6 +45,9 @@ fn bare_strategies_round_trip() { for bc in cubic_c2_variants() { round_trip(&bc); } + for c1 in cubic_c1_variants() { + round_trip(&c1); + } round_trip(&GridTransform::::log(Linear)); round_trip(&ValuesTransform::::log(Linear)); } @@ -61,6 +71,9 @@ macro_rules! enum_round_trip_test { for bc in cubic_c2_variants() { round_trip(&$Enum::from(bc)); } + for c1 in cubic_c1_variants() { + round_trip(&$Enum::from(c1)); + } let inner: Box<$Enum> = Box::new($Enum::::from(Linear)); round_trip(&$Enum::::from(GridTransform::log(inner.clone()))); round_trip(&$Enum::::from(ValuesTransform::log(inner))); From 51e3a05349ec8f479218d81c2fb03a1f9a1d7184 Mon Sep 17 00:00:00 2001 From: Kyle Carow Date: Sun, 23 Aug 2026 16:37:16 -0600 Subject: [PATCH 2/2] add tests --- src/interpolator/n/tests.rs | 12 ++++++ src/interpolator/one/tests.rs | 25 ++++++++++++ src/interpolator/three/tests.rs | 65 ++++++++++++++++++++++++++++++ src/interpolator/two/tests.rs | 71 +++++++++++++++++++++++++++++++++ 4 files changed, 173 insertions(+) diff --git a/src/interpolator/n/tests.rs b/src/interpolator/n/tests.rs index e61ec4f..56c8049 100644 --- a/src/interpolator/n/tests.rs +++ b/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. diff --git a/src/interpolator/one/tests.rs b/src/interpolator/one/tests.rs index b26d5e2..d49894a 100644 --- a/src/interpolator/one/tests.rs +++ b/src/interpolator/one/tests.rs @@ -317,6 +317,31 @@ fn test_cubic_c1_linear_exact() { 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 diff --git a/src/interpolator/three/tests.rs b/src/interpolator/three/tests.rs index f1a65c7..99c0320 100644 --- a/src/interpolator/three/tests.rs +++ b/src/interpolator/three/tests.rs @@ -415,6 +415,71 @@ fn test_cubic_c2_clamped_uses_given_derivative() { ); } +#[test] +fn test_cubic_c1_knot_exactness() { + // Values at all knots must be reproduced exactly regardless of data shape, under + // both cache modes. + fn f(x: f64, y: f64, z: f64) -> f64 { + x * x * y + y * y * z + z * z * x + } + let grid = [0., 1., 2., 3.]; + let values = Array3::from_shape_fn((4, 4, 4), |(i, j, k)| f(grid[i], grid[j], grid[k])); + for cache_mode in [CubicC1CacheMode::Full, CubicC1CacheMode::None] { + let interp = Interp3D::new( + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + values.clone(), + strategy::CubicC1::new().with_cache_mode(cache_mode), + Extrapolate::Error, + ) + .unwrap(); + for (i, &xi) in grid.iter().enumerate() { + for (j, &yj) in grid.iter().enumerate() { + for (k, &zk) in grid.iter().enumerate() { + assert_approx_eq!( + interp.interpolate(&[xi, yj, zk]).unwrap(), + values[[i, j, k]] + ); + } + } + } + } +} + +#[test] +fn test_cubic_c1_interior_accuracy() { + // Unlike `CubicC2`'s `NotAKnot`, `CubicC1`'s finite differences carry a real error + // term for genuinely nonlinear data, so this checks bounded accuracy against a + // known function, not exact reproduction, under both cache modes. `0.5` is a real, + // checked bound (max observed error ~0.24 at these points). + fn f(x: f64, y: f64, z: f64) -> f64 { + x * x * y + y * y * z + z * z * x + } + let grid = [0., 1., 2., 3.]; + let values = Array3::from_shape_fn((4, 4, 4), |(i, j, k)| f(grid[i], grid[j], grid[k])); + for cache_mode in [CubicC1CacheMode::Full, CubicC1CacheMode::None] { + let interp = Interp3D::new( + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + values.clone(), + strategy::CubicC1::new().with_cache_mode(cache_mode), + Extrapolate::Error, + ) + .unwrap(); + for &(x, y, z) in &[(0.5, 0.5, 0.5), (1.5, 2.5, 0.25), (2.25, 0.75, 1.5)] { + let got = interp.interpolate(&[x, y, z]).unwrap(); + let expected = f(x, y, z); + assert!( + (got - expected).abs() < 0.5, + "f({x}, {y}, {z}) = {expected}, got {got} (diff {})", + (got - expected).abs() + ); + } + } +} + #[test] fn test_cubic_c1_linear_exact() { // Linear data: finite differences recover the exact constant slope on each axis, so diff --git a/src/interpolator/two/tests.rs b/src/interpolator/two/tests.rs index d34a4e0..f22c8fd 100644 --- a/src/interpolator/two/tests.rs +++ b/src/interpolator/two/tests.rs @@ -443,6 +443,77 @@ fn test_cubic_c2_mixed_endpoints_scipy_oracle() { } } +#[test] +fn test_cubic_c1_knot_exactness() { + // Values at all knots must be reproduced exactly regardless of data shape, under + // both cache modes (also exercises `None`'s grid-point-exactness at the boundary, + // where the local window is clipped rather than the usual 4-point interior case). + let grid = array![0., 1., 2., 3.]; + let values = array![ + [0.5, 1.2, 0.3, 2.1], + [1.8, 0.4, 2.5, 1.1], + [0.9, 2.2, 1.4, 0.6], + [2.3, 1.0, 0.7, 1.9], + ]; + for cache_mode in [CubicC1CacheMode::Full, CubicC1CacheMode::None] { + let interp = Interp2D::new( + grid.clone(), + grid.clone(), + values.clone(), + strategy::CubicC1::new().with_cache_mode(cache_mode), + 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_cubic_c1_interior_accuracy() { + // Unlike `CubicC2`'s `NotAKnot` (which reproduces any degree-<=3 polynomial + // exactly), `CubicC1`'s finite-difference derivatives carry an O(h^2) error term + // for genuinely nonlinear data, so this checks bounded accuracy against a known + // function, not exact reproduction, under both cache modes. + fn f(x: f64, y: f64) -> f64 { + x * x * y + x * y * y + } + let grid = array![0., 1., 2., 3.]; + let values = array![ + [0., 0., 0., 0.], + [0., 2., 6., 12.], + [0., 6., 16., 30.], + [0., 12., 30., 54.], + ]; + for cache_mode in [CubicC1CacheMode::Full, CubicC1CacheMode::None] { + let interp = Interp2D::new( + grid.clone(), + grid.clone(), + values.clone(), + strategy::CubicC1::new().with_cache_mode(cache_mode), + Extrapolate::Error, + ) + .unwrap(); + for &(x, y) in &[(0.5, 0.5), (1.5, 2.5), (2.5, 1.5), (0.25, 2.75)] { + let got = interp.interpolate(&[x, y]).unwrap(); + let expected = f(x, y); + assert!( + (got - expected).abs() < 0.5, + "f({x}, {y}) = {expected}, got {got} (diff {})", + (got - expected).abs() + ); + } + } +} + #[test] fn test_cubic_c1_linear_exact() { // Linear data: finite differences recover the exact constant slope on each axis, so