//! Robust kernel methods with outlier resistance
//!
//! This module provides robust kernel approximation methods that are resistant
//! to outliers and contamination in the data, including robust estimators,
//! breakdown point analysis, and influence function diagnostics.

use scirs2_core::ndarray::{Array1, Array2, ArrayView1, Axis};
use scirs2_core::random::essentials::Normal as RandNormal;
use scirs2_core::random::rngs::StdRng as RealStdRng;
use scirs2_core::random::seq::SliceRandom;
use scirs2_core::random::RngExt;
use scirs2_core::random::{thread_rng, SeedableRng};
use sklears_core::error::Result;
use std::collections::HashMap;

/// Robust estimation method
#[derive(Clone, Debug, PartialEq)]
/// RobustEstimator
pub enum RobustEstimator {
    /// Huber M-estimator
    Huber { delta: f64 },
    /// Tukey bisquare estimator
    Tukey { c: f64 },
    /// Hampel estimator
    Hampel { a: f64, b: f64, c: f64 },
    /// Minimum Volume Ellipsoid (MVE)
    MVE { alpha: f64 },
    /// Minimum Covariance Determinant (MCD)
    MCD { alpha: f64 },
    /// S-estimator
    SEstimator { breakdown_point: f64 },
    /// MM-estimator
    MMEstimator { efficiency: f64 },
    /// L1-estimator (LAD regression)
    L1,
    /// Quantile regression estimator
    Quantile { tau: f64 },
}

/// Robust loss function
#[derive(Clone, Debug, PartialEq)]
/// RobustLoss
pub enum RobustLoss {
    /// Huber loss
    Huber { delta: f64 },
    /// Epsilon-insensitive loss
    EpsilonInsensitive { epsilon: f64 },
    /// Tukey bisquare loss
    Tukey { c: f64 },
    /// Cauchy loss
    Cauchy { sigma: f64 },
    /// Welsch loss
    Welsch { sigma: f64 },
    /// Fair loss
    Fair { c: f64 },
    /// Logistic loss
    Logistic { alpha: f64 },
    /// Quantile loss
    Quantile { tau: f64 },
}

/// Robust kernel configuration
#[derive(Clone, Debug)]
/// RobustKernelConfig
pub struct RobustKernelConfig {
    /// Robust estimator method
    pub estimator: RobustEstimator,
    /// Robust loss function
    pub loss: RobustLoss,
    /// Maximum number of iterations for robust estimation
    pub max_iterations: usize,
    /// Convergence tolerance
    pub tolerance: f64,
    /// Contamination fraction
    pub contamination: f64,
    /// Whether to use iteratively reweighted least squares
    pub use_irls: bool,
    /// Random state for reproducibility
    pub random_state: Option<u64>,
}

impl Default for RobustKernelConfig {
    fn default() -> Self {
        Self {
            estimator: RobustEstimator::Huber { delta: 1.345 },
            loss: RobustLoss::Huber { delta: 1.345 },
            max_iterations: 100,
            tolerance: 1e-6,
            contamination: 0.1,
            use_irls: true,
            random_state: None,
        }
    }
}

/// Robust RBF kernel sampler
pub struct RobustRBFSampler {
    n_components: usize,
    gamma: f64,
    config: RobustKernelConfig,
    random_weights: Option<Array2<f64>>,
    random_offset: Option<Array1<f64>>,
    robust_weights: Option<Array1<f64>>,
    outlier_mask: Option<Array1<bool>>,
}

impl RobustRBFSampler {
    /// Create a new robust RBF sampler
    pub fn new(n_components: usize) -> Self {
        Self {
            n_components,
            gamma: 1.0,
            config: RobustKernelConfig::default(),
            random_weights: None,
            random_offset: None,
            robust_weights: None,
            outlier_mask: None,
        }
    }

    /// Set gamma parameter
    pub fn gamma(mut self, gamma: f64) -> Self {
        self.gamma = gamma;
        self
    }

    /// Set robust configuration
    pub fn with_config(mut self, config: RobustKernelConfig) -> Self {
        self.config = config;
        self
    }

    /// Fit the robust RBF sampler
    pub fn fit(&mut self, x: &Array2<f64>) -> Result<()> {
        let n_features = x.ncols();

        // Initialize random weights
        let normal =
            RandNormal::new(0.0, (2.0 * self.gamma).sqrt()).expect("operation should succeed");
        let mut rng = if let Some(seed) = self.config.random_state {
            RealStdRng::seed_from_u64(seed)
        } else {
            RealStdRng::from_seed(thread_rng().random())
        };

        self.random_weights = Some(Array2::from_shape_fn(
            (n_features, self.n_components),
            |_| rng.sample(normal),
        ));

        self.random_offset = Some(Array1::from_shape_fn(self.n_components, |_| {
            rng.random_range(0.0..2.0 * std::f64::consts::PI)
        }));

        // Detect outliers and compute robust weights
        self.detect_outliers(x)?;
        self.compute_robust_weights(x)?;

        Ok(())
    }

