nautilus_core/

math.rs

// -------------------------------------------------------------------------------------------------
//  Copyright (C) 2015-2026 Nautech Systems Pty Ltd. All rights reserved.
//  https://nautechsystems.io
//
//  Licensed under the GNU Lesser General Public License Version 3.0 (the "License");
//  You may not use this file except in compliance with the License.
//  You may obtain a copy of the License at https://www.gnu.org/licenses/lgpl-3.0.en.html
//
//  Unless required by applicable law or agreed to in writing, software
//  distributed under the License is distributed on an "AS IS" BASIS,
//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//  See the License for the specific language governing permissions and
//  limitations under the License.
// -------------------------------------------------------------------------------------------------

//! Mathematical functions and interpolation utilities.
//!
//! This module provides essential mathematical operations for quantitative trading,
//! including linear and quadratic interpolation functions commonly used in financial
//! data processing and analysis.
//!
//! # Epsilon Values
//!
//! Two epsilon thresholds are used in this module:
//!
//! - **`f64::EPSILON * 2.0` (~4.44e-16), scaled by magnitude:** Used for detecting
//!   near-coincident abscissas in `linear_weight` and `quad_polynomial` to prevent
//!   division instability. The threshold scales with `max(1, |x1|, |x2|)` so that
//!   large-magnitude inputs whose difference is below their floating-point
//!   resolution are also rejected.
//!
//! - **`1e-8`:** Used in `quadratic_interpolation` for detecting exact sample points.
//!   This is an application-level threshold appropriate for typical financial data.

/// Macro for approximate floating-point equality comparison.
///
/// This macro compares two floating-point values with a specified epsilon tolerance,
/// providing a safe alternative to exact equality checks which can fail due to
/// floating-point precision issues.
///
/// # Usage
///
/// ```rust
/// use nautilus_core::approx_eq;
///
/// let a = 0.1 + 0.2;
/// let b = 0.3;
/// assert!(approx_eq!(f64, a, b, epsilon = 1e-10));
/// ```
#[macro_export]
macro_rules! approx_eq {
    ($type:ty, $left:expr, $right:expr, epsilon = $epsilon:expr) => {{
        let left_val: $type = $left;
        let right_val: $type = $right;
        (left_val - right_val).abs() < $epsilon
    }};
}

/// Calculates the interpolation weight between `x1` and `x2` for a value `x`.
///
/// The returned weight `w` satisfies `y = (1 - w) * y1 + w * y2` when
/// interpolating ordinates that correspond to abscissas `x1` and `x2`.
///
/// # Panics
///
/// Panics if:
/// - Any input is NaN or infinite.
/// - `x1` and `x2` are too close relative to their magnitude, which would
///   cause division by zero or numerical instability.
#[inline]
#[must_use]
pub fn linear_weight(x1: f64, x2: f64, x: f64) -> f64 {
    assert!(
        x1.is_finite() && x2.is_finite() && x.is_finite(),
        "All inputs must be finite: x1={x1}, x2={x2}, x={x}"
    );

    assert!(
        !too_close_for_interpolation(x1, x2),
        "`x1` ({x1}) and `x2` ({x2}) are too close for stable interpolation"
    );
    (x - x1) / (x2 - x1)
}

/// Performs linear interpolation using a weight factor.
///
/// Given ordinates `y1` and `y2` and a weight `x1_diff`, computes the
/// interpolated value using the formula: `y1 + x1_diff * (y2 - y1)`.
#[inline]
#[must_use]
pub fn linear_weighting(y1: f64, y2: f64, x1_diff: f64) -> f64 {
    x1_diff.mul_add(y2 - y1, y1)
}

/// Finds the position for interpolation in a sorted array.
///
/// Returns the index of the largest element in `xs` that is less than `x`,
/// clamped to the valid range `[0, xs.len() - 1]`.
///
/// # Edge Cases
///
/// - For empty arrays, returns 0
/// - For single-element arrays, always returns index 0, regardless of whether `x > xs[0]`
/// - For values below the minimum, returns 0
/// - For values at or above the maximum, returns `xs.len() - 1`
#[inline]
#[must_use]
pub fn pos_search(x: f64, xs: &[f64]) -> usize {
    if xs.is_empty() {
        return 0;
    }

    let n_elem = xs.len();
    let pos = xs.partition_point(|&val| val < x);
    pos.saturating_sub(1).min(n_elem - 1)
}

