scirs2-interpolate 0.4.2

//! Local Polynomial Regression Models
//!
//! This module provides implementations of local polynomial regression methods,
//! which fit polynomial models locally around each prediction point using
//! weighted least squares. These methods are particularly useful for nonparametric
//! regression and smoothing noisy data.
//!
//! The main difference from Moving Least Squares is the additional focus on
//! regression diagnostics, bandwidth selection, and statistical properties.

use scirs2_core::ndarray::ArrayStatCompat;
use scirs2_core::ndarray::{Array1, Array2, ArrayView1, ArrayView2};
use scirs2_core::numeric::{Float, FromPrimitive};
use std::fmt::Debug;
use std::marker::PhantomData;

use super::mls::{PolynomialBasis, WeightFunction};
use crate::error::{InterpolateError, InterpolateResult};
use statrs::statistics::Statistics;

/// Local polynomial regression model result
#[derive(Debug, Clone)]
pub struct RegressionResult<F: Float> {
    /// Fitted value at the query point
    pub value: F,

    /// Standard error of the fitted value
    pub std_error: F,

    /// Confidence interval (lower, upper) if confidence level was specified
    pub confidence_interval: Option<(F, F)>,

    /// Local polynomial coefficients
    pub coefficients: Array1<F>,

    /// Effective degrees of freedom used in this fit
    pub effective_df: F,

    /// Local R² (coefficient of determination)
    pub r_squared: F,
}

/// Configuration for local polynomial regression
#[derive(Debug, Clone)]
pub struct LocalPolynomialConfig<F: Float> {
    /// Bandwidth parameter (controls locality/smoothness)
    pub bandwidth: F,

    /// Weight function to use
    pub weight_fn: WeightFunction,

    /// Polynomial basis to use
    pub basis: PolynomialBasis,

    /// Confidence level for intervals (e.g., 0.95)
    pub confidence_level: Option<F>,

    /// Whether to compute robust standard errors
    pub robust_se: bool,

    /// Maximum number of points to use (for efficiency)
    pub max_points: Option<usize>,

    /// Small value to add to denominators for numerical stability
    pub epsilon: F,
}

impl<F: Float + FromPrimitive> Default for LocalPolynomialConfig<F> {
    fn default() -> Self {
        Self {
            bandwidth: F::from_f64(0.2).expect("Operation failed"),
            weight_fn: WeightFunction::Gaussian,
            basis: PolynomialBasis::Linear,
            confidence_level: None,
            robust_se: false,
            max_points: None,
            epsilon: F::from_f64(1e-10).expect("Operation failed"),
        }
    }
}

