quantrs2-device 0.1.3

//! Analysis utilities for dynamical decoupling

use scirs2_core::ndarray::{Array1, Array2};
use std::collections::HashMap;

use super::{config::NoiseType, performance::DDPerformanceAnalysis, sequences::DDSequence};
use crate::DeviceResult;

// SciRS2 dependencies with fallbacks
#[cfg(feature = "scirs2")]
use scirs2_stats::{mean, std};

#[cfg(not(feature = "scirs2"))]
use super::fallback_scirs2::{inv, mean, std, trace};
use scirs2_core::random::prelude::*;

/// Statistical analysis results for DD sequences
#[derive(Debug, Clone)]
pub struct DDStatisticalAnalysis {
    /// Basic statistical measures
    pub basic_statistics: BasicStatistics,
    /// Advanced statistical analysis
    pub advanced_statistics: AdvancedStatistics,
    /// Machine learning insights
    pub ml_insights: Option<MLInsights>,
    /// Uncertainty quantification
    pub uncertainty_analysis: UncertaintyAnalysis,
}

/// Basic statistical measures
#[derive(Debug, Clone)]
pub struct BasicStatistics {
    /// Sample mean
    pub mean: f64,
    /// Sample standard deviation
    pub std_deviation: f64,
    /// Sample variance
    pub variance: f64,
    /// Skewness
    pub skewness: f64,
    /// Kurtosis
    pub kurtosis: f64,
    /// Median
    pub median: f64,
    /// Interquartile range
    pub iqr: f64,
}

/// Advanced statistical analysis
#[derive(Debug, Clone)]
pub struct AdvancedStatistics {
    /// Multivariate analysis
    pub multivariate_analysis: MultivariateAnalysis,
    /// Time series analysis
    pub time_series_analysis: Option<TimeSeriesAnalysis>,
    /// Bayesian analysis
    pub bayesian_analysis: Option<BayesianAnalysis>,
    /// Non-parametric analysis
    pub non_parametric_analysis: NonParametricAnalysis,
}

/// Multivariate statistical analysis
#[derive(Debug, Clone)]
pub struct MultivariateAnalysis {
    /// Principal component analysis
    pub pca_results: PCAResults,
    /// Factor analysis
    pub factor_analysis: FactorAnalysisResults,
    /// Cluster analysis
    pub cluster_analysis: ClusterAnalysisResults,
    /// Discriminant analysis
    pub discriminant_analysis: DiscriminantAnalysisResults,
}

/// Principal Component Analysis results
#[derive(Debug, Clone)]
pub struct PCAResults {
    /// Principal components
    pub components: Array2<f64>,
    /// Explained variance ratio
    pub explained_variance_ratio: Array1<f64>,
    /// Cumulative explained variance
    pub cumulative_variance: Array1<f64>,
    /// Number of components to retain
    pub n_components_retain: usize,
}

/// Factor analysis results
#[derive(Debug, Clone)]
pub struct FactorAnalysisResults {
    /// Factor loadings
    pub loadings: Array2<f64>,
    /// Communalities
    pub communalities: Array1<f64>,
    /// Specific variances
    pub specific_variances: Array1<f64>,
    /// Factor scores
    pub factor_scores: Array2<f64>,
}

/// Cluster analysis results
#[derive(Debug, Clone)]
pub struct ClusterAnalysisResults {
    /// Cluster labels
    pub cluster_labels: Array1<i32>,
    /// Cluster centers
    pub cluster_centers: Array2<f64>,
    /// Silhouette scores
    pub silhouette_scores: Array1<f64>,
    /// Inertia (within-cluster sum of squares)
    pub inertia: f64,
}

/// Discriminant analysis results
#[derive(Debug, Clone)]
pub struct DiscriminantAnalysisResults {
    /// Linear discriminants
    pub linear_discriminants: Array2<f64>,
    /// Classification accuracy
    pub classification_accuracy: f64,
    /// Cross-validation score
    pub cv_score: f64,
    /// Feature importance
    pub feature_importance: Array1<f64>,
}

