quantrs2-anneal 0.1.3

//! Main solution clustering analyzer implementation

use scirs2_core::random::prelude::*;
use scirs2_core::random::ChaCha8Rng;
use scirs2_core::random::{Rng, SeedableRng};
use std::collections::{HashMap, VecDeque};
use std::time::{Duration, Instant};

use super::algorithms::{ClusteringAlgorithm, DistanceMetric, LinkageType};
use super::config::{ClusteringConfig, FeatureExtractionMethod};
use super::error::{ClusteringError, ClusteringResult};
use super::types::{
    AnalysisStatistics, ClusterQualityMetrics, ClusterStatistics, ClusteringPerformanceMetrics,
    ClusteringResults, ConnectivityAnalysis, ConvergenceAnalysis, CorrelationAnalysis,
    DifficultyLevel, DistributionAnalysis, DistributionType, EfficiencyMetrics, EnergyBasin,
    EnergyStatistics, FunnelAnalysis, LandscapeAnalysis, MultiModalityAnalysis,
    OptimizationRecommendation, OverallClusteringQuality, PlateauAnalysis, PriorityLevel,
    RecommendationType, RuggednessMetrics, ScalabilityMetrics, SolutionCluster, SolutionMetadata,
    SolutionPoint, StatisticalSummary,
};
use crate::simulator::AnnealingSolution;

/// Solution clustering analyzer
pub struct SolutionClusteringAnalyzer {
    /// Configuration
    config: ClusteringConfig,
    /// Cached distance matrices
    distance_cache: HashMap<String, Vec<Vec<f64>>>,
    /// Analysis statistics
    stats: AnalysisStatistics,
}

impl SolutionClusteringAnalyzer {
    /// Create a new solution clustering analyzer
    #[must_use]
    pub fn new(config: ClusteringConfig) -> Self {
        Self {
            config,
            distance_cache: HashMap::new(),
            stats: AnalysisStatistics {
                total_solutions: 0,
                total_time: Duration::from_secs(0),
                cache_hit_rate: 0.0,
                peak_memory: 0,
            },
        }
    }

    /// Analyze a collection of solutions
    pub fn analyze_solutions(
        &mut self,
        solutions: &[AnnealingSolution],
    ) -> ClusteringResult<ClusteringResults> {
        let start_time = Instant::now();

        // Convert solutions to solution points
        let solution_points = self.convert_solutions(solutions)?;

        // Extract features if needed
        let featured_points = self.extract_features(solution_points)?;

        // Perform clustering
        let clusters = self.perform_clustering(&featured_points)?;

        // Perform landscape analysis
        let landscape_analysis = self.analyze_landscape(&featured_points, &clusters)?;

        // Perform statistical analysis
        let statistical_summary = self.perform_statistical_analysis(&featured_points, &clusters)?;

        // Calculate overall quality metrics
        let overall_quality = self.calculate_overall_quality(&clusters, &featured_points)?;

        // Generate recommendations
        let recommendations =
            self.generate_recommendations(&clusters, &landscape_analysis, &statistical_summary)?;

        // Update statistics
        self.stats.total_solutions += solutions.len();
        self.stats.total_time += start_time.elapsed();

        Ok(ClusteringResults {
            clusters,
            algorithm: self.config.algorithm.clone(),
            distance_metric: self.config.distance_metric.clone(),
            overall_quality,
            landscape_analysis,
            statistical_summary,
            performance_metrics: ClusteringPerformanceMetrics {
                clustering_time: start_time.elapsed(),
                analysis_time: start_time.elapsed(),
                memory_usage: 0, // Simplified
                scalability_metrics: ScalabilityMetrics {
                    time_complexity: "O(n^2)".to_string(),
                    space_complexity: "O(n^2)".to_string(),
                    scaling_factor: 2.0,
                    parallelization_efficiency: 0.8,
                },
                efficiency_metrics: EfficiencyMetrics {
                    convergence_efficiency: 0.85,
                    resource_utilization: 0.75,
                    quality_time_ratio: 0.9,
                    robustness: 0.8,
                },
            },
            recommendations,
        })
    }

    /// Convert annealing solutions to solution points
    pub fn convert_solutions(
        &self,
        solutions: &[AnnealingSolution],
    ) -> ClusteringResult<Vec<SolutionPoint>> {
        let mut solution_points = Vec::new();

        for (i, solution) in solutions.iter().enumerate() {
            let mut metrics = HashMap::new();
            metrics.insert("energy".to_string(), solution.best_energy);
            metrics.insert("num_evaluations".to_string(), solution.total_sweeps as f64);

            solution_points.push(SolutionPoint {
                solution: solution.best_spins.clone(),
                energy: solution.best_energy,
                metrics,
                metadata: SolutionMetadata {
                    id: i,
                    source: "annealing".to_string(),
                    timestamp: Instant::now(),
                    iterations: solution.total_sweeps,
                    quality_rank: None,
                    is_feasible: true, // Simplified
                },
                features: None,
            });
        }

        Ok(solution_points)
    }