/// Local Polynomial Regression model
///
/// This model fits a polynomial of specified degree locally around
/// each prediction point using weighted least squares. The weights
/// depend on the distance between the prediction point and the data points.
///
/// The implementation includes:
/// - Multiple weight function options
/// - Polynomial basis options (constant, linear, quadratic)
/// - Standard error estimates
/// - Optional confidence intervals
/// - Local R² and effective degrees of freedom
/// - Support for robust standard errors
///
/// # Examples
///
/// ```
/// # #[cfg(feature = "linalg")]
/// # {
/// use scirs2_core::ndarray::{Array1, Array2, Axis};
/// use scirs2_interpolate::local::polynomial::{
///     LocalPolynomialRegression, LocalPolynomialConfig
/// };
/// use scirs2_interpolate::local::mls::{WeightFunction, PolynomialBasis};
///
/// // Create some 1D data with noise
/// let x = Array1::<f64>::linspace(0.0, 10.0, 50);
/// let mut y = Array1::<f64>::zeros(50);
/// for (i, x_val) in x.iter().enumerate() {
///     // y = sin(x) + some noise
///     y[i] = x_val.sin() + 0.1 * 0.3;
/// }
///
/// // Create a 2D array of points
/// let points = x.clone().insert_axis(Axis(1));
///
/// // Configure and create the local polynomial regression model
/// let config = LocalPolynomialConfig::<f64> {
///     bandwidth: 1.0,
///     weight_fn: WeightFunction::Gaussian,
///     basis: PolynomialBasis::Quadratic,
///     confidence_level: Some(0.95),
///     ..LocalPolynomialConfig::default()
/// };
///
/// let loess = LocalPolynomialRegression::with_config(
///     points,
///     y,
///     config
/// ).expect("Operation failed");
///
/// // Predict at a new point
/// let query = Array1::from_vec(vec![5.0]);
/// let result = loess.fit_at_point(&query.view()).expect("Operation failed");
///
/// // Access the fitted value
/// println!("Fitted value: {}", result.value);
///
/// // Access confidence interval
/// if let Some((lower, upper)) = result.confidence_interval {
///     println!("95% CI: ({}, {})", lower, upper);
/// }
/// # }
/// ```
#[derive(Debug, Clone)]
pub struct LocalPolynomialRegression<F>
where
    F: Float + FromPrimitive + Debug,
{
    /// Points coordinates (input locations)
    points: Array2<F>,

    /// Values at points
    values: Array1<F>,

    /// Configuration for the regression
    config: LocalPolynomialConfig<F>,

    /// Precomputed standard deviation of the response
    response_sd: F,

    /// Marker for generic type parameter
    _phantom: PhantomData<F>,
}

impl<F> LocalPolynomialRegression<F>
where
    F: Float + FromPrimitive + Debug,
{
    /// Get the precomputed standard deviation of the response
    pub fn response_sd(&self) -> F {
        self.response_sd
    }
}