/// Time series analysis
#[derive(Debug, Clone)]
pub struct TimeSeriesAnalysis {
    /// Trend analysis
    pub trend_analysis: TrendAnalysis,
    /// Seasonality analysis
    pub seasonality_analysis: SeasonalityAnalysis,
    /// Stationarity test results
    pub stationarity_tests: StationarityTests,
    /// Autocorrelation analysis
    pub autocorrelation: AutocorrelationAnalysis,
}

/// Trend analysis results
#[derive(Debug, Clone)]
pub struct TrendAnalysis {
    /// Linear trend slope
    pub linear_slope: f64,
    /// Trend significance
    pub trend_p_value: f64,
    /// R-squared
    pub r_squared: f64,
    /// Trend direction
    pub trend_direction: TrendDirection,
}

/// Trend direction
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum TrendDirection {
    Increasing,
    Decreasing,
    Stable,
    NonLinear,
}

/// Seasonality analysis
#[derive(Debug, Clone)]
pub struct SeasonalityAnalysis {
    /// Seasonal period
    pub seasonal_period: Option<usize>,
    /// Seasonal strength
    pub seasonal_strength: f64,
    /// Seasonal decomposition
    pub seasonal_decomposition: SeasonalDecomposition,
}

/// Seasonal decomposition
#[derive(Debug, Clone)]
pub struct SeasonalDecomposition {
    /// Trend component
    pub trend: Array1<f64>,
    /// Seasonal component
    pub seasonal: Array1<f64>,
    /// Residual component
    pub residual: Array1<f64>,
}

/// Stationarity test results
#[derive(Debug, Clone)]
pub struct StationarityTests {
    /// Augmented Dickey-Fuller test
    pub adf_test: StationarityTest,
    /// KPSS test
    pub kpss_test: StationarityTest,
    /// Phillips-Perron test
    pub pp_test: StationarityTest,
}

/// Individual stationarity test
#[derive(Debug, Clone)]
pub struct StationarityTest {
    /// Test statistic
    pub statistic: f64,
    /// p-value
    pub p_value: f64,
    /// Critical values
    pub critical_values: HashMap<String, f64>,
    /// Is stationary
    pub is_stationary: bool,
}

/// Autocorrelation analysis
#[derive(Debug, Clone)]
pub struct AutocorrelationAnalysis {
    /// Autocorrelation function
    pub acf: Array1<f64>,
    /// Partial autocorrelation function
    pub pacf: Array1<f64>,
    /// Significant lags
    pub significant_lags: Vec<usize>,
    /// Ljung-Box test result
    pub ljung_box_test: LjungBoxTest,
}

/// Ljung-Box test for autocorrelation
#[derive(Debug, Clone)]
pub struct LjungBoxTest {
    /// Test statistic
    pub statistic: f64,
    /// p-value
    pub p_value: f64,
    /// Degrees of freedom
    pub df: usize,
    /// Has significant autocorrelation
    pub has_autocorrelation: bool,
}

/// Bayesian analysis results
#[derive(Debug, Clone)]
pub struct BayesianAnalysis {
    /// Posterior distributions
    pub posterior_distributions: HashMap<String, PosteriorDistribution>,
    /// Credible intervals
    pub credible_intervals: HashMap<String, (f64, f64)>,
    /// Bayes factors
    pub bayes_factors: HashMap<String, f64>,
    /// Model evidence
    pub model_evidence: f64,
}

/// Posterior distribution
#[derive(Debug, Clone)]
pub struct PosteriorDistribution {
    /// Sample values
    pub samples: Array1<f64>,
    /// Mean
    pub mean: f64,
    /// Standard deviation
    pub std: f64,
    /// Highest density interval
    pub hdi: (f64, f64),
}

/// Non-parametric analysis
#[derive(Debug, Clone)]
pub struct NonParametricAnalysis {
    /// Rank-based statistics
    pub rank_statistics: RankStatistics,
    /// Permutation test results
    pub permutation_tests: PermutationTests,
    /// Bootstrap analysis
    pub bootstrap_analysis: BootstrapAnalysis,
    /// Kernel density estimation
    pub kde_analysis: KDEAnalysis,
}

