scirs2_interpolate/

statistical.rs

//! Advanced statistical interpolation methods
//!
//! This module provides statistical interpolation techniques that go beyond
//! deterministic interpolation, including:
//! - Bootstrap confidence intervals
//! - Bayesian interpolation with posterior distributions
//! - Quantile interpolation/regression
//! - Robust interpolation methods
//! - Stochastic interpolation for random fields

#![allow(clippy::too_many_arguments)]
#![allow(dead_code)]

use crate::error::{InterpolateError, InterpolateResult};
use scirs2_core::ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis};
use scirs2_core::numeric::{Float, FromPrimitive};
use scirs2_core::random::{rngs::StdRng, Rng, SeedableRng};
use scirs2_core::random::{Distribution, Normal, StandardNormal};
use statrs::statistics::Statistics;
use std::fmt::{Debug, Display};

/// Configuration for bootstrap confidence intervals
#[derive(Debug, Clone)]
pub struct BootstrapConfig {
    /// Number of bootstrap samples
    pub n_samples: usize,
    /// Confidence level (e.g., 0.95 for 95% CI)
    pub confidence_level: f64,
    /// Random seed for reproducibility
    pub seed: Option<u64>,
}

impl Default for BootstrapConfig {
    fn default() -> Self {
        Self {
            n_samples: 1000,
            confidence_level: 0.95,
            seed: None,
        }
    }
}

/// Result from bootstrap interpolation including confidence intervals
#[derive(Debug, Clone)]
pub struct BootstrapResult<T: Float> {
    /// Point estimate (median of bootstrap samples)
    pub estimate: Array1<T>,
    /// Lower confidence bound
    pub lower_bound: Array1<T>,
    /// Upper confidence bound
    pub upper_bound: Array1<T>,
    /// Standard error estimate
    pub std_error: Array1<T>,
}

/// Bootstrap interpolation with confidence intervals
///
/// This method performs interpolation with uncertainty quantification
/// using bootstrap resampling of the input data.
pub struct BootstrapInterpolator<T: Float> {
    /// Configuration for bootstrap
    config: BootstrapConfig,
    /// Base interpolator factory
    interpolator_factory:
        Box<dyn Fn(&ArrayView1<T>, &ArrayView1<T>) -> InterpolateResult<Box<dyn Fn(T) -> T>>>,
}

impl<T: Float + FromPrimitive + Debug + Display + std::iter::Sum> BootstrapInterpolator<T> {
    /// Create a new bootstrap interpolator
    pub fn new<F>(config: BootstrapConfig, interpolator_factory: F) -> Self
    where
        F: Fn(&ArrayView1<T>, &ArrayView1<T>) -> InterpolateResult<Box<dyn Fn(T) -> T>> + 'static,
    {
        Self {
            config,
            interpolator_factory: Box::new(interpolator_factory),
        }
    }

    /// Perform bootstrap interpolation at given points
    pub fn interpolate(
        &self,
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        xnew: &ArrayView1<T>,
    ) -> InterpolateResult<BootstrapResult<T>> {
        if x.len() != y.len() {
            return Err(InterpolateError::DimensionMismatch(
                "x and y must have the same length".to_string(),
            ));
        }

        let n = x.len();
        let m = xnew.len();
        let mut rng = match self.config.seed {
            Some(seed) => StdRng::seed_from_u64(seed),
            None => {
                let mut rng = scirs2_core::random::rng();
                StdRng::from_rng(&mut rng)
            }
        };

        // Storage for bootstrap samples
        let mut bootstrap_results = Array2::<T>::zeros((self.config.n_samples, m));

        // Perform bootstrap resampling
        for i in 0..self.config.n_samples {
            // Resample indices with replacement
            let indices: Vec<usize> = (0..n).map(|_| rng.gen_range(0..n)).collect();

            // Create resampled data
            let x_resampled = indices.iter().map(|&idx| x[idx]).collect::<Array1<_>>();
            let y_resampled = indices.iter().map(|&idx| y[idx]).collect::<Array1<_>>();

            // Create interpolator for this bootstrap sample
            let interpolator =
                (self.interpolator_factory)(&x_resampled.view(), &y_resampled.view())?;

            // Evaluate at _new points
            for (j, &x_val) in xnew.iter().enumerate() {
                bootstrap_results[[i, j]] = interpolator(x_val);
            }
        }

        // Calculate statistics
        let alpha = T::from(1.0 - self.config.confidence_level).unwrap();
        let lower_percentile = alpha / T::from(2.0).unwrap();
        let upper_percentile = T::one() - alpha / T::from(2.0).unwrap();

        let mut estimate = Array1::zeros(m);
        let mut lower_bound = Array1::zeros(m);
        let mut upper_bound = Array1::zeros(m);
        let mut std_error = Array1::zeros(m);

        for j in 0..m {
            let column = bootstrap_results.index_axis(Axis(1), j);
            let mut sorted_col = column.to_vec();
            sorted_col.sort_by(|a, b| a.partial_cmp(b).unwrap());

            // Median as point estimate
            let median_idx = self.config.n_samples / 2;
            estimate[j] = sorted_col[median_idx];

            // Confidence bounds
            let lower_idx = (lower_percentile * T::from(self.config.n_samples).unwrap())
                .to_usize()
                .unwrap();
            let upper_idx = (upper_percentile * T::from(self.config.n_samples).unwrap())
                .to_usize()
                .unwrap();
            lower_bound[j] = sorted_col[lower_idx];
            upper_bound[j] = sorted_col[upper_idx];

            // Standard error
            let mean = column.mean().unwrap();
            let variance = column
                .iter()
                .map(|&val| (val - mean) * (val - mean))
                .sum::<T>()
                / T::from(self.config.n_samples - 1).unwrap();
            std_error[j] = variance.sqrt();
        }

        Ok(BootstrapResult {
            estimate,
            lower_bound,
            upper_bound,
            std_error,
        })
    }
}

/// Configuration for Bayesian interpolation
pub struct BayesianConfig<T: Float> {
    /// Prior mean function
    pub prior_mean: Box<dyn Fn(T) -> T>,
    /// Prior variance
    pub prior_variance: T,
    /// Measurement noise variance
    pub noise_variance: T,
    /// RBF kernel length scale parameter
    pub length_scale: T,
    /// Number of posterior samples to draw
    pub n_posterior_samples: usize,
}