    /// Extract features from solution points
    fn extract_features(
        &self,
        mut solution_points: Vec<SolutionPoint>,
    ) -> ClusteringResult<Vec<SolutionPoint>> {
        match &self.config.feature_extraction {
            FeatureExtractionMethod::Raw => {
                for point in &mut solution_points {
                    point.features = Some(point.solution.iter().map(|&x| f64::from(x)).collect());
                }
            }
            FeatureExtractionMethod::EnergyBased => {
                for point in &mut solution_points {
                    let mut features = vec![point.energy];
                    features.extend(point.solution.iter().map(|&x| f64::from(x)));
                    point.features = Some(features);
                }
            }
            FeatureExtractionMethod::Structural => {
                for point in &mut solution_points {
                    let features = self.extract_structural_features(&point.solution);
                    point.features = Some(features);
                }
            }
            FeatureExtractionMethod::PCA { num_components } => {
                // Simplified PCA implementation
                let features = self.apply_pca(&solution_points, *num_components)?;
                for (point, feature_vec) in solution_points.iter_mut().zip(features.iter()) {
                    point.features = Some(feature_vec.clone());
                }
            }
            _ => {
                // Default to raw features
                for point in &mut solution_points {
                    point.features = Some(point.solution.iter().map(|&x| f64::from(x)).collect());
                }
            }
        }

        Ok(solution_points)
    }

    /// Extract structural features from a solution
    #[must_use]
    pub fn extract_structural_features(&self, solution: &[i8]) -> Vec<f64> {
        let mut features = Vec::new();

        // Basic structural features
        let num_ones = solution.iter().filter(|&&x| x == 1).count() as f64;
        let num_neg_ones = solution.iter().filter(|&&x| x == -1).count() as f64;

        features.push(num_ones);
        features.push(num_neg_ones);
        features.push(num_ones / solution.len() as f64); // Fraction of +1 spins

        // Consecutive patterns
        let mut consecutive_ones = 0;
        let mut consecutive_neg_ones = 0;
        let mut max_consecutive_ones = 0;
        let mut max_consecutive_neg_ones = 0;

        for &spin in solution {
            if spin == 1 {
                consecutive_ones += 1;
                consecutive_neg_ones = 0;
                max_consecutive_ones = max_consecutive_ones.max(consecutive_ones);
            } else {
                consecutive_neg_ones += 1;
                consecutive_ones = 0;
                max_consecutive_neg_ones = max_consecutive_neg_ones.max(consecutive_neg_ones);
            }
        }

        features.push(f64::from(max_consecutive_ones));
        features.push(f64::from(max_consecutive_neg_ones));

        // Transition count
        let transitions = solution
            .windows(2)
            .filter(|window| window[0] != window[1])
            .count() as f64;

        features.push(transitions);

        features
    }

    /// Apply PCA to solution points (simplified implementation)
    fn apply_pca(
        &self,
        solution_points: &[SolutionPoint],
        num_components: usize,
    ) -> ClusteringResult<Vec<Vec<f64>>> {
        if solution_points.is_empty() {
            return Ok(Vec::new());
        }

        let n = solution_points.len();
        let d = solution_points[0].solution.len();

        // Create data matrix
        let mut data = vec![vec![0.0; d]; n];
        for (i, point) in solution_points.iter().enumerate() {
            for (j, &spin) in point.solution.iter().enumerate() {
                data[i][j] = f64::from(spin);
            }
        }

        // Center the data
        let mut means = vec![0.0; d];
        for j in 0..d {
            means[j] = data.iter().map(|row| row[j]).sum::<f64>() / n as f64;
        }

        for i in 0..n {
            for j in 0..d {
                data[i][j] -= means[j];
            }
        }

        // Simplified PCA: just take first num_components dimensions
        let mut pca_data = Vec::new();
        for i in 0..n {
            let mut pca_row = Vec::new();
            for j in 0..num_components.min(d) {
                pca_row.push(data[i][j]);
            }
            pca_data.push(pca_row);
        }

        Ok(pca_data)
    }

    /// Perform clustering on solution points
    fn perform_clustering(
        &self,
        solution_points: &[SolutionPoint],
    ) -> ClusteringResult<Vec<SolutionCluster>> {
        match &self.config.algorithm {
            ClusteringAlgorithm::KMeans { k, max_iterations } => {
                self.kmeans_clustering(solution_points, *k, *max_iterations)
            }
            ClusteringAlgorithm::Hierarchical {
                linkage,
                distance_threshold,
            } => self.hierarchical_clustering(solution_points, linkage, *distance_threshold),
            ClusteringAlgorithm::DBSCAN { eps, min_samples } => {
                self.dbscan_clustering(solution_points, *eps, *min_samples)
            }
            _ => {
                // Default to k-means
                self.kmeans_clustering(solution_points, 5, 100)
            }
        }
    }

