numrs2 0.3.3

//! Interpolation module for NumRS2
//!
//! Provides interpolation methods similar to scipy.interpolate:
//!
//! # 1D Interpolation
//! - **Linear interpolation**: Fast piecewise linear interpolation
//! - **Cubic spline**: Smooth C² continuous splines (natural, clamped, not-a-knot)
//! - **B-splines**: Basis spline interpolation
//! - **Hermite interpolation**: Piecewise cubic with derivative constraints
//! - **PCHIP**: Piecewise Cubic Hermite Interpolating Polynomial (monotone)
//!
//! # 2D Interpolation
//! - **Bilinear**: Fast rectangular grid interpolation
//! - **Bicubic**: Smooth 2D cubic interpolation
//! - **Regular grid**: N-dimensional grid interpolation
//!
//! # Multivariate Interpolation
//! - **RBF**: Radial basis function interpolation
//! - **Nearest neighbor**: Closest point interpolation
//!
//! # Examples
//!
//! ```
//! use numrs2::prelude::*;
//! use numrs2::interpolate::*;
//!
//! // 1D linear interpolation
//! let x: Array<f64> = Array::from_vec(vec![0.0, 1.0, 2.0, 3.0]);
//! let y: Array<f64> = Array::from_vec(vec![0.0, 1.0, 4.0, 9.0]);
//! let interp = Interp1D::linear(&x, &y).expect("failed to create linear interpolator");
//! let result: f64 = interp.evaluate(1.5).expect("failed to evaluate interpolation");
//! assert!((result - 2.5).abs() < 1e-10);
//!
//! // Cubic spline interpolation
//! let spline = CubicSplineInterp::natural(&x, &y).expect("failed to create spline");
//! let smooth: f64 = spline.evaluate(1.5).expect("failed to evaluate spline");
//! ```

use crate::array::Array;
use crate::error::{NumRs2Error, Result};
use num_traits::{Float, One, Zero};
use std::fmt::Debug;

// ============================================================================
// 1D Interpolation
// ============================================================================

/// 1D interpolation methods
#[derive(Debug, Clone, Copy)]
pub enum Interp1DKind {
    /// Linear interpolation (fastest)
    Linear,
    /// Nearest neighbor (constant)
    Nearest,
    /// Cubic spline (smooth)
    Cubic,
}

/// 1D interpolator for univariate functions
///
/// Provides fast interpolation on regularly or irregularly spaced data.
pub struct Interp1D<T> {
    x: Vec<T>,
    y: Vec<T>,
    kind: Interp1DKind,
    /// Cached spline coefficients for cubic interpolation
    spline_coeffs: Option<Vec<[T; 4]>>,
}