/// Rank-based statistics
#[derive(Debug, Clone)]
pub struct RankStatistics {
    /// Spearman rank correlation
    pub spearman_correlation: f64,
    /// Kendall's tau
    pub kendall_tau: f64,
    /// Mann-Whitney U test
    pub mann_whitney_u: MannWhitneyTest,
    /// Wilcoxon signed-rank test
    pub wilcoxon_test: WilcoxonTest,
}

/// Mann-Whitney U test
#[derive(Debug, Clone)]
pub struct MannWhitneyTest {
    /// U statistic
    pub u_statistic: f64,
    /// p-value
    pub p_value: f64,
    /// Effect size
    pub effect_size: f64,
}

/// Wilcoxon signed-rank test
#[derive(Debug, Clone)]
pub struct WilcoxonTest {
    /// Test statistic
    pub statistic: f64,
    /// p-value
    pub p_value: f64,
    /// Effect size
    pub effect_size: f64,
}

/// Permutation test results
#[derive(Debug, Clone)]
pub struct PermutationTests {
    /// Permutation p-values
    pub p_values: HashMap<String, f64>,
    /// Effect sizes
    pub effect_sizes: HashMap<String, f64>,
    /// Number of permutations
    pub n_permutations: usize,
}

/// Bootstrap analysis
#[derive(Debug, Clone)]
pub struct BootstrapAnalysis {
    /// Bootstrap confidence intervals
    pub confidence_intervals: HashMap<String, (f64, f64)>,
    /// Bootstrap bias
    pub bias: HashMap<String, f64>,
    /// Bootstrap standard errors
    pub standard_errors: HashMap<String, f64>,
    /// Number of bootstrap samples
    pub n_bootstrap: usize,
}

/// Kernel density estimation analysis
#[derive(Debug, Clone)]
pub struct KDEAnalysis {
    /// Density estimates
    pub density_estimates: Array1<f64>,
    /// Bandwidth
    pub bandwidth: f64,
    /// Grid points
    pub grid_points: Array1<f64>,
    /// Kernel type
    pub kernel_type: String,
}

/// Machine learning insights
#[derive(Debug, Clone)]
pub struct MLInsights {
    /// Feature importance from random forest
    pub feature_importance: Array1<f64>,
    /// Anomaly detection results
    pub anomaly_detection: AnomalyDetectionResults,
    /// Predictive modeling results
    pub predictive_modeling: PredictiveModelingResults,
    /// Dimensionality reduction
    pub dimensionality_reduction: DimensionalityReduction,
}

/// Anomaly detection results
#[derive(Debug, Clone)]
pub struct AnomalyDetectionResults {
    /// Anomaly scores
    pub anomaly_scores: Array1<f64>,
    /// Anomaly threshold
    pub threshold: f64,
    /// Detected anomalies
    pub anomalies: Vec<usize>,
    /// Isolation forest results
    pub isolation_forest: IsolationForestResults,
}

/// Isolation forest results
#[derive(Debug, Clone)]
pub struct IsolationForestResults {
    /// Anomaly scores
    pub scores: Array1<f64>,
    /// Path lengths
    pub path_lengths: Array1<f64>,
    /// Contamination rate
    pub contamination: f64,
}

/// Predictive modeling results
#[derive(Debug, Clone)]
pub struct PredictiveModelingResults {
    /// Model performance metrics
    pub performance_metrics: HashMap<String, f64>,
    /// Cross-validation scores
    pub cv_scores: Array1<f64>,
    /// Learning curves
    pub learning_curves: LearningCurves,
    /// Model interpretability
    pub interpretability: ModelInterpretability,
}

/// Learning curves
#[derive(Debug, Clone)]
pub struct LearningCurves {
    /// Training sizes
    pub training_sizes: Array1<f64>,
    /// Training scores
    pub training_scores: Array1<f64>,
    /// Validation scores
    pub validation_scores: Array1<f64>,
}

