1use crate::error::OptimizeError;
8use crate::unconstrained::{minimize, Method, Options};
9use scirs2_core::ndarray::{Array1, Array2, ArrayView1};
10use std::collections::HashMap;
11
12#[derive(Debug, Clone)]
14pub struct ClusteringOptions {
15 pub distance_threshold: f64,
17 pub function_tolerance: f64,
19 pub max_clusters: Option<usize>,
21 pub algorithm: ClusteringAlgorithm,
23 pub normalize_coordinates: bool,
25 pub use_function_values: bool,
27 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#[derive(Debug, Clone, Copy)]
47pub enum ClusteringAlgorithm {
48 Hierarchical,
50 KMeans,
52 Density,
54 Threshold,
56}
57
58#[derive(Debug, Clone)]
60pub struct LocalMinimum<S> {
61 pub x: Array1<f64>,
63 pub f: f64,
65 pub fun_value: S,
67 pub nit: usize,
69 pub func_evals: usize,
71 pub success: bool,
73 pub start_point: Array1<f64>,
75 pub cluster_id: Option<usize>,
77 pub cluster_distance: Option<f64>,
79}
80
81#[derive(Debug, Clone)]
83pub struct ClusteringResult<S> {
84 pub minima: Vec<LocalMinimum<S>>,
86 pub centroids: Vec<ClusterCentroid>,
88 pub num_clusters: usize,
90 pub silhouette_score: Option<f64>,
92 pub wcss: f64,
94 pub best_minimum: Option<LocalMinimum<S>>,
96}
97
98#[derive(Debug, Clone)]
100pub struct ClusterCentroid {
101 pub x: Array1<f64>,
103 pub f_avg: f64,
105 pub f_min: f64,
107 pub size: usize,
109 pub radius: f64,
111}
112
113#[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 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 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 cluster_minima(&mut minima, &clustering_opts)?;
168
169 let centroids = compute_cluster_centroids(&minima)?;
171 let num_clusters = centroids.len();
172 let wcss = compute_wcss(&minima, ¢roids);
173 let silhouette_score = compute_silhouette_score(&minima);
174
175 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#[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#[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 let distances = compute_distance_matrix(minima, options);
232
233 let mut cluster_assignments = vec![None; n];
235 let mut next_cluster_id = n;
236
237 for (i, assignment) in cluster_assignments.iter_mut().enumerate().take(n) {
239 *assignment = Some(i);
240 }
241
242 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 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 if min_dist > options.distance_threshold {
263 break;
264 }
265
266 if let Some(max_clusters) = options.max_clusters {
268 if active_clusters.len() <= max_clusters {
269 break;
270 }
271 }
272
273 let (i, j) = merge_pair;
275 let mut merged_cluster = active_clusters[i].clone();
276 merged_cluster.extend(&active_clusters[j]);
277
278 for &point_idx in &merged_cluster {
280 cluster_assignments[point_idx] = Some(next_cluster_id);
281 }
282
283 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 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 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#[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 let k = if let Some(max_k) = options.max_clusters {
335 std::cmp::min(max_k, n)
336 } else {
337 std::cmp::min((n as f64).sqrt() as usize + 1, n)
339 };
340
341 if k >= n {
342 for (i, minimum) in minima.iter_mut().enumerate() {
344 minimum.cluster_id = Some(i);
345 }
346 return Ok(());
347 }
348
349 let features = extract_features(minima, options);
351 let dim = features.ncols();
352
353 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 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), ¢roids.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 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 let centroid_change = (¢roids - &new_centroids).mapv(|x| x.abs()).sum();
407
408 centroids = new_centroids;
409
410 if centroid_change < tolerance {
411 break;
412 }
413 }
414
415 for (i, minimum) in minima.iter_mut().enumerate() {
417 minimum.cluster_id = Some(assignments[i]);
418 }
419
420 Ok(())
421}
422
423#[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; 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 let neighbors: Vec<usize> = (0..n).filter(|&j| distances[[i, j]] <= eps).collect();
456
457 if neighbors.len() < min_pts {
458 continue;
460 }
461
462 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 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 for (i, minimum) in minima.iter_mut().enumerate() {
498 minimum.cluster_id = cluster_assignments[i];
499 }
500
501 Ok(())
502}
503
504#[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 cluster_assignments[i] = Some(cluster_id);
531
532 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 for (i, minimum) in minima.iter_mut().enumerate() {
549 minimum.cluster_id = cluster_assignments[i];
550 }
551
552 Ok(())
553}
554
555#[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#[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 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 let func_dist = (min1.f - min2.f).abs();
597
598 coord_dist + options.function_weight * func_dist
600}
601
602#[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#[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 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 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 if options.normalize_coordinates {
656 normalize_features(&mut features, coord_dim);
657 }
658
659 features
660}
661
662#[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 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#[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 let first_idx = 0; centroids.row_mut(0).assign(&features.row(first_idx));
697
698 for c in 1..k {
700 let mut distances = vec![f64::INFINITY; n];
701
702 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, ¢roid);
708 *distance = distance.min(dist);
709 }
710 }
711
712 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#[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#[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 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 for minimum in &cluster_minima {
767 centroid_x = ¢roid_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 let mut max_radius = 0.0;
778 for minimum in &cluster_minima {
779 let dist = (&minimum.x - ¢roid_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#[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 = ¢roids[cluster_id];
807 let dist = (&minimum.x - ¢roid.x).mapv(|x| x.powi(2)).sum();
808 wcss += dist;
809 }
810 }
811 }
812
813 wcss
814}
815
816#[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 let mut intra_sum = 0.0;
833 let mut intra_count = 0;
834
835 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 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#[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#[derive(Debug, Clone, Copy)]
901pub enum StartPointStrategy {
902 Random,
903 LatinHypercube,
904 Grid,
905 Sobol,
906}
907
908#[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 let t = (i * num_points + points.len()) as f64 / (num_points * dim) as f64;
919 let random_val = (t * 17.0).fract(); point[i] = low + random_val * (high - low);
921 }
922 points.push(point);
923 }
924
925 points
926}
927
928#[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 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; point[j] = low + segment * (high - low);
940 }
941 points.push(point);
942 }
943
944 points
945}
946
947#[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 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 _points.truncate(num_points);
988 _points
989}
990
991#[allow(dead_code)]
993fn generate_sobol_points(bounds: &[(f64, f64)], num_points: usize) -> Vec<Array1<f64>> {
994 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 let mut n = i + 1;
1003 let base = 2 + j; 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 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 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 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}