1use super::{utils, MultiObjectiveConfig, MultiObjectiveOptimizer};
19use crate::error::OptimizeError;
20use crate::multi_objective::crossover::{CrossoverOperator, SimulatedBinaryCrossover};
21use crate::multi_objective::mutation::{MutationOperator, PolynomialMutation};
22use crate::multi_objective::selection::{SelectionOperator, TournamentSelection};
23use crate::multi_objective::solutions::{MultiObjectiveResult, MultiObjectiveSolution, Population};
24use scirs2_core::ndarray::{Array1, ArrayView1};
25use scirs2_core::random::rngs::StdRng;
26use scirs2_core::random::{Rng, RngExt, SeedableRng};
27use std::cmp::Ordering;
28
29pub struct NSGAIII {
31 config: MultiObjectiveConfig,
32 n_objectives: usize,
33 n_variables: usize,
34 reference_points: Vec<Array1<f64>>,
36 population: Population,
37 generation: usize,
38 n_evaluations: usize,
39 rng: StdRng,
40 crossover: SimulatedBinaryCrossover,
41 mutation: PolynomialMutation,
42 selection: TournamentSelection,
43 ideal_point: Array1<f64>,
45 niche_count: Vec<usize>,
47 convergence_history: Vec<f64>,
48}
49
50impl NSGAIII {
51 pub fn new(population_size: usize, n_objectives: usize, n_variables: usize) -> Self {
53 let config = MultiObjectiveConfig {
54 population_size,
55 ..Default::default()
56 };
57 Self::with_config(config, n_objectives, n_variables, None)
58 }
59
60 pub fn with_config(
62 config: MultiObjectiveConfig,
63 n_objectives: usize,
64 n_variables: usize,
65 custom_reference_points: Option<Vec<Array1<f64>>>,
66 ) -> Self {
67 let seed = config.random_seed.unwrap_or_else(|| {
68 use std::time::{SystemTime, UNIX_EPOCH};
69 SystemTime::now()
70 .duration_since(UNIX_EPOCH)
71 .map(|d| d.as_secs())
72 .unwrap_or(42)
73 });
74
75 let rng = StdRng::seed_from_u64(seed);
76
77 let crossover =
78 SimulatedBinaryCrossover::new(config.crossover_eta, config.crossover_probability);
79 let mutation = PolynomialMutation::new(config.mutation_probability, config.mutation_eta);
80 let selection = TournamentSelection::new(2);
81
82 let reference_points = custom_reference_points.unwrap_or_else(|| {
84 let n_partitions = Self::auto_partitions(n_objectives, config.population_size);
85 utils::generate_das_dennis_points(n_objectives, n_partitions)
86 });
87
88 let niche_count = vec![0; reference_points.len()];
89
90 Self {
91 config,
92 n_objectives,
93 n_variables,
94 reference_points,
95 population: Population::new(),
96 generation: 0,
97 n_evaluations: 0,
98 rng,
99 crossover,
100 mutation,
101 selection,
102 ideal_point: Array1::from_elem(n_objectives, f64::INFINITY),
103 niche_count,
104 convergence_history: Vec::new(),
105 }
106 }
107
108 fn auto_partitions(n_objectives: usize, pop_size: usize) -> usize {
110 if n_objectives <= 2 {
113 return pop_size.max(4);
114 }
115 let mut p = 1;
117 loop {
118 let n_points = binomial_coefficient(p + n_objectives - 1, n_objectives - 1);
119 if n_points >= pop_size / 2 {
120 return p;
121 }
122 p += 1;
123 if p > 50 {
124 return p;
125 }
126 }
127 }
128
129 fn evaluate_individual<F>(
131 &mut self,
132 variables: &Array1<f64>,
133 objective_function: &F,
134 ) -> Result<Array1<f64>, OptimizeError>
135 where
136 F: Fn(&ArrayView1<f64>) -> Array1<f64> + Send + Sync,
137 {
138 self.n_evaluations += 1;
139
140 if let Some(max_evals) = self.config.max_evaluations {
141 if self.n_evaluations > max_evals {
142 return Err(OptimizeError::MaxEvaluationsReached);
143 }
144 }
145
146 let objectives = objective_function(&variables.view());
147 if objectives.len() != self.n_objectives {
148 return Err(OptimizeError::InvalidInput(format!(
149 "Expected {} objectives, got {}",
150 self.n_objectives,
151 objectives.len()
152 )));
153 }
154
155 Ok(objectives)
156 }
157
158 fn update_ideal_point(&mut self, solutions: &[MultiObjectiveSolution]) {
160 for sol in solutions {
161 for (i, &obj) in sol.objectives.iter().enumerate() {
162 if obj < self.ideal_point[i] {
163 self.ideal_point[i] = obj;
164 }
165 }
166 }
167 }
168
169 fn normalize_objectives(&self, solutions: &[MultiObjectiveSolution]) -> Vec<Array1<f64>> {
171 let n = solutions.len();
172 if n == 0 {
173 return vec![];
174 }
175
176 let translated: Vec<Array1<f64>> = solutions
178 .iter()
179 .map(|sol| &sol.objectives - &self.ideal_point)
180 .collect();
181
182 let mut extreme_points = Vec::with_capacity(self.n_objectives);
184 for j in 0..self.n_objectives {
185 let mut best_asf = f64::INFINITY;
186 let mut best_idx = 0;
187 for (i, t) in translated.iter().enumerate() {
188 let asf = (0..self.n_objectives)
190 .map(|k| {
191 if k == j {
192 t[k] } else {
194 t[k] * 1e6 }
196 })
197 .fold(f64::NEG_INFINITY, f64::max);
198 if asf < best_asf {
199 best_asf = asf;
200 best_idx = i;
201 }
202 }
203 extreme_points.push(translated[best_idx].clone());
204 }
205
206 let mut intercepts = Array1::from_elem(self.n_objectives, 1.0);
208 if extreme_points.len() == self.n_objectives {
209 for (j, ep) in extreme_points.iter().enumerate() {
212 let val = ep[j];
213 if val > 1e-10 {
214 intercepts[j] = val;
215 }
216 }
217 }
218
219 translated
221 .iter()
222 .map(|t| {
223 let mut normalized = Array1::zeros(self.n_objectives);
224 for j in 0..self.n_objectives {
225 normalized[j] = if intercepts[j] > 1e-10 {
226 t[j] / intercepts[j]
227 } else {
228 t[j]
229 };
230 }
231 normalized
232 })
233 .collect()
234 }
235
236 fn associate_to_reference_points(
239 &self,
240 normalized_objectives: &[Array1<f64>],
241 ) -> Vec<(usize, f64)> {
242 normalized_objectives
243 .iter()
244 .map(|obj| {
245 let mut min_dist = f64::INFINITY;
246 let mut min_idx = 0;
247
248 for (rp_idx, rp) in self.reference_points.iter().enumerate() {
249 let dist = perpendicular_distance(obj, rp);
250 if dist < min_dist {
251 min_dist = dist;
252 min_idx = rp_idx;
253 }
254 }
255
256 (min_idx, min_dist)
257 })
258 .collect()
259 }
260
261 fn niching_selection(
263 &mut self,
264 last_front_indices: &[usize],
265 all_solutions: &[MultiObjectiveSolution],
266 associations: &[(usize, f64)],
267 selected_so_far: &[usize],
268 remaining: usize,
269 ) -> Vec<usize> {
270 let mut niche_count = vec![0usize; self.reference_points.len()];
272 for &idx in selected_so_far {
273 let (rp_idx, _) = associations[idx];
274 niche_count[rp_idx] += 1;
275 }
276
277 let mut selected = Vec::with_capacity(remaining);
278 let mut available: Vec<usize> = last_front_indices.to_vec();
279
280 for _ in 0..remaining {
281 if available.is_empty() {
282 break;
283 }
284
285 let relevant_rps: Vec<usize> = available
288 .iter()
289 .map(|&idx| associations[idx].0)
290 .collect::<std::collections::HashSet<_>>()
291 .into_iter()
292 .collect();
293
294 if relevant_rps.is_empty() {
295 break;
296 }
297
298 let min_niche = relevant_rps
299 .iter()
300 .map(|&rp| niche_count[rp])
301 .min()
302 .unwrap_or(0);
303
304 let min_niche_rps: Vec<usize> = relevant_rps
306 .iter()
307 .filter(|&&rp| niche_count[rp] == min_niche)
308 .copied()
309 .collect();
310
311 let chosen_rp_idx = self.rng.random_range(0..min_niche_rps.len());
313 let chosen_rp = min_niche_rps[chosen_rp_idx];
314
315 let rp_members: Vec<usize> = available
317 .iter()
318 .filter(|&&idx| associations[idx].0 == chosen_rp)
319 .copied()
320 .collect();
321
322 if rp_members.is_empty() {
323 continue;
324 }
325
326 let chosen_member = if min_niche == 0 {
327 *rp_members
329 .iter()
330 .min_by(|&&a, &&b| {
331 associations[a]
332 .1
333 .partial_cmp(&associations[b].1)
334 .unwrap_or(Ordering::Equal)
335 })
336 .unwrap_or(&rp_members[0])
337 } else {
338 rp_members[self.rng.random_range(0..rp_members.len())]
340 };
341
342 selected.push(chosen_member);
343 niche_count[chosen_rp] += 1;
344 available.retain(|&idx| idx != chosen_member);
345 }
346
347 selected
348 }
349
350 fn create_offspring<F>(
352 &mut self,
353 objective_function: &F,
354 ) -> Result<Vec<MultiObjectiveSolution>, OptimizeError>
355 where
356 F: Fn(&ArrayView1<f64>) -> Array1<f64> + Send + Sync,
357 {
358 let mut offspring = Vec::new();
359 let pop_solutions = self.population.solutions().to_vec();
360
361 while offspring.len() < self.config.population_size {
362 let selected = self.selection.select(&pop_solutions, 2);
363 if selected.len() < 2 {
364 break;
365 }
366
367 let p1_vars = selected[0].variables.as_slice().unwrap_or(&[]);
368 let p2_vars = selected[1].variables.as_slice().unwrap_or(&[]);
369
370 let (mut c1_vars, mut c2_vars) = self.crossover.crossover(p1_vars, p2_vars);
371
372 let bounds: Vec<(f64, f64)> = if let Some((lower, upper)) = &self.config.bounds {
373 lower
374 .iter()
375 .zip(upper.iter())
376 .map(|(&l, &u)| (l, u))
377 .collect()
378 } else {
379 vec![(-1.0, 1.0); self.n_variables]
380 };
381
382 self.mutation.mutate(&mut c1_vars, &bounds);
383 self.mutation.mutate(&mut c2_vars, &bounds);
384
385 let c1_arr = Array1::from_vec(c1_vars);
386 let c1_obj = self.evaluate_individual(&c1_arr, objective_function)?;
387 offspring.push(MultiObjectiveSolution::new(c1_arr, c1_obj));
388
389 if offspring.len() < self.config.population_size {
390 let c2_arr = Array1::from_vec(c2_vars);
391 let c2_obj = self.evaluate_individual(&c2_arr, objective_function)?;
392 offspring.push(MultiObjectiveSolution::new(c2_arr, c2_obj));
393 }
394 }
395
396 Ok(offspring)
397 }
398
399 fn environmental_selection(
401 &mut self,
402 combined: Vec<MultiObjectiveSolution>,
403 ) -> Vec<MultiObjectiveSolution> {
404 let target_size = self.config.population_size;
405 if combined.len() <= target_size {
406 return combined;
407 }
408
409 let mut temp_pop = Population::from_solutions(combined.clone());
411 let fronts = temp_pop.non_dominated_sort();
412
413 let mut selected_indices: Vec<usize> = Vec::new();
415 let mut last_front_idx = 0;
416
417 for (fi, front) in fronts.iter().enumerate() {
418 if selected_indices.len() + front.len() <= target_size {
419 selected_indices.extend(front);
420 } else {
421 last_front_idx = fi;
422 break;
423 }
424 }
425
426 if selected_indices.len() >= target_size {
427 return selected_indices
428 .iter()
429 .take(target_size)
430 .map(|&i| combined[i].clone())
431 .collect();
432 }
433
434 let remaining = target_size - selected_indices.len();
436 let last_front = &fronts[last_front_idx];
437
438 self.update_ideal_point(&combined);
440
441 let all_considered: Vec<usize> = selected_indices
443 .iter()
444 .chain(last_front.iter())
445 .copied()
446 .collect();
447 let all_solutions: Vec<MultiObjectiveSolution> = all_considered
448 .iter()
449 .map(|&i| combined[i].clone())
450 .collect();
451
452 let normalized = self.normalize_objectives(&all_solutions);
453 let associations = self.associate_to_reference_points(&normalized);
454
455 let n_selected = selected_indices.len();
457 let last_front_local: Vec<usize> = (n_selected..all_solutions.len()).collect();
458 let selected_local: Vec<usize> = (0..n_selected).collect();
459
460 let niching_result = self.niching_selection(
461 &last_front_local,
462 &all_solutions,
463 &associations,
464 &selected_local,
465 remaining,
466 );
467
468 let mut result: Vec<MultiObjectiveSolution> = selected_indices
470 .iter()
471 .map(|&i| combined[i].clone())
472 .collect();
473
474 for local_idx in niching_result {
475 let global_idx = all_considered[local_idx];
476 result.push(combined[global_idx].clone());
477 }
478
479 result
480 }
481
482 fn calculate_metrics(&mut self) {
484 if let Some(ref_point) = &self.config.reference_point {
485 let pareto_front = self.population.extract_pareto_front();
486 let hv = utils::calculate_hypervolume(&pareto_front, ref_point);
487 self.convergence_history.push(hv);
488 }
489 }
490}
491
492impl MultiObjectiveOptimizer for NSGAIII {
493 fn optimize<F>(&mut self, objective_function: F) -> Result<MultiObjectiveResult, OptimizeError>
494 where
495 F: Fn(&ArrayView1<f64>) -> Array1<f64> + Send + Sync,
496 {
497 self.initialize_population()?;
498
499 let initial_vars = utils::generate_random_population(
501 self.config.population_size,
502 self.n_variables,
503 &self.config.bounds,
504 );
505
506 let mut initial_solutions = Vec::new();
507 for vars in initial_vars {
508 let objs = self.evaluate_individual(&vars, &objective_function)?;
509 initial_solutions.push(MultiObjectiveSolution::new(vars, objs));
510 }
511
512 self.update_ideal_point(&initial_solutions);
513 self.population = Population::from_solutions(initial_solutions);
514
515 while self.generation < self.config.max_generations {
517 if self.check_convergence() {
518 break;
519 }
520 self.evolve_generation(&objective_function)?;
521 }
522
523 let pareto_front = self.population.extract_pareto_front();
525 let hypervolume = self
526 .config
527 .reference_point
528 .as_ref()
529 .map(|rp| utils::calculate_hypervolume(&pareto_front, rp));
530
531 let mut result = MultiObjectiveResult::new(
532 pareto_front,
533 self.population.solutions().to_vec(),
534 self.n_evaluations,
535 self.generation,
536 );
537 result.hypervolume = hypervolume;
538 result.metrics.convergence_history = self.convergence_history.clone();
539 result.metrics.population_stats = self.population.calculate_statistics();
540
541 Ok(result)
542 }
543
544 fn evolve_generation<F>(&mut self, objective_function: &F) -> Result<(), OptimizeError>
545 where
546 F: Fn(&ArrayView1<f64>) -> Array1<f64> + Send + Sync,
547 {
548 let offspring = self.create_offspring(objective_function)?;
549
550 let mut combined = self.population.solutions().to_vec();
551 combined.extend(offspring);
552
553 let next_pop = self.environmental_selection(combined);
554 self.population = Population::from_solutions(next_pop);
555
556 self.generation += 1;
557 self.calculate_metrics();
558
559 Ok(())
560 }
561
562 fn initialize_population(&mut self) -> Result<(), OptimizeError> {
563 self.population.clear();
564 self.generation = 0;
565 self.n_evaluations = 0;
566 self.ideal_point = Array1::from_elem(self.n_objectives, f64::INFINITY);
567 self.niche_count = vec![0; self.reference_points.len()];
568 self.convergence_history.clear();
569 Ok(())
570 }
571
572 fn check_convergence(&self) -> bool {
573 if let Some(max_evals) = self.config.max_evaluations {
574 if self.n_evaluations >= max_evals {
575 return true;
576 }
577 }
578
579 if self.convergence_history.len() >= 10 {
580 let recent = &self.convergence_history[self.convergence_history.len() - 10..];
581 let max_hv = recent.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
582 let min_hv = recent.iter().fold(f64::INFINITY, |a, &b| a.min(b));
583 if (max_hv - min_hv) < self.config.tolerance {
584 return true;
585 }
586 }
587
588 false
589 }
590
591 fn get_population(&self) -> &Population {
592 &self.population
593 }
594
595 fn get_generation(&self) -> usize {
596 self.generation
597 }
598
599 fn get_evaluations(&self) -> usize {
600 self.n_evaluations
601 }
602
603 fn name(&self) -> &str {
604 "NSGA-III"
605 }
606}
607
608fn perpendicular_distance(point: &Array1<f64>, direction: &Array1<f64>) -> f64 {
610 let dir_norm_sq = direction.dot(direction);
611 if dir_norm_sq < 1e-30 {
612 return point.dot(point).sqrt();
613 }
614
615 let proj_scalar = point.dot(direction) / dir_norm_sq;
617 let proj = proj_scalar * direction;
618
619 let diff = point - &proj;
621 diff.dot(&diff).sqrt()
622}
623
624fn binomial_coefficient(n: usize, k: usize) -> usize {
626 if k > n {
627 return 0;
628 }
629 if k == 0 || k == n {
630 return 1;
631 }
632 let k = k.min(n - k);
633 let mut result = 1usize;
634 for i in 0..k {
635 result = result.saturating_mul(n - i);
636 result /= i + 1;
637 }
638 result
639}
640
641#[cfg(test)]
642mod tests {
643 use super::*;
644 use scirs2_core::ndarray::{array, s};
645
646 fn zdt1(x: &ArrayView1<f64>) -> Array1<f64> {
647 let f1 = x[0];
648 let g = 1.0 + 9.0 * x.slice(s![1..]).sum() / (x.len() - 1) as f64;
649 let f2 = g * (1.0 - (f1 / g).sqrt());
650 array![f1, f2]
651 }
652
653 fn dtlz1(x: &ArrayView1<f64>) -> Array1<f64> {
654 let n = x.len();
656 let k = n - 2; let g: f64 = (0..k)
658 .map(|i| {
659 let xi = x[n - k + i];
660 (xi - 0.5).powi(2) - (20.0 * std::f64::consts::PI * (xi - 0.5)).cos()
661 })
662 .sum::<f64>();
663 let g = 100.0 * (k as f64 + g);
664
665 let f1 = 0.5 * x[0] * x[1] * (1.0 + g);
666 let f2 = 0.5 * x[0] * (1.0 - x[1]) * (1.0 + g);
667 let f3 = 0.5 * (1.0 - x[0]) * (1.0 + g);
668 array![f1, f2, f3]
669 }
670
671 #[test]
672 fn test_nsga3_creation() {
673 let nsga3 = NSGAIII::new(100, 3, 5);
674 assert_eq!(nsga3.n_objectives, 3);
675 assert_eq!(nsga3.n_variables, 5);
676 assert!(!nsga3.reference_points.is_empty());
677 assert_eq!(nsga3.generation, 0);
678 }
679
680 #[test]
681 fn test_nsga3_with_config() {
682 let mut config = MultiObjectiveConfig::default();
683 config.population_size = 50;
684 config.max_generations = 10;
685 config.random_seed = Some(42);
686
687 let nsga3 = NSGAIII::with_config(config, 2, 3, None);
688 assert_eq!(nsga3.n_objectives, 2);
689 assert!(!nsga3.reference_points.is_empty());
690 }
691
692 #[test]
693 fn test_nsga3_custom_reference_points() {
694 let rps = vec![
695 array![1.0, 0.0, 0.0],
696 array![0.0, 1.0, 0.0],
697 array![0.0, 0.0, 1.0],
698 array![0.5, 0.5, 0.0],
699 array![0.5, 0.0, 0.5],
700 array![0.0, 0.5, 0.5],
701 array![0.333, 0.333, 0.334],
702 ];
703
704 let config = MultiObjectiveConfig {
705 population_size: 20,
706 max_generations: 5,
707 random_seed: Some(42),
708 ..Default::default()
709 };
710
711 let nsga3 = NSGAIII::with_config(config, 3, 5, Some(rps.clone()));
712 assert_eq!(nsga3.reference_points.len(), 7);
713 }
714
715 #[test]
716 fn test_nsga3_optimize_zdt1() {
717 let mut config = MultiObjectiveConfig::default();
718 config.max_generations = 10;
719 config.population_size = 20;
720 config.bounds = Some((Array1::zeros(3), Array1::ones(3)));
721 config.random_seed = Some(42);
722
723 let mut nsga3 = NSGAIII::with_config(config, 2, 3, None);
724 let result = nsga3.optimize(zdt1);
725
726 assert!(result.is_ok());
727 let res = result.expect("should succeed");
728 assert!(res.success);
729 assert!(!res.pareto_front.is_empty());
730 assert!(res.n_evaluations > 0);
731 }
732
733 #[test]
734 fn test_nsga3_optimize_dtlz1() {
735 let mut config = MultiObjectiveConfig::default();
736 config.max_generations = 10;
737 config.population_size = 20;
738 config.bounds = Some((Array1::zeros(5), Array1::ones(5)));
739 config.random_seed = Some(42);
740
741 let mut nsga3 = NSGAIII::with_config(config, 3, 5, None);
742 let result = nsga3.optimize(dtlz1);
743
744 assert!(result.is_ok());
745 let res = result.expect("should succeed");
746 assert!(res.success);
747 assert!(!res.pareto_front.is_empty());
748 }
749
750 #[test]
751 fn test_nsga3_max_evaluations() {
752 let mut config = MultiObjectiveConfig::default();
753 config.max_generations = 1000;
754 config.max_evaluations = Some(50);
755 config.population_size = 10;
756 config.bounds = Some((Array1::zeros(3), Array1::ones(3)));
757 config.random_seed = Some(42);
758
759 let mut nsga3 = NSGAIII::with_config(config, 2, 3, None);
760 let result = nsga3.optimize(zdt1);
761 assert!(result.is_ok());
762 let res = result.expect("should succeed");
763 assert!(res.n_evaluations <= 60); }
765
766 #[test]
767 fn test_perpendicular_distance() {
768 let point = array![1.0, 1.0];
769 let direction = array![1.0, 0.0];
770
771 let dist = perpendicular_distance(&point, &direction);
772 assert!(
773 (dist - 1.0).abs() < 1e-10,
774 "Distance should be 1.0, got {}",
775 dist
776 );
777
778 let point2 = array![2.0, 0.0];
780 let dist2 = perpendicular_distance(&point2, &direction);
781 assert!(dist2.abs() < 1e-10, "Distance should be 0.0, got {}", dist2);
782 }
783
784 #[test]
785 fn test_binomial_coefficient() {
786 assert_eq!(binomial_coefficient(5, 2), 10);
787 assert_eq!(binomial_coefficient(10, 3), 120);
788 assert_eq!(binomial_coefficient(4, 0), 1);
789 assert_eq!(binomial_coefficient(4, 4), 1);
790 assert_eq!(binomial_coefficient(0, 0), 1);
791 }
792
793 #[test]
794 fn test_nsga3_name() {
795 let nsga3 = NSGAIII::new(50, 2, 3);
796 assert_eq!(nsga3.name(), "NSGA-III");
797 }
798
799 #[test]
800 fn test_nsga3_convergence_check() {
801 let config = MultiObjectiveConfig {
802 tolerance: 1e-10,
803 max_generations: 2,
804 population_size: 10,
805 ..Default::default()
806 };
807 let nsga3 = NSGAIII::with_config(config, 2, 2, None);
808 assert!(!nsga3.check_convergence());
809 }
810}