impl<T: Float + FromPrimitive> Default for BayesianConfig<T> {
    fn default() -> Self {
        Self {
            prior_mean: Box::new(|_| T::zero()),
            prior_variance: T::one(),
            noise_variance: T::from(0.01).unwrap(),
            length_scale: T::one(),
            n_posterior_samples: 100,
        }
    }
}

impl<T: Float + FromPrimitive> BayesianConfig<T> {
    /// Set the RBF kernel length scale parameter
    pub fn with_length_scale(mut self, lengthscale: T) -> Self {
        self.length_scale = lengthscale;
        self
    }

    /// Set the prior variance
    pub fn with_prior_variance(mut self, variance: T) -> Self {
        self.prior_variance = variance;
        self
    }

    /// Set the noise variance
    pub fn with_noise_variance(mut self, variance: T) -> Self {
        self.noise_variance = variance;
        self
    }

    /// Set the number of posterior samples
    pub fn with_n_posterior_samples(mut self, nsamples: usize) -> Self {
        self.n_posterior_samples = nsamples;
        self
    }
}

/// Bayesian interpolation with full posterior distribution
///
/// This provides interpolation with full uncertainty quantification
/// through Bayesian inference.
pub struct BayesianInterpolator<T: Float> {
    config: BayesianConfig<T>,
    x_obs: Array1<T>,
    y_obs: Array1<T>,
}

impl<
        T: Float
            + FromPrimitive
            + Debug
            + Display
            + std::ops::AddAssign
            + std::ops::SubAssign
            + std::ops::MulAssign
            + std::ops::DivAssign
            + std::ops::RemAssign,
    > BayesianInterpolator<T>
{
    /// Create a new Bayesian interpolator
    pub fn new(
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        config: BayesianConfig<T>,
    ) -> InterpolateResult<Self> {
        if x.len() != y.len() {
            return Err(InterpolateError::DimensionMismatch(
                "x and y must have the same length".to_string(),
            ));
        }

        Ok(Self {
            config,
            x_obs: x.to_owned(),
            y_obs: y.to_owned(),
        })
    }

    /// Get posterior mean at given points using proper Gaussian process regression
    pub fn posterior_mean(&self, xnew: &ArrayView1<T>) -> InterpolateResult<Array1<T>> {
        let n = self.x_obs.len();
        let m = xnew.len();

        if n == 0 {
            return Err(InterpolateError::invalid_input(
                "No observed data points".to_string(),
            ));
        }

        // Compute covariance matrix K(X, X) + σ²I
        let mut k_xx = Array2::<T>::zeros((n, n));
        let length_scale = self.config.length_scale;

        // Build covariance matrix with RBF kernel
        for i in 0..n {
            for j in 0..n {
                let dist_sq = (self.x_obs[i] - self.x_obs[j]).powi(2);
                k_xx[[i, j]] = self.config.prior_variance
                    * (-dist_sq / (T::from(2.0).unwrap() * length_scale.powi(2))).exp();

                // Add noise variance to diagonal
                if i == j {
                    k_xx[[i, j]] += self.config.noise_variance;
                }
            }
        }

        // Solve the linear system K * weights = y_obs using Cholesky decomposition
        // This is more numerically stable than matrix inversion
        let weights = match self.solve_gp_system(&k_xx.view(), &self.y_obs.view()) {
            Ok(w) => w,
            Err(_) => {
                // Fallback to regularized system if Cholesky fails
                let regularization = T::from(1e-6).unwrap();
                for i in 0..n {
                    k_xx[[i, i]] += regularization;
                }
                self.solve_gp_system(&k_xx.view(), &self.y_obs.view())?
            }
        };

        // Compute cross-covariance K(X*, X)
        let mut k_star_x = Array2::<T>::zeros((m, n));
        for i in 0..m {
            for j in 0..n {
                let dist_sq = (xnew[i] - self.x_obs[j]).powi(2);
                k_star_x[[i, j]] = self.config.prior_variance
                    * (-dist_sq / (T::from(2.0).unwrap() * length_scale.powi(2))).exp();
            }
        }

        // Compute posterior mean: μ* = K(X*, X) * weights
        let mut mean = Array1::zeros(m);
        for i in 0..m {
            let mut sum = T::zero();
            for j in 0..n {
                sum += k_star_x[[i, j]] * weights[j];
            }
            // Add prior mean
            mean[i] = (self.config.prior_mean)(xnew[i]) + sum;
        }

        Ok(mean)
    }

    /// Solve the GP linear system using available numerical methods
    fn solve_gp_system(
        &self,
        k_matrix: &ArrayView2<T>,
        y_obs: &ArrayView1<T>,
    ) -> InterpolateResult<Array1<T>> {
        use crate::structured_matrix::solve_dense_system;

        // Try using the structured _matrix solver
        match solve_dense_system(k_matrix, y_obs) {
            Ok(solution) => Ok(solution),
            Err(_) => {
                // Additional fallback: use simple weighted average if _matrix is ill-conditioned
                let n = y_obs.len();
                let weights = Array1::from_elem(n, T::one() / T::from(n).unwrap());
                Ok(weights)
            }
        }
    }

    /// Draw samples from the posterior distribution
    pub fn posterior_samples(
        &self,
        xnew: &ArrayView1<T>,
        n_samples: usize,
    ) -> InterpolateResult<Array2<T>> {
        let mean = self.posterior_mean(xnew)?;
        let m = xnew.len();

        // For computational efficiency, we use a simplified approach that captures
        // the main posterior uncertainty while avoiding expensive matrix operations.
        // A full implementation would compute the posterior covariance matrix:
        // Σ* = K(X*, X*) - K(X*, X)[K(X, X) + σ²I]^(-1)K(X, X*)

        let mut samples = Array2::zeros((n_samples, m));
        let mut rng = scirs2_core::random::rng();

        // Compute approximate posterior variance at each point
        let length_scale = T::one();
        for j in 0..m {
            // Compute posterior variance as prior variance minus reduction from observations
            let mut reduction_factor = T::zero();
            let mut total_influence = T::zero();

            for i in 0..self.x_obs.len() {
                let dist_sq = (xnew[j] - self.x_obs[i]).powi(2);
                let influence = (-dist_sq / (T::from(2.0).unwrap() * length_scale.powi(2))).exp();
                total_influence += influence;
                reduction_factor += influence * influence;
            }

            // Approximate posterior variance
            let noise_ratio = self.config.noise_variance / self.config.prior_variance;
            let posterior_var = self.config.prior_variance
                * (T::one()
                    - reduction_factor / (total_influence + noise_ratio + T::from(1e-8).unwrap()));

            // Ensure positive variance
            let std_dev = posterior_var.max(T::from(1e-12).unwrap()).sqrt();

            // Draw _samples for this query point
            for i in 0..n_samples {
                if let Ok(normal) =
                    Normal::new(mean[j].to_f64().unwrap(), std_dev.to_f64().unwrap())
                {
                    samples[[i, j]] = T::from(normal.sample(&mut rng)).unwrap();
                } else {
                    samples[[i, j]] = mean[j];
                }
            }
        }

        Ok(samples)
    }
}

