rustyqlib/core/interpolation/
linear.rs

1//! Piecewise linear interpolation and the pillar bracketing shared by
2//! curves and smiles.
3
4/// Linear blend `a + w (b - a)`.
5pub fn lerp(a: f64, b: f64, w: f64) -> f64 {
6    a + w * (b - a)
7}
8
9/// Bracket `x` in the sorted pillar grid `xs`: returns `(idx, w)` where
10/// `xs[idx - 1] <= x <= xs[idx]` and `w` is the weight of the upper
11/// pillar. `x` is clamped to the grid, so `idx` is always in
12/// `[1, xs.len() - 1]` and `w` in `[0, 1]`.
13pub fn bracket(xs: &[f64], x: f64) -> (usize, f64) {
14    assert!(xs.len() >= 2, "need at least two pillars");
15    let n = xs.len();
16    if x <= xs[0] {
17        return (1, 0.0);
18    }
19    if x >= xs[n - 1] {
20        return (n - 1, 1.0);
21    }
22    let idx = xs.partition_point(|&xi| xi < x);
23    let (x0, x1) = (xs[idx - 1], xs[idx]);
24    (idx, (x - x0) / (x1 - x0))
25}
26
27/// Piecewise linear interpolation of `(xs, ys)` at `x`, flat beyond the
28/// ends. `xs` must be sorted strictly increasing.
29pub fn linear_interp(xs: &[f64], ys: &[f64], x: f64) -> f64 {
30    assert_eq!(xs.len(), ys.len());
31    let (idx, w) = bracket(xs, x);
32    lerp(ys[idx - 1], ys[idx], w)
33}
34
35/// [`linear_interp`] over `(x, y)` pairs sorted by `x` — the smile
36/// storage shape.
37pub fn interp_pairs(points: &[(f64, f64)], x: f64) -> f64 {
38    assert!(points.len() >= 2, "need at least two points");
39    let n = points.len();
40    if x <= points[0].0 {
41        return points[0].1;
42    }
43    if x >= points[n - 1].0 {
44        return points[n - 1].1;
45    }
46    let idx = points.partition_point(|&(xi, _)| xi < x);
47    let (x0, y0) = points[idx - 1];
48    let (x1, y1) = points[idx];
49    lerp(y0, y1, (x - x0) / (x1 - x0))
50}
51
52#[cfg(test)]
53mod tests {
54    use super::*;
55
56    #[test]
57    fn interpolates_and_extrapolates_flat() {
58        let xs = [1.0, 2.0, 4.0];
59        let ys = [10.0, 20.0, 40.0];
60        assert_eq!(linear_interp(&xs, &ys, 1.5), 15.0);
61        assert_eq!(linear_interp(&xs, &ys, 3.0), 30.0);
62        assert_eq!(linear_interp(&xs, &ys, 0.0), 10.0); // flat left
63        assert_eq!(linear_interp(&xs, &ys, 9.0), 40.0); // flat right
64        assert_eq!(interp_pairs(&[(1.0, 10.0), (2.0, 20.0)], 1.25), 12.5);
65    }
66
67    #[test]
68    fn bracket_clamps_to_the_grid() {
69        let xs = [0.0, 1.0, 3.0];
70        assert_eq!(bracket(&xs, -5.0), (1, 0.0));
71        assert_eq!(bracket(&xs, 5.0), (2, 1.0));
72        let (idx, w) = bracket(&xs, 2.0);
73        assert_eq!(idx, 2);
74        assert!((w - 0.5).abs() < 1e-15);
75    }
76}
rustyqlib/core/interpolation/linear.rs

rustyqlib/core/interpolation/
linear.rs
