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

nautilus_core/
math.rs
