trustformers_debug/model_diagnostics/

analytics.rs

//! Advanced analytics for model behavior analysis.
//!
//! This module provides sophisticated analytical capabilities including
//! clustering analysis, temporal dynamics monitoring, representation
//! stability assessment, and multi-dimensional data visualization.
// reason: debug/profiling scaffolding — structs are constructed and their fields/methods
// are retained for the data model, serialization completeness, and future consumers that
// do not yet read every member. Consolidated from many item-level #[allow(dead_code)].
#![allow(dead_code)]

use anyhow::Result;
use std::collections::{HashMap, VecDeque};

use super::types::{
    ActivationHeatmap, ClusteringResults, DriftInfo, HiddenStateAnalysis, LayerActivationStats,
    ModelPerformanceMetrics, RepresentationStability, TemporalDynamics,
};

/// Advanced analytics engine for model behavior analysis.
#[derive(Debug)]
pub struct AdvancedAnalytics {
    /// Analytics configuration
    config: AnalyticsConfig,
    /// Historical hidden state data
    hidden_states_history: VecDeque<HiddenStateData>,
    /// Performance correlation data
    performance_correlations: HashMap<String, CorrelationData>,
    /// Temporal analysis cache
    temporal_analysis_cache: TemporalAnalysisCache,
    /// Clustering analysis results
    clustering_results_cache: HashMap<String, ClusteringResults>,
}

/// Configuration for advanced analytics.
#[derive(Debug, Clone)]
pub struct AnalyticsConfig {
    /// Maximum number of historical samples to retain
    pub max_history_samples: usize,
    /// Minimum samples required for clustering analysis
    pub min_clustering_samples: usize,
    /// Number of clusters for k-means analysis
    pub default_num_clusters: usize,
    /// Window size for temporal analysis
    pub temporal_analysis_window: usize,
    /// Drift detection sensitivity
    pub drift_detection_sensitivity: f64,
    /// Correlation analysis threshold
    pub correlation_threshold: f64,
    /// Enable advanced visualizations
    pub enable_visualizations: bool,
}

/// Hidden state data for analysis.
#[derive(Debug, Clone)]
pub struct HiddenStateData {
    /// Layer name
    pub layer_name: String,
    /// Hidden state vectors
    pub hidden_states: Vec<Vec<f64>>,
    /// Corresponding labels or metadata
    pub labels: Option<Vec<String>>,
    /// Timestamp of collection
    pub timestamp: chrono::DateTime<chrono::Utc>,
    /// Training step when collected
    pub training_step: usize,
}

/// Correlation analysis data.
#[derive(Debug, Clone)]
pub struct CorrelationData {
    /// Metric name
    pub metric_name: String,
    /// Historical values
    pub values: VecDeque<f64>,
    /// Correlations with other metrics
    pub correlations: HashMap<String, f64>,
    /// Last update timestamp
    pub last_updated: chrono::DateTime<chrono::Utc>,
}

/// Temporal analysis cache for performance optimization.
#[derive(Debug, Clone)]
pub struct TemporalAnalysisCache {
    /// Cached drift detection results
    pub drift_results: HashMap<String, DriftInfo>,
    /// Cached temporal consistency scores
    pub consistency_scores: HashMap<String, f64>,
    /// Cached stability windows
    pub stability_windows: HashMap<String, Vec<(usize, usize)>>,
    /// Last analysis timestamp
    pub last_analysis: chrono::DateTime<chrono::Utc>,
}

/// Advanced clustering analysis parameters.
#[derive(Debug, Clone)]
pub struct ClusteringParameters {
    /// Number of clusters
    pub num_clusters: usize,
    /// Maximum iterations for k-means
    pub max_iterations: usize,
    /// Convergence tolerance
    pub tolerance: f64,
    /// Random seed for reproducibility
    pub random_seed: Option<u64>,
    /// Distance metric to use
    pub distance_metric: DistanceMetric,
}

/// Distance metrics for clustering.
#[derive(Debug, Clone)]
pub enum DistanceMetric {
    /// Euclidean distance
    Euclidean,
    /// Manhattan distance
    Manhattan,
    /// Cosine similarity
    Cosine,
    /// Minkowski distance with parameter p
    Minkowski { p: f64 },
}

/// Dimensionality reduction parameters.
#[derive(Debug, Clone)]
pub struct DimensionalityReductionParams {
    /// Target dimensions
    pub target_dimensions: usize,
    /// Reduction method
    pub method: ReductionMethod,
    /// Preserve variance ratio
    pub preserve_variance_ratio: f64,
}

/// Dimensionality reduction methods.
#[derive(Debug, Clone)]
pub enum ReductionMethod {
    /// Principal Component Analysis
    PCA,
    /// t-SNE
    TSNE { perplexity: f64 },
    /// UMAP
    UMAP { n_neighbors: usize, min_dist: f64 },
}

/// Visualization generation parameters.
#[derive(Debug, Clone)]
pub struct VisualizationParams {
    /// Output dimensions (width, height)
    pub dimensions: (usize, usize),
    /// Color scheme
    pub color_scheme: ColorScheme,
    /// Include annotations
    pub include_annotations: bool,
    /// Export format
    pub export_format: ExportFormat,
}

/// Color schemes for visualizations.
#[derive(Debug, Clone)]
pub enum ColorScheme {
    /// Viridis color scale
    Viridis,
    /// Plasma color scale
    Plasma,
    /// Inferno color scale
    Inferno,
    /// Custom color map
    Custom(Vec<(f64, f64, f64)>),
}

/// Export formats for visualizations.
#[derive(Debug, Clone)]
pub enum ExportFormat {
    /// PNG image
    PNG,
    /// SVG vector
    SVG,
    /// JSON data
    JSON,
    /// CSV data
    CSV,
}

