scirs2_interpolate/

utils.rs

//! Utility functions for interpolation
//!
//! This module provides helper functions for interpolation tasks.

use crate::error::{InterpolateError, InterpolateResult};
use scirs2_core::ndarray::{Array1, ArrayView1};
use scirs2_core::numeric::{Float, FromPrimitive};
use scirs2_core::safe_ops::{safe_divide, safe_sqrt};
use std::fmt::{Debug, Display};

/// Compute the error estimate for interpolation
///
/// This function performs leave-one-out cross-validation to estimate
/// the interpolation error. It removes each data point in turn, fits
/// the interpolation to the remaining data, and measures the prediction
/// error at the removed point.
///
/// # Arguments
///
/// * `x` - Original x coordinates
/// * `y` - Original y values  
/// * `interp_fn` - Function that performs interpolation at given x values
///
/// # Returns
///
/// Root mean square error (RMSE) from leave-one-out cross-validation
///
/// # Examples
///
/// ```rust
/// use scirs2_core::ndarray::{Array1, ArrayView1};
/// use scirs2_interpolate::utils::error_estimate;
/// use scirs2_interpolate::error::InterpolateResult;
///
/// // Sample data with some noise
/// let x = Array1::from_vec(vec![0.0, 1.0, 2.0, 3.0, 4.0]);
/// let y = Array1::from_vec(vec![0.1, 0.9, 2.1, 2.9, 4.1]);
///
/// // Define a simple linear interpolation function
/// let linear_interp = |x_train: &ArrayView1<f64>, y_train: &ArrayView1<f64>, x_test: &ArrayView1<f64>| -> InterpolateResult<Array1<f64>> {
///     // Simplified linear interpolation implementation
///     let mut result = Array1::zeros(x_test.len());
///     for (i, &x_val) in x_test.iter().enumerate() {
///         // Find nearest neighbors and interpolate
///         if x_train.len() >= 2 {
///             result[i] = x_val; // Simplified: just return x for y=x function
///         }
///     }
///     Ok(result)
/// };
///
/// let rmse = error_estimate(&x.view(), &y.view(), linear_interp).expect("Operation failed");
/// println!("Cross-validation RMSE: {}", rmse);
/// ```
#[allow(dead_code)]
pub fn error_estimate<F, Func>(
    x: &ArrayView1<F>,
    y: &ArrayView1<F>,
    interp_fn: Func,
) -> InterpolateResult<F>
where
    F: Float + FromPrimitive + Debug + Display,
    Func: Fn(&ArrayView1<F>, &ArrayView1<F>, &ArrayView1<F>) -> InterpolateResult<Array1<F>>,
{
    if x.len() != y.len() {
        return Err(InterpolateError::invalid_input(
            "x and y arrays must have the same length",
        ));
    }

    if x.len() < 3 {
        return Err(InterpolateError::insufficient_points(
            3,
            x.len(),
            "interpolation error estimation",
        ));
    }

    let mut sum_squared_error = F::zero();
    let n = x.len();

    for i in 0..n {
        // Create leave-one-out dataset
        let mut x_loo = Vec::with_capacity(n - 1);
        let mut y_loo = Vec::with_capacity(n - 1);

        for j in 0..n {
            if i != j {
                x_loo.push(x[j]);
                y_loo.push(y[j]);
            }
        }

        let x_loo_array = Array1::from_vec(x_loo);
        let y_loo_array = Array1::from_vec(y_loo);

        // Predict at the left-out point
        let x_test = Array1::from_vec(vec![x[i]]);
        let y_pred = interp_fn(&x_loo_array.view(), &y_loo_array.view(), &x_test.view())?;

        // Compute squared error
        let error = y_pred[0] - y[i];
        sum_squared_error = sum_squared_error + error * error;
    }

    // Return RMSE
    let n_f = F::from_usize(n).ok_or_else(|| {
        InterpolateError::ComputationError(
            "Failed to convert array length to float type".to_string(),
        )
    })?;

    let variance = safe_divide(sum_squared_error, n_f).map_err(|_| {
        InterpolateError::ComputationError("Division by zero in RMSE calculation".to_string())
    })?;

    let rmse = safe_sqrt(variance).map_err(|_| {
        InterpolateError::ComputationError(
            "Square root of negative value in RMSE calculation".to_string(),
        )
    })?;

    Ok(rmse)
}