    /// K-means clustering implementation
    pub fn kmeans_clustering(
        &self,
        solution_points: &[SolutionPoint],
        k: usize,
        max_iterations: usize,
    ) -> ClusteringResult<Vec<SolutionCluster>> {
        if solution_points.len() < k {
            return Err(ClusteringError::InsufficientData {
                required: k,
                actual: solution_points.len(),
            });
        }

        let n = solution_points.len();
        let features = solution_points
            .iter()
            .map(|p| {
                p.features.as_ref().ok_or_else(|| {
                    ClusteringError::DataError("Solution point missing features".to_string())
                })
            })
            .collect::<Result<Vec<_>, _>>()?;
        let d = features[0].len();

        // Initialize centroids randomly
        let mut rng = match self.config.seed {
            Some(seed) => ChaCha8Rng::seed_from_u64(seed),
            None => ChaCha8Rng::seed_from_u64(thread_rng().random()),
        };

        let mut centroids = Vec::new();
        for _ in 0..k {
            let mut centroid = Vec::new();
            for _ in 0..d {
                centroid.push(rng.random_range(-1.0..1.0));
            }
            centroids.push(centroid);
        }

        let mut assignments = vec![0; n];

        // K-means iterations
        for _iteration in 0..max_iterations {
            let mut changed = false;

            // Assign points to closest centroids
            for (i, feature_vec) in features.iter().enumerate() {
                let mut best_cluster = 0;
                let mut best_distance = f64::INFINITY;

                for (j, centroid) in centroids.iter().enumerate() {
                    let distance = self.calculate_distance(feature_vec, centroid)?;
                    if distance < best_distance {
                        best_distance = distance;
                        best_cluster = j;
                    }
                }

                if assignments[i] != best_cluster {
                    assignments[i] = best_cluster;
                    changed = true;
                }
            }

            // Update centroids
            for j in 0..k {
                let cluster_points: Vec<_> = features
                    .iter()
                    .enumerate()
                    .filter(|(i, _)| assignments[*i] == j)
                    .map(|(_, features)| *features)
                    .collect();

                if !cluster_points.is_empty() {
                    for dim in 0..d {
                        centroids[j][dim] =
                            cluster_points.iter().map(|point| point[dim]).sum::<f64>()
                                / cluster_points.len() as f64;
                    }
                }
            }

            if !changed {
                break;
            }
        }

        // Create clusters
        let mut clusters = Vec::new();
        for cluster_id in 0..k {
            let cluster_solutions: Vec<_> = solution_points
                .iter()
                .enumerate()
                .filter(|(i, _)| assignments[*i] == cluster_id)
                .map(|(_, point)| point.clone())
                .collect();

            if !cluster_solutions.is_empty() {
                let statistics = self.calculate_cluster_statistics(&cluster_solutions);
                let quality_metrics = self
                    .calculate_cluster_quality_metrics(&cluster_solutions, &centroids[cluster_id]);

                clusters.push(SolutionCluster {
                    id: cluster_id,
                    solutions: cluster_solutions,
                    centroid: centroids[cluster_id].clone(),
                    representative: None, // Will be set later
                    statistics,
                    quality_metrics,
                });
            }
        }

        Ok(clusters)
    }

    /// Hierarchical clustering implementation (simplified)
    fn hierarchical_clustering(
        &self,
        solution_points: &[SolutionPoint],
        _linkage: &LinkageType,
        distance_threshold: f64,
    ) -> ClusteringResult<Vec<SolutionCluster>> {
        // Simplified implementation using single linkage
        let n = solution_points.len();
        let mut clusters: Vec<Vec<usize>> = (0..n).map(|i| vec![i]).collect();

        while clusters.len() > 1 {
            let mut min_distance = f64::INFINITY;
            let mut merge_indices = (0, 1);

            // Find closest clusters
            for i in 0..clusters.len() {
                for j in (i + 1)..clusters.len() {
                    let distance = self.calculate_cluster_distance(
                        &clusters[i],
                        &clusters[j],
                        solution_points,
                    )?;
                    if distance < min_distance {
                        min_distance = distance;
                        merge_indices = (i, j);
                    }
                }
            }

            if min_distance > distance_threshold {
                break;
            }

            // Merge clusters
            let (i, j) = merge_indices;
            let mut merged_cluster = clusters[i].clone();
            merged_cluster.extend_from_slice(&clusters[j]);

            // Remove original clusters and add merged cluster
            if i < j {
                clusters.remove(j);
                clusters.remove(i);
            } else {
                clusters.remove(i);
                clusters.remove(j);
            }
            clusters.push(merged_cluster);
        }

        // Convert to SolutionCluster format
        let mut result_clusters = Vec::new();
        for (cluster_id, cluster_indices) in clusters.iter().enumerate() {
            let cluster_solutions: Vec<_> = cluster_indices
                .iter()
                .map(|&i| solution_points[i].clone())
                .collect();

            if !cluster_solutions.is_empty() {
                let centroid = self.calculate_centroid(&cluster_solutions);
                let statistics = self.calculate_cluster_statistics(&cluster_solutions);
                let quality_metrics =
                    self.calculate_cluster_quality_metrics(&cluster_solutions, &centroid);

                result_clusters.push(SolutionCluster {
                    id: cluster_id,
                    solutions: cluster_solutions,
                    centroid,
                    representative: None,
                    statistics,
                    quality_metrics,
                });
            }
        }

        Ok(result_clusters)
    }