/// Evaluates the quadratic Lagrange polynomial defined by three points.
///
/// Given points `(x0, y0)`, `(x1, y1)`, `(x2, y2)` this returns *P(x)* where
/// *P* is the unique polynomial of degree ≤ 2 passing through the three
/// points.
///
/// # Panics
///
/// Panics if:
/// - Any input is NaN or infinite.
/// - Any two abscissas are too close relative to their magnitude, which would
///   cause division by zero or numerical instability.
#[inline]
#[must_use]
pub fn quad_polynomial(x: f64, x0: f64, x1: f64, x2: f64, y0: f64, y1: f64, y2: f64) -> f64 {
    assert!(
        x.is_finite()
            && x0.is_finite()
            && x1.is_finite()
            && x2.is_finite()
            && y0.is_finite()
            && y1.is_finite()
            && y2.is_finite(),
        "All inputs must be finite: x={x}, x0={x0}, x1={x1}, x2={x2}, y0={y0}, y1={y1}, y2={y2}"
    );

    // Protect against coincident x values that would lead to division by zero
    assert!(
        !too_close_for_interpolation(x0, x1)
            && !too_close_for_interpolation(x0, x2)
            && !too_close_for_interpolation(x1, x2),
        "Abscissas are too close for stable interpolation: x0={x0}, x1={x1}, x2={x2}"
    );

    y0 * (x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2))
        + y1 * (x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2))
        + y2 * (x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1))
}

/// Performs quadratic interpolation for the point `x` given vectors of abscissas `xs` and ordinates `ys`.
///
/// # Panics
///
/// Panics if `xs.len() < 3` or `xs.len() != ys.len()`.
#[must_use]
pub fn quadratic_interpolation(x: f64, xs: &[f64], ys: &[f64]) -> f64 {
    let n_elem = xs.len();
    let epsilon = 1e-8;

    assert!(
        n_elem >= 3,
        "Need at least 3 points for quadratic interpolation"
    );
    assert_eq!(xs.len(), ys.len(), "xs and ys must have the same length");

    if x <= xs[0] {
        return ys[0];
    }

    if x >= xs[n_elem - 1] {
        return ys[n_elem - 1];
    }

    let pos = pos_search(x, xs);

    if (xs[pos] - x).abs() < epsilon {
        return ys[pos];
    }

    if pos == 0 {
        return quad_polynomial(x, xs[0], xs[1], xs[2], ys[0], ys[1], ys[2]);
    }

    if pos == n_elem - 2 {
        return quad_polynomial(
            x,
            xs[n_elem - 3],
            xs[n_elem - 2],
            xs[n_elem - 1],
            ys[n_elem - 3],
            ys[n_elem - 2],
            ys[n_elem - 1],
        );
    }

    let w = linear_weight(xs[pos], xs[pos + 1], x);

    linear_weighting(
        quad_polynomial(
            x,
            xs[pos - 1],
            xs[pos],
            xs[pos + 1],
            ys[pos - 1],
            ys[pos],
            ys[pos + 1],
        ),
        quad_polynomial(
            x,
            xs[pos],
            xs[pos + 1],
            xs[pos + 2],
            ys[pos],
            ys[pos + 1],
            ys[pos + 2],
        ),
        w,
    )
}

/// Returns `true` if `x1` and `x2` are too close, relative to their magnitude,
/// for numerically stable interpolation.
///
/// The threshold scales with `max(1, |x1|, |x2|)`, giving an absolute floor of
/// `2 * f64::EPSILON` for small values and a relative bound for large values.
/// [`linear_weight`] and [`quad_polynomial`] panic when this returns `true` for
/// their abscissas.
#[inline]
#[must_use]
pub fn too_close_for_interpolation(x1: f64, x2: f64) -> bool {
    const EPSILON: f64 = f64::EPSILON * 2.0; // ~4.44e-16
    let scale = 1.0_f64.max(x1.abs()).max(x2.abs());
    (x2 - x1).abs() < EPSILON * scale
}

#[cfg(test)]
mod tests {
    use rstest::*;

    use super::*;

    #[rstest]
    #[case(0.0, 10.0, 5.0, 0.5)]
    #[case(1.0, 3.0, 2.0, 0.5)]
    #[case(0.0, 1.0, 0.25, 0.25)]
    #[case(0.0, 1.0, 0.75, 0.75)]
    fn test_linear_weight_valid_cases(
        #[case] x1: f64,
        #[case] x2: f64,
        #[case] x: f64,
        #[case] expected: f64,
    ) {
        let result = linear_weight(x1, x2, x);
        assert!(
            approx_eq!(f64, result, expected, epsilon = 1e-10),
            "Expected {expected}, was {result}"
        );
    }

    #[rstest]
    #[should_panic(expected = "too close for stable interpolation")]
    fn test_linear_weight_zero_divisor() {
        let _ = linear_weight(1.0, 1.0, 0.5);
    }