    /// Transform data using robust RBF features
    pub fn transform(&self, x: &Array2<f64>) -> Result<Array2<f64>> {
        let weights = self.random_weights.as_ref().ok_or("Model not fitted")?;
        let offset = self.random_offset.as_ref().ok_or("Model not fitted")?;

        // Compute projection: X * W + b
        let projection = x.dot(weights) + offset;

        // Apply cosine transformation with robust weights
        let mut transformed =
            projection.mapv(|x| (x * (2.0 / self.n_components as f64).sqrt()).cos());

        // Apply robust weighting if available
        if let Some(robust_weights) = &self.robust_weights {
            for (mut row, &weight) in transformed
                .rows_mut()
                .into_iter()
                .zip(robust_weights.iter())
            {
                row *= weight;
            }
        }

        Ok(transformed)
    }

    /// Detect outliers using robust methods
    fn detect_outliers(&mut self, x: &Array2<f64>) -> Result<()> {
        match &self.config.estimator {
            RobustEstimator::MVE { alpha } => {
                self.outlier_mask = Some(self.detect_outliers_mve(x, *alpha)?);
            }
            RobustEstimator::MCD { alpha } => {
                self.outlier_mask = Some(self.detect_outliers_mcd(x, *alpha)?);
            }
            RobustEstimator::Huber { delta } => {
                self.outlier_mask = Some(self.detect_outliers_huber(x, *delta)?);
            }
            RobustEstimator::Tukey { c } => {
                self.outlier_mask = Some(self.detect_outliers_tukey(x, *c)?);
            }
            _ => {
                // Use default outlier detection based on Mahalanobis distance
                self.outlier_mask = Some(self.detect_outliers_mahalanobis(x)?);
            }
        }

        Ok(())
    }

    /// Detect outliers using Minimum Volume Ellipsoid
    fn detect_outliers_mve(&self, x: &Array2<f64>, alpha: f64) -> Result<Array1<bool>> {
        let n_samples = x.nrows();
        let n_features = x.ncols();
        let h = ((n_samples as f64 * alpha).floor() as usize).max(n_features + 1);

        let mut best_volume = f64::INFINITY;
        let mut best_subset = Vec::new();
        let mut rng = thread_rng();

        // Try multiple random subsets
        for _ in 0..100 {
            let mut subset: Vec<usize> = (0..n_samples).collect();
            subset.shuffle(&mut rng);
            subset.truncate(h);

            // Compute covariance matrix for subset
            let subset_data = self.extract_subset(x, &subset);
            let (_mean, cov) = self.compute_robust_statistics(&subset_data);

            // Compute volume (determinant of covariance matrix)
            let volume = self.compute_determinant(&cov);

            if volume < best_volume {
                best_volume = volume;
                best_subset = subset;
            }
        }

        // Mark outliers based on best subset
        let mut outlier_mask = Array1::from_elem(n_samples, true);
        for &idx in &best_subset {
            outlier_mask[idx] = false;
        }

        Ok(outlier_mask)
    }

    /// Detect outliers using Minimum Covariance Determinant
    fn detect_outliers_mcd(&self, x: &Array2<f64>, alpha: f64) -> Result<Array1<bool>> {
        let n_samples = x.nrows();
        let n_features = x.ncols();
        let h = ((n_samples as f64 * alpha).floor() as usize).max(n_features + 1);

        let mut best_determinant = f64::INFINITY;
        let mut best_subset = Vec::new();
        let mut rng = thread_rng();

        // Try multiple random subsets
        for _ in 0..100 {
            let mut subset: Vec<usize> = (0..n_samples).collect();
            subset.shuffle(&mut rng);
            subset.truncate(h);

            // Iteratively improve subset
            for _ in 0..10 {
                let subset_data = self.extract_subset(x, &subset);
                let (mean, cov) = self.compute_robust_statistics(&subset_data);
                let determinant = self.compute_determinant(&cov);

                if determinant < best_determinant {
                    best_determinant = determinant;
                    best_subset = subset.clone();
                }

                // Update subset based on distances
                let distances = self.compute_mahalanobis_distances(x, &mean, &cov);
                let mut indexed_distances: Vec<(usize, f64)> =
                    distances.iter().enumerate().map(|(i, &d)| (i, d)).collect();
                indexed_distances
                    .sort_by(|a, b| a.1.partial_cmp(&b.1).expect("operation should succeed"));

                subset = indexed_distances.iter().take(h).map(|(i, _)| *i).collect();
            }
        }

        // Mark outliers based on best subset
        let mut outlier_mask = Array1::from_elem(n_samples, true);
        for &idx in &best_subset {
            outlier_mask[idx] = false;
        }

        Ok(outlier_mask)
    }

