From 5c3e7d093a3e2cd7ceb9f8aaa1bf25627894ae63 Mon Sep 17 00:00:00 2001 From: Kyle Carow Date: Mon, 3 Aug 2026 01:00:55 -0600 Subject: [PATCH 1/2] initial cubic spline implementation --- src/interpolator/n/strategies.rs | 20 ++++++ src/interpolator/n/tests.rs | 42 +++++++++++ src/interpolator/one/strategies.rs | 104 +++++++++++++++++++++++++++ src/interpolator/one/tests.rs | 54 ++++++++++++++ src/interpolator/three/strategies.rs | 24 +++++++ src/interpolator/three/tests.rs | 55 ++++++++++++++ src/interpolator/two/strategies.rs | 24 +++++++ src/interpolator/two/tests.rs | 49 +++++++++++++ src/strategy/mod.rs | 52 ++++++++++++++ src/strategy/traits.rs | 94 ++++++++++++++++++++++++ 10 files changed, 518 insertions(+) diff --git a/src/interpolator/n/strategies.rs b/src/interpolator/n/strategies.rs index 4dff2ab..8dcf512 100644 --- a/src/interpolator/n/strategies.rs +++ b/src/interpolator/n/strategies.rs @@ -225,3 +225,23 @@ where false } } + +impl StrategyND for CubicSpline +where + D: Data + RawDataClone + Clone, + D::Elem: Float + Debug, +{ + fn interpolate( + &self, + data: &InterpDataND, + point: &[D::Elem], + ) -> Result { + let grids: Vec> = data.grid.iter().map(|g| g.view()).collect(); + Ok(spline_eval_nd_recursive(&grids, data.values.view(), point)) + } + + /// Returns `true`: the boundary cubic polynomials extend naturally. + fn allow_extrapolate(&self) -> bool { + true + } +} diff --git a/src/interpolator/n/tests.rs b/src/interpolator/n/tests.rs index 0845bf0..d5e3db1 100644 --- a/src/interpolator/n/tests.rs +++ b/src/interpolator/n/tests.rs @@ -1,5 +1,47 @@ use super::*; +#[test] +fn test_cubic_spline() { + // f(x, y) = 2x + y: linear, reproduced exactly by any spline + let interp = InterpND::new( + vec![array![0., 1., 2.], array![0., 1., 2.]], + array![[0., 1., 2.], [2., 3., 4.], [4., 5., 6.]].into_dyn(), + strategy::CubicSpline::new(), + Extrapolate::Enable, + ) + .unwrap(); + assert_approx_eq!(interp.interpolate(&[0.5, 0.5]).unwrap(), 1.5); + assert_approx_eq!(interp.interpolate(&[1., 1.]).unwrap(), 3.); + assert_approx_eq!(interp.interpolate(&[3., 1.]).unwrap(), 7.); // extrapolation +} + +#[test] +fn test_cubic_spline_knot_exactness() { + let interp = InterpND::new( + vec![array![0., 1., 2., 3.], array![0., 1., 2., 3.]], + array![ + [0., 1., 4., 9.], + [1., 2., 5., 10.], + [4., 5., 8., 13.], + [9., 10., 13., 18.], + ] + .into_dyn(), + strategy::CubicSpline::new(), + Extrapolate::Error, + ) + .unwrap(); + for i in 0..4usize { + for j in 0..4usize { + let xi = interp.data.grid[0][i]; + let yj = interp.data.grid[1][j]; + assert_approx_eq!( + interp.interpolate(&[xi, yj]).unwrap(), + interp.data.values[[i, j]] + ); + } + } +} + #[test] fn test_linear_0d() { let interp = InterpND::new( diff --git a/src/interpolator/one/strategies.rs b/src/interpolator/one/strategies.rs index 625e155..42bd87b 100644 --- a/src/interpolator/one/strategies.rs +++ b/src/interpolator/one/strategies.rs @@ -114,3 +114,107 @@ where false } } + +impl Strategy1D for CubicSpline +where + D: Data + RawDataClone + Clone, + D::Elem: Float + Debug, +{ + /// Computes and caches the per-interval spline coefficients. + /// + /// Uses a natural cubic spline (zero second derivative at endpoints). + /// Solves the resulting tridiagonal system via the Thomas algorithm. + fn init(&mut self, data: &InterpData1D) -> Result<(), ValidateError> { + let x = data.grid[0].view(); + let y = data.values.view(); + let n = x.len() - 1; // number of intervals + + let two = D::Elem::one() + D::Elem::one(); + let six = two + two + two; + + // h[i] = x[i+1] - x[i] + let h: Vec = (0..n).map(|i| x[i + 1] - x[i]).collect(); + + // Solve for interior second derivatives M[1..n-1] (natural BCs: M[0] = M[n] = 0). + // System size m = n-1; for n=1 (two points) the system is empty and the spline + // degenerates to linear interpolation. + let m = n - 1; + let mut m_inner = vec![D::Elem::zero(); m]; + + if m > 0 { + // RHS for equation k (unknown = M[k+1]): + // 6 * ((y[k+2]-y[k+1])/h[k+1] - (y[k+1]-y[k])/h[k]) + let rhs: Vec = (0..m) + .map(|k| six * ((y[k + 2] - y[k + 1]) / h[k + 1] - (y[k + 1] - y[k]) / h[k])) + .collect(); + + // Thomas algorithm (tridiagonal solve). + // Sub-diagonal: h[k], main: 2*(h[k]+h[k+1]), super: h[k+1]. + let mut cp = vec![D::Elem::zero(); m]; // modified super-diagonal + let mut dp = vec![D::Elem::zero(); m]; // modified RHS + + let w0 = two * (h[0] + h[1]); + cp[0] = if m > 1 { h[1] / w0 } else { D::Elem::zero() }; + dp[0] = rhs[0] / w0; + + for k in 1..m { + let w = two * (h[k] + h[k + 1]) - h[k] * cp[k - 1]; + cp[k] = if k < m - 1 { + h[k + 1] / w + } else { + D::Elem::zero() + }; + dp[k] = (rhs[k] - h[k] * dp[k - 1]) / w; + } + + m_inner[m - 1] = dp[m - 1]; + for k in (0..m - 1).rev() { + m_inner[k] = dp[k] - cp[k] * m_inner[k + 1]; + } + } + + // M[i]: M[0]=0, M[1..n-1]=m_inner[0..m-1], M[n]=0. + // For interval i, M[i] = m_inner[i-1] (0 at boundaries). + let m_at = |i: usize| -> D::Elem { + if i == 0 || i == n { + D::Elem::zero() + } else { + m_inner[i - 1] + } + }; + + // S_i(x) = y[i] + b[i]*dx + c[i]*dx^2 + d[i]*dx^3 where dx = x - x[i] + self.b = (0..n) + .map(|i| (y[i + 1] - y[i]) / h[i] - h[i] * (two * m_at(i) + m_at(i + 1)) / six) + .collect(); + self.c = (0..n).map(|i| m_at(i) / two).collect(); + self.d = (0..n) + .map(|i| (m_at(i + 1) - m_at(i)) / (six * h[i])) + .collect(); + + Ok(()) + } + + fn interpolate( + &self, + data: &InterpData1D, + point: &[D::Elem; 1], + ) -> Result { + let grid = data.grid[0].view(); + // For extrapolation, clamp to the boundary interval. + let i = if &point[0] < grid.first().unwrap() { + 0 + } else if &point[0] > grid.last().unwrap() { + grid.len() - 2 + } else { + find_nearest_index(grid, &point[0]) + }; + let dx = point[0] - grid[i]; + Ok(data.values[i] + self.b[i] * dx + self.c[i] * dx * dx + self.d[i] * dx * dx * dx) + } + + /// Returns `true`: the boundary cubic polynomials extend naturally. + fn allow_extrapolate(&self) -> bool { + true + } +} diff --git a/src/interpolator/one/tests.rs b/src/interpolator/one/tests.rs index 4ba2c76..bf2fa14 100644 --- a/src/interpolator/one/tests.rs +++ b/src/interpolator/one/tests.rs @@ -1,5 +1,59 @@ use super::*; +#[test] +fn test_cubic_spline() { + // Linear data: any spline reproduces it exactly (all second derivatives = 0) + let interp = Interp1D::new( + array![0., 1., 2., 3.], + array![1., 3., 5., 7.], // f(x) = 2x + 1 + strategy::CubicSpline::new(), + Extrapolate::Enable, + ) + .unwrap(); + // Knot values + let x = interp.data.grid[0].clone(); + for (i, xi) in x.iter().enumerate() { + assert_approx_eq!(interp.interpolate(&[*xi]).unwrap(), interp.data.values[i]); + } + // Midpoints + assert_approx_eq!(interp.interpolate(&[0.5]).unwrap(), 2.0); + assert_approx_eq!(interp.interpolate(&[1.5]).unwrap(), 4.0); + assert_approx_eq!(interp.interpolate(&[2.5]).unwrap(), 6.0); + // Extrapolation via boundary polynomials + assert_approx_eq!(interp.interpolate(&[-1.0]).unwrap(), -1.0); + assert_approx_eq!(interp.interpolate(&[4.0]).unwrap(), 9.0); +} + +#[test] +fn test_cubic_spline_knot_exactness() { + // Values at all knots must be reproduced exactly regardless of data shape + let interp = Interp1D::new( + array![0., 1., 2., 3., 4.], + array![0., 1., 4., 9., 16.], // f(x) = x^2 + strategy::CubicSpline::new(), + Extrapolate::Error, + ) + .unwrap(); + let x = interp.data.grid[0].clone(); + for (i, xi) in x.iter().enumerate() { + assert_approx_eq!(interp.interpolate(&[*xi]).unwrap(), interp.data.values[i]); + } +} + +#[test] +fn test_cubic_spline_two_points() { + // Degenerate case: 2 points → degenerates to linear interpolation + let interp = Interp1D::new( + array![0., 1.], + array![0., 2.], + strategy::CubicSpline::new(), + Extrapolate::Enable, + ) + .unwrap(); + assert_approx_eq!(interp.interpolate(&[0.5]).unwrap(), 1.0); + assert_approx_eq!(interp.interpolate(&[2.0]).unwrap(), 4.0); // extrapolation +} + #[test] fn test_invalid_args() { let interp = Interp1D::new( diff --git a/src/interpolator/three/strategies.rs b/src/interpolator/three/strategies.rs index f1a9c97..1587080 100644 --- a/src/interpolator/three/strategies.rs +++ b/src/interpolator/three/strategies.rs @@ -241,3 +241,27 @@ where false } } + +impl Strategy3D for CubicSpline +where + D: Data + RawDataClone + Clone, + D::Elem: Float + Debug, +{ + fn interpolate( + &self, + data: &InterpData3D, + point: &[D::Elem; 3], + ) -> Result { + let grids: Vec> = data.grid.iter().map(|g| g.view()).collect(); + Ok(spline_eval_nd_recursive( + &grids, + data.values.view().into_dyn(), + point, + )) + } + + /// Returns `true`: the boundary cubic polynomials extend naturally. + fn allow_extrapolate(&self) -> bool { + true + } +} diff --git a/src/interpolator/three/tests.rs b/src/interpolator/three/tests.rs index e4f5b74..16c70b0 100644 --- a/src/interpolator/three/tests.rs +++ b/src/interpolator/three/tests.rs @@ -1,5 +1,60 @@ use super::*; +#[test] +fn test_cubic_spline() { + // f(x, y, z) = x + 2y + 3z: linear, reproduced exactly by any spline + let interp = Interp3D::new( + array![0., 1., 2.], + array![0., 1., 2.], + array![0., 1., 2.], + array![ + [[0., 3., 6.], [2., 5., 8.], [4., 7., 10.]], + [[1., 4., 7.], [3., 6., 9.], [5., 8., 11.]], + [[2., 5., 8.], [4., 7., 10.], [6., 9., 12.]], + ], + strategy::CubicSpline::new(), + Extrapolate::Enable, + ) + .unwrap(); + // Knots + assert_approx_eq!(interp.interpolate(&[1., 1., 1.]).unwrap(), 6.); + // Midpoints + assert_approx_eq!(interp.interpolate(&[0.5, 0.5, 0.5]).unwrap(), 3.); + assert_approx_eq!(interp.interpolate(&[1.0, 0.5, 1.0]).unwrap(), 5.); + // Extrapolation + assert_approx_eq!(interp.interpolate(&[3., 1., 1.]).unwrap(), 8.); +} + +#[test] +fn test_cubic_spline_knot_exactness() { + let interp = Interp3D::new( + array![0., 1., 2.], + array![0., 1., 2.], + array![0., 1., 2.], + array![ + [[0., 1., 2.], [3., 4., 5.], [6., 7., 8.]], + [[9., 10., 11.], [12., 13., 14.], [15., 16., 17.]], + [[18., 19., 20.], [21., 22., 23.], [24., 25., 26.]], + ], + strategy::CubicSpline::new(), + Extrapolate::Error, + ) + .unwrap(); + let x = interp.data.grid[0].clone(); + let y = interp.data.grid[1].clone(); + let z = interp.data.grid[2].clone(); + for (i, xi) in x.iter().enumerate() { + for (j, yj) in y.iter().enumerate() { + for (k, zk) in z.iter().enumerate() { + assert_approx_eq!( + interp.interpolate(&[*xi, *yj, *zk]).unwrap(), + interp.data.values[[i, j, k]] + ); + } + } + } +} + #[test] fn test_linear() { let interp = Interp3D::new( diff --git a/src/interpolator/two/strategies.rs b/src/interpolator/two/strategies.rs index bbecc58..9838575 100644 --- a/src/interpolator/two/strategies.rs +++ b/src/interpolator/two/strategies.rs @@ -172,3 +172,27 @@ where false } } + +impl Strategy2D for CubicSpline +where + D: Data + RawDataClone + Clone, + D::Elem: Float + Debug, +{ + fn interpolate( + &self, + data: &InterpData2D, + point: &[D::Elem; 2], + ) -> Result { + let grids: Vec> = data.grid.iter().map(|g| g.view()).collect(); + Ok(spline_eval_nd_recursive( + &grids, + data.values.view().into_dyn(), + point, + )) + } + + /// Returns `true`: the boundary cubic polynomials extend naturally. + fn allow_extrapolate(&self) -> bool { + true + } +} diff --git a/src/interpolator/two/tests.rs b/src/interpolator/two/tests.rs index 594a4b9..3882cd6 100644 --- a/src/interpolator/two/tests.rs +++ b/src/interpolator/two/tests.rs @@ -1,5 +1,54 @@ use super::*; +#[test] +fn test_cubic_spline() { + // f(x, y) = 2x + y: linear in both dims, reproduced exactly by any spline + let interp = Interp2D::new( + array![0., 1., 2.], + array![0., 1., 2.], + array![[0., 1., 2.], [2., 3., 4.], [4., 5., 6.]], + strategy::CubicSpline::new(), + Extrapolate::Enable, + ) + .unwrap(); + // Knots + assert_approx_eq!(interp.interpolate(&[1., 1.]).unwrap(), 3.); + assert_approx_eq!(interp.interpolate(&[2., 0.]).unwrap(), 4.); + // Midpoints + assert_approx_eq!(interp.interpolate(&[0.5, 0.5]).unwrap(), 1.5); + assert_approx_eq!(interp.interpolate(&[1.5, 1.0]).unwrap(), 4.); + // Extrapolation + assert_approx_eq!(interp.interpolate(&[3., 1.]).unwrap(), 7.); + assert_approx_eq!(interp.interpolate(&[1., 3.]).unwrap(), 5.); +} + +#[test] +fn test_cubic_spline_knot_exactness() { + let interp = Interp2D::new( + array![0., 1., 2., 3.], + array![0., 1., 2., 3.], + array![ + [0., 1., 4., 9.], + [1., 2., 5., 10.], + [4., 5., 8., 13.], + [9., 10., 13., 18.], + ], // f(x, y) = x^2 + y + strategy::CubicSpline::new(), + Extrapolate::Error, + ) + .unwrap(); + let x = interp.data.grid[0].clone(); + let y = interp.data.grid[1].clone(); + for (i, xi) in x.iter().enumerate() { + for (j, yj) in y.iter().enumerate() { + assert_approx_eq!( + interp.interpolate(&[*xi, *yj]).unwrap(), + interp.data.values[[i, j]] + ); + } + } +} + #[test] fn test_linear() { let interp = Interp2D::new( diff --git a/src/strategy/mod.rs b/src/strategy/mod.rs index 63f4d2f..8ef2979 100644 --- a/src/strategy/mod.rs +++ b/src/strategy/mod.rs @@ -161,6 +161,58 @@ mod step_serde { } } +/// Natural cubic spline interpolation (). +/// +/// Constructs a C² piecewise cubic polynomial through all data points using +/// natural boundary conditions (zero second derivative at the endpoints). +/// Coefficients are precomputed in [`Strategy1D::init`], called automatically +/// by [`Interp1D::new`] and [`Interp1D::set_strategy`]. +/// +/// Supports [`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: Interp1DOwned = Interp1D::new( +/// array![0., 1., 2., 3.], +/// array![1., 3., 5., 7.], +/// strategy::CubicSpline::new(), +/// 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))] +pub struct CubicSpline { + pub(crate) b: Vec, + pub(crate) c: Vec, + pub(crate) d: Vec, +} + +impl Default for CubicSpline { + fn default() -> Self { + Self { + b: Vec::new(), + c: Vec::new(), + d: Vec::new(), + } + } +} + +impl CubicSpline { + /// Create a new cubic spline strategy with no precomputed coefficients. + /// Coefficients are computed automatically when passed to [`Interp1D::new`]. + pub fn new() -> Self { + Self::default() + } +} + #[cfg(test)] mod tests { #[allow(unused_imports)] diff --git a/src/strategy/traits.rs b/src/strategy/traits.rs index b339163..328b60b 100644 --- a/src/strategy/traits.rs +++ b/src/strategy/traits.rs @@ -298,6 +298,100 @@ where clone_trait_object!( StrategyND); +/// Solve for natural cubic spline coefficients and evaluate at `point` in one pass. +/// +/// Implements the Thomas algorithm for the tridiagonal second-derivative system, +/// then evaluates the resulting cubic polynomial. Used by the 2-D, 3-D, and N-D +/// [`CubicSpline`](crate::strategy::CubicSpline) strategy implementations, which +/// cannot precompute coefficients across all slices. +/// +/// `x` and `y` must each have length ≥ 2 (guaranteed by [`InterpData`](crate::data::InterpData) validation). +pub(crate) fn spline_eval_1d(x: ArrayView1, y: ArrayView1, point: T) -> T { + 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 m = n - 1; + let mut m_inner = vec![T::zero(); m]; + + if m > 0 { + let rhs: Vec = (0..m) + .map(|k| six * ((y[k + 2] - y[k + 1]) / h[k + 1] - (y[k + 1] - y[k]) / h[k])) + .collect(); + + let mut cp = vec![T::zero(); m]; + let mut dp = vec![T::zero(); m]; + + let w0 = two * (h[0] + h[1]); + cp[0] = if m > 1 { h[1] / w0 } else { T::zero() }; + dp[0] = rhs[0] / w0; + + for k in 1..m { + let w = two * (h[k] + h[k + 1]) - h[k] * cp[k - 1]; + cp[k] = if k < m - 1 { h[k + 1] / w } else { T::zero() }; + dp[k] = (rhs[k] - h[k] * dp[k - 1]) / w; + } + + m_inner[m - 1] = dp[m - 1]; + for k in (0..m - 1).rev() { + m_inner[k] = dp[k] - cp[k] * m_inner[k + 1]; + } + } + + let m_at = |i: usize| -> T { + if i == 0 || i == n { + T::zero() + } else { + m_inner[i - 1] + } + }; + + let i = if point < *x.first().unwrap() { + 0 + } else if point > *x.last().unwrap() { + n - 1 + } else { + find_nearest_index(x, &point) + }; + + let dx = point - x[i]; + let b = (y[i + 1] - y[i]) / h[i] - h[i] * (two * m_at(i) + m_at(i + 1)) / six; + let c = m_at(i) / two; + let d = (m_at(i + 1) - m_at(i)) / (six * h[i]); + y[i] + b * dx + c * dx * dx + d * dx * dx * dx +} + +/// Recursively evaluates an N-D natural cubic spline via sequential 1-D slicing. +/// +/// For each grid point along the first axis, recursively evaluates the remaining +/// (N-1)-D spline, producing a 1-D array of values. A final [`spline_eval_1d`] call +/// interpolates those values along the first axis at `point[0]`. +/// +/// `grids[k].len() == values.len_of(Axis(k))` must hold at every recursion level, +/// which is guaranteed by [`InterpData`](crate::data::InterpData) validation. +pub(crate) fn spline_eval_nd_recursive( + grids: &[ArrayView1], + values: ArrayViewD, + point: &[T], +) -> T { + debug_assert_eq!(grids.len(), point.len()); + + if grids.len() == 1 { + let y = values + .into_dimensionality::() + .expect("internal: 1-D base case reached with non-1-D values array"); + return spline_eval_1d(grids[0], y, point[0]); + } + + let n = grids[0].len(); + let g: Vec = (0..n) + .map(|i| spline_eval_nd_recursive(&grids[1..], values.index_axis(Axis(0), i), &point[1..])) + .collect(); + + spline_eval_1d(grids[0], ArrayView1::from(&g), point[0]) +} + impl StrategyND for Box> where D: Data + RawDataClone + Clone, From 01984a3639560589ed997cda3e491bb9e54d28aa Mon Sep 17 00:00:00 2001 From: Kyle Carow Date: Mon, 3 Aug 2026 11:27:32 -0600 Subject: [PATCH 2/2] incremental progress --- src/interpolator/n/strategies.rs | 2 +- src/interpolator/n/tests.rs | 4 +- src/interpolator/one/strategies.rs | 87 ++-------- src/interpolator/one/tests.rs | 9 +- src/interpolator/three/strategies.rs | 5 +- src/interpolator/three/tests.rs | 4 +- src/interpolator/two/strategies.rs | 5 +- src/interpolator/two/tests.rs | 4 +- src/lib.rs | 2 + src/strategy/mod.rs | 92 +++++++++-- src/strategy/spline.rs | 236 +++++++++++++++++++++++++++ src/strategy/traits.rs | 216 ------------------------ src/strategy/utils.rs | 125 ++++++++++++++ 13 files changed, 469 insertions(+), 322 deletions(-) create mode 100644 src/strategy/spline.rs create mode 100644 src/strategy/utils.rs diff --git a/src/interpolator/n/strategies.rs b/src/interpolator/n/strategies.rs index 8dcf512..0eac0d8 100644 --- a/src/interpolator/n/strategies.rs +++ b/src/interpolator/n/strategies.rs @@ -237,7 +237,7 @@ where point: &[D::Elem], ) -> Result { let grids: Vec> = data.grid.iter().map(|g| g.view()).collect(); - Ok(spline_eval_nd_recursive(&grids, data.values.view(), point)) + spline_eval_nd_recursive(&grids, data.values.view(), point, &self.boundary_conditions) } /// Returns `true`: the boundary cubic polynomials extend naturally. diff --git a/src/interpolator/n/tests.rs b/src/interpolator/n/tests.rs index d5e3db1..7809524 100644 --- a/src/interpolator/n/tests.rs +++ b/src/interpolator/n/tests.rs @@ -6,7 +6,7 @@ fn test_cubic_spline() { let interp = InterpND::new( vec![array![0., 1., 2.], array![0., 1., 2.]], array![[0., 1., 2.], [2., 3., 4.], [4., 5., 6.]].into_dyn(), - strategy::CubicSpline::new(), + strategy::CubicSpline::natural(), Extrapolate::Enable, ) .unwrap(); @@ -26,7 +26,7 @@ fn test_cubic_spline_knot_exactness() { [9., 10., 13., 18.], ] .into_dyn(), - strategy::CubicSpline::new(), + strategy::CubicSpline::natural(), Extrapolate::Error, ) .unwrap(); diff --git a/src/interpolator/one/strategies.rs b/src/interpolator/one/strategies.rs index 42bd87b..1632a29 100644 --- a/src/interpolator/one/strategies.rs +++ b/src/interpolator/one/strategies.rs @@ -120,78 +120,10 @@ where D: Data + RawDataClone + Clone, D::Elem: Float + Debug, { - /// Computes and caches the per-interval spline coefficients. - /// - /// Uses a natural cubic spline (zero second derivative at endpoints). - /// Solves the resulting tridiagonal system via the Thomas algorithm. + /// Computes and caches `M[0..=n]` using [`compute_m`] for the configured BC. fn init(&mut self, data: &InterpData1D) -> Result<(), ValidateError> { - let x = data.grid[0].view(); - let y = data.values.view(); - let n = x.len() - 1; // number of intervals - - let two = D::Elem::one() + D::Elem::one(); - let six = two + two + two; - - // h[i] = x[i+1] - x[i] - let h: Vec = (0..n).map(|i| x[i + 1] - x[i]).collect(); - - // Solve for interior second derivatives M[1..n-1] (natural BCs: M[0] = M[n] = 0). - // System size m = n-1; for n=1 (two points) the system is empty and the spline - // degenerates to linear interpolation. - let m = n - 1; - let mut m_inner = vec![D::Elem::zero(); m]; - - if m > 0 { - // RHS for equation k (unknown = M[k+1]): - // 6 * ((y[k+2]-y[k+1])/h[k+1] - (y[k+1]-y[k])/h[k]) - let rhs: Vec = (0..m) - .map(|k| six * ((y[k + 2] - y[k + 1]) / h[k + 1] - (y[k + 1] - y[k]) / h[k])) - .collect(); - - // Thomas algorithm (tridiagonal solve). - // Sub-diagonal: h[k], main: 2*(h[k]+h[k+1]), super: h[k+1]. - let mut cp = vec![D::Elem::zero(); m]; // modified super-diagonal - let mut dp = vec![D::Elem::zero(); m]; // modified RHS - - let w0 = two * (h[0] + h[1]); - cp[0] = if m > 1 { h[1] / w0 } else { D::Elem::zero() }; - dp[0] = rhs[0] / w0; - - for k in 1..m { - let w = two * (h[k] + h[k + 1]) - h[k] * cp[k - 1]; - cp[k] = if k < m - 1 { - h[k + 1] / w - } else { - D::Elem::zero() - }; - dp[k] = (rhs[k] - h[k] * dp[k - 1]) / w; - } - - m_inner[m - 1] = dp[m - 1]; - for k in (0..m - 1).rev() { - m_inner[k] = dp[k] - cp[k] * m_inner[k + 1]; - } - } - - // M[i]: M[0]=0, M[1..n-1]=m_inner[0..m-1], M[n]=0. - // For interval i, M[i] = m_inner[i-1] (0 at boundaries). - let m_at = |i: usize| -> D::Elem { - if i == 0 || i == n { - D::Elem::zero() - } else { - m_inner[i - 1] - } - }; - - // S_i(x) = y[i] + b[i]*dx + c[i]*dx^2 + d[i]*dx^3 where dx = x - x[i] - self.b = (0..n) - .map(|i| (y[i + 1] - y[i]) / h[i] - h[i] * (two * m_at(i) + m_at(i + 1)) / six) - .collect(); - self.c = (0..n).map(|i| m_at(i) / two).collect(); - self.d = (0..n) - .map(|i| (m_at(i + 1) - m_at(i)) / (six * h[i])) - .collect(); - + let new_m = compute_m(data.grid[0].view(), data.values.view(), self.bc_for_dim(0))?; + self.m = new_m; Ok(()) } @@ -209,8 +141,17 @@ where } else { find_nearest_index(grid, &point[0]) }; - let dx = point[0] - grid[i]; - Ok(data.values[i] + self.b[i] * dx + self.c[i] * dx * dx + self.d[i] * dx * dx * dx) + let two = D::Elem::one() + D::Elem::one(); + let six = two + two + two; + let h = grid[i + 1] - grid[i]; + let dx = point[0] - grid[i]; // t - x[i] + let dx_r = h - dx; // x[i+1] - t + let six_h = six * h; + let h2_over_six = h * h / six; + Ok(self.m[i] * dx_r * dx_r * dx_r / six_h + + self.m[i + 1] * dx * dx * dx / six_h + + (data.values[i] - self.m[i] * h2_over_six) * dx_r / h + + (data.values[i + 1] - self.m[i + 1] * h2_over_six) * dx / h) } /// Returns `true`: the boundary cubic polynomials extend naturally. diff --git a/src/interpolator/one/tests.rs b/src/interpolator/one/tests.rs index bf2fa14..5d2ef50 100644 --- a/src/interpolator/one/tests.rs +++ b/src/interpolator/one/tests.rs @@ -6,7 +6,7 @@ fn test_cubic_spline() { let interp = Interp1D::new( array![0., 1., 2., 3.], array![1., 3., 5., 7.], // f(x) = 2x + 1 - strategy::CubicSpline::new(), + strategy::CubicSpline::not_a_knot(), Extrapolate::Enable, ) .unwrap(); @@ -30,7 +30,7 @@ fn test_cubic_spline_knot_exactness() { let interp = Interp1D::new( array![0., 1., 2., 3., 4.], array![0., 1., 4., 9., 16.], // f(x) = x^2 - strategy::CubicSpline::new(), + strategy::CubicSpline::not_a_knot(), Extrapolate::Error, ) .unwrap(); @@ -42,11 +42,12 @@ fn test_cubic_spline_knot_exactness() { #[test] fn test_cubic_spline_two_points() { - // Degenerate case: 2 points → degenerates to linear interpolation + // Degenerate case: 2 points → degenerates to linear interpolation. + // Uses Natural BC since NotAKnot requires ≥ 3 points. let interp = Interp1D::new( array![0., 1.], array![0., 2.], - strategy::CubicSpline::new(), + strategy::CubicSpline::natural(), Extrapolate::Enable, ) .unwrap(); diff --git a/src/interpolator/three/strategies.rs b/src/interpolator/three/strategies.rs index 1587080..a86df88 100644 --- a/src/interpolator/three/strategies.rs +++ b/src/interpolator/three/strategies.rs @@ -253,11 +253,12 @@ where point: &[D::Elem; 3], ) -> Result { let grids: Vec> = data.grid.iter().map(|g| g.view()).collect(); - Ok(spline_eval_nd_recursive( + spline_eval_nd_recursive( &grids, data.values.view().into_dyn(), point, - )) + &self.boundary_conditions, + ) } /// Returns `true`: the boundary cubic polynomials extend naturally. diff --git a/src/interpolator/three/tests.rs b/src/interpolator/three/tests.rs index 16c70b0..76af1cd 100644 --- a/src/interpolator/three/tests.rs +++ b/src/interpolator/three/tests.rs @@ -12,7 +12,7 @@ fn test_cubic_spline() { [[1., 4., 7.], [3., 6., 9.], [5., 8., 11.]], [[2., 5., 8.], [4., 7., 10.], [6., 9., 12.]], ], - strategy::CubicSpline::new(), + strategy::CubicSpline::natural(), Extrapolate::Enable, ) .unwrap(); @@ -36,7 +36,7 @@ fn test_cubic_spline_knot_exactness() { [[9., 10., 11.], [12., 13., 14.], [15., 16., 17.]], [[18., 19., 20.], [21., 22., 23.], [24., 25., 26.]], ], - strategy::CubicSpline::new(), + strategy::CubicSpline::natural(), Extrapolate::Error, ) .unwrap(); diff --git a/src/interpolator/two/strategies.rs b/src/interpolator/two/strategies.rs index 9838575..1984b1d 100644 --- a/src/interpolator/two/strategies.rs +++ b/src/interpolator/two/strategies.rs @@ -184,11 +184,12 @@ where point: &[D::Elem; 2], ) -> Result { let grids: Vec> = data.grid.iter().map(|g| g.view()).collect(); - Ok(spline_eval_nd_recursive( + spline_eval_nd_recursive( &grids, data.values.view().into_dyn(), point, - )) + &self.boundary_conditions, + ) } /// Returns `true`: the boundary cubic polynomials extend naturally. diff --git a/src/interpolator/two/tests.rs b/src/interpolator/two/tests.rs index 3882cd6..cac838d 100644 --- a/src/interpolator/two/tests.rs +++ b/src/interpolator/two/tests.rs @@ -7,7 +7,7 @@ fn test_cubic_spline() { array![0., 1., 2.], array![0., 1., 2.], array![[0., 1., 2.], [2., 3., 4.], [4., 5., 6.]], - strategy::CubicSpline::new(), + strategy::CubicSpline::natural(), Extrapolate::Enable, ) .unwrap(); @@ -33,7 +33,7 @@ fn test_cubic_spline_knot_exactness() { [4., 5., 8., 13.], [9., 10., 13., 18.], ], // f(x, y) = x^2 + y - strategy::CubicSpline::new(), + strategy::CubicSpline::not_a_knot(), Extrapolate::Error, ) .unwrap(); diff --git a/src/lib.rs b/src/lib.rs index bd78834..b8512fc 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -42,7 +42,9 @@ pub use interpolator::data; pub(crate) use interpolator::data::*; pub(crate) use error::*; +pub(crate) use strategy::spline::*; pub(crate) use strategy::traits::*; +pub(crate) use strategy::utils::*; pub(crate) use std::fmt::Debug; diff --git a/src/strategy/mod.rs b/src/strategy/mod.rs index 8ef2979..4594953 100644 --- a/src/strategy/mod.rs +++ b/src/strategy/mod.rs @@ -3,7 +3,9 @@ use super::*; pub mod enums; +pub(crate) mod spline; pub mod traits; +pub(crate) mod utils; /// Linear interpolation: #[derive(Debug, Clone, PartialEq)] @@ -161,10 +163,10 @@ mod step_serde { } } -/// Natural cubic spline interpolation (). +/// Cubic spline interpolation (). /// -/// Constructs a C² piecewise cubic polynomial through all data points using -/// natural boundary conditions (zero second derivative at the endpoints). +/// Constructs a C² piecewise cubic polynomial through all data points. +/// The boundary condition is set by [`boundary_conditions`](CubicSpline::boundary_conditions). /// Coefficients are precomputed in [`Strategy1D::init`], called automatically /// by [`Interp1D::new`] and [`Interp1D::set_strategy`]. /// @@ -180,7 +182,7 @@ mod step_serde { /// let interp: Interp1DOwned = Interp1D::new( /// array![0., 1., 2., 3.], /// array![1., 3., 5., 7.], -/// strategy::CubicSpline::new(), +/// strategy::CubicSpline::not_a_knot(), /// Extrapolate::Enable, /// ) /// .unwrap(); @@ -190,26 +192,80 @@ mod step_serde { #[derive(Debug, Clone, PartialEq)] #[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] pub struct CubicSpline { - pub(crate) b: Vec, - pub(crate) c: Vec, - pub(crate) d: Vec, + /// Boundary conditions, one per dimension or a single entry broadcast to all. + pub boundary_conditions: Vec>, + /// Second derivatives `M[i] = S''(x_i)` at each grid point, length `n + 1` + /// for `n` intervals. Populated by [`Strategy1D::init`]; boundary values + /// are determined by [`boundary_conditions`](CubicBC). + /// + /// Not included in the serialized form. Call [`Interpolator::validate`] after + /// deserializing a 1-D interpolator to recompute these coefficients before use. + #[cfg_attr(feature = "serde", serde(skip))] + pub(crate) m: Vec, } -impl Default for CubicSpline { - fn default() -> Self { +/// Boundary conditions for [`CubicSpline`]. +#[derive(Debug, Clone, PartialEq)] +#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] +pub enum CubicBC { + /// C³ continuity at the second and penultimate knots; no extra input required. + /// Generally gives better accuracy than [`Natural`](CubicBC::Natural) + /// for smooth functions. + NotAKnot, + /// Zero second derivative at both endpoints. + Natural, + /// Specified first derivative at both endpoints. + Clamped { + /// First derivative at the left (lower) endpoint. + left: T, + /// First derivative at the right (upper) endpoint. + right: T, + }, + /// First and second derivatives match at both endpoints. + /// Requires `values[0] == values[n]`. + Periodic, +} + +impl CubicSpline { + /// Returns the boundary condition for the given dimension. + /// A single-entry vec is broadcast to all dimensions. + pub(crate) fn bc_for_dim(&self, dim: usize) -> &CubicBC { + let bcs = &self.boundary_conditions; + &bcs[if bcs.len() == 1 { 0 } else { dim }] + } + + /// 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 { - b: Vec::new(), - c: Vec::new(), - d: Vec::new(), + boundary_conditions: vec![CubicBC::NotAKnot], + m: Vec::new(), } } -} -impl CubicSpline { - /// Create a new cubic spline strategy with no precomputed coefficients. - /// Coefficients are computed automatically when passed to [`Interp1D::new`]. - pub fn new() -> Self { - Self::default() + /// Create a cubic spline with natural (zero second derivative at endpoints) BCs. + pub fn natural() -> Self { + Self { + boundary_conditions: vec![CubicBC::Natural], + m: Vec::new(), + } + } + + /// Create a cubic spline with specified first derivatives at both endpoints. + pub fn clamped(left: T, right: T) -> Self { + Self { + boundary_conditions: vec![CubicBC::Clamped { left, right }], + m: Vec::new(), + } + } + + /// Create a cubic spline with periodic boundary conditions. + /// Requires `values[0] == values[n]`. + pub fn periodic() -> Self { + Self { + boundary_conditions: vec![CubicBC::Periodic], + m: Vec::new(), + } } } diff --git a/src/strategy/spline.rs b/src/strategy/spline.rs new file mode 100644 index 0000000..8fba0bf --- /dev/null +++ b/src/strategy/spline.rs @@ -0,0 +1,236 @@ +//! Cubic spline algorithms shared across all dimensionalities. + +use super::*; + +/// 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 [`Strategy1D::init`] (stores in `self.m`) and [`spline_eval_1d`] (on the fly). +pub(crate) fn compute_m( + x: ArrayView1, + y: ArrayView1, + bc: &CubicBC, +) -> Result, ValidateError> { + 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(); + + Ok(match bc { + CubicBC::NotAKnot => { + if n < 3 { + return Err(ValidateError::Other( + "CubicSpline with NotAKnot requires at least 4 data points".into(), + )); + } + let mut sub: Vec = h[1..n - 2].to_vec(); + sub.push(h[n - 2] * h[n - 2] - h[n - 1] * h[n - 1]); + let mut sup = vec![h[1] * h[1] - h[0] * h[0]]; + sup.extend_from_slice(&h[2..n - 1]); + let mut diag = vec![(h[0] + h[1]) * (h[0] + two * h[1])]; + for k in 1..n - 2 { + diag.push(two * (h[k] + h[k + 1])); + } + diag.push((h[n - 2] + h[n - 1]) * (two * h[n - 2] + h[n - 1])); + let mut rhs = vec![h[1] * u[0]]; + for k in 1..n - 2 { + rhs.push(u[k]); + } + rhs.push(h[n - 2] * u[n - 2]); + let inner = thomas(&sub, &diag, &sup, &rhs); + let m0 = ((h[0] + h[1]) * inner[0] - h[0] * inner[1]) / h[1]; + let mn = ((h[n - 2] + h[n - 1]) * inner[n - 2] - h[n - 1] * inner[n - 3]) / h[n - 2]; + let mut m = vec![m0]; + m.extend(inner); + m.push(mn); + m + } + CubicBC::Natural => { + let mut sub = h[..n.saturating_sub(1)].to_vec(); + sub.push(T::zero()); + let mut diag = vec![T::one()]; + for k in 0..n.saturating_sub(1) { + diag.push(two * (h[k] + h[k + 1])); + } + diag.push(T::one()); + let mut sup = vec![T::zero()]; + sup.extend_from_slice(&h[1..n]); + let mut rhs = vec![T::zero()]; + rhs.extend_from_slice(&u); + rhs.push(T::zero()); + thomas(&sub, &diag, &sup, &rhs) + } + CubicBC::Clamped { left, right } => { + let (l, r) = (*left, *right); + let mut diag = vec![two * h[0]]; + for k in 0..n.saturating_sub(1) { + diag.push(two * (h[k] + h[k + 1])); + } + diag.push(two * h[n - 1]); + let mut rhs = vec![six * (slopes[0] - l)]; + rhs.extend_from_slice(&u); + rhs.push(six * (r - slopes[n - 1])); + thomas(&h, &diag, &h, &rhs) + } + CubicBC::Periodic => { + if y[0] != y[n] { + return Err(ValidateError::Other( + "CubicSpline with Periodic BC requires values[0] == values[n]".into(), + )); + } + 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 eval_spline_from_m( + x: ArrayView1, + y: ArrayView1, + m: &[T], + point: T, +) -> T { + let n = x.len() - 1; + let two = T::one() + T::one(); + let six = two + two + two; + let i = if point < *x.first().unwrap() { + 0 + } else if point > *x.last().unwrap() { + n - 1 + } else { + find_nearest_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 +} + +/// Computes and evaluates a 1-D cubic spline at `point` in one pass, respecting `bc`. +/// +/// Used by the sequential N-D path in [`spline_eval_nd_recursive`]. +pub(crate) fn spline_eval_1d( + x: ArrayView1, + y: ArrayView1, + point: T, + bc: &CubicBC, +) -> Result { + let m = compute_m(x, y, bc).map_err(|e| InterpolateError::Other(e.to_string()))?; + Ok(eval_spline_from_m(x, y, &m, point)) +} + +/// Recursively evaluates an N-D cubic spline via sequential 1-D slicing, respecting `bcs`. +/// +/// `bcs` length 1 broadcasts to all dimensions; length N applies per-dimension. +pub(crate) fn spline_eval_nd_recursive( + grids: &[ArrayView1], + values: ArrayViewD, + point: &[T], + bcs: &[CubicBC], +) -> Result { + debug_assert_eq!(grids.len(), point.len()); + debug_assert!(!bcs.is_empty()); + + let current_bc = &bcs[0]; + let next_bcs = if bcs.len() > 1 { &bcs[1..] } else { bcs }; + + if grids.len() == 1 { + let y = values + .into_dimensionality::() + .map_err(|_| InterpolateError::Other( + "internal: non-1-D values at 1-D base case, grids.len() == values.ndim() invariant broken".into() + ))?; + return spline_eval_1d(grids[0], y, point[0], current_bc); + } + + let n = grids[0].len(); + let g: Vec = (0..n) + .map(|i| { + spline_eval_nd_recursive( + &grids[1..], + values.index_axis(Axis(0), i), + &point[1..], + next_bcs, + ) + }) + .collect::, _>>()?; + + spline_eval_1d(grids[0], ArrayView1::from(&g), point[0], current_bc) +} diff --git a/src/strategy/traits.rs b/src/strategy/traits.rs index 328b60b..248a2a7 100644 --- a/src/strategy/traits.rs +++ b/src/strategy/traits.rs @@ -2,128 +2,6 @@ use super::*; -/// Find nearest index in `arr` left of `target` -/// -/// This method contains code from RouteE Compass, another open-source NLR-developed tool -/// -/// -pub fn find_nearest_index(arr: ArrayView1, target: &T) -> usize { - if target == arr.last().unwrap() { - return arr.len() - 2; - } - - let mut low = 0; - let mut high = arr.len() - 1; - - while low < high { - let mid = low + (high - low) / 2; - - if &arr[mid] >= target { - high = mid; - } else { - low = mid + 1; - } - } - - if low > 0 && &arr[low] >= target { - low - 1 - } else { - low - } -} - -/// Returns the step index for `point` in `grid` using the given [`StepDirection`]. -/// -/// Handles all exact grid-point edge cases that arise from [`find_nearest_index`]'s -/// interval semantics (returning the lower bracket rather than the exact position). -pub(crate) fn step_index( - dir: StepDirection, - grid: ArrayView1, - point: &T, -) -> usize { - match dir { - StepDirection::Lower => { - let x_l = find_nearest_index(grid, point); - // find_nearest_index returns i where grid[i] < point <= grid[i+1] for interior - // matches, so an exact match at grid[i+1] gives i instead of i+1. Correct both: - if point == grid.last().unwrap() { - grid.len() - 1 - } else if *point == grid[x_l + 1] { - x_l + 1 - } else { - x_l - } - } - StepDirection::Upper => { - // find_nearest_index returns 0 when point == grid[0], giving x_l+1 = 1 - // which would skip values[0]. Handle the first-element case explicitly: - if point == grid.first().unwrap() { - 0 - } else { - find_nearest_index(grid, point) + 1 - } - } - } -} - -/// Returns the exact grid index if `point` lies on `grid[lower]` or `grid[lower+1]`, else `None`. -/// -/// Used to short-circuit interpolation when a query point coincides with a grid coordinate. -pub(crate) fn exact_index( - grid: ArrayView1, - lower: usize, - point: &T, -) -> Option { - if grid[lower] == *point { - Some(lower) - } else if grid[lower + 1] == *point { - Some(lower + 1) - } else { - None - } -} - -/// Computes the lower bracket index for a uniformly-spaced grid in O(1). -/// -/// Equivalent to [`find_nearest_index`] but replaces binary search with direct arithmetic. -/// Only valid when the grid spacing is uniform — validate with [`check_uniform_grid`] first. -pub(crate) fn uniform_lower_index(grid0: T, step: T, n: usize, point: T) -> usize { - let t = (point - grid0) / step; - if t < T::zero() { - 0 - } else { - t.floor().to_usize().unwrap_or(0).min(n - 2) - } -} - -/// Validates that `grid` is uniformly spaced within floating-point tolerance. -/// -/// Uses a relative tolerance of 1024 × ε to accommodate accumulated floating-point rounding -/// error in grids constructed from repeated arithmetic. -pub(crate) fn check_uniform_grid( - grid: ArrayView1, - dim: usize, -) -> Result<(), ValidateError> { - let step = grid[1] - grid[0]; - // 1024 * epsilon via 10 doublings — avoids numeric literal casting - let tolerance = { - let mut tol = T::epsilon(); - for _ in 0..10 { - tol = tol + tol; - } - step.abs() * tol - }; - for i in 1..grid.len() - 1 { - let gap = grid[i + 1] - grid[i]; - if (gap - step).abs() > tolerance { - return Err(ValidateError::Other(format!( - "LinearUniform: grid[{dim}] is not uniformly spaced (gap at index {i})" - ))); - } - } - Ok(()) -} - /// 1-D interpolation strategy. pub trait Strategy1D: Debug + DynClone where @@ -298,100 +176,6 @@ where clone_trait_object!( StrategyND); -/// Solve for natural cubic spline coefficients and evaluate at `point` in one pass. -/// -/// Implements the Thomas algorithm for the tridiagonal second-derivative system, -/// then evaluates the resulting cubic polynomial. Used by the 2-D, 3-D, and N-D -/// [`CubicSpline`](crate::strategy::CubicSpline) strategy implementations, which -/// cannot precompute coefficients across all slices. -/// -/// `x` and `y` must each have length ≥ 2 (guaranteed by [`InterpData`](crate::data::InterpData) validation). -pub(crate) fn spline_eval_1d(x: ArrayView1, y: ArrayView1, point: T) -> T { - 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 m = n - 1; - let mut m_inner = vec![T::zero(); m]; - - if m > 0 { - let rhs: Vec = (0..m) - .map(|k| six * ((y[k + 2] - y[k + 1]) / h[k + 1] - (y[k + 1] - y[k]) / h[k])) - .collect(); - - let mut cp = vec![T::zero(); m]; - let mut dp = vec![T::zero(); m]; - - let w0 = two * (h[0] + h[1]); - cp[0] = if m > 1 { h[1] / w0 } else { T::zero() }; - dp[0] = rhs[0] / w0; - - for k in 1..m { - let w = two * (h[k] + h[k + 1]) - h[k] * cp[k - 1]; - cp[k] = if k < m - 1 { h[k + 1] / w } else { T::zero() }; - dp[k] = (rhs[k] - h[k] * dp[k - 1]) / w; - } - - m_inner[m - 1] = dp[m - 1]; - for k in (0..m - 1).rev() { - m_inner[k] = dp[k] - cp[k] * m_inner[k + 1]; - } - } - - let m_at = |i: usize| -> T { - if i == 0 || i == n { - T::zero() - } else { - m_inner[i - 1] - } - }; - - let i = if point < *x.first().unwrap() { - 0 - } else if point > *x.last().unwrap() { - n - 1 - } else { - find_nearest_index(x, &point) - }; - - let dx = point - x[i]; - let b = (y[i + 1] - y[i]) / h[i] - h[i] * (two * m_at(i) + m_at(i + 1)) / six; - let c = m_at(i) / two; - let d = (m_at(i + 1) - m_at(i)) / (six * h[i]); - y[i] + b * dx + c * dx * dx + d * dx * dx * dx -} - -/// Recursively evaluates an N-D natural cubic spline via sequential 1-D slicing. -/// -/// For each grid point along the first axis, recursively evaluates the remaining -/// (N-1)-D spline, producing a 1-D array of values. A final [`spline_eval_1d`] call -/// interpolates those values along the first axis at `point[0]`. -/// -/// `grids[k].len() == values.len_of(Axis(k))` must hold at every recursion level, -/// which is guaranteed by [`InterpData`](crate::data::InterpData) validation. -pub(crate) fn spline_eval_nd_recursive( - grids: &[ArrayView1], - values: ArrayViewD, - point: &[T], -) -> T { - debug_assert_eq!(grids.len(), point.len()); - - if grids.len() == 1 { - let y = values - .into_dimensionality::() - .expect("internal: 1-D base case reached with non-1-D values array"); - return spline_eval_1d(grids[0], y, point[0]); - } - - let n = grids[0].len(); - let g: Vec = (0..n) - .map(|i| spline_eval_nd_recursive(&grids[1..], values.index_axis(Axis(0), i), &point[1..])) - .collect(); - - spline_eval_1d(grids[0], ArrayView1::from(&g), point[0]) -} - impl StrategyND for Box> where D: Data + RawDataClone + Clone, diff --git a/src/strategy/utils.rs b/src/strategy/utils.rs new file mode 100644 index 0000000..817a062 --- /dev/null +++ b/src/strategy/utils.rs @@ -0,0 +1,125 @@ +//! Shared index and grid utilities for interpolation strategies. + +use super::*; + +/// Find nearest index in `arr` left of `target` +/// +/// This method contains code from RouteE Compass, another open-source NLR-developed tool +/// +/// +pub fn find_nearest_index(arr: ArrayView1, target: &T) -> usize { + if target == arr.last().unwrap() { + return arr.len() - 2; + } + + let mut low = 0; + let mut high = arr.len() - 1; + + while low < high { + let mid = low + (high - low) / 2; + + if &arr[mid] >= target { + high = mid; + } else { + low = mid + 1; + } + } + + if low > 0 && &arr[low] >= target { + low - 1 + } else { + low + } +} + +/// Returns the step index for `point` in `grid` using the given [`StepDirection`]. +/// +/// Handles all exact grid-point edge cases that arise from [`find_nearest_index`]'s +/// interval semantics (returning the lower bracket rather than the exact position). +pub(crate) fn step_index( + dir: StepDirection, + grid: ArrayView1, + point: &T, +) -> usize { + match dir { + StepDirection::Lower => { + let x_l = find_nearest_index(grid, point); + // find_nearest_index returns i where grid[i] < point <= grid[i+1] for interior + // matches, so an exact match at grid[i+1] gives i instead of i+1. Correct both: + if point == grid.last().unwrap() { + grid.len() - 1 + } else if *point == grid[x_l + 1] { + x_l + 1 + } else { + x_l + } + } + StepDirection::Upper => { + // find_nearest_index returns 0 when point == grid[0], giving x_l+1 = 1 + // which would skip values[0]. Handle the first-element case explicitly: + if point == grid.first().unwrap() { + 0 + } else { + find_nearest_index(grid, point) + 1 + } + } + } +} + +/// Returns the exact grid index if `point` lies on `grid[lower]` or `grid[lower+1]`, else `None`. +/// +/// Used to short-circuit interpolation when a query point coincides with a grid coordinate. +pub(crate) fn exact_index( + grid: ArrayView1, + lower: usize, + point: &T, +) -> Option { + if grid[lower] == *point { + Some(lower) + } else if grid[lower + 1] == *point { + Some(lower + 1) + } else { + None + } +} + +/// Computes the lower bracket index for a uniformly-spaced grid in O(1). +/// +/// Equivalent to [`find_nearest_index`] but replaces binary search with direct arithmetic. +/// Only valid when the grid spacing is uniform — validate with [`check_uniform_grid`] first. +pub(crate) fn uniform_lower_index(grid0: T, step: T, n: usize, point: T) -> usize { + let t = (point - grid0) / step; + if t < T::zero() { + 0 + } else { + t.floor().to_usize().unwrap_or(0).min(n - 2) + } +} + +/// Validates that `grid` is uniformly spaced within floating-point tolerance. +/// +/// Uses a relative tolerance of 1024 × ε to accommodate accumulated floating-point rounding +/// error in grids constructed from repeated arithmetic. +pub(crate) fn check_uniform_grid( + grid: ArrayView1, + dim: usize, +) -> Result<(), ValidateError> { + let step = grid[1] - grid[0]; + // 1024 * epsilon via 10 doublings — avoids numeric literal casting + let tolerance = { + let mut tol = T::epsilon(); + for _ in 0..10 { + tol = tol + tol; + } + step.abs() * tol + }; + for i in 1..grid.len() - 1 { + let gap = grid[i + 1] - grid[i]; + if (gap - step).abs() > tolerance { + return Err(ValidateError::Other(format!( + "LinearUniform: grid[{dim}] is not uniformly spaced (gap at index {i})" + ))); + } + } + Ok(()) +}