/// Find optimal interpolation parameters
///
/// # Arguments
///
/// * `x` - Original x coordinates
/// * `y` - Original y values
/// * `param_values` - Array of parameter values to try
/// * `interp_fn_builder` - Function that builds an interpolation function with a parameter
///
/// # Returns
///
/// The parameter value that minimizes the cross-validation error
#[allow(dead_code)]
pub fn optimize_parameter<F, Func, BuilderFunc>(
    x: &ArrayView1<F>,
    y: &ArrayView1<F>,
    param_values: &ArrayView1<F>,
    interp_fn_builder: BuilderFunc,
) -> InterpolateResult<F>
where
    F: Float + FromPrimitive + Debug + Display,
    Func: Fn(&ArrayView1<F>, &ArrayView1<F>, &ArrayView1<F>) -> InterpolateResult<Array1<F>>,
    BuilderFunc: Fn(F) -> Func,
{
    if param_values.is_empty() {
        return Err(InterpolateError::invalid_input(
            "at least one parameter value must be provided",
        ));
    }

    let mut best_param = param_values[0];
    let mut min_error = F::infinity();

    for &param in param_values.iter() {
        let interp_fn = interp_fn_builder(param);
        let error = error_estimate(x, y, interp_fn)?;

        if error < min_error {
            min_error = error;
            best_param = param;
        }
    }

    Ok(best_param)
}

/// Differentiate an interpolated function using finite differences
///
/// This function computes the derivative of an interpolated function
/// at a given point using central finite differences. This is useful
/// when you have an interpolant but need its derivative.
///
/// # Arguments
///
/// * `x` - Point at which to evaluate the derivative
/// * `h` - Step size for the finite difference (smaller = more accurate, but numerical issues)
/// * `evalfn` - Function that evaluates the interpolant at a point
///
/// # Returns
///
/// The approximate derivative of the interpolant at x
///
/// # Examples
///
/// ```rust
/// use scirs2_interpolate::utils::differentiate;
/// use scirs2_interpolate::error::InterpolateResult;
///
/// // Example: differentiate f(x) = x^3 at x = 2
/// // Expected derivative: f'(2) = 3 * 2^2 = 12
/// let cubic_fn = |x: f64| -> InterpolateResult<f64> {
///     Ok(x * x * x)
/// };
///
/// let derivative_at_2 = differentiate(2.0, 0.001, cubic_fn).expect("Operation failed");
/// assert!((derivative_at_2 - 12.0).abs() < 0.01); // Should be close to 12
///
/// // Example: differentiate sin(x) at x = π/2  
/// // Expected derivative: cos(π/2) = 0
/// let sin_fn = |x: f64| -> InterpolateResult<f64> {
///     Ok(x.sin())
/// };
///
/// let derivative_at_pi_2 = differentiate(std::f64::consts::PI / 2.0, 0.0001, sin_fn).expect("Operation failed");
/// assert!(derivative_at_pi_2.abs() < 0.01); // Should be close to 0
/// ```
#[allow(dead_code)]
pub fn differentiate<F, Func>(x: F, h: F, evalfn: Func) -> InterpolateResult<F>
where
    F: Float + FromPrimitive + Debug + Display,
    Func: Fn(F) -> InterpolateResult<F>,
{
    // Use central difference for better accuracy
    let f_plus = evalfn(x + h)?;
    let f_minus = evalfn(x - h)?;

    let two = F::from_f64(2.0).ok_or_else(|| {
        InterpolateError::ComputationError(
            "Failed to convert constant 2.0 to float type".to_string(),
        )
    })?;

    let denominator = two * h;
    let derivative = safe_divide(f_plus - f_minus, denominator).map_err(|_| {
        InterpolateError::ComputationError(
            "Division by zero in finite difference calculation (step size too small)".to_string(),
        )
    })?;

    Ok(derivative)
}

