rustyqlib/core/interpolation/
akima.rs

1//! Akima's spline: cubic Hermite interpolation with slopes chosen from
2//! local weighted secant differences. Purely local (a moved point only
3//! affects its neighborhood) and far less prone to the wide oscillations
4//! a global cubic spline develops around outliers and flat runs.
5
6use super::pchip::hermite;
7use crate::core::errors::RustyQLibError;
8
9/// An Akima spline interpolant.
10#[derive(Debug, Clone)]
11pub struct Akima {
12    xs: Vec<f64>,
13    ys: Vec<f64>,
14    d: Vec<f64>,
15}
16
17impl Akima {
18    pub fn new(xs: &[f64], ys: &[f64]) -> Result<Self, RustyQLibError> {
19        let n = xs.len();
20        if n < 2 || ys.len() != n {
21            return Err(RustyQLibError::invalid_input("akima", "need at least two knots with matching y values"));
22        }
23        if xs.windows(2).any(|w| w[1] <= w[0]) {
24            return Err(RustyQLibError::invalid_input("akima", "knots must be strictly increasing"));
25        }
26        // secants with Akima's quadratic extension at both ends:
27        // ext[k] corresponds to delta_{k-2} for knot arithmetic below
28        let mut ext = Vec::with_capacity(n + 3);
29        ext.resize(2, 0.0);
30        for i in 0..n - 1 {
31            ext.push((ys[i + 1] - ys[i]) / (xs[i + 1] - xs[i]));
32        }
33        let m = ext.len();
34        ext.push(2.0 * ext[m - 1] - ext[m - 2]);
35        ext.push(2.0 * ext[m] - ext[m - 1]);
36        ext[1] = 2.0 * ext[2] - ext[3];
37        ext[0] = 2.0 * ext[1] - ext[2];
38
39        let d: Vec<f64> = (0..n)
40            .map(|i| {
41                // slopes around knot i: ext[i..i+4] = delta_{i-2..i+1}
42                let w1 = (ext[i + 3] - ext[i + 2]).abs();
43                let w2 = (ext[i + 1] - ext[i]).abs();
44                if w1 + w2 > 1e-300 {
45                    (w1 * ext[i + 1] + w2 * ext[i + 2]) / (w1 + w2)
46                } else {
47                    0.5 * (ext[i + 1] + ext[i + 2])
48                }
49            })
50            .collect();
51        Ok(Akima { xs: xs.to_vec(), ys: ys.to_vec(), d })
52    }
53
54    /// Interpolant value at `x` (linear extrapolation with the end slope).
55    pub fn eval(&self, x: f64) -> f64 {
56        let n = self.xs.len();
57        if x <= self.xs[0] {
58            return self.ys[0] + self.d[0] * (x - self.xs[0]);
59        }
60        if x >= self.xs[n - 1] {
61            return self.ys[n - 1] + self.d[n - 1] * (x - self.xs[n - 1]);
62        }
63        let i = self.xs[1..n - 1].partition_point(|&xi| xi < x);
64        hermite(
65            x,
66            self.xs[i],
67            self.xs[i + 1],
68            self.ys[i],
69            self.ys[i + 1],
70            self.d[i],
71            self.d[i + 1],
72        )
73    }
74}
75
76#[cfg(test)]
77mod tests {
78    use super::super::cubic_spline::{BoundaryCondition, CubicSpline};
79    use super::*;
80
81    #[test]
82    fn interpolates_the_knots() {
83        let xs = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
84        let ys = [0.0, 0.5, 2.0, 1.5, 1.0, 3.0];
85        let a = Akima::new(&xs, &ys).unwrap();
86        for (x, y) in xs.iter().zip(&ys) {
87            assert!((a.eval(*x) - y).abs() < 1e-13);
88        }
89    }
90
91    #[test]
92    fn flat_runs_stay_flat_where_a_spline_rings() {
93        // Akima's classic showcase: a flat run next to a jump — the
94        // global spline rings along the flat section, Akima does not
95        let xs = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
96        let ys = [0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0];
97        let akima = Akima::new(&xs, &ys).unwrap();
98        let spline = CubicSpline::new(&xs, &ys, BoundaryCondition::Natural).unwrap();
99
100        // on the flat left section the Akima interpolant is exactly zero
101        let mut spline_rings = false;
102        for i in 0..=25 {
103            let x = i as f64 * 0.1; // [0, 2.5]
104            assert!(akima.eval(x).abs() < 1e-12, "akima rings at {x}");
105            if spline.eval(x).abs() > 1e-3 {
106                spline_rings = true;
107            }
108        }
109        assert!(spline_rings, "cubic spline unexpectedly local");
110    }
111
112    #[test]
113    fn straight_line_data_reproduces_the_line() {
114        let xs = [0.0, 1.0, 3.0, 6.0];
115        let ys: Vec<f64> = xs.iter().map(|x| 2.0 * x - 1.0).collect();
116        let a = Akima::new(&xs, &ys).unwrap();
117        for i in 0..=70 {
118            let x = i as f64 * 0.1;
119            assert!((a.eval(x) - (2.0 * x - 1.0)).abs() < 1e-12, "x = {x}");
120        }
121    }
122}
rustyqlib/core/interpolation/akima.rs

rustyqlib/core/interpolation/
akima.rs
