quantsupport 0.1.5

Rust library for 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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

use crate::ad::{dual::DualFwd, scalar::Scalar};

/// Numeric behavior required for bilinear interpolation.
pub trait BilinearValue: Scalar {
    /// 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 DualFwd {
    fn lerp(a: Self, b: Self, t: f64) -> Self {
        (a + (b - a) * t).into()
    }
}

/// Input point for bilinear interpolation.
#[derive(Clone)]
pub struct BilinearPoint<T: Scalar> {
    /// 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;
        }

        // Queries landing exactly on the last pillar clamp to the last cell.
        let clamp_upper = |grid: &[f64], v: f64, idx: usize| {
            if idx == grid.len() && (v - grid[grid.len() - 1]) <= f64::EPSILON * v.abs().max(1.0) {
                idx - 1
            } else {
                idx
            }
        };
        let ix = clamp_upper(&xs, x, xs.partition_point(|v| *v <= x));
        let iy = clamp_upper(&ys, y, 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::dual::DualFwd,
        math::interpolation::bilinear::{BilinearInterpolator, BilinearPoint},
    };

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