    /// DBSCAN clustering implementation (simplified)
    fn dbscan_clustering(
        &self,
        solution_points: &[SolutionPoint],
        eps: f64,
        min_samples: usize,
    ) -> ClusteringResult<Vec<SolutionCluster>> {
        let n = solution_points.len();
        let mut labels = vec![-1i32; n]; // -1 = noise, 0+ = cluster id
        let mut cluster_id = 0;

        for i in 0..n {
            if labels[i] != -1 {
                continue; // Already processed
            }

            let neighbors = self.find_neighbors(i, solution_points, eps)?;

            if neighbors.len() < min_samples {
                labels[i] = -1; // Mark as noise
                continue;
            }

            // Start new cluster
            labels[i] = cluster_id;
            let mut queue = VecDeque::from(neighbors);

            while let Some(j) = queue.pop_front() {
                if labels[j] == -1 {
                    labels[j] = cluster_id; // Change noise to border point
                } else if labels[j] != -1 {
                    continue; // Already in a cluster
                }

                labels[j] = cluster_id;
                let j_neighbors = self.find_neighbors(j, solution_points, eps)?;

                if j_neighbors.len() >= min_samples {
                    for &neighbor in &j_neighbors {
                        if labels[neighbor] == -1 || labels[neighbor] == cluster_id {
                            queue.push_back(neighbor);
                        }
                    }
                }
            }

            cluster_id += 1;
        }

        // Convert to SolutionCluster format
        let mut result_clusters = Vec::new();
        for cid in 0..cluster_id {
            let cluster_solutions: Vec<_> = solution_points
                .iter()
                .enumerate()
                .filter(|(i, _)| labels[*i] == cid)
                .map(|(_, point)| point.clone())
                .collect();

            if !cluster_solutions.is_empty() {
                let centroid = self.calculate_centroid(&cluster_solutions);
                let statistics = self.calculate_cluster_statistics(&cluster_solutions);
                let quality_metrics =
                    self.calculate_cluster_quality_metrics(&cluster_solutions, &centroid);

                result_clusters.push(SolutionCluster {
                    id: cid as usize,
                    solutions: cluster_solutions,
                    centroid,
                    representative: None,
                    statistics,
                    quality_metrics,
                });
            }
        }

        Ok(result_clusters)
    }

    /// Find neighbors within eps distance
    fn find_neighbors(
        &self,
        point_idx: usize,
        solution_points: &[SolutionPoint],
        eps: f64,
    ) -> ClusteringResult<Vec<usize>> {
        let mut neighbors = Vec::new();
        let point_features = solution_points[point_idx]
            .features
            .as_ref()
            .ok_or_else(|| {
                ClusteringError::DataError("Solution point missing features".to_string())
            })?;

        for (i, other_point) in solution_points.iter().enumerate() {
            if i != point_idx {
                let other_features = other_point.features.as_ref().ok_or_else(|| {
                    ClusteringError::DataError("Solution point missing features".to_string())
                })?;
                let distance = self.calculate_distance(point_features, other_features)?;
                if distance <= eps {
                    neighbors.push(i);
                }
            }
        }

        Ok(neighbors)
    }

    /// Calculate distance between cluster indices
    fn calculate_cluster_distance(
        &self,
        cluster1: &[usize],
        cluster2: &[usize],
        solution_points: &[SolutionPoint],
    ) -> ClusteringResult<f64> {
        let mut min_distance = f64::INFINITY;

        for &i in cluster1 {
            for &j in cluster2 {
                let features1 = solution_points[i].features.as_ref().ok_or_else(|| {
                    ClusteringError::DataError("Solution point missing features".to_string())
                })?;
                let features2 = solution_points[j].features.as_ref().ok_or_else(|| {
                    ClusteringError::DataError("Solution point missing features".to_string())
                })?;
                let distance = self.calculate_distance(features1, features2)?;
                min_distance = min_distance.min(distance);
            }
        }

        Ok(min_distance)
    }

