quantrs2-device 0.1.3

//! Statistical analysis components

use super::super::results::*;
use crate::DeviceResult;
use std::collections::HashMap;

/// Statistical analyzer for measurement data
pub struct StatisticalAnalyzer {
    // Configuration and state
}

impl StatisticalAnalyzer {
    /// Create new statistical analyzer
    pub const fn new() -> Self {
        Self {}
    }

    /// Perform statistical analysis
    pub fn analyze(
        &self,
        latencies: &[f64],
        confidences: &[f64],
    ) -> DeviceResult<StatisticalAnalysisResults> {
        let descriptive_stats = self.calculate_descriptive_statistics(latencies, confidences)?;
        let hypothesis_tests = self.perform_hypothesis_tests(latencies, confidences)?;
        let confidence_intervals = self.calculate_confidence_intervals(latencies, confidences)?;
        let effect_sizes = self.calculate_effect_sizes(latencies, confidences)?;

        Ok(StatisticalAnalysisResults {
            descriptive_stats,
            hypothesis_tests,
            confidence_intervals,
            effect_sizes,
        })
    }

    /// Calculate descriptive statistics
    fn calculate_descriptive_statistics(
        &self,
        latencies: &[f64],
        confidences: &[f64],
    ) -> DeviceResult<DescriptiveStatistics> {
        if latencies.is_empty() {
            return Ok(DescriptiveStatistics::default());
        }

        let mean_latency = latencies.iter().sum::<f64>() / latencies.len() as f64;
        let variance = latencies
            .iter()
            .map(|&x| (x - mean_latency).powi(2))
            .sum::<f64>()
            / latencies.len() as f64;
        let std_latency = variance.sqrt();

        // Calculate median
        let mut sorted_latencies = latencies.to_vec();
        sorted_latencies.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        let median_latency = if sorted_latencies.len() % 2 == 0 {
            let mid = sorted_latencies.len() / 2;
            f64::midpoint(sorted_latencies[mid - 1], sorted_latencies[mid])
        } else {
            sorted_latencies[sorted_latencies.len() / 2]
        };

        // Calculate percentiles
        let latency_percentiles = vec![
            self.percentile(&sorted_latencies, 25.0),
            self.percentile(&sorted_latencies, 75.0),
            self.percentile(&sorted_latencies, 95.0),
            self.percentile(&sorted_latencies, 99.0),
        ];

        let success_rate_stats = self.calculate_success_rate_stats(confidences)?;
        let error_rate_distribution = self.calculate_error_rate_distribution(confidences)?;

        Ok(DescriptiveStatistics {
            mean_latency,
            std_latency,
            median_latency,
            latency_percentiles,
            success_rate_stats,
            error_rate_distribution,
        })
    }

    /// Calculate percentile
    fn percentile(&self, sorted_data: &[f64], percentile: f64) -> f64 {
        if sorted_data.is_empty() {
            return 0.0;
        }

        let index = (percentile / 100.0) * (sorted_data.len() - 1) as f64;
        let lower = index.floor() as usize;
        let upper = index.ceil() as usize;

        if lower == upper {
            sorted_data[lower]
        } else {
            (index - lower as f64)
                .mul_add(sorted_data[upper] - sorted_data[lower], sorted_data[lower])
        }
    }

    /// Calculate success rate statistics
    fn calculate_success_rate_stats(
        &self,
        confidences: &[f64],
    ) -> DeviceResult<MeasurementSuccessStats> {
        if confidences.is_empty() {
            return Ok(MeasurementSuccessStats::default());
        }

        let high_confidence_count = confidences.iter().filter(|&&c| c > 0.95).count();
        let overall_success_rate = high_confidence_count as f64 / confidences.len() as f64;

        // Simple confidence interval calculation
        let n = confidences.len() as f64;
        let p = overall_success_rate;
        let se = ((p * (1.0 - p)) / n).sqrt();
        let z = 1.96; // 95% confidence
        let success_rate_ci = (p - z * se, p + z * se);

        Ok(MeasurementSuccessStats {
            overall_success_rate,
            per_qubit_success_rate: HashMap::new(),
            temporal_success_rate: vec![],
            success_rate_ci,
        })
    }

