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| 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    // Add last rectangle
45    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
55/// Monte Carlo approximation for hypervolume
56fn 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
105/// Calculate the Inverted Generational Distance (IGD)
106pub 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
126/// Calculate the Generational Distance (GD)
127pub 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
147/// Calculate spacing indicator
148pub 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}