1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
//! Worked example: `BoundedCmaInject` with Levenberg-Marquardt inner.
//!
//! Booth-as-residuals on a tight `[-1, 1]²` box — the unconstrained
//! minimum `(1, 3)` is outside the box, so the bound-active constrained
//! optimum sits on the box corner `(1, 1)`. Bounded CMA-ES does global
//! exploration under the adaptive BoundPenalty (Hansen / pycma); each
//! generation's best `k` candidates are polished by an unconstrained LM
//! inner working off the same residual + Jacobian; the refined points
//! are Mahalanobis-clipped (Hansen 2011 eq. 4) and injected back into
//! the population.
//!
//! Run: `cargo test --test example_bounded_cma_inject_lm --features
//! nalgebra -- --nocapture`. The `--nocapture` flag is what lets the
//! progress prints land in your terminal.
#![cfg(feature = "nalgebra")]
use basin::problems::BoothBoxedResiduals;
use basin::{BasicPopulationState, BoundedCmaEs, BoundedCmaInject, Executor, LevenbergMarquardt};
use nalgebra::{DMatrix, DVector};
#[test]
fn example_bounded_cma_inject_lm_on_booth_corner() {
// -----------------------------------------------------------------
// 1. Problem: Booth as residuals (so LM has a Jacobian to chew on)
// with bounds `[-1, 1]²` that exclude the unconstrained min
// `(1, 3)`. The corner `(1, 1)` is the constrained optimum.
// -----------------------------------------------------------------
let lower = DVector::from_vec(vec![-1.0, -1.0]);
let upper = DVector::from_vec(vec![1.0, 1.0]);
let problem = BoothBoxedResiduals::<DVector<f64>>::new(lower, upper);
// -----------------------------------------------------------------
// 2. Outer: bounded CMA-ES with default population size.
// Start at the opposite corner so CMA has to drift across the
// box; σ = 0.3 keeps the initial distribution inside.
// -----------------------------------------------------------------
let m0 = DVector::from_vec(vec![-0.5, -0.5]);
let lambda = BoundedCmaEs::<DVector<f64>, DMatrix<f64>>::default_lambda(2);
let cma = BoundedCmaEs::<DVector<f64>, DMatrix<f64>>::new(m0, 0.3, 42);
// -----------------------------------------------------------------
// 3. Memetic wrapper: top `k = 1` candidate per generation gets
// polished by LM. LM's tol_grad = 1e-8 default terminates the
// inner cleanly once ‖Jᵀr‖_∞ ≤ 1e-8 (2–3 iters on this
// quadratic), so we don't need to cap inner iters explicitly.
// -----------------------------------------------------------------
let solver = BoundedCmaInject::with_inner_solver(cma, LevenbergMarquardt::new()).with_k(1);
// -----------------------------------------------------------------
// 4. Drive.
// -----------------------------------------------------------------
let state = BasicPopulationState::<DVector<f64>>::with_size(lambda);
let result = Executor::new(problem, solver, state).max_iter(100).run();
let p = result.param();
eprintln!();
eprintln!("=== BoundedCmaInject + Levenberg-Marquardt on BoothBoxedResiduals ===");
eprintln!("bounds: [-1, 1]²");
eprintln!("unconstrained min: (1, 3) (outside the box)");
eprintln!("constrained min: (1, 1) (corner of the box)");
eprintln!();
eprintln!("final iterate: ({:.10}, {:.10})", p[0], p[1]);
eprintln!("final cost: {:e}", result.cost());
eprintln!("outer iters: {}", result.iter());
eprintln!(
"cost evals: {} (rolls up LM residual + Jacobian calls)",
result.cost_evals()
);
eprintln!("termination: {:?}", result.reason);
eprintln!();
// -----------------------------------------------------------------
// 5. Sanity check: we should be at the corner (1, 1) within the
// BoundPenalty's tolerance. Inner LM doesn't see the bounds (it's
// unconstrained), so the polished iterates may temporarily step
// toward the unconstrained min (1, 3); BoundPenalty repairs them
// back into feasibility on the next generation.
// -----------------------------------------------------------------
let err = (p[0] - 1.0).abs().max((p[1] - 1.0).abs());
assert!(
err <= 5e-3,
"expected ≈ (1, 1) within 5e-3, got ({}, {}) — err = {}",
p[0],
p[1],
err
);
}