rustyqlib/core/interpolation/
pchip.rs

1//! PCHIP: Fritsch-Carlson monotone piecewise cubic Hermite
2//! interpolation. Shape-preserving — the interpolant is monotone
3//! wherever the data are, so it never overshoots or invents wiggles the
4//! way an unconstrained cubic spline can. The safe choice for curves
5//! that must stay monotone (discount factors, CDFs) or non-negative.
6use crate::core::errors::RustyQLibError;
7
8/// A monotonicity-preserving cubic Hermite interpolant.
9#[derive(Debug, Clone)]
10pub struct Pchip {
11    xs: Vec<f64>,
12    ys: Vec<f64>,
13    /// Knot slopes after the Fritsch-Carlson limiter.
14    d: Vec<f64>,
15}
16
17impl Pchip {
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("pchip", "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("pchip", "knots must be strictly increasing"));
25        }
26        let h: Vec<f64> = xs.windows(2).map(|w| w[1] - w[0]).collect();
27        let delta: Vec<f64> =
28            (0..n - 1).map(|i| (ys[i + 1] - ys[i]) / h[i]).collect();
29
30        let mut d = vec![0.0; n];
31        if n == 2 {
32            d[0] = delta[0];
33            d[1] = delta[0];
34        } else {
35            // interior: weighted harmonic mean when the secants agree in
36            // sign, zero otherwise (this is what preserves monotonicity)
37            for i in 1..n - 1 {
38                if delta[i - 1] * delta[i] > 0.0 {
39                    let w1 = 2.0 * h[i] + h[i - 1];
40                    let w2 = h[i] + 2.0 * h[i - 1];
41                    d[i] = (w1 + w2) / (w1 / delta[i - 1] + w2 / delta[i]);
42                }
43            }
44            d[0] = end_slope(h[0], h[1], delta[0], delta[1]);
45            d[n - 1] = end_slope(h[n - 2], h[n - 3], delta[n - 2], delta[n - 3]);
46        }
47        Ok(Pchip { xs: xs.to_vec(), ys: ys.to_vec(), d })
48    }
49
50    /// Interpolant value at `x` (linear extrapolation with the end slope).
51    pub fn eval(&self, x: f64) -> f64 {
52        let n = self.xs.len();
53        if x <= self.xs[0] {
54            return self.ys[0] + self.d[0] * (x - self.xs[0]);
55        }
56        if x >= self.xs[n - 1] {
57            return self.ys[n - 1] + self.d[n - 1] * (x - self.xs[n - 1]);
58        }
59        let i = self.xs[1..n - 1].partition_point(|&xi| xi < x);
60        hermite(
61            x,
62            self.xs[i],
63            self.xs[i + 1],
64            self.ys[i],
65            self.ys[i + 1],
66            self.d[i],
67            self.d[i + 1],
68        )
69    }
70}
71
72/// Three-point end slope with the Fritsch-Carlson clips.
73fn end_slope(h0: f64, h1: f64, delta0: f64, delta1: f64) -> f64 {
74    let d = ((2.0 * h0 + h1) * delta0 - h0 * delta1) / (h0 + h1);
75    if d * delta0 <= 0.0 {
76        0.0
77    } else if delta0 * delta1 < 0.0 && d.abs() > 3.0 * delta0.abs() {
78        3.0 * delta0
79    } else {
80        d
81    }
82}
83
84/// Cubic Hermite basis evaluation on `[x0, x1]`.
85pub(crate) fn hermite(x: f64, x0: f64, x1: f64, y0: f64, y1: f64, d0: f64, d1: f64) -> f64 {
86    let h = x1 - x0;
87    let t = (x - x0) / h;
88    let t2 = t * t;
89    let t3 = t2 * t;
90    y0 * (2.0 * t3 - 3.0 * t2 + 1.0)
91        + y1 * (-2.0 * t3 + 3.0 * t2)
92        + d0 * h * (t3 - 2.0 * t2 + t)
93        + d1 * h * (t3 - t2)
94}
95
96#[cfg(test)]
97mod tests {
98    use super::super::cubic_spline::{BoundaryCondition, CubicSpline};
99    use super::*;
100
101    /// A monotone step-like data set that makes free splines overshoot.
102    fn step_data() -> (Vec<f64>, Vec<f64>) {
103        (
104            vec![0.0, 1.0, 2.0, 3.0, 4.0, 5.0],
105            vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0],
106        )
107    }
108
109    #[test]
110    fn interpolates_the_knots() {
111        let (xs, ys) = step_data();
112        let p = Pchip::new(&xs, &ys).unwrap();
113        for (x, y) in xs.iter().zip(&ys) {
114            assert!((p.eval(*x) - y).abs() < 1e-14);
115        }
116    }
117
118    #[test]
119    fn preserves_monotonicity_where_a_cubic_spline_overshoots() {
120        let (xs, ys) = step_data();
121        let pchip = Pchip::new(&xs, &ys).unwrap();
122        let spline = CubicSpline::new(&xs, &ys, BoundaryCondition::Natural).unwrap();
123
124        let mut prev = f64::NEG_INFINITY;
125        let mut spline_overshoots = false;
126        for i in 0..=500 {
127            let x = i as f64 * 0.01;
128            let v = pchip.eval(x);
129            // monotone and inside the data range
130            assert!(v >= prev - 1e-12, "pchip not monotone at {x}");
131            assert!((-1e-12..=1.0 + 1e-12).contains(&v), "pchip overshoots at {x}");
132            prev = v;
133            let s = spline.eval(x);
134            if !(0.0..=1.0).contains(&s) {
135                spline_overshoots = true;
136            }
137        }
138        // and the comparison is meaningful: the free spline DOES overshoot
139        assert!(spline_overshoots, "cubic spline unexpectedly shape-preserving");
140    }
141
142    #[test]
143    fn flat_data_stays_exactly_flat() {
144        let xs = [0.0, 1.0, 2.0, 3.0];
145        let ys = [5.0, 5.0, 5.0, 5.0];
146        let p = Pchip::new(&xs, &ys).unwrap();
147        for i in 0..=30 {
148            assert!((p.eval(i as f64 * 0.1) - 5.0).abs() < 1e-12);
149        }
150    }
151}
rustyqlib/core/interpolation/pchip.rs

rustyqlib/core/interpolation/
pchip.rs