/// Statistical analysis results.
#[derive(Debug, Clone, serde::Serialize, serde::Deserialize, Default)]
pub struct StatisticalAnalysis {
    /// Mean values
    pub means: Vec<f64>,
    /// Standard deviations
    pub std_devs: Vec<f64>,
    /// Correlation matrix
    pub correlation_matrix: Vec<Vec<f64>>,
    /// Principal components
    pub principal_components: Vec<Vec<f64>>,
    /// Explained variance ratios
    pub explained_variance_ratios: Vec<f64>,
    /// Statistical significance tests
    pub significance_tests: Vec<SignificanceTest>,
}

/// Statistical significance test result.
#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
pub struct SignificanceTest {
    /// Test name
    pub test_name: String,
    /// Test statistic
    pub statistic: f64,
    /// P-value
    pub p_value: f64,
    /// Degrees of freedom
    pub degrees_of_freedom: Option<usize>,
    /// Confidence interval
    pub confidence_interval: Option<(f64, f64)>,
}

/// Anomaly detection results.
#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
pub struct AnomalyDetectionResults {
    /// Detected anomalies
    pub anomalies: Vec<Anomaly>,
    /// Anomaly scores for all data points
    pub anomaly_scores: Vec<f64>,
    /// Detection threshold used
    pub threshold: f64,
    /// Detection method
    pub method: AnomalyDetectionMethod,
}

/// Individual anomaly information.
#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
pub struct Anomaly {
    /// Index in dataset
    pub index: usize,
    /// Anomaly score
    pub score: f64,
    /// Timestamp
    pub timestamp: chrono::DateTime<chrono::Utc>,
    /// Additional context
    pub context: HashMap<String, String>,
}

/// Anomaly detection methods.
#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
pub enum AnomalyDetectionMethod {
    /// Isolation Forest
    IsolationForest { n_trees: usize },
    /// Local Outlier Factor
    LocalOutlierFactor { n_neighbors: usize },
    /// One-Class SVM
    OneClassSVM { nu: f64 },
    /// Statistical threshold
    StatisticalThreshold { n_std: f64 },
}

impl Default for AnalyticsConfig {
    fn default() -> Self {
        Self {
            max_history_samples: 10000,
            min_clustering_samples: 50,
            default_num_clusters: 8,
            temporal_analysis_window: 100,
            drift_detection_sensitivity: 0.05,
            correlation_threshold: 0.7,
            enable_visualizations: true,
        }
    }
}

impl Default for ClusteringParameters {
    fn default() -> Self {
        Self {
            num_clusters: 8,
            max_iterations: 100,
            tolerance: 1e-4,
            random_seed: Some(42),
            distance_metric: DistanceMetric::Euclidean,
        }
    }
}

impl AdvancedAnalytics {
    /// Create a new advanced analytics engine.
    pub fn new() -> Self {
        Self {
            config: AnalyticsConfig::default(),
            hidden_states_history: VecDeque::new(),
            performance_correlations: HashMap::new(),
            temporal_analysis_cache: TemporalAnalysisCache::new(),
            clustering_results_cache: HashMap::new(),
        }
    }

    /// Create analytics engine with custom configuration.
    pub fn with_config(config: AnalyticsConfig) -> Self {
        Self {
            config,
            hidden_states_history: VecDeque::new(),
            performance_correlations: HashMap::new(),
            temporal_analysis_cache: TemporalAnalysisCache::new(),
            clustering_results_cache: HashMap::new(),
        }
    }

    /// Record hidden state data for analysis.
    pub fn record_hidden_states(&mut self, hidden_states: HiddenStateData) {
        self.hidden_states_history.push_back(hidden_states);

        while self.hidden_states_history.len() > self.config.max_history_samples {
            self.hidden_states_history.pop_front();
        }
    }

    /// Record performance metrics for correlation analysis.
    pub fn record_performance_metrics(&mut self, metrics: &ModelPerformanceMetrics) {
        self.update_correlation_data("loss", metrics.loss);
        self.update_correlation_data("throughput", metrics.throughput_samples_per_sec);
        self.update_correlation_data("memory_usage", metrics.memory_usage_mb);

        if let Some(accuracy) = metrics.accuracy {
            self.update_correlation_data("accuracy", accuracy);
        }

        if let Some(gpu_util) = metrics.gpu_utilization {
            self.update_correlation_data("gpu_utilization", gpu_util);
        }
    }

    /// Perform comprehensive hidden state analysis.
    pub fn analyze_hidden_states(&self, layer_name: &str) -> Result<HiddenStateAnalysis> {
        let layer_data: Vec<_> = self
            .hidden_states_history
            .iter()
            .filter(|data| data.layer_name == layer_name)
            .collect();

        if layer_data.is_empty() {
            return Err(anyhow::anyhow!(
                "No hidden state data available for layer: {}",
                layer_name
            ));
        }

        // Extract all hidden states for this layer
        let all_states: Vec<Vec<f64>> =
            layer_data.iter().flat_map(|data| data.hidden_states.iter()).cloned().collect();

        if all_states.is_empty() {
            return Err(anyhow::anyhow!(
                "No hidden states found for layer: {}",
                layer_name
            ));
        }

        let dimensionality = all_states[0].len();

        // Perform clustering analysis
        let clustering_results = self.perform_clustering_analysis(&all_states)?;

        // Analyze temporal dynamics
        let temporal_dynamics = self.analyze_temporal_dynamics(&layer_data)?;

        // Assess representation stability
        let representation_stability = self.assess_representation_stability(&all_states)?;

        // Calculate information content
        let information_content = self.calculate_information_content(&all_states)?;

        Ok(HiddenStateAnalysis {
            dimensionality,
            information_content,
            clustering_results,
            temporal_dynamics,
            representation_stability,
        })
    }