/// Quantile interpolation/regression
///
/// Interpolates specific quantiles of the response distribution
pub struct QuantileInterpolator<T: Float> {
    /// Quantile to interpolate (e.g., 0.5 for median)
    quantile: T,
    /// Bandwidth for local quantile estimation
    bandwidth: T,
}

impl<T: Float + FromPrimitive + Debug + Display> QuantileInterpolator<T>
where
    T: std::iter::Sum<T> + for<'a> std::iter::Sum<&'a T>,
{
    /// Create a new quantile interpolator
    pub fn new(quantile: T, bandwidth: T) -> InterpolateResult<Self> {
        if quantile <= T::zero() || quantile >= T::one() {
            return Err(InterpolateError::InvalidValue(
                "Quantile must be between 0 and 1".to_string(),
            ));
        }

        Ok(Self {
            quantile,
            bandwidth,
        })
    }

    /// Interpolate quantile at given points
    pub fn interpolate(
        &self,
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        xnew: &ArrayView1<T>,
    ) -> InterpolateResult<Array1<T>> {
        if x.len() != y.len() {
            return Err(InterpolateError::DimensionMismatch(
                "x and y must have the same length".to_string(),
            ));
        }

        let n = x.len();
        let m = xnew.len();
        let mut result = Array1::zeros(m);

        // Local quantile regression
        for j in 0..m {
            let x_target = xnew[j];

            // Compute weights based on distance
            let mut weights = Vec::with_capacity(n);
            let mut weighted_values = Vec::with_capacity(n);

            for i in 0..n {
                let dist = (x[i] - x_target).abs() / self.bandwidth;
                let weight = if dist < T::one() {
                    (T::one() - dist * dist * dist).powi(3) // Tricube kernel
                } else {
                    T::zero()
                };

                if weight > T::epsilon() {
                    weights.push(weight);
                    weighted_values.push((y[i], weight));
                }
            }

            if weighted_values.is_empty() {
                // No nearby points, use nearest neighbor
                let nearest_idx = x
                    .iter()
                    .enumerate()
                    .min_by_key(|(_, &xi)| ((xi - x_target).abs().to_f64().unwrap() * 1e6) as i64)
                    .map(|(i_, _)| i_)
                    .unwrap();
                result[j] = y[nearest_idx];
            } else {
                // Sort by value
                weighted_values.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());

                // Find weighted quantile
                let total_weight: T = weights.iter().sum();
                let target_weight = self.quantile * total_weight;

                let mut cumulative_weight = T::zero();
                for (val, weight) in weighted_values {
                    cumulative_weight = cumulative_weight + weight;
                    if cumulative_weight >= target_weight {
                        result[j] = val;
                        break;
                    }
                }
            }
        }

        Ok(result)
    }
}

/// Robust interpolation methods resistant to outliers
pub struct RobustInterpolator<T: Float> {
    /// Tuning constant for robustness
    tuning_constant: T,
    /// Maximum iterations for iterative reweighting
    max_iterations: usize,
    /// Convergence tolerance
    tolerance: T,
}

impl<T: Float + FromPrimitive + Debug + Display> RobustInterpolator<T> {
    /// Create a new robust interpolator using M-estimation
    pub fn new(_tuningconstant: T) -> Self {
        Self {
            tuning_constant: _tuningconstant,
            max_iterations: 100,
            tolerance: T::from(1e-6).unwrap(),
        }
    }

    /// Perform robust interpolation using iteratively reweighted least squares
    pub fn interpolate(
        &self,
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        xnew: &ArrayView1<T>,
    ) -> InterpolateResult<Array1<T>> {
        // Use local polynomial regression with robust weights
        let n = x.len();
        let m = xnew.len();
        let mut result = Array1::zeros(m);

        for j in 0..m {
            let x_target = xnew[j];

            // Initial weights (uniform)
            let mut weights = vec![T::one(); n];
            let mut prev_estimate = T::zero();

            // Iteratively reweighted least squares
            for _iter in 0..self.max_iterations {
                // Weighted linear regression
                let mut sum_w = T::zero();
                let mut sum_wx = T::zero();
                let mut sum_wy = T::zero();
                let mut sum_wxx = T::zero();
                let mut sum_wxy = T::zero();

                for i in 0..n {
                    let w = weights[i];
                    let dx = x[i] - x_target;
                    sum_w = sum_w + w;
                    sum_wx = sum_wx + w * dx;
                    sum_wy = sum_wy + w * y[i];
                    sum_wxx = sum_wxx + w * dx * dx;
                    sum_wxy = sum_wxy + w * dx * y[i];
                }

                // Solve for coefficients
                let det = sum_w * sum_wxx - sum_wx * sum_wx;
                let estimate = if det.abs() > T::epsilon() {
                    (sum_wxx * sum_wy - sum_wx * sum_wxy) / det
                } else {
                    sum_wy / sum_w
                };

                // Check convergence
                if (estimate - prev_estimate).abs() < self.tolerance {
                    result[j] = estimate;
                    break;
                }
                prev_estimate = estimate;

                // Update weights using Huber's psi function
                for i in 0..n {
                    let residual = y[i] - estimate;
                    let scaled_residual = residual / self.tuning_constant;

                    weights[i] = if scaled_residual.abs() <= T::one() {
                        T::one()
                    } else {
                        T::one() / scaled_residual.abs()
                    };
                }
            }

            result[j] = prev_estimate;
        }

        Ok(result)
    }
}

/// Stochastic interpolation for random fields
///
/// Provides interpolation that preserves the stochastic properties
/// of the underlying random field.
pub struct StochasticInterpolator<T: Float> {
    /// Correlation length scale
    correlation_length: T,
    /// Field variance
    field_variance: T,
    /// Number of realizations to generate
    n_realizations: usize,
}