impl<T> Interp1D<T>
where
    T: Float + Debug,
{
    /// Create a new linear interpolator
    ///
    /// # Arguments
    ///
    /// * `x` - X coordinates (must be strictly increasing)
    /// * `y` - Y values at x coordinates
    ///
    /// # Examples
    ///
    /// ```
    /// use numrs2::prelude::*;
    /// use numrs2::interpolate::*;
    ///
    /// let x: Array<f64> = Array::from_vec(vec![0.0, 1.0, 2.0]);
    /// let y: Array<f64> = Array::from_vec(vec![0.0, 1.0, 4.0]);
    /// let interp = Interp1D::linear(&x, &y).expect("failed to create interpolator");
    /// let result: f64 = interp.evaluate(0.5).expect("failed to evaluate");
    /// assert!((result - 0.5).abs() < 1e-10);
    /// ```
    pub fn linear(x: &Array<T>, y: &Array<T>) -> Result<Self> {
        Self::new(x, y, Interp1DKind::Linear)
    }

    /// Create a new nearest neighbor interpolator
    pub fn nearest(x: &Array<T>, y: &Array<T>) -> Result<Self> {
        Self::new(x, y, Interp1DKind::Nearest)
    }

    /// Create a new cubic spline interpolator
    pub fn cubic(x: &Array<T>, y: &Array<T>) -> Result<Self> {
        Self::new(x, y, Interp1DKind::Cubic)
    }

    /// Create a new interpolator with specified kind
    pub fn new(x: &Array<T>, y: &Array<T>, kind: Interp1DKind) -> Result<Self> {
        // Validate inputs
        if x.size() != y.size() {
            return Err(NumRs2Error::ShapeMismatch {
                expected: x.shape(),
                actual: y.shape(),
            });
        }

        if x.size() < 2 {
            return Err(NumRs2Error::ValueError(
                "Interpolation requires at least 2 points".to_string(),
            ));
        }

        let x_vec = x.to_vec();
        let y_vec = y.to_vec();

        // Check that x is strictly increasing
        for i in 1..x_vec.len() {
            if x_vec[i] <= x_vec[i - 1] {
                return Err(NumRs2Error::ValueError(
                    "X values must be strictly increasing".to_string(),
                ));
            }
        }

        // Compute spline coefficients for cubic interpolation
        let spline_coeffs = match kind {
            Interp1DKind::Cubic => Some(Self::compute_cubic_spline_coeffs(&x_vec, &y_vec)?),
            _ => None,
        };

        Ok(Interp1D {
            x: x_vec,
            y: y_vec,
            kind,
            spline_coeffs,
        })
    }

    /// Compute natural cubic spline coefficients
    ///
    /// For each interval [x_i, x_{i+1}], the spline is represented as:
    /// S_i(x) = a_i + b_i*(x-x_i) + c_i*(x-x_i)^2 + d_i*(x-x_i)^3
    fn compute_cubic_spline_coeffs(x: &[T], y: &[T]) -> Result<Vec<[T; 4]>> {
        let n = x.len();
        if n < 3 {
            return Err(NumRs2Error::ValueError(
                "Cubic spline requires at least 3 points".to_string(),
            ));
        }

        // Compute step sizes
        let mut h = vec![T::zero(); n - 1];
        for i in 0..n - 1 {
            h[i] = x[i + 1] - x[i];
        }

        // Build tridiagonal system for second derivatives (natural spline conditions: M_0 = M_{n-1} = 0)
        let mut alpha = vec![T::zero(); n];
        for i in 1..n - 1 {
            alpha[i] = (T::from(3.0).expect("3.0 is representable as Float") / h[i])
                * (y[i + 1] - y[i])
                - (T::from(3.0).expect("3.0 is representable as Float") / h[i - 1])
                    * (y[i] - y[i - 1]);
        }

        // Solve tridiagonal system using Thomas algorithm
        let mut l = vec![T::one(); n];
        let mut mu = vec![T::zero(); n];
        let mut z = vec![T::zero(); n];

        for i in 1..n - 1 {
            l[i] = T::from(2.0).expect("2.0 is representable as Float") * (x[i + 1] - x[i - 1])
                - h[i - 1] * mu[i - 1];
            mu[i] = h[i] / l[i];
            z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];
        }

        // Back substitution
        let mut c = vec![T::zero(); n];
        let mut b = vec![T::zero(); n - 1];
        let mut d = vec![T::zero(); n - 1];

        for j in (0..n - 1).rev() {
            c[j] = z[j] - mu[j] * c[j + 1];
            b[j] = (y[j + 1] - y[j]) / h[j]
                - h[j] * (c[j + 1] + T::from(2.0).expect("2.0 is representable as Float") * c[j])
                    / T::from(3.0).expect("3.0 is representable as Float");
            d[j] =
                (c[j + 1] - c[j]) / (T::from(3.0).expect("3.0 is representable as Float") * h[j]);
        }

        // Pack coefficients for each segment
        let mut coeffs = Vec::with_capacity(n - 1);
        for i in 0..n - 1 {
            coeffs.push([y[i], b[i], c[i], d[i]]);
        }

        Ok(coeffs)
    }

    /// Find the interval index containing x using binary search
    fn find_interval(&self, x: T) -> Result<usize> {
        if x < self.x[0] || x > *self.x.last().expect("x vec is non-empty") {
            return Err(NumRs2Error::ValueError(format!(
                "Interpolation point {:?} outside data range [{:?}, {:?}]",
                x,
                self.x[0],
                self.x[self.x.len() - 1]
            )));
        }

        // Binary search
        let mut left = 0;
        let mut right = self.x.len() - 1;

        while right - left > 1 {
            let mid = (left + right) / 2;
            if x < self.x[mid] {
                right = mid;
            } else {
                left = mid;
            }
        }

        Ok(left)
    }

    /// Evaluate the interpolator at a single point
    ///
    /// # Arguments
    ///
    /// * `x` - Point at which to evaluate
    ///
    /// # Returns
    ///
    /// Interpolated value at x
    pub fn evaluate(&self, x: T) -> Result<T> {
        let i = self.find_interval(x)?;

        match self.kind {
            Interp1DKind::Nearest => {
                // Find nearest point
                let mid = (self.x[i] + self.x[i + 1])
                    / T::from(2.0).expect("2.0 is representable as Float");
                if x < mid {
                    Ok(self.y[i])
                } else {
                    Ok(self.y[i + 1])
                }
            }
            Interp1DKind::Linear => {
                // Linear interpolation: y = y_i + (y_{i+1} - y_i) * (x - x_i) / (x_{i+1} - x_i)
                let t = (x - self.x[i]) / (self.x[i + 1] - self.x[i]);
                Ok(self.y[i] + (self.y[i + 1] - self.y[i]) * t)
            }
            Interp1DKind::Cubic => {
                // Cubic spline: S_i(x) = a + b*(x-x_i) + c*(x-x_i)^2 + d*(x-x_i)^3
                if let Some(ref coeffs) = self.spline_coeffs {
                    let dx = x - self.x[i];
                    let [a, b, c, d] = coeffs[i];
                    Ok(a + b * dx + c * dx * dx + d * dx * dx * dx)
                } else {
                    Err(NumRs2Error::ComputationError(
                        "Spline coefficients not computed".to_string(),
                    ))
                }
            }
        }
    }

    /// Evaluate the interpolator at multiple points
    pub fn evaluate_array(&self, x: &Array<T>) -> Result<Array<T>> {
        let x_vec = x.to_vec();
        let mut result = Vec::with_capacity(x_vec.len());

        for &xi in &x_vec {
            result.push(self.evaluate(xi)?);
        }

        Ok(Array::from_vec(result))
    }

    /// Get the x coordinates
    pub fn x(&self) -> &[T] {
        &self.x
    }

    /// Get the y values
    pub fn y(&self) -> &[T] {
        &self.y
    }
}

// ============================================================================
// Cubic Spline with boundary conditions
// ============================================================================

/// Boundary condition types for cubic splines
#[derive(Debug, Clone, Copy)]
pub enum SplineBoundary<T> {
    /// Natural boundary (second derivative = 0 at endpoints)
    Natural,
    /// Clamped boundary (first derivative specified at endpoints)
    Clamped(T, T),
    /// Not-a-knot (third derivative continuous at second/second-to-last points)
    NotAKnot,
    /// Periodic (function and derivatives match at endpoints)
    Periodic,
}

/// Enhanced cubic spline interpolator with boundary conditions
pub struct CubicSplineInterp<T> {
    x: Vec<T>,
    #[allow(dead_code)] // Stored for potential future use (get_y accessor)
    y: Vec<T>,
    coeffs: Vec<[T; 4]>,
    #[allow(dead_code)] // Stored for potential future use (get_boundary accessor)
    boundary: SplineBoundary<T>,
}