    /// Perform clustering analysis on hidden states.
    pub fn perform_clustering_analysis(&self, data: &[Vec<f64>]) -> Result<ClusteringResults> {
        if data.len() < self.config.min_clustering_samples {
            return Err(anyhow::anyhow!("Insufficient data for clustering analysis"));
        }

        let params = ClusteringParameters::default();
        let num_clusters = params.num_clusters.min(data.len() / 2);

        // Simple k-means clustering implementation
        let mut cluster_centers = self.initialize_cluster_centers(data, num_clusters)?;
        let mut cluster_assignments = vec![0; data.len()];

        for _iteration in 0..params.max_iterations {
            // Assign points to nearest cluster
            let mut new_assignments = vec![0; data.len()];
            for (i, point) in data.iter().enumerate() {
                let mut best_distance = f64::INFINITY;
                let mut best_cluster = 0;

                for (j, center) in cluster_centers.iter().enumerate() {
                    let distance =
                        self.calculate_distance(point, center, &params.distance_metric)?;
                    if distance < best_distance {
                        best_distance = distance;
                        best_cluster = j;
                    }
                }
                new_assignments[i] = best_cluster;
            }

            // Check for convergence
            if new_assignments == cluster_assignments {
                break;
            }
            cluster_assignments = new_assignments;

            // Update cluster centers
            cluster_centers =
                self.update_cluster_centers(data, &cluster_assignments, num_clusters)?;
        }

        // Calculate silhouette score
        let silhouette_score =
            self.calculate_silhouette_score(data, &cluster_assignments, &cluster_centers)?;

        // Calculate inertia
        let inertia = self.calculate_inertia(data, &cluster_assignments, &cluster_centers)?;

        Ok(ClusteringResults {
            num_clusters,
            cluster_centers,
            cluster_assignments,
            silhouette_score,
            inertia,
        })
    }

    /// Analyze temporal dynamics of hidden states.
    pub fn analyze_temporal_dynamics(
        &self,
        layer_data: &[&HiddenStateData],
    ) -> Result<TemporalDynamics> {
        if layer_data.len() < 2 {
            return Err(anyhow::anyhow!("Insufficient temporal data"));
        }

        // Calculate temporal consistency
        let temporal_consistency = self.calculate_temporal_consistency(layer_data)?;

        // Calculate change rate
        let change_rate = self.calculate_change_rate(layer_data)?;

        // Identify stability windows
        let stability_windows = self.identify_stability_windows(layer_data)?;

        // Detect distribution drift
        let drift_detection = self.detect_distribution_drift(layer_data)?;

        Ok(TemporalDynamics {
            temporal_consistency,
            change_rate,
            stability_windows,
            drift_detection,
        })
    }

    /// Assess representation stability.
    pub fn assess_representation_stability(
        &self,
        hidden_states: &[Vec<f64>],
    ) -> Result<RepresentationStability> {
        if hidden_states.is_empty() {
            return Err(anyhow::anyhow!("No hidden states provided"));
        }

        // Calculate overall stability score
        let stability_score = self.calculate_stability_score(hidden_states)?;

        let variance_across_batches = self.calculate_batch_variance(hidden_states)?;

        // Calculate consistency measure
        let consistency_measure = self.calculate_consistency_measure(hidden_states)?;

        let robustness_to_noise = self.assess_noise_robustness(hidden_states)?;

        Ok(RepresentationStability {
            stability_score,
            variance_across_batches,
            consistency_measure,
            robustness_to_noise,
        })
    }

    /// Generate activation heatmap for visualization.
    pub fn generate_activation_heatmap(
        &self,
        layer_stats: &[LayerActivationStats],
    ) -> Result<ActivationHeatmap> {
        if layer_stats.is_empty() {
            return Err(anyhow::anyhow!("No layer statistics provided"));
        }

        // Create heatmap data based on layer statistics
        let mut data = Vec::new();
        let mut min_val = f64::INFINITY;
        let mut max_val = f64::NEG_INFINITY;

        for stats in layer_stats {
            let row = vec![
                stats.mean_activation,
                stats.std_activation,
                stats.min_activation,
                stats.max_activation,
                stats.dead_neurons_ratio,
                stats.saturated_neurons_ratio,
                stats.sparsity,
            ];

            for &val in &row {
                min_val = min_val.min(val);
                max_val = max_val.max(val);
            }

            data.push(row);
        }

        let dimensions = (data.len(), data.first().map_or(0, |row| row.len()));

        Ok(ActivationHeatmap {
            data,
            dimensions,
            value_range: (min_val, max_val),
            interpretation: "Activation statistics heatmap showing layer behavior patterns"
                .to_string(),
        })
    }

    /// Perform anomaly detection on performance metrics.
    pub fn detect_performance_anomalies(&self) -> Result<AnomalyDetectionResults> {
        // Extract performance data
        let mut all_values = Vec::new();
        for correlation_data in self.performance_correlations.values() {
            all_values.extend(correlation_data.values.iter().cloned());
        }

        if all_values.is_empty() {
            return Err(anyhow::anyhow!("No performance data available"));
        }

        // Use statistical threshold method for anomaly detection
        let method = AnomalyDetectionMethod::StatisticalThreshold { n_std: 2.0 };

        let mean = all_values.iter().sum::<f64>() / all_values.len() as f64;
        let variance =
            all_values.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / all_values.len() as f64;
        let std_dev = variance.sqrt();

        let threshold = mean + 2.0 * std_dev;

        let mut anomalies = Vec::new();
        let mut anomaly_scores = Vec::new();

        for (i, &value) in all_values.iter().enumerate() {
            let score = (value - mean).abs() / std_dev;
            anomaly_scores.push(score);

            if value > threshold {
                anomalies.push(Anomaly {
                    index: i,
                    score,
                    timestamp: chrono::Utc::now(),
                    context: HashMap::new(),
                });
            }
        }

        Ok(AnomalyDetectionResults {
            anomalies,
            anomaly_scores,
            threshold,
            method,
        })
    }

    /// Calculate correlation matrix for all metrics.
    pub fn calculate_correlation_matrix(&self) -> Result<Vec<Vec<f64>>> {
        let metric_names: Vec<_> = self.performance_correlations.keys().cloned().collect();
        let n_metrics = metric_names.len();

        if n_metrics == 0 {
            return Err(anyhow::anyhow!(
                "No metrics available for correlation analysis"
            ));
        }

        let mut correlation_matrix = vec![vec![0.0; n_metrics]; n_metrics];

        for (i, metric1) in metric_names.iter().enumerate() {
            for (j, metric2) in metric_names.iter().enumerate() {
                if i == j {
                    correlation_matrix[i][j] = 1.0;
                } else {
                    let correlation = self.calculate_correlation(metric1, metric2)?;
                    correlation_matrix[i][j] = correlation;
                }
            }
        }

        Ok(correlation_matrix)
    }

