numra_interp/
lib.rs

1#![allow(clippy::needless_range_loop)]
2#![allow(clippy::manual_range_contains)]
3
4//! Interpolation and splines for Numra.
5//!
6//! This crate provides 1D interpolation methods:
7//!
8//! - **Linear** ([`Linear`]): Piecewise linear interpolation
9//! - **Cubic spline** ([`CubicSpline`]): Natural, clamped, and not-a-knot variants
10//! - **PCHIP** ([`Pchip`]): Monotonicity-preserving piecewise cubic Hermite
11//! - **Akima** ([`Akima`]): Robust piecewise cubic (tolerant of outliers)
12//! - **Barycentric Lagrange** ([`BarycentricLagrange`]): Stable polynomial interpolation
13//!
14//! Author: Moussa Leblouba
15//! Date: 9 February 2026
16//! Modified: 2 May 2026
17
18pub mod akima;
19pub mod cubic_spline;
20pub mod error;
21pub mod hermite;
22pub mod linear;
23pub mod polynomial;
24
25pub use akima::Akima;
26pub use cubic_spline::CubicSpline;
27pub use error::InterpError;
28pub use hermite::Pchip;
29pub use linear::Linear;
30pub use polynomial::BarycentricLagrange;
31
32use numra_core::Scalar;
33
34/// Trait for 1D interpolants.
35pub trait Interpolant<S: Scalar> {
36    /// Evaluate the interpolant at `x`.
37    /// Extrapolates using the first/last segment if `x` is outside the data range.
38    fn interpolate(&self, x: S) -> S;
39
40    /// Evaluate at multiple points.
41    fn interpolate_vec(&self, xs: &[S]) -> Vec<S> {
42        xs.iter().map(|&x| self.interpolate(x)).collect()
43    }
44
45    /// First derivative at `x`, if available.
46    fn derivative(&self, x: S) -> Option<S> {
47        let _ = x;
48        None
49    }
50
51    /// Definite integral from `a` to `b`, if available.
52    fn integrate(&self, a: S, b: S) -> Option<S> {
53        let _ = (a, b);
54        None
55    }
56}
57
58/// Interpolation method selector for [`interp1d`].
59pub enum Interp1dMethod {
60    Linear,
61    CubicSpline,
62    Pchip,
63    Akima,
64}
65
66/// Convenience: create an interpolant from data using the specified method.
67pub fn interp1d<S: Scalar>(
68    x: &[S],
69    y: &[S],
70    method: Interp1dMethod,
71) -> Result<Box<dyn Interpolant<S>>, InterpError> {
72    match method {
73        Interp1dMethod::Linear => Ok(Box::new(Linear::new(x, y)?)),
74        Interp1dMethod::CubicSpline => Ok(Box::new(CubicSpline::natural(x, y)?)),
75        Interp1dMethod::Pchip => Ok(Box::new(Pchip::new(x, y)?)),
76        Interp1dMethod::Akima => Ok(Box::new(Akima::new(x, y)?)),
77    }
78}
79
80// ============================================================================
81// Internal utilities
82// ============================================================================
83
84/// Binary search for the interval containing `x`.
85/// Returns `i` such that `xs[i] <= x < xs[i+1]`, clamped to valid range.
86pub(crate) fn search_sorted<S: Scalar>(xs: &[S], x: S) -> usize {
87    let n = xs.len();
88    if n < 2 {
89        return 0;
90    }
91    if x.to_f64() <= xs[0].to_f64() {
92        return 0;
93    }
94    if x.to_f64() >= xs[n - 1].to_f64() {
95        return n - 2;
96    }
97    let mut lo = 0;
98    let mut hi = n - 1;
99    while lo < hi - 1 {
100        let mid = (lo + hi) / 2;
101        if x < xs[mid] {
102            hi = mid;
103        } else {
104            lo = mid;
105        }
106    }
107    lo
108}
109
110/// Validate interpolation data: matching lengths, sufficient points, sorted.
111pub(crate) fn validate_data<S: Scalar>(
112    x: &[S],
113    y: &[S],
114    min_points: usize,
115) -> Result<(), InterpError> {
116    if x.len() != y.len() {
117        return Err(InterpError::DimensionMismatch {
118            x_len: x.len(),
119            y_len: y.len(),
120        });
121    }
122    if x.len() < min_points {
123        return Err(InterpError::TooFewPoints {
124            got: x.len(),
125            min: min_points,
126        });
127    }
128    for i in 1..x.len() {
129        if x[i].to_f64() <= x[i - 1].to_f64() {
130            return Err(InterpError::UnsortedData);
131        }
132    }
133    Ok(())
134}
135
136/// Evaluate piecewise cubic: a[i] + b[i]*t + c[i]*t^2 + d[i]*t^3 where t = x - x_knots[i].
137pub(crate) fn eval_piecewise_cubic<S: Scalar>(
138    x_knots: &[S],
139    a: &[S],
140    b: &[S],
141    c: &[S],
142    d: &[S],
143    x: S,
144) -> S {
145    let i = search_sorted(x_knots, x);
146    let t = x - x_knots[i];
147    a[i] + t * (b[i] + t * (c[i] + t * d[i]))
148}
149
150/// Derivative of piecewise cubic: b[i] + 2*c[i]*t + 3*d[i]*t^2.
151pub(crate) fn eval_piecewise_cubic_deriv<S: Scalar>(
152    x_knots: &[S],
153    b: &[S],
154    c: &[S],
155    d: &[S],
156    x: S,
157) -> S {
158    let i = search_sorted(x_knots, x);
159    let t = x - x_knots[i];
160    b[i] + t * (S::TWO * c[i] + S::from_f64(3.0) * t * d[i])
161}
162
163/// Integrate piecewise cubic from `lo` to `hi`.
164pub(crate) fn integrate_piecewise_cubic<S: Scalar>(
165    x_knots: &[S],
166    a: &[S],
167    b: &[S],
168    c: &[S],
169    d: &[S],
170    lo: S,
171    hi: S,
172) -> S {
173    if hi.to_f64() <= lo.to_f64() {
174        return S::ZERO;
175    }
176    let i_lo = search_sorted(x_knots, lo);
177    let i_hi = search_sorted(x_knots, hi);
178
179    let antideriv = |idx: usize, t: S| -> S {
180        let t2 = t * t;
181        let t3 = t2 * t;
182        let t4 = t3 * t;
183        a[idx] * t
184            + b[idx] * t2 * S::HALF
185            + c[idx] * t3 / S::from_f64(3.0)
186            + d[idx] * t4 / S::from_f64(4.0)
187    };
188
189    let mut result = S::ZERO;
190    for i in i_lo..=i_hi {
191        let t_lo = if i == i_lo { lo - x_knots[i] } else { S::ZERO };
192        let t_hi = if i == i_hi {
193            hi - x_knots[i]
194        } else {
195            x_knots[i + 1] - x_knots[i]
196        };
197        result += antideriv(i, t_hi) - antideriv(i, t_lo);
198    }
199    result
200}
numra_interp/lib.rs

numra_interp/
lib.rs
