scirs2-interpolate 0.4.0

//! Trait implementations for cubic splines
//!
//! This module contains the implementations of various interpolation traits
//! for the CubicSpline struct, providing compatibility with the broader
//! interpolation framework.

use crate::error::InterpolateResult;
use crate::traits::{InterpolationFloat, SplineInterpolator, Interpolator};
use scirs2_core::ndarray::{ArrayView1, ArrayView2, Array1, Array2};
use super::core::CubicSpline;

impl<F> SplineInterpolator<F> for CubicSpline<F>
where
    F: InterpolationFloat,
{
    /// Compute derivatives at multiple query points
    ///
    /// This method provides a batch interface for computing derivatives
    /// of a specified order at multiple points.
    ///
    /// # Arguments
    ///
    /// * `querypoints` - 2D array of query points (must have 1 column for 1D splines)
    /// * `order` - Order of derivative to compute
    ///
    /// # Returns
    ///
    /// Vector of derivative values at each query point
    ///
    /// # Errors
    ///
    /// Returns an error if:
    /// - Query points have more than 1 column (cubic splines are 1D only)
    /// - Any query point is outside the interpolation domain
    /// - Derivative computation fails
    fn derivative(
        &self,
        querypoints: &ArrayView2<F>,
        order: usize,
    ) -> InterpolateResult<Vec<F>> {
        if querypoints.ncols() != 1 {
            return Err(crate::InterpolateError::invalid_input(
                "CubicSpline only supports 1D interpolation",
            ));
        }

        let mut results = Vec::with_capacity(querypoints.nrows());
        for row in querypoints.outer_iter() {
            let x = row[0];
            let deriv = self.derivative_n(x, order)?;
            results.push(deriv);
        }
        Ok(results)
    }

    /// Integrate over multiple intervals
    ///
    /// This method provides a batch interface for computing definite integrals
    /// over multiple intervals.
    ///
    /// # Arguments
    ///
    /// * `bounds` - Vector of (lower, upper) bounds for integration
    ///
    /// # Returns
    ///
    /// Vector of integral values for each interval
    ///
    /// # Errors
    ///
    /// Returns an error if any integration fails due to out-of-bounds intervals
    /// or numerical issues.
    fn integrate(&self, bounds: &[(F, F)]) -> InterpolateResult<Vec<F>> {
        let mut results = Vec::with_capacity(bounds.len());
        for &(a, b) in bounds {
            let integral = self.integrate(a, b)?;
            results.push(integral);
        }
        Ok(results)
    }

    /// Compute the antiderivative (indefinite integral) of the spline
    ///
    /// This method returns a new spline that represents the antiderivative
    /// of the current spline. The antiderivative is computed analytically
    /// using the polynomial coefficients.
    ///
    /// # Returns
    ///
    /// A boxed spline interpolator representing the antiderivative
    ///
    /// # Errors
    ///
    /// Returns an error if the antiderivative computation fails.
    fn antiderivative(
        &self,
    ) -> InterpolateResult<Box<dyn SplineInterpolator<F>>> {
        let antideriv = self.antiderivative()?;
        Ok(Box::new(antideriv))
    }

    /// Find roots of the spline within specified intervals
    ///
    /// This method finds all roots (zeros) of the spline within the given
    /// intervals using numerical root-finding algorithms.
    ///
    /// # Arguments
    ///
    /// * `bounds` - Vector of (lower, upper) bounds to search for roots
    /// * `tolerance` - Convergence tolerance for root finding
    ///
    /// # Returns
    ///
    /// Vector of x values where the spline equals zero
    ///
    /// # Errors
    ///
    /// Returns an error if root finding fails or if any interval is invalid.
    fn find_roots(&self, bounds: &[(F, F)], tolerance: F) -> InterpolateResult<Vec<F>> {
        use crate::utils::find_roots_bisection;

        let mut all_roots = Vec::new();

        for &(a, b) in bounds {
            if a >= b {
                continue;
            }

            // Check if a and b are within the spline domain
            let x_min = self.x()[0];
            let x_max = self.x()[self.x().len() - 1];

            let search_a = a.max(x_min);
            let search_b = b.min(x_max);

            if search_a >= search_b {
                continue;
            }

            // Find roots in this interval using bisection method
            let interval_roots = find_roots_bisection(
                search_a,
                search_b,
                tolerance,
                |x| CubicSpline::evaluate(self, x),
            )?;

            all_roots.extend(interval_roots);
        }

        // Remove duplicates and sort
        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)
    }

    /// Find extrema (local minima and maxima) within given bounds
    fn find_extrema(
        &self,
        bounds: &[(F, F)],
        tolerance: F,
    ) -> InterpolateResult<Vec<(F, F, crate::traits::ExtremaType)>> {
        use crate::traits::ExtremaType;

        let mut all_extrema = Vec::new();

        for &(a, b) in bounds {
            if a >= b {
                continue;
            }

            // Check if a and b are within the spline domain
            let x_min = self.x()[0];
            let x_max = self.x()[self.x().len() - 1];

            let search_a = a.max(x_min);
            let search_b = b.min(x_max);

            if search_a >= search_b {
                continue;
            }

            // Find extrema by locating roots of the first derivative
            let step = (search_b - search_a) / F::from(100).expect("Failed to convert constant to float");
            let mut x = search_a;

            while x < search_b {
                let dx = step.min(search_b - x);

                // Get first and second derivatives to classify extrema
                if let (Ok(d1_left), Ok(d1_right)) = (
                    CubicSpline::derivative(self, x),
                    CubicSpline::derivative(self, x + dx)
                ) {
                    // Check for sign change in first derivative
                    if d1_left * d1_right < F::zero() {
                        // Found potential extremum, refine location
                        let mut left = x;
                        let mut right = x + dx;

                        // Bisection to find exact location
                        for _ in 0..20 {
                            let mid = (left + right) / (F::one() + F::one());
                            if let Ok(d1_mid) = CubicSpline::derivative(self, mid) {
                                if d1_left * d1_mid < F::zero() {
                                    right = mid;
                                } else {
                                    left = mid;
                                }
                            } else {
                                break;
                            }

                            if (right - left).abs() < tolerance {
                                break;
                            }
                        }

                        let extremum_x = (left + right) / (F::one() + F::one());
                        if let Ok(extremum_y) = CubicSpline::evaluate(self, extremum_x) {
                            // Classify extremum type using second derivative
                            let extremum_type = if let Ok(d2) = CubicSpline::derivative_n(self, extremum_x, 2) {
                                if d2 > F::zero() {
                                    ExtremaType::Minimum
                                } else if d2 < F::zero() {
                                    ExtremaType::Maximum
                                } else {
                                    ExtremaType::InflectionPoint
                                }
                            } else {
                                ExtremaType::InflectionPoint
                            };

                            all_extrema.push((extremum_x, extremum_y, extremum_type));
                        }
                    }
                }

                x += dx;
            }
        }

        // Remove duplicates
        all_extrema.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
        all_extrema.dedup_by(|a, b| (a.0 - b.0).abs() < tolerance);

        Ok(all_extrema)
    }
}