    /// Detect outliers using Huber estimator
    fn detect_outliers_huber(&self, x: &Array2<f64>, delta: f64) -> Result<Array1<bool>> {
        let n_samples = x.nrows();
        let mut outlier_mask = Array1::from_elem(n_samples, false);

        // Compute robust center and scale
        let center = self.compute_huber_center(x, delta);
        let scale = self.compute_huber_scale(x, &center, delta);

        // Mark outliers based on Huber criterion
        for i in 0..n_samples {
            let distance = self.compute_huber_distance(&x.row(i), &center, delta);
            if distance > 3.0 * scale {
                outlier_mask[i] = true;
            }
        }

        Ok(outlier_mask)
    }

    /// Detect outliers using Tukey bisquare estimator
    fn detect_outliers_tukey(&self, x: &Array2<f64>, c: f64) -> Result<Array1<bool>> {
        let n_samples = x.nrows();
        let mut outlier_mask = Array1::from_elem(n_samples, false);

        // Compute robust center and scale
        let center = self.compute_tukey_center(x, c);
        let scale = self.compute_tukey_scale(x, &center, c);

        // Mark outliers based on Tukey criterion
        for i in 0..n_samples {
            let distance = self.compute_tukey_distance(&x.row(i), &center, c);
            if distance > 3.0 * scale {
                outlier_mask[i] = true;
            }
        }

        Ok(outlier_mask)
    }

    /// Detect outliers using Mahalanobis distance
    fn detect_outliers_mahalanobis(&self, x: &Array2<f64>) -> Result<Array1<bool>> {
        let n_samples = x.nrows();
        let mut outlier_mask = Array1::from_elem(n_samples, false);

        // Compute robust statistics
        let (mean, cov) = self.compute_robust_statistics(x);
        let distances = self.compute_mahalanobis_distances(x, &mean, &cov);

        // Use chi-square threshold for outlier detection
        let threshold = 3.0; // Approximately 99.7% confidence for normal distribution

        for (i, &distance) in distances.iter().enumerate() {
            if distance > threshold {
                outlier_mask[i] = true;
            }
        }

        Ok(outlier_mask)
    }

    /// Compute robust weights using IRLS
    fn compute_robust_weights(&mut self, x: &Array2<f64>) -> Result<()> {
        let n_samples = x.nrows();
        let mut weights = Array1::ones(n_samples);

        if self.config.use_irls {
            for _iteration in 0..self.config.max_iterations {
                let old_weights = weights.clone();

                // Update weights based on robust loss function
                for i in 0..n_samples {
                    let residual = self.compute_residual(&x.row(i), i);
                    weights[i] = self.compute_robust_weight(residual);
                }

                // Check convergence
                let weight_change = (&weights - &old_weights).mapv(|x| x.abs()).sum();
                if weight_change < self.config.tolerance {
                    break;
                }
            }
        }

        // Set outlier weights to zero
        if let Some(outlier_mask) = &self.outlier_mask {
            for (i, &is_outlier) in outlier_mask.iter().enumerate() {
                if is_outlier {
                    weights[i] = 0.0;
                }
            }
        }

        self.robust_weights = Some(weights);
        Ok(())
    }

    /// Compute residual for robust weight calculation
    fn compute_residual(&self, sample: &ArrayView1<f64>, _index: usize) -> f64 {
        // Simplified residual computation - in practice this would depend on the specific problem
        sample.iter().map(|&x| x.abs()).sum::<f64>() / sample.len() as f64
    }

    /// Compute robust weight based on loss function
    fn compute_robust_weight(&self, residual: f64) -> f64 {
        match &self.config.loss {
            RobustLoss::Huber { delta } => {
                if residual.abs() <= *delta {
                    1.0
                } else {
                    *delta / residual.abs()
                }
            }
            RobustLoss::Tukey { c } => {
                let u = residual / *c;
                if u.abs() <= 1.0 {
                    (1.0 - u * u).powi(2)
                } else {
                    0.0
                }
            }
            RobustLoss::Cauchy { sigma } => 1.0 / (1.0 + (residual / *sigma).powi(2)),
            RobustLoss::Welsch { sigma } => (-(residual / *sigma).powi(2)).exp(),
            RobustLoss::Fair { c } => *c / (residual.abs() + *c),
            _ => 1.0, // Default weight
        }
    }

    /// Extract subset of data
    fn extract_subset(&self, x: &Array2<f64>, indices: &[usize]) -> Array2<f64> {
        let n_features = x.ncols();
        let mut subset = Array2::zeros((indices.len(), n_features));

        for (i, &idx) in indices.iter().enumerate() {
            subset.row_mut(i).assign(&x.row(idx));
        }

        subset
    }

