rustyqlib/core/interpolation/
cubic_spline.rs

1//! C2 cubic spline interpolation with the three classical boundary
2//! conditions.
3//!
4//! The spline is built on the knots' second derivatives `M_i`, solved
5//! from the tridiagonal continuity equations plus two boundary rows:
6//!
7//! - **Natural**: `M = 0` at both ends (zero curvature — the default
8//!   when nothing is known about the ends, slightly flattens there);
9//! - **Clamped**: first derivatives prescribed at both ends (use when
10//!   end slopes are known, e.g. from a model);
11//! - **Not-a-Knot**: third-derivative continuity across the first and
12//!   last interior knots — the best all-round accuracy without extra
13//!   information; reproduces a single cubic polynomial exactly.
14//!
15//! Evaluation outside the knot range extrapolates linearly with the end
16//! slope (curvature is not continued — safer for financial data).
17
18use crate::core::optimization::numerics::solve_dense;
19use crate::core::errors::RustyQLibError;
20
21/// Boundary condition for [`CubicSpline`].
22#[derive(Debug, Clone, Copy, PartialEq)]
23pub enum BoundaryCondition {
24    Natural,
25    /// Prescribed end slopes `(start, end)`.
26    Clamped { start_slope: f64, end_slope: f64 },
27    NotAKnot,
28}
29
30/// A cubic spline through `(xs, ys)` knots.
31#[derive(Debug, Clone)]
32pub struct CubicSpline {
33    xs: Vec<f64>,
34    ys: Vec<f64>,
35    /// Second derivatives at the knots.
36    m: Vec<f64>,
37}
38
39impl CubicSpline {
40    /// Build a spline through the knots (`xs` strictly increasing,
41    /// at least two points; Not-a-Knot needs at least four and falls
42    /// back to Natural below that).
43    pub fn new(xs: &[f64], ys: &[f64], bc: BoundaryCondition) -> Result<Self, RustyQLibError> {
44        let n = xs.len();
45        if n < 2 || ys.len() != n {
46            return Err(RustyQLibError::invalid_input("cubic_spline", "need at least two knots with matching y values"));
47        }
48        if xs.windows(2).any(|w| w[1] <= w[0]) {
49            return Err(RustyQLibError::invalid_input("cubic_spline", "knots must be strictly increasing"));
50        }
51        if n == 2 {
52            // a segment: linear, zero curvature
53            return Ok(CubicSpline { xs: xs.to_vec(), ys: ys.to_vec(), m: vec![0.0; 2] });
54        }
55        let bc = match bc {
56            BoundaryCondition::NotAKnot if n < 4 => BoundaryCondition::Natural,
57            other => other,
58        };
59
60        // dense (n x n) system in the second derivatives M; knot counts
61        // are small (pillars), so O(n^3) is irrelevant next to clarity
62        let h: Vec<f64> = xs.windows(2).map(|w| w[1] - w[0]).collect();
63        let mut a = vec![vec![0.0; n]; n];
64        let mut b = vec![0.0; n];
65        for i in 1..n - 1 {
66            a[i][i - 1] = h[i - 1] / 6.0;
67            a[i][i] = (h[i - 1] + h[i]) / 3.0;
68            a[i][i + 1] = h[i] / 6.0;
69            b[i] = (ys[i + 1] - ys[i]) / h[i] - (ys[i] - ys[i - 1]) / h[i - 1];
70        }
71        match bc {
72            BoundaryCondition::Natural => {
73                a[0][0] = 1.0;
74                a[n - 1][n - 1] = 1.0;
75            }
76            BoundaryCondition::Clamped { start_slope, end_slope } => {
77                a[0][0] = h[0] / 3.0;
78                a[0][1] = h[0] / 6.0;
79                b[0] = (ys[1] - ys[0]) / h[0] - start_slope;
80                a[n - 1][n - 2] = h[n - 2] / 6.0;
81                a[n - 1][n - 1] = h[n - 2] / 3.0;
82                b[n - 1] = end_slope - (ys[n - 1] - ys[n - 2]) / h[n - 2];
83            }
84            BoundaryCondition::NotAKnot => {
85                // third-derivative continuity at the second and
86                // second-to-last knots
87                a[0][0] = h[1];
88                a[0][1] = -(h[0] + h[1]);
89                a[0][2] = h[0];
90                a[n - 1][n - 3] = h[n - 2];
91                a[n - 1][n - 2] = -(h[n - 3] + h[n - 2]);
92                a[n - 1][n - 1] = h[n - 3];
93            }
94        }
95        let m = solve_dense(&mut a, &mut b).ok_or(RustyQLibError::invalid_input("cubic_spline", "singular spline system"))?;
96        Ok(CubicSpline { xs: xs.to_vec(), ys: ys.to_vec(), m })
97    }
98
99    fn segment(&self, x: f64) -> usize {
100        let n = self.xs.len();
101        self.xs[1..n - 1].partition_point(|&xi| xi < x)
102    }
103
104    /// Spline value at `x` (linear extrapolation beyond the knots).
105    pub fn eval(&self, x: f64) -> f64 {
106        let n = self.xs.len();
107        if x <= self.xs[0] {
108            return self.ys[0] + self.derivative(self.xs[0]) * (x - self.xs[0]);
109        }
110        if x >= self.xs[n - 1] {
111            return self.ys[n - 1] + self.derivative(self.xs[n - 1]) * (x - self.xs[n - 1]);
112        }
113        let i = self.segment(x);
114        let h = self.xs[i + 1] - self.xs[i];
115        let (dl, dr) = (self.xs[i + 1] - x, x - self.xs[i]);
116        self.m[i] * dl * dl * dl / (6.0 * h)
117            + self.m[i + 1] * dr * dr * dr / (6.0 * h)
118            + (self.ys[i] / h - self.m[i] * h / 6.0) * dl
119            + (self.ys[i + 1] / h - self.m[i + 1] * h / 6.0) * dr
120    }
121
122    /// First derivative at `x` (constant beyond the knots).
123    pub fn derivative(&self, x: f64) -> f64 {
124        let n = self.xs.len();
125        let x = x.clamp(self.xs[0], self.xs[n - 1]);
126        let i = self.segment(x).min(n - 2);
127        let h = self.xs[i + 1] - self.xs[i];
128        let (dl, dr) = (self.xs[i + 1] - x, x - self.xs[i]);
129        -self.m[i] * dl * dl / (2.0 * h)
130            + self.m[i + 1] * dr * dr / (2.0 * h)
131            + (self.ys[i + 1] - self.ys[i]) / h
132            - (self.m[i + 1] - self.m[i]) * h / 6.0
133    }
134
135    /// Second derivative at `x` (zero beyond the knots).
136    pub fn second_derivative(&self, x: f64) -> f64 {
137        let n = self.xs.len();
138        if x < self.xs[0] || x > self.xs[n - 1] {
139            return 0.0;
140        }
141        let i = self.segment(x).min(n - 2);
142        let h = self.xs[i + 1] - self.xs[i];
143        self.m[i] * (self.xs[i + 1] - x) / h + self.m[i + 1] * (x - self.xs[i]) / h
144    }
145}
146
147#[cfg(test)]
148mod tests {
149    use super::*;
150
151    fn knots() -> (Vec<f64>, Vec<f64>) {
152        let xs: Vec<f64> = (0..7).map(|i| i as f64 * 0.5).collect();
153        let ys: Vec<f64> = xs.iter().map(|x| x.sin()).collect();
154        (xs, ys)
155    }
156
157    #[test]
158    fn all_boundary_conditions_interpolate_the_knots() {
159        let (xs, ys) = knots();
160        for bc in [
161            BoundaryCondition::Natural,
162            BoundaryCondition::Clamped { start_slope: 1.0, end_slope: 3.0_f64.cos() },
163            BoundaryCondition::NotAKnot,
164        ] {
165            let s = CubicSpline::new(&xs, &ys, bc).unwrap();
166            for (x, y) in xs.iter().zip(&ys) {
167                assert!((s.eval(*x) - y).abs() < 1e-12, "{bc:?} at {x}");
168            }
169        }
170    }
171
172    #[test]
173    fn not_a_knot_reproduces_a_cubic_exactly() {
174        // the defining property: for data from one cubic polynomial the
175        // not-a-knot spline IS that polynomial
176        let p = |x: f64| 2.0 - x + 3.0 * x * x - 0.5 * x * x * x;
177        let xs: Vec<f64> = (0..6).map(|i| i as f64).collect();
178        let ys: Vec<f64> = xs.iter().map(|&x| p(x)).collect();
179        let s = CubicSpline::new(&xs, &ys, BoundaryCondition::NotAKnot).unwrap();
180        for i in 0..=50 {
181            let x = i as f64 * 0.1;
182            assert!((s.eval(x) - p(x)).abs() < 1e-9, "x = {x}");
183        }
184    }
185
186    #[test]
187    fn natural_ends_have_zero_curvature() {
188        let (xs, ys) = knots();
189        let s = CubicSpline::new(&xs, &ys, BoundaryCondition::Natural).unwrap();
190        assert!(s.second_derivative(xs[0]).abs() < 1e-10);
191        assert!(s.second_derivative(*xs.last().unwrap()).abs() < 1e-10);
192    }
193
194    #[test]
195    fn clamped_ends_match_the_prescribed_slopes() {
196        let (xs, ys) = knots();
197        let bc = BoundaryCondition::Clamped { start_slope: 1.0, end_slope: 3.0_f64.cos() };
198        let s = CubicSpline::new(&xs, &ys, bc).unwrap();
199        assert!((s.derivative(0.0) - 1.0).abs() < 1e-10);
200        assert!((s.derivative(3.0) - 3.0_f64.cos()).abs() < 1e-10);
201        // clamping with the true sin slopes beats natural near the ends
202        let natural = CubicSpline::new(&xs, &ys, BoundaryCondition::Natural).unwrap();
203        let x = 0.1;
204        assert!((s.eval(x) - x.sin()).abs() < (natural.eval(x) - x.sin()).abs());
205    }
206
207    #[test]
208    fn rejects_bad_knots() {
209        assert!(CubicSpline::new(&[0.0, 0.0, 1.0], &[1.0, 2.0, 3.0], BoundaryCondition::Natural)
210            .is_err());
211        assert!(CubicSpline::new(&[0.0], &[1.0], BoundaryCondition::Natural).is_err());
212    }
213}
rustyqlib/core/interpolation/cubic_spline.rs

rustyqlib/core/interpolation/
cubic_spline.rs