    /// Calculate error rate distribution
    fn calculate_error_rate_distribution(
        &self,
        confidences: &[f64],
    ) -> DeviceResult<ErrorRateDistribution> {
        let error_rates: Vec<f64> = confidences.iter().map(|&c| 1.0 - c).collect();

        // Simple histogram (10 bins)
        let mut histogram = vec![(0.0, 0); 10];
        for &error_rate in &error_rates {
            let bin = ((error_rate * 10.0).floor() as usize).min(9);
            histogram[bin].1 += 1;
            histogram[bin].0 = (bin as f64 + 0.5) / 10.0;
        }

        Ok(ErrorRateDistribution {
            histogram,
            best_fit_distribution: "normal".to_string(),
            distribution_parameters: vec![0.0, 1.0],
            goodness_of_fit: 0.95,
        })
    }

    /// Perform hypothesis tests
    fn perform_hypothesis_tests(
        &self,
        latencies: &[f64],
        confidences: &[f64],
    ) -> DeviceResult<HypothesisTestResults> {
        let mut independence_tests = HashMap::new();
        let mut stationarity_tests = HashMap::new();
        let mut normality_tests = HashMap::new();
        let mut comparison_tests = HashMap::new();

        // Mock tests for now
        independence_tests.insert(
            "latency_independence".to_string(),
            StatisticalTest {
                statistic: 1.23,
                p_value: 0.15,
                critical_value: 1.96,
                is_significant: false,
                effect_size: Some(0.1),
            },
        );

        Ok(HypothesisTestResults {
            independence_tests,
            stationarity_tests,
            normality_tests,
            comparison_tests,
        })
    }

    /// Calculate confidence intervals
    fn calculate_confidence_intervals(
        &self,
        latencies: &[f64],
        confidences: &[f64],
    ) -> DeviceResult<ConfidenceIntervals> {
        let mut mean_intervals = HashMap::new();
        let mut bootstrap_intervals = HashMap::new();
        let mut prediction_intervals = HashMap::new();

        if !latencies.is_empty() {
            let mean = latencies.iter().sum::<f64>() / latencies.len() as f64;
            let se = self.calculate_standard_error(latencies);
            let margin = 1.96 * se; // 95% confidence

            mean_intervals.insert("latency".to_string(), (mean - margin, mean + margin));
        }

        Ok(ConfidenceIntervals {
            confidence_level: 0.95,
            mean_intervals,
            bootstrap_intervals,
            prediction_intervals,
        })
    }

    /// Calculate standard error
    fn calculate_standard_error(&self, data: &[f64]) -> f64 {
        if data.len() < 2 {
            return 0.0;
        }

        let mean = data.iter().sum::<f64>() / data.len() as f64;
        let variance =
            data.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / (data.len() - 1) as f64;

        (variance / data.len() as f64).sqrt()
    }

    /// Calculate effect sizes
    fn calculate_effect_sizes(
        &self,
        latencies: &[f64],
        confidences: &[f64],
    ) -> DeviceResult<EffectSizeAnalysis> {
        let mut cohens_d = HashMap::new();
        let mut correlations = HashMap::new();
        let mut r_squared = HashMap::new();
        let mut practical_significance = HashMap::new();

        // Simple correlation between latency and confidence
        if latencies.len() == confidences.len() && !latencies.is_empty() {
            let correlation = self.calculate_correlation(latencies, confidences);
            correlations.insert("latency_confidence".to_string(), correlation);
            r_squared.insert("latency_confidence".to_string(), correlation.powi(2));
            practical_significance
                .insert("latency_confidence".to_string(), correlation.abs() > 0.3);
        }

        Ok(EffectSizeAnalysis {
            cohens_d,
            correlations,
            r_squared,
            practical_significance,
        })
    }

    /// Calculate Pearson correlation
    fn calculate_correlation(&self, x: &[f64], y: &[f64]) -> f64 {
        if x.len() != y.len() || x.len() < 2 {
            return 0.0;
        }

        let n = x.len() as f64;
        let mean_x = x.iter().sum::<f64>() / n;
        let mean_y = y.iter().sum::<f64>() / n;

        let numerator: f64 = x
            .iter()
            .zip(y.iter())
            .map(|(&xi, &yi)| (xi - mean_x) * (yi - mean_y))
            .sum();

        let sum_sq_x: f64 = x.iter().map(|&xi| (xi - mean_x).powi(2)).sum();
        let sum_sq_y: f64 = y.iter().map(|&yi| (yi - mean_y).powi(2)).sum();

        let denominator = (sum_sq_x * sum_sq_y).sqrt();

        if denominator > 1e-10 {
            numerator / denominator
        } else {
            0.0
        }
    }
}

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