    /// Perform statistical analysis on collected data.
    pub fn perform_statistical_analysis(&self) -> Result<StatisticalAnalysis> {
        if self.performance_correlations.is_empty() {
            return Err(anyhow::anyhow!(
                "No data available for statistical analysis"
            ));
        }

        // Calculate means and standard deviations
        let mut means = Vec::new();
        let mut std_devs = Vec::new();

        for correlation_data in self.performance_correlations.values() {
            let values: Vec<f64> = correlation_data.values.iter().cloned().collect();
            if !values.is_empty() {
                let mean = values.iter().sum::<f64>() / values.len() as f64;
                let variance =
                    values.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / values.len() as f64;
                let std_dev = variance.sqrt();

                means.push(mean);
                std_devs.push(std_dev);
            }
        }

        // Calculate correlation matrix
        let correlation_matrix = self.calculate_correlation_matrix()?;

        // REAL principal components: the eigen-decomposition of the
        // correlation matrix (PCA on standardised variables), largest
        // eigenvalue first. `principal_components[k]` is the k-th loading
        // vector and `explained_variance_ratios[k]` its share of the total
        // variance.
        //
        // These used to be an all-ones matrix and a flat `1/n` vector, i.e.
        // "every variable loads equally on every component and each component
        // explains the same share", published as a PCA result.
        let (principal_components, explained_variance_ratios) =
            Self::principal_components_of(&correlation_matrix);

        // Hypothesis tests need a stated null and paired samples to test it
        // against; `perform_statistical_analysis` receives neither. This used
        // to carry one "Sample t-test" with statistic 1.0, p-value 0.05 and a
        // (0.0, 1.0) confidence interval -- constants regardless of the data.
        let significance_tests = Vec::new();

        Ok(StatisticalAnalysis {
            means,
            std_devs,
            correlation_matrix,
            principal_components,
            explained_variance_ratios,
            significance_tests,
        })
    }

    /// Eigen-decompose a symmetric correlation matrix into principal
    /// components and explained-variance ratios, largest first.
    ///
    /// Returns empty vectors for an empty or non-square input.
    fn principal_components_of(correlation_matrix: &[Vec<f64>]) -> (Vec<Vec<f64>>, Vec<f64>) {
        let n = correlation_matrix.len();
        if n == 0 || correlation_matrix.iter().any(|row| row.len() != n) {
            return (Vec::new(), Vec::new());
        }
        if correlation_matrix.iter().flatten().any(|v| !v.is_finite()) {
            return (Vec::new(), Vec::new());
        }

        let matrix = nalgebra::DMatrix::<f64>::from_fn(n, n, |i, j| correlation_matrix[i][j]);
        // The correlation matrix is symmetric by construction, so the
        // symmetric (Jacobi) eigensolver applies and returns real eigenvalues.
        let eigen = nalgebra::linalg::SymmetricEigen::new(matrix);

        let mut order: Vec<usize> = (0..n).collect();
        order.sort_by(|&a, &b| {
            eigen.eigenvalues[b]
                .partial_cmp(&eigen.eigenvalues[a])
                .unwrap_or(std::cmp::Ordering::Equal)
        });

        // Negative eigenvalues can only appear as round-off on a
        // positive-semidefinite correlation matrix; clamp them at zero.
        let total: f64 = eigen.eigenvalues.iter().map(|v| v.max(0.0)).sum();
        let components: Vec<Vec<f64>> = order
            .iter()
            .map(|&k| eigen.eigenvectors.column(k).iter().copied().collect())
            .collect();
        let ratios: Vec<f64> = order
            .iter()
            .map(|&k| if total > 0.0 { eigen.eigenvalues[k].max(0.0) / total } else { 0.0 })
            .collect();

        (components, ratios)
    }

    /// Generate comprehensive analytics report.
    pub fn generate_analytics_report(&self) -> Result<AnalyticsReport> {
        let correlation_matrix = self.calculate_correlation_matrix().unwrap_or_default();
        let statistical_analysis = self.perform_statistical_analysis().unwrap_or_default();
        let anomaly_detection = self.detect_performance_anomalies().unwrap_or_default();

        // Analyze each layer if data is available
        let mut layer_analyses = HashMap::new();
        let unique_layers: std::collections::HashSet<String> =
            self.hidden_states_history.iter().map(|data| data.layer_name.clone()).collect();

        for layer_name in unique_layers {
            if let Ok(analysis) = self.analyze_hidden_states(&layer_name) {
                layer_analyses.insert(layer_name, analysis);
            }
        }

        Ok(AnalyticsReport {
            correlation_matrix,
            statistical_analysis,
            layer_analyses,
            anomaly_detection,
            temporal_summary: self.generate_temporal_summary(),
            recommendations: self.generate_analytics_recommendations(),
        })
    }

    // Helper methods

    /// Update correlation data for a metric.
    fn update_correlation_data(&mut self, metric_name: &str, value: f64) {
        let correlation_data = self
            .performance_correlations
            .entry(metric_name.to_string())
            .or_insert_with(|| CorrelationData {
                metric_name: metric_name.to_string(),
                values: VecDeque::new(),
                correlations: HashMap::new(),
                last_updated: chrono::Utc::now(),
            });

        correlation_data.values.push_back(value);
        correlation_data.last_updated = chrono::Utc::now();

        // Limit history size
        while correlation_data.values.len() > self.config.max_history_samples {
            correlation_data.values.pop_front();
        }
    }

