nautilus_core/
math.rs

1// -------------------------------------------------------------------------------------------------
2//  Copyright (C) 2015-2025 Posei Systems Pty Ltd. All rights reserved.
3//  https://poseitrader.io
4//
5//  Licensed under the GNU Lesser General Public License Version 3.0 (the "License");
6//  You may not use this file except in compliance with the License.
7//  You may obtain a copy of the License at https://www.gnu.org/licenses/lgpl-3.0.en.html
8//
9//  Unless required by applicable law or agreed to in writing, software
10//  distributed under the License is distributed on an "AS IS" BASIS,
11//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12//  See the License for the specific language governing permissions and
13//  limitations under the License.
14// -------------------------------------------------------------------------------------------------
15
16/// Calculates the interpolation weight between `x1` and `x2` for a value `x`.
17///
18/// The returned weight `w` satisfies `y = (1 - w) * y1 + w * y2` when
19/// interpolating ordinates that correspond to abscissas `x1` and `x2`.
20///
21/// # Panics
22///
23/// Panics if `x1 == x2` because the denominator becomes zero.
24#[inline]
25#[must_use]
26pub fn linear_weight(x1: f64, x2: f64, x: f64) -> f64 {
27    assert!(
28        x1 != x2,
29        "`x1` and `x2` must differ to compute a linear weight"
30    );
31    (x - x1) / (x2 - x1)
32}
33
34#[inline]
35#[must_use]
36pub fn linear_weighting(y1: f64, y2: f64, x1_diff: f64) -> f64 {
37    x1_diff.mul_add(y2 - y1, y1)
38}
39
40#[inline]
41#[must_use]
42pub fn pos_search(x: f64, xs: &[f64]) -> usize {
43    let n_elem = xs.len();
44    let pos = xs.partition_point(|&val| val < x);
45    std::cmp::min(std::cmp::max(pos.saturating_sub(1), 0), n_elem - 1)
46}
47
48/// Evaluates the quadratic Lagrange polynomial defined by three points.
49///
50/// Given points `(x0, y0)`, `(x1, y1)`, `(x2, y2)` this returns *P(x)* where
51/// *P* is the unique polynomial of degree ≤ 2 passing through the three
52/// points.
53///
54/// # Panics
55///
56/// Panics if any two abscissas are identical because the interpolation
57/// coefficients would involve division by zero.
58#[inline]
59#[must_use]
60pub fn quad_polynomial(x: f64, x0: f64, x1: f64, x2: f64, y0: f64, y1: f64, y2: f64) -> f64 {
61    // Protect against coincident x values that would lead to division by zero
62    assert!(
63        x0 != x1 && x0 != x2 && x1 != x2,
64        "Abscissas must be distinct"
65    );
66
67    y0 * (x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2))
68        + y1 * (x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2))
69        + y2 * (x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1))
70}
71
72/// Performs quadratic interpolation for the point `x` given vectors of abscissas `xs` and ordinates `ys`.
73///
74/// # Panics
75///
76/// Panics if `xs.len() < 3` or `xs.len() != ys.len()`.
77#[must_use]
78pub fn quadratic_interpolation(x: f64, xs: &[f64], ys: &[f64]) -> f64 {
79    let n_elem = xs.len();
80    let epsilon = 1e-8;
81
82    assert!(
83        (n_elem >= 3),
84        "Need at least 3 points for quadratic interpolation"
85    );
86
87    if x <= xs[0] {
88        return ys[0];
89    }
90
91    if x >= xs[n_elem - 1] {
92        return ys[n_elem - 1];
93    }
94
95    let pos = pos_search(x, xs);
96
97    if (xs[pos] - x).abs() < epsilon {
98        return ys[pos];
99    }
100
101    if pos == 0 {
102        return quad_polynomial(x, xs[0], xs[1], xs[2], ys[0], ys[1], ys[2]);
103    }
104
105    if pos == n_elem - 2 {
106        return quad_polynomial(
107            x,
108            xs[n_elem - 3],
109            xs[n_elem - 2],
110            xs[n_elem - 1],
111            ys[n_elem - 3],
112            ys[n_elem - 2],
113            ys[n_elem - 1],
114        );
115    }
116
117    let w = linear_weight(xs[pos], xs[pos + 1], x);
118
119    linear_weighting(
120        quad_polynomial(
121            x,
122            xs[pos - 1],
123            xs[pos],
124            xs[pos + 1],
125            ys[pos - 1],
126            ys[pos],
127            ys[pos + 1],
128        ),
129        quad_polynomial(
130            x,
131            xs[pos],
132            xs[pos + 1],
133            xs[pos + 2],
134            ys[pos],
135            ys[pos + 1],
136            ys[pos + 2],
137        ),
138        w,
139    )
140}
141
142////////////////////////////////////////////////////////////////////////////////
143// Tests
144////////////////////////////////////////////////////////////////////////////////
145#[cfg(test)]
146mod tests {
147    use rstest::*;
148
149    use super::*;
150
151    #[rstest]
152    #[should_panic(expected = "must differ to compute a linear weight")]
153    fn test_linear_weight_zero_divisor() {
154        let _ = linear_weight(1.0, 1.0, 0.5);
155    }
156
157    #[rstest]
158    #[should_panic(expected = "Abscissas must be distinct")]
159    fn test_quad_polynomial_duplicate_x() {
160        let _ = quad_polynomial(0.5, 1.0, 1.0, 2.0, 0.0, 1.0, 4.0);
161    }
162}
nautilus_core/math.rs

nautilus_core/
math.rs