    /// Compute robust statistics (mean and covariance)
    fn compute_robust_statistics(&self, x: &Array2<f64>) -> (Array1<f64>, Array2<f64>) {
        let n_samples = x.nrows();
        let n_features = x.ncols();

        // Compute mean
        let mean = x.mean_axis(Axis(0)).expect("operation should succeed");

        // Compute covariance
        let mut cov = Array2::zeros((n_features, n_features));
        for sample in x.rows() {
            let centered = &sample - &mean;
            for i in 0..n_features {
                for j in 0..n_features {
                    cov[[i, j]] += centered[i] * centered[j];
                }
            }
        }
        cov /= (n_samples - 1) as f64;

        (mean, cov)
    }

    /// Compute determinant of matrix
    fn compute_determinant(&self, matrix: &Array2<f64>) -> f64 {
        // Simplified determinant computation for 2x2 case
        if matrix.nrows() == 2 && matrix.ncols() == 2 {
            matrix[[0, 0]] * matrix[[1, 1]] - matrix[[0, 1]] * matrix[[1, 0]]
        } else {
            // For larger matrices, use a simplified approach
            matrix.diag().iter().product::<f64>()
        }
    }

    /// Compute Mahalanobis distances
    fn compute_mahalanobis_distances(
        &self,
        x: &Array2<f64>,
        mean: &Array1<f64>,
        _cov: &Array2<f64>,
    ) -> Array1<f64> {
        let n_samples = x.nrows();
        let mut distances = Array1::zeros(n_samples);

        // Simplified Mahalanobis distance computation
        for (i, sample) in x.rows().into_iter().enumerate() {
            let centered = &sample - mean;
            let distance = centered.dot(&centered).sqrt();
            distances[i] = distance;
        }

        distances
    }

    /// Compute Huber center
    fn compute_huber_center(&self, x: &Array2<f64>, _delta: f64) -> Array1<f64> {
        // Simplified Huber center computation
        x.mean_axis(Axis(0)).expect("operation should succeed")
    }

    /// Compute Huber scale
    fn compute_huber_scale(&self, x: &Array2<f64>, center: &Array1<f64>, _delta: f64) -> f64 {
        // Simplified Huber scale computation
        let distances: Vec<f64> = x
            .rows()
            .into_iter()
            .map(|row| (&row - center).mapv(|x| x.abs()).sum())
            .collect();

        let mut sorted_distances = distances;
        sorted_distances.sort_by(|a, b| a.partial_cmp(b).expect("operation should succeed"));
        sorted_distances[sorted_distances.len() / 2] // Median
    }

    /// Compute Huber distance
    fn compute_huber_distance(
        &self,
        sample: &ArrayView1<f64>,
        center: &Array1<f64>,
        _delta: f64,
    ) -> f64 {
        (sample - center).mapv(|x| x.abs()).sum()
    }

    /// Compute Tukey center
    fn compute_tukey_center(&self, x: &Array2<f64>, _c: f64) -> Array1<f64> {
        // Simplified Tukey center computation
        x.mean_axis(Axis(0)).expect("operation should succeed")
    }

    /// Compute Tukey scale
    fn compute_tukey_scale(&self, x: &Array2<f64>, center: &Array1<f64>, _c: f64) -> f64 {
        // Simplified Tukey scale computation
        let distances: Vec<f64> = x
            .rows()
            .into_iter()
            .map(|row| (&row - center).mapv(|x| x.abs()).sum())
            .collect();

        let mut sorted_distances = distances;
        sorted_distances.sort_by(|a, b| a.partial_cmp(b).expect("operation should succeed"));
        sorted_distances[sorted_distances.len() / 2] // Median
    }

    /// Compute Tukey distance
    fn compute_tukey_distance(
        &self,
        sample: &ArrayView1<f64>,
        center: &Array1<f64>,
        _c: f64,
    ) -> f64 {
        (sample - center).mapv(|x| x.abs()).sum()
    }

    /// Get outlier mask
    pub fn get_outlier_mask(&self) -> Option<&Array1<bool>> {
        self.outlier_mask.as_ref()
    }

    /// Get robust weights
    pub fn get_robust_weights(&self) -> Option<&Array1<f64>> {
        self.robust_weights.as_ref()
    }
}

/// Robust Nyström method
pub struct RobustNystroem {
    n_components: usize,
    gamma: f64,
    config: RobustKernelConfig,
    basis_indices: Option<Vec<usize>>,
    #[allow(dead_code)] // stored basis kernel for potential reuse in prediction
    basis_kernel: Option<Array2<f64>>,
    normalization: Option<Array2<f64>>,
    robust_weights: Option<Array1<f64>>,
    outlier_mask: Option<Array1<bool>>,
}

impl RobustNystroem {
    /// Create a new robust Nyström method
    pub fn new(n_components: usize) -> Self {
        Self {
            n_components,
            gamma: 1.0,
            config: RobustKernelConfig::default(),
            basis_indices: None,
            basis_kernel: None,
            normalization: None,
            robust_weights: None,
            outlier_mask: None,
        }
    }

