hpr_analysis/optimize/benchmark.rs
1//! Test functions with known minima, for checking an optimizer.
2//!
3//! Each is a test function of CMA-ES's authors, from N. Hansen, S. D. Müller and P. Koumoutsakos,
4//! "Reducing the time complexity of the derandomized evolution strategy with covariance matrix
5//! adaptation (CMA-ES)", *Evolutionary Computation* 11(1), 1–18 (2003),
6//! <https://doi.org/10.1162/106365603321828970>, Table 1 (printed p. 7), with the variables
7//! numbered from 0 here. Their minima are known exactly, so a run's best point and value can be
8//! held to them; that paper runs each to `f = 10⁻¹⁰`, as the tests do.
9//!
10//! | Function | What it tests | Minimum |
11//! |---|---|---|
12//! | [`sphere`] | the step size alone | 0 at `x = 0` |
13//! | [`ellipsoid`] | coefficients from 1 to 10⁶, so the distribution must grow 1,000 times longer one way than the other | 0 at `x = 0` |
14//! | [`rotated_ellipsoid`] | the same along axes that aren't the variables' | 0 at `x = 0` |
15//! | [`rosenbrock`] | following a long, curved valley | 0 at `x = 1` |
16//!
17//! [`constrained`] holds three problems whose minima lie on their constraints' edges.
18//! [`mixed`] holds three whose second half of variables take only whole numbers.
19
20/// The ellipsoid's condition number: the ratio of its largest curvature to its smallest.
21pub const ELLIPSOID_CONDITION: f64 = 1e6;
22
23/// The sphere, `f(x) = Σ xᵢ²`.
24pub fn sphere(x: &[f64]) -> f64 {
25 x.iter().map(|xi| xi * xi).sum()
26}
27
28/// The ellipsoid, `f(x) = Σ 10^(6 i/(n − 1)) xᵢ²` for `i = 0 … n − 1`: its curvatures run from 1
29/// to [`ELLIPSOID_CONDITION`]. With one variable it is the sphere.
30pub fn ellipsoid(x: &[f64]) -> f64 {
31 let n = x.len();
32 if n < 2 {
33 return sphere(x);
34 }
35 let denominator = (n - 1) as f64;
36 x.iter()
37 .enumerate()
38 .map(|(i, xi)| ELLIPSOID_CONDITION.powf(i as f64 / denominator) * xi * xi)
39 .sum()
40}
41
42/// The ellipsoid turned so that its axes aren't the variables': `ellipsoid(H x)`, with `H` the
43/// reflection `I − 2 v vᵀ/(vᵀv)`, `vᵢ = i + 1`. `H` is orthogonal, so the minimum is still 0 at
44/// `x = 0`, but no variable can be stepped alone to reach it. CMA-ES is invariant under rotations
45/// and reflections of the variables, so it should take about as many evaluations as on
46/// [`ellipsoid`].
47pub fn rotated_ellipsoid(x: &[f64]) -> f64 {
48 let vv: f64 = (1..=x.len()).map(|i| (i * i) as f64).sum();
49 let vx: f64 = x
50 .iter()
51 .enumerate()
52 .map(|(i, xi)| (i + 1) as f64 * xi)
53 .sum();
54 let factor = 2.0 * vx / vv;
55 let hx: Vec<f64> = x
56 .iter()
57 .enumerate()
58 .map(|(i, xi)| xi - factor * (i + 1) as f64)
59 .collect();
60 ellipsoid(&hx)
61}
62
63/// Rosenbrock's function, `f(x) = Σ [100 (xᵢ² − xᵢ₊₁)² + (1 − xᵢ)²]` for `i = 0 … n − 2`. Its
64/// global minimum is 0 at `x = 1`. "For higher dimension, even function 8 f_Rosen has a local
65/// minimum near y = (−1, 1, …, 1)ᵀ" (N. Hansen and A. Ostermeier, *Evolutionary Computation* 9(2), 159–195,
66/// 2001, <https://doi.org/10.1162/106365601750190398>, footnote 18), and CMA-ES sometimes misses
67/// the global one: in 1 to 3 of 20 runs at 4 to 16 variables in S. Kern, N. Hansen and
68/// P. Koumoutsakos, "Local meta-models for optimization using evolution strategies", *PPSN IX*
69/// (2006), Table 3.
70pub fn rosenbrock(x: &[f64]) -> f64 {
71 x.windows(2)
72 .map(|w| 100.0 * (w[0] * w[0] - w[1]).powi(2) + (1.0 - w[0]).powi(2))
73 .sum()
74}
75
76#[cfg(test)]
77mod tests {
78 use super::*;
79
80 #[test]
81 fn minima_are_zero_where_stated() {
82 let zero = [0.0; 10];
83 let one = [1.0; 10];
84 assert_eq!(sphere(&zero), 0.0);
85 assert_eq!(ellipsoid(&zero), 0.0);
86 assert_eq!(rotated_ellipsoid(&zero), 0.0);
87 assert_eq!(rosenbrock(&one), 0.0);
88 }
89
90 #[test]
91 fn values_at_a_point_by_hand() {
92 // Two variables: the ellipsoid's curvatures are 1 and 1e6.
93 assert_eq!(ellipsoid(&[1.0, 1.0]), 1.0 + 1e6);
94 assert_eq!(ellipsoid(&[3.0]), 9.0);
95 // Rosenbrock at the origin: one term per pair, each 100·0 + 1.
96 assert_eq!(rosenbrock(&[0.0; 10]), 9.0);
97 // H with v = (1, 2) maps e₁ to e₁ − (2/5)(1, 2) = (0.6, −0.8).
98 let turned = rotated_ellipsoid(&[1.0, 0.0]);
99 let by_hand = 0.6f64.powi(2) + 1e6 * 0.8f64.powi(2);
100 assert!((turned - by_hand).abs() <= 1e-9, "{turned} vs {by_hand}");
101 }
102
103 /// The reflection keeps lengths: on the sphere's terms (all curvatures 1) the turned and
104 /// unturned values agree, checked through a one-variable-at-a-time identity.
105 #[test]
106 fn reflection_keeps_lengths() {
107 let x = [0.3, -1.2, 2.5, 0.7];
108 let vv = 30.0;
109 let vx: f64 = x
110 .iter()
111 .enumerate()
112 .map(|(i, xi)| (i + 1) as f64 * xi)
113 .sum();
114 let hx: Vec<f64> = x
115 .iter()
116 .enumerate()
117 .map(|(i, xi)| xi - 2.0 * vx / vv * (i + 1) as f64)
118 .collect();
119 assert!((sphere(&hx) - sphere(&x)).abs() <= 1e-14);
120 }
121}
122
123/// Constrained test problems with known minima on their constraints' edges, for
124/// [`Run::tell_constrained`](super::cmaes::Run::tell_constrained). Each gives an
125/// [`Evaluation`](super::Evaluation) with constraints written `g(x) ≤ 0`.
126///
127/// | Problem | Constraint | Minimum |
128/// |---|---|---|
129/// | [`constrained::sphere_above`] | `x₀ ≥ 1` | 1 at `x = (1, 0, …, 0)` |
130/// | [`constrained::tangent`] | `Σ xᵢ ≥ n` | `n` at `x = 1` |
131/// | [`constrained::g06`] | two circles | `(x₀* − 10)³ + (x₁* − 20)³` at their crossing |
132pub mod constrained {
133 use crate::optimize::Evaluation;
134
135 /// The sphere `Σ xᵢ²` with `x₀ ≥ 1`, as `g = 1 − x₀ ≤ 0`. Its minimum is 1, at
136 /// `x = (1, 0, …, 0)`: the constraint binds, as the sphere's own minimum breaks it.
137 pub fn sphere_above(x: &[f64]) -> Evaluation {
138 let g = 1.0 - x.first().copied().unwrap_or(0.0);
139 Evaluation::constrained(super::sphere(x), &[g])
140 }
141
142 /// The tangent problem: the sphere with `Σ xᵢ ≥ n`, as `g = n − Σ xᵢ ≤ 0`, a constraint along
143 /// no variable's axis. By symmetry and Lagrange's condition `2 xᵢ = λ`, the minimum is `n`, at
144 /// `x = 1`.
145 pub fn tangent(x: &[f64]) -> Evaluation {
146 let n = x.len() as f64;
147 let g = n - x.iter().sum::<f64>();
148 Evaluation::constrained(super::sphere(x), &[g])
149 }
150
151 /// Problem g06 of the CEC 2006 constrained benchmark (J. J. Liang et al., "Problem
152 /// definitions and evaluation criteria for the CEC 2006 special session on constrained
153 /// real-parameter optimization", Nanyang Technological University (2006)), two variables:
154 /// `f = (x₀ − 10)³ + (x₁ − 20)³`, with `g₁ = 100 − (x₀ − 5)² − (x₁ − 5)² ≤ 0` (outside one
155 /// circle) and `g₂ = (x₀ − 6)² + (x₁ − 5)² − 82.81 ≤ 0` (inside another), and bounds
156 /// `13 ≤ x₀ ≤ 100`, `0 ≤ x₁ ≤ 100` for the caller to set. The minimum is where the circles
157 /// cross: subtracting the two edges gives `2 x₀ − 11 = 17.19`, so [`G06_X0`], and
158 /// [`g06_x1`] below the centers.
159 pub fn g06(x: &[f64]) -> Evaluation {
160 let (a, b) = (
161 x.first().copied().unwrap_or(0.0),
162 x.get(1).copied().unwrap_or(0.0),
163 );
164 let f = (a - 10.0).powi(3) + (b - 20.0).powi(3);
165 let g1 = 100.0 - (a - 5.0).powi(2) - (b - 5.0).powi(2);
166 let g2 = (a - 6.0).powi(2) + (b - 5.0).powi(2) - 82.81;
167 Evaluation::constrained(f, &[g1, g2])
168 }
169
170 /// g06's minimum's first coordinate, `x₀* = 14.095`.
171 pub const G06_X0: f64 = 14.095;
172
173 /// g06's minimum's second coordinate, `x₁* = 5 − √(100 − (x₀* − 5)²)`.
174 pub fn g06_x1() -> f64 {
175 5.0 - (100.0 - (G06_X0 - 5.0).powi(2)).sqrt()
176 }
177}
178
179/// Test functions of continuous and integer variables together: those of R. Hamano, S. Saito,
180/// M. Nomura and S. Shirakawa, "CMA-ES with Margin: Lower-Bounding Marginal Probability for
181/// Mixed-Integer Black-Box Optimization", GECCO 2022, <https://arxiv.org/abs/2205.13482>, §5.1
182/// (p. 7). The first `⌊n/2⌋` variables are continuous and the rest integer
183/// ([`Variable::integer`](super::Variable::integer)); each function is given the values the
184/// optimizer encodes, whole numbers in the integer variables.
185///
186/// | Function | Integer variables | Minimum |
187/// |---|---|---|
188/// | [`sphere_int`](mixed::sphere_int) | from −10 to 10 | 0 at `x = 0` |
189/// | [`ellipsoid_int`](mixed::ellipsoid_int) | from −10 to 10, with the largest coefficients | 0 at `x = 0` |
190/// | [`sphere_one_max`](mixed::sphere_one_max) | 0 or 1 | 0 at continuous 0, integer 1 |
191pub mod mixed {
192 use super::{ellipsoid, sphere};
193
194 /// SphereInt, `f(x) = Σ xᵢ²` over every variable: [`sphere`].
195 pub fn sphere_int(x: &[f64]) -> f64 {
196 sphere(x)
197 }
198
199 /// EllipsoidInt, `f(x) = Σ (1000^(i/(n − 1)) xᵢ)²`: [`ellipsoid`], whose coefficients are the
200 /// same, so the integer variables, the last, have the largest.
201 pub fn ellipsoid_int(x: &[f64]) -> f64 {
202 ellipsoid(x)
203 }
204
205 /// SphereOneMax, `f(x) = Σ xᵢ² + n_b − Σ x_k`, the first sum over the first `⌊n/2⌋`
206 /// variables (continuous) and the second over the `n_b` others (0 or 1).
207 pub fn sphere_one_max(x: &[f64]) -> f64 {
208 let (continuous, binary) = x.split_at(x.len() / 2);
209 // Cast: at most 200 variables.
210 sphere(continuous) + binary.len() as f64 - binary.iter().sum::<f64>()
211 }
212}
213
214/// Test functions for global optimization with few evaluations, for [`ego`](super::ego): the
215/// Branin and Hartmann functions of L. C. W. Dixon and G. P. Szegö, "The global optimisation
216/// problem: an introduction", in *Towards Global Optimisation 2*, North-Holland, 1–15 (1978),
217/// on which D. R. Jones, M. Schonlau and W. J. Welch, "Efficient global optimization of
218/// expensive black-box functions", *Journal of Global Optimization* 13, 455–492 (1998),
219/// <https://doi.org/10.1023/A:1008306431147>, run EGO (their Table 1). The constants are as
220/// S. Surjanovic and D. Bingham's Virtual Library of Simulation Experiments lists them,
221/// <https://www.sfu.ca/~ssurjano/optimization.html>. Branin has three minima, all of the same
222/// value; the Hartmann functions have local minima above their least value. A missing variable
223/// is taken as 0.
224///
225/// | Function | Variables | Minimum |
226/// |---|---|---|
227/// | [`branin`](global::branin) | `x₀` in `[−5, 10]`, `x₁` in `[0, 15]` | `5/(4π) ≈ 0.397887` at three points ([`BRANIN_MINIMA`](global::BRANIN_MINIMA)) |
228/// | [`hartmann3`](global::hartmann3) | 3, each in `[0, 1]` | `≈ −3.86278` ([`HARTMANN3_MINIMUM`](global::HARTMANN3_MINIMUM)) |
229/// | [`hartmann6`](global::hartmann6) | 6, each in `[0, 1]` | `≈ −3.32237` ([`HARTMANN6_MINIMUM`](global::HARTMANN6_MINIMUM)) |
230pub mod global {
231 use std::f64::consts::PI;
232
233 /// Branin's minimum, `s t = 10/(8π) = 5/(4π)`: the squared term vanishes and `cos x₀ = −1`.
234 pub const BRANIN_MINIMUM: f64 = 5.0 / (4.0 * PI);
235
236 /// Branin's three minima: `x₀ = −π, π, 3π`, where `cos x₀ = −1`, and `x₁` the root of its
237 /// squared term, `b x₀² − c x₀ + r` (12.275, 2.275 and 2.475).
238 pub const BRANIN_MINIMA: [[f64; 2]; 3] = [[-PI, 12.275], [PI, 2.275], [3.0 * PI, 2.475]];
239
240 /// Branin's function of two variables, `x₀` in `[−5, 10]` and `x₁` in `[0, 15]`:
241 /// `f = (x₁ − b x₀² + c x₀ − r)² + s (1 − t) cos x₀ + s`, with `b = 5.1/(4π²)`, `c = 5/π`,
242 /// `r = 6`, `s = 10` and `t = 1/(8π)`.
243 pub fn branin(x: &[f64]) -> f64 {
244 let (b, c, r, s, t) = (5.1 / (4.0 * PI * PI), 5.0 / PI, 6.0, 10.0, 1.0 / (8.0 * PI));
245 let (x0, x1) = (
246 x.first().copied().unwrap_or(0.0),
247 x.get(1).copied().unwrap_or(0.0),
248 );
249 let square = x1 - b * x0 * x0 + c * x0 - r;
250 square * square + s * (1.0 - t) * x0.cos() + s
251 }
252
253 /// The Hartmann 3 function's least value, as printed to six figures.
254 pub const HARTMANN3_MINIMUM: f64 = -3.86278;
255
256 /// The Hartmann 6 function's least value, as printed to six figures.
257 pub const HARTMANN6_MINIMUM: f64 = -3.32237;
258
259 /// The Hartmann functions' weights `αᵢ`.
260 const ALPHA: [f64; 4] = [1.0, 1.2, 3.0, 3.2];
261
262 /// `f = −Σᵢ αᵢ exp(−Σⱼ aᵢⱼ (xⱼ − pᵢⱼ)²)`, the Hartmann functions' form.
263 fn hartmann<const N: usize>(x: &[f64], a: &[[f64; N]; 4], p: &[[f64; N]; 4]) -> f64 {
264 -(0..4)
265 .map(|i| {
266 let exponent: f64 = (0..N)
267 .map(|j| {
268 let d = x.get(j).copied().unwrap_or(0.0) - p[i][j];
269 a[i][j] * d * d
270 })
271 .sum();
272 ALPHA[i] * (-exponent).exp()
273 })
274 .sum::<f64>()
275 }
276
277 /// The Hartmann 3 function, each variable in `[0, 1]`; its least value is about −3.86278,
278 /// near `(0.114614, 0.555649, 0.852547)`. The last center's first coordinate is 0.0381 here,
279 /// as the library lists it; some listings print 0.03815, the value that quoted point belongs
280 /// to. The two least values differ by about 2×10⁻⁶ (both −3.86278 to six figures), as that
281 /// coordinate's weight, 0.1, is small.
282 pub fn hartmann3(x: &[f64]) -> f64 {
283 const A: [[f64; 3]; 4] = [
284 [3.0, 10.0, 30.0],
285 [0.1, 10.0, 35.0],
286 [3.0, 10.0, 30.0],
287 [0.1, 10.0, 35.0],
288 ];
289 const P: [[f64; 3]; 4] = [
290 [0.3689, 0.1170, 0.2673],
291 [0.4699, 0.4387, 0.7470],
292 [0.1091, 0.8732, 0.5547],
293 [0.0381, 0.5743, 0.8828],
294 ];
295 hartmann(x, &A, &P)
296 }
297
298 /// The Hartmann 6 function, each variable in `[0, 1]`; its least value is about −3.32237,
299 /// near `(0.20169, 0.150011, 0.476874, 0.275332, 0.311652, 0.6573)`.
300 pub fn hartmann6(x: &[f64]) -> f64 {
301 const A: [[f64; 6]; 4] = [
302 [10.0, 3.0, 17.0, 3.5, 1.7, 8.0],
303 [0.05, 10.0, 17.0, 0.1, 8.0, 14.0],
304 [3.0, 3.5, 1.7, 10.0, 17.0, 8.0],
305 [17.0, 8.0, 0.05, 10.0, 0.1, 14.0],
306 ];
307 const P: [[f64; 6]; 4] = [
308 [0.1312, 0.1696, 0.5569, 0.0124, 0.8283, 0.5886],
309 [0.2329, 0.4135, 0.8307, 0.3736, 0.1004, 0.9991],
310 [0.2348, 0.1451, 0.3522, 0.2883, 0.3047, 0.6650],
311 [0.4047, 0.8828, 0.8732, 0.5743, 0.1091, 0.0381],
312 ];
313 hartmann(x, &A, &P)
314 }
315}
316
317/// Two-goal test problems with known Pareto fronts, for [`nsga2`](super::nsga2): ZDT1, ZDT2 and
318/// ZDT3 of E. Zitzler, K. Deb and L. Thiele, "Comparison of multiobjective evolutionary
319/// algorithms: empirical results", *Evolutionary Computation* 8(2), 173–195 (2000),
320/// <https://doi.org/10.1162/106365600568202>, §4, eqs. (7)–(9) (pp. 177–178), each of `n` variables in `[0, 1]`
321/// (the paper's `n = 30`):
322///
323/// `f₁ = x₀`, `g = 1 + 9 Σᵢ₌₁ⁿ⁻¹ xᵢ / (n − 1)`, `f₂ = g h(f₁, g)`.
324///
325/// The front is `g = 1`, every variable but the first zero, where `f₂ = h(f₁, 1)`:
326///
327/// | Problem | `h(f₁, g)` | The front |
328/// |---|---|---|
329/// | ZDT1 | `1 − √(f₁/g)` | `f₂ = 1 − √f₁`, convex, `f₁` from 0 to 1 |
330/// | ZDT2 | `1 − (f₁/g)²` | `f₂ = 1 − f₁²`, concave |
331/// | ZDT3 | `1 − √(f₁/g) − (f₁/g) sin(10π f₁)` | five separate pieces of `f₂ = 1 − √f₁ − f₁ sin(10π f₁)` ([`zdt::ZDT3_PIECES`]) |
332pub mod zdt {
333 /// One of the three problems.
334 #[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
335 #[non_exhaustive]
336 pub enum Zdt {
337 /// ZDT1: a convex front.
338 One,
339 /// ZDT2: a concave front.
340 Two,
341 /// ZDT3: a front in five pieces.
342 Three,
343 }
344
345 /// The `f₁` ranges of ZDT3's five pieces of front. Along `f₂ = 1 − √f₁ − f₁ sin(10π f₁)`,
346 /// a point is on the front only if `f₂` is below its value at every smaller `f₁`: each piece
347 /// ends at a local minimum of the curve (where its slope is zero), and the next starts where
348 /// the curve, falling again, first drops below that minimum. These were solved to 40 digits
349 /// (mpmath 1.3.0's `findroot`), here rounded to the nearest `f64`, and agree with the
350 /// ten-digit ranges published for ZDT3; the tests check each end's equation.
351 pub const ZDT3_PIECES: [(f64, f64); 5] = [
352 (0.0, 0.083_001_534_926_911_63),
353 (0.182_228_728_029_399_77, 0.257_762_363_387_830_2),
354 (0.409_313_674_808_656_8, 0.453_882_104_088_830_2),
355 (0.618_396_794_439_265_8, 0.652_511_703_804_662_5),
356 (0.823_331_798_326_632_7, 0.851_832_865_436_413_9),
357 ];
358
359 /// The whole range of `f₁`, ZDT1's and ZDT2's one piece of front.
360 const WHOLE: [(f64, f64); 1] = [(0.0, 1.0)];
361
362 /// `g = 1 + 9 Σᵢ₌₁ⁿ⁻¹ xᵢ / (n − 1)`; 1 for a single variable.
363 fn g(x: &[f64]) -> f64 {
364 let rest = x.get(1..).unwrap_or(&[]);
365 if rest.is_empty() {
366 return 1.0;
367 }
368 // Cast: at most 200 variables.
369 1.0 + 9.0 * rest.iter().sum::<f64>() / rest.len() as f64
370 }
371
372 impl Zdt {
373 /// `h(f₁, g)`.
374 fn h(self, f1: f64, g: f64) -> f64 {
375 let r = f1 / g;
376 match self {
377 Self::One => 1.0 - r.sqrt(),
378 Self::Two => 1.0 - r * r,
379 Self::Three => 1.0 - r.sqrt() - r * (10.0 * std::f64::consts::PI * f1).sin(),
380 }
381 }
382
383 /// The problem's two goals at `x`, `[f₁, f₂]`; `x₀` is taken as 0 if `x` is empty.
384 pub fn evaluate(self, x: &[f64]) -> Vec<f64> {
385 let f1 = x.first().copied().unwrap_or(0.0);
386 let g = g(x);
387 vec![f1, g * self.h(f1, g)]
388 }
389
390 /// The curve the front lies on, `f₂ = h(f₁, 1)`, at `f₁`.
391 pub fn front_f2(self, f1: f64) -> f64 {
392 self.h(f1, 1.0)
393 }
394
395 /// The `f₁` ranges of the front's pieces: one for ZDT1 and ZDT2, five for ZDT3.
396 pub fn pieces(self) -> &'static [(f64, f64)] {
397 match self {
398 Self::One | Self::Two => &WHOLE,
399 Self::Three => &ZDT3_PIECES,
400 }
401 }
402
403 /// `points` points of the front, `f₁` evenly spaced from 0 to 1, those outside ZDT3's
404 /// pieces left out: a reference for
405 /// [`inverted_generational_distance`](crate::optimize::nsga2::inverted_generational_distance).
406 pub fn reference(self, points: usize) -> Vec<Vec<f64>> {
407 // Cast: a count of points, exact in f64 below 2⁵³.
408 let last = points.saturating_sub(1).max(1) as f64;
409 (0..points)
410 .map(|i| i as f64 / last)
411 .filter(|&f1| {
412 self.pieces()
413 .iter()
414 .any(|&(lo, hi)| (lo..=hi).contains(&f1))
415 })
416 .map(|f1| vec![f1, self.front_f2(f1)])
417 .collect()
418 }
419
420 /// The Euclidean distance from goals `f = [f₁, f₂]` to the front, the curve itself
421 /// rather than points along it.
422 ///
423 /// Any point of the front bounds it, so the nearest point lies within `f₁ ± d₀` of `f₁`,
424 /// where `d₀` is the distance to the front at `f₁` itself (or a piece's nearer end):
425 /// that window of each piece is searched on a grid of 400 steps in `s` (`f₁ = s²` for
426 /// ZDT1 and ZDT3, whose slope is infinite at 0; `f₁ = s` for ZDT2), and the best step
427 /// refined by golden-section search over its two neighbours.
428 pub fn distance_to_front(self, f: &[f64]) -> f64 {
429 const STEPS: usize = 400;
430 let (a, b) = (
431 f.first().copied().unwrap_or(0.0),
432 f.get(1).copied().unwrap_or(0.0),
433 );
434 let distance = |f1: f64| (a - f1).hypot(b - self.front_f2(f1));
435 let squared = self != Self::Two;
436 let to_f1 = |s: f64| if squared { s * s } else { s };
437 let to_s = |f1: f64| if squared { f1.sqrt() } else { f1 };
438 let d0 = self
439 .pieces()
440 .iter()
441 .map(|&(lo, hi)| distance(a.clamp(lo, hi)))
442 .fold(f64::INFINITY, f64::min);
443 let mut best = d0;
444 for &(lo, hi) in self.pieces() {
445 let (left, right) = ((a - d0).max(lo), (a + d0).min(hi));
446 if left > right {
447 continue;
448 }
449 let (s0, s1) = (to_s(left), to_s(right));
450 // Cast: STEPS is small and exact in f64.
451 let at = |k: usize| s0 + (s1 - s0) * k as f64 / STEPS as f64;
452 let on = |s: f64| distance(to_f1(s));
453 let (mut k_best, mut d_best) = (0, f64::INFINITY);
454 for k in 0..=STEPS {
455 let d = on(at(k));
456 if d < d_best {
457 (k_best, d_best) = (k, d);
458 }
459 }
460 let (mut lo_s, mut hi_s) =
461 (at(k_best.saturating_sub(1)), at((k_best + 1).min(STEPS)));
462 let ratio = (5.0_f64.sqrt() - 1.0) / 2.0;
463 for _ in 0..100 {
464 let p = hi_s - ratio * (hi_s - lo_s);
465 let q = lo_s + ratio * (hi_s - lo_s);
466 if on(p) <= on(q) {
467 hi_s = q;
468 } else {
469 lo_s = p;
470 }
471 }
472 best = best.min(d_best).min(on(0.5 * (lo_s + hi_s)));
473 }
474 best
475 }
476 }
477
478 #[cfg(test)]
479 mod tests {
480 use super::*;
481
482 use std::f64::consts::PI;
483
484 /// The slope of ZDT3's front curve, `dh/df₁ = −1/(2√f₁) − sin(10π f₁) − 10π f₁ cos(10π f₁)`.
485 fn slope(f1: f64) -> f64 {
486 -0.5 / f1.sqrt() - (10.0 * PI * f1).sin() - 10.0 * PI * f1 * (10.0 * PI * f1).cos()
487 }
488
489 /// Each piece of ZDT3's front ends where the curve's slope is zero, and the next starts
490 /// at the same height, to rounding; the ends agree with the ten-digit ranges pymoo
491 /// 0.6.2's `zdt.py` prints, but for its second start, 0.182228780, a digit short of
492 /// 0.1822287280.
493 #[test]
494 fn zdt3_pieces_meet_their_equations() {
495 let h = |f1| Zdt::Three.front_f2(f1);
496 for (k, &(lo, hi)) in ZDT3_PIECES.iter().enumerate() {
497 assert!(
498 slope(hi).abs() <= 1e-12,
499 "piece {k} end: slope {}",
500 slope(hi)
501 );
502 if k > 0 {
503 let previous = ZDT3_PIECES[k - 1].1;
504 assert!((h(lo) - h(previous)).abs() <= 1e-15, "piece {k} start");
505 assert!(slope(lo) < 0.0);
506 }
507 }
508 let published = [
509 (0.0, 0.083_001_534_9),
510 (0.182_228_728_0, 0.257_762_363_4),
511 (0.409_313_674_8, 0.453_882_104_1),
512 (0.618_396_794_4, 0.652_511_703_8),
513 (0.823_331_798_3, 0.851_832_865_4),
514 ];
515 for (&(lo, hi), (plo, phi)) in ZDT3_PIECES.iter().zip(published) {
516 assert!((lo - plo).abs() <= 5e-11 && (hi - phi).abs() <= 5e-11);
517 }
518 }
519
520 /// On a scan of 200,001 points of the curve and the pieces' ends, a point is a new
521 /// lowest `f₂` exactly when it is inside a piece. (A piece's start ties the previous end's
522 /// `f₂`, so is itself dominated: the ends only set the lowest.)
523 #[test]
524 fn zdt3_pieces_are_the_curves_undominated_points() {
525 let ends: Vec<f64> = ZDT3_PIECES.iter().flat_map(|&(lo, hi)| [lo, hi]).collect();
526 let mut scan: Vec<f64> = (0..=200_000).map(|i| f64::from(i) / 200_000.0).collect();
527 scan.extend(&ends);
528 scan.sort_by(f64::total_cmp);
529 let mut lowest = f64::INFINITY;
530 for f1 in scan {
531 let f2 = Zdt::Three.front_f2(f1);
532 if !ends.contains(&f1) {
533 let inside = ZDT3_PIECES.iter().any(|&(lo, hi)| (lo..=hi).contains(&f1));
534 assert_eq!(f2 < lowest, inside, "f₁ = {f1}");
535 }
536 lowest = lowest.min(f2);
537 }
538 }
539
540 #[test]
541 fn the_front_is_g_equal_to_one() {
542 let mut x = vec![0.0; 30];
543 for f1 in [0.0, 0.04, 0.2, 0.5, 0.83, 1.0] {
544 x[0] = f1;
545 for p in [Zdt::One, Zdt::Two, Zdt::Three] {
546 let f = p.evaluate(&x);
547 assert_eq!(f, vec![f1, p.front_f2(f1)]);
548 }
549 }
550 x[5] = 0.29;
551 // g = 1 + 9 × 0.29 / 29 = 1.09.
552 let f = Zdt::Two.evaluate(&x);
553 let r = 1.0 / 1.09;
554 assert!((f[1] - 1.09 * (1.0 - r * r)).abs() <= 1e-15);
555 }
556
557 /// The distance to the front: zero on it; by hand for a point straight out from ZDT2's
558 /// curve along its normal; and never more than a brute-force scan of 10⁶ points of it,
559 /// nor less than that scan by more than its spacing could hide.
560 #[test]
561 fn distance_to_the_front() {
562 for p in [Zdt::One, Zdt::Two, Zdt::Three] {
563 for f1 in [0.0, 0.01, 0.25, 0.42, 0.65, 0.83] {
564 if p.pieces().iter().any(|&(lo, hi)| (lo..=hi).contains(&f1)) {
565 assert!(p.distance_to_front(&[f1, p.front_f2(f1)]) <= 1e-15);
566 }
567 }
568 }
569 // ZDT2 at f₁ = 0.5: the curve's normal is (2 f₁, 1)/√(1 + 4 f₁²) = (1, 1)/√2.
570 let d = 0.01;
571 let out = [0.5 + d / 2f64.sqrt(), 0.75 + d / 2f64.sqrt()];
572 assert!((Zdt::Two.distance_to_front(&out) - d).abs() <= 1e-12);
573 let points = [
574 [0.3, 0.5],
575 [0.001, 1.2],
576 [0.12, 0.9],
577 [0.3, 0.3],
578 [0.55, -0.1],
579 [0.9, -0.5],
580 [1.2, 0.1],
581 [-0.1, 1.1],
582 ];
583 for p in [Zdt::One, Zdt::Two, Zdt::Three] {
584 for f in points {
585 let mut scan = f64::INFINITY;
586 for &(lo, hi) in p.pieces() {
587 for i in 0..=1_000_000 {
588 let f1 = lo + (hi - lo) * f64::from(i) / 1e6;
589 scan = scan.min((f[0] - f1).hypot(f[1] - p.front_f2(f1)));
590 }
591 }
592 let d = p.distance_to_front(&f);
593 assert!(
594 d <= scan + 1e-15,
595 "{p:?} {f:?}: {d} above the scan's {scan}"
596 );
597 assert!(
598 scan - d <= 1e-6,
599 "{p:?} {f:?}: {d} far below the scan's {scan}"
600 );
601 }
602 }
603 }
604
605 #[test]
606 fn reference_points_lie_on_the_front() {
607 let r = Zdt::One.reference(1001);
608 assert_eq!(r.len(), 1001);
609 assert_eq!(r[500], vec![0.5, 1.0 - 0.5f64.sqrt()]);
610 let r3 = Zdt::Three.reference(1001);
611 assert!(r3.iter().all(|f| Zdt::Three.distance_to_front(f) <= 1e-15));
612 // The pieces' f₁ lengths sum to 0.2657, so about a quarter of the points.
613 assert_eq!(r3.len(), 265);
614 }
615 }
616}