impl<T: Float + FromPrimitive + Debug + Display> StochasticInterpolator<T> {
    /// Create a new stochastic interpolator
    pub fn new(correlation_length: T, field_variance: T, n_realizations: usize) -> Self {
        Self {
            correlation_length,
            field_variance,
            n_realizations,
        }
    }

    /// Generate stochastic realizations of the interpolated field
    pub fn interpolate_realizations(
        &self,
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        xnew: &ArrayView1<T>,
    ) -> InterpolateResult<Array2<T>> {
        let n = x.len();
        let m = xnew.len();
        let mut realizations = Array2::zeros((self.n_realizations, m));

        let mut rng = scirs2_core::random::rng();

        for r in 0..self.n_realizations {
            // Generate a realization using conditional simulation
            for j in 0..m {
                let x_target = xnew[j];

                // Kriging interpolation with added noise
                let mut weighted_sum = T::zero();
                let mut weight_sum = T::zero();

                for i in 0..n {
                    let dist = (x[i] - x_target).abs() / self.correlation_length;
                    let weight = (-dist * dist).exp();
                    weighted_sum = weighted_sum + weight * y[i];
                    weight_sum = weight_sum + weight;
                }

                let mean = if weight_sum > T::epsilon() {
                    weighted_sum / weight_sum
                } else {
                    T::zero()
                };

                // Add stochastic component
                let std_dev =
                    (self.field_variance * (T::one() - weight_sum / T::from(n).unwrap())).sqrt();
                let normal_sample: f64 = StandardNormal.sample(&mut rng);
                let noise: T = T::from(normal_sample).unwrap() * std_dev;

                realizations[[r, j]] = mean + noise;
            }
        }

        Ok(realizations)
    }

    /// Get mean and variance of the stochastic interpolation
    pub fn interpolate_statistics(
        &self,
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        xnew: &ArrayView1<T>,
    ) -> InterpolateResult<(Array1<T>, Array1<T>)> {
        let realizations = self.interpolate_realizations(x, y, xnew)?;

        let mean = realizations.mean_axis(Axis(0)).unwrap();
        let variance = realizations.var_axis(Axis(0), T::from(1.0).unwrap());

        Ok((mean, variance))
    }
}

/// Factory functions for creating statistical interpolators
/// Create a bootstrap interpolator with linear base interpolation
#[allow(dead_code)]
pub fn make_bootstrap_linear_interpolator<
    T: Float + FromPrimitive + Debug + Display + 'static + std::iter::Sum,
>(
    config: BootstrapConfig,
) -> BootstrapInterpolator<T> {
    BootstrapInterpolator::new(config, |x, y| {
        // Create a simple linear interpolator
        let x_owned = x.to_owned();
        let y_owned = y.to_owned();
        Ok(Box::new(move |xnew| {
            // Simple linear interpolation
            if xnew <= x_owned[0] {
                y_owned[0]
            } else if xnew >= x_owned[x_owned.len() - 1] {
                y_owned[y_owned.len() - 1]
            } else {
                // Find surrounding points
                let mut i = 0;
                for j in 1..x_owned.len() {
                    if xnew <= x_owned[j] {
                        i = j - 1;
                        break;
                    }
                }

                let alpha = (xnew - x_owned[i]) / (x_owned[i + 1] - x_owned[i]);
                y_owned[i] * (T::one() - alpha) + y_owned[i + 1] * alpha
            }
        }))
    })
}

/// Create a Bayesian interpolator with default configuration
#[allow(dead_code)]
pub fn make_bayesian_interpolator<T: crate::traits::InterpolationFloat>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
) -> InterpolateResult<BayesianInterpolator<T>> {
    BayesianInterpolator::new(x, y, BayesianConfig::default())
}

/// Create a median (0.5 quantile) interpolator
#[allow(dead_code)]
pub fn make_median_interpolator<T>(bandwidth: T) -> InterpolateResult<QuantileInterpolator<T>>
where
    T: Float + FromPrimitive + Debug + Display + std::iter::Sum<T> + for<'a> std::iter::Sum<&'a T>,
{
    QuantileInterpolator::new(T::from(0.5).unwrap(), bandwidth)
}

/// Create a robust interpolator with default Huber tuning
#[allow(dead_code)]
pub fn make_robust_interpolator<T: crate::traits::InterpolationFloat>() -> RobustInterpolator<T> {
    RobustInterpolator::new(T::from(1.345).unwrap()) // Huber's recommended value
}

/// Create a stochastic interpolator with default parameters
#[allow(dead_code)]
pub fn make_stochastic_interpolator<T: crate::traits::InterpolationFloat>(
    correlation_length: T,
) -> StochasticInterpolator<T> {
    StochasticInterpolator::new(correlation_length, T::one(), 100)
}

/// Ensemble interpolation combining multiple methods
///
/// Provides interpolation using an ensemble of different methods
/// to improve robustness and uncertainty quantification.
pub struct EnsembleInterpolator<T: Float> {
    /// Weight for each interpolation method
    weights: Array1<T>,
    /// Interpolation methods in the ensemble
    methods: Vec<
        Box<dyn Fn(&ArrayView1<T>, &ArrayView1<T>, &ArrayView1<T>) -> InterpolateResult<Array1<T>>>,
    >,
    /// Whether to normalize weights
    normalize_weights: bool,
}

impl<T: crate::traits::InterpolationFloat> EnsembleInterpolator<T> {
    /// Create a new ensemble interpolator
    pub fn new() -> Self {
        Self {
            weights: Array1::zeros(0),
            methods: Vec::new(),
            normalize_weights: true,
        }
    }

    /// Add a linear interpolation method to the ensemble
    pub fn add_linear_method(mut self, weight: T) -> Self {
        self.weights = if self.weights.is_empty() {
            Array1::from_vec(vec![weight])
        } else {
            let mut new_weights = self.weights.to_vec();
            new_weights.push(weight);
            Array1::from_vec(new_weights)
        };

        self.methods.push(Box::new(|x, y, xnew| {
            let mut result = Array1::zeros(xnew.len());
            for (i, &x_val) in xnew.iter().enumerate() {
                // Linear interpolation
                if x_val <= x[0] {
                    result[i] = y[0];
                } else if x_val >= x[x.len() - 1] {
                    result[i] = y[y.len() - 1];
                } else {
                    // Find surrounding points
                    for j in 1..x.len() {
                        if x_val <= x[j] {
                            let alpha = (x_val - x[j - 1]) / (x[j] - x[j - 1]);
                            result[i] = y[j - 1] * (T::one() - alpha) + y[j] * alpha;
                            break;
                        }
                    }
                }
            }
            Ok(result)
        }));
        self
    }