/// Model interpretability
#[derive(Debug, Clone)]
pub struct ModelInterpretability {
    /// SHAP values
    pub shap_values: Array2<f64>,
    /// Feature attributions
    pub feature_attributions: Array1<f64>,
    /// Partial dependence plots
    pub partial_dependence: HashMap<String, (Array1<f64>, Array1<f64>)>,
}

/// Dimensionality reduction results
#[derive(Debug, Clone)]
pub struct DimensionalityReduction {
    /// t-SNE results
    pub tsne_results: TSNEResults,
    /// UMAP results
    pub umap_results: UMAPResults,
    /// Manifold learning
    pub manifold_learning: ManifoldLearningResults,
}

/// t-SNE results
#[derive(Debug, Clone)]
pub struct TSNEResults {
    /// Embedded coordinates
    pub embedding: Array2<f64>,
    /// Perplexity
    pub perplexity: f64,
    /// KL divergence
    pub kl_divergence: f64,
}

/// UMAP results
#[derive(Debug, Clone)]
pub struct UMAPResults {
    /// Embedded coordinates
    pub embedding: Array2<f64>,
    /// Number of neighbors
    pub n_neighbors: usize,
    /// Minimum distance
    pub min_dist: f64,
}

/// Manifold learning results
#[derive(Debug, Clone)]
pub struct ManifoldLearningResults {
    /// Intrinsic dimensionality estimate
    pub intrinsic_dimension: usize,
    /// Local linearity scores
    pub local_linearity: Array1<f64>,
    /// Manifold embedding
    pub embedding: Array2<f64>,
}

/// Uncertainty quantification
#[derive(Debug, Clone)]
pub struct UncertaintyAnalysis {
    /// Aleatory uncertainty (inherent randomness)
    pub aleatory_uncertainty: f64,
    /// Epistemic uncertainty (model uncertainty)
    pub epistemic_uncertainty: f64,
    /// Total uncertainty
    pub total_uncertainty: f64,
    /// Uncertainty sources
    pub uncertainty_sources: HashMap<String, f64>,
    /// Sensitivity analysis
    pub sensitivity_analysis: SensitivityAnalysis,
}

/// Sensitivity analysis
#[derive(Debug, Clone)]
pub struct SensitivityAnalysis {
    /// First-order sensitivity indices
    pub first_order_indices: Array1<f64>,
    /// Total sensitivity indices
    pub total_indices: Array1<f64>,
    /// Interaction effects
    pub interaction_effects: Array2<f64>,
    /// Morris screening results
    pub morris_screening: MorrisScreeningResults,
}

/// Morris screening results
#[derive(Debug, Clone)]
pub struct MorrisScreeningResults {
    /// Elementary effects
    pub elementary_effects: Array2<f64>,
    /// Means of elementary effects
    pub mu: Array1<f64>,
    /// Standard deviations of elementary effects
    pub sigma: Array1<f64>,
    /// Modified means
    pub mu_star: Array1<f64>,
}

/// DD statistical analyzer
pub struct DDStatisticalAnalyzer;

impl DDStatisticalAnalyzer {
    /// Perform comprehensive statistical analysis
    pub fn perform_statistical_analysis(
        sequence: &DDSequence,
        performance_analysis: &DDPerformanceAnalysis,
    ) -> DeviceResult<DDStatisticalAnalysis> {
        // Create sample data for analysis
        let sample_data = Self::generate_sample_data(sequence, performance_analysis)?;

        let basic_statistics = Self::calculate_basic_statistics(&sample_data)?;
        let advanced_statistics = Self::perform_advanced_analysis(&sample_data)?;
        let ml_insights = Self::extract_ml_insights(&sample_data)?;
        let uncertainty_analysis = Self::quantify_uncertainty(&sample_data)?;

        Ok(DDStatisticalAnalysis {
            basic_statistics,
            advanced_statistics,
            ml_insights: Some(ml_insights),
            uncertainty_analysis,
        })
    }

