numra_interp/
akima.rs

1//! Akima interpolation.
2//!
3//! A piecewise cubic method that is robust to outliers. Uses a weighted
4//! average of adjacent secants to determine slopes.
5//!
6//! Author: Moussa Leblouba
7//! Date: 9 February 2026
8//! Modified: 2 May 2026
9
10use numra_core::Scalar;
11
12use crate::cubic_spline::coefficients_from_slopes;
13use crate::error::InterpError;
14use crate::{eval_piecewise_cubic, eval_piecewise_cubic_deriv, integrate_piecewise_cubic};
15use crate::{validate_data, Interpolant};
16
17/// Akima interpolant.
18///
19/// Uses Akima's method for computing slopes, which is robust to outliers
20/// in the data. Requires at least 5 data points.
21pub struct Akima<S: Scalar> {
22    x: Vec<S>,
23    a: Vec<S>,
24    b: Vec<S>,
25    c: Vec<S>,
26    d: Vec<S>,
27}
28
29impl<S: Scalar> Akima<S> {
30    /// Create an Akima interpolant from data points.
31    ///
32    /// # Errors
33    ///
34    /// Returns error if fewer than 5 points or data is not sorted.
35    pub fn new(x: &[S], y: &[S]) -> Result<Self, InterpError> {
36        validate_data(x, y, 5)?;
37        let n = x.len();
38        let slopes = compute_akima_slopes(x, y);
39        let (a, b, c, d) = coefficients_from_slopes(x, y, &slopes);
40
41        Ok(Self {
42            x: x[..n].to_vec(),
43            a,
44            b,
45            c,
46            d,
47        })
48    }
49}
50
51impl<S: Scalar> Interpolant<S> for Akima<S> {
52    fn interpolate(&self, x: S) -> S {
53        eval_piecewise_cubic(&self.x, &self.a, &self.b, &self.c, &self.d, x)
54    }
55
56    fn derivative(&self, x: S) -> Option<S> {
57        Some(eval_piecewise_cubic_deriv(
58            &self.x, &self.b, &self.c, &self.d, x,
59        ))
60    }
61
62    fn integrate(&self, a: S, b: S) -> Option<S> {
63        Some(integrate_piecewise_cubic(
64            &self.x, &self.a, &self.b, &self.c, &self.d, a, b,
65        ))
66    }
67}
68
69/// Compute Akima slopes.
70fn compute_akima_slopes<S: Scalar>(x: &[S], y: &[S]) -> Vec<S> {
71    let n = x.len();
72    let n_seg = n - 1;
73
74    // Compute secants s[0..n_seg-1]
75    let s: Vec<S> = (0..n_seg)
76        .map(|i| (y[i + 1] - y[i]) / (x[i + 1] - x[i]))
77        .collect();
78
79    // Extend secant array with 2 values on each side (Akima's extrapolation)
80    // s_ext[-2], s_ext[-1], s_ext[0], ..., s_ext[n_seg-1], s_ext[n_seg], s_ext[n_seg+1]
81    let sm2 = S::from_f64(3.0) * s[0] - S::TWO * s[1];
82    let sm1 = S::TWO * s[0] - s[1];
83    let sp1 = S::TWO * s[n_seg - 1] - s[n_seg - 2];
84    let sp2 = S::from_f64(3.0) * s[n_seg - 1] - S::TWO * s[n_seg - 2];
85
86    // Build extended secant array indexed 0..n_seg+3 corresponding to original -2..n_seg+1
87    let mut s_ext = Vec::with_capacity(n_seg + 4);
88    s_ext.push(sm2);
89    s_ext.push(sm1);
90    for &si in &s {
91        s_ext.push(si);
92    }
93    s_ext.push(sp1);
94    s_ext.push(sp2);
95    // Now s_ext[k] corresponds to original index k-2
96
97    let mut slopes = Vec::with_capacity(n);
98    for i in 0..n {
99        // Map to s_ext indices: s_{i-2}..s_{i+1} → s_ext[i]..s_ext[i+3]
100        let w1 = (s_ext[i + 3] - s_ext[i + 2]).abs();
101        let w2 = (s_ext[i + 1] - s_ext[i]).abs();
102        let denom = w1 + w2;
103        if denom > S::EPSILON * S::from_f64(100.0) {
104            slopes.push((w1 * s_ext[i + 1] + w2 * s_ext[i + 2]) / denom);
105        } else {
106            slopes.push((s_ext[i + 1] + s_ext[i + 2]) * S::HALF);
107        }
108    }
109
110    slopes
111}
112
113#[cfg(test)]
114mod tests {
115    use super::*;
116    use approx::assert_relative_eq;
117
118    #[test]
119    fn test_akima_at_knots() {
120        let x = vec![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
121        let y = vec![0.0, 1.0, 0.5, 1.5, 1.0, 2.0];
122        let a = Akima::new(&x, &y).unwrap();
123        for (&xi, &yi) in x.iter().zip(y.iter()) {
124            assert_relative_eq!(a.interpolate(xi), yi, epsilon = 1e-12);
125        }
126    }
127
128    #[test]
129    fn test_akima_smooth() {
130        // Akima on sin(x)
131        let n = 20;
132        let x: Vec<f64> = (0..n)
133            .map(|i| i as f64 * core::f64::consts::PI * 2.0 / (n - 1) as f64)
134            .collect();
135        let y: Vec<f64> = x.iter().map(|&xi| xi.sin()).collect();
136        let a = Akima::new(&x, &y).unwrap();
137
138        let test_x = 1.0;
139        let err = (a.interpolate(test_x) - test_x.sin()).abs();
140        assert!(err < 1e-3, "Akima error too large: {}", err);
141    }
142
143    #[test]
144    fn test_akima_outlier_robustness() {
145        // Data with an outlier at x=2
146        let x = vec![0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
147        let y = vec![0.0, 0.0, 10.0, 0.0, 0.0, 0.0, 0.0]; // outlier at x=2
148
149        let a = Akima::new(&x, &y).unwrap();
150
151        // Akima should not overshoot excessively away from the outlier
152        let y_at_4_5 = a.interpolate(4.5);
153        assert!(
154            y_at_4_5.abs() < 1.0,
155            "Akima overshoots away from outlier: {}",
156            y_at_4_5
157        );
158    }
159
160    #[test]
161    fn test_akima_linear_data() {
162        // Linear data should give linear interpolation
163        let x = vec![0.0, 1.0, 2.0, 3.0, 4.0];
164        let y = vec![0.0, 1.0, 2.0, 3.0, 4.0];
165        let a = Akima::new(&x, &y).unwrap();
166        assert_relative_eq!(a.interpolate(1.5), 1.5, epsilon = 1e-12);
167        assert_relative_eq!(a.interpolate(2.7), 2.7, epsilon = 1e-12);
168    }
169
170    #[test]
171    fn test_akima_derivative() {
172        let x = vec![0.0, 1.0, 2.0, 3.0, 4.0];
173        let y = vec![0.0, 1.0, 2.0, 3.0, 4.0];
174        let a = Akima::new(&x, &y).unwrap();
175        // Derivative of linear function should be 1
176        assert_relative_eq!(a.derivative(2.0).unwrap(), 1.0, epsilon = 1e-10);
177    }
178
179    #[test]
180    fn test_akima_errors() {
181        assert!(Akima::<f64>::new(&[0.0, 1.0, 2.0, 3.0], &[0.0, 1.0, 2.0, 3.0]).is_err());
182    }
183}
numra_interp/akima.rs

numra_interp/
akima.rs