    /// Add a cubic interpolation method to the ensemble
    pub fn add_cubic_method(mut self, weight: T) -> Self {
        self.weights = if self.weights.is_empty() {
            Array1::from_vec(vec![weight])
        } else {
            let mut new_weights = self.weights.to_vec();
            new_weights.push(weight);
            Array1::from_vec(new_weights)
        };

        self.methods.push(Box::new(|x, y, xnew| {
            // Cubic spline interpolation using natural boundary conditions
            use crate::spline::CubicSpline;

            // Need at least 3 points for cubic spline
            if x.len() < 3 {
                return Err(InterpolateError::invalid_input(
                    "Cubic spline requires at least 3 data points".to_string(),
                ));
            }

            // Create cubic spline with natural boundary conditions
            let spline = CubicSpline::new(x, y)?;

            // Evaluate at all query points
            let mut result = Array1::zeros(xnew.len());
            for (i, &x_val) in xnew.iter().enumerate() {
                // Handle extrapolation by clamping to boundary values
                if x_val < x[0] {
                    result[i] = y[0];
                } else if x_val > x[x.len() - 1] {
                    result[i] = y[y.len() - 1];
                } else {
                    // Evaluate cubic spline within the valid range
                    result[i] = spline.evaluate(x_val)?;
                }
            }
            Ok(result)
        }));
        self
    }

    /// Perform ensemble interpolation
    pub fn interpolate(
        &self,
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        xnew: &ArrayView1<T>,
    ) -> InterpolateResult<Array1<T>> {
        if self.methods.is_empty() {
            return Err(InterpolateError::InvalidState(
                "No interpolation methods in ensemble".to_string(),
            ));
        }

        let mut weighted_results = Array1::zeros(xnew.len());
        let mut total_weight = T::zero();

        for (i, method) in self.methods.iter().enumerate() {
            let result = method(x, y, xnew)?;
            let weight = self.weights[i];

            for j in 0..xnew.len() {
                weighted_results[j] += weight * result[j];
            }
            total_weight += weight;
        }

        // Normalize if requested
        if self.normalize_weights && total_weight > T::zero() {
            for val in weighted_results.iter_mut() {
                *val /= total_weight;
            }
        }

        Ok(weighted_results)
    }

    /// Get ensemble variance (measure of uncertainty)
    pub fn interpolate_with_variance(
        &self,
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        xnew: &ArrayView1<T>,
    ) -> InterpolateResult<(Array1<T>, Array1<T>)> {
        if self.methods.is_empty() {
            return Err(InterpolateError::InvalidState(
                "No interpolation methods in ensemble".to_string(),
            ));
        }

        let mut all_results = Vec::new();

        // Collect results from all methods
        for method in self.methods.iter() {
            let result = method(x, y, xnew)?;
            all_results.push(result);
        }

        // Compute weighted mean
        let mut weighted_mean = Array1::zeros(xnew.len());
        let mut total_weight = T::zero();

        for (i, result) in all_results.iter().enumerate() {
            let weight = self.weights[i];
            for j in 0..xnew.len() {
                weighted_mean[j] += weight * result[j];
            }
            total_weight += weight;
        }

        if total_weight > T::zero() {
            for val in weighted_mean.iter_mut() {
                *val /= total_weight;
            }
        }

        // Compute weighted variance
        let mut variance = Array1::zeros(xnew.len());
        if all_results.len() > 1 {
            for (i, result) in all_results.iter().enumerate() {
                let weight = self.weights[i];
                for j in 0..xnew.len() {
                    let diff = result[j] - weighted_mean[j];
                    variance[j] += weight * diff * diff;
                }
            }

            if total_weight > T::zero() {
                for val in variance.iter_mut() {
                    *val /= total_weight;
                }
            }
        }

        Ok((weighted_mean, variance))
    }
}

impl<T: crate::traits::InterpolationFloat> Default for EnsembleInterpolator<T> {
    fn default() -> Self {
        Self::new()
    }
}

/// Cross-validation based uncertainty estimation
///
/// Provides uncertainty estimates using leave-one-out cross-validation
/// or k-fold cross-validation.
pub struct CrossValidationUncertainty {
    /// Number of folds for k-fold CV (if 0, use leave-one-out)
    k_folds: usize,
    /// Random seed for fold assignment
    seed: Option<u64>,
}

impl CrossValidationUncertainty {
    /// Create a new cross-validation uncertainty estimator
    pub fn new(_kfolds: usize) -> Self {
        Self {
            k_folds: _kfolds,
            seed: None,
        }
    }

    /// Set random seed for reproducible fold assignment
    pub fn with_seed(mut self, seed: u64) -> Self {
        self.seed = Some(seed);
        self
    }

    /// Estimate uncertainty using cross-validation
    pub fn estimate_uncertainty<T, F>(
        &self,
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        xnew: &ArrayView1<T>,
        interpolator_factory: F,
    ) -> InterpolateResult<(Array1<T>, Array1<T>)>
    where
        T: Clone
            + Copy
            + scirs2_core::numeric::Float
            + scirs2_core::numeric::FromPrimitive
            + std::iter::Sum,
        F: Fn(&ArrayView1<T>, &ArrayView1<T>, &ArrayView1<T>) -> InterpolateResult<Array1<T>>,
    {
        let n = x.len();
        let _m = xnew.len();

        if self.k_folds == 0 || self.k_folds >= n {
            // Leave-one-out cross-validation
            self.leave_one_out_uncertainty(x, y, xnew, interpolator_factory)
        } else {
            // K-fold cross-validation
            self.k_fold_uncertainty(x, y, xnew, interpolator_factory)
        }
    }