    /// Initialize cluster centers using k-means++ algorithm.
    fn initialize_cluster_centers(
        &self,
        data: &[Vec<f64>],
        num_clusters: usize,
    ) -> Result<Vec<Vec<f64>>> {
        if data.is_empty() || num_clusters == 0 {
            return Err(anyhow::anyhow!("Invalid input for cluster initialization"));
        }

        let mut centers = Vec::new();
        let _dimensions = data[0].len();

        // Deterministic seeding: take the first point rather than a random
        // one, so repeated runs over the same data cluster identically. (This
        // is the only departure from textbook k-means++, which samples the
        // first centre uniformly at random.)
        centers.push(data[0].clone());

        // Remaining centres by the k-means++ D^2 rule.
        for _ in 1..num_clusters {
            if centers.len() >= data.len() {
                break;
            }

            let mut best_distance = 0.0;
            let mut best_point = data[0].clone();

            for point in data {
                let mut min_distance = f64::INFINITY;
                for center in &centers {
                    let distance =
                        self.calculate_distance(point, center, &DistanceMetric::Euclidean)?;
                    min_distance = min_distance.min(distance);
                }

                if min_distance > best_distance {
                    best_distance = min_distance;
                    best_point = point.clone();
                }
            }

            centers.push(best_point);
        }

        Ok(centers)
    }

    /// Calculate distance between two points.
    fn calculate_distance(
        &self,
        point1: &[f64],
        point2: &[f64],
        metric: &DistanceMetric,
    ) -> Result<f64> {
        if point1.len() != point2.len() {
            return Err(anyhow::anyhow!("Points must have same dimensionality"));
        }

        match metric {
            DistanceMetric::Euclidean => {
                let sum_squared =
                    point1.iter().zip(point2.iter()).map(|(a, b)| (a - b).powi(2)).sum::<f64>();
                Ok(sum_squared.sqrt())
            },
            DistanceMetric::Manhattan => {
                let sum_abs =
                    point1.iter().zip(point2.iter()).map(|(a, b)| (a - b).abs()).sum::<f64>();
                Ok(sum_abs)
            },
            DistanceMetric::Cosine => {
                let dot_product = point1.iter().zip(point2.iter()).map(|(a, b)| a * b).sum::<f64>();
                let norm1 = point1.iter().map(|x| x.powi(2)).sum::<f64>().sqrt();
                let norm2 = point2.iter().map(|x| x.powi(2)).sum::<f64>().sqrt();

                if norm1 == 0.0 || norm2 == 0.0 {
                    Ok(1.0)
                } else {
                    Ok(1.0 - (dot_product / (norm1 * norm2)))
                }
            },
            DistanceMetric::Minkowski { p } => {
                let sum_powered = point1
                    .iter()
                    .zip(point2.iter())
                    .map(|(a, b)| (a - b).abs().powf(*p))
                    .sum::<f64>();
                Ok(sum_powered.powf(1.0 / p))
            },
        }
    }

    /// Update cluster centers based on current assignments.
    fn update_cluster_centers(
        &self,
        data: &[Vec<f64>],
        assignments: &[usize],
        num_clusters: usize,
    ) -> Result<Vec<Vec<f64>>> {
        let dimensions = data[0].len();
        let mut new_centers = vec![vec![0.0; dimensions]; num_clusters];
        let mut cluster_counts = vec![0; num_clusters];

        // Sum points in each cluster
        for (point, &cluster_id) in data.iter().zip(assignments.iter()) {
            if cluster_id < num_clusters {
                for (i, &value) in point.iter().enumerate() {
                    new_centers[cluster_id][i] += value;
                }
                cluster_counts[cluster_id] += 1;
            }
        }

        // Average to get new centers
        for (cluster_id, count) in cluster_counts.iter().enumerate() {
            if *count > 0 {
                for value in &mut new_centers[cluster_id] {
                    *value /= *count as f64;
                }
            }
        }

        Ok(new_centers)
    }

    /// Calculate silhouette score for clustering quality.
    fn calculate_silhouette_score(
        &self,
        data: &[Vec<f64>],
        assignments: &[usize],
        centers: &[Vec<f64>],
    ) -> Result<f64> {
        if data.is_empty() {
            return Ok(0.0);
        }

        let mut total_score = 0.0;
        let mut valid_points = 0;

        for (i, point) in data.iter().enumerate() {
            let cluster_id = assignments[i];

            // Calculate average distance to points in same cluster (a)
            let mut same_cluster_distances = Vec::new();
            for (j, other_point) in data.iter().enumerate() {
                if i != j && assignments[j] == cluster_id {
                    let distance =
                        self.calculate_distance(point, other_point, &DistanceMetric::Euclidean)?;
                    same_cluster_distances.push(distance);
                }
            }

            let a = if same_cluster_distances.is_empty() {
                0.0
            } else {
                same_cluster_distances.iter().sum::<f64>() / same_cluster_distances.len() as f64
            };

            // Calculate minimum average distance to points in other clusters (b)
            let mut min_other_cluster_distance = f64::INFINITY;
            for (other_cluster_id, _) in centers.iter().enumerate() {
                if other_cluster_id != cluster_id {
                    let mut other_cluster_distances = Vec::new();
                    for (j, other_point) in data.iter().enumerate() {
                        if assignments[j] == other_cluster_id {
                            let distance = self.calculate_distance(
                                point,
                                other_point,
                                &DistanceMetric::Euclidean,
                            )?;
                            other_cluster_distances.push(distance);
                        }
                    }

                    if !other_cluster_distances.is_empty() {
                        let avg_distance = other_cluster_distances.iter().sum::<f64>()
                            / other_cluster_distances.len() as f64;
                        min_other_cluster_distance = min_other_cluster_distance.min(avg_distance);
                    }
                }
            }

            let b = min_other_cluster_distance;

            if a < b {
                total_score += (b - a) / b;
            } else if a > b {
                total_score += (b - a) / a;
            }
            // If a == b, silhouette score is 0 (no contribution)

            valid_points += 1;
        }

        Ok(if valid_points > 0 { total_score / valid_points as f64 } else { 0.0 })
    }

