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//! Mathematical functions and interpolation utilities.
17//!
18//! This module provides essential mathematical operations for quantitative trading,
19//! including linear and quadratic interpolation functions commonly used in financial
20//! data processing and analysis.
21
22/// Calculates the interpolation weight between `x1` and `x2` for a value `x`.
23///
24/// The returned weight `w` satisfies `y = (1 - w) * y1 + w * y2` when
25/// interpolating ordinates that correspond to abscissas `x1` and `x2`.
26///
27/// # Panics
28///
29/// Panics if `x1 == x2` because the denominator becomes zero.
30#[inline]
31#[must_use]
32pub fn linear_weight(x1: f64, x2: f64, x: f64) -> f64 {
33    assert!(
34        x1 != x2,
35        "`x1` and `x2` must differ to compute a linear weight"
36    );
37    (x - x1) / (x2 - x1)
38}
39
40/// Performs linear interpolation using a weight factor.
41///
42/// Given ordinates `y1` and `y2` and a weight `x1_diff`, computes the
43/// interpolated value using the formula: `y1 + x1_diff * (y2 - y1)`.
44#[inline]
45#[must_use]
46pub fn linear_weighting(y1: f64, y2: f64, x1_diff: f64) -> f64 {
47    x1_diff.mul_add(y2 - y1, y1)
48}
49
50/// Finds the position for interpolation in a sorted array.
51///
52/// Returns the index of the largest element in `xs` that is less than `x`,
53/// clamped to the valid range `[0, xs.len() - 1]`.
54#[inline]
55#[must_use]
56pub fn pos_search(x: f64, xs: &[f64]) -> usize {
57    let n_elem = xs.len();
58    let pos = xs.partition_point(|&val| val < x);
59    std::cmp::min(std::cmp::max(pos.saturating_sub(1), 0), n_elem - 1)
60}
61
62/// Evaluates the quadratic Lagrange polynomial defined by three points.
63///
64/// Given points `(x0, y0)`, `(x1, y1)`, `(x2, y2)` this returns *P(x)* where
65/// *P* is the unique polynomial of degree ≤ 2 passing through the three
66/// points.
67///
68/// # Panics
69///
70/// Panics if any two abscissas are identical because the interpolation
71/// coefficients would involve division by zero.
72#[inline]
73#[must_use]
74pub fn quad_polynomial(x: f64, x0: f64, x1: f64, x2: f64, y0: f64, y1: f64, y2: f64) -> f64 {
75    // Protect against coincident x values that would lead to division by zero
76    assert!(
77        x0 != x1 && x0 != x2 && x1 != x2,
78        "Abscissas must be distinct"
79    );
80
81    y0 * (x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2))
82        + y1 * (x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2))
83        + y2 * (x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1))
84}
85
86/// Performs quadratic interpolation for the point `x` given vectors of abscissas `xs` and ordinates `ys`.
87///
88/// # Panics
89///
90/// Panics if `xs.len() < 3` or `xs.len() != ys.len()`.
91#[must_use]
92pub fn quadratic_interpolation(x: f64, xs: &[f64], ys: &[f64]) -> f64 {
93    let n_elem = xs.len();
94    let epsilon = 1e-8;
95
96    assert!(
97        (n_elem >= 3),
98        "Need at least 3 points for quadratic interpolation"
99    );
100    assert_eq!(xs.len(), ys.len(), "xs and ys must have the same length");
101
102    if x <= xs[0] {
103        return ys[0];
104    }
105
106    if x >= xs[n_elem - 1] {
107        return ys[n_elem - 1];
108    }
109
110    let pos = pos_search(x, xs);
111
112    if (xs[pos] - x).abs() < epsilon {
113        return ys[pos];
114    }
115
116    if pos == 0 {
117        return quad_polynomial(x, xs[0], xs[1], xs[2], ys[0], ys[1], ys[2]);
118    }
119
120    if pos == n_elem - 2 {
121        return quad_polynomial(
122            x,
123            xs[n_elem - 3],
124            xs[n_elem - 2],
125            xs[n_elem - 1],
126            ys[n_elem - 3],
127            ys[n_elem - 2],
128            ys[n_elem - 1],
129        );
130    }
131
132    let w = linear_weight(xs[pos], xs[pos + 1], x);
133
134    linear_weighting(
135        quad_polynomial(
136            x,
137            xs[pos - 1],
138            xs[pos],
139            xs[pos + 1],
140            ys[pos - 1],
141            ys[pos],
142            ys[pos + 1],
143        ),
144        quad_polynomial(
145            x,
146            xs[pos],
147            xs[pos + 1],
148            xs[pos + 2],
149            ys[pos],
150            ys[pos + 1],
151            ys[pos + 2],
152        ),
153        w,
154    )
155}
156
157////////////////////////////////////////////////////////////////////////////////
158// Tests
159////////////////////////////////////////////////////////////////////////////////
160#[cfg(test)]
161mod tests {
162    use rstest::*;
163
164    use super::*;
165
166    #[rstest]
167    #[should_panic(expected = "must differ to compute a linear weight")]
168    fn test_linear_weight_zero_divisor() {
169        let _ = linear_weight(1.0, 1.0, 0.5);
170    }
171
172    #[rstest]
173    #[should_panic(expected = "Abscissas must be distinct")]
174    fn test_quad_polynomial_duplicate_x() {
175        let _ = quad_polynomial(0.5, 1.0, 1.0, 2.0, 0.0, 1.0, 4.0);
176    }
177}
nautilus_core/math.rs

nautilus_core/
math.rs