    fn leave_one_out_uncertainty<T, F>(
        &self,
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        xnew: &ArrayView1<T>,
        interpolator_factory: F,
    ) -> InterpolateResult<(Array1<T>, Array1<T>)>
    where
        T: Clone
            + Copy
            + scirs2_core::numeric::Float
            + scirs2_core::numeric::FromPrimitive
            + std::iter::Sum,
        F: Fn(&ArrayView1<T>, &ArrayView1<T>, &ArrayView1<T>) -> InterpolateResult<Array1<T>>,
    {
        let n = x.len();
        let m = xnew.len();
        let mut predictions = Array2::zeros((n, m));

        // Leave-one-out cross-validation
        for i in 0..n {
            // Create training set without point i
            let mut x_train = Vec::new();
            let mut y_train = Vec::new();

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

            let x_train_array = Array1::from_vec(x_train);
            let y_train_array = Array1::from_vec(y_train);

            // Train on reduced dataset and predict
            let pred = interpolator_factory(&x_train_array.view(), &y_train_array.view(), xnew)?;
            for j in 0..m {
                predictions[[i, j]] = pred[j];
            }
        }

        // Compute mean and variance of predictions
        let mut mean = Array1::zeros(m);
        let mut variance = Array1::zeros(m);

        for j in 0..m {
            let col = predictions.column(j);
            let sum: T = col.iter().copied().sum();
            mean[j] = sum / T::from(n).unwrap();

            let var_sum: T = col
                .iter()
                .map(|&val| (val - mean[j]) * (val - mean[j]))
                .sum();
            variance[j] = var_sum / T::from(n - 1).unwrap();
        }

        Ok((mean, variance))
    }

    fn k_fold_uncertainty<T, F>(
        &self,
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        xnew: &ArrayView1<T>,
        interpolator_factory: F,
    ) -> InterpolateResult<(Array1<T>, Array1<T>)>
    where
        T: Clone
            + Copy
            + scirs2_core::numeric::Float
            + scirs2_core::numeric::FromPrimitive
            + std::iter::Sum,
        F: Fn(&ArrayView1<T>, &ArrayView1<T>, &ArrayView1<T>) -> InterpolateResult<Array1<T>>,
    {
        let n = x.len();
        let m = xnew.len();
        let fold_size = n / self.k_folds;
        let mut predictions = Vec::new();

        let mut rng = match self.seed {
            Some(seed) => StdRng::seed_from_u64(seed),
            None => {
                let mut rng = scirs2_core::random::rng();
                StdRng::from_rng(&mut rng)
            }
        };

        // Create shuffled indices
        let mut indices: Vec<usize> = (0..n).collect();
        use scirs2_core::random::seq::SliceRandom;
        indices.shuffle(&mut rng);

        // K-fold cross-validation
        for fold in 0..self.k_folds {
            let start_idx = fold * fold_size;
            let end_idx = if fold == self.k_folds - 1 {
                n
            } else {
                (fold + 1) * fold_size
            };

            // Create training set (excluding current fold)
            let mut x_train = Vec::new();
            let mut y_train = Vec::new();

            for &idx in &indices[..start_idx] {
                x_train.push(x[idx]);
                y_train.push(y[idx]);
            }
            for &idx in &indices[end_idx..] {
                x_train.push(x[idx]);
                y_train.push(y[idx]);
            }

            let x_train_array = Array1::from_vec(x_train);
            let y_train_array = Array1::from_vec(y_train);

            // Train and predict
            let pred = interpolator_factory(&x_train_array.view(), &y_train_array.view(), xnew)?;
            predictions.push(pred);
        }

        // Compute statistics across folds
        let mut mean = Array1::zeros(m);
        let mut variance = Array1::zeros(m);

        for j in 0..m {
            let values: Vec<T> = predictions.iter().map(|pred| pred[j]).collect();
            let sum: T = values.iter().copied().sum();
            mean[j] = sum / T::from(self.k_folds).unwrap();

            let var_sum: T = values
                .iter()
                .map(|&val| (val - mean[j]) * (val - mean[j]))
                .sum();
            variance[j] = var_sum / T::from(self.k_folds - 1).unwrap();
        }

        Ok((mean, variance))
    }
}

/// Create an ensemble interpolator with linear and cubic methods
#[allow(dead_code)]
pub fn make_ensemble_interpolator<
    T: Float
        + FromPrimitive
        + Debug
        + Display
        + Copy
        + std::iter::Sum
        + crate::traits::InterpolationFloat,
>() -> EnsembleInterpolator<T> {
    EnsembleInterpolator::new()
        .add_linear_method(T::from(0.6).unwrap())
        .add_cubic_method(T::from(0.4).unwrap())
}

/// Create a cross-validation uncertainty estimator with leave-one-out
#[allow(dead_code)]
pub fn make_loocv_uncertainty() -> CrossValidationUncertainty {
    CrossValidationUncertainty::new(0) // 0 means leave-one-out
}

/// Create a cross-validation uncertainty estimator with k-fold CV
#[allow(dead_code)]
pub fn make_kfold_uncertainty(k: usize) -> CrossValidationUncertainty {
    CrossValidationUncertainty::new(k)
}

/// Isotonic (monotonic) regression interpolator
///
/// Performs interpolation while maintaining monotonicity constraints.
/// This is useful for dose-response relationships and other applications
/// where the underlying relationship must be monotonic.
#[derive(Debug, Clone)]
pub struct IsotonicInterpolator<T: Float> {
    /// Fitted isotonic values at training points
    fitted_values: Array1<T>,
    /// Training x coordinates (sorted)
    x_data: Array1<T>,
    /// Whether interpolation should be increasing (true) or decreasing (false)
    increasing: bool,
}

impl<T: Float + FromPrimitive + Debug + Display + Copy + std::iter::Sum> IsotonicInterpolator<T> {
    /// Create a new isotonic interpolator
    pub fn new(x: &ArrayView1<T>, y: &ArrayView1<T>, increasing: bool) -> InterpolateResult<Self> {
        if x.len() != y.len() {
            return Err(InterpolateError::DimensionMismatch(
                "x and y must have the same length".to_string(),
            ));
        }

        if x.len() < 2 {
            return Err(InterpolateError::invalid_input(
                "Need at least 2 points for isotonic regression".to_string(),
            ));
        }

        // Sort by x values
        let mut indices: Vec<usize> = (0..x.len()).collect();
        indices.sort_by(|&i, &j| x[i].partial_cmp(&x[j]).unwrap());

        let x_sorted: Array1<T> = indices.iter().map(|&i| x[i]).collect();
        let y_sorted: Array1<T> = indices.iter().map(|&i| y[i]).collect();

        // Apply pool-adjacent-violators algorithm
        let fitted_values = Self::pool_adjacent_violators(&y_sorted.view(), increasing)?;

        Ok(Self {
            fitted_values,
            x_data: x_sorted,
            increasing,
        })
    }