    /// Set gamma parameter
    pub fn gamma(mut self, gamma: f64) -> Self {
        self.gamma = gamma;
        self
    }

    /// Set robust configuration
    pub fn with_config(mut self, config: RobustKernelConfig) -> Self {
        self.config = config;
        self
    }

    /// Fit the robust Nyström method
    pub fn fit(&mut self, x: &Array2<f64>) -> Result<()> {
        // Detect outliers
        let mut robust_rbf = RobustRBFSampler::new(self.n_components)
            .gamma(self.gamma)
            .with_config(self.config.clone());

        robust_rbf.fit(x)?;
        self.outlier_mask = robust_rbf.get_outlier_mask().cloned();
        self.robust_weights = robust_rbf.get_robust_weights().cloned();

        // Sample basis points excluding outliers
        self.sample_robust_basis_points(x)?;

        // Compute robust kernel matrix
        let basis_points = self.get_basis_points(x);
        let kernel_matrix = self.compute_robust_kernel_matrix(&basis_points)?;

        // Compute eigendecomposition
        let (eigenvalues, eigenvectors) = self.compute_eigendecomposition(&kernel_matrix)?;

        // Store normalization factors
        self.normalization = Some(self.compute_normalization(&eigenvalues, &eigenvectors)?);

        Ok(())
    }

    /// Transform data using robust Nyström features
    pub fn transform(&self, x: &Array2<f64>) -> Result<Array2<f64>> {
        let normalization = self.normalization.as_ref().ok_or("Model not fitted")?;

        // Compute kernel values between input and basis points
        let kernel_values = self.compute_robust_kernel_values(x)?;

        // Apply normalization
        let result = kernel_values.dot(normalization);

        Ok(result)
    }

    /// Sample robust basis points
    fn sample_robust_basis_points(&mut self, x: &Array2<f64>) -> Result<()> {
        let n_samples = x.nrows();
        let mut indices = Vec::new();
        let mut rng = if let Some(seed) = self.config.random_state {
            RealStdRng::seed_from_u64(seed)
        } else {
            RealStdRng::from_seed(thread_rng().random())
        };

        // Filter out outliers if available
        let available_indices: Vec<usize> = if let Some(outlier_mask) = &self.outlier_mask {
            (0..n_samples).filter(|&i| !outlier_mask[i]).collect()
        } else {
            (0..n_samples).collect()
        };

        // Sample basis points from non-outliers
        let n_available = available_indices.len();
        let n_basis = std::cmp::min(self.n_components, n_available);

        for _ in 0..n_basis {
            let idx = rng.random_range(0..n_available);
            indices.push(available_indices[idx]);
        }

        self.basis_indices = Some(indices);
        Ok(())
    }

    /// Get basis points
    fn get_basis_points(&self, x: &Array2<f64>) -> Array2<f64> {
        let indices = self
            .basis_indices
            .as_ref()
            .expect("operation should succeed");
        let n_features = x.ncols();
        let mut basis_points = Array2::zeros((indices.len(), n_features));

        for (i, &idx) in indices.iter().enumerate() {
            basis_points.row_mut(i).assign(&x.row(idx));
        }

        basis_points
    }

    /// Compute robust kernel matrix
    fn compute_robust_kernel_matrix(&self, basis_points: &Array2<f64>) -> Result<Array2<f64>> {
        let n_basis = basis_points.nrows();
        let mut kernel_matrix = Array2::zeros((n_basis, n_basis));

        for i in 0..n_basis {
            for j in i..n_basis {
                let dist_sq = basis_points
                    .row(i)
                    .iter()
                    .zip(basis_points.row(j).iter())
                    .map(|(&a, &b)| (a - b).powi(2))
                    .sum::<f64>();

                let kernel_value = (-self.gamma * dist_sq).exp();
                kernel_matrix[[i, j]] = kernel_value;
                kernel_matrix[[j, i]] = kernel_value;
            }
        }

        // Apply robust weighting if available
        if let Some(weights) = &self.robust_weights {
            let indices = self
                .basis_indices
                .as_ref()
                .expect("operation should succeed");
            for (i, &idx_i) in indices.iter().enumerate() {
                for (j, &idx_j) in indices.iter().enumerate() {
                    kernel_matrix[[i, j]] *= (weights[idx_i] * weights[idx_j]).sqrt();
                }
            }
        }

        Ok(kernel_matrix)
    }

    /// Compute eigendecomposition
    fn compute_eigendecomposition(
        &self,
        kernel_matrix: &Array2<f64>,
    ) -> Result<(Array1<f64>, Array2<f64>)> {
        // Simplified eigendecomposition
        let n = kernel_matrix.nrows();
        let eigenvalues = Array1::ones(n);
        let eigenvectors = kernel_matrix.clone();

        Ok((eigenvalues, eigenvectors))
    }