    /// Calculate inertia (sum of squared distances to centroids).
    fn calculate_inertia(
        &self,
        data: &[Vec<f64>],
        assignments: &[usize],
        centers: &[Vec<f64>],
    ) -> Result<f64> {
        let mut inertia = 0.0;

        for (point, &cluster_id) in data.iter().zip(assignments.iter()) {
            if cluster_id < centers.len() {
                let distance = self.calculate_distance(
                    point,
                    &centers[cluster_id],
                    &DistanceMetric::Euclidean,
                )?;
                inertia += distance.powi(2);
            }
        }

        Ok(inertia)
    }

    /// Calculate information content of hidden states.
    fn calculate_information_content(&self, hidden_states: &[Vec<f64>]) -> Result<f64> {
        if hidden_states.is_empty() {
            return Ok(0.0);
        }

        let dimensions = hidden_states[0].len();
        let mut total_variance = 0.0;

        for dim in 0..dimensions {
            let values: Vec<f64> = hidden_states.iter().map(|state| state[dim]).collect();
            if values.len() > 1 {
                let mean = values.iter().sum::<f64>() / values.len() as f64;
                let variance = values.iter().map(|x| (x - mean).powi(2)).sum::<f64>()
                    / (values.len() - 1) as f64;
                total_variance += variance;
            }
        }

        // Information content as normalized variance
        Ok(total_variance / dimensions as f64)
    }

    /// Calculate temporal consistency of hidden states.
    fn calculate_temporal_consistency(&self, layer_data: &[&HiddenStateData]) -> Result<f64> {
        if layer_data.len() < 2 {
            return Ok(1.0);
        }

        let mut consistency_scores = Vec::new();

        for i in 1..layer_data.len() {
            let prev_states = &layer_data[i - 1].hidden_states;
            let curr_states = &layer_data[i].hidden_states;

            if !prev_states.is_empty() && !curr_states.is_empty() {
                // Simple consistency measure based on mean state similarity
                let prev_mean = self.calculate_mean_state(prev_states);
                let curr_mean = self.calculate_mean_state(curr_states);

                if prev_mean.len() == curr_mean.len() {
                    let distance = self.calculate_distance(
                        &prev_mean,
                        &curr_mean,
                        &DistanceMetric::Euclidean,
                    )?;
                    consistency_scores.push(1.0 / (1.0 + distance));
                }
            }
        }

        Ok(if consistency_scores.is_empty() {
            1.0
        } else {
            consistency_scores.iter().sum::<f64>() / consistency_scores.len() as f64
        })
    }

    /// Calculate mean state from a collection of states.
    fn calculate_mean_state(&self, states: &[Vec<f64>]) -> Vec<f64> {
        if states.is_empty() {
            return Vec::new();
        }

        let dimensions = states[0].len();
        let mut mean_state = vec![0.0; dimensions];

        for state in states {
            for (i, &value) in state.iter().enumerate() {
                if i < dimensions {
                    mean_state[i] += value;
                }
            }
        }

        for value in &mut mean_state {
            *value /= states.len() as f64;
        }

        mean_state
    }

    /// Calculate change rate between temporal samples.
    fn calculate_change_rate(&self, layer_data: &[&HiddenStateData]) -> Result<f64> {
        if layer_data.len() < 2 {
            return Ok(0.0);
        }

        let mut total_change = 0.0;
        let mut valid_comparisons = 0;

        for i in 1..layer_data.len() {
            let prev_mean = self.calculate_mean_state(&layer_data[i - 1].hidden_states);
            let curr_mean = self.calculate_mean_state(&layer_data[i].hidden_states);

            if !prev_mean.is_empty() && !curr_mean.is_empty() && prev_mean.len() == curr_mean.len()
            {
                let change =
                    self.calculate_distance(&prev_mean, &curr_mean, &DistanceMetric::Euclidean)?;
                total_change += change;
                valid_comparisons += 1;
            }
        }

        Ok(if valid_comparisons > 0 {
            total_change / valid_comparisons as f64
        } else {
            0.0
        })
    }

    /// Identify stability windows in temporal data.
    fn identify_stability_windows(
        &self,
        layer_data: &[&HiddenStateData],
    ) -> Result<Vec<(usize, usize)>> {
        if layer_data.len() < 3 {
            return Ok(Vec::new());
        }

        let mut stability_windows = Vec::new();
        let mut window_start = 0;
        let stability_threshold = 0.1; // Configurable threshold

        for i in 1..layer_data.len() {
            let prev_mean = self.calculate_mean_state(&layer_data[i - 1].hidden_states);
            let curr_mean = self.calculate_mean_state(&layer_data[i].hidden_states);

            if !prev_mean.is_empty() && !curr_mean.is_empty() && prev_mean.len() == curr_mean.len()
            {
                let change = self
                    .calculate_distance(&prev_mean, &curr_mean, &DistanceMetric::Euclidean)
                    .unwrap_or(f64::INFINITY);

                if change > stability_threshold {
                    // End of stability window
                    if i - window_start > 2 {
                        stability_windows.push((window_start, i - 1));
                    }
                    window_start = i;
                }
            }
        }

        // Handle final window
        if layer_data.len() - window_start > 2 {
            stability_windows.push((window_start, layer_data.len() - 1));
        }

        Ok(stability_windows)
    }

    /// Detect distribution drift in temporal data.
    fn detect_distribution_drift(&self, layer_data: &[&HiddenStateData]) -> Result<DriftInfo> {
        if layer_data.len() < self.config.temporal_analysis_window {
            return Ok(DriftInfo {
                drift_detected: false,
                drift_magnitude: 0.0,
                drift_direction: "unknown".to_string(),
                onset_step: None,
            });
        }

        let window_size = self.config.temporal_analysis_window;
        let mid_point = layer_data.len() / 2;

        // Compare early and late windows
        let early_data = &layer_data[0..window_size.min(mid_point)];
        let late_data = &layer_data[mid_point.max(layer_data.len() - window_size)..];

        let early_mean = self.calculate_aggregated_mean(early_data);
        let late_mean = self.calculate_aggregated_mean(late_data);

        if early_mean.len() == late_mean.len() && !early_mean.is_empty() {
            let drift_magnitude =
                self.calculate_distance(&early_mean, &late_mean, &DistanceMetric::Euclidean)?;
            let drift_detected = drift_magnitude > self.config.drift_detection_sensitivity;

            Ok(DriftInfo {
                drift_detected,
                drift_magnitude,
                drift_direction: if drift_detected {
                    "forward".to_string()
                } else {
                    "stable".to_string()
                },
                onset_step: if drift_detected { Some(mid_point) } else { None },
            })
        } else {
            Ok(DriftInfo {
                drift_detected: false,
                drift_magnitude: 0.0,
                drift_direction: "unknown".to_string(),
                onset_step: None,
            })
        }
    }