/// Integrate an interpolated function using Simpson's rule
///
/// This function computes the definite integral of an interpolated function
/// over a specified interval using composite Simpson's rule. This is useful
/// for computing areas under interpolated curves or other integral quantities.
///
/// # Arguments
///
/// * `a` - Lower bound of integration
/// * `b` - Upper bound of integration  
/// * `n` - Number of intervals for the quadrature (must be even and >= 2)
/// * `evalfn` - Function that evaluates the interpolant at a point
///
/// # Returns
///
/// The approximate definite integral of the interpolant from a to b
///
/// # Examples
///
/// ```rust
/// use scirs2_interpolate::utils::integrate;
/// use scirs2_interpolate::error::InterpolateResult;
///
/// // Example: integrate f(x) = x^2 from 0 to 2
/// // Expected result: ∫₀² x² dx = [x³/3]₀² = 8/3 ≈ 2.667
/// let quadratic_fn = |x: f64| -> InterpolateResult<f64> {
///     Ok(x * x)
/// };
///
/// let integral = integrate(0.0, 2.0, 100, quadratic_fn).expect("Operation failed");
/// assert!((integral - 8.0/3.0).abs() < 0.001);
///
/// // Example: integrate sin(x) from 0 to π
/// // Expected result: ∫₀^π sin(x) dx = [-cos(x)]₀^π = 2
/// let sin_fn = |x: f64| -> InterpolateResult<f64> {
///     Ok(x.sin())
/// };
///
/// let integral_sin = integrate(0.0, std::f64::consts::PI, 200, sin_fn).expect("Operation failed");
/// assert!((integral_sin - 2.0).abs() < 0.001);
/// ```
#[allow(dead_code)]
pub fn integrate<F, Func>(a: F, b: F, n: usize, evalfn: Func) -> InterpolateResult<F>
where
    F: Float + FromPrimitive + Debug + Display,
    Func: Fn(F) -> InterpolateResult<F>,
{
    if a > b {
        return integrate(b, a, n, evalfn).map(|result| -result);
    }

    // Use composite Simpson's rule for integration
    if n < 2 {
        return Err(InterpolateError::InvalidValue(
            "number of intervals must be at least 2".to_string(),
        ));
    }

    if !n.is_multiple_of(2) {
        return Err(InterpolateError::InvalidValue(
            "number of intervals must be even".to_string(),
        ));
    }

    let n_f = F::from_usize(n).ok_or_else(|| {
        InterpolateError::ComputationError(
            "Failed to convert number of intervals to float type".to_string(),
        )
    })?;

    let h = safe_divide(b - a, n_f).map_err(|_| {
        InterpolateError::ComputationError(
            "Division by zero in step size calculation (zero intervals)".to_string(),
        )
    })?;

    let mut sum = evalfn(a)? + evalfn(b)?;

    // Even-indexed points (except endpoints)
    let two = F::from_f64(2.0).ok_or_else(|| {
        InterpolateError::ComputationError(
            "Failed to convert constant 2.0 to float type".to_string(),
        )
    })?;

    for i in 1..n {
        if i % 2 == 0 {
            let i_f = F::from_usize(i).ok_or_else(|| {
                InterpolateError::ComputationError(
                    "Failed to convert index to float type".to_string(),
                )
            })?;
            let x_i = a + i_f * h;
            sum = sum + two * evalfn(x_i)?;
        }
    }

    // Odd-indexed points
    let four = F::from_f64(4.0).ok_or_else(|| {
        InterpolateError::ComputationError(
            "Failed to convert constant 4.0 to float type".to_string(),
        )
    })?;

    for i in 1..n {
        if i % 2 == 1 {
            let i_f = F::from_usize(i).ok_or_else(|| {
                InterpolateError::ComputationError(
                    "Failed to convert index to float type".to_string(),
                )
            })?;
            let x_i = a + i_f * h;
            sum = sum + four * evalfn(x_i)?;
        }
    }

    let three = F::from_f64(3.0).ok_or_else(|| {
        InterpolateError::ComputationError(
            "Failed to convert constant 3.0 to float type".to_string(),
        )
    })?;

    let integral = safe_divide(h * sum, three).map_err(|_| {
        InterpolateError::ComputationError(
            "Division by zero in Simpson's rule calculation".to_string(),
        )
    })?;

    Ok(integral)
}