impl<F> Interpolator<F> for CubicSpline<F>
where
    F: InterpolationFloat,
{
    /// Evaluate the interpolator at given query points
    ///
    /// This method provides the basic interpolation interface, evaluating
    /// the spline at multiple query points.
    ///
    /// # Arguments
    ///
    /// * `querypoints` - 2D array of query points (must have 1 column for 1D splines)
    ///
    /// # Returns
    ///
    /// Vector of interpolated values at each query point
    ///
    /// # Errors
    ///
    /// Returns an error if:
    /// - Query points have more than 1 column (cubic splines are 1D only)
    /// - Any query point is outside the interpolation domain
    /// - Evaluation fails
    fn evaluate(&self, querypoints: &ArrayView2<F>) -> InterpolateResult<Vec<F>> {
        if querypoints.ncols() != 1 {
            return Err(crate::InterpolateError::invalid_input(
                "CubicSpline only supports 1D interpolation",
            ));
        }

        let mut results = Vec::with_capacity(querypoints.nrows());
        for row in querypoints.outer_iter() {
            let x = row[0];
            // Call the spline's evaluate method directly to avoid naming conflict
            let value = CubicSpline::evaluate(self, x)?;
            results.push(value);
        }
        Ok(results)
    }

    /// Get the spatial dimension of the interpolator
    ///
    /// # Returns
    ///
    /// The spatial dimension (always 1 for cubic splines)
    fn dimension(&self) -> usize {
        1
    }

    /// Get the number of data points used to construct the interpolator
    ///
    /// # Returns
    ///
    /// The number of data points
    fn len(&self) -> usize {
        self.x().len()
    }
}