    /// Generate sample data for analysis
    fn generate_sample_data(
        _sequence: &DDSequence,
        performance_analysis: &DDPerformanceAnalysis,
    ) -> DeviceResult<Array2<f64>> {
        // Create synthetic dataset based on performance metrics
        let n_samples = 100;
        let n_features = performance_analysis.metrics.len();
        let mut data = Array2::zeros((n_samples, n_features));

        // Fill with simulated data
        for i in 0..n_samples {
            for j in 0..n_features {
                data[[i, j]] = thread_rng().random::<f64>();
            }
        }

        Ok(data)
    }

    /// Calculate basic statistics
    fn calculate_basic_statistics(data: &Array2<f64>) -> DeviceResult<BasicStatistics> {
        let flat_data = data.iter().copied().collect::<Vec<f64>>();
        let n = flat_data.len() as f64;

        let mean = flat_data.iter().sum::<f64>() / n;
        let variance = flat_data.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / (n - 1.0);
        let std_deviation = variance.sqrt();

        // Calculate higher moments
        let skewness = flat_data
            .iter()
            .map(|x| ((x - mean) / std_deviation).powi(3))
            .sum::<f64>()
            / n;

        let kurtosis = flat_data
            .iter()
            .map(|x| ((x - mean) / std_deviation).powi(4))
            .sum::<f64>()
            / n
            - 3.0;

        // Calculate median and IQR
        let mut sorted_data = flat_data;
        sorted_data.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

        let median = if sorted_data.len() % 2 == 0 {
            f64::midpoint(
                sorted_data[sorted_data.len() / 2 - 1],
                sorted_data[sorted_data.len() / 2],
            )
        } else {
            sorted_data[sorted_data.len() / 2]
        };

        let q1 = sorted_data[sorted_data.len() / 4];
        let q3 = sorted_data[3 * sorted_data.len() / 4];
        let iqr = q3 - q1;

        Ok(BasicStatistics {
            mean,
            std_deviation,
            variance,
            skewness,
            kurtosis,
            median,
            iqr,
        })
    }

    /// Perform advanced statistical analysis (simplified)
    fn perform_advanced_analysis(_data: &Array2<f64>) -> DeviceResult<AdvancedStatistics> {
        // Simplified implementations
        let multivariate_analysis = MultivariateAnalysis {
            pca_results: PCAResults {
                components: Array2::eye(2),
                explained_variance_ratio: Array1::from_vec(vec![0.8, 0.2]),
                cumulative_variance: Array1::from_vec(vec![0.8, 1.0]),
                n_components_retain: 2,
            },
            factor_analysis: FactorAnalysisResults {
                loadings: Array2::eye(2),
                communalities: Array1::from_vec(vec![0.8, 0.9]),
                specific_variances: Array1::from_vec(vec![0.2, 0.1]),
                factor_scores: Array2::zeros((10, 2)),
            },
            cluster_analysis: ClusterAnalysisResults {
                cluster_labels: Array1::zeros(10),
                cluster_centers: Array2::zeros((3, 2)),
                silhouette_scores: Array1::from_vec(vec![0.8, 0.7, 0.9]),
                inertia: 10.0,
            },
            discriminant_analysis: DiscriminantAnalysisResults {
                linear_discriminants: Array2::eye(2),
                classification_accuracy: 0.95,
                cv_score: 0.93,
                feature_importance: Array1::from_vec(vec![0.7, 0.3]),
            },
        };

        Ok(AdvancedStatistics {
            multivariate_analysis,
            time_series_analysis: None,
            bayesian_analysis: None,
            non_parametric_analysis: NonParametricAnalysis {
                rank_statistics: RankStatistics {
                    spearman_correlation: 0.8,
                    kendall_tau: 0.7,
                    mann_whitney_u: MannWhitneyTest {
                        u_statistic: 50.0,
                        p_value: 0.05,
                        effect_size: 0.5,
                    },
                    wilcoxon_test: WilcoxonTest {
                        statistic: 25.0,
                        p_value: 0.03,
                        effect_size: 0.6,
                    },
                },
                permutation_tests: PermutationTests {
                    p_values: HashMap::new(),
                    effect_sizes: HashMap::new(),
                    n_permutations: 1000,
                },
                bootstrap_analysis: BootstrapAnalysis {
                    confidence_intervals: HashMap::new(),
                    bias: HashMap::new(),
                    standard_errors: HashMap::new(),
                    n_bootstrap: 1000,
                },
                kde_analysis: KDEAnalysis {
                    density_estimates: Array1::zeros(100),
                    bandwidth: 0.1,
                    grid_points: Array1::zeros(100),
                    kernel_type: "gaussian".to_string(),
                },
            },
        })
    }