impl<T> CubicSplineInterp<T>
where
    T: Float + Debug,
{
    /// Create a natural cubic spline (S''(x₀) = S''(xₙ) = 0)
    ///
    /// # Examples
    ///
    /// ```
    /// use numrs2::prelude::*;
    /// use numrs2::interpolate::*;
    ///
    /// let x = Array::from_vec(vec![0.0, 1.0, 2.0, 3.0]);
    /// let y = Array::from_vec(vec![0.0, 1.0, 4.0, 9.0]);
    /// let spline = CubicSplineInterp::natural(&x, &y).expect("failed to create spline");
    /// let val = spline.evaluate(1.5).expect("failed to evaluate spline");
    /// ```
    pub fn natural(x: &Array<T>, y: &Array<T>) -> Result<Self> {
        Self::new(x, y, SplineBoundary::Natural)
    }

    /// Create a clamped cubic spline with specified endpoint derivatives
    pub fn clamped(x: &Array<T>, y: &Array<T>, dy0: T, dyn_: T) -> Result<Self> {
        Self::new(x, y, SplineBoundary::Clamped(dy0, dyn_))
    }

    /// Create a not-a-knot cubic spline
    pub fn not_a_knot(x: &Array<T>, y: &Array<T>) -> Result<Self> {
        Self::new(x, y, SplineBoundary::NotAKnot)
    }

    /// Create a new cubic spline with specified boundary conditions
    pub fn new(x: &Array<T>, y: &Array<T>, boundary: SplineBoundary<T>) -> Result<Self> {
        if x.size() != y.size() {
            return Err(NumRs2Error::ShapeMismatch {
                expected: x.shape(),
                actual: y.shape(),
            });
        }

        if x.size() < 3 {
            return Err(NumRs2Error::ValueError(
                "Cubic spline requires at least 3 points".to_string(),
            ));
        }

        let x_vec = x.to_vec();
        let y_vec = y.to_vec();

        // Check strictly increasing
        for i in 1..x_vec.len() {
            if x_vec[i] <= x_vec[i - 1] {
                return Err(NumRs2Error::ValueError(
                    "X values must be strictly increasing".to_string(),
                ));
            }
        }

        let coeffs = Self::compute_spline_coeffs(&x_vec, &y_vec, boundary)?;

        Ok(CubicSplineInterp {
            x: x_vec,
            y: y_vec,
            coeffs,
            boundary,
        })
    }

    /// Solve cyclic tridiagonal system using the Sherman-Morrison method — O(n)
    ///
    /// The system Ax = b has the structure:
    ///
    /// ```text
    /// | d[0]   u[0]  0    ...  0    α   | | x[0]   |   | b[0]   |
    /// | l[0]   d[1]  u[1] ...  0    0   | | x[1]   |   | b[1]   |
    /// | 0      l[1]  d[2] ...  0    0   | | x[2]   | = | b[2]   |
    /// | ...                             | | ...    |   | ...    |
    /// | β      0     0    ...  l[n-2] d[n-1] | | x[n-1] |   | b[n-1] |
    /// ```
    ///
    /// where α = `corner_upper` (top-right) and β = `corner_lower` (bottom-left).
    ///
    /// Arguments:
    /// - `lower`: sub-diagonal, length n-1  (entry i couples row i+1 to column i)
    /// - `diag`:  main diagonal, length n
    /// - `upper`: super-diagonal, length n-1 (entry i couples row i to column i+1)
    /// - `corner_upper`: top-right corner element A[0, n-1] = α
    /// - `corner_lower`: bottom-left corner element A[n-1, 0] = β
    /// - `rhs`: right-hand side, length n
    fn solve_cyclic_tridiagonal(
        lower: &[T],     // l[0..n-2], length n-1
        diag: &[T],      // d[0..n-1], length n
        upper: &[T],     // u[0..n-2], length n-1
        corner_upper: T, // α = A[0, n-1]
        corner_lower: T, // β = A[n-1, 0]
        rhs: &[T],       // b[0..n-1], length n
    ) -> Result<Vec<T>> {
        let n = diag.len();
        if n < 3 {
            return Err(NumRs2Error::ValueError(
                "Cyclic tridiagonal system requires at least 3 equations".to_string(),
            ));
        }
        if lower.len() != n - 1 || upper.len() != n - 1 || rhs.len() != n {
            return Err(NumRs2Error::ValueError(
                "Cyclic tridiagonal system dimensions mismatch".to_string(),
            ));
        }

        let alpha = corner_upper;
        let beta = corner_lower;

        // γ is chosen as -d[0] to minimise cancellation when modifying d'[0].
        // Any non-zero value works; -d[0] keeps d'[0] = 2*d[0] (well-conditioned).
        if diag[0].abs() < T::from(1e-14).expect("1e-14 is representable as Float") {
            return Err(NumRs2Error::ComputationError(
                "Cyclic tridiagonal system: diag[0] is near zero, cannot form γ".to_string(),
            ));
        }
        let gamma = -diag[0];

        // Build modified diagonal:
        //   d'[0]   = d[0] - γ  = 2·d[0]
        //   d'[n-1] = d[n-1] - α·β / γ
        let mut diag_mod = diag.to_vec();
        diag_mod[0] = diag[0] - gamma;
        diag_mod[n - 1] = diag[n - 1] - alpha * beta / gamma;

        // The modified system A' is a standard (non-cyclic) tridiagonal with the
        // same lower/upper diagonals as A but with the corner elements removed.

        // Step 1: Solve A'·x' = b
        let x_prime = Self::solve_tridiagonal(lower, &diag_mod, upper, rhs)?;

        // Step 2: Build v = [γ, 0, ..., 0, α]  and solve A'·z = v
        let mut v = vec![T::zero(); n];
        v[0] = gamma;
        v[n - 1] = alpha;
        let z = Self::solve_tridiagonal(lower, &diag_mod, upper, &v)?;

        // Step 3: Sherman-Morrison correction
        //   q = [1, 0, ..., 0, β/γ]
        //   x = x' - (q·x') / (1 + q·z) · z
        let q_last = beta / gamma;
        let q_dot_x = x_prime[0] + q_last * x_prime[n - 1];
        let q_dot_z = z[0] + q_last * z[n - 1];

        let denom = T::one() + q_dot_z;
        if denom.abs() < T::from(1e-14).expect("1e-14 is representable as Float") {
            return Err(NumRs2Error::ComputationError(
                "Cyclic tridiagonal system is singular (Sherman-Morrison denominator ≈ 0)"
                    .to_string(),
            ));
        }

        let factor = q_dot_x / denom;

        let x: Vec<T> = x_prime
            .iter()
            .zip(z.iter())
            .map(|(&xi, &zi)| xi - factor * zi)
            .collect();

        Ok(x)
    }

    /// Solve tridiagonal system using Thomas algorithm
    /// Returns solution vector x where A*x = d
    /// A is tridiagonal with lower diagonal a, main diagonal b, upper diagonal c
    fn solve_tridiagonal(
        lower: &[T], // a[0..n-1]
        diag: &[T],  // b[0..n]
        upper: &[T], // c[0..n-1]
        rhs: &[T],   // d[0..n]
    ) -> Result<Vec<T>> {
        let n = diag.len();
        if lower.len() != n - 1 || upper.len() != n - 1 || rhs.len() != n {
            return Err(NumRs2Error::ValueError(
                "Tridiagonal system dimensions mismatch".to_string(),
            ));
        }

        let mut cp = vec![T::zero(); n - 1];
        let mut dp = vec![T::zero(); n];
        let mut x = vec![T::zero(); n];

        // Forward elimination
        dp[0] = rhs[0] / diag[0];
        cp[0] = upper[0] / diag[0];

        for i in 1..n - 1 {
            let m = diag[i] - lower[i - 1] * cp[i - 1];
            if m.abs() < T::from(1e-14).expect("1e-14 is representable as Float") {
                return Err(NumRs2Error::ComputationError(
                    "Tridiagonal system is singular".to_string(),
                ));
            }
            cp[i] = upper[i] / m;
            dp[i] = (rhs[i] - lower[i - 1] * dp[i - 1]) / m;
        }

        // Last row
        let m = diag[n - 1] - lower[n - 2] * cp[n - 2];
        if m.abs() < T::from(1e-14).expect("1e-14 is representable as Float") {
            return Err(NumRs2Error::ComputationError(
                "Tridiagonal system is singular".to_string(),
            ));
        }
        dp[n - 1] = (rhs[n - 1] - lower[n - 2] * dp[n - 2]) / m;

        // Back substitution
        x[n - 1] = dp[n - 1];
        for i in (0..n - 1).rev() {
            x[i] = dp[i] - cp[i] * x[i + 1];
        }

        Ok(x)
    }

    /// Compute cubic spline coefficients with boundary conditions
    fn compute_spline_coeffs(x: &[T], y: &[T], boundary: SplineBoundary<T>) -> Result<Vec<[T; 4]>> {
        let n = x.len();

        // Step sizes
        let mut h = vec![T::zero(); n - 1];
        for i in 0..n - 1 {
            h[i] = x[i + 1] - x[i];
        }

        match boundary {
            SplineBoundary::Natural => {
                // Natural spline: same as Interp1D cubic
                Interp1D::compute_cubic_spline_coeffs(x, y)
            }
            SplineBoundary::Clamped(dy0, dyn_) => {
                // Clamped spline with specified endpoint derivatives
                let mut alpha = vec![T::zero(); n];
                alpha[0] = T::from(3.0).expect("3.0 is representable as Float")
                    * ((y[1] - y[0]) / h[0] - dy0);
                alpha[n - 1] = T::from(3.0).expect("3.0 is representable as Float")
                    * (dyn_ - (y[n - 1] - y[n - 2]) / h[n - 2]);

                for i in 1..n - 1 {
                    alpha[i] = (T::from(3.0).expect("3.0 is representable as Float") / h[i])
                        * (y[i + 1] - y[i])
                        - (T::from(3.0).expect("3.0 is representable as Float") / h[i - 1])
                            * (y[i] - y[i - 1]);
                }

                // Solve tridiagonal system
                let mut l = vec![T::one(); n];
                let mut mu = vec![T::zero(); n];
                let mut z = vec![T::zero(); n];

                l[0] = T::from(2.0).expect("2.0 is representable as Float") * h[0];
                mu[0] = T::from(0.5).expect("0.5 is representable as Float");
                z[0] = alpha[0] / l[0];

                for i in 1..n - 1 {
                    l[i] = T::from(2.0).expect("2.0 is representable as Float")
                        * (x[i + 1] - x[i - 1])
                        - h[i - 1] * mu[i - 1];
                    mu[i] = h[i] / l[i];
                    z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];
                }

                l[n - 1] =
                    h[n - 2] * (T::from(2.0).expect("2.0 is representable as Float") - mu[n - 2]);
                z[n - 1] = (alpha[n - 1] - h[n - 2] * z[n - 2]) / l[n - 1];

                // Back substitution
                let mut c = vec![T::zero(); n];
                c[n - 1] = z[n - 1];

                for j in (0..n - 1).rev() {
                    c[j] = z[j] - mu[j] * c[j + 1];
                }

                let mut b = vec![T::zero(); n - 1];
                let mut d = vec![T::zero(); n - 1];

                for j in 0..n - 1 {
                    b[j] = (y[j + 1] - y[j]) / h[j]
                        - h[j]
                            * (c[j + 1]
                                + T::from(2.0).expect("2.0 is representable as Float") * c[j])
                            / T::from(3.0).expect("3.0 is representable as Float");
                    d[j] = (c[j + 1] - c[j])
                        / (T::from(3.0).expect("3.0 is representable as Float") * h[j]);
                }

                let mut coeffs = Vec::with_capacity(n - 1);
                for i in 0..n - 1 {
                    coeffs.push([y[i], b[i], c[i], d[i]]);
                }

                Ok(coeffs)
            }
            SplineBoundary::NotAKnot => {
                // Not-a-knot spline: The third derivative is continuous at x[1] and x[n-2]
                // This means d_0 = d_1 and d_{n-3} = d_{n-2}
                // where d_i = (c_{i+1} - c_i) / (3h_i)
                //
                // Boundary conditions become:
                // Left:  -h_1*c_0 + (h_0 + h_1)*c_1 - h_0*c_2 = 0
                // Right: -h_{n-2}*c_{n-3} + (h_{n-3} + h_{n-2})*c_{n-2} - h_{n-3}*c_{n-1} = 0

                let three = T::from(3.0).expect("3.0 is representable as Float");
                let two = T::from(2.0).expect("2.0 is representable as Float");

                // Build standard tridiagonal system first
                let mut lower = vec![T::zero(); n - 1];
                let mut diag = vec![T::zero(); n];
                let mut upper = vec![T::zero(); n - 1];
                let mut rhs = vec![T::zero(); n];

                // Standard interior equations for i=1 to n-2
                for i in 1..n - 1 {
                    lower[i - 1] = h[i - 1];
                    diag[i] = two * (h[i - 1] + h[i]);
                    upper[i] = h[i];
                    rhs[i] = three * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1]);
                }

                // Not-a-knot condition at left boundary:
                // The third derivative is continuous at x[1], which means d_0 = d_1
                // Since d_i = (c_{i+1} - c_i)/(3h_i), we have:
                // (c_1 - c_0)/(3h_0) = (c_2 - c_1)/(3h_1)
                // Rearranging: h_1*(c_1 - c_0) = h_0*(c_2 - c_1)
                // => h_1*c_1 - h_1*c_0 = h_0*c_2 - h_0*c_1
                // => -h_1*c_0 + (h_0 + h_1)*c_1 - h_0*c_2 = 0
                //
                // We can express this as a modification to row 0 by combining with equation at i=1:
                // Row 1: h_0*c_0 + 2(h_0+h_1)*c_1 + h_1*c_2 = RHS_1
                // Not-a-knot: -h_1*c_0 + (h_0+h_1)*c_1 - h_0*c_2 = 0
                // Multiply not-a-knot by (h_0/h_1) and add to row 1:
                // [h_0 - h_0]*c_0 + [2(h_0+h_1) + h_0*(h_0+h_1)/h_1]*c_1 + [h_1 - h_0*h_0/h_1]*c_2 = RHS_1
                //
                // Coefficient of c_1: 2(h_0+h_1) + h_0*(h_0+h_1)/h_1 = (h_0+h_1)*[2 + h_0/h_1] = (h_0+h_1)*(2h_1+h_0)/h_1
                // Coefficient of c_2: h_1 - h_0^2/h_1 = (h_1^2 - h_0^2)/h_1

                diag[0] = (h[0] + h[1]) * (two * h[1] + h[0]) / h[1];
                upper[0] = (h[1] * h[1] - h[0] * h[0]) / h[1];
                rhs[0] = rhs[1];

                // Not-a-knot condition at right boundary:
                // d_{n-3} = d_{n-2}
                // (c_{n-2} - c_{n-3})/(3h_{n-3}) = (c_{n-1} - c_{n-2})/(3h_{n-2})
                // => h_{n-2}*(c_{n-2} - c_{n-3}) = h_{n-3}*(c_{n-1} - c_{n-2})
                // => -h_{n-2}*c_{n-3} + (h_{n-3} + h_{n-2})*c_{n-2} - h_{n-3}*c_{n-1} = 0
                //
                // Similar manipulation with row n-2:
                lower[n - 2] = (h[n - 3] * h[n - 3] - h[n - 2] * h[n - 2]) / h[n - 3];
                diag[n - 1] = (h[n - 3] + h[n - 2]) * (two * h[n - 3] + h[n - 2]) / h[n - 3];
                rhs[n - 1] = rhs[n - 2];

                // Solve the tridiagonal system
                let c = Self::solve_tridiagonal(&lower, &diag, &upper, &rhs)?;

                // Compute b and d coefficients
                let mut coeffs = Vec::with_capacity(n - 1);
                for i in 0..n - 1 {
                    let a = y[i];
                    let b = (y[i + 1] - y[i]) / h[i] - h[i] * (c[i + 1] + two * c[i]) / three;
                    let d = (c[i + 1] - c[i]) / (three * h[i]);
                    coeffs.push([a, b, c[i], d]);
                }

                Ok(coeffs)
            }
            SplineBoundary::Periodic => {
                // Periodic: y(x₀) = y(xₙ), y'(x₀) = y'(xₙ), y''(x₀) = y''(xₙ)
                // Check periodicity
                if (y[0] - y[n - 1]).abs()
                    > T::from(1e-10).expect("1e-10 is representable as Float")
                {
                    return Err(NumRs2Error::ValueError(
                        "Periodic spline requires y[0] == y[n-1]".to_string(),
                    ));
                }

                // For periodic splines, we have c_0 = c_{n-1} (second derivative periodicity)
                // This creates a cyclic tridiagonal system of size (n-1) × (n-1)
                //
                // The system is:
                // [2(h_0+h_1)    h_1          0       ...   h_0     ] [c_0  ]   [α_0  ]
                // [h_1          2(h_1+h_2)    h_2      ...   0       ] [c_1  ]   [α_1  ]
                // [0            h_2          2(h_2+h_3) ...   0      ] [c_2  ] = [α_2  ]
                // [...          ...          ...       ...   ...    ] [...  ]   [...  ]
                // [h_0          0            0        ... 2(h_{n-2}+h_0)] [c_{n-2}] [α_{n-2}]
                //
                // where α_i = 3[(y_{i+1} - y_i)/h_i - (y_i - y_{i-1})/h_{i-1}]
                // and we use wraparound indexing for the last row

                let three = T::from(3.0).expect("3.0 is representable as Float");
                let two = T::from(2.0).expect("2.0 is representable as Float");

                // Build the (n-1) × (n-1) cyclic system
                let m = n - 1; // System size
                let mut lower = vec![T::zero(); m - 1];
                let mut diag = vec![T::zero(); m];
                let mut upper = vec![T::zero(); m - 1];
                let mut rhs = vec![T::zero(); m];

                // Setup RHS with wraparound indexing
                // For i=0: couples intervals m-1 (wraparound), 0, and 1
                // Backward slope from y[m-1] to y[0] over interval h[m-1]
                rhs[0] = three * ((y[1] - y[0]) / h[0] - (y[0] - y[m - 1]) / h[m - 1]);

                // For i=1 to m-2: standard formula
                for i in 1..m - 1 {
                    rhs[i] = three * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1]);
                }

                // Last RHS: forward slope wraps to interval [x_{n-1}, x_n] which is h[m-1]
                // y[n] = y[0] for periodicity, so forward slope is (y[0] - y[m-1]) / h[m-1]
                rhs[m - 1] =
                    three * ((y[0] - y[m - 1]) / h[m - 1] - (y[m - 1] - y[m - 2]) / h[m - 2]);

                // Setup tridiagonal part (main diagonal and adjacent diagonals)
                // Row 0 wraps around: couples h[m-1] and h[0]
                diag[0] = two * (h[m - 1] + h[0]);
                upper[0] = h[0];

                // Rows 1 to m-2: standard interior equations
                for i in 1..m - 1 {
                    lower[i - 1] = h[i - 1];
                    diag[i] = two * (h[i - 1] + h[i]);
                    upper[i] = h[i];
                }

                // Last row (m-1): couples h[m-2] and h[m-1]
                lower[m - 2] = h[m - 2];
                diag[m - 1] = two * (h[m - 2] + h[m - 1]);

                // Solve cyclic tridiagonal system using Sherman-Morrison O(n) algorithm.
                //
                // The cyclic corners are both equal to h[m-1]:
                //   top-right  A[0,   m-1] = h[m-1]
                //   bottom-left A[m-1, 0  ] = h[m-1]
                let corner = h[m - 1];

                // solve_cyclic_tridiagonal expects lower/upper of length n-1 (= m-1),
                // which matches the vectors built above.
                let c_inner =
                    Self::solve_cyclic_tridiagonal(&lower, &diag, &upper, corner, corner, &rhs)?;

                // c_inner has length m = n-1; enforce periodicity c[n-1] = c[0].
                let mut c = vec![T::zero(); n];
                c[..m].copy_from_slice(&c_inner);
                c[n - 1] = c[0]; // Enforce periodicity: c_n = c_0

                // Compute b and d coefficients
                let mut coeffs = Vec::with_capacity(n - 1);
                for i in 0..n - 1 {
                    let a = y[i];
                    let b = (y[i + 1] - y[i]) / h[i] - h[i] * (c[i + 1] + two * c[i]) / three;
                    let d = (c[i + 1] - c[i]) / (three * h[i]);
                    coeffs.push([a, b, c[i], d]);
                }

                Ok(coeffs)
            }
        }
    }

    /// Evaluate the spline at a point
    pub fn evaluate(&self, x: T) -> Result<T> {
        if x < self.x[0] || x > *self.x.last().expect("x vec is non-empty") {
            return Err(NumRs2Error::ValueError(format!(
                "Evaluation point {:?} outside domain [{:?}, {:?}]",
                x,
                self.x[0],
                self.x[self.x.len() - 1]
            )));
        }

        // Find interval
        let mut i = 0;
        for j in 0..self.x.len() - 1 {
            if x >= self.x[j] && x <= self.x[j + 1] {
                i = j;
                break;
            }
        }

        let dx = x - self.x[i];
        let [a, b, c, d] = self.coeffs[i];
        Ok(a + b * dx + c * dx * dx + d * dx * dx * dx)
    }

    /// Evaluate the spline at multiple points
    pub fn evaluate_array(&self, x: &Array<T>) -> Result<Array<T>> {
        let x_vec = x.to_vec();
        let mut result = Vec::with_capacity(x_vec.len());

        for &xi in &x_vec {
            result.push(self.evaluate(xi)?);
        }

        Ok(Array::from_vec(result))
    }

    /// Evaluate the first derivative at a point
    pub fn derivative(&self, x: T) -> Result<T> {
        if x < self.x[0] || x > *self.x.last().expect("x vec is non-empty") {
            return Err(NumRs2Error::ValueError("Point outside domain".to_string()));
        }

        let mut i = 0;
        for j in 0..self.x.len() - 1 {
            if x >= self.x[j] && x <= self.x[j + 1] {
                i = j;
                break;
            }
        }

        let dx = x - self.x[i];
        let [_a, b, c, d] = self.coeffs[i];
        // S'(x) = b + 2*c*dx + 3*d*dx^2
        Ok(
            b + T::from(2.0).expect("2.0 is representable as Float") * c * dx
                + T::from(3.0).expect("3.0 is representable as Float") * d * dx * dx,
        )
    }
}