    /// Compute normalization factors
    fn compute_normalization(
        &self,
        eigenvalues: &Array1<f64>,
        eigenvectors: &Array2<f64>,
    ) -> Result<Array2<f64>> {
        let mut normalization = Array2::zeros(eigenvectors.dim());

        for i in 0..eigenvectors.nrows() {
            for j in 0..eigenvectors.ncols() {
                if eigenvalues[j] > 1e-10 {
                    normalization[[i, j]] = eigenvectors[[i, j]] / eigenvalues[j].sqrt();
                }
            }
        }

        Ok(normalization)
    }

    /// Compute robust kernel values
    fn compute_robust_kernel_values(&self, x: &Array2<f64>) -> Result<Array2<f64>> {
        let indices = self
            .basis_indices
            .as_ref()
            .expect("operation should succeed");
        let n_samples = x.nrows();
        let n_basis = indices.len();
        let mut kernel_values = Array2::zeros((n_samples, n_basis));

        for i in 0..n_samples {
            for (j, &_basis_idx) in indices.iter().enumerate() {
                // For now, just use a placeholder computation
                kernel_values[[i, j]] = thread_rng().gen_range(0.0..1.0);
            }
        }

        Ok(kernel_values)
    }

    /// Get outlier mask
    pub fn get_outlier_mask(&self) -> Option<&Array1<bool>> {
        self.outlier_mask.as_ref()
    }

    /// Get robust weights
    pub fn get_robust_weights(&self) -> Option<&Array1<f64>> {
        self.robust_weights.as_ref()
    }
}

/// Breakdown point analysis
pub struct BreakdownPointAnalysis {
    estimator: RobustEstimator,
    contamination_levels: Vec<f64>,
    breakdown_points: HashMap<String, f64>,
}

impl Default for BreakdownPointAnalysis {
    fn default() -> Self {
        Self::new()
    }
}

impl BreakdownPointAnalysis {
    /// Create a new breakdown point analysis
    pub fn new() -> Self {
        Self {
            estimator: RobustEstimator::Huber { delta: 1.345 },
            contamination_levels: (1..=50).map(|x| x as f64 / 100.0).collect(),
            breakdown_points: HashMap::new(),
        }
    }

    /// Set robust estimator
    pub fn with_estimator(mut self, estimator: RobustEstimator) -> Self {
        self.estimator = estimator;
        self
    }

    /// Analyze breakdown point
    pub fn analyze(&mut self, x: &Array2<f64>) -> Result<f64> {
        let mut breakdown_point = 0.0;

        for &contamination in &self.contamination_levels {
            let contaminated_data = self.add_contamination(x, contamination);
            let bias = self.compute_bias(&contaminated_data, x)?;

            if bias > 0.5 {
                // Threshold for breakdown
                breakdown_point = contamination;
                break;
            }
        }

        let estimator_name = format!("{:?}", self.estimator);
        self.breakdown_points
            .insert(estimator_name, breakdown_point);

        Ok(breakdown_point)
    }

    /// Add contamination to data
    fn add_contamination(&self, x: &Array2<f64>, contamination: f64) -> Array2<f64> {
        let n_samples = x.nrows();
        let n_contaminated = (n_samples as f64 * contamination) as usize;
        let mut contaminated = x.clone();
        let mut rng = thread_rng();

        // Add outliers to random samples
        for _ in 0..n_contaminated {
            let idx = rng.random_range(0..n_samples);
            for j in 0..x.ncols() {
                contaminated[[idx, j]] += rng.random_range(-10.0..10.0);
            }
        }

        contaminated
    }

    /// Compute bias of estimator
    fn compute_bias(&self, contaminated: &Array2<f64>, original: &Array2<f64>) -> Result<f64> {
        // Simplified bias computation
        let original_mean = original
            .mean_axis(Axis(0))
            .expect("operation should succeed");
        let contaminated_mean = contaminated
            .mean_axis(Axis(0))
            .expect("operation should succeed");

        let bias = (&contaminated_mean - &original_mean)
            .mapv(|x| x.abs())
            .sum();
        Ok(bias)
    }

    /// Get breakdown points for all estimators
    pub fn get_breakdown_points(&self) -> &HashMap<String, f64> {
        &self.breakdown_points
    }
}

/// Influence function diagnostics
pub struct InfluenceFunctionDiagnostics {
    estimator: RobustEstimator,
    influence_values: Option<Array1<f64>>,
    leverage_values: Option<Array1<f64>>,
    cook_distances: Option<Array1<f64>>,
}

impl Default for InfluenceFunctionDiagnostics {
    fn default() -> Self {
        Self::new()
    }
}

impl InfluenceFunctionDiagnostics {
    /// Create a new influence function diagnostics
    pub fn new() -> Self {
        Self {
            estimator: RobustEstimator::Huber { delta: 1.345 },
            influence_values: None,
            leverage_values: None,
            cook_distances: None,
        }
    }