    /// Pool-adjacent-violators algorithm for isotonic regression
    fn pool_adjacent_violators(
        y: &ArrayView1<T>,
        increasing: bool,
    ) -> InterpolateResult<Array1<T>> {
        let n = y.len();
        let mut fitted = y.to_owned();
        let mut weights = Array1::<T>::ones(n);

        loop {
            let mut changed = false;

            for i in 0..n - 1 {
                let violates = if increasing {
                    fitted[i] > fitted[i + 1]
                } else {
                    fitted[i] < fitted[i + 1]
                };

                if violates {
                    // Pool adjacent blocks
                    let total_weight = weights[i] + weights[i + 1];
                    let weighted_mean =
                        (fitted[i] * weights[i] + fitted[i + 1] * weights[i + 1]) / total_weight;

                    fitted[i] = weighted_mean;
                    fitted[i + 1] = weighted_mean;
                    weights[i] = total_weight;
                    weights[i + 1] = total_weight;

                    changed = true;
                }
            }

            if !changed {
                break;
            }
        }

        Ok(fitted)
    }

    /// Interpolate at new points
    pub fn interpolate(&self, xnew: &ArrayView1<T>) -> InterpolateResult<Array1<T>> {
        let mut result = Array1::zeros(xnew.len());

        for (i, &x) in xnew.iter().enumerate() {
            // Find position in sorted data
            let idx = match self
                .x_data
                .as_slice()
                .unwrap()
                .binary_search_by(|&probe| probe.partial_cmp(&x).unwrap())
            {
                Ok(exact_idx) => {
                    result[i] = self.fitted_values[exact_idx];
                    continue;
                }
                Err(insert_idx) => insert_idx,
            };

            // Linear interpolation between adjacent fitted values
            if idx == 0 {
                result[i] = self.fitted_values[0];
            } else if idx >= self.x_data.len() {
                result[i] = self.fitted_values[self.x_data.len() - 1];
            } else {
                let x0 = self.x_data[idx - 1];
                let x1 = self.x_data[idx];
                let y0 = self.fitted_values[idx - 1];
                let y1 = self.fitted_values[idx];

                let t = (x - x0) / (x1 - x0);
                result[i] = y0 + t * (y1 - y0);
            }
        }

        Ok(result)
    }
}

/// Kernel Density Estimation (KDE) based interpolator
///
/// Uses kernel density estimation to create smooth interpolations
/// based on probability density functions.
#[derive(Debug, Clone)]
pub struct KDEInterpolator<T: Float> {
    /// Training data points
    x_data: Array1<T>,
    y_data: Array1<T>,
    /// Kernel bandwidth
    bandwidth: T,
    /// Kernel type
    kernel_type: KDEKernel,
}

/// Kernel types for KDE interpolation
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum KDEKernel {
    /// Gaussian (normal) kernel
    Gaussian,
    /// Epanechnikov kernel (more efficient)
    Epanechnikov,
    /// Triangular kernel
    Triangular,
    /// Uniform (box) kernel
    Uniform,
}

impl<T: Float + FromPrimitive + Debug + Display + Copy> KDEInterpolator<T> {
    /// Create a new KDE interpolator
    pub fn new(
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        bandwidth: T,
        kernel_type: KDEKernel,
    ) -> InterpolateResult<Self> {
        if x.len() != y.len() {
            return Err(InterpolateError::DimensionMismatch(
                "x and y must have the same length".to_string(),
            ));
        }

        if bandwidth <= T::zero() {
            return Err(InterpolateError::invalid_input(
                "Bandwidth must be positive".to_string(),
            ));
        }

        Ok(Self {
            x_data: x.to_owned(),
            y_data: y.to_owned(),
            bandwidth,
            kernel_type,
        })
    }

    /// Kernel function evaluation
    fn kernel(&self, u: T) -> T {
        match self.kernel_type {
            KDEKernel::Gaussian => {
                let pi = T::from(std::f64::consts::PI).unwrap();
                let two = T::from(2.0).unwrap();
                let exp_arg = -u * u / two;
                exp_arg.exp() / (two * pi).sqrt()
            }
            KDEKernel::Epanechnikov => {
                if u.abs() <= T::one() {
                    let three_fourths = T::from(0.75).unwrap();
                    three_fourths * (T::one() - u * u)
                } else {
                    T::zero()
                }
            }
            KDEKernel::Triangular => {
                if u.abs() <= T::one() {
                    T::one() - u.abs()
                } else {
                    T::zero()
                }
            }
            KDEKernel::Uniform => {
                if u.abs() <= T::one() {
                    T::from(0.5).unwrap()
                } else {
                    T::zero()
                }
            }
        }
    }

    /// Interpolate at new points using KDE
    pub fn interpolate(&self, xnew: &ArrayView1<T>) -> InterpolateResult<Array1<T>> {
        let mut result = Array1::zeros(xnew.len());

        for (i, &x) in xnew.iter().enumerate() {
            let mut weighted_sum = T::zero();
            let mut weight_sum = T::zero();

            for j in 0..self.x_data.len() {
                let u = (x - self.x_data[j]) / self.bandwidth;
                let kernel_weight = self.kernel(u);

                weighted_sum = weighted_sum + kernel_weight * self.y_data[j];
                weight_sum = weight_sum + kernel_weight;
            }

            if weight_sum > T::zero() {
                result[i] = weighted_sum / weight_sum;
            } else {
                // Fallback to nearest neighbor
                let mut min_dist = T::infinity();
                let mut nearest_y = self.y_data[0];

                for j in 0..self.x_data.len() {
                    let dist = (x - self.x_data[j]).abs();
                    if dist < min_dist {
                        min_dist = dist;
                        nearest_y = self.y_data[j];
                    }
                }

                result[i] = nearest_y;
            }
        }

        Ok(result)
    }
}

/// Empirical Bayes interpolator
///
/// Uses empirical Bayes methods for shrinkage-based interpolation.
/// Particularly useful when dealing with multiple related functions
/// or when prior information is available.
#[derive(Debug, Clone)]
pub struct EmpiricalBayesInterpolator<T: Float> {
    /// Training data
    x_data: Array1<T>,
    y_data: Array1<T>,
    /// Shrinkage parameters
    shrinkage_factor: T,
    /// Prior mean function
    prior_mean: T,
    /// Noise variance estimate
    noise_variance: T,
}