/// Find roots using bisection method
///
/// This function uses the bisection method to find all roots of a function within a given interval.
///
/// # Arguments
///
/// * `a` - Left boundary of search interval
/// * `b` - Right boundary of search interval  
/// * `tolerance` - Tolerance for root finding accuracy
/// * `evalfn` - Function to evaluate
///
/// # Returns
///
/// Vector of roots found in the interval
///
#[allow(dead_code)]
pub fn find_roots_bisection<F, Func>(
    a: F,
    b: F,
    tolerance: F,
    evalfn: Func,
) -> InterpolateResult<Vec<F>>
where
    F: Float + FromPrimitive + Debug + Display,
    Func: Fn(F) -> InterpolateResult<F>,
{
    let mut roots = Vec::new();

    if a >= b {
        return Ok(roots);
    }

    // Evaluate at endpoints
    let fa = evalfn(a)?;
    let fb = evalfn(b)?;

    // If either endpoint is close to zero, it's a root
    if fa.abs() < tolerance {
        roots.push(a);
    }
    if fb.abs() < tolerance && (b - a).abs() > tolerance {
        roots.push(b);
    }

    // If signs are the same, no root in interval by intermediate value theorem
    if fa * fb > F::zero() {
        return Ok(roots);
    }

    // Binary search for root
    let mut left = a;
    let mut right = b;
    let mut f_left = fa;
    let mut _f_right = fb;

    while (right - left).abs() > tolerance {
        let mid = left + (right - left) / F::from_f64(2.0).expect("Operation failed");
        let f_mid = evalfn(mid)?;

        if f_mid.abs() < tolerance {
            roots.push(mid);
            break;
        }

        if f_left * f_mid < F::zero() {
            right = mid;
            _f_right = f_mid;
        } else {
            left = mid;
            f_left = f_mid;
        }
    }

    // If we didn't find exact root, add the midpoint
    if roots.is_empty() {
        let root = left + (right - left) / F::from_f64(2.0).expect("Operation failed");
        let f_root = evalfn(root)?;
        if f_root.abs() < tolerance * F::from_f64(10.0).expect("Operation failed") {
            roots.push(root);
        }
    }

    Ok(roots)
}

/// Find multiple roots by subdividing interval
///
/// This function subdivides the interval and searches for roots in each subdivision.
///
/// # Arguments
///
/// * `a` - Left boundary of search interval
/// * `b` - Right boundary of search interval
/// * `tolerance` - Tolerance for root finding accuracy
/// * `subdivisions` - Number of subdivisions to search
/// * `evalfn` - Function to evaluate
///
/// # Returns
///
/// Vector of roots found in the interval
///
#[allow(dead_code)]
pub fn find_multiple_roots<F, Func>(
    a: F,
    b: F,
    tolerance: F,
    subdivisions: usize,
    evalfn: Func,
) -> InterpolateResult<Vec<F>>
where
    F: Float + FromPrimitive + Debug + Display,
    Func: Fn(F) -> InterpolateResult<F> + Copy,
{
    let mut all_roots = Vec::new();

    if subdivisions == 0 {
        return Ok(all_roots);
    }

    let step = (b - a) / F::from_usize(subdivisions).expect("Operation failed");

    for i in 0..subdivisions {
        let left = a + F::from_usize(i).expect("Operation failed") * step;
        let right = a + F::from_usize(i + 1).expect("Operation failed") * step;

        match find_roots_bisection(left, right, tolerance, evalfn) {
            Ok(mut roots) => all_roots.append(&mut roots),
            Err(_) => continue,
        }
    }

    // Sort and remove duplicates
    all_roots.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    all_roots.dedup_by(|a, b| (*a - *b).abs() < tolerance);

    Ok(all_roots)
}

