Shekel's foxholes
The problem
Shekel's foxholes, the fifth function of De Jong's test bed, is a function of two variables to minimize:
1 / f(x) = 1/500 + Σⱼ₌₁²⁵ 1 / (j + (x₁ − a₁ⱼ)⁶ + (x₂ − a₂ⱼ)⁶), x₁, x₂ in [−65.536, 65.536]
where the 25 points aⱼ are the grid (−32, −16, 0, 16, 32)², x₁ varying first: a₁ = (−32, −32), a₂ = (−16, −32), …, a₂₅ = (32, 32).
It comes from De Jong's thesis (1975, appendix A.6), read in the scan that George Mason University's EC lab published, where it is test function F5, "synthesized as suggested by Shekel (1971)", with cⱼ = j, K = 500, this grid and these bounds; De Jong gives the minimum as ≅ 1. Yao, Liu and Lin (1999, f14) restate it the same way. Shekel's own paper couldn't be read.
The best known minimum is 0.9980038377944502 at (−31.97833483565697, −31.978334837300796), in the
first hole: Newton's method, computed to 40 digits and rounded. The other holes pull it a little
towards the middle of the box, and the value at (−32, −32) is 0.9980038388186489. The other holes'
minima are about their j: 1.99203, 2.98211, 3.96825, 4.95049, … up to 23.8 for the last. It's not
proven global, and genoxide's problems::ShekelFoxholes gives it as a best known value.
What makes it hard
Away from the holes, all 25 terms are tiny and f is nearly 500: half of the box is above 499.95. Each hole is a flat-bottomed well, since the sixth powers are nearly 0 within about 1 of its center and very large beyond: a random point falls below 10, in one of the ten deepest holes, with a probability of 0.36%, and below 2, in the deepest, with 0.026%. Between the holes, the plane gives no slope that points to a deeper one.
So a search that finds a hole and settles in it is stuck: to find a deeper one, it has to sample the plane again, 16 away. And the holes are many and of all depths: the deepest is only a little deeper than the next, 0.998 against 1.992.
Representation
A Real genome of 2 genes, each in [−65.536, 65.536]: the point (x₁, x₂) itself. The fitness is f,
to minimize. The function, its bounds and its best known minimum are genoxide's
problems::ShekelFoxholes.
Algorithm
Four algorithms, each from seeds 1 to 30, each run stopping once its value is within 1e-8 of the best known minimum, or after 10,000 evaluations:
- CMA-ES (Hansen and Ostermeier, 2001, Evolutionary Computation 9(2): 159-195), with genoxide's defaults: a population of 6 and a step size of 0.3 of each gene's range, 39, from a random start;
- CMA-ES with IPOP restarts (Auger and Hansen, 2005, IEEE CEC 2005: 1769-1776): a run that has converged starts again from a random point with twice the population;
- particle swarm optimization (Kennedy and Eberhart, 1995), 40 particles with Clerc and Kennedy's constriction coefficients and a global topology;
- a real-coded genetic algorithm (De Jong's were binary): a population of 50, tournaments of 3, simulated binary crossover (Deb and Agrawal, 1995) with η = 15 and polynomial mutation with η = 20 at a rate of 1/2 per gene.
Output
The first line gives the best known minimum, the seeds and the budget. Then a row per algorithm: how
many of the 30 runs reach the best known minimum and how many don't, and the median and largest
number of evaluations of the runs that reach it. In Python, run evaluates the function in Rust, so
both versions print the same table.
The page's plot shows the 30 runs of CMA-ES without restarts, each at its best point so far, over the function's contour, with the best known minimum marked: they end in holes all over the grid. A curve gives the best and the median run's error, on a logarithmic axis.
The project page plays this run back.
Good results
A good result reaches the best known minimum in every run. The genetic algorithm does, after a median of 2,049 evaluations and at most 5,081, and so does PSO, after a median of 2,140 and at most 4,120: their populations of 50 and 40 sample the plane widely before they gather.
CMA-ES never does without restarts: the best points of its runs lie in 17 different holes, one of them the deepest (seed 27, at 0.99894, on its side), 4 in the third, at (0, −32), 4 in the one in the middle, at (0, 0), and the others elsewhere. Its population of 6 narrows down on a hole that its first samples find. With IPOP restarts, only one run reaches the minimum, after 5,280 evaluations: the restarts sample the plane again, and 18 runs have a best point in the deepest hole, down to 0.9980038389, but aren't within 1e-8 of its bottom when the budget ends.
Known optimum: 0.99800 at (−31.97833, −31.97833) (best known)
Source: examples/shekel_foxholes
Interactive run: tachsin.gr/projects/genoxide/examples/shekel-foxholes
cargo run --release --example shekel_foxholes
//! Shekel's foxholes: minimize De Jong's fifth function, a plane with 25 narrow holes of different
//! depths, with CMA-ES without and with IPOP restarts, particle swarm optimization and a genetic
//! algorithm, from 30 seeds each.
//!
//! The table counts the runs that reach the deepest hole's minimum, to within 1e-8, and the
//! evaluations they take. The function, its bounds and its minimum come from genoxide's
//! `problems::ShekelFoxholes`.
//!
//! With `GENOXIDE_TRACE=<file>`, it also writes a trace of its runs for the plot on the example's
//! page, with `trace.rs`.
//!
//! ```text
//! cargo run --release --example shekel_foxholes
//! ```
mod trace;
use genoxide::observer::Snapshot;
use genoxide::prelude::*;
use genoxide::problems::{Problem, ShekelFoxholes};
const SEEDS: u64 = 30;
const BUDGET: u64 = 10_000;
// a run stops once its error to the best known minimum is at most this
const ERROR: f64 = 1e-8;
// the algorithms of the table, in its order
const ALGORITHMS: [Algorithm; 4] = [
Algorithm::Cmaes,
Algorithm::CmaesIpop,
Algorithm::Pso,
Algorithm::Ga,
];
// the algorithm whose runs the trace records
const TRACED: Algorithm = Algorithm::Cmaes;
#[derive(Clone, Copy, PartialEq, Eq)]
enum Algorithm {
Cmaes,
CmaesIpop,
Pso,
Ga,
}
impl Algorithm {
fn name(self) -> &'static str {
match self {
Algorithm::Cmaes => "CMA-ES",
Algorithm::CmaesIpop => "CMA-ES with IPOP",
Algorithm::Pso => "PSO",
Algorithm::Ga => "GA",
}
}
}
fn main() -> Result<()> {
let problem = ShekelFoxholes;
let target = problem.optimum().expect("known").value() + ERROR;
println!(
"Shekel's foxholes: best known minimum 0.99800 near (-32, -32), {SEEDS} seeds, {BUDGET} \
evaluations at most per run"
);
println!("runs at min elsewhere evaluations: median largest");
// with GENOXIDE_TRACE=<file>, a trace of the runs of CMA-ES for the plot on the example's page
let mut trace = trace::Trace::from_env();
for algorithm in ALGORITHMS {
let (mut reached, mut elsewhere) = (0, 0);
// the evaluations of the runs that reach the target
let mut evaluations = Vec::new();
for seed in 1..=SEEDS {
let traced = algorithm == TRACED;
let outcome = run(algorithm, seed, target, |snapshot| {
if traced {
trace.record(snapshot);
}
})?;
if outcome.stop_reason() == StopReason::Target {
reached += 1;
evaluations.push(outcome.evaluations());
} else {
elsewhere += 1;
}
}
let largest = evaluations.iter().max().copied().unwrap_or(0);
println!(
"{:<16} {reached:>6} {elsewhere:>9} {:>19} {largest:>7}",
algorithm.name(),
median(&mut evaluations)
);
}
println!("evaluations: of the runs that reach the best known minimum, to within 1e-8");
trace.write();
Ok(())
}
// a run of `algorithm` from `seed`, until its best is at most `target` or it has used BUDGET
// evaluations, which calls `record` after each generation
fn run(
algorithm: Algorithm,
seed: u64,
target: f64,
record: impl FnMut(&Snapshot<'_, Reals>),
) -> Result<Outcome<Reals>> {
let problem = ShekelFoxholes;
let real = problem.representation();
let stop = Stop::target(target).or(Stop::evaluations(BUDGET));
match algorithm {
Algorithm::Cmaes => {
let cmaes = Cmaes::builder(real).minimize().seed(seed).build()?;
Engine::new(cmaes, problem)
.stop_when(stop)
.on_generation(record)
.run()
}
Algorithm::CmaesIpop => {
let cmaes = Cmaes::builder(real)
.restarts(cmaes::Restarts::Ipop)
.minimize()
.seed(seed)
.build()?;
Engine::new(cmaes, problem)
.stop_when(stop)
.on_generation(record)
.run()
}
Algorithm::Pso => {
let pso = Pso::builder(real)
.population_size(40)
.minimize()
.seed(seed)
.build()?;
Engine::new(pso, problem)
.stop_when(stop)
.on_generation(record)
.run()
}
Algorithm::Ga => {
let rate = 1.0 / real.bounds().len() as f64;
let ga = Ga::builder(real)
.population_size(50)
.select(Tournament::new(3)?)
.crossover(SimulatedBinaryCrossover::new(15.0)?)
.mutate(PolynomialMutation::per_gene(rate, 20.0)?)
.minimize()
.seed(seed)
.build()?;
Engine::new(ga, problem)
.stop_when(stop)
.on_generation(record)
.run()
}
}
}
// the median of the evaluations, rounded down; 0 without any
fn median(evaluations: &mut [u64]) -> u64 {
evaluations.sort_unstable();
let middle = evaluations.len() / 2;
match evaluations.len() {
0 => 0,
n if n % 2 == 1 => evaluations[middle],
_ => (evaluations[middle - 1] + evaluations[middle]) / 2,
}
}
python examples/shekel_foxholes/main.py
"""Shekel's foxholes: minimize De Jong's fifth function, a plane with 25 narrow holes of different
depths, with CMA-ES without and with IPOP restarts, particle swarm optimization and a genetic
algorithm, from 30 seeds each.
The table counts the runs that reach the deepest hole's minimum, to within 1e-8, and the evaluations
they take. The function, its bounds and its minimum come from genoxide's `problems::ShekelFoxholes`.
With ``GENOXIDE_TRACE=<file>``, it also writes a trace of its runs for the plot on the example's
page, with trace.py.
python examples/shekel_foxholes/main.py
"""
import genoxide as gx
from trace import Trace
SEEDS = 30
BUDGET = 10_000
# a run stops once its error to the best known minimum is at most this
ERROR = 1e-8
# the algorithms of the table, in its order
ALGORITHMS = ["CMA-ES", "CMA-ES with IPOP", "PSO", "GA"]
# the algorithm whose runs the trace records
TRACED = "CMA-ES"
def median(evaluations):
"""The median of the evaluations, rounded down; 0 without any."""
evaluations = sorted(evaluations)
middle = len(evaluations) // 2
if not evaluations:
return 0
if len(evaluations) % 2:
return evaluations[middle]
return (evaluations[middle - 1] + evaluations[middle]) // 2
def build(name, genome, seed):
"""The algorithm called ``name``, on ``genome``, from ``seed``."""
algorithms = {
"CMA-ES": lambda: gx.Cmaes(genome, objective="minimize", seed=seed),
"CMA-ES with IPOP": lambda: gx.Cmaes(
genome, restarts="ipop", objective="minimize", seed=seed
),
"PSO": lambda: gx.Pso(genome, population_size=40, objective="minimize", seed=seed),
"GA": lambda: gx.Ga(
genome,
population_size=50,
select=gx.Tournament(3),
crossover=gx.SimulatedBinaryCrossover(15.0),
mutation=gx.PolynomialMutation(20.0, rate=1 / 2),
objective="minimize",
seed=seed,
),
}
return algorithms[name]()
problem = gx.problems.ShekelFoxholes()
target = problem.optimum.value + ERROR
print(
f"Shekel's foxholes: best known minimum 0.99800 near (-32, -32), {SEEDS} seeds, {BUDGET} "
"evaluations at most per run"
)
print("runs at min elsewhere evaluations: median largest")
# with GENOXIDE_TRACE=<file>, a trace of the runs of CMA-ES for the plot on the example's page
trace = Trace(problem)
for name in ALGORITHMS:
reached, elsewhere = 0, 0
# the evaluations of the runs that reach the target
evaluations = []
for seed in range(1, SEEDS + 1):
traced = name == TRACED
result = build(name, problem.genome, seed).run(
problem,
target=target,
evaluations=BUDGET,
on_generation=trace.on_generation if traced else None,
)
if result.stop_reason == "target":
reached += 1
evaluations.append(result.evaluations)
else:
elsewhere += 1
largest = max(evaluations, default=0)
print(f"{name:<16} {reached:>6} {elsewhere:>9} {median(evaluations):>19} {largest:>7}")
print("evaluations: of the runs that reach the best known minimum, to within 1e-8")
trace.write()
What it prints, from a seeded run:
Shekel's foxholes: best known minimum 0.99800 near (-32, -32), 30 seeds, 10000 evaluations at most per run
runs at min elsewhere evaluations: median largest
CMA-ES 0 30 0 0
CMA-ES with IPOP 1 29 5280 5280
PSO 30 0 2140 4120
GA 30 0 2049 5081
evaluations: of the runs that reach the best known minimum, to within 1e-8