impl<T: Float + FromPrimitive + Debug + Display + Copy + std::iter::Sum>
    EmpiricalBayesInterpolator<T>
{
    /// Create a new empirical Bayes interpolator
    pub fn new(x: &ArrayView1<T>, y: &ArrayView1<T>) -> InterpolateResult<Self> {
        if x.len() != y.len() {
            return Err(InterpolateError::DimensionMismatch(
                "x and y must have the same length".to_string(),
            ));
        }

        if x.len() < 3 {
            return Err(InterpolateError::invalid_input(
                "Need at least 3 points for empirical Bayes".to_string(),
            ));
        }

        // Estimate prior mean (overall mean)
        let prior_mean = y.iter().copied().sum::<T>() / T::from(y.len()).unwrap();

        // Estimate noise variance using residuals
        let residuals: Array1<T> = y.iter().map(|&yi| yi - prior_mean).collect();
        let noise_variance =
            residuals.iter().map(|&r| r * r).sum::<T>() / T::from(residuals.len() - 1).unwrap();

        // Compute shrinkage factor using James-Stein type estimator
        let signal_variance = noise_variance.max(T::from(1e-10).unwrap());
        let shrinkage_factor = noise_variance / (noise_variance + signal_variance);

        Ok(Self {
            x_data: x.to_owned(),
            y_data: y.to_owned(),
            shrinkage_factor,
            prior_mean,
            noise_variance,
        })
    }

    /// Create empirical Bayes interpolator with custom prior
    pub fn with_prior(
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        prior_mean: T,
        shrinkage_factor: T,
    ) -> InterpolateResult<Self> {
        if x.len() != y.len() {
            return Err(InterpolateError::DimensionMismatch(
                "x and y must have the same length".to_string(),
            ));
        }

        let residuals: Array1<T> = y.iter().map(|&yi| yi - prior_mean).collect();
        let noise_variance =
            residuals.iter().map(|&r| r * r).sum::<T>() / T::from(residuals.len().max(1)).unwrap();

        Ok(Self {
            x_data: x.to_owned(),
            y_data: y.to_owned(),
            shrinkage_factor,
            prior_mean,
            noise_variance,
        })
    }

    /// Interpolate using empirical Bayes shrinkage
    pub fn interpolate(&self, xnew: &ArrayView1<T>) -> InterpolateResult<Array1<T>> {
        let mut result = Array1::zeros(xnew.len());

        for (i, &x) in xnew.iter().enumerate() {
            // Find nearest neighbors for local estimation
            let mut distances: Vec<(T, usize)> = self
                .x_data
                .iter()
                .enumerate()
                .map(|(j, &xi)| ((x - xi).abs(), j))
                .collect();
            distances.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());

            // Use k nearest neighbors (k = 3 or n/2, whichever is smaller)
            let k = (3_usize).min(self.x_data.len() / 2).max(1);
            let mut local_mean = T::zero();
            let mut total_weight = T::zero();

            for &(dist, j) in distances.iter().take(k) {
                let weight = if dist == T::zero() {
                    T::one()
                } else {
                    T::one() / (T::one() + dist)
                };
                local_mean = local_mean + weight * self.y_data[j];
                total_weight = total_weight + weight;
            }

            if total_weight > T::zero() {
                local_mean = local_mean / total_weight;
            } else {
                local_mean = self.prior_mean;
            }

            // Apply empirical Bayes shrinkage
            let shrunk_estimate = (T::one() - self.shrinkage_factor) * local_mean
                + self.shrinkage_factor * self.prior_mean;

            result[i] = shrunk_estimate;
        }

        Ok(result)
    }

    /// Get shrinkage factor
    pub fn get_shrinkage_factor(&self) -> T {
        self.shrinkage_factor
    }

    /// Get prior mean
    pub fn get_prior_mean(&self) -> T {
        self.prior_mean
    }

    /// Get noise variance estimate
    pub fn get_noise_variance(&self) -> T {
        self.noise_variance
    }
}

/// Convenience function to create an isotonic interpolator (increasing)
#[allow(dead_code)]
pub fn make_isotonic_interpolator<
    T: Float + FromPrimitive + Debug + Display + Copy + std::iter::Sum,
>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
) -> InterpolateResult<IsotonicInterpolator<T>> {
    IsotonicInterpolator::new(x, y, true)
}

/// Convenience function to create a decreasing isotonic interpolator
#[allow(dead_code)]
pub fn make_decreasing_isotonic_interpolator<
    T: Float + FromPrimitive + Debug + Display + Copy + std::iter::Sum,
>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
) -> InterpolateResult<IsotonicInterpolator<T>> {
    IsotonicInterpolator::new(x, y, false)
}

/// Convenience function to create a KDE interpolator with Gaussian kernel
#[allow(dead_code)]
pub fn make_kde_interpolator<T: crate::traits::InterpolationFloat + Copy>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
    bandwidth: T,
) -> InterpolateResult<KDEInterpolator<T>> {
    KDEInterpolator::new(x, y, bandwidth, KDEKernel::Gaussian)
}

/// Convenience function to create a KDE interpolator with automatic bandwidth selection
#[allow(dead_code)]
pub fn make_auto_kde_interpolator<
    T: Float + FromPrimitive + Debug + Display + Copy + std::iter::Sum,
>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
) -> InterpolateResult<KDEInterpolator<T>> {
    // Scott's rule for bandwidth selection
    let n = T::from(x.len()).unwrap();
    let x_std = {
        let mean = x.iter().copied().sum::<T>() / n;
        let variance = x.iter().map(|&xi| (xi - mean) * (xi - mean)).sum::<T>() / (n - T::one());
        variance.sqrt()
    };

    let bandwidth = x_std * n.powf(-T::from(0.2).unwrap()); // n^(-1/5)
    KDEInterpolator::new(x, y, bandwidth, KDEKernel::Gaussian)
}

/// Convenience function to create an empirical Bayes interpolator
#[allow(dead_code)]
pub fn make_empirical_bayes_interpolator<
    T: Float + FromPrimitive + Debug + Display + Copy + std::iter::Sum,
>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
) -> InterpolateResult<EmpiricalBayesInterpolator<T>> {
    EmpiricalBayesInterpolator::new(x, y)
}