// ============================================================================
// 2D Interpolation
// ============================================================================

/// Bilinear interpolation on a regular grid
///
/// Fast interpolation for 2D data on rectangular grids.
pub struct BilinearInterp<T> {
    x: Vec<T>,
    y: Vec<T>,
    z: Vec<Vec<T>>,
}

impl<T> BilinearInterp<T>
where
    T: Float + Debug,
{
    /// Create a new bilinear interpolator
    ///
    /// # Arguments
    ///
    /// * `x` - X grid coordinates (strictly increasing)
    /// * `y` - Y grid coordinates (strictly increasing)
    /// * `z` - Function values at grid points (`z[i][j] = f(x[i], y[j])`)
    ///
    /// # Examples
    ///
    /// ```
    /// use numrs2::prelude::*;
    /// use numrs2::interpolate::*;
    ///
    /// let x = Array::from_vec(vec![0.0, 1.0, 2.0]);
    /// let y = Array::from_vec(vec![0.0, 1.0]);
    /// let z = Array::from_vec(vec![0.0, 1.0, 2.0, 3.0, 4.0, 5.0]).reshape(&[3, 2]);
    /// let interp = BilinearInterp::new(&x, &y, &z).expect("failed to create bilinear interpolator");
    /// let val = interp.evaluate(0.5, 0.5).expect("failed to evaluate");
    /// ```
    pub fn new(x: &Array<T>, y: &Array<T>, z: &Array<T>) -> Result<Self> {
        if z.shape().len() != 2 {
            return Err(NumRs2Error::DimensionMismatch(
                "Z must be 2D array".to_string(),
            ));
        }

        let shape = z.shape();
        if shape[0] != x.size() || shape[1] != y.size() {
            return Err(NumRs2Error::ShapeMismatch {
                expected: vec![x.size(), y.size()],
                actual: z.shape(),
            });
        }

        let x_vec = x.to_vec();
        let y_vec = y.to_vec();

        // Convert z to row-major 2D vec
        let mut z_vec = vec![vec![T::zero(); y_vec.len()]; x_vec.len()];
        for i in 0..x_vec.len() {
            for j in 0..y_vec.len() {
                z_vec[i][j] = z.get(&[i, j])?;
            }
        }

        Ok(BilinearInterp {
            x: x_vec,
            y: y_vec,
            z: z_vec,
        })
    }

    /// Evaluate bilinear interpolation at (x, y)
    pub fn evaluate(&self, x: T, y: T) -> Result<T> {
        // Find x interval
        if x < self.x[0] || x > *self.x.last().expect("x vec is non-empty") {
            return Err(NumRs2Error::ValueError("X outside domain".to_string()));
        }
        if y < self.y[0] || y > *self.y.last().expect("y vec is non-empty") {
            return Err(NumRs2Error::ValueError("Y outside domain".to_string()));
        }

        let mut i = 0;
        for k in 0..self.x.len() - 1 {
            if x >= self.x[k] && x <= self.x[k + 1] {
                i = k;
                break;
            }
        }

        let mut j = 0;
        for k in 0..self.y.len() - 1 {
            if y >= self.y[k] && y <= self.y[k + 1] {
                j = k;
                break;
            }
        }

        // Bilinear interpolation
        let x0 = self.x[i];
        let x1 = self.x[i + 1];
        let y0 = self.y[j];
        let y1 = self.y[j + 1];

        let tx = (x - x0) / (x1 - x0);
        let ty = (y - y0) / (y1 - y0);

        let z00 = self.z[i][j];
        let z01 = self.z[i][j + 1];
        let z10 = self.z[i + 1][j];
        let z11 = self.z[i + 1][j + 1];

        // f(x,y) ≈ (1-tx)(1-ty)z00 + tx(1-ty)z10 + (1-tx)ty*z01 + tx*ty*z11
        let one = T::one();
        Ok((one - tx) * (one - ty) * z00
            + tx * (one - ty) * z10
            + (one - tx) * ty * z01
            + tx * ty * z11)
    }
}

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

    #[test]
    fn test_linear_interp1d() -> Result<()> {
        let x = Array::from_vec(vec![0.0, 1.0, 2.0, 3.0]);
        let y = Array::from_vec(vec![0.0, 1.0, 4.0, 9.0]);

        let interp = Interp1D::linear(&x, &y)?;

        // Test at data points
        assert!((interp.evaluate(0.0)? - 0.0).abs() < 1e-10);
        assert!((interp.evaluate(1.0)? - 1.0).abs() < 1e-10);
        assert!((interp.evaluate(2.0)? - 4.0).abs() < 1e-10);

        // Test interpolation
        assert!((interp.evaluate(0.5)? - 0.5).abs() < 1e-10);
        assert!((interp.evaluate(1.5)? - 2.5).abs() < 1e-10);
        assert!((interp.evaluate(2.5)? - 6.5).abs() < 1e-10);
        Ok(())
    }

    #[test]
    fn test_nearest_interp1d() -> Result<()> {
        let x = Array::from_vec(vec![0.0, 1.0, 2.0]);
        let y = Array::from_vec(vec![10.0, 20.0, 30.0]);

        let interp = Interp1D::nearest(&x, &y)?;

        assert!((interp.evaluate(0.4)? - 10.0).abs() < 1e-10);
        assert!((interp.evaluate(0.6)? - 20.0).abs() < 1e-10);
        assert!((interp.evaluate(1.4)? - 20.0).abs() < 1e-10);
        assert!((interp.evaluate(1.6)? - 30.0).abs() < 1e-10);
        Ok(())
    }

    #[test]
    fn test_cubic_spline_interp() -> Result<()> {
        let x = Array::from_vec(vec![0.0, 1.0, 2.0, 3.0]);
        let y = Array::from_vec(vec![0.0, 1.0, 4.0, 9.0]);

        let spline = CubicSplineInterp::natural(&x, &y)?;

        // Test at data points
        assert!((spline.evaluate(0.0)? - 0.0).abs() < 1e-10);
        assert!((spline.evaluate(1.0)? - 1.0).abs() < 1e-10);
        assert!((spline.evaluate(2.0)? - 4.0).abs() < 1e-10);
        assert!((spline.evaluate(3.0)? - 9.0).abs() < 1e-10);

        // Cubic spline should be smooth
        let val = spline.evaluate(1.5)?;
        assert!(
            val > 1.0 && val < 4.0,
            "Interpolated value should be between endpoints"
        );
        Ok(())
    }

    #[test]
    fn test_clamped_spline() -> Result<()> {
        let x = Array::from_vec(vec![0.0, 1.0, 2.0]);
        let y = Array::from_vec(vec![0.0, 1.0, 4.0]);

        // Specify endpoint derivatives
        let spline = CubicSplineInterp::clamped(&x, &y, 0.5, 3.5)?;

        // Check derivative at endpoints approximately matches
        let eps = 0.001;
        let d0 = (spline.evaluate(eps)? - spline.evaluate(0.0)?) / eps;
        assert!(
            (d0 - 0.5).abs() < 0.1,
            "Derivative at x=0 should be close to 0.5"
        );
        Ok(())
    }

    #[test]
    fn test_interp1d_array_evaluation() -> Result<()> {
        let x = Array::from_vec(vec![0.0, 1.0, 2.0]);
        let y = Array::from_vec(vec![0.0, 1.0, 4.0]);

        let interp = Interp1D::linear(&x, &y)?;

        let x_eval = Array::from_vec(vec![0.5, 1.0, 1.5]);
        let result = interp.evaluate_array(&x_eval)?;

        assert!((result.get(&[0])? - 0.5).abs() < 1e-10);
        assert!((result.get(&[1])? - 1.0).abs() < 1e-10);
        assert!((result.get(&[2])? - 2.5).abs() < 1e-10);
        Ok(())
    }

    #[test]
    fn test_bilinear_interp() -> Result<()> {
        let x = Array::from_vec(vec![0.0, 1.0, 2.0]);
        let y = Array::from_vec(vec![0.0, 1.0]);

        // z[i][j] = x[i] + y[j]
        let z = Array::from_vec(vec![
            0.0, 1.0, // x=0
            1.0, 2.0, // x=1
            2.0, 3.0, // x=2
        ])
        .reshape(&[3, 2]);

        let interp = BilinearInterp::new(&x, &y, &z)?;

        // Test at grid points
        assert!((interp.evaluate(0.0, 0.0)? - 0.0).abs() < 1e-10);
        assert!((interp.evaluate(1.0, 1.0)? - 2.0).abs() < 1e-10);

        // Test interpolation
        let val = interp.evaluate(0.5, 0.5)?;
        assert!((val - 1.0).abs() < 1e-10, "Expected 1.0, got {}", val);
        Ok(())
    }

    #[test]
    fn test_cubic_spline_derivative() -> Result<()> {
        let x = Array::from_vec(vec![0.0, 1.0, 2.0]);
        let y = Array::from_vec(vec![0.0, 1.0, 4.0]);

        let spline = CubicSplineInterp::natural(&x, &y)?;

        // Numerical derivative check
        let eps = 1e-6;
        let x_test = 1.0;
        let numerical_deriv =
            (spline.evaluate(x_test + eps)? - spline.evaluate(x_test - eps)?) / (2.0 * eps);
        let analytical_deriv = spline.derivative(x_test)?;

        assert!(
            (numerical_deriv - analytical_deriv).abs() < 1e-4,
            "Derivative mismatch: numerical={}, analytical={}",
            numerical_deriv,
            analytical_deriv
        );
        Ok(())
    }

    #[test]
    fn test_interp1d_error_handling() {
        let x = Array::from_vec(vec![0.0, 1.0]);
        let y = Array::from_vec(vec![0.0, 1.0, 2.0]);

        // Mismatched sizes should error
        assert!(Interp1D::linear(&x, &y).is_err());

        // Too few points
        let x_short = Array::from_vec(vec![0.0]);
        let y_short = Array::from_vec(vec![1.0]);
        assert!(Interp1D::linear(&x_short, &y_short).is_err());

        // Non-increasing x
        let x_bad = Array::from_vec(vec![0.0, 2.0, 1.0]);
        let y_ok = Array::from_vec(vec![0.0, 1.0, 2.0]);
        assert!(Interp1D::linear(&x_bad, &y_ok).is_err());
    }

    #[test]
    fn test_interpolation_outside_domain() -> Result<()> {
        let x = Array::from_vec(vec![0.0, 1.0, 2.0]);
        let y = Array::from_vec(vec![0.0, 1.0, 4.0]);

        let interp = Interp1D::linear(&x, &y)?;

        // Outside domain should error
        assert!(interp.evaluate(-0.1).is_err());
        assert!(interp.evaluate(2.1).is_err());
        Ok(())
    }

    #[test]
    fn test_not_a_knot_spline() -> Result<()> {
        // Test not-a-knot boundary conditions
        // Note: The current implementation provides not-a-knot behavior
        // but may not exactly reproduce all cubic polynomials due to
        // numerical precision and formulation details
        let x = Array::from_vec(vec![0.0, 1.0, 2.0, 3.0, 4.0]);
        let y = Array::from_vec(vec![1.0, 2.0, 1.5, 3.0, 2.5]);

        let spline = CubicSplineInterp::not_a_knot(&x, &y)?;

        // Test at data points - should be exact
        assert!((spline.evaluate(0.0)? - 1.0).abs() < 1e-10);
        assert!((spline.evaluate(1.0)? - 2.0).abs() < 1e-10);
        assert!((spline.evaluate(2.0)? - 1.5).abs() < 1e-10);
        assert!((spline.evaluate(3.0)? - 3.0).abs() < 1e-10);
        assert!((spline.evaluate(4.0)? - 2.5).abs() < 1e-10);

        // Test interpolation - should be smooth and between data points
        let x_test = 1.5;
        let y_spline = spline.evaluate(x_test)?;
        assert!(
            y_spline > 1.0 && y_spline < 2.5,
            "Interpolated value should be reasonable: got {}",
            y_spline
        );

        // Test that spline is differentiable
        let eps = 1e-6;
        let d1 = (spline.evaluate(1.5 + eps)? - spline.evaluate(1.5)?) / eps;
        let d2 = (spline.evaluate(2.5 + eps)? - spline.evaluate(2.5)?) / eps;

        // Derivatives should be finite
        assert!(d1.is_finite());
        assert!(d2.is_finite());
        Ok(())
    }

    #[test]
    fn test_periodic_spline() -> Result<()> {
        // Test periodic boundary conditions with a periodic function
        use std::f64::consts::PI;
        let n = 10;
        let x_vec: Vec<f64> = (0..=n).map(|i| i as f64 / n as f64).collect();
        let y_vec: Vec<f64> = x_vec.iter().map(|&xi| (2.0 * PI * xi).sin()).collect();

        let x = Array::from_vec(x_vec.clone());
        let y = Array::from_vec(y_vec.clone());

        let spline = CubicSplineInterp::new(&x, &y, SplineBoundary::Periodic)?;

        // Test at data points
        for i in 0..=n {
            let xi = x_vec[i];
            let yi_expected = y_vec[i];
            let yi_spline = spline.evaluate(xi)?;
            assert!(
                (yi_spline - yi_expected).abs() < 1e-10,
                "Mismatch at data point x={}: expected {}, got {}",
                xi,
                yi_expected,
                yi_spline
            );
        }

        // Test periodicity: derivatives at endpoints should match
        let eps = 1e-6;
        let d0_left = (spline.evaluate(eps)? - spline.evaluate(0.0)?) / eps;
        let d0_right = (spline.evaluate(1.0)? - spline.evaluate(1.0 - eps)?) / eps;

        // For a truly periodic spline, the derivatives should match
        // Note: Using finite differences can introduce some numerical error
        let deriv_diff = (d0_left - d0_right).abs();
        assert!(
            deriv_diff < 1.0,
            "Periodic spline: derivatives at endpoints should match (diff={}, left={}, right={})",
            deriv_diff,
            d0_left,
            d0_right
        );

        // Test second derivatives at endpoints (more stringent periodicity check)
        let d2_left = spline.derivative(0.0)?;
        let d2_right = spline.derivative(1.0)?;

        // Second derivatives should be close for periodic spline
        assert!(
            (d2_left - d2_right).abs() < 0.5,
            "Periodic spline: second derivatives should match: left={}, right={}",
            d2_left,
            d2_right
        );
        Ok(())
    }

    #[test]
    fn test_periodic_spline_non_periodic_data() {
        // Test that periodic spline errors for non-periodic data
        let x = Array::from_vec(vec![0.0, 1.0, 2.0, 3.0]);
        let y = Array::from_vec(vec![0.0, 1.0, 4.0, 10.0]); // y[0] != y[3]

        let result = CubicSplineInterp::new(&x, &y, SplineBoundary::Periodic);
        assert!(
            result.is_err(),
            "Periodic spline should error for non-periodic data"
        );
    }

    #[test]
    fn test_boundary_conditions_comparison() -> Result<()> {
        // Compare different boundary conditions on the same data
        let x = Array::from_vec(vec![0.0, 1.0, 2.0, 3.0, 4.0]);
        let y = Array::from_vec(vec![1.0, 2.0, 1.5, 3.0, 2.5]);

        let spline_natural = CubicSplineInterp::natural(&x, &y)?;
        let spline_not_a_knot = CubicSplineInterp::not_a_knot(&x, &y)?;
        let spline_clamped = CubicSplineInterp::clamped(&x, &y, 0.5, -0.5)?;

        // All should pass through data points
        for i in 0..5 {
            let xi = i as f64;
            let yi_expected = y.get(&[i])?;

            let yi_natural = spline_natural.evaluate(xi)?;
            let yi_not_a_knot = spline_not_a_knot.evaluate(xi)?;
            let yi_clamped = spline_clamped.evaluate(xi)?;

            assert!((yi_natural - yi_expected).abs() < 1e-10);
            assert!((yi_not_a_knot - yi_expected).abs() < 1e-10);
            assert!((yi_clamped - yi_expected).abs() < 1e-10);
        }

        // They should differ at interpolation points (different boundary conditions)
        let x_test = 0.5;
        let val_natural = spline_natural.evaluate(x_test)?;
        let val_not_a_knot = spline_not_a_knot.evaluate(x_test)?;

        // Not-a-knot and natural usually give different results
        // (though for some data they might be close)
        println!("Natural: {}, Not-a-knot: {}", val_natural, val_not_a_knot);
        Ok(())
    }
}