#[cfg(test)]
mod tests {
    use super::*;
    // interpolation functions
    use scirs2_core::ndarray::array;

    #[test]
    fn test_error_estimate() {
        use crate::{cubic_interpolate, linear_interpolate};

        // Linear data y = 2x + 1 (not y == x, so a hypothetical "just echo
        // the input back" bug in error_estimate/interp_fn couldn't
        // coincidentally pass this check).
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0];
        let y = array![1.0, 3.0, 5.0, 7.0, 9.0];

        // Leave-one-out CV error for linear interpolation on perfectly
        // linear data should be ~0 (up to floating point): removing any one
        // point still leaves >= 4 points on the same line, and linear
        // interpolation/extrapolation along that line recovers the held-out
        // point exactly.
        let error =
            error_estimate(&x.view(), &y.view(), linear_interpolate).expect("Operation failed");
        assert!(
            error < 1e-10,
            "linear LOO-CV RMSE should be ~0, got {error}"
        );

        // `cubic_interpolate`'s `Interp1d::Cubic` (see
        // `interp1d::mod::cubic_interp`) is a uniform-parametrization
        // Catmull-Rom spline: its tangent estimate at an interior point p1
        // is `(p2 - p0)`, which only approximates the true derivative when
        // p1 is spaced evenly between its neighbors p0 and p2. Leave-one-out
        // removes one x value at a time, which makes the *remaining* points
        // non-uniformly spaced (e.g. removing x=1 from [0,1,2,3,4] leaves
        // gaps 2,1,1) -- so even though the underlying y=2x+1 is perfectly
        // linear, this specific (uniform-spacing-assuming) formula no longer
        // reproduces it exactly. This is a genuine, deterministic property of
        // the simplified Catmull-Rom implementation (matches the wide 0.3
        // tolerance already documented on `cubic_interpolate`'s own doctest),
        // not a bug to paper over with a near-zero tolerance: the measured
        // RMSE here is a reproducible ~0.0791, so 0.1 gives real headroom
        // while still catching a genuine regression.
        let error =
            error_estimate(&x.view(), &y.view(), cubic_interpolate).expect("Operation failed");
        assert!(
            error < 0.1,
            "cubic LOO-CV RMSE should stay within the known Catmull-Rom \
             non-uniform-spacing error budget, got {error}"
        );
    }

    #[test]
    fn test_differentiate() {
        // Function: f(x) = x^2
        let f = |x: f64| -> InterpolateResult<f64> { Ok(x * x) };

        // At x=2, f'(x) = 2x = 4
        let derivative = differentiate(2.0, 0.001, f).expect("Operation failed");
        assert!((derivative - 4.0).abs() < 1e-5);

        // At x=3, f'(x) = 2x = 6
        let derivative = differentiate(3.0, 0.001, f).expect("Operation failed");
        assert!((derivative - 6.0).abs() < 1e-5);
    }

    #[test]
    fn test_integrate() {
        // Function: f(x) = x^2
        // Integral from 0 to 1: x^3/3 = 1/3
        let f = |x: f64| -> InterpolateResult<f64> { Ok(x * x) };

        let integral = integrate(0.0, 1.0, 100, f).expect("Operation failed");
        assert!((integral - 1.0 / 3.0).abs() < 1e-5);

        // Integral from 0 to 2: x^3/3 = 8/3
        let integral = integrate(0.0, 2.0, 100, f).expect("Operation failed");
        assert!((integral - 8.0 / 3.0).abs() < 1e-5);
    }
}