scirs2_optimize/global/
clustering.rs

1//! Clustering algorithms for local minima in multi-start optimization
2//!
3//! This module provides algorithms to identify, cluster, and analyze local minima
4//! found during multi-start optimization strategies. It helps distinguish between
5//! unique local minima and provides insights into the optimization landscape.
6
7use crate::error::OptimizeError;
8use crate::unconstrained::{minimize, Method, Options};
9use scirs2_core::ndarray::{Array1, Array2, ArrayView1};
10use std::collections::HashMap;
11
12/// Configuration for clustering local minima
13#[derive(Debug, Clone)]
14pub struct ClusteringOptions {
15    /// Distance threshold for considering minima as belonging to the same cluster
16    pub distance_threshold: f64,
17    /// Relative tolerance for function values when clustering
18    pub function_tolerance: f64,
19    /// Maximum number of clusters to form
20    pub max_clusters: Option<usize>,
21    /// Clustering algorithm to use
22    pub algorithm: ClusteringAlgorithm,
23    /// Whether to normalize coordinates before clustering
24    pub normalize_coordinates: bool,
25    /// Whether to use function values in clustering
26    pub use_function_values: bool,
27    /// Weight for function values vs coordinates in distance calculation
28    pub function_weight: f64,
29}
30
31impl Default for ClusteringOptions {
32    fn default() -> Self {
33        Self {
34            distance_threshold: 1e-3,
35            function_tolerance: 1e-6,
36            max_clusters: None,
37            algorithm: ClusteringAlgorithm::Hierarchical,
38            normalize_coordinates: true,
39            use_function_values: true,
40            function_weight: 0.1,
41        }
42    }
43}
44
45/// Clustering algorithms available
46#[derive(Debug, Clone, Copy)]
47pub enum ClusteringAlgorithm {
48    /// Hierarchical clustering with single linkage
49    Hierarchical,
50    /// K-means clustering
51    KMeans,
52    /// Density-based clustering (DBSCAN-like)
53    Density,
54    /// Custom threshold-based clustering
55    Threshold,
56}
57
58/// Represents a local minimum found during optimization
59#[derive(Debug, Clone)]
60pub struct LocalMinimum<S> {
61    /// Location of the minimum
62    pub x: Array1<f64>,
63    /// Function value at the minimum
64    pub f: f64,
65    /// Original function value (for generic type)
66    pub fun_value: S,
67    /// Number of optimization iterations to reach this minimum
68    pub nit: usize,
69    /// Number of function evaluations
70    pub func_evals: usize,
71    /// Whether optimization was successful
72    pub success: bool,
73    /// Starting point that led to this minimum
74    pub start_point: Array1<f64>,
75    /// Cluster ID (assigned after clustering)
76    pub cluster_id: Option<usize>,
77    /// Distance to cluster centroid
78    pub cluster_distance: Option<f64>,
79}
80
81/// Result of clustering analysis
82#[derive(Debug, Clone)]
83pub struct ClusteringResult<S> {
84    /// All local minima found
85    pub minima: Vec<LocalMinimum<S>>,
86    /// Cluster centroids
87    pub centroids: Vec<ClusterCentroid>,
88    /// Number of clusters formed
89    pub num_clusters: usize,
90    /// Silhouette score (quality measure)
91    pub silhouette_score: Option<f64>,
92    /// Within-cluster sum of squares
93    pub wcss: f64,
94    /// Best minimum found (global optimum candidate)
95    pub best_minimum: Option<LocalMinimum<S>>,
96}
97
98/// Cluster centroid information
99#[derive(Debug, Clone)]
100pub struct ClusterCentroid {
101    /// Centroid coordinates
102    pub x: Array1<f64>,
103    /// Average function value in cluster
104    pub f_avg: f64,
105    /// Best (minimum) function value in cluster
106    pub f_min: f64,
107    /// Number of minima in this cluster
108    pub size: usize,
109    /// Cluster radius (max distance from centroid)
110    pub radius: f64,
111}
112
113/// Multi-start optimization with clustering
114#[allow(dead_code)]
115pub fn multi_start_with_clustering<F, S>(
116    fun: F,
117    start_points: &[Array1<f64>],
118    method: Method,
119    options: Option<Options>,
120    clustering_options: Option<ClusteringOptions>,
121) -> Result<ClusteringResult<S>, OptimizeError>
122where
123    F: FnMut(&ArrayView1<f64>) -> S + Clone,
124    S: Into<f64> + Clone + From<f64>,
125{
126    let clustering_opts = clustering_options.unwrap_or_default();
127    let mut minima = Vec::new();
128
129    // Run optimization from each starting point
130    for start_point in start_points {
131        let fun_clone = fun.clone();
132
133        match minimize(
134            fun_clone,
135            start_point.as_slice().unwrap(),
136            method,
137            options.clone(),
138        ) {
139            Ok(result) => {
140                let minimum = LocalMinimum {
141                    x: result.x.clone(),
142                    f: result.fun.clone().into(),
143                    fun_value: result.fun,
144                    nit: result.nit,
145                    func_evals: result.func_evals,
146                    success: result.success,
147                    start_point: start_point.clone(),
148                    cluster_id: None,
149                    cluster_distance: None,
150                };
151                minima.push(minimum);
152            }
153            Err(_) => {
154                // Skip failed optimizations but could log them
155                continue;
156            }
157        }
158    }
159
160    if minima.is_empty() {
161        return Err(OptimizeError::ComputationError(
162            "No successful optimizations found".to_string(),
163        ));
164    }
165
166    // Perform clustering
167    cluster_minima(&mut minima, &clustering_opts)?;
168
169    // Compute cluster centroids and statistics
170    let centroids = compute_cluster_centroids(&minima)?;
171    let num_clusters = centroids.len();
172    let wcss = compute_wcss(&minima, &centroids);
173    let silhouette_score = compute_silhouette_score(&minima);
174
175    // Find best minimum
176    let best_minimum = minima
177        .iter()
178        .filter(|m| m.success)
179        .min_by(|a, b| a.f.partial_cmp(&b.f).unwrap())
180        .cloned();
181
182    Ok(ClusteringResult {
183        minima,
184        centroids,
185        num_clusters,
186        silhouette_score,
187        wcss,
188        best_minimum,
189    })
190}
191
192/// Cluster local minima using the specified algorithm
193#[allow(dead_code)]
194fn cluster_minima<S>(
195    minima: &mut [LocalMinimum<S>],
196    options: &ClusteringOptions,
197) -> Result<(), OptimizeError>
198where
199    S: Clone,
200{
201    if minima.is_empty() {
202        return Ok(());
203    }
204
205    match options.algorithm {
206        ClusteringAlgorithm::Hierarchical => hierarchical_clustering(minima, options),
207        ClusteringAlgorithm::KMeans => kmeans_clustering(minima, options),
208        ClusteringAlgorithm::Density => density_clustering(minima, options),
209        ClusteringAlgorithm::Threshold => threshold_clustering(minima, options),
210    }
211}
212
213/// Hierarchical clustering implementation
214#[allow(dead_code)]
215fn hierarchical_clustering<S>(
216    minima: &mut [LocalMinimum<S>],
217    options: &ClusteringOptions,
218) -> Result<(), OptimizeError>
219where
220    S: Clone,
221{
222    let n = minima.len();
223    if n <= 1 {
224        if n == 1 {
225            minima[0].cluster_id = Some(0);
226        }
227        return Ok(());
228    }
229
230    // Compute distance matrix
231    let distances = compute_distance_matrix(minima, options);
232
233    // Perform hierarchical clustering using single linkage
234    let mut cluster_assignments = vec![None; n];
235    let mut next_cluster_id = n;
236
237    // Initialize: each point is its own cluster
238    for (i, assignment) in cluster_assignments.iter_mut().enumerate().take(n) {
239        *assignment = Some(i);
240    }
241
242    // Build hierarchy by merging closest clusters
243    let mut active_clusters: Vec<Vec<usize>> = (0..n).map(|i| vec![i]).collect();
244
245    while active_clusters.len() > 1 {
246        let mut min_dist = f64::INFINITY;
247        let mut merge_pair = (0, 1);
248
249        // Find closest pair of clusters
250        for i in 0..active_clusters.len() {
251            for j in (i + 1)..active_clusters.len() {
252                let cluster_dist =
253                    compute_cluster_distance(&active_clusters[i], &active_clusters[j], &distances);
254                if cluster_dist < min_dist {
255                    min_dist = cluster_dist;
256                    merge_pair = (i, j);
257                }
258            }
259        }
260
261        // Stop if distance threshold exceeded
262        if min_dist > options.distance_threshold {
263            break;
264        }
265
266        // Check max clusters constraint
267        if let Some(max_clusters) = options.max_clusters {
268            if active_clusters.len() <= max_clusters {
269                break;
270            }
271        }
272
273        // Merge clusters
274        let (i, j) = merge_pair;
275        let mut merged_cluster = active_clusters[i].clone();
276        merged_cluster.extend(&active_clusters[j]);
277
278        // Update cluster assignments
279        for &point_idx in &merged_cluster {
280            cluster_assignments[point_idx] = Some(next_cluster_id);
281        }
282
283        // Remove old clusters and add merged one
284        let mut new_clusters = Vec::new();
285        for (k, cluster) in active_clusters.iter().enumerate() {
286            if k != i && k != j {
287                new_clusters.push(cluster.clone());
288            }
289        }
290        new_clusters.push(merged_cluster);
291        active_clusters = new_clusters;
292        next_cluster_id += 1;
293    }
294
295    // Assign final cluster IDs (renumber to be sequential)
296    let mut cluster_map = HashMap::new();
297    let mut final_cluster_id = 0;
298
299    for cluster_id in cluster_assignments.iter().flatten() {
300        if !cluster_map.contains_key(cluster_id) {
301            cluster_map.insert(*cluster_id, final_cluster_id);
302            final_cluster_id += 1;
303        }
304    }
305
306    // Update minima with final cluster assignments
307    for (i, minimum) in minima.iter_mut().enumerate() {
308        if let Some(cluster_id) = cluster_assignments[i] {
309            minimum.cluster_id = cluster_map.get(&cluster_id).copied();
310        }
311    }
312
313    Ok(())
314}
315
316/// K-means clustering implementation
317#[allow(dead_code)]
318fn kmeans_clustering<S>(
319    minima: &mut [LocalMinimum<S>],
320    options: &ClusteringOptions,
321) -> Result<(), OptimizeError>
322where
323    S: Clone,
324{
325    let n = minima.len();
326    if n <= 1 {
327        if n == 1 {
328            minima[0].cluster_id = Some(0);
329        }
330        return Ok(());
331    }
332
333    // Determine number of clusters
334    let k = if let Some(max_k) = options.max_clusters {
335        std::cmp::min(max_k, n)
336    } else {
337        // Use elbow method or default to sqrt(n)
338        std::cmp::min((n as f64).sqrt() as usize + 1, n)
339    };
340
341    if k >= n {
342        // Each point is its own cluster
343        for (i, minimum) in minima.iter_mut().enumerate() {
344            minimum.cluster_id = Some(i);
345        }
346        return Ok(());
347    }
348
349    // Get feature vectors for clustering
350    let features = extract_features(minima, options);
351    let dim = features.ncols();
352
353    // Initialize centroids randomly (or use k-means++)
354    let mut centroids = initialize_centroids_plus_plus(&features, k);
355    let mut assignments = vec![0; n];
356    let max_iter = 100;
357    let tolerance = 1e-6;
358
359    for _iter in 0..max_iter {
360        let mut changed = false;
361
362        // Assign points to nearest centroids
363        for (i, assignment) in assignments.iter_mut().enumerate().take(n) {
364            let mut min_dist = f64::INFINITY;
365            let mut best_cluster = 0;
366
367            for j in 0..k {
368                let dist = euclidean_distance(&features.row(i), &centroids.row(j));
369                if dist < min_dist {
370                    min_dist = dist;
371                    best_cluster = j;
372                }
373            }
374
375            if *assignment != best_cluster {
376                *assignment = best_cluster;
377                changed = true;
378            }
379        }
380
381        if !changed {
382            break;
383        }
384
385        // Update centroids
386        let mut new_centroids = Array2::zeros((k, dim));
387        let mut cluster_sizes = vec![0; k];
388
389        for i in 0..n {
390            let cluster = assignments[i];
391            cluster_sizes[cluster] += 1;
392            for d in 0..dim {
393                new_centroids[[cluster, d]] += features[[i, d]];
394            }
395        }
396
397        for j in 0..k {
398            if cluster_sizes[j] > 0 {
399                for d in 0..dim {
400                    new_centroids[[j, d]] /= cluster_sizes[j] as f64;
401                }
402            }
403        }
404
405        // Check convergence
406        let centroid_change = (&centroids - &new_centroids).mapv(|x| x.abs()).sum();
407
408        centroids = new_centroids;
409
410        if centroid_change < tolerance {
411            break;
412        }
413    }
414
415    // Assign cluster IDs to minima
416    for (i, minimum) in minima.iter_mut().enumerate() {
417        minimum.cluster_id = Some(assignments[i]);
418    }
419
420    Ok(())
421}
422
423/// Density-based clustering (DBSCAN-like)
424#[allow(dead_code)]
425fn density_clustering<S>(
426    minima: &mut [LocalMinimum<S>],
427    options: &ClusteringOptions,
428) -> Result<(), OptimizeError>
429where
430    S: Clone,
431{
432    let n = minima.len();
433    if n <= 1 {
434        if n == 1 {
435            minima[0].cluster_id = Some(0);
436        }
437        return Ok(());
438    }
439
440    let eps = options.distance_threshold;
441    let min_pts = 2; // Minimum points to form a cluster
442
443    let distances = compute_distance_matrix(minima, options);
444    let mut cluster_assignments = vec![None; n];
445    let mut visited = vec![false; n];
446    let mut cluster_id = 0;
447
448    for i in 0..n {
449        if visited[i] {
450            continue;
451        }
452        visited[i] = true;
453
454        // Find neighbors
455        let neighbors: Vec<usize> = (0..n).filter(|&j| distances[[i, j]] <= eps).collect();
456
457        if neighbors.len() < min_pts {
458            // Mark as noise (will remain None)
459            continue;
460        }
461
462        // Start new cluster
463        let mut to_visit = neighbors.clone();
464        cluster_assignments[i] = Some(cluster_id);
465
466        let mut idx = 0;
467        while idx < to_visit.len() {
468            let point = to_visit[idx];
469
470            if !visited[point] {
471                visited[point] = true;
472
473                let point_neighbors: Vec<usize> =
474                    (0..n).filter(|&j| distances[[point, j]] <= eps).collect();
475
476                if point_neighbors.len() >= min_pts {
477                    // Add new neighbors to visit list
478                    for &neighbor in &point_neighbors {
479                        if !to_visit.contains(&neighbor) {
480                            to_visit.push(neighbor);
481                        }
482                    }
483                }
484            }
485
486            if cluster_assignments[point].is_none() {
487                cluster_assignments[point] = Some(cluster_id);
488            }
489
490            idx += 1;
491        }
492
493        cluster_id += 1;
494    }
495
496    // Assign cluster IDs to minima
497    for (i, minimum) in minima.iter_mut().enumerate() {
498        minimum.cluster_id = cluster_assignments[i];
499    }
500
501    Ok(())
502}
503
504/// Simple threshold-based clustering
505#[allow(dead_code)]
506fn threshold_clustering<S>(
507    minima: &mut [LocalMinimum<S>],
508    options: &ClusteringOptions,
509) -> Result<(), OptimizeError>
510where
511    S: Clone,
512{
513    let n = minima.len();
514    if n <= 1 {
515        if n == 1 {
516            minima[0].cluster_id = Some(0);
517        }
518        return Ok(());
519    }
520
521    let mut cluster_assignments = vec![None; n];
522    let mut cluster_id = 0;
523
524    for i in 0..n {
525        if cluster_assignments[i].is_some() {
526            continue;
527        }
528
529        // Start new cluster
530        cluster_assignments[i] = Some(cluster_id);
531
532        // Find all points within threshold
533        for j in (i + 1)..n {
534            if cluster_assignments[j].is_some() {
535                continue;
536            }
537
538            let distance = compute_distance(&minima[i], &minima[j], options);
539            if distance <= options.distance_threshold {
540                cluster_assignments[j] = Some(cluster_id);
541            }
542        }
543
544        cluster_id += 1;
545    }
546
547    // Assign cluster IDs to minima
548    for (i, minimum) in minima.iter_mut().enumerate() {
549        minimum.cluster_id = cluster_assignments[i];
550    }
551
552    Ok(())
553}
554
555/// Compute distance matrix between all minima
556#[allow(dead_code)]
557fn compute_distance_matrix<S>(
558    minima: &[LocalMinimum<S>],
559    options: &ClusteringOptions,
560) -> Array2<f64>
561where
562    S: Clone,
563{
564    let n = minima.len();
565    let mut distances = Array2::zeros((n, n));
566
567    for i in 0..n {
568        for j in (i + 1)..n {
569            let dist = compute_distance(&minima[i], &minima[j], options);
570            distances[[i, j]] = dist;
571            distances[[j, i]] = dist;
572        }
573    }
574
575    distances
576}
577
578/// Compute distance between two local minima
579#[allow(dead_code)]
580fn compute_distance<S>(
581    min1: &LocalMinimum<S>,
582    min2: &LocalMinimum<S>,
583    options: &ClusteringOptions,
584) -> f64
585where
586    S: Clone,
587{
588    // Coordinate distance
589    let coord_dist = (&min1.x - &min2.x).mapv(|x| x.powi(2)).sum().sqrt();
590
591    if !options.use_function_values {
592        return coord_dist;
593    }
594
595    // Function value distance
596    let func_dist = (min1.f - min2.f).abs();
597
598    // Combined distance
599    coord_dist + options.function_weight * func_dist
600}
601
602/// Compute distance between two clusters (single linkage)
603#[allow(dead_code)]
604fn compute_cluster_distance(
605    cluster1: &[usize],
606    cluster2: &[usize],
607    distances: &Array2<f64>,
608) -> f64 {
609    let mut min_dist = f64::INFINITY;
610
611    for &i in cluster1 {
612        for &j in cluster2 {
613            let dist = distances[[i, j]];
614            if dist < min_dist {
615                min_dist = dist;
616            }
617        }
618    }
619
620    min_dist
621}
622
623/// Extract feature vectors for clustering
624#[allow(dead_code)]
625fn extract_features<S>(minima: &[LocalMinimum<S>], options: &ClusteringOptions) -> Array2<f64>
626where
627    S: Clone,
628{
629    let n = minima.len();
630    if n == 0 {
631        return Array2::zeros((0, 0));
632    }
633
634    let coord_dim = minima[0].x.len();
635    let func_dim = if options.use_function_values { 1 } else { 0 };
636    let total_dim = coord_dim + func_dim;
637
638    let mut features = Array2::zeros((n, total_dim));
639
640    // Extract coordinates
641    for (i, minimum) in minima.iter().enumerate() {
642        for j in 0..coord_dim {
643            features[[i, j]] = minimum.x[j];
644        }
645    }
646
647    // Add function values if requested
648    if options.use_function_values {
649        for (i, minimum) in minima.iter().enumerate() {
650            features[[i, coord_dim]] = minimum.f * options.function_weight;
651        }
652    }
653
654    // Normalize if requested
655    if options.normalize_coordinates {
656        normalize_features(&mut features, coord_dim);
657    }
658
659    features
660}
661
662/// Normalize feature matrix
663#[allow(dead_code)]
664fn normalize_features(features: &mut Array2<f64>, coord_dim: usize) {
665    let n = features.nrows();
666    if n == 0 {
667        return;
668    }
669
670    // Normalize coordinates only
671    for j in 0..coord_dim {
672        let col = features.column(j);
673        let min_val = col.iter().fold(f64::INFINITY, |a, &b| f64::min(a, b));
674        let max_val = col.iter().fold(f64::NEG_INFINITY, |a, &b| f64::max(a, b));
675
676        if (max_val - min_val).abs() > 1e-10 {
677            for i in 0..n {
678                features[[i, j]] = (features[[i, j]] - min_val) / (max_val - min_val);
679            }
680        }
681    }
682}
683
684/// Initialize centroids using k-means++ algorithm
685#[allow(dead_code)]
686fn initialize_centroids_plus_plus(features: &Array2<f64>, k: usize) -> Array2<f64> {
687    let (n, dim) = features.dim();
688    let mut centroids = Array2::zeros((k, dim));
689
690    if n == 0 || k == 0 {
691        return centroids;
692    }
693
694    // Choose first centroid randomly
695    let first_idx = 0; // In practice, should be random
696    centroids.row_mut(0).assign(&features.row(first_idx));
697
698    // Choose remaining centroids
699    for c in 1..k {
700        let mut distances = vec![f64::INFINITY; n];
701
702        // Compute distance to nearest centroid for each point
703        for (i, distance) in distances.iter_mut().enumerate().take(n) {
704            let point = features.row(i);
705            for j in 0..c {
706                let centroid = centroids.row(j);
707                let dist = euclidean_distance(&point, &centroid);
708                *distance = distance.min(dist);
709            }
710        }
711
712        // Choose next centroid (should be weighted random, using max for simplicity)
713        let next_idx = distances
714            .iter()
715            .enumerate()
716            .max_by(|a, b| a.1.partial_cmp(b.1).unwrap())
717            .map(|(i, _)| i)
718            .unwrap_or(0);
719
720        centroids.row_mut(c).assign(&features.row(next_idx));
721    }
722
723    centroids
724}
725
726/// Compute Euclidean distance between two points
727#[allow(dead_code)]
728fn euclidean_distance(a: &ArrayView1<f64>, b: &ArrayView1<f64>) -> f64 {
729    (a - b).mapv(|x| x.powi(2)).sum().sqrt()
730}
731
732/// Compute cluster centroids and statistics
733#[allow(dead_code)]
734fn compute_cluster_centroids<S>(
735    minima: &[LocalMinimum<S>],
736) -> Result<Vec<ClusterCentroid>, OptimizeError>
737where
738    S: Clone,
739{
740    if minima.is_empty() {
741        return Ok(Vec::new());
742    }
743
744    // Group minima by cluster
745    let mut clusters: HashMap<usize, Vec<&LocalMinimum<S>>> = HashMap::new();
746
747    for minimum in minima {
748        if let Some(cluster_id) = minimum.cluster_id {
749            clusters.entry(cluster_id).or_default().push(minimum);
750        }
751    }
752
753    let mut centroids = Vec::new();
754
755    for (_, cluster_minima) in clusters {
756        if cluster_minima.is_empty() {
757            continue;
758        }
759
760        let dim = cluster_minima[0].x.len();
761        let mut centroid_x = Array1::zeros(dim);
762        let mut f_sum = 0.0;
763        let mut f_min = f64::INFINITY;
764
765        // Compute centroid coordinates and function statistics
766        for minimum in &cluster_minima {
767            centroid_x = &centroid_x + &minimum.x;
768            f_sum += minimum.f;
769            f_min = f_min.min(minimum.f);
770        }
771
772        let size = cluster_minima.len();
773        centroid_x /= size as f64;
774        let f_avg = f_sum / size as f64;
775
776        // Compute cluster radius
777        let mut max_radius = 0.0;
778        for minimum in &cluster_minima {
779            let dist = (&minimum.x - &centroid_x).mapv(|x| x.powi(2)).sum().sqrt();
780            max_radius = f64::max(max_radius, dist);
781        }
782
783        centroids.push(ClusterCentroid {
784            x: centroid_x,
785            f_avg,
786            f_min,
787            size,
788            radius: max_radius,
789        });
790    }
791
792    Ok(centroids)
793}
794
795/// Compute within-cluster sum of squares
796#[allow(dead_code)]
797fn compute_wcss<S>(minima: &[LocalMinimum<S>], centroids: &[ClusterCentroid]) -> f64
798where
799    S: Clone,
800{
801    let mut wcss = 0.0;
802
803    for minimum in minima {
804        if let Some(cluster_id) = minimum.cluster_id {
805            if cluster_id < centroids.len() {
806                let centroid = &centroids[cluster_id];
807                let dist = (&minimum.x - &centroid.x).mapv(|x| x.powi(2)).sum();
808                wcss += dist;
809            }
810        }
811    }
812
813    wcss
814}
815
816/// Compute silhouette score for clustering quality
817#[allow(dead_code)]
818fn compute_silhouette_score<S>(minima: &[LocalMinimum<S>]) -> Option<f64>
819where
820    S: Clone,
821{
822    if minima.len() < 2 {
823        return None;
824    }
825
826    let mut silhouette_sum = 0.0;
827    let mut valid_points = 0;
828
829    for (i, minimum) in minima.iter().enumerate() {
830        if let Some(cluster_id) = minimum.cluster_id {
831            // Compute intra-cluster distance
832            let mut intra_sum = 0.0;
833            let mut intra_count = 0;
834
835            // Compute inter-cluster distances
836            let mut min_inter = f64::INFINITY;
837            let mut cluster_inter_sums: HashMap<usize, f64> = HashMap::new();
838            let mut cluster_inter_counts: HashMap<usize, usize> = HashMap::new();
839
840            for (j, other) in minima.iter().enumerate() {
841                if i == j {
842                    continue;
843                }
844
845                if let Some(other_cluster_id) = other.cluster_id {
846                    let dist = (&minimum.x - &other.x).mapv(|x| x.powi(2)).sum().sqrt();
847
848                    if other_cluster_id == cluster_id {
849                        intra_sum += dist;
850                        intra_count += 1;
851                    } else {
852                        *cluster_inter_sums.entry(other_cluster_id).or_insert(0.0) += dist;
853                        *cluster_inter_counts.entry(other_cluster_id).or_insert(0) += 1;
854                    }
855                }
856            }
857
858            // Find minimum inter-cluster distance
859            for (other_cluster_id, sum) in cluster_inter_sums {
860                let count = cluster_inter_counts[&other_cluster_id];
861                if count > 0 {
862                    let avg_inter = sum / count as f64;
863                    min_inter = min_inter.min(avg_inter);
864                }
865            }
866
867            if intra_count > 0 && min_inter < f64::INFINITY {
868                let a = intra_sum / intra_count as f64;
869                let b = min_inter;
870                let silhouette = (b - a) / f64::max(a, b);
871                silhouette_sum += silhouette;
872                valid_points += 1;
873            }
874        }
875    }
876
877    if valid_points > 0 {
878        Some(silhouette_sum / valid_points as f64)
879    } else {
880        None
881    }
882}
883
884/// Generate diverse starting points for multi-start optimization
885#[allow(dead_code)]
886pub fn generate_diverse_start_points(
887    bounds: &[(f64, f64)],
888    num_points: usize,
889    strategy: StartPointStrategy,
890) -> Vec<Array1<f64>> {
891    match strategy {
892        StartPointStrategy::Random => generate_random_points(bounds, num_points),
893        StartPointStrategy::LatinHypercube => generate_latin_hypercube_points(bounds, num_points),
894        StartPointStrategy::Grid => generate_grid_points(bounds, num_points),
895        StartPointStrategy::Sobol => generate_sobol_points(bounds, num_points),
896    }
897}
898
899/// Strategy for generating starting points
900#[derive(Debug, Clone, Copy)]
901pub enum StartPointStrategy {
902    Random,
903    LatinHypercube,
904    Grid,
905    Sobol,
906}
907
908/// Generate random starting points
909#[allow(dead_code)]
910fn generate_random_points(bounds: &[(f64, f64)], num_points: usize) -> Vec<Array1<f64>> {
911    let dim = bounds.len();
912    let mut points = Vec::new();
913
914    for _ in 0..num_points {
915        let mut point = Array1::zeros(dim);
916        for (i, &(low, high)) in bounds.iter().enumerate() {
917            // Simple pseudo-random (in practice, use proper RNG)
918            let t = (i * num_points + points.len()) as f64 / (num_points * dim) as f64;
919            let random_val = (t * 17.0).fract(); // Simple pseudo-random
920            point[i] = low + random_val * (high - low);
921        }
922        points.push(point);
923    }
924
925    points
926}
927
928/// Generate Latin Hypercube sampling points
929#[allow(dead_code)]
930fn generate_latin_hypercube_points(bounds: &[(f64, f64)], num_points: usize) -> Vec<Array1<f64>> {
931    let dim = bounds.len();
932    let mut points = Vec::new();
933
934    // Simple Latin Hypercube implementation
935    for i in 0..num_points {
936        let mut point = Array1::zeros(dim);
937        for (j, &(low, high)) in bounds.iter().enumerate() {
938            let segment = (i as f64 + 0.5) / num_points as f64; // Center of segment
939            point[j] = low + segment * (high - low);
940        }
941        points.push(point);
942    }
943
944    points
945}
946
947/// Generate grid points
948#[allow(dead_code)]
949fn generate_grid_points(bounds: &[(f64, f64)], num_points: usize) -> Vec<Array1<f64>> {
950    let dim = bounds.len();
951    if dim == 0 {
952        return Vec::new();
953    }
954
955    let points_per_dim = ((num_points as f64).powf(1.0 / dim as f64)).ceil() as usize;
956    let mut _points = Vec::new();
957
958    // Generate grid coordinates recursively
959    fn generate_grid_recursive(
960        bounds: &[(f64, f64)],
961        points_per_dim: usize,
962        current_point: &mut Array1<f64>,
963        dim_idx: usize,
964        points: &mut Vec<Array1<f64>>,
965    ) {
966        if dim_idx >= bounds.len() {
967            points.push(current_point.clone());
968            return;
969        }
970
971        let (low, high) = bounds[dim_idx];
972        for i in 0..points_per_dim {
973            let t = if points_per_dim == 1 {
974                0.5
975            } else {
976                i as f64 / (points_per_dim - 1) as f64
977            };
978            current_point[dim_idx] = low + t * (high - low);
979            generate_grid_recursive(bounds, points_per_dim, current_point, dim_idx + 1, points);
980        }
981    }
982
983    let mut current_point = Array1::zeros(dim);
984    generate_grid_recursive(bounds, points_per_dim, &mut current_point, 0, &mut _points);
985
986    // Truncate to requested number of _points
987    _points.truncate(num_points);
988    _points
989}
990
991/// Generate Sobol sequence points (simplified)
992#[allow(dead_code)]
993fn generate_sobol_points(bounds: &[(f64, f64)], num_points: usize) -> Vec<Array1<f64>> {
994    // Simplified Sobol sequence (in practice, use proper implementation)
995    let dim = bounds.len();
996    let mut points = Vec::new();
997
998    for i in 0..num_points {
999        let mut point = Array1::zeros(dim);
1000        for (j, &(low, high)) in bounds.iter().enumerate() {
1001            // Van der Corput sequence for dimension j
1002            let mut n = i + 1;
1003            let base = 2 + j; // Different base for each dimension
1004            let mut result = 0.0;
1005            let mut f = 1.0 / base as f64;
1006
1007            while n > 0 {
1008                result += (n % base) as f64 * f;
1009                n /= base;
1010                f /= base as f64;
1011            }
1012
1013            point[j] = low + result * (high - low);
1014        }
1015        points.push(point);
1016    }
1017
1018    points
1019}
1020
1021#[cfg(test)]
1022mod tests {
1023    use super::*;
1024    use approx::assert_abs_diff_eq;
1025
1026    #[test]
1027    fn test_simple_clustering() {
1028        // Create some test minima
1029        let mut minima = vec![
1030            LocalMinimum {
1031                x: Array1::from_vec(vec![0.0, 0.0]),
1032                f: 1.0,
1033                fun_value: 1.0,
1034                nit: 10,
1035                func_evals: 20,
1036                success: true,
1037                start_point: Array1::from_vec(vec![1.0, 1.0]),
1038                cluster_id: None,
1039                cluster_distance: None,
1040            },
1041            LocalMinimum {
1042                x: Array1::from_vec(vec![0.1, 0.1]),
1043                f: 1.1,
1044                fun_value: 1.1,
1045                nit: 12,
1046                func_evals: 24,
1047                success: true,
1048                start_point: Array1::from_vec(vec![1.1, 1.1]),
1049                cluster_id: None,
1050                cluster_distance: None,
1051            },
1052            LocalMinimum {
1053                x: Array1::from_vec(vec![5.0, 5.0]),
1054                f: 2.0,
1055                fun_value: 2.0,
1056                nit: 15,
1057                func_evals: 30,
1058                success: true,
1059                start_point: Array1::from_vec(vec![5.5, 5.5]),
1060                cluster_id: None,
1061                cluster_distance: None,
1062            },
1063        ];
1064
1065        let options = ClusteringOptions {
1066            distance_threshold: 1.0,
1067            algorithm: ClusteringAlgorithm::Threshold,
1068            ..Default::default()
1069        };
1070
1071        threshold_clustering(&mut minima, &options).unwrap();
1072
1073        // First two should be in same cluster, third in different cluster
1074        assert_eq!(minima[0].cluster_id, minima[1].cluster_id);
1075        assert_ne!(minima[0].cluster_id, minima[2].cluster_id);
1076    }
1077
1078    #[test]
1079    fn test_distance_computation() {
1080        let min1 = LocalMinimum {
1081            x: Array1::from_vec(vec![0.0, 0.0]),
1082            f: 1.0,
1083            fun_value: 1.0,
1084            nit: 10,
1085            func_evals: 20,
1086            success: true,
1087            start_point: Array1::from_vec(vec![1.0, 1.0]),
1088            cluster_id: None,
1089            cluster_distance: None,
1090        };
1091
1092        let min2 = LocalMinimum {
1093            x: Array1::from_vec(vec![3.0, 4.0]),
1094            f: 2.0,
1095            fun_value: 2.0,
1096            nit: 12,
1097            func_evals: 24,
1098            success: true,
1099            start_point: Array1::from_vec(vec![3.5, 4.5]),
1100            cluster_id: None,
1101            cluster_distance: None,
1102        };
1103
1104        let options = ClusteringOptions::default();
1105        let distance = compute_distance(&min1, &min2, &options);
1106
1107        // Euclidean distance is 5.0, plus function weight * |1.0 - 2.0| = 5.0 + 0.1 * 1.0 = 5.1
1108        assert_abs_diff_eq!(distance, 5.1, epsilon = 1e-10);
1109    }
1110
1111    #[test]
1112    fn test_start_point_generation() {
1113        let bounds = vec![(0.0, 10.0), (-5.0, 5.0)];
1114        let num_points = 5;
1115
1116        let random_points =
1117            generate_diverse_start_points(&bounds, num_points, StartPointStrategy::Random);
1118        assert_eq!(random_points.len(), num_points);
1119
1120        for point in &random_points {
1121            assert!(point[0] >= 0.0 && point[0] <= 10.0);
1122            assert!(point[1] >= -5.0 && point[1] <= 5.0);
1123        }
1124
1125        let grid_points =
1126            generate_diverse_start_points(&bounds, num_points, StartPointStrategy::Grid);
1127        assert_eq!(grid_points.len(), num_points);
1128
1129        for point in &grid_points {
1130            assert!(point[0] >= 0.0 && point[0] <= 10.0);
1131            assert!(point[1] >= -5.0 && point[1] <= 5.0);
1132        }
1133    }
1134}