    /// Calculate distance between two feature vectors
    pub fn calculate_distance(
        &self,
        features1: &[f64],
        features2: &[f64],
    ) -> ClusteringResult<f64> {
        if features1.len() != features2.len() {
            return Err(ClusteringError::DimensionMismatch {
                expected: features1.len(),
                actual: features2.len(),
            });
        }

        match self.config.distance_metric {
            DistanceMetric::Euclidean => {
                let sum_sq: f64 = features1
                    .iter()
                    .zip(features2.iter())
                    .map(|(a, b)| (a - b).powi(2))
                    .sum();
                Ok(sum_sq.sqrt())
            }
            DistanceMetric::Manhattan => {
                let sum_abs: f64 = features1
                    .iter()
                    .zip(features2.iter())
                    .map(|(a, b)| (a - b).abs())
                    .sum();
                Ok(sum_abs)
            }
            DistanceMetric::Hamming => {
                let diff_count = features1
                    .iter()
                    .zip(features2.iter())
                    .filter(|(a, b)| (*a - *b).abs() > 1e-10)
                    .count();
                Ok(diff_count as f64)
            }
            DistanceMetric::Cosine => {
                let dot_product: f64 = features1
                    .iter()
                    .zip(features2.iter())
                    .map(|(a, b)| a * b)
                    .sum();

                let norm1: f64 = features1.iter().map(|x| x * x).sum::<f64>().sqrt();
                let norm2: f64 = features2.iter().map(|x| x * x).sum::<f64>().sqrt();

                if norm1 > 1e-10 && norm2 > 1e-10 {
                    Ok(1.0 - dot_product / (norm1 * norm2))
                } else {
                    Ok(1.0)
                }
            }
            _ => {
                // Default to Euclidean
                let sum_sq: f64 = features1
                    .iter()
                    .zip(features2.iter())
                    .map(|(a, b)| (a - b).powi(2))
                    .sum();
                Ok(sum_sq.sqrt())
            }
        }
    }

    /// Calculate centroid of a cluster
    fn calculate_centroid(&self, cluster_solutions: &[SolutionPoint]) -> Vec<f64> {
        if cluster_solutions.is_empty() {
            return Vec::new();
        }

        let features_dim = cluster_solutions[0]
            .features
            .as_ref()
            .expect("solution point should have features for centroid calculation")
            .len();
        let mut centroid = vec![0.0; features_dim];

        for solution in cluster_solutions {
            let features = solution
                .features
                .as_ref()
                .expect("solution point should have features for centroid calculation");
            for (i, &value) in features.iter().enumerate() {
                centroid[i] += value;
            }
        }

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

        centroid
    }