    /// Calculate aggregated mean across multiple data samples.
    fn calculate_aggregated_mean(&self, layer_data: &[&HiddenStateData]) -> Vec<f64> {
        let all_states: Vec<Vec<f64>> =
            layer_data.iter().flat_map(|data| data.hidden_states.iter()).cloned().collect();

        self.calculate_mean_state(&all_states)
    }

    /// Calculate stability score for representation.
    fn calculate_stability_score(&self, hidden_states: &[Vec<f64>]) -> Result<f64> {
        if hidden_states.len() < 2 {
            return Ok(1.0);
        }

        let mut stability_scores = Vec::new();
        let window_size = (hidden_states.len() / 10).max(2);

        for i in window_size..hidden_states.len() {
            let current_window = &hidden_states[i - window_size..i];
            let mean_current = self.calculate_mean_state(current_window);

            if i >= 2 * window_size {
                let prev_window = &hidden_states[i - 2 * window_size..i - window_size];
                let mean_prev = self.calculate_mean_state(prev_window);

                if mean_current.len() == mean_prev.len() && !mean_current.is_empty() {
                    let distance = self.calculate_distance(
                        &mean_current,
                        &mean_prev,
                        &DistanceMetric::Euclidean,
                    )?;
                    stability_scores.push(1.0 / (1.0 + distance));
                }
            }
        }

        Ok(if stability_scores.is_empty() {
            1.0
        } else {
            stability_scores.iter().sum::<f64>() / stability_scores.len() as f64
        })
    }

    /// Mean per-dimension sample variance across the supplied hidden states.
    fn calculate_batch_variance(&self, hidden_states: &[Vec<f64>]) -> Result<f64> {
        if hidden_states.is_empty() {
            return Ok(0.0);
        }

        let dimensions = hidden_states[0].len();
        let mut total_variance = 0.0;

        for dim in 0..dimensions {
            let values: Vec<f64> = hidden_states.iter().map(|state| state[dim]).collect();
            if values.len() > 1 {
                let mean = values.iter().sum::<f64>() / values.len() as f64;
                let variance = values.iter().map(|x| (x - mean).powi(2)).sum::<f64>()
                    / (values.len() - 1) as f64;
                total_variance += variance;
            }
        }

        Ok(total_variance / dimensions as f64)
    }

    /// Calculate consistency measure for representation.
    fn calculate_consistency_measure(&self, hidden_states: &[Vec<f64>]) -> Result<f64> {
        if hidden_states.len() < 2 {
            return Ok(1.0);
        }

        // Calculate pairwise similarities and return average
        let mut similarities = Vec::new();
        let sample_size = hidden_states.len().min(100); // Limit for performance

        for i in 0..sample_size {
            for j in (i + 1)..sample_size {
                let distance = self.calculate_distance(
                    &hidden_states[i],
                    &hidden_states[j],
                    &DistanceMetric::Cosine,
                )?;
                similarities.push(1.0 - distance); // Convert distance to similarity
            }
        }

        Ok(if similarities.is_empty() {
            1.0
        } else {
            similarities.iter().sum::<f64>() / similarities.len() as f64
        })
    }

    /// Variance-derived noise-robustness PROXY, `1 / (1 + variance)`, in
    /// `(0, 1]`.
    ///
    /// It is a monotone transform of the real
    /// [`Self::calculate_batch_variance`], not a measurement of robustness:
    /// nothing is perturbed and the model is never re-evaluated. Tightly
    /// clustered representations score near 1 and widely spread ones near 0,
    /// which is a heuristic for -- not evidence of -- noise robustness.
    fn assess_noise_robustness(&self, hidden_states: &[Vec<f64>]) -> Result<f64> {
        self.calculate_batch_variance(hidden_states).map(|variance| {
            // High variance might indicate low robustness to noise
            1.0 / (1.0 + variance)
        })
    }

    /// Calculate correlation between two metrics.
    fn calculate_correlation(&self, metric1: &str, metric2: &str) -> Result<f64> {
        let data1 = self
            .performance_correlations
            .get(metric1)
            .ok_or_else(|| anyhow::anyhow!("Metric {} not found", metric1))?;

        let data2 = self
            .performance_correlations
            .get(metric2)
            .ok_or_else(|| anyhow::anyhow!("Metric {} not found", metric2))?;

        let values1: Vec<f64> = data1.values.iter().cloned().collect();
        let values2: Vec<f64> = data2.values.iter().cloned().collect();

        if values1.len() != values2.len() || values1.is_empty() {
            return Ok(0.0);
        }

        let mean1 = values1.iter().sum::<f64>() / values1.len() as f64;
        let mean2 = values2.iter().sum::<f64>() / values2.len() as f64;

        let numerator: f64 = values1
            .iter()
            .zip(values2.iter())
            .map(|(x1, x2)| (x1 - mean1) * (x2 - mean2))
            .sum();

        let var1: f64 = values1.iter().map(|x| (x - mean1).powi(2)).sum();
        let var2: f64 = values2.iter().map(|x| (x - mean2).powi(2)).sum();

        let denominator = (var1 * var2).sqrt();

        Ok(if denominator == 0.0 { 0.0 } else { numerator / denominator })
    }

