scirs2_optimize/multi_objective/solutions/
mod.rs

1//! Multi-objective solution representations and utilities
2//!
3//! This module provides core data structures for multi-objective optimization:
4//! - Solution representation with objectives and constraints
5//! - Population management utilities
6//! - Pareto front analysis
7//! - Solution comparison and ranking
8
9use scirs2_core::ndarray::{Array1, ArrayView1};
10use serde::{Deserialize, Serialize};
11use std::cmp::Ordering;
12use std::collections::HashMap;
13
14/// Type alias for backward compatibility
15pub type Solution = MultiObjectiveSolution;
16
17/// Represents a solution in multi-objective optimization
18#[derive(Debug, Clone, Serialize, Deserialize)]
19pub struct MultiObjectiveSolution {
20    /// Decision variables
21    pub variables: Array1<f64>,
22    /// Objective function values
23    pub objectives: Array1<f64>,
24    /// Constraint violation (0.0 if feasible)
25    pub constraint_violation: f64,
26    /// Rank in the population (for NSGA-II/III)
27    pub rank: usize,
28    /// Crowding distance (for NSGA-II)
29    pub crowding_distance: f64,
30    /// Additional metadata
31    pub metadata: HashMap<String, f64>,
32}
33
34impl MultiObjectiveSolution {
35    /// Create a new solution
36    pub fn new(variables: Array1<f64>, objectives: Array1<f64>) -> Self {
37        Self {
38            variables,
39            objectives,
40            constraint_violation: 0.0,
41            rank: 0,
42            crowding_distance: 0.0,
43            metadata: HashMap::new(),
44        }
45    }
46
47    /// Create a new solution with constraint violation
48    pub fn new_with_constraints(
49        variables: Array1<f64>,
50        objectives: Array1<f64>,
51        constraint_violation: f64,
52    ) -> Self {
53        Self {
54            variables,
55            objectives,
56            constraint_violation,
57            rank: 0,
58            crowding_distance: 0.0,
59            metadata: HashMap::new(),
60        }
61    }
62
63    /// Check if this solution dominates another
64    pub fn dominates(&self, other: &Self) -> bool {
65        // First check constraint violations
66        if self.constraint_violation < other.constraint_violation {
67            return true;
68        }
69        if self.constraint_violation > other.constraint_violation {
70            return false;
71        }
72
73        // If both are equally feasible/infeasible, compare objectives
74        let mut at_least_one_better = false;
75
76        for (obj1, obj2) in self.objectives.iter().zip(other.objectives.iter()) {
77            if obj1 > obj2 {
78                return false; // Assuming minimization
79            }
80            if obj1 < obj2 {
81                at_least_one_better = true;
82            }
83        }
84
85        at_least_one_better
86    }
87
88    /// Check if this solution is dominated by another
89    pub fn is_dominated_by(&self, other: &Self) -> bool {
90        other.dominates(self)
91    }
92
93    /// Check if solutions are non-dominated with respect to each other
94    pub fn is_non_dominated_with(&self, other: &Self) -> bool {
95        !self.dominates(other) && !other.dominates(self)
96    }
97
98    /// Calculate distance to another solution in objective space
99    pub fn objective_distance(&self, other: &Self) -> f64 {
100        self.objectives
101            .iter()
102            .zip(other.objectives.iter())
103            .map(|(a, b)| (a - b).powi(2))
104            .sum::<f64>()
105            .sqrt()
106    }
107
108    /// Calculate distance to another solution in variable space
109    pub fn variable_distance(&self, other: &Self) -> f64 {
110        self.variables
111            .iter()
112            .zip(other.variables.iter())
113            .map(|(a, b)| (a - b).powi(2))
114            .sum::<f64>()
115            .sqrt()
116    }
117
118    /// Get number of objectives
119    pub fn n_objectives(&self) -> usize {
120        self.objectives.len()
121    }
122
123    /// Get number of variables
124    pub fn n_variables(&self) -> usize {
125        self.variables.len()
126    }
127
128    /// Check if solution is feasible
129    pub fn is_feasible(&self) -> bool {
130        self.constraint_violation <= 1e-10
131    }
132
133    /// Set metadata value
134    pub fn set_metadata(&mut self, key: String, value: f64) {
135        self.metadata.insert(key, value);
136    }
137
138    /// Get metadata value
139    pub fn get_metadata(&self, key: &str) -> Option<f64> {
140        self.metadata.get(key).copied()
141    }
142
143    /// Clone with new objectives (useful for constraint handling)
144    pub fn with_objectives(&self, objectives: Array1<f64>) -> Self {
145        Self {
146            variables: self.variables.clone(),
147            objectives,
148            constraint_violation: self.constraint_violation,
149            rank: self.rank,
150            crowding_distance: self.crowding_distance,
151            metadata: self.metadata.clone(),
152        }
153    }
154
155    /// Clone with new variables
156    pub fn with_variables(&self, variables: Array1<f64>) -> Self {
157        Self {
158            variables,
159            objectives: self.objectives.clone(),
160            constraint_violation: self.constraint_violation,
161            rank: self.rank,
162            crowding_distance: self.crowding_distance,
163            metadata: self.metadata.clone(),
164        }
165    }
166}
167
168impl PartialEq for MultiObjectiveSolution {
169    fn eq(&self, other: &Self) -> bool {
170        self.variables == other.variables && self.objectives == other.objectives
171    }
172}
173
174impl PartialOrd for MultiObjectiveSolution {
175    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
176        // Compare by rank first, then by crowding distance
177        match self.rank.cmp(&other.rank) {
178            Ordering::Equal => {
179                // Higher crowding distance is better (reversed comparison)
180                other.crowding_distance.partial_cmp(&self.crowding_distance)
181            }
182            other_order => Some(other_order),
183        }
184    }
185}
186
187/// Result of multi-objective optimization
188#[derive(Debug, Clone, Serialize, Deserialize)]
189pub struct MultiObjectiveResult {
190    /// Pareto front solutions
191    pub pareto_front: Vec<MultiObjectiveSolution>,
192    /// All final population solutions
193    pub population: Vec<MultiObjectiveSolution>,
194    /// Number of function evaluations
195    pub n_evaluations: usize,
196    /// Number of generations/iterations
197    pub n_generations: usize,
198    /// Success flag
199    pub success: bool,
200    /// Convergence message
201    pub message: String,
202    /// Hypervolume indicator (if reference point provided)
203    pub hypervolume: Option<f64>,
204    /// Additional metrics
205    pub metrics: OptimizationMetrics,
206}
207
208impl MultiObjectiveResult {
209    /// Create a new result
210    pub fn new(
211        pareto_front: Vec<MultiObjectiveSolution>,
212        population: Vec<MultiObjectiveSolution>,
213        n_evaluations: usize,
214        n_generations: usize,
215    ) -> Self {
216        Self {
217            pareto_front,
218            population,
219            n_evaluations,
220            n_generations,
221            success: true,
222            message: "Optimization completed successfully".to_string(),
223            hypervolume: None,
224            metrics: OptimizationMetrics::default(),
225        }
226    }
227
228    /// Create a failed result
229    pub fn failed(message: String, n_evaluations: usize, n_generations: usize) -> Self {
230        Self {
231            pareto_front: vec![],
232            population: vec![],
233            n_evaluations,
234            n_generations,
235            success: false,
236            message,
237            hypervolume: None,
238            metrics: OptimizationMetrics::default(),
239        }
240    }
241
242    /// Get best solution for a specific objective
243    pub fn best_for_objective(&self, objective_index: usize) -> Option<&MultiObjectiveSolution> {
244        self.pareto_front.iter().min_by(|a, b| {
245            a.objectives[objective_index]
246                .partial_cmp(&b.objectives[objective_index])
247                .unwrap_or(Ordering::Equal)
248        })
249    }
250
251    /// Get solution closest to ideal point
252    pub fn closest_to_ideal(&self, ideal_point: &Array1<f64>) -> Option<&MultiObjectiveSolution> {
253        self.pareto_front.iter().min_by(|a, b| {
254            let dist_a = a
255                .objectives
256                .iter()
257                .zip(ideal_point.iter())
258                .map(|(obj, ideal)| (obj - ideal).powi(2))
259                .sum::<f64>();
260            let dist_b = b
261                .objectives
262                .iter()
263                .zip(ideal_point.iter())
264                .map(|(obj, ideal)| (obj - ideal).powi(2))
265                .sum::<f64>();
266            dist_a.partial_cmp(&dist_b).unwrap_or(Ordering::Equal)
267        })
268    }
269
270    /// Get number of solutions in Pareto front
271    pub fn pareto_front_size(&self) -> usize {
272        self.pareto_front.len()
273    }
274
275    /// Check if result contains feasible solutions
276    pub fn has_feasible_solutions(&self) -> bool {
277        self.pareto_front.iter().any(|sol| sol.is_feasible())
278    }
279
280    /// Get all feasible solutions
281    pub fn feasible_solutions(&self) -> Vec<&MultiObjectiveSolution> {
282        self.pareto_front
283            .iter()
284            .filter(|sol| sol.is_feasible())
285            .collect()
286    }
287}
288
289/// Optimization metrics for tracking performance
290#[derive(Debug, Clone, Serialize, Deserialize, Default)]
291pub struct OptimizationMetrics {
292    /// Convergence history (hypervolume over generations)
293    pub convergence_history: Vec<f64>,
294    /// Diversity metrics over generations
295    pub diversity_history: Vec<f64>,
296    /// Average objective values over generations
297    pub average_objectives: Vec<Array1<f64>>,
298    /// Best objective values over generations
299    pub best_objectives: Vec<Array1<f64>>,
300    /// Population statistics
301    pub population_stats: PopulationStatistics,
302}
303
304/// Population statistics
305#[derive(Debug, Clone, Serialize, Deserialize)]
306pub struct PopulationStatistics {
307    /// Mean objective values
308    pub mean_objectives: Array1<f64>,
309    /// Standard deviation of objectives
310    pub std_objectives: Array1<f64>,
311    /// Minimum objective values
312    pub min_objectives: Array1<f64>,
313    /// Maximum objective values
314    pub max_objectives: Array1<f64>,
315    /// Feasibility ratio
316    pub feasibility_ratio: f64,
317    /// Average constraint violation
318    pub avg_constraint_violation: f64,
319}
320
321impl Default for PopulationStatistics {
322    fn default() -> Self {
323        Self {
324            mean_objectives: Array1::zeros(0),
325            std_objectives: Array1::zeros(0),
326            min_objectives: Array1::zeros(0),
327            max_objectives: Array1::zeros(0),
328            feasibility_ratio: 0.0,
329            avg_constraint_violation: 0.0,
330        }
331    }
332}
333
334/// Population management utilities
335#[derive(Debug, Clone)]
336pub struct Population {
337    solutions: Vec<MultiObjectiveSolution>,
338}
339
340impl Population {
341    /// Create a new empty population
342    pub fn new() -> Self {
343        Self { solutions: vec![] }
344    }
345
346    /// Create a new population with specified capacity
347    pub fn with_capacity(
348        population_size: usize,
349        _n_objectives: usize,
350        _n_variables: usize,
351    ) -> Self {
352        Self {
353            solutions: Vec::with_capacity(population_size),
354        }
355    }
356
357    /// Create population from solutions
358    pub fn from_solutions(solutions: Vec<MultiObjectiveSolution>) -> Self {
359        Self { solutions }
360    }
361
362    /// Add a solution to the population
363    pub fn add(&mut self, solution: MultiObjectiveSolution) {
364        self.solutions.push(solution);
365    }
366
367    /// Get all solutions
368    pub fn solutions(&self) -> &[MultiObjectiveSolution] {
369        &self.solutions
370    }
371
372    /// Get mutable reference to solutions
373    pub fn solutions_mut(&mut self) -> &mut Vec<MultiObjectiveSolution> {
374        &mut self.solutions
375    }
376
377    /// Get population size
378    pub fn size(&self) -> usize {
379        self.solutions.len()
380    }
381
382    /// Check if population is empty
383    pub fn is_empty(&self) -> bool {
384        self.solutions.is_empty()
385    }
386
387    /// Clear the population
388    pub fn clear(&mut self) {
389        self.solutions.clear();
390    }
391
392    /// Extract Pareto front from population
393    pub fn extract_pareto_front(&self) -> Vec<MultiObjectiveSolution> {
394        let mut pareto_front: Vec<MultiObjectiveSolution> = Vec::new();
395
396        for candidate in &self.solutions {
397            let mut is_dominated = false;
398
399            // Check if candidate is dominated by any existing solution in Pareto front
400            for existing in &pareto_front {
401                if existing.dominates(candidate) {
402                    is_dominated = true;
403                    break;
404                }
405            }
406
407            if !is_dominated {
408                // Remove solutions from Pareto front that are dominated by candidate
409                pareto_front.retain(|existing| !candidate.dominates(existing));
410                pareto_front.push(candidate.clone());
411            }
412        }
413
414        pareto_front
415    }
416
417    /// Perform non-dominated sorting (for NSGA-II)
418    pub fn non_dominated_sort(&mut self) -> Vec<Vec<usize>> {
419        let n = self.solutions.len();
420        let mut fronts = Vec::new();
421        let mut domination_counts = vec![0; n];
422        let mut dominated_solutions: Vec<Vec<usize>> = vec![Vec::new(); n];
423
424        // Calculate domination relationships
425        for i in 0..n {
426            for j in 0..n {
427                if i != j {
428                    if self.solutions[i].dominates(&self.solutions[j]) {
429                        dominated_solutions[i].push(j);
430                    } else if self.solutions[j].dominates(&self.solutions[i]) {
431                        domination_counts[i] += 1;
432                    }
433                }
434            }
435        }
436
437        // Find first front
438        let mut current_front = Vec::new();
439        for i in 0..n {
440            if domination_counts[i] == 0 {
441                self.solutions[i].rank = 0;
442                current_front.push(i);
443            }
444        }
445
446        let mut front_number = 0;
447
448        while !current_front.is_empty() {
449            fronts.push(current_front.clone());
450            let mut next_front = Vec::new();
451
452            for &i in &current_front {
453                for &j in &dominated_solutions[i] {
454                    domination_counts[j] -= 1;
455                    if domination_counts[j] == 0 {
456                        self.solutions[j].rank = front_number + 1;
457                        next_front.push(j);
458                    }
459                }
460            }
461
462            current_front = next_front;
463            front_number += 1;
464        }
465
466        fronts
467    }
468
469    /// Calculate crowding distances for solutions in the same front
470    pub fn calculate_crowding_distances(&mut self, front_indices: &[usize]) {
471        let front_size = front_indices.len();
472        if front_size <= 2 {
473            // Set infinite distance for boundary solutions
474            for &i in front_indices {
475                self.solutions[i].crowding_distance = f64::INFINITY;
476            }
477            return;
478        }
479
480        let n_objectives = self.solutions[front_indices[0]].n_objectives();
481
482        // Initialize distances to zero
483        for &i in front_indices {
484            self.solutions[i].crowding_distance = 0.0;
485        }
486
487        // Calculate distance for each objective
488        for obj in 0..n_objectives {
489            // Sort by current objective
490            let mut sorted_indices = front_indices.to_vec();
491            sorted_indices.sort_by(|&a, &b| {
492                self.solutions[a].objectives[obj]
493                    .partial_cmp(&self.solutions[b].objectives[obj])
494                    .unwrap_or(Ordering::Equal)
495            });
496
497            // Set boundary solutions to infinite distance
498            self.solutions[sorted_indices[0]].crowding_distance = f64::INFINITY;
499            self.solutions[sorted_indices[front_size - 1]].crowding_distance = f64::INFINITY;
500
501            // Calculate objective range
502            let obj_min = self.solutions[sorted_indices[0]].objectives[obj];
503            let obj_max = self.solutions[sorted_indices[front_size - 1]].objectives[obj];
504            let obj_range = obj_max - obj_min;
505
506            if obj_range > 0.0 {
507                // Calculate distances for intermediate solutions
508                for i in 1..front_size - 1 {
509                    if self.solutions[sorted_indices[i]].crowding_distance != f64::INFINITY {
510                        let distance = (self.solutions[sorted_indices[i + 1]].objectives[obj]
511                            - self.solutions[sorted_indices[i - 1]].objectives[obj])
512                            / obj_range;
513                        self.solutions[sorted_indices[i]].crowding_distance += distance;
514                    }
515                }
516            }
517        }
518    }
519
520    /// Select best solutions using NSGA-II crowded comparison
521    pub fn select_best(&mut self, target_size: usize) -> Vec<MultiObjectiveSolution> {
522        if self.solutions.len() <= target_size {
523            return self.solutions.clone();
524        }
525
526        // Perform non-dominated sorting
527        let fronts = self.non_dominated_sort();
528
529        // Calculate crowding distances for each front
530        for front in &fronts {
531            self.calculate_crowding_distances(front);
532        }
533
534        let mut selected = Vec::new();
535
536        // Add complete fronts until we reach capacity
537        for front in &fronts {
538            if selected.len() + front.len() <= target_size {
539                for &i in front {
540                    selected.push(self.solutions[i].clone());
541                }
542            } else {
543                // Add partial front based on crowding distance
544                let remaining = target_size - selected.len();
545                let mut front_solutions: Vec<_> =
546                    front.iter().map(|&i| &self.solutions[i]).collect();
547
548                // Sort by crowding distance (descending)
549                front_solutions.sort_by(|a, b| {
550                    b.crowding_distance
551                        .partial_cmp(&a.crowding_distance)
552                        .unwrap_or(Ordering::Equal)
553                });
554
555                for solution in front_solutions.iter().take(remaining) {
556                    selected.push((*solution).clone());
557                }
558                break;
559            }
560        }
561
562        selected
563    }
564
565    /// Calculate population statistics
566    pub fn calculate_statistics(&self) -> PopulationStatistics {
567        if self.solutions.is_empty() {
568            return PopulationStatistics::default();
569        }
570
571        let n_objectives = self.solutions[0].n_objectives();
572        let n_solutions = self.solutions.len();
573
574        let mut mean_objectives = Array1::zeros(n_objectives);
575        let mut min_objectives = Array1::from_elem(n_objectives, f64::INFINITY);
576        let mut max_objectives = Array1::from_elem(n_objectives, f64::NEG_INFINITY);
577
578        let mut feasible_count = 0;
579        let mut total_violation = 0.0;
580
581        // Calculate means, mins, maxs
582        for solution in &self.solutions {
583            for (i, &obj) in solution.objectives.iter().enumerate() {
584                mean_objectives[i] += obj;
585                min_objectives[i] = min_objectives[i].min(obj);
586                max_objectives[i] = max_objectives[i].max(obj);
587            }
588
589            if solution.is_feasible() {
590                feasible_count += 1;
591            }
592            total_violation += solution.constraint_violation;
593        }
594
595        mean_objectives /= n_solutions as f64;
596
597        // Calculate standard deviations
598        let mut std_objectives = Array1::zeros(n_objectives);
599        for solution in &self.solutions {
600            for (i, &obj) in solution.objectives.iter().enumerate() {
601                let diff = obj - mean_objectives[i];
602                std_objectives[i] += diff * diff;
603            }
604        }
605        std_objectives = std_objectives.mapv(|x: f64| (x / n_solutions as f64).sqrt());
606
607        PopulationStatistics {
608            mean_objectives,
609            std_objectives,
610            min_objectives,
611            max_objectives,
612            feasibility_ratio: feasible_count as f64 / n_solutions as f64,
613            avg_constraint_violation: total_violation / n_solutions as f64,
614        }
615    }
616}
617
618impl Default for Population {
619    fn default() -> Self {
620        Self::new()
621    }
622}
623
624#[cfg(test)]
625mod tests {
626    use super::*;
627    use scirs2_core::ndarray::array;
628
629    #[test]
630    fn test_solution_creation() {
631        let variables = array![1.0, 2.0, 3.0];
632        let objectives = array![4.0, 5.0];
633        let solution = MultiObjectiveSolution::new(variables.clone(), objectives.clone());
634
635        assert_eq!(solution.variables, variables);
636        assert_eq!(solution.objectives, objectives);
637        assert_eq!(solution.constraint_violation, 0.0);
638        assert!(solution.is_feasible());
639    }
640
641    #[test]
642    fn test_domination() {
643        let sol1 = MultiObjectiveSolution::new(array![1.0], array![1.0, 2.0]);
644        let sol2 = MultiObjectiveSolution::new(array![2.0], array![2.0, 3.0]);
645        let sol3 = MultiObjectiveSolution::new(array![3.0], array![0.5, 3.5]);
646
647        assert!(sol1.dominates(&sol2)); // sol1 is better in both objectives
648        assert!(!sol2.dominates(&sol1));
649        assert!(sol1.is_non_dominated_with(&sol3)); // Neither dominates the other
650    }
651
652    #[test]
653    fn test_constraint_domination() {
654        let feasible =
655            MultiObjectiveSolution::new_with_constraints(array![1.0], array![2.0, 3.0], 0.0);
656        let infeasible =
657            MultiObjectiveSolution::new_with_constraints(array![2.0], array![1.0, 2.0], 1.0);
658
659        assert!(feasible.dominates(&infeasible)); // Feasible dominates infeasible
660        assert!(!infeasible.dominates(&feasible));
661    }
662
663    #[test]
664    fn test_population_pareto_front() {
665        let mut population = Population::new();
666        population.add(MultiObjectiveSolution::new(array![1.0], array![1.0, 3.0]));
667        population.add(MultiObjectiveSolution::new(array![2.0], array![2.0, 2.0]));
668        population.add(MultiObjectiveSolution::new(array![3.0], array![3.0, 1.0]));
669        population.add(MultiObjectiveSolution::new(array![4.0], array![2.5, 2.5])); // Dominated
670
671        let pareto_front = population.extract_pareto_front();
672        assert_eq!(pareto_front.len(), 3); // First three solutions are non-dominated
673    }
674
675    #[test]
676    fn test_non_dominated_sorting() {
677        let mut population = Population::new();
678        population.add(MultiObjectiveSolution::new(array![1.0], array![1.0, 3.0]));
679        population.add(MultiObjectiveSolution::new(array![2.0], array![3.0, 1.0]));
680        population.add(MultiObjectiveSolution::new(array![3.0], array![2.0, 2.0]));
681        population.add(MultiObjectiveSolution::new(array![4.0], array![4.0, 4.0])); // Dominated
682
683        let fronts = population.non_dominated_sort();
684
685        assert_eq!(fronts.len(), 2); // Two fronts
686        assert_eq!(fronts[0].len(), 3); // First front has 3 solutions
687        assert_eq!(fronts[1].len(), 1); // Second front has 1 solution
688
689        // Check ranks
690        for &i in &fronts[0] {
691            assert_eq!(population.solutions[i].rank, 0);
692        }
693        for &i in &fronts[1] {
694            assert_eq!(population.solutions[i].rank, 1);
695        }
696    }
697
698    #[test]
699    fn test_crowding_distance() {
700        let mut population = Population::new();
701        population.add(MultiObjectiveSolution::new(array![1.0], array![1.0, 3.0]));
702        population.add(MultiObjectiveSolution::new(array![2.0], array![2.0, 2.0]));
703        population.add(MultiObjectiveSolution::new(array![3.0], array![3.0, 1.0]));
704
705        let front = vec![0, 1, 2];
706        population.calculate_crowding_distances(&front);
707
708        // Boundary solutions should have infinite distance
709        assert_eq!(population.solutions[0].crowding_distance, f64::INFINITY);
710        assert_eq!(population.solutions[2].crowding_distance, f64::INFINITY);
711        // Middle solution should have finite distance
712        assert!(population.solutions[1].crowding_distance.is_finite());
713        assert!(population.solutions[1].crowding_distance > 0.0);
714    }
715
716    #[test]
717    fn test_population_statistics() {
718        let mut population = Population::new();
719        population.add(MultiObjectiveSolution::new(array![1.0], array![1.0, 3.0]));
720        population.add(MultiObjectiveSolution::new(array![2.0], array![2.0, 2.0]));
721        population.add(MultiObjectiveSolution::new(array![3.0], array![3.0, 1.0]));
722
723        let stats = population.calculate_statistics();
724
725        assert_eq!(stats.mean_objectives[0], 2.0); // (1+2+3)/3
726        assert_eq!(stats.mean_objectives[1], 2.0); // (3+2+1)/3
727        assert_eq!(stats.min_objectives[0], 1.0);
728        assert_eq!(stats.max_objectives[0], 3.0);
729        assert_eq!(stats.feasibility_ratio, 1.0); // All feasible
730    }
731}
scirs2_optimize/multi_objective/solutions/mod.rs

scirs2_optimize/multi_objective/solutions/
mod.rs