    /// Calculate cluster statistics
    fn calculate_cluster_statistics(
        &self,
        cluster_solutions: &[SolutionPoint],
    ) -> ClusterStatistics {
        if cluster_solutions.is_empty() {
            return ClusterStatistics {
                size: 0,
                mean_energy: 0.0,
                energy_std: 0.0,
                min_energy: 0.0,
                max_energy: 0.0,
                intra_cluster_distance: 0.0,
                diameter: 0.0,
                density: 0.0,
            };
        }

        let energies: Vec<f64> = cluster_solutions.iter().map(|s| s.energy).collect();
        let mean_energy = energies.iter().sum::<f64>() / energies.len() as f64;
        let variance = energies
            .iter()
            .map(|e| (e - mean_energy).powi(2))
            .sum::<f64>()
            / energies.len() as f64;
        let energy_std = variance.sqrt();

        let min_energy = energies.iter().fold(f64::INFINITY, |a, &b| a.min(b));
        let max_energy = energies.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));

        // Calculate intra-cluster distance and diameter
        let mut total_distance = 0.0;
        let mut max_distance = 0.0f64;
        let mut distance_count = 0;

        for i in 0..cluster_solutions.len() {
            for j in (i + 1)..cluster_solutions.len() {
                if let (Some(features1), Some(features2)) = (
                    cluster_solutions[i].features.as_ref(),
                    cluster_solutions[j].features.as_ref(),
                ) {
                    if let Ok(distance) = self.calculate_distance(features1, features2) {
                        total_distance += distance;
                        max_distance = max_distance.max(distance);
                        distance_count += 1;
                    }
                }
            }
        }

        let intra_cluster_distance = if distance_count > 0 {
            total_distance / f64::from(distance_count)
        } else {
            0.0
        };

        ClusterStatistics {
            size: cluster_solutions.len(),
            mean_energy,
            energy_std,
            min_energy,
            max_energy,
            intra_cluster_distance,
            diameter: max_distance,
            density: if max_distance > 0.0 {
                cluster_solutions.len() as f64 / max_distance
            } else {
                0.0
            },
        }
    }

    /// Calculate cluster quality metrics
    fn calculate_cluster_quality_metrics(
        &self,
        cluster_solutions: &[SolutionPoint],
        centroid: &[f64],
    ) -> ClusterQualityMetrics {
        let mut inertia = 0.0;

        for solution in cluster_solutions {
            if let Some(features) = solution.features.as_ref() {
                if let Ok(distance) = self.calculate_distance(features, centroid) {
                    inertia += distance * distance;
                }
            }
        }

        ClusterQualityMetrics {
            silhouette_coefficient: 0.5, // Simplified
            inertia,
            calinski_harabasz_index: 1.0, // Simplified
            davies_bouldin_index: 1.0,    // Simplified
            stability: 0.8,               // Simplified
        }
    }

    /// Analyze the solution landscape
    fn analyze_landscape(
        &self,
        solution_points: &[SolutionPoint],
        clusters: &[SolutionCluster],
    ) -> ClusteringResult<LandscapeAnalysis> {
        let energy_statistics = self.calculate_energy_statistics(solution_points);
        let basins = self.detect_energy_basins(solution_points, clusters);
        let connectivity = self.analyze_connectivity(solution_points);
        let multi_modality = self.analyze_multi_modality(solution_points);
        let ruggedness = self.calculate_ruggedness_metrics(solution_points);
        let funnel_analysis = self.analyze_funnel_structure(solution_points, clusters);

        Ok(LandscapeAnalysis {
            energy_statistics,
            basins,
            connectivity,
            multi_modality,
            ruggedness,
            funnel_analysis,
        })
    }

    /// Calculate energy statistics
    #[must_use]
    pub fn calculate_energy_statistics(
        &self,
        solution_points: &[SolutionPoint],
    ) -> EnergyStatistics {
        let energies: Vec<f64> = solution_points.iter().map(|s| s.energy).collect();

        if energies.is_empty() {
            return EnergyStatistics {
                mean: 0.0,
                std_dev: 0.0,
                min: 0.0,
                max: 0.0,
                percentiles: Vec::new(),
                skewness: 0.0,
                kurtosis: 0.0,
                num_distinct_energies: 0,
            };
        }

        let mut sorted_energies = energies.clone();
        sorted_energies.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

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

        let min = sorted_energies[0];
        let max = sorted_energies[sorted_energies.len() - 1];

        // Calculate percentiles
        let percentiles = vec![
            sorted_energies[sorted_energies.len() * 25 / 100],
            sorted_energies[sorted_energies.len() * 50 / 100],
            sorted_energies[sorted_energies.len() * 75 / 100],
        ];

        // Calculate skewness and kurtosis (simplified)
        let skewness = if std_dev > 1e-10 {
            energies
                .iter()
                .map(|e| ((e - mean) / std_dev).powi(3))
                .sum::<f64>()
                / energies.len() as f64
        } else {
            0.0
        };

        let kurtosis = if std_dev > 1e-10 {
            energies
                .iter()
                .map(|e| ((e - mean) / std_dev).powi(4))
                .sum::<f64>()
                / energies.len() as f64
                - 3.0
        } else {
            0.0
        };

        let mut sorted_energies = energies;
        sorted_energies.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        sorted_energies.dedup_by(|a, b| (*a - *b).abs() < 1e-10);
        let num_distinct_energies = sorted_energies.len();

        EnergyStatistics {
            mean,
            std_dev,
            min,
            max,
            percentiles,
            skewness,
            kurtosis,
            num_distinct_energies,
        }
    }

    /// Detect energy basins in the landscape
    fn detect_energy_basins(
        &self,
        _solution_points: &[SolutionPoint],
        clusters: &[SolutionCluster],
    ) -> Vec<EnergyBasin> {
        let mut basins = Vec::new();

        for (basin_id, cluster) in clusters.iter().enumerate() {
            let energies: Vec<f64> = cluster.solutions.iter().map(|s| s.energy).collect();

            if !energies.is_empty() {
                let min_energy = energies.iter().fold(f64::INFINITY, |a, &b| a.min(b));
                let max_energy = energies.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));

                basins.push(EnergyBasin {
                    id: basin_id,
                    solutions: cluster.solutions.iter().map(|s| s.metadata.id).collect(),
                    min_energy,
                    size: cluster.solutions.len(),
                    depth: 0.0, // Would need global minimum to calculate
                    width: max_energy - min_energy,
                    escape_barrier: 0.0, // Simplified
                });
            }
        }

        basins
    }

    /// Analyze connectivity of the solution landscape
    const fn analyze_connectivity(
        &self,
        solution_points: &[SolutionPoint],
    ) -> ConnectivityAnalysis {
        // Simplified connectivity analysis
        ConnectivityAnalysis {
            num_components: 1, // Simplified
            largest_component_size: solution_points.len(),
            average_path_length: 2.0,    // Simplified
            clustering_coefficient: 0.3, // Simplified
            diameter: 10,                // Simplified
        }
    }

    /// Analyze multi-modality of the landscape
    fn analyze_multi_modality(&self, solution_points: &[SolutionPoint]) -> MultiModalityAnalysis {
        // Simplified multi-modality analysis
        let energies: Vec<f64> = solution_points.iter().map(|s| s.energy).collect();

        MultiModalityAnalysis {
            num_modes: 2, // Simplified
            mode_energies: if energies.len() >= 2 {
                vec![energies[0], energies[energies.len() / 2]]
            } else {
                energies
            },
            mode_strengths: vec![0.6, 0.4], // Simplified
            inter_mode_distances: vec![vec![0.0, 1.0], vec![1.0, 0.0]], // Simplified
            multi_modality_index: 0.3,      // Simplified
        }
    }

    /// Calculate ruggedness metrics
    fn calculate_ruggedness_metrics(&self, solution_points: &[SolutionPoint]) -> RuggednessMetrics {
        // Simplified ruggedness calculation
        RuggednessMetrics {
            autocorrelation: vec![1.0, 0.8, 0.6, 0.4, 0.2], // Simplified
            ruggedness_coefficient: 0.5,                    // Simplified
            num_local_optima: solution_points.len() / 10,   // Simplified
            epistasis: 0.3,                                 // Simplified
            neutrality: 0.1,                                // Simplified
        }
    }

    /// Analyze funnel structure
    fn analyze_funnel_structure(
        &self,
        _solution_points: &[SolutionPoint],
        clusters: &[SolutionCluster],
    ) -> FunnelAnalysis {
        // Simplified funnel analysis
        FunnelAnalysis {
            num_funnels: clusters.len(),
            funnel_depths: clusters.iter().map(|c| c.statistics.energy_std).collect(),
            funnel_widths: clusters.iter().map(|c| c.statistics.diameter).collect(),
            global_funnel: Some(0), // Simplified
            competition_index: 0.5, // Simplified
        }
    }

    /// Perform statistical analysis
    fn perform_statistical_analysis(
        &self,
        solution_points: &[SolutionPoint],
        clusters: &[SolutionCluster],
    ) -> ClusteringResult<StatisticalSummary> {
        let cluster_size_distribution = clusters.iter().map(|c| c.statistics.size).collect();

        let energy_distribution = DistributionAnalysis {
            distribution_type: DistributionType::Normal, // Simplified
            parameters: HashMap::from([("mean".to_string(), 0.0), ("std".to_string(), 1.0)]),
            goodness_of_fit: 0.8,
            confidence_intervals: vec![(0.1, 0.9)],
        };

        let convergence_analysis = ConvergenceAnalysis {
            trajectory_clusters: Vec::new(), // Simplified
            convergence_rates: vec![0.1, 0.2, 0.15],
            plateau_analysis: PlateauAnalysis {
                num_plateaus: 2,
                plateau_durations: vec![10, 15],
                plateau_energies: vec![-1.0, -0.5],
                escape_probabilities: vec![0.3, 0.7],
            },
            premature_convergence: false,
            diversity_evolution: vec![1.0, 0.8, 0.6, 0.4, 0.2],
        };

        let correlation_analysis = CorrelationAnalysis {
            variable_correlations: vec![
                vec![1.0; solution_points[0].solution.len()];
                solution_points[0].solution.len()
            ],
            energy_correlations: vec![0.1; solution_points[0].solution.len()],
            significant_correlations: Vec::new(),
            correlation_patterns: Vec::new(),
        };

        let outliers = Vec::new(); // Simplified

        Ok(StatisticalSummary {
            cluster_size_distribution,
            energy_distribution,
            convergence_analysis,
            correlation_analysis,
            outliers,
        })
    }

    /// Calculate overall clustering quality
    fn calculate_overall_quality(
        &self,
        clusters: &[SolutionCluster],
        solution_points: &[SolutionPoint],
    ) -> ClusteringResult<OverallClusteringQuality> {
        let silhouette_score = clusters
            .iter()
            .map(|c| c.quality_metrics.silhouette_coefficient)
            .sum::<f64>()
            / clusters.len() as f64;

        let inter_cluster_separation = self.calculate_inter_cluster_separation(clusters)?;
        let cluster_cohesion = self.calculate_cluster_cohesion(clusters);

        Ok(OverallClusteringQuality {
            silhouette_score,
            adjusted_rand_index: None,
            normalized_mutual_information: None,
            inter_cluster_separation,
            cluster_cohesion,
            num_clusters: clusters.len(),
            optimal_num_clusters: self.estimate_optimal_clusters(solution_points)?,
        })
    }

    /// Calculate inter-cluster separation
    fn calculate_inter_cluster_separation(
        &self,
        clusters: &[SolutionCluster],
    ) -> ClusteringResult<f64> {
        if clusters.len() < 2 {
            return Ok(0.0);
        }

        let mut total_separation = 0.0;
        let mut count = 0;

        for i in 0..clusters.len() {
            for j in (i + 1)..clusters.len() {
                let distance =
                    self.calculate_distance(&clusters[i].centroid, &clusters[j].centroid)?;
                total_separation += distance;
                count += 1;
            }
        }

        Ok(total_separation / f64::from(count))
    }

    /// Calculate cluster cohesion
    fn calculate_cluster_cohesion(&self, clusters: &[SolutionCluster]) -> f64 {
        if clusters.is_empty() {
            return 0.0;
        }

        clusters
            .iter()
            .map(|c| 1.0 / (1.0 + c.statistics.intra_cluster_distance))
            .sum::<f64>()
            / clusters.len() as f64
    }

    /// Estimate optimal number of clusters
    fn estimate_optimal_clusters(
        &self,
        solution_points: &[SolutionPoint],
    ) -> ClusteringResult<usize> {
        // Simplified elbow method
        let max_k = solution_points.len().min(10);
        let mut inertias = Vec::new();

        for k in 1..=max_k {
            if let Ok(clusters) = self.kmeans_clustering(solution_points, k, 50) {
                let total_inertia: f64 = clusters.iter().map(|c| c.quality_metrics.inertia).sum();
                inertias.push(total_inertia);
            }
        }

        // Find elbow (simplified)
        let optimal_k = if inertias.len() >= 3 {
            let mut max_diff = 0.0;
            let mut optimal = 1;

            for i in 1..inertias.len() - 1 {
                let diff = 2.0f64.mul_add(-inertias[i], inertias[i - 1]) + inertias[i + 1];
                if diff > max_diff {
                    max_diff = diff;
                    optimal = i + 1;
                }
            }
            optimal
        } else {
            inertias.len()
        };

        Ok(optimal_k)
    }

    /// Generate optimization recommendations
    fn generate_recommendations(
        &self,
        clusters: &[SolutionCluster],
        landscape_analysis: &LandscapeAnalysis,
        _statistical_summary: &StatisticalSummary,
    ) -> ClusteringResult<Vec<OptimizationRecommendation>> {
        let mut recommendations = Vec::new();

        // Recommendation based on cluster quality
        if clusters
            .iter()
            .any(|c| c.quality_metrics.silhouette_coefficient < 0.3)
        {
            recommendations.push(OptimizationRecommendation {
                recommendation_type: RecommendationType::ParameterTuning,
                description: "Low cluster quality detected. Consider tuning annealing parameters or using different initialization strategies.".to_string(),
                expected_improvement: 0.2,
                difficulty: DifficultyLevel::Easy,
                priority: PriorityLevel::High,
                evidence: vec!["Low silhouette coefficients in multiple clusters".to_string()],
            });
        }

        // Recommendation based on energy landscape
        if landscape_analysis.multi_modality.num_modes > 3 {
            recommendations.push(OptimizationRecommendation {
                recommendation_type: RecommendationType::MultiStart,
                description: "Multiple modes detected in energy landscape. Consider using multi-start optimization or parallel runs.".to_string(),
                expected_improvement: 0.3,
                difficulty: DifficultyLevel::Moderate,
                priority: PriorityLevel::Medium,
                evidence: vec![format!("{} modes detected", landscape_analysis.multi_modality.num_modes)],
            });
        }

        // Recommendation based on cluster sizes
        let cluster_sizes: Vec<usize> = clusters.iter().map(|c| c.statistics.size).collect();
        let size_variance = cluster_sizes
            .iter()
            .map(|&size| {
                (size as f64
                    - cluster_sizes.iter().sum::<usize>() as f64 / cluster_sizes.len() as f64)
                    .powi(2)
            })
            .sum::<f64>()
            / cluster_sizes.len() as f64;

        if size_variance > 100.0 {
            recommendations.push(OptimizationRecommendation {
                recommendation_type: RecommendationType::AlgorithmModification,
                description: "Highly unbalanced cluster sizes suggest potential convergence issues. Consider adjusting cooling schedule or using adaptive algorithms.".to_string(),
                expected_improvement: 0.15,
                difficulty: DifficultyLevel::Moderate,
                priority: PriorityLevel::Medium,
                evidence: vec![format!("Cluster size variance: {:.2}", size_variance)],
            });
        }

        Ok(recommendations)
    }
}