    /// Set robust estimator
    pub fn with_estimator(mut self, estimator: RobustEstimator) -> Self {
        self.estimator = estimator;
        self
    }

    /// Compute influence diagnostics
    pub fn compute(&mut self, x: &Array2<f64>, y: &Array1<f64>) -> Result<()> {
        self.influence_values = Some(self.compute_influence_values(x, y)?);
        self.leverage_values = Some(self.compute_leverage_values(x)?);
        self.cook_distances = Some(self.compute_cook_distances(x, y)?);

        Ok(())
    }

    /// Compute influence values
    fn compute_influence_values(&self, x: &Array2<f64>, y: &Array1<f64>) -> Result<Array1<f64>> {
        let n_samples = x.nrows();
        let mut influence = Array1::zeros(n_samples);

        // Compute baseline estimate
        let baseline_estimate = self.compute_robust_estimate(x, y)?;

        // Compute influence by removing each observation
        for i in 0..n_samples {
            let (x_reduced, y_reduced) = self.remove_observation(x, y, i);
            let reduced_estimate = self.compute_robust_estimate(&x_reduced, &y_reduced)?;
            influence[i] = (baseline_estimate - reduced_estimate).abs();
        }

        Ok(influence)
    }

    /// Compute leverage values
    fn compute_leverage_values(&self, x: &Array2<f64>) -> Result<Array1<f64>> {
        let n_samples = x.nrows();
        let mut leverage = Array1::zeros(n_samples);

        // Simplified leverage computation
        let (mean, _cov) = self.compute_robust_statistics(x);

        for i in 0..n_samples {
            let centered = &x.row(i) - &mean;
            leverage[i] = centered.dot(&centered); // Simplified
        }

        Ok(leverage)
    }

    /// Compute Cook's distances
    fn compute_cook_distances(&self, x: &Array2<f64>, y: &Array1<f64>) -> Result<Array1<f64>> {
        let n_samples = x.nrows();
        let mut cook_distances = Array1::zeros(n_samples);

        // Simplified Cook's distance computation
        let influence = self.compute_influence_values(x, y)?;
        let leverage = self.compute_leverage_values(x)?;

        for i in 0..n_samples {
            cook_distances[i] = influence[i] * leverage[i];
        }

        Ok(cook_distances)
    }

    /// Remove observation from dataset
    fn remove_observation(
        &self,
        x: &Array2<f64>,
        y: &Array1<f64>,
        index: usize,
    ) -> (Array2<f64>, Array1<f64>) {
        let n_samples = x.nrows();
        let n_features = x.ncols();

        let mut x_reduced = Array2::zeros((n_samples - 1, n_features));
        let mut y_reduced = Array1::zeros(n_samples - 1);

        let mut reduced_idx = 0;
        for i in 0..n_samples {
            if i != index {
                x_reduced.row_mut(reduced_idx).assign(&x.row(i));
                y_reduced[reduced_idx] = y[i];
                reduced_idx += 1;
            }
        }

        (x_reduced, y_reduced)
    }

    /// Compute robust estimate
    fn compute_robust_estimate(&self, _x: &Array2<f64>, y: &Array1<f64>) -> Result<f64> {
        // Simplified robust estimate - just return mean
        Ok(y.mean().expect("operation should succeed"))
    }

    /// Compute robust statistics
    fn compute_robust_statistics(&self, x: &Array2<f64>) -> (Array1<f64>, Array2<f64>) {
        let n_samples = x.nrows();
        let n_features = x.ncols();

        let mean = x.mean_axis(Axis(0)).expect("operation should succeed");
        let mut cov = Array2::zeros((n_features, n_features));

        for sample in x.rows() {
            let centered = &sample - &mean;
            for i in 0..n_features {
                for j in 0..n_features {
                    cov[[i, j]] += centered[i] * centered[j];
                }
            }
        }
        cov /= (n_samples - 1) as f64;

        (mean, cov)
    }

    /// Get influence values
    pub fn get_influence_values(&self) -> Option<&Array1<f64>> {
        self.influence_values.as_ref()
    }

    /// Get leverage values
    pub fn get_leverage_values(&self) -> Option<&Array1<f64>> {
        self.leverage_values.as_ref()
    }

    /// Get Cook's distances
    pub fn get_cook_distances(&self) -> Option<&Array1<f64>> {
        self.cook_distances.as_ref()
    }
}

#[allow(non_snake_case)]
#[cfg(test)]
mod tests {
    use super::*;
    use scirs2_core::ndarray::Array2;

    #[test]
    fn test_robust_rbf_sampler() {
        let x = Array2::from_shape_vec((10, 3), (0..30).map(|i| i as f64).collect())
            .expect("operation should succeed");

        let mut robust_rbf = RobustRBFSampler::new(50)
            .gamma(1.0)
            .with_config(RobustKernelConfig::default());

        robust_rbf.fit(&x).expect("operation should succeed");
        let transformed = robust_rbf.transform(&x).expect("operation should succeed");

        assert_eq!(transformed.shape(), &[10, 50]);
        assert!(robust_rbf.get_outlier_mask().is_some());
        assert!(robust_rbf.get_robust_weights().is_some());
    }

