nautilus_core/
math.rs

1// -------------------------------------------------------------------------------------------------
2//  Copyright (C) 2015-2025 Nautech Systems Pty Ltd. All rights reserved.
3//  https://nautechsystems.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/// Macro for approximate floating-point equality comparison.
23///
24/// This macro compares two floating-point values with a specified epsilon tolerance,
25/// providing a safe alternative to exact equality checks which can fail due to
26/// floating-point precision issues.
27///
28/// # Usage
29///
30/// ```rust
31/// use nautilus_core::approx_eq;
32///
33/// let a = 0.1 + 0.2;
34/// let b = 0.3;
35/// assert!(approx_eq!(f64, a, b, epsilon = 1e-10));
36/// ```
37#[macro_export]
38macro_rules! approx_eq {
39    ($type:ty, $left:expr, $right:expr, epsilon = $epsilon:expr) => {{
40        let left_val: $type = $left;
41        let right_val: $type = $right;
42        (left_val - right_val).abs() < $epsilon
43    }};
44    ($type:ty, $left:expr, $right:expr, epsilon = $epsilon:expr, ulps = $ulps:expr) => {{
45        let left_val: $type = $left;
46        let right_val: $type = $right;
47        // For compatibility, we use epsilon comparison and ignore ulps
48        (left_val - right_val).abs() < $epsilon
49    }};
50}
51
52/// Calculates the interpolation weight between `x1` and `x2` for a value `x`.
53///
54/// The returned weight `w` satisfies `y = (1 - w) * y1 + w * y2` when
55/// interpolating ordinates that correspond to abscissas `x1` and `x2`.
56///
57/// # Panics
58///
59/// Panics if `x1` and `x2` are too close (within machine epsilon), which would
60/// cause division by zero or numerical instability. Uses `f64::EPSILON` * 2.0 to
61/// account for floating-point rounding in the difference computation.
62#[inline]
63#[must_use]
64pub fn linear_weight(x1: f64, x2: f64, x: f64) -> f64 {
65    const EPSILON: f64 = f64::EPSILON * 2.0; // ~4.44e-16
66    let diff = (x2 - x1).abs();
67    assert!(
68        diff >= EPSILON,
69        "`x1` ({x1}) and `x2` ({x2}) are too close for stable interpolation (diff: {diff}, min: {EPSILON})"
70    );
71    (x - x1) / (x2 - x1)
72}
73
74/// Performs linear interpolation using a weight factor.
75///
76/// Given ordinates `y1` and `y2` and a weight `x1_diff`, computes the
77/// interpolated value using the formula: `y1 + x1_diff * (y2 - y1)`.
78#[inline]
79#[must_use]
80pub fn linear_weighting(y1: f64, y2: f64, x1_diff: f64) -> f64 {
81    x1_diff.mul_add(y2 - y1, y1)
82}
83
84/// Finds the position for interpolation in a sorted array.
85///
86/// Returns the index of the largest element in `xs` that is less than `x`,
87/// clamped to the valid range `[0, xs.len() - 1]`.
88///
89/// # Edge Cases
90///
91/// - For empty arrays, returns 0
92/// - For single-element arrays, always returns index 0, regardless of whether `x > xs[0]`
93/// - For values below the minimum, returns 0
94/// - For values at or above the maximum, returns `xs.len() - 1`
95///
96/// # Examples
97///
98/// ```
99/// use nautilus_core::math::pos_search;
100///
101/// // Normal case: find position between elements
102/// assert_eq!(pos_search(2.5, &[1.0, 2.0, 3.0, 4.0]), 1); // Between indices 1 and 2
103///
104/// // Single element: always returns 0
105/// assert_eq!(pos_search(5.0, &[10.0]), 0); // Even though x < xs[0]
106/// assert_eq!(pos_search(15.0, &[10.0]), 0); // Even though x > xs[0]
107/// ```
108#[inline]
109#[must_use]
110pub fn pos_search(x: f64, xs: &[f64]) -> usize {
111    if xs.is_empty() {
112        return 0;
113    }
114
115    let n_elem = xs.len();
116    let pos = xs.partition_point(|&val| val < x);
117    std::cmp::min(std::cmp::max(pos.saturating_sub(1), 0), n_elem - 1)
118}
119
120/// Evaluates the quadratic Lagrange polynomial defined by three points.
121///
122/// Given points `(x0, y0)`, `(x1, y1)`, `(x2, y2)` this returns *P(x)* where
123/// *P* is the unique polynomial of degree ≤ 2 passing through the three
124/// points.
125///
126/// # Panics
127///
128/// Panics if any two abscissas are too close (within machine epsilon), which would
129/// cause division by zero or numerical instability in the interpolation.
130#[inline]
131#[must_use]
132pub fn quad_polynomial(x: f64, x0: f64, x1: f64, x2: f64, y0: f64, y1: f64, y2: f64) -> f64 {
133    const EPSILON: f64 = f64::EPSILON * 2.0; // ~4.44e-16
134
135    // Protect against coincident x values that would lead to division by zero
136    let diff_01 = (x0 - x1).abs();
137    let diff_02 = (x0 - x2).abs();
138    let diff_12 = (x1 - x2).abs();
139
140    assert!(
141        diff_01 >= EPSILON && diff_02 >= EPSILON && diff_12 >= EPSILON,
142        "Abscissas are too close for stable interpolation: x0={x0}, x1={x1}, x2={x2} (diffs: {diff_01:.2e}, {diff_02:.2e}, {diff_12:.2e}, min: {EPSILON})"
143    );
144
145    y0 * (x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2))
146        + y1 * (x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2))
147        + y2 * (x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1))
148}
149
150/// Performs quadratic interpolation for the point `x` given vectors of abscissas `xs` and ordinates `ys`.
151///
152/// # Panics
153///
154/// Panics if `xs.len() < 3` or `xs.len() != ys.len()`.
155#[must_use]
156pub fn quadratic_interpolation(x: f64, xs: &[f64], ys: &[f64]) -> f64 {
157    let n_elem = xs.len();
158    let epsilon = 1e-8;
159
160    assert!(
161        n_elem >= 3,
162        "Need at least 3 points for quadratic interpolation"
163    );
164    assert_eq!(xs.len(), ys.len(), "xs and ys must have the same length");
165
166    if x <= xs[0] {
167        return ys[0];
168    }
169
170    if x >= xs[n_elem - 1] {
171        return ys[n_elem - 1];
172    }
173
174    let pos = pos_search(x, xs);
175
176    if (xs[pos] - x).abs() < epsilon {
177        return ys[pos];
178    }
179
180    if pos == 0 {
181        return quad_polynomial(x, xs[0], xs[1], xs[2], ys[0], ys[1], ys[2]);
182    }
183
184    if pos == n_elem - 2 {
185        return quad_polynomial(
186            x,
187            xs[n_elem - 3],
188            xs[n_elem - 2],
189            xs[n_elem - 1],
190            ys[n_elem - 3],
191            ys[n_elem - 2],
192            ys[n_elem - 1],
193        );
194    }
195
196    let w = linear_weight(xs[pos], xs[pos + 1], x);
197
198    linear_weighting(
199        quad_polynomial(
200            x,
201            xs[pos - 1],
202            xs[pos],
203            xs[pos + 1],
204            ys[pos - 1],
205            ys[pos],
206            ys[pos + 1],
207        ),
208        quad_polynomial(
209            x,
210            xs[pos],
211            xs[pos + 1],
212            xs[pos + 2],
213            ys[pos],
214            ys[pos + 1],
215            ys[pos + 2],
216        ),
217        w,
218    )
219}
220
221////////////////////////////////////////////////////////////////////////////////
222// Tests
223////////////////////////////////////////////////////////////////////////////////
224#[cfg(test)]
225mod tests {
226    use rstest::*;
227
228    use super::*;
229
230    #[rstest]
231    #[case(0.0, 10.0, 5.0, 0.5)]
232    #[case(1.0, 3.0, 2.0, 0.5)]
233    #[case(0.0, 1.0, 0.25, 0.25)]
234    #[case(0.0, 1.0, 0.75, 0.75)]
235    fn test_linear_weight_valid_cases(
236        #[case] x1: f64,
237        #[case] x2: f64,
238        #[case] x: f64,
239        #[case] expected: f64,
240    ) {
241        let result = linear_weight(x1, x2, x);
242        assert!(
243            approx_eq!(f64, result, expected, epsilon = 1e-10),
244            "Expected {expected}, was {result}"
245        );
246    }
247
248    #[rstest]
249    #[should_panic(expected = "too close for stable interpolation")]
250    fn test_linear_weight_zero_divisor() {
251        let _ = linear_weight(1.0, 1.0, 0.5);
252    }
253
254    #[rstest]
255    #[should_panic(expected = "too close for stable interpolation")]
256    fn test_linear_weight_near_equal_values() {
257        // Values within machine epsilon should be rejected
258        let _ = linear_weight(1.0, 1.0 + f64::EPSILON, 0.5);
259    }
260
261    #[rstest]
262    fn test_linear_weight_with_small_differences() {
263        // High-resolution data (e.g., nanosecond timestamps as seconds) should work
264        let result = linear_weight(0.0, 1e-12, 5e-13);
265        assert!(result.is_finite());
266        assert!((result - 0.5).abs() < 1e-10); // Should be approximately 0.5
267    }
268
269    #[rstest]
270    fn test_linear_weight_just_above_epsilon() {
271        // Values differing by more than machine epsilon should work
272        let result = linear_weight(1.0, 1.0 + 1e-9, 1.0 + 5e-10);
273        // Should not panic and return a reasonable value
274        assert!(result.is_finite());
275    }
276
277    #[rstest]
278    #[case(1.0, 3.0, 0.5, 2.0)]
279    #[case(10.0, 20.0, 0.25, 12.5)]
280    #[case(0.0, 10.0, 0.0, 0.0)]
281    #[case(0.0, 10.0, 1.0, 10.0)]
282    fn test_linear_weighting(
283        #[case] y1: f64,
284        #[case] y2: f64,
285        #[case] weight: f64,
286        #[case] expected: f64,
287    ) {
288        let result = linear_weighting(y1, y2, weight);
289        assert!(
290            approx_eq!(f64, result, expected, epsilon = 1e-10),
291            "Expected {expected}, was {result}"
292        );
293    }
294
295    #[rstest]
296    #[case(5.0, &[1.0, 2.0, 3.0, 4.0, 6.0, 7.0], 3)]
297    #[case(1.5, &[1.0, 2.0, 3.0, 4.0], 0)]
298    #[case(0.5, &[1.0, 2.0, 3.0, 4.0], 0)]
299    #[case(4.5, &[1.0, 2.0, 3.0, 4.0], 3)]
300    #[case(2.0, &[1.0, 2.0, 3.0, 4.0], 0)]
301    fn test_pos_search(#[case] x: f64, #[case] xs: &[f64], #[case] expected: usize) {
302        let result = pos_search(x, xs);
303        assert_eq!(result, expected);
304    }
305
306    #[rstest]
307    fn test_pos_search_edge_cases() {
308        // Single element array
309        let result = pos_search(5.0, &[10.0]);
310        assert_eq!(result, 0);
311
312        // Value at exact boundary
313        let result = pos_search(3.0, &[1.0, 2.0, 3.0, 4.0]);
314        assert_eq!(result, 1); // Index of largest element < 3.0 is index 1 (value 2.0)
315
316        // Two element array
317        let result = pos_search(1.5, &[1.0, 2.0]);
318        assert_eq!(result, 0);
319    }
320
321    #[rstest]
322    fn test_pos_search_empty_slice() {
323        let empty: &[f64] = &[];
324        assert_eq!(pos_search(42.0, empty), 0);
325    }
326
327    #[rstest]
328    fn test_quad_polynomial_linear_case() {
329        // Test with three collinear points - should behave like linear interpolation
330        let result = quad_polynomial(1.5, 1.0, 2.0, 3.0, 1.0, 2.0, 3.0);
331        assert!(approx_eq!(f64, result, 1.5, epsilon = 1e-10));
332    }
333
334    #[rstest]
335    fn test_quad_polynomial_parabola() {
336        // Test with a simple parabola y = x^2
337        // Points: (0,0), (1,1), (2,4)
338        let result = quad_polynomial(1.5, 0.0, 1.0, 2.0, 0.0, 1.0, 4.0);
339        let expected = 1.5 * 1.5; // Should be 2.25
340        assert!(approx_eq!(f64, result, expected, epsilon = 1e-10));
341    }
342
343    #[rstest]
344    #[should_panic(expected = "too close for stable interpolation")]
345    fn test_quad_polynomial_duplicate_x() {
346        let _ = quad_polynomial(0.5, 1.0, 1.0, 2.0, 0.0, 1.0, 4.0);
347    }
348
349    #[rstest]
350    #[should_panic(expected = "too close for stable interpolation")]
351    fn test_quad_polynomial_near_equal_x_values() {
352        // x0 and x1 differ by only machine epsilon
353        let _ = quad_polynomial(0.5, 1.0, 1.0 + f64::EPSILON, 2.0, 0.0, 1.0, 4.0);
354    }
355
356    #[rstest]
357    fn test_quad_polynomial_with_small_differences() {
358        // High-resolution data should work (e.g., 1e-12 spacing)
359        let result = quad_polynomial(5e-13, 0.0, 1e-12, 2e-12, 0.0, 1.0, 4.0);
360        assert!(result.is_finite());
361    }
362
363    #[rstest]
364    fn test_quad_polynomial_just_above_epsilon() {
365        // Values differing by more than machine epsilon should work
366        let result = quad_polynomial(0.5, 0.0, 1.0 + 1e-9, 2.0, 0.0, 1.0, 4.0);
367        // Should not panic and return a reasonable value
368        assert!(result.is_finite());
369    }
370
371    #[rstest]
372    fn test_quadratic_interpolation_boundary_conditions() {
373        let xs = vec![1.0, 2.0, 3.0, 4.0, 5.0];
374        let ys = vec![1.0, 4.0, 9.0, 16.0, 25.0]; // y = x^2
375
376        // Test below minimum
377        let result = quadratic_interpolation(0.5, &xs, &ys);
378        assert_eq!(result, ys[0]);
379
380        // Test above maximum
381        let result = quadratic_interpolation(6.0, &xs, &ys);
382        assert_eq!(result, ys[4]);
383    }
384
385    #[rstest]
386    fn test_quadratic_interpolation_exact_points() {
387        let xs = vec![1.0, 2.0, 3.0, 4.0, 5.0];
388        let ys = vec![1.0, 4.0, 9.0, 16.0, 25.0];
389
390        // Test exact points
391        for (i, &x) in xs.iter().enumerate() {
392            let result = quadratic_interpolation(x, &xs, &ys);
393            assert!(approx_eq!(f64, result, ys[i], epsilon = 1e-6));
394        }
395    }
396
397    #[rstest]
398    fn test_quadratic_interpolation_intermediate_values() {
399        let xs = vec![1.0, 2.0, 3.0, 4.0, 5.0];
400        let ys = vec![1.0, 4.0, 9.0, 16.0, 25.0]; // y = x^2
401
402        // Test interpolation between points
403        let result = quadratic_interpolation(2.5, &xs, &ys);
404        let expected = 2.5 * 2.5; // Should be close to 6.25
405        assert!((result - expected).abs() < 0.1); // Allow some interpolation error
406    }
407
408    #[rstest]
409    #[should_panic(expected = "Need at least 3 points")]
410    fn test_quadratic_interpolation_insufficient_points() {
411        let xs = vec![1.0, 2.0];
412        let ys = vec![1.0, 4.0];
413        let _ = quadratic_interpolation(1.5, &xs, &ys);
414    }
415
416    #[rstest]
417    #[should_panic(expected = "xs and ys must have the same length")]
418    fn test_quadratic_interpolation_mismatched_lengths() {
419        let xs = vec![1.0, 2.0, 3.0];
420        let ys = vec![1.0, 4.0];
421        let _ = quadratic_interpolation(1.5, &xs, &ys);
422    }
423}
nautilus_core/math.rs

nautilus_core/
math.rs