    /// Generate temporal summary.
    fn generate_temporal_summary(&self) -> String {
        format!(
            "Temporal analysis: {} hidden state samples collected across {} layers. \
            Average stability observed with {} correlation metrics tracked.",
            self.hidden_states_history.len(),
            self.hidden_states_history
                .iter()
                .map(|data| &data.layer_name)
                .collect::<std::collections::HashSet<_>>()
                .len(),
            self.performance_correlations.len()
        )
    }

    /// Generate analytics recommendations.
    fn generate_analytics_recommendations(&self) -> Vec<String> {
        let mut recommendations = Vec::new();

        if self.performance_correlations.len() < 3 {
            recommendations.push(
                "Collect more performance metrics for comprehensive correlation analysis"
                    .to_string(),
            );
        }

        if self.hidden_states_history.len() < 50 {
            recommendations
                .push("Increase hidden state sampling for better temporal analysis".to_string());
        }

        recommendations
            .push("Consider implementing automated anomaly detection alerts".to_string());
        recommendations.push("Enable advanced visualization for better insights".to_string());

        recommendations
    }
}

/// Comprehensive analytics report.
#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
pub struct AnalyticsReport {
    /// Correlation matrix between metrics
    pub correlation_matrix: Vec<Vec<f64>>,
    /// Statistical analysis results
    pub statistical_analysis: StatisticalAnalysis,
    /// Layer-specific analyses
    pub layer_analyses: HashMap<String, HiddenStateAnalysis>,
    /// Anomaly detection results
    pub anomaly_detection: AnomalyDetectionResults,
    /// Temporal analysis summary
    pub temporal_summary: String,
    /// Analytics recommendations
    pub recommendations: Vec<String>,
}

impl TemporalAnalysisCache {
    /// Create a new temporal analysis cache.
    fn new() -> Self {
        Self {
            drift_results: HashMap::new(),
            consistency_scores: HashMap::new(),
            stability_windows: HashMap::new(),
            last_analysis: chrono::Utc::now(),
        }
    }
}

impl Default for AnomalyDetectionResults {
    fn default() -> Self {
        Self {
            anomalies: Vec::new(),
            anomaly_scores: Vec::new(),
            threshold: 0.0,
            method: AnomalyDetectionMethod::StatisticalThreshold { n_std: 2.0 },
        }
    }
}

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

#[cfg(test)]
mod tests {
    use super::*;

    // ---- Wave 6c debug-sweep2: real PCA ----------------------------------

    #[test]
    fn principal_components_of_the_identity_are_unit_and_equally_weighted() {
        let identity = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
        let (components, ratios) = AdvancedAnalytics::principal_components_of(&identity);
        assert_eq!(components.len(), 2);
        // Uncorrelated unit-variance variables split the variance evenly.
        for ratio in &ratios {
            assert!((ratio - 0.5).abs() < 1e-9, "expected 0.5, got {ratio}");
        }
    }

    #[test]
    fn principal_components_of_perfectly_correlated_variables_collapse_to_one() {
        // Two perfectly correlated variables: PC1 explains everything.
        let correlated = vec![vec![1.0, 1.0], vec![1.0, 1.0]];
        let (components, ratios) = AdvancedAnalytics::principal_components_of(&correlated);
        assert_eq!(components.len(), 2);
        assert!(
            (ratios[0] - 1.0).abs() < 1e-9,
            "PC1 must explain all the variance, got {}",
            ratios[0]
        );
        assert!(
            ratios[1].abs() < 1e-9,
            "PC2 must explain none, got {}",
            ratios[1]
        );
        // The old placeholder returned a flat 1/n for every component, so this
        // pair would both have been 0.5.
        assert!(
            (ratios[0] - ratios[1]).abs() > 0.5,
            "the ratios must actually differ"
        );
    }

    #[test]
    fn principal_components_are_ordered_by_explained_variance() {
        let matrix = vec![vec![1.0, 0.8], vec![0.8, 1.0]];
        let (_, ratios) = AdvancedAnalytics::principal_components_of(&matrix);
        assert!(
            ratios[0] >= ratios[1],
            "components must be sorted descending: {ratios:?}"
        );
        assert!(
            (ratios.iter().sum::<f64>() - 1.0).abs() < 1e-9,
            "ratios must sum to 1"
        );
        // Eigenvalues of [[1,0.8],[0.8,1]] are 1.8 and 0.2 => 0.9 / 0.1.
        assert!((ratios[0] - 0.9).abs() < 1e-9, "{ratios:?}");
    }

    #[test]
    fn principal_components_reject_a_malformed_matrix() {
        assert_eq!(
            AdvancedAnalytics::principal_components_of(&[]),
            (Vec::new(), Vec::new())
        );
        let ragged = vec![vec![1.0, 0.0], vec![0.0]];
        assert_eq!(
            AdvancedAnalytics::principal_components_of(&ragged),
            (Vec::new(), Vec::new())
        );
    }

    #[test]
    fn test_advanced_analytics_creation() {
        let analytics = AdvancedAnalytics::new();
        assert_eq!(analytics.hidden_states_history.len(), 0);
        assert_eq!(analytics.performance_correlations.len(), 0);
    }

    #[test]
    fn test_distance_calculation() {
        let analytics = AdvancedAnalytics::new();
        let point1 = vec![1.0, 2.0, 3.0];
        let point2 = vec![4.0, 5.0, 6.0];

        let distance = analytics
            .calculate_distance(&point1, &point2, &DistanceMetric::Euclidean)
            .expect("operation failed in test");
        assert!(distance > 0.0);
    }

    #[test]
    fn test_clustering_parameters() {
        let params = ClusteringParameters::default();
        assert_eq!(params.num_clusters, 8);
        assert_eq!(params.max_iterations, 100);
    }

    #[test]
    fn test_correlation_calculation() {
        let mut analytics = AdvancedAnalytics::new();

        // Add some test data
        analytics.update_correlation_data("metric1", 1.0);
        analytics.update_correlation_data("metric1", 2.0);
        analytics.update_correlation_data("metric2", 3.0);
        analytics.update_correlation_data("metric2", 4.0);

        assert_eq!(analytics.performance_correlations.len(), 2);
    }
}