    #[rstest]
    #[should_panic(expected = "too close for stable interpolation")]
    fn test_linear_weight_near_equal_values() {
        // Values within machine epsilon should be rejected
        let _ = linear_weight(1.0, 1.0 + f64::EPSILON, 0.5);
    }

    #[rstest]
    fn test_linear_weight_with_small_differences() {
        // High-resolution data (e.g., nanosecond timestamps as seconds) should work
        let result = linear_weight(0.0, 1e-12, 5e-13);
        assert!(result.is_finite());
        assert!((result - 0.5).abs() < 1e-10); // Should be approximately 0.5
    }

    #[rstest]
    fn test_linear_weight_just_above_epsilon() {
        // Values differing by more than machine epsilon should work
        let result = linear_weight(1.0, 1.0 + 1e-9, 1.0 + 5e-10);
        // Should not panic and return a reasonable value
        assert!(result.is_finite());
    }

    #[rstest]
    #[should_panic(expected = "too close for stable interpolation")]
    fn test_linear_weight_large_magnitude_near_equal_panics() {
        // At magnitude 1e16 the f64 resolution is 2.0, so a spacing of 2.0 is
        // below the scaled threshold (~4.44) and must be rejected
        let _ = linear_weight(1e16, 1e16 + 2.0, 1e16 + 1.0);
    }

    #[rstest]
    fn test_linear_weight_large_magnitude_adequate_spacing() {
        let result = linear_weight(1e16, 1e16 + 1e10, 1e16 + 5e9);
        assert!(result.is_finite());
        assert!((result - 0.5).abs() < 1e-6);
    }

    #[rstest]
    #[case(1.0, 3.0, 0.5, 2.0)]
    #[case(10.0, 20.0, 0.25, 12.5)]
    #[case(0.0, 10.0, 0.0, 0.0)]
    #[case(0.0, 10.0, 1.0, 10.0)]
    fn test_linear_weighting(
        #[case] y1: f64,
        #[case] y2: f64,
        #[case] weight: f64,
        #[case] expected: f64,
    ) {
        let result = linear_weighting(y1, y2, weight);
        assert!(
            approx_eq!(f64, result, expected, epsilon = 1e-10),
            "Expected {expected}, was {result}"
        );
    }

    #[rstest]
    #[case(5.0, &[1.0, 2.0, 3.0, 4.0, 6.0, 7.0], 3)]
    #[case(1.5, &[1.0, 2.0, 3.0, 4.0], 0)]
    #[case(0.5, &[1.0, 2.0, 3.0, 4.0], 0)]
    #[case(4.5, &[1.0, 2.0, 3.0, 4.0], 3)]
    #[case(2.0, &[1.0, 2.0, 3.0, 4.0], 0)]
    fn test_pos_search(#[case] x: f64, #[case] xs: &[f64], #[case] expected: usize) {
        let result = pos_search(x, xs);
        assert_eq!(result, expected);
    }

    #[rstest]
    fn test_pos_search_edge_cases() {
        // Single element array
        let result = pos_search(5.0, &[10.0]);
        assert_eq!(result, 0);

        // Value at exact boundary
        let result = pos_search(3.0, &[1.0, 2.0, 3.0, 4.0]);
        assert_eq!(result, 1); // Index of largest element < 3.0 is index 1 (value 2.0)

        // Two element array
        let result = pos_search(1.5, &[1.0, 2.0]);
        assert_eq!(result, 0);
    }

    #[rstest]
    fn test_pos_search_empty_slice() {
        let empty: &[f64] = &[];
        assert_eq!(pos_search(42.0, empty), 0);
    }

    #[rstest]
    fn test_quad_polynomial_linear_case() {
        // Test with three collinear points - should behave like linear interpolation
        let result = quad_polynomial(1.5, 1.0, 2.0, 3.0, 1.0, 2.0, 3.0);
        assert!(approx_eq!(f64, result, 1.5, epsilon = 1e-10));
    }

    #[rstest]
    fn test_quad_polynomial_parabola() {
        // Test with a simple parabola y = x^2
        // Points: (0,0), (1,1), (2,4)
        let result = quad_polynomial(1.5, 0.0, 1.0, 2.0, 0.0, 1.0, 4.0);
        let expected = 1.5 * 1.5; // Should be 2.25
        assert!(approx_eq!(f64, result, expected, epsilon = 1e-10));
    }

    #[rstest]
    #[should_panic(expected = "too close for stable interpolation")]
    fn test_quad_polynomial_duplicate_x() {
        let _ = quad_polynomial(0.5, 1.0, 1.0, 2.0, 0.0, 1.0, 4.0);
    }