impl<F> LocalPolynomialRegression<F>
where
    F: Float + FromPrimitive + Debug + 'static,
{
    /// Create a new LocalPolynomialRegression with default configuration
    ///
    /// # Arguments
    ///
    /// * `points` - Point coordinates with shape (n_points, n_dims)
    /// * `values` - Values at each point with shape (n_points,)
    /// * `bandwidth` - Bandwidth parameter controlling locality (larger = smoother)
    ///
    /// # Returns
    ///
    /// A new LocalPolynomialRegression model
    pub fn new(points: Array2<F>, values: Array1<F>, bandwidth: F) -> InterpolateResult<Self> {
        let config = LocalPolynomialConfig {
            bandwidth,
            ..LocalPolynomialConfig::default()
        };

        Self::with_config(points, values, config)
    }

    /// Create a new LocalPolynomialRegression with custom configuration
    ///
    /// # Arguments
    ///
    /// * `points` - Point coordinates with shape (n_points, n_dims)
    /// * `values` - Values at each point with shape (n_points,)
    /// * `config` - Configuration for the regression
    ///
    /// # Returns
    ///
    /// A new LocalPolynomialRegression model
    pub fn with_config(
        points: Array2<F>,
        values: Array1<F>,
        config: LocalPolynomialConfig<F>,
    ) -> InterpolateResult<Self> {
        // Validate inputs
        if points.shape()[0] != values.len() {
            return Err(InterpolateError::DimensionMismatch(
                "Number of points must match number of values".to_string(),
            ));
        }

        if points.shape()[0] < 2 {
            return Err(InterpolateError::InvalidValue(
                "At least 2 points are required for local polynomial regression".to_string(),
            ));
        }

        if config.bandwidth <= F::zero() {
            return Err(InterpolateError::InvalidValue(
                "Bandwidth parameter must be positive".to_string(),
            ));
        }

        // Precompute standard deviation of the response for standardization
        let mean = values.sum() / F::from_usize(values.len()).expect("Operation failed");
        let sum_squared_dev = values.fold(F::zero(), |acc, &v| acc + (v - mean).powi(2));
        let variance = sum_squared_dev / F::from_usize(values.len() - 1).expect("Operation failed");
        let response_sd = variance.sqrt();

        Ok(Self {
            points,
            values,
            config,
            response_sd,
            _phantom: PhantomData,
        })
    }

    /// Fit the local polynomial regression at a given point
    ///
    /// # Arguments
    ///
    /// * `x` - Query point with shape (n_dims,)
    ///
    /// # Returns
    ///
    /// Regression result at the query point
    pub fn fit_at_point(&self, x: &ArrayView1<F>) -> InterpolateResult<RegressionResult<F>> {
        // Check dimensions
        if x.len() != self.points.shape()[1] {
            return Err(InterpolateError::DimensionMismatch(
                "Query point dimension must match training points".to_string(),
            ));
        }

        // Get points to use for local fit
        let (indices, distances) = self.find_relevant_points(x)?;

        if indices.is_empty() {
            return Err(InterpolateError::invalid_input(
                "No points found within effective range".to_string(),
            ));
        }

        // Compute weights
        let weights = self.compute_weights(&distances)?;

        // Create basis functions for these points
        let local_points = self.extract_local_points(&indices);
        let basis_functions = self.create_basis_functions(&local_points, x)?;

        // Get local values
        let local_values = self.extract_local_values(&indices);

        // Perform weighted least squares fit
        let result = self.fit_weighted_least_squares(
            &local_points,
            &local_values,
            x,
            &weights,
            &basis_functions,
        )?;

        Ok(result)
    }

    /// Fit the local polynomial regression at multiple points
    ///
    /// # Arguments
    ///
    /// * `points` - Query points with shape (n_points, n_dims)
    ///
    /// # Returns
    ///
    /// Array of fitted values at the query points
    pub fn fit_multiple(&self, points: &ArrayView2<F>) -> InterpolateResult<Array1<F>> {
        // Check dimensions
        if points.shape()[1] != self.points.shape()[1] {
            return Err(InterpolateError::DimensionMismatch(
                "Query points dimension must match training points".to_string(),
            ));
        }

        let n_points = points.shape()[0];
        let mut results = Array1::zeros(n_points);

        // Fit at each point
        for i in 0..n_points {
            let point = points.slice(scirs2_core::ndarray::s![i, ..]);
            let result = self.fit_at_point(&point)?;
            results[i] = result.value;
        }

        Ok(results)
    }

    /// Find points within relevant range of the query point
    ///
    /// Returns indices of points to use and their distances to the query point
    fn find_relevant_points(&self, x: &ArrayView1<F>) -> InterpolateResult<(Vec<usize>, Vec<F>)> {
        let n_points = self.points.shape()[0];
        let n_dims = self.points.shape()[1];

        // Compute squared distances
        let mut distances = Vec::with_capacity(n_points);
        for i in 0..n_points {
            let mut d_squared = F::zero();
            for j in 0..n_dims {
                let diff = x[j] - self.points[[i, j]];
                d_squared = d_squared + diff * diff;
            }
            let dist = d_squared.sqrt();
            distances.push((i, dist));
        }

        // Sort by distance
        distances.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));

        // Apply max_points limit if specified
        let limit = match self.config.max_points {
            Some(limit) => std::cmp::min(limit, n_points),
            None => n_points,
        };

        // Filter out points with zero weight (if using compactly supported weight function)
        let effective_radius = match self.config.weight_fn {
            WeightFunction::WendlandC2 | WeightFunction::CubicSpline => self.config.bandwidth,
            _ => F::infinity(),
        };

        let mut indices = Vec::new();
        let mut dist_values = Vec::new();

        for &(idx, dist) in distances.iter().take(limit) {
            if dist <= effective_radius {
                indices.push(idx);
                dist_values.push(dist);
            }
        }

        // Ensure we have enough points for the basis
        let n_dims = self.points.shape()[1];
        let min_points = match self.config.basis {
            PolynomialBasis::Constant => 1,
            PolynomialBasis::Linear => n_dims + 1,
            PolynomialBasis::Quadratic => ((n_dims + 1) * (n_dims + 2)) / 2,
        };

        if indices.len() < min_points {
            // If not enough points with compact support, take the closest ones
            indices = distances
                .iter()
                .take(min_points)
                .map(|&(idx, _)| idx)
                .collect();
            dist_values = distances
                .iter()
                .take(min_points)
                .map(|&(_, dist)| dist)
                .collect();
        }

        Ok((indices, dist_values))
    }

    /// Compute weights based on distances and the weight function
    fn compute_weights(&self, distances: &[F]) -> InterpolateResult<Array1<F>> {
        let n = distances.len();
        let mut weights = Array1::zeros(n);

        for (i, &d) in distances.iter().enumerate() {
            // Normalize distance by bandwidth
            let r = d / self.config.bandwidth;

            // Compute weight based on the chosen weight function
            let weight = match self.config.weight_fn {
                WeightFunction::Gaussian => (-r * r).exp(),
                WeightFunction::WendlandC2 => {
                    if r < F::one() {
                        let t = F::one() - r;
                        let factor = F::from_f64(4.0).expect("Operation failed") * r + F::one();
                        t.powi(4) * factor
                    } else {
                        F::zero()
                    }
                }
                WeightFunction::InverseDistance => F::one() / (self.config.epsilon + r * r),
                WeightFunction::CubicSpline => {
                    if r < F::from_f64(1.0 / 3.0).expect("Operation failed") {
                        let r2 = r * r;
                        let r3 = r2 * r;
                        F::from_f64(2.0 / 3.0).expect("Operation failed")
                            - F::from_f64(9.0).expect("Operation failed") * r2
                            + F::from_f64(19.0).expect("Operation failed") * r3
                    } else if r < F::one() {
                        let t = F::from_f64(2.0).expect("Operation failed")
                            - F::from_f64(3.0).expect("Operation failed") * r;
                        F::from_f64(1.0 / 3.0).expect("Operation failed") * t.powi(3)
                    } else {
                        F::zero()
                    }
                }
            };

            weights[i] = weight;
        }

        // Normalize weights to sum to 1 for numerical stability
        let sum = weights.sum();
        if sum > F::zero() {
            weights.mapv_inplace(|w| w / sum);
        } else {
            // If all weights are zero (shouldn't happen), use equal weights
            weights.fill(F::from_f64(1.0 / n as f64).expect("Operation failed"));
        }

        Ok(weights)
    }

    /// Extract local points based on indices
    fn extract_local_points(&self, indices: &[usize]) -> Array2<F> {
        let n_points = indices.len();
        let n_dims = self.points.shape()[1];
        let mut local_points = Array2::zeros((n_points, n_dims));

        for (i, &idx) in indices.iter().enumerate() {
            let row = self.points.row(idx);
            local_points.row_mut(i).assign(&row);
        }

        local_points
    }

    /// Extract local values based on indices
    fn extract_local_values(&self, indices: &[usize]) -> Array1<F> {
        let mut local_values = Array1::zeros(indices.len());

        for (i, &idx) in indices.iter().enumerate() {
            local_values[i] = self.values[idx];
        }

        local_values
    }

    /// Create basis functions for the given local points
    fn create_basis_functions(
        &self,
        local_points: &Array2<F>,
        x: &ArrayView1<F>,
    ) -> InterpolateResult<Array2<F>> {
        let n_points = local_points.shape()[0];
        let n_dims = local_points.shape()[1];

        // Determine number of basis functions
        let n_basis = match self.config.basis {
            PolynomialBasis::Constant => 1,
            PolynomialBasis::Linear => n_dims + 1,
            PolynomialBasis::Quadratic => ((n_dims + 1) * (n_dims + 2)) / 2,
        };

        let mut basis = Array2::zeros((n_points, n_basis));

        // Fill basis functions for each point
        for i in 0..n_points {
            let point = local_points.row(i);
            let mut col = 0;

            // Constant term
            basis[[i, col]] = F::one();
            col += 1;

            if self.config.basis == PolynomialBasis::Linear
                || self.config.basis == PolynomialBasis::Quadratic
            {
                // Linear terms
                for j in 0..n_dims {
                    basis[[i, col]] = point[j] - x[j]; // Centered at query point for numerical stability
                    col += 1;
                }
            }

            if self.config.basis == PolynomialBasis::Quadratic {
                // Quadratic terms
                for j in 0..n_dims {
                    for k in j..n_dims {
                        let term_j = point[j] - x[j];
                        let term_k = point[k] - x[k];
                        basis[[i, col]] = term_j * term_k;
                        col += 1;
                    }
                }
            }
        }

        Ok(basis)
    }

    /// Perform weighted least squares fit and compute diagnostics
    #[allow(clippy::too_many_lines)]
    fn fit_weighted_least_squares(
        &self,
        local_points: &Array2<F>,
        local_values: &Array1<F>,
        _x: &ArrayView1<F>,
        weights: &Array1<F>,
        basis: &Array2<F>,
    ) -> InterpolateResult<RegressionResult<F>> {
        let n_points = local_points.shape()[0];
        let n_basis = basis.shape()[1];

        // Create the weighted basis matrix and target vector
        let mut w_basis = Array2::zeros((n_points, n_basis));
        let mut w_values = Array1::zeros(n_points);

        for i in 0..n_points {
            let sqrt_w = weights[i].sqrt();
            for j in 0..n_basis {
                w_basis[[i, j]] = basis[[i, j]] * sqrt_w;
            }
            w_values[i] = local_values[i] * sqrt_w;
        }

        // Solve the least squares problem: (X'WX)β = X'Wy
        #[cfg(feature = "linalg")]
        let xtx = w_basis.t().dot(&w_basis);
        #[cfg(not(feature = "linalg"))]
        let _xtx = w_basis.t().dot(&w_basis);
        let xty = w_basis.t().dot(&w_values);

        // Compute the hat matrix diagonal (leverage values)

        // Solve the system for coefficients
        #[cfg(feature = "linalg")]
        let coefficients = {
            use scirs2_linalg::solve;
            let xtx_f64 = xtx.mapv(|x| x.to_f64().expect("Operation failed"));
            let xty_f64 = xty.mapv(|x| x.to_f64().expect("Operation failed"));
            match solve(&xtx_f64.view(), &xty_f64.view(), None) {
                Ok(c) => c.mapv(|_x| F::from_f64(_x).expect("Operation failed")),
                Err(_) => {
                    // Fallback: use local weighted mean for numerical stability
                    let mut mean = F::zero();
                    let mut sum_weights = F::zero();

                    for i in 0..n_points {
                        mean = mean + weights[i] * local_values[i];
                        sum_weights = sum_weights + weights[i];
                    }

                    if sum_weights > F::zero() {
                        mean = mean / sum_weights;
                    } else {
                        mean = local_values.mean_or(F::zero());
                    }

                    let mut result = Array1::zeros(n_basis);
                    result[0] = mean;
                    result
                }
            }
        };

        #[cfg(not(feature = "linalg"))]
        let coefficients = {
            // Fallback implementation when linalg is not available
            // Simple diagonal approximation
            let mut result = Array1::zeros(xty.len());
            // For simple approximation, just use the weighted mean
            let mut mean = F::zero();
            let mut sum_weights = F::zero();

            for i in 0..n_points {
                mean = mean + weights[i] * local_values[i];
                sum_weights = sum_weights + weights[i];
            }

            if sum_weights > F::zero() {
                mean = mean / sum_weights;
            } else {
                mean = local_values.mean_or(F::zero());
            }

            result[0] = mean;
            result
        };

        // The fitted value at the query point is just the constant term coefficient
        // since all other basis functions (linear, quadratic) evaluate to 0 at the query point
        let fitted_value = coefficients[0];

        // H = X(X'WX)^(-1)X'W
        // We only need the diagonal of H for standard errors

        // Try to compute the inverse of X'WX
        #[cfg(feature = "linalg")]
        let xtx_inv = {
            use scirs2_linalg::inv;
            let xtx_f64 = xtx.mapv(|x| x.to_f64().expect("Operation failed"));
            match inv(&xtx_f64.view(), None) {
                Ok(inv) => inv.mapv(|_x| F::from_f64(_x).expect("Operation failed")),
                Err(_) => {
                    // If inversion fails, return a simpler result without diagnostics
                    return Ok(RegressionResult {
                        value: fitted_value,
                        std_error: F::zero(),
                        confidence_interval: None,
                        coefficients,
                        effective_df: F::from_f64(1.0).expect("Operation failed"),
                        r_squared: F::zero(),
                    });
                }
            }
        };

        #[cfg(not(feature = "linalg"))]
        {
            // Without linalg, return a simpler result
            Ok(RegressionResult {
                value: fitted_value,
                std_error: F::zero(),
                confidence_interval: None,
                coefficients,
                effective_df: F::from_f64(1.0).expect("Operation failed"),
                r_squared: F::zero(),
            })
        }

        // The following code only runs when linalg feature is enabled
        #[cfg(feature = "linalg")]
        {
            // Compute fitted _values for all local _points
            let fitted_local = basis.dot(&coefficients);

            // Compute residuals
            let residuals = local_values - &fitted_local;

            // Compute sum of squared residuals
            let ssr = residuals
                .iter()
                .zip(weights.iter())
                .fold(F::zero(), |acc, (&r, &w)| acc + w * r * r);

            // Compute total sum of squares
            let mean = local_values
                .iter()
                .zip(weights.iter())
                .fold(F::zero(), |acc, (&y, &w)| acc + w * y)
                / weights.sum();

            let sst = local_values
                .iter()
                .zip(weights.iter())
                .fold(F::zero(), |acc, (&y, &w)| acc + w * (y - mean) * (y - mean));

            // Compute R-squared
            let r_squared = if sst > F::zero() {
                F::one() - (ssr / sst)
            } else {
                F::zero()
            };

            // Compute leverage _values (diagonal of hat matrix)
            let mut leverage = Array1::zeros(n_points);
            for i in 0..n_points {
                let w_row = w_basis.row(i);
                let h_ii = w_row.dot(&xtx_inv.dot(&w_row));
                leverage[i] = h_ii;
            }

            // Compute effective degrees of freedom (sum of leverage values)
            let effective_df = leverage.sum();

            // Compute standard error of the fitted value

            // Get the first row of (X'WX)^(-1) which corresponds to the intercept
            let xtx_inv_row1 = xtx_inv.row(0);

            // Compute variance of residuals (MSE)
            let n_effective = F::from_usize(n_points).expect("Operation failed") - effective_df;
            let mse = if n_effective > F::zero() {
                ssr / n_effective
            } else {
                F::zero()
            };

            // Standard error depends on the type (robust or regular)
            let std_error = if self.config.robust_se {
                // Compute robust HC3 standard errors
                let mut sum_squared_weighted_residuals = F::zero();

                for i in 0..n_points {
                    // Adjust residuals by leverage: r_i / (1 - h_ii)
                    let adjusted_residual = if leverage[i] < F::one() {
                        residuals[i] / (F::one() - leverage[i])
                    } else {
                        residuals[i]
                    };

                    // Sum of squared weighted residuals
                    let weighted_residual = basis[[i, 0]] * adjusted_residual;
                    sum_squared_weighted_residuals =
                        sum_squared_weighted_residuals + weighted_residual * weighted_residual;
                }

                // Robust variance estimate
                (xtx_inv_row1[0] * sum_squared_weighted_residuals).sqrt()
            } else {
                // Regular standard error
                (xtx_inv_row1[0] * mse).sqrt()
            };

            // Compute confidence interval if requested
            let confidence_interval = self.config.confidence_level.map(|level| {
                // Get t-critical value (approximate with normal distribution for simplicity)
                let alpha = F::one() - level;
                let half_alpha = alpha / F::from_f64(2.0).expect("Operation failed");

                // Approximate critical value assuming a normal distribution
                // More accurate would be a t-distribution with n_effective degrees of freedom
                let z_critical = if half_alpha <= F::from_f64(0.001).expect("Operation failed") {
                    F::from_f64(3.09).expect("Operation failed") // ~99.9% CI
                } else if half_alpha <= F::from_f64(0.005).expect("Operation failed") {
                    F::from_f64(2.81).expect("Operation failed") // ~99% CI
                } else if half_alpha <= F::from_f64(0.01).expect("Operation failed") {
                    F::from_f64(2.58).expect("Operation failed") // ~98% CI
                } else if half_alpha <= F::from_f64(0.025).expect("Operation failed") {
                    F::from_f64(1.96).expect("Operation failed") // ~95% CI
                } else if half_alpha <= F::from_f64(0.05).expect("Operation failed") {
                    F::from_f64(1.645).expect("Operation failed") // ~90% CI
                } else {
                    F::from_f64(1.28).expect("Operation failed") // ~80% CI
                };

                let margin = z_critical * std_error;
                (fitted_value - margin, fitted_value + margin)
            });

            Ok(RegressionResult {
                value: fitted_value,
                std_error,
                confidence_interval,
                coefficients,
                effective_df,
                r_squared,
            })
        }
    }

    /// Get the configuration used by this local polynomial regression
    pub fn config(&self) -> &LocalPolynomialConfig<F> {
        &self.config
    }

    /// Get the points used by this local polynomial regression
    pub fn points(&self) -> &Array2<F> {
        &self.points
    }

    /// Get the values used by this local polynomial regression
    pub fn values(&self) -> &Array1<F> {
        &self.values
    }

    /// Compute the optimal bandwidth using leave-one-out cross-validation
    ///
    /// # Arguments
    ///
    /// * `bandwidths` - Array of bandwidths to evaluate
    ///
    /// # Returns
    ///
    /// The bandwidth that minimizes prediction error in cross-validation
    pub fn select_bandwidth(&self, bandwidths: &[F]) -> InterpolateResult<F> {
        if bandwidths.is_empty() {
            return Err(InterpolateError::InvalidValue(
                "Bandwidths array cannot be empty".to_string(),
            ));
        }

        let n_points = self.points.shape()[0];
        let mut best_bandwidth = bandwidths[0];
        let mut min_error = F::infinity();

        for &bandwidth in bandwidths {
            if bandwidth <= F::zero() {
                continue;
            }

            // Create a new model with this bandwidth
            let config = LocalPolynomialConfig {
                bandwidth,
                ..self.config.clone()
            };

            let model = Self::with_config(self.points.clone(), self.values.clone(), config)?;

            // Perform leave-one-out cross-validation
            let mut total_squared_error = F::zero();

            for i in 0..n_points {
                let point = self.points.row(i).to_owned();

                // Make prediction at this point (the model will automatically
                // not use the point itself since it's being left out)
                let result = model.fit_at_point(&point.view())?;

                // Compute squared error
                let error = result.value - self.values[i];
                total_squared_error = total_squared_error + error * error;
            }

            // Compute mean squared error
            let mse = total_squared_error / F::from_usize(n_points).expect("Operation failed");

            // Update best bandwidth if this one is better
            if mse < min_error {
                min_error = mse;
                best_bandwidth = bandwidth;
            }
        }

        Ok(best_bandwidth)
    }
}

