scirs2_optimize/multi_objective/indicators/
mod.rs1use crate::multi_objective::solutions::Solution;
6
7pub 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 if n_objectives == 2 {
17 hypervolume_2d(pareto_front, reference_point)
18 } else {
19 hypervolume_monte_carlo(pareto_front, reference_point, 10000)
21 }
22}
23
24fn 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| a.objectives[0].partial_cmp(&b.objectives[0]).unwrap());
28
29 let mut volume = 0.0;
30 let mut prev_x = 0.0;
31
32 for solution in &sorted_front {
33 let x = solution.objectives[0];
34 let y = solution.objectives[1];
35
36 if x < reference_point[0] && y < reference_point[1] {
37 let width = x - prev_x;
38 let height = reference_point[1] - y;
39 volume += width * height;
40 prev_x = x;
41 }
42 }
43
44 if prev_x < reference_point[0] {
46 let last_y = sorted_front.last().map(|s| s.objectives[1]).unwrap_or(0.0);
47 if last_y < reference_point[1] {
48 volume += (reference_point[0] - prev_x) * (reference_point[1] - last_y);
49 }
50 }
51
52 volume
53}
54
55fn hypervolume_monte_carlo(
57 pareto_front: &[Solution],
58 reference_point: &[f64],
59 n_samples: usize,
60) -> f64 {
61 use scirs2_core::random::Rng;
62 let mut rng = scirs2_core::random::rng();
63 let n_objectives = reference_point.len();
64
65 let mut count = 0;
66
67 for _ in 0..n_samples {
68 let point: Vec<f64> = (0..n_objectives)
69 .map(|i| rng.random::<f64>() * reference_point[i])
70 .collect();
71
72 if is_dominated_by_front(&point, pareto_front) {
73 count += 1;
74 }
75 }
76
77 let total_volume: f64 = reference_point.iter().product();
78 total_volume * (count as f64 / n_samples as f64)
79}
80
81fn is_dominated_by_front(point: &[f64], pareto_front: &[Solution]) -> bool {
82 for solution in pareto_front {
83 if dominates(solution.objectives.as_slice().unwrap(), point) {
84 return true;
85 }
86 }
87 false
88}
89
90fn dominates(a: &[f64], b: &[f64]) -> bool {
91 let mut at_least_one_better = false;
92
93 for i in 0..a.len() {
94 if a[i] > b[i] {
95 return false;
96 }
97 if a[i] < b[i] {
98 at_least_one_better = true;
99 }
100 }
101
102 at_least_one_better
103}
104
105pub fn igd(pareto_front: &[Solution], true_pareto_front: &[Vec<f64>]) -> f64 {
107 if true_pareto_front.is_empty() {
108 return f64::INFINITY;
109 }
110
111 let sum: f64 = true_pareto_front
112 .iter()
113 .map(|true_point| {
114 pareto_front
115 .iter()
116 .map(|solution| {
117 euclidean_distance(solution.objectives.as_slice().unwrap(), true_point)
118 })
119 .fold(f64::INFINITY, |a, b| a.min(b))
120 })
121 .sum();
122
123 sum / true_pareto_front.len() as f64
124}
125
126pub fn generational_distance(pareto_front: &[Solution], true_pareto_front: &[Vec<f64>]) -> f64 {
128 if pareto_front.is_empty() {
129 return f64::INFINITY;
130 }
131
132 let sum: f64 = pareto_front
133 .iter()
134 .map(|solution| {
135 true_pareto_front
136 .iter()
137 .map(|true_point| {
138 euclidean_distance(solution.objectives.as_slice().unwrap(), true_point)
139 })
140 .fold(f64::INFINITY, |a, b| a.min(b))
141 })
142 .sum();
143
144 sum / pareto_front.len() as f64
145}
146
147pub fn spacing(pareto_front: &[Solution]) -> f64 {
149 if pareto_front.len() < 2 {
150 return 0.0;
151 }
152
153 let distances: Vec<f64> = pareto_front
154 .iter()
155 .map(|sol| {
156 pareto_front
157 .iter()
158 .filter(|other| !std::ptr::eq(*other, sol))
159 .map(|other| {
160 euclidean_distance(
161 sol.objectives.as_slice().unwrap(),
162 other.objectives.as_slice().unwrap(),
163 )
164 })
165 .fold(f64::INFINITY, |a, b| a.min(b))
166 })
167 .collect();
168
169 let mean_distance = distances.iter().sum::<f64>() / distances.len() as f64;
170
171 let variance = distances
172 .iter()
173 .map(|d| (d - mean_distance).powi(2))
174 .sum::<f64>()
175 / distances.len() as f64;
176
177 variance.sqrt()
178}
179
180fn euclidean_distance(a: &[f64], b: &[f64]) -> f64 {
181 a.iter()
182 .zip(b.iter())
183 .map(|(x, y)| (x - y).powi(2))
184 .sum::<f64>()
185 .sqrt()
186}