    #[rstest]
    #[should_panic(expected = "too close for stable interpolation")]
    fn test_quad_polynomial_near_equal_x_values() {
        // x0 and x1 differ by only machine epsilon
        let _ = quad_polynomial(0.5, 1.0, 1.0 + f64::EPSILON, 2.0, 0.0, 1.0, 4.0);
    }

    #[rstest]
    fn test_quad_polynomial_with_small_differences() {
        // High-resolution data should work (e.g., 1e-12 spacing)
        let result = quad_polynomial(5e-13, 0.0, 1e-12, 2e-12, 0.0, 1.0, 4.0);
        assert!(result.is_finite());
    }

    #[rstest]
    fn test_quad_polynomial_just_above_epsilon() {
        // Values differing by more than machine epsilon should work
        let result = quad_polynomial(0.5, 0.0, 1.0 + 1e-9, 2.0, 0.0, 1.0, 4.0);
        // Should not panic and return a reasonable value
        assert!(result.is_finite());
    }

    #[rstest]
    #[should_panic(expected = "too close for stable interpolation")]
    fn test_quad_polynomial_large_magnitude_near_equal_panics() {
        // x0 and x1 are adjacent representable values at magnitude 1e16
        let _ = quad_polynomial(1e16, 1e16, 1e16 + 2.0, 2e16, 0.0, 1.0, 4.0);
    }

    #[rstest]
    // `allow` not `expect`: nightly clippy does not fire `float_cmp` inside `assert_eq!`
    #[allow(
        clippy::float_cmp,
        reason = "boundary inputs must return the exact boundary ys value"
    )]
    fn test_quadratic_interpolation_boundary_conditions() {
        let xs = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let ys = vec![1.0, 4.0, 9.0, 16.0, 25.0]; // y = x^2

        // Test below minimum
        let result = quadratic_interpolation(0.5, &xs, &ys);
        assert_eq!(result, ys[0]);

        // Test above maximum
        let result = quadratic_interpolation(6.0, &xs, &ys);
        assert_eq!(result, ys[4]);
    }

    #[rstest]
    fn test_quadratic_interpolation_exact_points() {
        let xs = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let ys = vec![1.0, 4.0, 9.0, 16.0, 25.0];

        // Test exact points
        for (i, &x) in xs.iter().enumerate() {
            let result = quadratic_interpolation(x, &xs, &ys);
            assert!(approx_eq!(f64, result, ys[i], epsilon = 1e-6));
        }
    }

    #[rstest]
    fn test_quadratic_interpolation_intermediate_values() {
        let xs = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let ys = vec![1.0, 4.0, 9.0, 16.0, 25.0]; // y = x^2

        // Test interpolation between points
        let result = quadratic_interpolation(2.5, &xs, &ys);
        let expected = 2.5 * 2.5; // Should be close to 6.25
        assert!((result - expected).abs() < 0.1); // Allow some interpolation error
    }

    #[rstest]
    #[should_panic(expected = "Need at least 3 points")]
    fn test_quadratic_interpolation_insufficient_points() {
        let xs = vec![1.0, 2.0];
        let ys = vec![1.0, 4.0];
        let _ = quadratic_interpolation(1.5, &xs, &ys);
    }

    #[rstest]
    #[should_panic(expected = "xs and ys must have the same length")]
    fn test_quadratic_interpolation_mismatched_lengths() {
        let xs = vec![1.0, 2.0, 3.0];
        let ys = vec![1.0, 4.0];
        let _ = quadratic_interpolation(1.5, &xs, &ys);
    }

    #[rstest]
    #[case(f64::NAN, 0.0, 1.0)]
    #[case(0.0, f64::NAN, 1.0)]
    #[case(0.0, 1.0, f64::NAN)]
    #[case(f64::INFINITY, 0.0, 1.0)]
    #[case(0.0, f64::NEG_INFINITY, 1.0)]
    #[should_panic(expected = "All inputs must be finite")]
    fn test_linear_weight_non_finite_panics(#[case] x1: f64, #[case] x2: f64, #[case] x: f64) {
        let _ = linear_weight(x1, x2, x);
    }

    #[rstest]
    #[should_panic(expected = "All inputs must be finite")]
    fn test_quad_polynomial_nan_panics() {
        let _ = quad_polynomial(f64::NAN, 0.0, 1.0, 2.0, 0.0, 1.0, 4.0);
    }

    #[rstest]
    #[should_panic(expected = "All inputs must be finite")]
    fn test_quad_polynomial_infinity_panics() {
        let _ = quad_polynomial(0.5, f64::INFINITY, 1.0, 2.0, 0.0, 1.0, 4.0);
    }
}