/// Create a local polynomial regression model (LOESS)
///
/// This is a convenience function to create a local polynomial regression
/// model with commonly used default settings (Gaussian weights, linear fit).
///
/// # Arguments
///
/// * `points` - Point coordinates with shape (n_points, n_dims)
/// * `values` - Values at each point with shape (n_points,)
/// * `bandwidth` - Bandwidth parameter controlling locality (larger = smoother)
///
/// # Returns
///
/// A LocalPolynomialRegression model with common default settings
#[allow(dead_code)]
pub fn make_loess<F>(
    points: Array2<F>,
    values: Array1<F>,
    bandwidth: F,
) -> InterpolateResult<LocalPolynomialRegression<F>>
where
    F: Float + FromPrimitive + Debug + 'static,
{
    LocalPolynomialRegression::new(points, values, bandwidth)
}

/// Create a robust local polynomial regression model
///
/// This model uses robust standard errors and is less sensitive to outliers.
///
/// # Arguments
///
/// * `points` - Point coordinates with shape (n_points, n_dims)
/// * `values` - Values at each point with shape (n_points,)
/// * `bandwidth` - Bandwidth parameter controlling locality (larger = smoother)
/// * `confidence_level` - Confidence level for intervals (e.g., 0.95)
///
/// # Returns
///
/// A LocalPolynomialRegression model with robust error estimates
#[allow(dead_code)]
pub fn make_robust_loess<F>(
    points: Array2<F>,
    values: Array1<F>,
    bandwidth: F,
    confidence_level: F,
) -> InterpolateResult<LocalPolynomialRegression<F>>
where
    F: Float + FromPrimitive + Debug + 'static,
{
    let config = LocalPolynomialConfig {
        bandwidth,
        weight_fn: WeightFunction::Gaussian,
        basis: PolynomialBasis::Linear,
        confidence_level: Some(confidence_level),
        robust_se: true,
        max_points: None,
        epsilon: F::from_f64(1e-10).expect("Operation failed"),
    };

    LocalPolynomialRegression::with_config(points, values, config)
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_abs_diff_eq;
    use scirs2_core::ndarray::{array, Axis};

    #[test]
    fn test_local_polynomial_regression() {
        // Simple test with 1D data
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0].insert_axis(Axis(1));
        let y = array![0.0, 1.0, 4.0, 9.0, 16.0, 25.0]; // y = x^2

        let loess = LocalPolynomialRegression::new(
            x.clone(),
            y.clone(),
            1.5, // bandwidth
        )
        .expect("Operation failed");

        // Test at x = 2.5 (should be close to 6.25)
        let query = array![2.5];
        let result = loess.fit_at_point(&query.view()).expect("Operation failed");

        assert_abs_diff_eq!(result.value, 6.25, epsilon = 1.5);
    }

    #[test]
    fn test_confidence_intervals() {
        // Simple test with 1D data
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0].insert_axis(Axis(1));
        let y = array![0.0, 1.0, 4.0, 9.0, 16.0, 25.0]; // y = x^2

        let config = LocalPolynomialConfig {
            bandwidth: 1.5,
            confidence_level: Some(0.95),
            ..LocalPolynomialConfig::default()
        };

        let loess = LocalPolynomialRegression::with_config(x.clone(), y.clone(), config)
            .expect("Operation failed");

        // Test at x = 2.5
        let query = array![2.5];
        let result = loess.fit_at_point(&query.view()).expect("Operation failed");

        // Confidence interval should exist and contain the true value (6.25)
        // Note: confidence intervals require the linalg feature
        #[cfg(feature = "linalg")]
        {
            let (lower, upper) = result.confidence_interval.expect("Operation failed");
            assert!(lower < 6.25);
            assert!(upper > 6.25);
        }
        #[cfg(not(feature = "linalg"))]
        {
            assert!(result.confidence_interval.is_none());
        }
    }
}