scirs2_optimize/multi_objective/indicators/
mod.rs

1//! Performance indicators for multi-objective optimization
2//!
3//! Metrics to evaluate the quality of Pareto fronts.
4
5use crate::multi_objective::solutions::Solution;
6
7/// Calculate the hypervolume indicator
8pub fn hypervolume(pareto_front: &[Solution], reference_point: &[f64]) -> f64 {
9    if pareto_front.is_empty() {
10        return 0.0;
11    }
12
13    let n_objectives = pareto_front[0].objectives.len();
14
15    // Simple 2D hypervolume calculation
16    if n_objectives == 2 {
17        hypervolume_2d(pareto_front, reference_point)
18    } else {
19        // TODO: Implement WFG algorithm for higher dimensions
20        hypervolume_monte_carlo(pareto_front, reference_point, 10000)
21    }
22}
23
24/// Calculate 2D hypervolume
25fn hypervolume_2d(pareto_front: &[Solution], reference_point: &[f64]) -> f64 {
26    let mut sorted_front = pareto_front.to_vec();
27    sorted_front.sort_by(|a, b| {
28        a.objectives[0]
29            .partial_cmp(&b.objectives[0])
30            .expect("Operation failed")
31    });
32
33    let mut volume = 0.0;
34    let mut prev_x = 0.0;
35
36    for solution in &sorted_front {
37        let x = solution.objectives[0];
38        let y = solution.objectives[1];
39
40        if x < reference_point[0] && y < reference_point[1] {
41            let width = x - prev_x;
42            let height = reference_point[1] - y;
43            volume += width * height;
44            prev_x = x;
45        }
46    }
47
48    // Add last rectangle
49    if prev_x < reference_point[0] {
50        let last_y = sorted_front.last().map(|s| s.objectives[1]).unwrap_or(0.0);
51        if last_y < reference_point[1] {
52            volume += (reference_point[0] - prev_x) * (reference_point[1] - last_y);
53        }
54    }
55
56    volume
57}
58
59/// Monte Carlo approximation for hypervolume
60fn hypervolume_monte_carlo(
61    pareto_front: &[Solution],
62    reference_point: &[f64],
63    n_samples: usize,
64) -> f64 {
65    use scirs2_core::random::Rng;
66    let mut rng = scirs2_core::random::rng();
67    let n_objectives = reference_point.len();
68
69    let mut count = 0;
70
71    for _ in 0..n_samples {
72        let point: Vec<f64> = (0..n_objectives)
73            .map(|i| rng.random::<f64>() * reference_point[i])
74            .collect();
75
76        if is_dominated_by_front(&point, pareto_front) {
77            count += 1;
78        }
79    }
80
81    let total_volume: f64 = reference_point.iter().product();
82    total_volume * (count as f64 / n_samples as f64)
83}
84
85fn is_dominated_by_front(point: &[f64], pareto_front: &[Solution]) -> bool {
86    for solution in pareto_front {
87        if dominates(
88            solution.objectives.as_slice().expect("Operation failed"),
89            point,
90        ) {
91            return true;
92        }
93    }
94    false
95}
96
97fn dominates(a: &[f64], b: &[f64]) -> bool {
98    let mut at_least_one_better = false;
99
100    for i in 0..a.len() {
101        if a[i] > b[i] {
102            return false;
103        }
104        if a[i] < b[i] {
105            at_least_one_better = true;
106        }
107    }
108
109    at_least_one_better
110}
111
112/// Calculate the Inverted Generational Distance (IGD)
113pub fn igd(pareto_front: &[Solution], true_pareto_front: &[Vec<f64>]) -> f64 {
114    if true_pareto_front.is_empty() {
115        return f64::INFINITY;
116    }
117
118    let sum: f64 = true_pareto_front
119        .iter()
120        .map(|true_point| {
121            pareto_front
122                .iter()
123                .map(|solution| {
124                    euclidean_distance(
125                        solution.objectives.as_slice().expect("Operation failed"),
126                        true_point,
127                    )
128                })
129                .fold(f64::INFINITY, |a, b| a.min(b))
130        })
131        .sum();
132
133    sum / true_pareto_front.len() as f64
134}
135
136/// Calculate the Generational Distance (GD)
137pub fn generational_distance(pareto_front: &[Solution], true_pareto_front: &[Vec<f64>]) -> f64 {
138    if pareto_front.is_empty() {
139        return f64::INFINITY;
140    }
141
142    let sum: f64 = pareto_front
143        .iter()
144        .map(|solution| {
145            true_pareto_front
146                .iter()
147                .map(|true_point| {
148                    euclidean_distance(
149                        solution.objectives.as_slice().expect("Operation failed"),
150                        true_point,
151                    )
152                })
153                .fold(f64::INFINITY, |a, b| a.min(b))
154        })
155        .sum();
156
157    sum / pareto_front.len() as f64
158}
159
160/// Calculate spacing indicator
161pub fn spacing(pareto_front: &[Solution]) -> f64 {
162    if pareto_front.len() < 2 {
163        return 0.0;
164    }
165
166    let distances: Vec<f64> = pareto_front
167        .iter()
168        .map(|sol| {
169            pareto_front
170                .iter()
171                .filter(|other| !std::ptr::eq(*other, sol))
172                .map(|other| {
173                    euclidean_distance(
174                        sol.objectives.as_slice().expect("Operation failed"),
175                        other.objectives.as_slice().expect("Operation failed"),
176                    )
177                })
178                .fold(f64::INFINITY, |a, b| a.min(b))
179        })
180        .collect();
181
182    let mean_distance = distances.iter().sum::<f64>() / distances.len() as f64;
183
184    let variance = distances
185        .iter()
186        .map(|d| (d - mean_distance).powi(2))
187        .sum::<f64>()
188        / distances.len() as f64;
189
190    variance.sqrt()
191}
192
193fn euclidean_distance(a: &[f64], b: &[f64]) -> f64 {
194    a.iter()
195        .zip(b.iter())
196        .map(|(x, y)| (x - y).powi(2))
197        .sum::<f64>()
198        .sqrt()
199}
scirs2_optimize/multi_objective/indicators/mod.rs

scirs2_optimize/multi_objective/indicators/
mod.rs