    /// Extract ML insights (simplified)
    fn extract_ml_insights(_data: &Array2<f64>) -> DeviceResult<MLInsights> {
        Ok(MLInsights {
            feature_importance: Array1::from_vec(vec![0.5, 0.3, 0.2]),
            anomaly_detection: AnomalyDetectionResults {
                anomaly_scores: Array1::zeros(100),
                threshold: 0.5,
                anomalies: vec![5, 15, 87],
                isolation_forest: IsolationForestResults {
                    scores: Array1::zeros(100),
                    path_lengths: Array1::zeros(100),
                    contamination: 0.1,
                },
            },
            predictive_modeling: PredictiveModelingResults {
                performance_metrics: HashMap::new(),
                cv_scores: Array1::from_vec(vec![0.9, 0.85, 0.92, 0.88, 0.90]),
                learning_curves: LearningCurves {
                    training_sizes: Array1::from_vec(vec![10.0, 25.0, 50.0, 75.0, 100.0]),
                    training_scores: Array1::from_vec(vec![0.8, 0.85, 0.9, 0.92, 0.93]),
                    validation_scores: Array1::from_vec(vec![0.75, 0.82, 0.87, 0.89, 0.90]),
                },
                interpretability: ModelInterpretability {
                    shap_values: Array2::zeros((100, 3)),
                    feature_attributions: Array1::from_vec(vec![0.4, 0.35, 0.25]),
                    partial_dependence: HashMap::new(),
                },
            },
            dimensionality_reduction: DimensionalityReduction {
                tsne_results: TSNEResults {
                    embedding: Array2::zeros((100, 2)),
                    perplexity: 30.0,
                    kl_divergence: 1.5,
                },
                umap_results: UMAPResults {
                    embedding: Array2::zeros((100, 2)),
                    n_neighbors: 15,
                    min_dist: 0.1,
                },
                manifold_learning: ManifoldLearningResults {
                    intrinsic_dimension: 2,
                    local_linearity: Array1::zeros(100),
                    embedding: Array2::zeros((100, 2)),
                },
            },
        })
    }

    /// Quantify uncertainty (simplified)
    fn quantify_uncertainty(_data: &Array2<f64>) -> DeviceResult<UncertaintyAnalysis> {
        let mut uncertainty_sources = HashMap::new();
        uncertainty_sources.insert("measurement_noise".to_string(), 0.3);
        uncertainty_sources.insert("model_uncertainty".to_string(), 0.2);
        uncertainty_sources.insert("parameter_uncertainty".to_string(), 0.1);

        Ok(UncertaintyAnalysis {
            aleatory_uncertainty: 0.3,
            epistemic_uncertainty: 0.2,
            total_uncertainty: 0.5,
            uncertainty_sources,
            sensitivity_analysis: SensitivityAnalysis {
                first_order_indices: Array1::from_vec(vec![0.4, 0.3, 0.2, 0.1]),
                total_indices: Array1::from_vec(vec![0.5, 0.4, 0.25, 0.15]),
                interaction_effects: Array2::zeros((4, 4)),
                morris_screening: MorrisScreeningResults {
                    elementary_effects: Array2::zeros((100, 4)),
                    mu: Array1::from_vec(vec![0.2, 0.15, 0.1, 0.05]),
                    sigma: Array1::from_vec(vec![0.1, 0.08, 0.06, 0.04]),
                    mu_star: Array1::from_vec(vec![0.25, 0.18, 0.12, 0.08]),
                },
            },
        })
    }
}