quantsupport 0.1.0

Rust library for fixed-income, derivative pricing and risk analytics.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

use crate::ad::adreal::{ADReal, IsReal};

/// Numeric behavior required for bilinear interpolation.
pub trait BilinearValue: IsReal {
    /// Linearly interpolates between `a` and `b` using scalar weight `t`.
    fn lerp(a: Self, b: Self, t: f64) -> Self;
}

impl BilinearValue for f64 {
    fn lerp(a: Self, b: Self, t: f64) -> Self {
        (b - a).mul_add(t, a)
    }
}

impl BilinearValue for ADReal {
    fn lerp(a: Self, b: Self, t: f64) -> Self {
        (a + (b - a) * t).into()
    }
}

/// # `BilinearPoint`
/// Input point for bilinear interpolation.
#[derive(Clone)]
pub struct BilinearPoint<T: IsReal> {
    /// X coordinate.
    pub x: f64,
    /// Y coordinate.
    pub y: f64,
    /// Point value.
    pub value: T,
}

/// # `BilinearInterpolator`
/// Bilinear interpolation implementation independent from volatility containers.
pub struct BilinearInterpolator;

impl BilinearInterpolator {
    /// Interpolate at (`x`,`y`) from rectangular points.
    #[must_use]
    pub fn interpolate<T: BilinearValue>(
        x: f64,
        y: f64,
        mut points: Vec<BilinearPoint<T>>,
    ) -> Option<T> {
        if points.is_empty() {
            return None;
        }

        points.sort_by(|a, b| a.x.total_cmp(&b.x).then_with(|| a.y.total_cmp(&b.y)));

        let mut xs = points.iter().map(|p| p.x).collect::<Vec<_>>();
        xs.sort_by(f64::total_cmp);
        xs.dedup_by(|a, b| (*a - *b).abs() < f64::EPSILON);

        let mut ys = points.iter().map(|p| p.y).collect::<Vec<_>>();
        ys.sort_by(f64::total_cmp);
        ys.dedup_by(|a, b| (*a - *b).abs() < f64::EPSILON);

        if xs.len() < 2 || ys.len() < 2 {
            return None;
        }

        let ix = xs.partition_point(|v| *v <= x);
        let iy = ys.partition_point(|v| *v <= y);
        if ix == 0 || ix >= xs.len() || iy == 0 || iy >= ys.len() {
            return None;
        }

        let (x0, x1) = (xs[ix - 1], xs[ix]);
        let (y0, y1) = (ys[iy - 1], ys[iy]);

        let lookup = |lx: f64, ly: f64| {
            points
                .iter()
                .find(|p| (p.x - lx).abs() < f64::EPSILON && (p.y - ly).abs() < f64::EPSILON)
                .cloned()
        };

        let p00 = lookup(x0, y0)?;
        let p10 = lookup(x1, y0)?;
        let p01 = lookup(x0, y1)?;
        let p11 = lookup(x1, y1)?;

        let tx = if (x1 - x0).abs() < f64::EPSILON {
            0.0
        } else {
            (x - x0) / (x1 - x0)
        };
        let ty = if (y1 - y0).abs() < f64::EPSILON {
            0.0
        } else {
            (y - y0) / (y1 - y0)
        };

        let v0 = T::lerp(p00.value, p10.value, tx);
        let v1 = T::lerp(p01.value, p11.value, tx);
        Some(T::lerp(v0, v1, ty))
    }
}

#[cfg(test)]
mod tests {
    use crate::{
        ad::adreal::{ADReal, IsReal},
        math::interpolation::bilinear::{BilinearInterpolator, BilinearPoint},
    };

    #[test]
    fn bilinear_interpolates_center() {
        let points = vec![
            BilinearPoint { x: 0.0, y: 0.0, value: ADReal::from(1.0) },
            BilinearPoint { x: 1.0, y: 0.0, value: ADReal::from(3.0) },
            BilinearPoint { x: 0.0, y: 1.0, value: ADReal::from(2.0) },
            BilinearPoint { x: 1.0, y: 1.0, value: ADReal::from(4.0) },
        ];
        let out = BilinearInterpolator::interpolate(0.5, 0.5, points).expect("center interpolation");
        assert!((out.value() - 2.5).abs() < 1e-12);
    }
}