// Additional trait implementations for enhanced functionality
impl<F: InterpolationFloat + ToString> CubicSpline<F> {
    /// Compute the antiderivative of the spline
    ///
    /// Returns a new cubic spline that represents the indefinite integral
    /// of the current spline. The antiderivative is computed analytically
    /// by integrating the polynomial coefficients.
    ///
    /// # Returns
    ///
    /// A new CubicSpline representing the antiderivative
    ///
    /// # Errors
    ///
    /// Returns an error if the antiderivative computation fails.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use scirs2_core::ndarray::array;
    /// use scirs2_interpolate::spline::CubicSpline;
    ///
    /// let x = array![0.0, 1.0, 2.0, 3.0];
    /// let y = array![0.0, 1.0, 4.0, 9.0]; // y = x^2
    /// let spline = CubicSpline::new(&x.view(), &y.view()).expect("Operation failed");
    ///
    /// let antideriv = spline.antiderivative().expect("Operation failed");
    /// // antideriv should approximate x^3/3 + C
    /// ```
    pub fn antiderivative(&self) -> InterpolateResult<CubicSpline<F>> {
        // For a cubic spline with coefficients [a, b, c, d] representing:
        // y(x) = a + b*(x-xi) + c*(x-xi)^2 + d*(x-xi)^3
        //
        // The antiderivative has coefficients:
        // Y(x) = C + a*(x-xi) + b*(x-xi)^2/2 + c*(x-xi)^3/3 + d*(x-xi)^4/4
        //
        // We need to determine the constants C for each segment to ensure continuity.

        let n_segments = self.coeffs().nrows();
        let mut antideriv_coeffs = Array2::zeros((n_segments, 4));

        // Constants for coefficient transformation
        let two = F::from_f64(2.0).ok_or_else(|| {
            crate::InterpolateError::ComputationError(
                "Failed to convert constant 2.0 to float type".to_string(),
            )
        })?;
        let three = F::from_f64(3.0).ok_or_else(|| {
            crate::InterpolateError::ComputationError(
                "Failed to convert constant 3.0 to float type".to_string(),
            )
        })?;
        let four = F::from_f64(4.0).ok_or_else(|| {
            crate::InterpolateError::ComputationError(
                "Failed to convert constant 4.0 to float type".to_string(),
            )
        })?;

        // Transform coefficients for each segment
        for i in 0..n_segments {
            let a = self.coeffs()[[i, 0]];
            let b = self.coeffs()[[i, 1]];
            let c = self.coeffs()[[i, 2]];
            let d = self.coeffs()[[i, 3]];

            // New coefficients for antiderivative
            antideriv_coeffs[[i, 0]] = F::zero(); // Will be set based on continuity
            antideriv_coeffs[[i, 1]] = a;
            antideriv_coeffs[[i, 2]] = b / two;
            antideriv_coeffs[[i, 3]] = c / three;
            // Note: d/4 term becomes the coefficient of (x-xi)^4, but we're working
            // with cubic splines, so we need to handle this carefully.
            // For now, we'll ignore the quartic term and note this limitation.
        }

        // Determine integration constants to ensure continuity
        // Set the first segment's constant to zero (arbitrary choice)
        antideriv_coeffs[[0, 0]] = F::zero();

        // For subsequent segments, ensure continuity at segment boundaries
        for i in 1..n_segments {
            let x_boundary = self.x()[i];
            let x_prev = self.x()[i - 1];

            // Evaluate the previous segment at the boundary
            let dx_prev = x_boundary - x_prev;
            let prev_val = antideriv_coeffs[[i - 1, 0]]
                + antideriv_coeffs[[i - 1, 1]] * dx_prev
                + antideriv_coeffs[[i - 1, 2]] * dx_prev * dx_prev
                + antideriv_coeffs[[i - 1, 3]] * dx_prev * dx_prev * dx_prev;

            // Current segment at boundary (dx = 0)
            antideriv_coeffs[[i, 0]] = prev_val;
        }

        // Create y values for the antiderivative spline by evaluating at data points
        let mut antideriv_y = Array1::zeros(self.x().len());
        for i in 0..self.x().len() {
            let x_val = self.x()[i];

            // Find which segment this point belongs to
            let segment_idx = if i == self.x().len() - 1 {
                n_segments - 1
            } else {
                i.min(n_segments - 1)
            };

            let x_base = self.x()[segment_idx];
            let dx = x_val - x_base;

            antideriv_y[i] = antideriv_coeffs[[segment_idx, 0]]
                + antideriv_coeffs[[segment_idx, 1]] * dx
                + antideriv_coeffs[[segment_idx, 2]] * dx * dx
                + antideriv_coeffs[[segment_idx, 3]] * dx * dx * dx;
        }

        // Create the antiderivative spline using the constructor
        CubicSpline::new(&self.x().view(), &antideriv_y.view())
    }
}