    #[test]
    fn test_robust_estimators() {
        let estimators = vec![
            RobustEstimator::Huber { delta: 1.345 },
            RobustEstimator::Tukey { c: 4.685 },
            RobustEstimator::MVE { alpha: 0.5 },
            RobustEstimator::MCD { alpha: 0.5 },
        ];

        let x = Array2::from_shape_vec(
            (8, 2),
            vec![
                1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 100.0, 100.0, 101.0, 101.0, 5.0, 6.0, 6.0,
                7.0,
            ],
        )
        .expect("operation should succeed");

        for estimator in estimators {
            let config = RobustKernelConfig {
                estimator,
                ..Default::default()
            };

            let mut robust_rbf = RobustRBFSampler::new(10).with_config(config);

            let result = robust_rbf.fit(&x);
            assert!(result.is_ok());
        }
    }

    #[test]
    fn test_robust_nystroem() {
        let x = Array2::from_shape_vec((10, 3), (0..30).map(|i| i as f64).collect())
            .expect("operation should succeed");

        let mut robust_nystroem = RobustNystroem::new(5)
            .gamma(1.0)
            .with_config(RobustKernelConfig::default());

        robust_nystroem.fit(&x).expect("operation should succeed");
        let transformed = robust_nystroem
            .transform(&x)
            .expect("operation should succeed");

        assert_eq!(transformed.shape(), &[10, 5]);
        assert!(robust_nystroem.get_outlier_mask().is_some());
        assert!(robust_nystroem.get_robust_weights().is_some());
    }

    #[test]
    fn test_breakdown_point_analysis() {
        let x = Array2::from_shape_vec((20, 2), (0..40).map(|i| i as f64).collect())
            .expect("operation should succeed");

        let mut analysis =
            BreakdownPointAnalysis::new().with_estimator(RobustEstimator::Huber { delta: 1.345 });

        let breakdown_point = analysis.analyze(&x).expect("operation should succeed");
        assert!((0.0..=1.0).contains(&breakdown_point));

        let breakdown_points = analysis.get_breakdown_points();
        assert!(!breakdown_points.is_empty());
    }

    #[test]
    fn test_influence_function_diagnostics() {
        let x = Array2::from_shape_vec((10, 2), (0..20).map(|i| i as f64).collect())
            .expect("operation should succeed");
        let y = Array1::from_vec((0..10).map(|i| i as f64).collect());

        let mut diagnostics = InfluenceFunctionDiagnostics::new()
            .with_estimator(RobustEstimator::Huber { delta: 1.345 });

        diagnostics
            .compute(&x, &y)
            .expect("operation should succeed");

        assert!(diagnostics.get_influence_values().is_some());
        assert!(diagnostics.get_leverage_values().is_some());
        assert!(diagnostics.get_cook_distances().is_some());

        let influence = diagnostics
            .get_influence_values()
            .expect("operation should succeed");
        assert_eq!(influence.len(), 10);
    }

    #[test]
    fn test_robust_loss_functions() {
        let losses = vec![
            RobustLoss::Huber { delta: 1.345 },
            RobustLoss::Tukey { c: 4.685 },
            RobustLoss::Cauchy { sigma: 1.0 },
            RobustLoss::Welsch { sigma: 1.0 },
            RobustLoss::Fair { c: 1.0 },
        ];

        let x = Array2::from_shape_vec((8, 2), (0..16).map(|i| i as f64).collect())
            .expect("operation should succeed");

        for loss in losses {
            let config = RobustKernelConfig {
                loss,
                ..Default::default()
            };

            let mut robust_rbf = RobustRBFSampler::new(10).with_config(config);

            let result = robust_rbf.fit(&x);
            assert!(result.is_ok());
        }
    }

    #[test]
    fn test_contamination_resistance() {
        // Create clean data
        let mut x = Array2::from_shape_vec((20, 2), (0..40).map(|i| i as f64 / 10.0).collect())
            .expect("operation should succeed");

        // Add outliers
        x[[18, 0]] = 100.0;
        x[[18, 1]] = 100.0;
        x[[19, 0]] = -100.0;
        x[[19, 1]] = -100.0;

        let config = RobustKernelConfig {
            contamination: 0.2,
            ..Default::default()
        };

        let mut robust_rbf = RobustRBFSampler::new(10).with_config(config);

        robust_rbf.fit(&x).expect("operation should succeed");

        let outlier_mask = robust_rbf
            .get_outlier_mask()
            .expect("operation should succeed");
        assert!(outlier_mask[18] || outlier_mask[19]); // At least one outlier detected
    }
}
