diff --git a/AGENTS.md b/AGENTS.md index f32b0f0b..8c61512c 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -136,10 +136,10 @@ Stops: `Stop::target(score)` (at least as good), `generations(n)`, `evaluations( - Maximize is the default; use `.minimize()`, don't negate. - `None`, `Fitness::invalid()` and NaN are invalid: worse than everything. - **Constraints:** return `(score, violation)`, 0 when feasible, adding up `constraint::at_most(value, limit)`, `at_least`, `equal(value, target, tolerance)`. Deb's rules: feasible beats infeasible, then score or violation decides. Select with `Tournament` or `Rank`: roulette and SUS give infeasible solutions no weight. `Penalty::new(weight)?.fitness(objective, score, violation)` is a static penalty instead. -- **Test problems:** `problems::{Sphere, AxisParallelEllipsoid, Schwefel1_2, Rastrigin, Rosenbrock, Ackley, Griewank, Schwefel2_26, Levy, Zakharov, StyblinskiTang, Michalewicz, Schwefel2_21, Schwefel2_22, DixonPrice, Trid, Powell, SumOfDifferentPowers, Step, Quartic, Penalized1, Penalized2, HighConditionedElliptic, BentCigar, Discus, DifferentPowers, BucheRastrigin, NonContinuousRastrigin, Weierstrass, Katsuura, HappyCat, HgBat, SchafferF7, RotatedHyperEllipsoid}::new(n)` (`Powell` takes a multiple of 4; `Quartic::noisy(n)` adds noise drawn from the genome) and `problems::{Himmelblau, Branin, GoldsteinPrice, SixHumpCamel, Hartmann3, Hartmann6, Shekel5, Shekel7, Shekel10, Easom, Eggholder, SchafferF6, Beale, Booth, Matyas, Bohachevsky1, Bohachevsky2, Bohachevsky3, ThreeHumpCamel, Langermann, ShekelFoxholes, Kowalik}` are fitness functions for `Engine::new(algorithm, problem)`, all minimized. The `problems::Problem` trait gives `representation()` (the bounds), `optimum()` (`value()`, `solutions()`), `reference()`; `problems::all()` lists them as `Box`. `problems::Shifted::new(problem, seed)` and `problems::Rotated::new(problem, seed)` make CEC/BBOB-style instances of any of them (e.g. CEC 2005's F10: `Rotated::new(Shifted::new(Rastrigin::new(n), seed), seed)`), keeping the optimum's value. Constrained, with fitness `(score, violation)` and `constraints(&x)` (`g <= 0`, then `h = 0`): `problems::cec2006::{G01, …, G24}` (equalities met within `EQUALITY_TOLERANCE` = 1e-4; `with_tolerance(δ)` for the problems with equalities, e.g. `G03::with_tolerance(δ)`), and `problems::engineering::{WeldedBeam, WeldedBeamRagsdell, PressureVessel, TensionCompressionSpring, SpeedReducer, ThreeBarTruss, CantileverBeam, CarSideImpact}`. `PressureVessel` and `SpeedReducer` round their discrete genes when evaluated; `design(&x)` gives the rounded design. `engineering::GearTrain` has an `Integer` genome and isn't in `all()`. `Optimum::is_proven()` is false for a best known value. +- **Test problems:** `problems::{Sphere, AxisParallelEllipsoid, Schwefel1_2, Rastrigin, Rosenbrock, Ackley, Griewank, Schwefel2_26, Levy, Zakharov, StyblinskiTang, Michalewicz, Schwefel2_21, Schwefel2_22, DixonPrice, Trid, Powell, SumOfDifferentPowers, Step, Quartic, Penalized1, Penalized2, HighConditionedElliptic, BentCigar, Discus, DifferentPowers, BucheRastrigin, NonContinuousRastrigin, Weierstrass, Katsuura, HappyCat, HgBat, SchafferF7, RotatedHyperEllipsoid}::new(n)` (`Powell` takes a multiple of 4; `Quartic::noisy(n)` adds noise drawn from the genome) and `problems::{Himmelblau, Branin, GoldsteinPrice, SixHumpCamel, Hartmann3, Hartmann6, Shekel5, Shekel7, Shekel10, Easom, Eggholder, SchafferF6, Beale, Booth, Matyas, Bohachevsky1, Bohachevsky2, Bohachevsky3, ThreeHumpCamel, Langermann, ShekelFoxholes, Kowalik}` are fitness functions for `Engine::new(algorithm, problem)`, all minimized. The `problems::Problem` trait gives `representation()` (the bounds), `optimum()` (`value()`, `solutions()`), `reference()`; `problems::all()` lists them as `Box`. `problems::Shifted::new(problem, seed)` and `problems::Rotated::new(problem, seed)` make CEC/BBOB-style instances of any of them (e.g. CEC 2005's F10: `Rotated::new(Shifted::new(Rastrigin::new(n), seed), seed)`), keeping the optimum's value and what the wrapped problem provides: its gradient (by the chain rule through the rotation), constraint values and their Jacobian. Constrained, with fitness `(score, violation)` and `constraints(&x)` (`g <= 0`, then `h = 0`): `problems::cec2006::{G01, …, G24}` (equalities met within `EQUALITY_TOLERANCE` = 1e-4; `with_tolerance(δ)` for the problems with equalities, e.g. `G03::with_tolerance(δ)`), and `problems::engineering::{WeldedBeam, WeldedBeamRagsdell, PressureVessel, TensionCompressionSpring, SpeedReducer, ThreeBarTruss, CantileverBeam, CarSideImpact}`. `PressureVessel` and `SpeedReducer` round their discrete genes when evaluated; `design(&x)` gives the rounded design. `engineering::GearTrain` has an `Integer` genome and isn't in `all()`. `Optimum::is_proven()` is false for a best known value. - **Extras:** return `Evaluated::new(value, info)` (`value` any of the above, `info` any `Send + Sync + 'static` type, e.g. a struct with a penalty's terms) to keep what the fitness function computed. Read it by type: `outcome.best_info::()`, `snapshot.info::(genome)` and `snapshot.best_info::()` in `.on_generation`, `hall_of_fame.info::(genome)`, and in `MultiEngine` `snapshot.info` and `outcome.info(genome)` for the front; `None` for another type. Never used by the search. Kept by genome for the population, the discarded and the best (copies share it); not in checkpoints. - **Gradients:** `Differentiable(|x: &Reals, gradient: &mut [f64]| value)`, for gradient-based methods; see [Gradients](#gradients-supplying-them). -- **Constraint values:** `Constrained::new(m, |x: &Reals, g: &mut [f64]| score)` writes the values of m constraints `gᵢ(x) <= 0`; `Constrained::differentiable(m, |x, gradient, g, jacobian| score)` also the gradient and the Jacobian (`jacobian[i * n + j]` = ∂gᵢ/∂xⱼ). Either is `(score, Σ max(0, gᵢ))` for any algorithm, and gives the values one by one to those that use them (`Mma`). The CEC 2006 problems with inequalities only and the engineering problems give their values (`problem.provides().inequalities`), not their gradients. +- **Constraint values:** `Constrained::new(m, |x: &Reals, g: &mut [f64]| score)` writes the values of m constraints `gᵢ(x) <= 0`; `Constrained::differentiable(m, |x, gradient, g, jacobian| score)` also the gradient and the Jacobian (`jacobian[i * n + j]` = ∂gᵢ/∂xⱼ). Either is `(score, Σ max(0, gᵢ))` for any algorithm, and gives the values one by one to those that use them (`Mma`). The CEC 2006 problems with inequalities only and the engineering problems give their values (`problem.provides().inequalities`), shifted or rotated too, not their gradients. - **Batch:** `Batch(|genomes: &[&G]| -> Vec)` scores a generation in one call, in order (SIMD, GPU, remote), in `Engine` or `MultiEngine`; the slice can be empty. A wrong count is `Error::FitnessCount`. See `examples/gpu` (wgpu). ## Templates @@ -909,7 +909,7 @@ How fitness functions give gradients to gradient-based methods ([L-BFGS-B](#l-bf - `Differentiable(|x: &Reals, gradient: &mut [f64]| value)` writes the gradient (zeroed, one value per gene) and returns the value; any algorithm takes it as a plain fitness function. `Batch(Differentiable(|xs: &[&Reals], gradients: &mut [f64]| values))`: flat, row-major, a row per genome. - A `FitnessFunction` declares it with `fn provides(&self) -> Provided { Provided::GRADIENT }` and writes it in `fn evaluate_with(&self, x, extras: &mut Extras<'_>)` when `extras.gradient()` is `Some` (`genoxide::engine::{Extras, Provided}`); the value must be `evaluate`'s, to the bit. -- The smooth test problems supply theirs: every classic function in `problems` but `Eggholder`, `Schwefel2_21` and `Schwefel2_22`, and not yet the CEC- and BBOB-style ones (`SumOfDifferentPowers` to `RotatedHyperEllipsoid`) or the `Shifted` and `Rotated` wrappers; `problem.provides().gradient`. +- The test problems supply theirs: every classic function in `problems` but `Eggholder`, `Schwefel2_21`, `Schwefel2_22`, `Step`, `NonContinuousRastrigin`, `Katsuura` and `Quartic::noisy`, whose derivative is undefined or 0 on sets of positive measure (at an isolated cone or cusp, such as `Ackley`'s or `HappyCat`'s, the kinked term's gradient is 0); `Shifted` and `Rotated` pass the wrapped problem's on (`Mᵀ ∇f` for a rotation); `problem.provides().gradient`. - `gradient::check(&function, &x)?` compares a supplied gradient with central differences: `.largest()` about 1e-10 when right, `.worst_gene()`. - An algorithm's `gradient::Gradients` setting: `Auto` (default: supplied if provided, else forward differences, n evaluations per gradient, up to `gradient::AUTO_LIMIT` = 10⁴ genes), `Supplied` (an error at the start of a run without one), `Forward { step: None }`, `Central { step: None }` (2n per gradient, more accurate). Finite differences count towards `Stop::evaluations`. - A NaN in a gradient follows the `NanPolicy`: invalid fitness, or `Error::NanFitness`. diff --git a/docs/features.md b/docs/features.md index edcca1d5..03fbbc5e 100644 --- a/docs/features.md +++ b/docs/features.md @@ -73,7 +73,7 @@ What genoxide has on main; [docs.rs](https://docs.rs/genoxide) documents the lat - **Constraints:** the values of `Constrained` fitness functions and of the constrained test problems, a Gaussian process per constraint, and the acquisition weighed by the probability of feasibility (Gardner et al. 2014); before a feasible point, that probability alone. - **Integer genomes:** the genes rounded inside the kernel (Garrido-Merchán and Hernández-Lobato 2020), the acquisition maximized on the lattice by a hill climb (every point of a small lattice at once), the run finished once the lattice is evaluated. - **Gaussian processes** (`model::gp`, unstable for one release): regression with a constant mean, an ARD Matérn 5/2 or squared exponential kernel, no noise by default (the model interpolates a deterministic function) or learned noise, inputs scaled to the unit cube and outputs standardized; hyperparameters by maximum marginal likelihood (Rasmussen and Williams, eq. 2.30 and 5.9) with genoxide's L-BFGS-B from fixed and seeded random starts; the posterior mean and variance with their gradients; Cholesky factorizations with a growing jitter. In Python as `gx.model.gp`. -- **Test problems** (`problems`): Sphere, the axis-parallel ellipsoid, Schwefel 1.2 and 2.26, Rastrigin, Rosenbrock, Ackley, Griewank, Levy, Zakharov, Styblinski-Tang, Michalewicz, Himmelblau, Branin, Goldstein-Price, the six-hump camel, Hartmann's functions in 3 and 6 dimensions, Shekel's with 5, 7 and 10 wells, Easom, the eggholder, Schaffer's F6, Schwefel 2.21 and 2.22, Dixon-Price, Trid, Powell's singular function, Beale, Booth, Matyas, Bohachevsky's three functions, the three-hump camel, Langermann, Shekel's foxholes, Kowalik, the sum of different powers, the step function, the quartic (with or without noise), Yao, Liu and Lin's two penalized functions, the high-conditioned elliptic, the bent cigar, the discus, BBOB's different powers, Büche-Rastrigin, the non-continuous Rastrigin, Weierstrass, Katsuura, HappyCat, HGBat, Schaffer's F7 and the rotated hyper-ellipsoid, each with its bounds, known optimum (or best known, for those found numerically) and reference, in Rust and Python, and the analytic gradient of the smooth ones up to Kowalik and Powell (all but the eggholder and Schwefel 2.21 and 2.22, which aren't differentiable; not yet the CEC- and BBOB-style ones); and shift and rotation wrappers, generated from a seed, for CEC- and BBOB-style instances of any of them. +- **Test problems** (`problems`): Sphere, the axis-parallel ellipsoid, Schwefel 1.2 and 2.26, Rastrigin, Rosenbrock, Ackley, Griewank, Levy, Zakharov, Styblinski-Tang, Michalewicz, Himmelblau, Branin, Goldstein-Price, the six-hump camel, Hartmann's functions in 3 and 6 dimensions, Shekel's with 5, 7 and 10 wells, Easom, the eggholder, Schaffer's F6, Schwefel 2.21 and 2.22, Dixon-Price, Trid, Powell's singular function, Beale, Booth, Matyas, Bohachevsky's three functions, the three-hump camel, Langermann, Shekel's foxholes, Kowalik, the sum of different powers, the step function, the quartic (with or without noise), Yao, Liu and Lin's two penalized functions, the high-conditioned elliptic, the bent cigar, the discus, BBOB's different powers, Büche-Rastrigin, the non-continuous Rastrigin, Weierstrass, Katsuura, HappyCat, HGBat, Schaffer's F7 and the rotated hyper-ellipsoid, each with its bounds, known optimum (or best known, for those found numerically) and reference, in Rust and Python, and the analytic gradient of all but seven (the eggholder, Schwefel 2.21 and 2.22, the step function, the non-continuous Rastrigin, Katsuura and the noisy quartic, whose derivative is undefined or 0 on sets of positive measure); and shift and rotation wrappers, generated from a seed, for CEC- and BBOB-style instances of any of them, which pass on the wrapped problem's gradient (by the chain rule), constraint values and Jacobian. - **Constrained test problems:** CEC 2006's g01-g24 (`problems::cec2006`), and the engineering design problems (`problems::engineering`): the welded beam in two forms, the pressure vessel, the tension/compression spring, the speed reducer, the gear train (integer), the three-bar truss, the cantilever beam and the car side impact, each with its optimum or best known solution and references, in Rust and Python. Those with inequalities only give their constraints' values one by one to the algorithms that use them. ## Multi-objective diff --git a/docs/problems-plan.md b/docs/problems-plan.md index 12dec314..72acd0d4 100644 --- a/docs/problems-plan.md +++ b/docs/problems-plan.md @@ -207,7 +207,15 @@ as (1, …, 1) (it's (−1, …, −1)) and its table I drops the square of f13' report prints Schaffer F7 with sin for BBOB's sin²; Katsuura's minima include the bounds ±5, where a PSO that stops particles at the bounds lands on them; at G06's minimum, shifted, rounding in x − o leaves an active constraint violated by 6e-14. The minima are proven from the formulas (every term -at least 0), but the noisy quartic's, which isn't known. +at least 0), but the noisy quartic's, which isn't known. Gradients (#416): all but the step +function, the non-continuous Rastrigin (flat steps), Katsuura (kinks 2⁻³³ apart) and the noisy +quartic (noise that jumps between any two genomes) supply their analytic gradient, 0 for the +kinked term at the cones and cusps of measure 0 (different powers, HappyCat, HGBat, Schaffer F7); +Büche-Rastrigin is differentiable at 0 despite T_osz, its terms being O(x²) there. Weierstrass's +is checked on its partial sums: the computed phase 2π 3ᵏ (x + 0.5) of its highest terms rounds by +2·10⁻⁶ rad, which no differences of the function resolve. `Shifted` passes the wrapped problem's +extras on at `x − o`; `Rotated` its gradient as `Mᵀ ∇f`, its constraint values, and its +Jacobian as `J M`. **Hand-computed test values (VC):** Rosenbrock (−1.2, 1) = 24.2; Powell (3, −1, 0, 1) = 215; Goldstein-Price (0, −1) = 3 and its three local minima's values from the paper; Beale (3, 0.5) = diff --git a/python/genoxide/problems/__init__.py b/python/genoxide/problems/__init__.py index 88b99168..e1374777 100644 --- a/python/genoxide/problems/__init__.py +++ b/python/genoxide/problems/__init__.py @@ -57,6 +57,13 @@ cmaes = gx.Cmaes(problem.genome, restarts="ipop", objective="minimize", seed=1) result = cmaes.run(problem, target=problem.optimum.value + 1e-8, evaluations=200_000) +The gradient methods (:class:`genoxide.Lbfgsb`, :class:`genoxide.FirstOrder`, :class:`genoxide.Mma`) +take a problem's analytic gradient, computed in Rust. Every classic function has one but the +eggholder, Schwefel 2.21 and 2.22, the step function, the non-continuous Rastrigin, Katsuura and +the noisy quartic, whose derivative is undefined or 0 on sets of positive measure. The wrappers +pass the problem's gradient on (turned by the rotation), and a constrained problem's constraint +values, which :class:`genoxide.Bo` models. + All problems here are minimized, on :class:`genoxide.Real` genomes except the gear train's :class:`genoxide.Integer` and :class:`Zdt5`'s :class:`genoxide.Binary`. Each class's docstring gives @@ -1507,7 +1514,8 @@ class Shifted(Problem[Real]): value, its solutions shifted (those that leave the box are dropped). That holds when the problem's minimum is its minimum over all of ℝⁿ, as for the functions that CEC and BBOB shift, not for one whose minimum is only the lowest in its box, such as :class:`Schwefel2_26`. The - same seed gives the same shift as Rust's ``problems::Shifted`` on every platform. + problem's analytic gradient and constraint values pass through, those at ``x − o``. The same + seed gives the same shift as Rust's ``problems::Shifted`` on every platform. """ problem: Problem[Real] @@ -1535,8 +1543,9 @@ class Rotated(Problem[Real]): turns about the optimum, as CEC 2005 and BBOB do, so the optimum stays in place with its value. Rotating a :class:`Shifted` problem gives CEC 2005's shifted rotated functions, such as its F10, ``Rotated(Shifted(Rastrigin(n), seed), seed)``. The bounds, the name and the - reference are the problem's. The same seed gives the same matrix as Rust's - ``problems::Rotated`` on every platform. + reference are the problem's. The problem's analytic gradient passes through by the chain rule, + ``Mᵀ ∇f(c + M (x − c))``, and so do the constraint values at ``c + M (x − c)``. The same seed + gives the same matrix as Rust's ``problems::Rotated`` on every platform. """ problem: Problem[Real] diff --git a/python/src/problems.rs b/python/src/problems.rs index b572fa3b..ee0ebb55 100644 --- a/python/src/problems.rs +++ b/python/src/problems.rs @@ -3,7 +3,7 @@ //! evaluated in Rust. use crate::run::{WithObjectives, with_objectives}; -use genoxide::engine::{FitnessFunction, IntoFitness}; +use genoxide::engine::{Extras, FitnessFunction, IntoFitness, Provided}; use genoxide::genome::{Binary, Bits, Integer, Integers, Real, Reals, Representation}; use genoxide::multi::problems::{self as multi, DynMultiProblem, MultiProblem, try_boxed}; use genoxide::multi::{IntoScores, MultiFitnessFunction, Scores}; @@ -1813,6 +1813,15 @@ impl FitnessFunction for Wrapped { fn evaluate(&self, genome: &Reals) -> Fitness { self.0.evaluate(genome) } + + // the wrapped problem's gradient and constraint values, which the wrappers pass on + fn provides(&self) -> Provided { + self.0.provides() + } + + fn evaluate_with(&self, genome: &Reals, extras: &mut Extras<'_>) -> Fitness { + self.0.evaluate_with(genome, extras) + } } impl problems::Problem for Wrapped { diff --git a/python/tests/test_bo.py b/python/tests/test_bo.py index 606be9f2..c61b92b6 100644 --- a/python/tests/test_bo.py +++ b/python/tests/test_bo.py @@ -355,6 +355,11 @@ def test_wrong_constraints_are_errors(): bo.run(toy, constraints=2, batch=True, evaluations=10) with pytest.raises(ValueError, match="own constraints"): bo.run(gx.problems.cec2006.G24(), constraints=2, evaluations=10) + # and so does a shifted or rotated one: the wrappers pass the values on + g24 = gx.problems.cec2006.G24() + for wrapped in [gx.problems.Shifted(g24, seed=1), gx.problems.Rotated(g24, seed=1)]: + with pytest.raises(ValueError, match="own constraints"): + bo.run(wrapped, constraints=2, evaluations=10) ucb = gx.Bo(gx.Real((0.0, 1.0), length=2), acquisition=gx.UpperConfidenceBound(2.0)) with pytest.raises(ValueError, match="upper confidence bound"): ucb.run(toy, constraints=2, evaluations=10) diff --git a/python/tests/test_lbfgsb.py b/python/tests/test_lbfgsb.py index 10459243..85db85de 100644 --- a/python/tests/test_lbfgsb.py +++ b/python/tests/test_lbfgsb.py @@ -46,6 +46,36 @@ def test_a_problem_converges_with_its_gradient_in_rust(): assert result.evaluations < 1_000 +def test_a_shifted_and_rotated_problem_keeps_its_gradient_in_rust(): + elliptic = gx.problems.HighConditionedElliptic(10) + problem = gx.problems.Rotated(gx.problems.Shifted(elliptic, seed=1), seed=1) + result = gx.Lbfgsb( + problem.genome, + gradients="supplied", + gradient_tolerance=1e-9, + function_tolerance=0.0, + objective="minimize", + seed=1, + ).run(problem, evaluations=10_000) + assert result.stop_reason == "converged" + assert result.best_fitness < 1e-16 + assert np.allclose(result.best_genome, problem.optimum.solutions[0], rtol=0.0, atol=1e-8) + # one evaluation per iteration or line-search trial, the gradient with it + assert result.evaluations < 1_000 + + +def test_a_problem_without_a_gradient_has_none_through_the_wrappers(): + for problem in [ + gx.problems.Shifted(gx.problems.Katsuura(5), seed=1), + gx.problems.Rotated(gx.problems.Step(5), seed=1), + gx.problems.Rotated(gx.problems.Shifted(gx.problems.Quartic(5, noisy=True), seed=1), seed=1), + ]: + with pytest.raises(ValueError, match="gradient"): + gx.Lbfgsb(problem.genome, gradients="supplied", objective="minimize").run( + problem, evaluations=100 + ) + + def test_every_way_to_give_the_gradient_takes_the_same_path(): def both(x): return rosenbrock(x), rosenbrock_gradient(x) diff --git a/src/problems.rs b/src/problems.rs index 8180b8f2..0d4f1b7b 100644 --- a/src/problems.rs +++ b/src/problems.rs @@ -89,14 +89,26 @@ //! `Shifted::new(Rastrigin::new(n), seed)` and `Rotated::new(Shifted::new(Rastrigin::new(n), seed), seed)`, //! with genoxide's own shift and matrix rather than the report's data files. //! -//! The classic functions of the table up to [`Kowalik`] but [`Eggholder`], [`Schwefel2_21`] and -//! [`Schwefel2_22`], which aren't differentiable everywhere, supply their analytic gradient to the -//! algorithms that want one (see [`gradient`](crate::gradient)): [`FitnessFunction::provides`] -//! says so, and [`FitnessFunction::evaluate_with`] computes it, with [`math`](crate::math)'s -//! functions. Ackley's has a cone at the origin, where its gradient is taken as 0, and Schwefel -//! 2.26's second derivative is unbounded at 0. The functions after [`Kowalik`], from -//! [`SumOfDifferentPowers`] on, and the [`Shifted`] and [`Rotated`] wrappers don't supply one -//! yet: algorithms estimate it by finite differences. +//! The classic functions of the table supply their analytic gradient to the algorithms that want +//! one (see [`gradient`](crate::gradient)): [`FitnessFunction::provides`] says so, and +//! [`FitnessFunction::evaluate_with`] computes it, with [`math`](crate::math)'s functions. All do +//! but six, whose derivative is undefined or 0 on sets of positive measure: [`Eggholder`], +//! [`Schwefel2_21`] and [`Schwefel2_22`] (absolute values and a maximum), [`Step`] and +//! [`NonContinuousRastrigin`] (flat steps), and [`Katsuura`] (kinks 2⁻³³ apart), and the noisy +//! [`Quartic`], whose noise jumps between any two genomes; for those, algorithms estimate a +//! gradient by finite differences. Where a function with a gradient has no derivative, on a set +//! of measure 0, the gradient of the term with the kink is taken as 0: at the cones of +//! [`Ackley`] and [`DifferentPowers`] at the origin, the cusps of [`HappyCat`] and [`HgBat`] where +//! their first term is 0, and those of [`SchafferF7`] where a pair of genes is 0. Schwefel 2.26's +//! second derivative is unbounded at 0, and [`Weierstrass`]'s gradient has terms that vary on +//! scales of 10⁻¹⁰. +//! +//! The [`Shifted`] and [`Rotated`] wrappers pass on what the wrapped problem provides: its +//! gradient (at `x − o` for a shift, and `Mᵀ ∇f(c + M (x − c))` for a rotation, by the chain +//! rule), and a constrained problem's constraint values and their Jacobian (the rows times `M` +//! for a rotation): a shifted or rotated [`G06`](cec2006::G06) still gives its constraints' +//! values to [`Bo`](crate::algorithm::Bo), and a problem that gives their Jacobian too, to +//! [`Mma`](crate::algorithm::Mma). //! //! Two submodules hold constrained problems, whose fitness is `(score, violation)`: //! @@ -393,8 +405,8 @@ pub trait DynProblem: Send + Sync { fn constraints(&self, genome: &Reals) -> Constraints; /// What the problem gives besides the fitness, as [`FitnessFunction::provides`]: the - /// gradient, for the smooth classic functions, and the constraints' values, for the - /// problems with inequalities only. Nothing by default. + /// gradient, for the classic functions that are differentiable almost everywhere, and the + /// constraints' values, for the problems with inequalities only. Nothing by default. fn provides(&self) -> Provided { Provided::NOTHING } diff --git a/src/problems/classic.rs b/src/problems/classic.rs index 5340782b..89c77207 100644 --- a/src/problems/classic.rs +++ b/src/problems/classic.rs @@ -2546,6 +2546,9 @@ scalable!( /// /// Bounds [−1, 1]ⁿ; minimum 0 at the origin; 30 dimensions by default. /// + /// Supplies its analytic gradient, `(i + 1) |xᵢ|^i sign(xᵢ)`, through + /// [`evaluate_with`](FitnessFunction::evaluate_with). + /// /// Its origin is unknown: definition and bounds as Molga and Smutnicki (2005, section 2.8) /// give them. Not yet checked against an original /// ([#168](https://github.com/tachsin/genoxide/issues/168)). @@ -2558,6 +2561,8 @@ scalable!( impl FitnessFunction for SumOfDifferentPowers { type Output = f64; + gradient!(gradients::sum_of_different_powers); + fn evaluate(&self, x: &Reals) -> f64 { x.iter() .enumerate() @@ -2580,6 +2585,8 @@ scalable!( /// Bounds [−100, 100]ⁿ; minimum 0 on the whole cube [−0.5, 0.5)ⁿ, here at the origin; 30 /// dimensions by default. /// + /// It supplies no gradient: its derivative is 0 on the steps and undefined at their edges. + /// /// Yao, X., Liu, Y. and Lin, G. (1999). Evolutionary programming made faster. *IEEE /// Transactions on Evolutionary Computation* 3(2): 82-102, function f6 (table I and the /// appendix, read): its definition, bounds and dimension. De Jong's (1975) F3, which it is @@ -2617,12 +2624,18 @@ scalable_problem!( /// Bounds [−1.28, 1.28]ⁿ; minimum 0 at the origin, without noise; 30 dimensions by default. The /// function is flat near the minimum: at 0.01 from it in every gene, it's below 10⁻⁵. /// +/// Without noise, it supplies its analytic gradient, `4 i xᵢ³`, through +/// [`evaluate_with`](FitnessFunction::evaluate_with); with noise, none (see below). +/// /// Yao, Liu and Lin (1999, f7) add a uniform random number in [0, 1) to each evaluation, so that /// an algorithm can't use differences smaller than the noise. A fitness function is deterministic /// in genoxide (a copy of a genome inherits its fitness), so [`noisy`](Quartic::noisy) draws the /// noise from a generator seeded with the genome's bits: the same genome always gets the same /// noise, and two genomes, however close, independent noises. Its minimum isn't known (it's the -/// smallest noise near the origin), so [`optimum`](Problem::optimum) is `None`. +/// smallest noise near the origin), so [`optimum`](Problem::optimum) is `None`. Nor has it a +/// gradient: its value jumps between any two genomes, so it has no derivative anywhere, and the +/// noiseless part's gradient would lead a method to that function's minimum, ignoring the noise +/// that the function is there to test. /// /// De Jong, K. A. (1975). *An Analysis of the Behavior of a Class of Genetic Adaptive Systems.* /// PhD thesis, University of Michigan, function F4, with Gaussian noise, in 30 dimensions on @@ -2713,6 +2726,33 @@ impl FitnessFunction for Quartic { quartic } } + + /// The gradient without noise; nothing with noise, whose value jumps between any two + /// genomes, however close, so that it has no derivative anywhere. + fn provides(&self) -> Provided { + if self.noisy { + Provided::NOTHING + } else { + Provided::GRADIENT + } + } + + /// The value at `x`, as [`evaluate`](FitnessFunction::evaluate), and, without noise, its + /// analytic gradient if it's wanted. + /// + /// # Panics + /// + /// If the gradient doesn't have a value per gene of `x`. + fn evaluate_with(&self, x: &Reals, extras: &mut Extras<'_>) -> f64 { + if !self.noisy + && let Some(gradient) = extras.gradient() + { + assert_eq!(gradient.len(), x.len(), "a gradient has a value per gene"); + gradient.fill(0.0); + gradients::quartic(x, gradient); + } + self.evaluate(x) + } } impl Problem for Quartic { @@ -2767,6 +2807,10 @@ scalable!( /// Bounds [−50, 50]ⁿ; minimum 0 at (−1, …, −1), where every yᵢ is 1 (every term is at least /// 0); 30 dimensions by default. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with): + /// the penalty's slope is 0 at ±10, where it starts, so the function is differentiable + /// everywhere. + /// /// Yao, X., Liu, Y. and Lin, G. (1999). Evolutionary programming made faster. *IEEE /// Transactions on Evolutionary Computation* 3(2): 82-102, function f12 (table I and the /// appendix, read), whose appendix misprints the minimizer as (1, …, 1). The function is @@ -2780,6 +2824,8 @@ scalable!( impl FitnessFunction for Penalized1 { type Output = f64; + gradient!(gradients::penalized_1); + fn evaluate(&self, x: &Reals) -> f64 { let n = x.len(); let y: Vec = x.iter().map(|xi| 1.0 + (xi + 1.0) / 4.0).collect(); @@ -2813,6 +2859,10 @@ scalable!( /// Bounds [−50, 50]ⁿ; minimum 0 at (1, …, 1) (every term is at least 0); 30 dimensions by /// default. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with): + /// the penalty's slope is 0 at ±5, where it starts, so the function is differentiable + /// everywhere. + /// /// Yao, X., Liu, Y. and Lin, G. (1999). Evolutionary programming made faster. *IEEE /// Transactions on Evolutionary Computation* 3(2): 82-102, function f13 (table I and the /// appendix, read). Table I prints the last term's `(xₙ − 1)` without its square, which the @@ -2827,6 +2877,8 @@ scalable!( impl FitnessFunction for Penalized2 { type Output = f64; + gradient!(gradients::penalized_2); + fn evaluate(&self, x: &Reals) -> f64 { let (Some(&first), Some(&last)) = (x.first(), x.last()) else { return 0.0; @@ -2851,7 +2903,7 @@ scalable_problem!( // 10⁶ to the power (i − 1) / (n − 1), for gene i from 0: from 1 for the first gene to 10⁶ for // the last -fn conditioning(i: usize, n: usize) -> f64 { +pub(super) fn conditioning(i: usize, n: usize) -> f64 { math::powf(1e6, i as f64 / (n - 1) as f64) } @@ -2861,6 +2913,8 @@ scalable!( /// /// Bounds [−100, 100]ⁿ; minimum 0 at the origin; at least 2 dimensions, 30 by default. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with). + /// /// Suganthan, P. N., Hansen, N., Liang, J. J., Deb, K., Chen, Y.-P., Auger, A. and Tiwari, S. /// (2005). *Problem Definitions and Evaluation Criteria for the CEC 2005 Special Session on /// Real-Parameter Optimization*, function F3 (read), shifted and rotated there: its @@ -2877,6 +2931,8 @@ scalable!( impl FitnessFunction for HighConditionedElliptic { type Output = f64; + gradient!(gradients::high_conditioned_elliptic); + fn evaluate(&self, x: &Reals) -> f64 { let n = x.len(); x.iter() @@ -2900,6 +2956,8 @@ scalable!( /// /// Bounds [−100, 100]ⁿ; minimum 0 at the origin; at least 2 dimensions, 30 by default. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with). + /// /// Hansen, N., Finck, S., Ros, R. and Auger, A. (2009). *Real-Parameter Black-Box /// Optimization Benchmarking 2009: Noiseless Functions Definitions.* INRIA research report /// RR-6829, function f12 (read), which composes it with an asymmetric transformation and two @@ -2916,6 +2974,8 @@ scalable!( impl FitnessFunction for BentCigar { type Output = f64; + gradient!(gradients::bent_cigar); + fn evaluate(&self, x: &Reals) -> f64 { let Some((first, rest)) = x.split_first() else { return 0.0; @@ -2938,6 +2998,8 @@ scalable!( /// /// Bounds [−100, 100]ⁿ; minimum 0 at the origin; at least 2 dimensions, 30 by default. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with). + /// /// Hansen, N., Finck, S., Ros, R. and Auger, A. (2009). *Real-Parameter Black-Box /// Optimization Benchmarking 2009: Noiseless Functions Definitions.* INRIA research report /// RR-6829, function f11 (read), which composes it with an oscillation and a rotation. This @@ -2953,6 +3015,8 @@ scalable!( impl FitnessFunction for Discus { type Output = f64; + gradient!(gradients::discus); + fn evaluate(&self, x: &Reals) -> f64 { let Some((first, rest)) = x.split_first() else { return 0.0; @@ -2975,6 +3039,10 @@ scalable!( /// /// Bounds [−5, 5]ⁿ; minimum 0 at the origin; at least 2 dimensions, 30 by default. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with). + /// The square root makes a cone at the origin, where the gradient is taken as 0, as + /// [`Ackley`]'s. + /// /// Hansen, N., Finck, S., Ros, R. and Auger, A. (2009). *Real-Parameter Black-Box /// Optimization Benchmarking 2009: Noiseless Functions Definitions.* INRIA research report /// RR-6829, function f14 (read): its definition, without the rotation, and its search @@ -2988,6 +3056,8 @@ scalable!( impl FitnessFunction for DifferentPowers { type Output = f64; + gradient!(gradients::different_powers); + fn evaluate(&self, x: &Reals) -> f64 { let n = x.len(); let exponent = |i: usize| 2.0 + 4.0 * i as f64 / (n.max(2) - 1) as f64; @@ -3009,7 +3079,7 @@ scalable_problem!( // BBOB's oscillation T_osz of one value: the identity, but for small smooth wiggles that scale // with the value -fn oscillation(x: f64) -> f64 { +pub(super) fn oscillation(x: f64) -> f64 { if x == 0.0 { return 0.0; } @@ -3031,6 +3101,10 @@ scalable!( /// /// Bounds [−5, 5]ⁿ; minimum 0 at the origin; at least 2 dimensions, 30 by default. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with). + /// T_osz has no derivative at 0, but each gene's term is O(x²) there, so the function's + /// derivative is 0 at 0 and it's differentiable everywhere. + /// /// Hansen, N., Finck, S., Ros, R. and Auger, A. (2009). *Real-Parameter Black-Box /// Optimization Benchmarking 2009: Noiseless Functions Definitions.* INRIA research report /// RR-6829, function f4 (read): its definition, with its optimum at the origin and no offset @@ -3044,6 +3118,8 @@ scalable!( impl FitnessFunction for BucheRastrigin { type Output = f64; + gradient!(gradients::buche_rastrigin); + fn evaluate(&self, x: &Reals) -> f64 { let n = x.len(); let mut cosines = 0.0; @@ -3081,6 +3157,9 @@ scalable!( /// Bounds [−5.12, 5.12]ⁿ; minimum 0 at the origin; 30 dimensions by default. `round` rounds /// halves away from 0, as MATLAB's does. /// + /// It supplies no gradient: away from the origin's (−½, ½), its derivative is 0 on the flat + /// steps and undefined at their jumps. + /// /// Liang, J. J., Qin, A. K., Suganthan, P. N. and Baskar, S. (2006). Comprehensive learning /// particle swarm optimizer for global optimization of multimodal functions. *IEEE /// Transactions on Evolutionary Computation* 10(3): 281-295, function f7 and table II @@ -3119,14 +3198,15 @@ scalable_problem!( url: "https://doi.org/10.1109/TEVC.2005.857610", ); -// Weierstrass's a, b and k_max -const WEIERSTRASS_A: f64 = 0.5; -const WEIERSTRASS_B: f64 = 3.0; -const WEIERSTRASS_TERMS: i32 = 21; +// Weierstrass's a, b and k_max + 1, the number of terms +pub(super) const WEIERSTRASS_A: f64 = 0.5; +pub(super) const WEIERSTRASS_B: f64 = 3.0; +pub(super) const WEIERSTRASS_TERMS: i32 = 21; -// Σₖ aᵏ cos(2π bᵏ (x + 0.5)), k from 0 to k_max -fn weierstrass_sum(x: f64) -> f64 { - (0..WEIERSTRASS_TERMS) +// Σₖ aᵏ cos(2π bᵏ (x + 0.5)), k from 0 to `terms` − 1: k_max + 1 terms for the function, fewer +// for the tests of its gradient +pub(super) fn weierstrass_sum(x: f64, terms: i32) -> f64 { + (0..terms) .map(|k| { let (ak, bk) = (math::powi(WEIERSTRASS_A, k), math::powi(WEIERSTRASS_B, k)); ak * math::cos(2.0 * PI * bk * (x + 0.5)) @@ -3144,6 +3224,10 @@ scalable!( /// 30 dimensions by default. Each gene's sum is at least −Σ aᵏ, reached where every cosine /// is −1, at the integers, and the second term is −n Σ aᵏ, since every bᵏ is odd. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with), + /// the exact derivative of the finite sum. Its highest terms vary on scales of 10⁻¹⁰, finer + /// than any finite difference of the computed function resolves. + /// /// Suganthan, P. N., Hansen, N., Liang, J. J., Deb, K., Chen, Y.-P., Auger, A. and Tiwari, S. /// (2005). *Problem Definitions and Evaluation Criteria for the CEC 2005 Special Session on /// Real-Parameter Optimization*, function F11 (read), shifted and rotated there: its @@ -3159,10 +3243,15 @@ scalable!( impl FitnessFunction for Weierstrass { type Output = f64; + gradient!(gradients::weierstrass); + fn evaluate(&self, x: &Reals) -> f64 { // the same sum at 0, so that the minimum is exactly 0 - let offset = x.len() as f64 * weierstrass_sum(0.0); - x.iter().map(|&xi| weierstrass_sum(xi)).sum::() - offset + let offset = x.len() as f64 * weierstrass_sum(0.0, WEIERSTRASS_TERMS); + x.iter() + .map(|&xi| weierstrass_sum(xi, WEIERSTRASS_TERMS)) + .sum::() + - offset } } @@ -3185,6 +3274,10 @@ scalable!( /// of 1/2 (where every term of the inner sums is 0): 21ⁿ global minima in the box. 30 /// dimensions by default. /// + /// It supplies no gradient: the function has kinks 2⁻³³ apart in each gene, where its inner + /// sums' slopes jump by up to 2, so its derivative changes on scales that no search step + /// resolves. + /// /// Hansen, N., Finck, S., Ros, R. and Auger, A. (2009). *Real-Parameter Black-Box /// Optimization Benchmarking 2009: Noiseless Functions Definitions.* INRIA research report /// RR-6829, function f23 (read), "based on the idea" of Katsuura, H. (1991). Continuous @@ -3239,6 +3332,10 @@ scalable!( /// Bounds [−5, 5]ⁿ; minimum 0 at (−1, …, −1), the only one: the second part is /// `Σ (xᵢ + 1)² / (2n)`, 0 only there, where the first is 0 too. 30 dimensions by default. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with). + /// The first term has a cusp on the sphere `Σ xᵢ² = n`, through the minimum, where its gradient + /// is taken as 0. + /// /// Beyer, H.-G. and Finck, S. (2012). HappyCat: a simple function class where well-known /// direct search algorithms do fail. *Parallel Problem Solving from Nature, PPSN XII*, LNCS /// 7491: 367-376, which couldn't be read. Its function has a parameter α that shapes the @@ -3264,6 +3361,8 @@ fn squares_and_sum(x: &Reals) -> (f64, f64) { impl FitnessFunction for HappyCat { type Output = f64; + gradient!(gradients::happy_cat); + fn evaluate(&self, x: &Reals) -> f64 { let n = x.len() as f64; let (squares, sum) = squares_and_sum(x); @@ -3289,6 +3388,10 @@ scalable!( /// Bounds [−5, 5]ⁿ; minimum 0 at (−1, …, −1), the only one: the second part is /// `Σ (xᵢ + 1)² / (2n)`, 0 only there, where the first is 0 too. 30 dimensions by default. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with). + /// The first term has a cusp where `(Σ xᵢ²)² = (Σ xᵢ)²`, through the minimum, where its + /// gradient is taken as 0. + /// /// Liang, J. J., Qu, B. Y. and Suganthan, P. N. (2013). *Problem Definitions and Evaluation /// Criteria for the CEC 2014 Special Session and Competition on Single Objective /// Real-Parameter Numerical Optimization.* Technical report 201311, Zhengzhou University and @@ -3304,6 +3407,8 @@ scalable!( impl FitnessFunction for HgBat { type Output = f64; + gradient!(gradients::hg_bat); + fn evaluate(&self, x: &Reals) -> f64 { let n = x.len() as f64; let (squares, sum) = squares_and_sum(x); @@ -3326,6 +3431,10 @@ scalable!( /// /// Bounds [−100, 100]ⁿ; minimum 0 at the origin; at least 2 dimensions, 30 by default. /// + /// Supplies its analytic gradient through [`evaluate_with`](FitnessFunction::evaluate_with). A + /// pair's `√sᵢ` has a cusp where both its genes are 0, at the minimum among others, where the + /// pair's term of the gradient is taken as 0. + /// /// Schaffer, J. D., Caruana, R. A., Eshelman, L. J. and Das, R. (1989). A study of control /// parameters affecting online performance of genetic algorithms for function optimization. /// *Proceedings of the Third International Conference on Genetic Algorithms*, Morgan @@ -3344,6 +3453,8 @@ scalable!( impl FitnessFunction for SchafferF7 { type Output = f64; + gradient!(gradients::schaffer_f7); + fn evaluate(&self, x: &Reals) -> f64 { let pairs = x.len().saturating_sub(1).max(1) as f64; let sum: f64 = x @@ -3378,6 +3489,9 @@ scalable!( /// /// Bounds [−65.536, 65.536]ⁿ; minimum 0 at the origin; 30 dimensions by default. /// + /// Supplies its analytic gradient, `2 (n − j + 1) xⱼ`, through + /// [`evaluate_with`](FitnessFunction::evaluate_with). + /// /// Its origin is unknown: definition and bounds as Molga and Smutnicki (2005, section 2.3, /// read) give them. Not yet checked against an original /// ([#168](https://github.com/tachsin/genoxide/issues/168)). @@ -3390,6 +3504,8 @@ scalable!( impl FitnessFunction for RotatedHyperEllipsoid { type Output = f64; + gradient!(gradients::rotated_hyper_ellipsoid); + fn evaluate(&self, x: &Reals) -> f64 { let mut prefix = 0.0; let mut sum = 0.0; diff --git a/src/problems/gradients.rs b/src/problems/gradients.rs index b4b2e94a..fa44cfc0 100644 --- a/src/problems/gradients.rs +++ b/src/problems/gradients.rs @@ -1,10 +1,13 @@ -//! The analytic gradients of the smooth classic functions: `∂f / ∂xᵢ` into `gradient[i]`, with -//! `gradient` as long as `x` and zeroed. Each is the derivative of the formula in its problem's -//! docs, written out by hand, and tested against central differences. +//! The analytic gradients of the classic functions that are differentiable almost everywhere: +//! `∂f / ∂xᵢ` into `gradient[i]`, with `gradient` as long as `x` and zeroed. Each is the +//! derivative of the formula in its problem's docs, written out by hand, and tested against +//! central differences. Where a function has no derivative, on a set of measure 0 (a cone's or a +//! cusp's apex), the gradient of the term with the kink is taken as 0, as Ackley's at its cone. use super::classic::{ FOXHOLES, HARTMANN_3_A, HARTMANN_3_P, HARTMANN_6_A, HARTMANN_6_P, HARTMANN_C, KOWALIK_A, KOWALIK_B_INVERSE, LANGERMANN_A, LANGERMANN_C, MICHALEWICZ_M, SHEKEL_A, SHEKEL_C, + WEIERSTRASS_A, WEIERSTRASS_B, WEIERSTRASS_TERMS, conditioning, oscillation, }; use crate::math; use std::f64::consts::PI; @@ -409,5 +412,247 @@ pub(super) fn kowalik(x: &[f64], gradient: &mut [f64]) { } } +// ---- the CEC- and BBOB-style functions ----------------------------------------------------------- + +// Σ |xᵢ|^(i+1) (i from 1): (i + 1) |xᵢ|^i sign(xᵢ), 0 at 0 +pub(super) fn sum_of_different_powers(x: &[f64], gradient: &mut [f64]) { + for (i, (g, &xi)) in gradient.iter_mut().zip(x).enumerate() { + let power = i as i32 + 2; + *g = f64::from(power) * math::powi(xi.abs(), power - 1) * xi.signum(); + } +} + +// Σ i xᵢ⁴: 4 i xᵢ³ +pub(super) fn quartic(x: &[f64], gradient: &mut [f64]) { + for (i, (g, &xi)) in gradient.iter_mut().zip(x).enumerate() { + *g = 4.0 * (i + 1) as f64 * math::powi(xi, 3); + } +} + +// the slope of Yao, Liu and Lin's penalty k (|x| − a)^m outside [−a, a]: k m (|x| − a)^(m − 1) +// sign(x), and 0 inside, where the penalty meets 0 with that slope +fn penalty_slope(x: f64, a: f64, k: f64, m: i32) -> f64 { + if x.abs() > a { + k * f64::from(m) * math::powi(x.abs() - a, m - 1) * x.signum() + } else { + 0.0 + } +} + +// (π / n) L + Σ u(xᵢ, 10, 100, 4), with yᵢ = 1 + (xᵢ + 1) / 4 and L = 10 sin²(πy₁) +// + Σ (yᵢ − 1)² (1 + 10 sin²(πyᵢ₊₁)) + (yₙ − 1)²: the derivatives of L in y, with +// d sin²(πy) / dy = π sin 2πy, times (π / n) dy/dx = π / (4n) +pub(super) fn penalized_1(x: &[f64], gradient: &mut [f64]) { + let n = x.len(); + if n == 0 { + return; + } + let y = |i: usize| 1.0 + (x[i] + 1.0) / 4.0; + gradient[0] += 10.0 * PI * math::sin(2.0 * PI * y(0)); + for i in 0..n - 1 { + let (yi, next) = (y(i), y(i + 1)); + gradient[i] += 2.0 * (yi - 1.0) * (1.0 + 10.0 * math::powi(math::sin(PI * next), 2)); + gradient[i + 1] += math::powi(yi - 1.0, 2) * 10.0 * PI * math::sin(2.0 * PI * next); + } + gradient[n - 1] += 2.0 * (y(n - 1) - 1.0); + let scale = PI / (4.0 * n as f64); + for (g, &xi) in gradient.iter_mut().zip(x) { + *g = scale * *g + penalty_slope(xi, 10.0, 100.0, 4); + } +} + +// 0.1 {sin²(3πx₁) + Σ (xᵢ − 1)² (1 + sin²(3πxᵢ₊₁)) + (xₙ − 1)² (1 + sin²(2πxₙ))} +// + Σ u(xᵢ, 5, 100, 4), with d sin²(cx) / dx = c sin 2cx +pub(super) fn penalized_2(x: &[f64], gradient: &mut [f64]) { + let n = x.len(); + if n == 0 { + return; + } + let three_pi = 3.0 * PI; + gradient[0] += three_pi * math::sin(2.0 * three_pi * x[0]); + for i in 0..n - 1 { + let (xi, next) = (x[i], x[i + 1]); + gradient[i] += 2.0 * (xi - 1.0) * (1.0 + math::powi(math::sin(three_pi * next), 2)); + gradient[i + 1] += math::powi(xi - 1.0, 2) * three_pi * math::sin(2.0 * three_pi * next); + } + let last = x[n - 1]; + let (sin, cos) = math::sin_cos(2.0 * PI * last); + gradient[n - 1] += + 2.0 * (last - 1.0) * (1.0 + sin * sin) + math::powi(last - 1.0, 2) * 4.0 * PI * sin * cos; + for (g, &xi) in gradient.iter_mut().zip(x) { + *g = 0.1 * *g + penalty_slope(xi, 5.0, 100.0, 4); + } +} + +// Σ cᵢ xᵢ², cᵢ = (10⁶)^((i−1)/(n−1)): 2 cᵢ xᵢ +pub(super) fn high_conditioned_elliptic(x: &[f64], gradient: &mut [f64]) { + let n = x.len(); + for (i, (g, xi)) in gradient.iter_mut().zip(x).enumerate() { + *g = 2.0 * conditioning(i, n) * xi; + } +} + +// x₁² + 10⁶ Σᵢ₌₂ⁿ xᵢ² +pub(super) fn bent_cigar(x: &[f64], gradient: &mut [f64]) { + for (i, (g, xi)) in gradient.iter_mut().zip(x).enumerate() { + *g = if i == 0 { 2.0 * xi } else { 2e6 * xi }; + } +} + +// 10⁶ x₁² + Σᵢ₌₂ⁿ xᵢ² +pub(super) fn discus(x: &[f64], gradient: &mut [f64]) { + for (i, (g, xi)) in gradient.iter_mut().zip(x).enumerate() { + *g = if i == 0 { 2e6 * xi } else { 2.0 * xi }; + } +} + +// √S, S = Σ |xᵢ|^pᵢ, pᵢ = 2 + 4 (i−1)/(n−1): pᵢ |xᵢ|^(pᵢ−1) sign(xᵢ) / (2√S). At the origin, the +// only point where S is 0, √S has a cone (along the first axis, it's |x₁|): its gradient is +// taken as 0 there +pub(super) fn different_powers(x: &[f64], gradient: &mut [f64]) { + let n = x.len(); + let exponent = |i: usize| 2.0 + 4.0 * i as f64 / (n.max(2) - 1) as f64; + let sum: f64 = x + .iter() + .enumerate() + .map(|(i, xi)| math::powf(xi.abs(), exponent(i))) + .sum(); + if sum <= 0.0 { + return; + } + let scale = 0.5 / sum.sqrt(); + for (i, (g, &xi)) in gradient.iter_mut().zip(x).enumerate() { + let p = exponent(i); + *g = scale * p * math::powf(xi.abs(), p - 1.0) * xi.signum(); + } +} + +// 10 (n − Σ cos 2πzᵢ) + Σ zᵢ² + 100 Σ max(0, |xᵢ| − 5)², zᵢ = sᵢ T_osz(xᵢ): +// (20π sin 2πzᵢ + 2zᵢ) sᵢ T'(xᵢ) + 200 max(0, |xᵢ| − 5) sign(xᵢ). With x̂ = ln |x| and +// T(x) = sign(x) exp(x̂ + w(x̂)), w(x̂) = 0.049 (sin c₁x̂ + sin c₂x̂), T'(x) = T(x) (1 + w'(x̂)) / x +// = exp(w(x̂)) (1 + w'(x̂)), between 0.11 and 2.1: T is increasing. At 0, where T has no +// derivative (T(x) / x oscillates as x goes to 0), the term 10 (1 − cos 2πz) + z² is O(z²) = +// O(x²) on either side, so its derivative is 0 +pub(super) fn buche_rastrigin(x: &[f64], gradient: &mut [f64]) { + let n = x.len(); + for (i, (g, &xi)) in gradient.iter_mut().zip(x).enumerate() { + let penalty = 200.0 * (xi.abs() - 5.0).max(0.0) * xi.signum(); + if xi == 0.0 { + *g = penalty; + continue; + } + let oscillated = oscillation(xi); + let mut scale = math::powf(10.0, 0.5 * i as f64 / (n.max(2) - 1) as f64); + if oscillated > 0.0 && i % 2 == 0 { + scale *= 10.0; + } + let z = scale * oscillated; + let logarithm = math::ln(xi.abs()); + let (c1, c2) = if xi > 0.0 { (10.0, 7.9) } else { (5.5, 3.1) }; + let (sin1, cos1) = math::sin_cos(c1 * logarithm); + let (sin2, cos2) = math::sin_cos(c2 * logarithm); + let wiggle = 0.049 * (sin1 + sin2); + let slope = 0.049 * (c1 * cos1 + c2 * cos2); + let derivative = math::exp(wiggle) * (1.0 + slope); + let outer = 20.0 * PI * math::sin(2.0 * PI * z) + 2.0 * z; + *g = outer * scale * derivative + penalty; + } +} + +// the derivative of Weierstrass's sum of its first `terms` terms, Σₖ aᵏ cos(2π bᵏ (x + 0.5)): +// −Σₖ aᵏ 2π bᵏ sin(2π bᵏ (x + 0.5)) = Σₖ aᵏ 2π bᵏ sin(2π bᵏ x), since every bᵏ is odd and +// sin(θ + π bᵏ) = −sin θ. The phase from x, not x + 0.5, is rounded the less the nearer x is to +// the minimum at 0, where every term is 0 +pub(super) fn weierstrass_slope(x: f64, terms: i32) -> f64 { + (0..terms) + .map(|k| { + let (ak, bk) = (math::powi(WEIERSTRASS_A, k), math::powi(WEIERSTRASS_B, k)); + ak * 2.0 * PI * bk * math::sin(2.0 * PI * bk * x) + }) + .sum() +} + +// Σᵢ Σₖ aᵏ cos(2π bᵏ (xᵢ + 0.5)) − n Σₖ aᵏ cos(π bᵏ): each gene's sum's derivative +pub(super) fn weierstrass(x: &[f64], gradient: &mut [f64]) { + for (g, &xi) in gradient.iter_mut().zip(x) { + *g = weierstrass_slope(xi, WEIERSTRASS_TERMS); + } +} + +// Σ xᵢ² and Σ xᵢ, in the order of the functions' own sums +fn squares_and_sum(x: &[f64]) -> (f64, f64) { + x.iter().fold((0.0, 0.0), |(squares, sum), xi| { + (squares + xi * xi, sum + xi) + }) +} + +// |S − n|^(1/4) + (S / 2 + T) / n + 1/2, S = Σ xᵢ², T = Σ xᵢ: (1/4) |S − n|^(−3/4) sign(S − n) 2xᵢ +// = xᵢ |S − n|^(1/4) / (2 (S − n)), plus (xᵢ + 1) / n. On the sphere S = n, where the first term +// has a cusp (one-sided slopes of −∞ and +∞ across it), its gradient is taken as 0 +pub(super) fn happy_cat(x: &[f64], gradient: &mut [f64]) { + let n = x.len() as f64; + let (squares, _) = squares_and_sum(x); + let difference = squares - n; + let groove = if difference != 0.0 { + 0.5 * math::powf(difference.abs(), 0.25) / difference + } else { + 0.0 + }; + for (g, &xi) in gradient.iter_mut().zip(x) { + *g = groove * xi + (xi + 1.0) / n; + } +} + +// |S² − T²|^(1/2) + (S / 2 + T) / n + 1/2, S = Σ xᵢ², T = Σ xᵢ: with U = S² − T², +// sign(U) (4S xᵢ − 2T) / (2 √|U|) = (2S xᵢ − T) √|U| / U, plus (xᵢ + 1) / n. Where U is 0, where +// the first term has a cusp, its gradient is taken as 0 +pub(super) fn hg_bat(x: &[f64], gradient: &mut [f64]) { + let n = x.len() as f64; + let (squares, sum) = squares_and_sum(x); + let difference = squares * squares - sum * sum; + let groove = if difference != 0.0 { + difference.abs().sqrt() / difference + } else { + 0.0 + }; + for (g, &xi) in gradient.iter_mut().zip(x) { + *g = groove * (2.0 * squares * xi - sum) + (xi + 1.0) / n; + } +} + +// (S / P)², S = Σᵢ h(sᵢ), P = n − 1 pairs, sᵢ = √(xᵢ² + xᵢ₊₁²), h(s) = √s (1 + sin²(50 s^(1/5))): +// 2 S / P² Σᵢ h'(sᵢ) ∂sᵢ/∂x, with h'(s) = (1 + sin²θ) / (2√s) + 10 sin(2θ) s^(−3/10), +// θ = 50 s^(1/5), and ∂sᵢ/∂xⱼ = xⱼ / sᵢ for the pair's two genes. Where a pair's genes are both +// 0, its √s has a cusp, and its term of the gradient is taken as 0 +pub(super) fn schaffer_f7(x: &[f64], gradient: &mut [f64]) { + let pairs = x.len().saturating_sub(1).max(1) as f64; + let mut sum = 0.0; + for (i, pair) in x.windows(2).enumerate() { + let s = (pair[0] * pair[0] + pair[1] * pair[1]).sqrt(); + let root = s.sqrt(); + let fifth = math::powf(s, 0.2); + let (sin, cos) = math::sin_cos(50.0 * fifth); + sum += root * (1.0 + sin * sin); + if s > 0.0 { + // 10 sin(2θ) s^(−3/10) = 20 sin θ cos θ s^(1/5) / √s + let slope = (1.0 + sin * sin) / (2.0 * root) + 20.0 * sin * cos * fifth / root; + gradient[i] += slope * pair[0] / s; + gradient[i + 1] += slope * pair[1] / s; + } + } + let outer = 2.0 * sum / (pairs * pairs); + for g in gradient.iter_mut() { + *g *= outer; + } +} + +// Σⱼ (n − j + 1) xⱼ² (j from 1): 2 (n − j + 1) xⱼ +pub(super) fn rotated_hyper_ellipsoid(x: &[f64], gradient: &mut [f64]) { + let n = x.len(); + for (j, (g, xi)) in gradient.iter_mut().zip(x).enumerate() { + *g = 2.0 * (n - j) as f64 * xi; + } +} + #[cfg(test)] mod tests; diff --git a/src/problems/gradients/tests.rs b/src/problems/gradients/tests.rs index 6cff5445..76c539bc 100644 --- a/src/problems/gradients/tests.rs +++ b/src/problems/gradients/tests.rs @@ -1,12 +1,20 @@ -use crate::StreamRng; -use crate::engine::Extras; +use crate::engine::{Extras, FitnessFunction, Provided}; use crate::genome::{Reals, Representation}; use crate::gradient; -use crate::problems::{self, DynProblem}; +use crate::problems::classic::{WEIERSTRASS_B, WEIERSTRASS_TERMS, weierstrass_sum}; +use crate::problems::{ + self, BentCigar, BucheRastrigin, DifferentPowers, Discus, DynProblem, HappyCat, HgBat, + HighConditionedElliptic, Katsuura, NonContinuousRastrigin, Penalized1, Penalized2, Problem, + Quartic, RotatedHyperEllipsoid, SchafferF7, Step, SumOfDifferentPowers, Weierstrass, +}; +use crate::{Fitness, StreamRng}; +use std::f64::consts::PI; -// the problems with an analytic gradient: every classic function but the three that aren't -// differentiable (a maximum, absolute values) -const DIFFERENTIABLE: [&str; 36] = [ +// the problems with an analytic gradient: every classic function but those whose derivative is +// undefined or 0 on sets of positive measure, or meaningless: a maximum and absolute values (the +// eggholder, Schwefel 2.21 and 2.22), flat steps (the step function, the non-continuous +// Rastrigin) and kinks every 2⁻³³ (Katsuura) +const DIFFERENTIABLE: [&str; 50] = [ "Sphere", "AxisParallelEllipsoid", "Schwefel1_2", @@ -43,8 +51,27 @@ const DIFFERENTIABLE: [&str; 36] = [ "Langermann", "ShekelFoxholes", "Kowalik", + "SumOfDifferentPowers", + "Quartic", + "Penalized1", + "Penalized2", + "HighConditionedElliptic", + "BentCigar", + "Discus", + "DifferentPowers", + "BucheRastrigin", + "Weierstrass", + "HappyCat", + "HgBat", + "SchafferF7", + "RotatedHyperEllipsoid", ]; +// the problems whose minimum is on a kink, where central differences straddle it and the +// gradient is 0 by convention: Büche-Rastrigin's (whose T_osz has no derivative at 0), HappyCat's +// and HGBat's (cusps, on a sphere and two) and Schaffer F7's (a cusp of each pair) +const KINKED_AT_THE_OPTIMUM: [&str; 4] = ["BucheRastrigin", "HappyCat", "HgBat", "SchafferF7"]; + // the analytic gradient of `problem` at `x` against central differences: the value is the plain // evaluation's, to the bit, and the gradient's error is within 1e-5 relative (absolute below 1), // plus the error of the differences themselves: their rounding error, about ε |f| / h, which @@ -101,7 +128,8 @@ fn the_smooth_problems_have_gradients() { fn gradients_match_central_differences_at_random_points() { let mut rng = StreamRng::seed_from_u64(3); for problem in problems::all() { - if !problem.provides().gradient { + // Weierstrass's highest terms are beyond differences (see its own test below) + if !problem.provides().gradient || problem.name() == "Weierstrass" { continue; } let real = problem.real(); @@ -124,9 +152,13 @@ fn gradients_match_central_differences_near_the_optimum() { } let real = problem.real(); let optimum = problem.optimum().expect("known"); + // Weierstrass's highest terms are beyond differences (see its own test below) + let differences = problem.name() != "Weierstrass"; let mut largest: f64 = 0.0; for solution in optimum.solutions() { - largest = largest.max(error(problem.as_ref(), solution)); + if differences && !KINKED_AT_THE_OPTIMUM.contains(&problem.name()) { + largest = largest.max(error(problem.as_ref(), solution)); + } // the optima are interior minima, rounded to the precision of f64: a gradient of // about 1e-14, but Michalewicz's 2e-6, whose minimizers are found by golden-section // search on terms as steep as sin²⁰ @@ -134,14 +166,17 @@ fn gradients_match_central_differences_near_the_optimum() { problem.evaluate_with(solution, &mut Extras::with_gradient(&mut gradient)); let norm = gradient.iter().fold(0.0, |norm: f64, g| norm.max(g.abs())); assert!(norm <= 1e-5, "{}: {norm}", problem.name()); - for _ in 0..20 { - // a thousandth of each range away, inside the bounds + for _ in 0..if differences { 20 } else { 0 } { + // about a thousandth of each range away, inside the bounds let x: Reals = solution .iter() .zip(real.bounds()) .map(|(&xi, range)| { let width = range.end() - range.start(); - let offset = (rng.unit_f64() * 2.0 - 1.0) * 1e-3 * width; + // from half to all of a thousandth: away from a kink at the optimum, where + // central differences straddle it + let u = rng.unit_f64() * 2.0 - 1.0; + let offset = (0.5 + 0.5 * u.abs()) * u.signum() * 1e-3 * width; (xi + offset).clamp(*range.start(), *range.end()) }) .collect(); @@ -151,3 +186,217 @@ fn gradients_match_central_differences_near_the_optimum() { assert!(largest <= 1.0, "{}: {largest}", problem.name()); } } + +// a problem behind `DynProblem` as a fitness function, for `gradient::check` +struct Function<'a>(&'a dyn DynProblem); + +impl FitnessFunction for Function<'_> { + type Output = Fitness; + + fn evaluate(&self, x: &Reals) -> Fitness { + self.0.evaluate(x) + } + + fn provides(&self) -> Provided { + self.0.provides() + } + + fn evaluate_with(&self, x: &Reals, extras: &mut Extras<'_>) -> Fitness { + self.0.evaluate_with(x, extras) + } +} + +// the CEC- and BBOB-style functions with a gradient that `gradient::check`'s differences resolve: +// all but Büche-Rastrigin and Weierstrass, whose ripples need finer steps (below) +fn checked() -> Vec> { + vec![ + problems::boxed(SumOfDifferentPowers::default()), + problems::boxed(Quartic::default()), + problems::boxed(Penalized1::default()), + problems::boxed(Penalized2::default()), + problems::boxed(HighConditionedElliptic::default()), + problems::boxed(BentCigar::default()), + problems::boxed(Discus::default()), + problems::boxed(DifferentPowers::default()), + problems::boxed(HappyCat::default()), + problems::boxed(HgBat::default()), + problems::boxed(SchafferF7::default()), + problems::boxed(RotatedHyperEllipsoid::default()), + ] +} + +// `gradient::check` at 200 random points of the box, which are away from the kinks (a set of +// measure 0) almost surely: every gene within 5e-8, relative (absolute below 1), of central +// differences, beyond the differences' own rounding error, about ε |f| / h, which is the larger +// where the value dwarfs a gene's slope (the bent cigar's first gene, the penalties' far genes) +#[test] +fn the_cec_and_bbob_style_gradients_pass_gradient_check() { + let mut rng = StreamRng::seed_from_u64(5); + for problem in checked() { + let real = problem.real(); + let mut largest: f64 = 0.0; + for _ in 0..200 { + let x = real.random_genome(&mut rng); + let check = gradient::check(&Function(problem.as_ref()), &x).expect("a gradient"); + let value = problem.evaluate(&x).score().expect("valid").abs().max(1.0); + for (gene, &error) in check.errors().iter().enumerate() { + let h = gradient::CENTRAL_STEP * x[gene].abs().max(1.0); + let scale = check.supplied()[gene] + .abs() + .max(check.estimated()[gene].abs()) + .max(1.0); + let rounding = 1e2 * f64::EPSILON * value / (h * scale); + largest = largest.max(error - rounding); + } + } + assert!(largest <= 5e-8, "{}: {largest}", problem.name()); + } +} + +// the slope of `f` at `x` by Richardson's extrapolation of central differences with the steps h +// and h / 2, whose error is O(h⁴) rather than O(h²) +fn richardson(f: impl Fn(f64) -> f64, x: f64, h: f64) -> f64 { + let central = |h: f64| { + let (above, below) = (x + h, x - h); + (f(above) - f(below)) / (above - below) + }; + (4.0 * central(h / 2.0) - central(h)) / 3.0 +} + +// Büche-Rastrigin's rings are up to 2000 radians per unit of a gene, and its oscillation T_osz +// wiggles at 10 / |x| near 0: central differences at gradient::check's steps are off by 1e-6 to +// 1e-4 there. A gene's slope, with the others fixed, by Richardson's extrapolation at steps of +// 2⁻²⁰ is within 1e-8 (relative, absolute below 1), beyond the rounding error +#[test] +fn the_buche_rastrigin_gradient_matches_finer_differences() { + let problem = BucheRastrigin::new(10); + let real = problem.representation(); + let mut rng = StreamRng::seed_from_u64(6); + let mut largest: f64 = 0.0; + for _ in 0..200 { + let x = real.random_genome(&mut rng); + let mut gradient = vec![0.0; x.len()]; + let value = problem.evaluate_with(&x, &mut Extras::with_gradient(&mut gradient)); + assert_eq!(value.to_bits(), problem.evaluate(&x).to_bits()); + for (gene, &g) in gradient.iter().enumerate() { + let h = exp2_floor(x[gene].abs().max(1.0)) * 2f64.powi(-20); + let along = |xi: f64| { + let mut point = x.clone(); + point[gene] = xi; + problem.evaluate(&point) + }; + let estimate = richardson(along, x[gene], h); + let scale = g.abs().max(estimate.abs()).max(1.0); + let rounding = 1e2 * f64::EPSILON * value.abs().max(1.0) / h; + largest = largest.max(((g - estimate).abs() - rounding) / scale); + } + } + assert!(largest <= 1e-8, "{largest}"); +} + +// Weierstrass's sum has terms up to 3²⁰ times the first's frequency, and the computed phase +// 2π 3ᵏ (x + 0.5) rounds by up to about ε 2π 3ᵏ, 2·10⁻⁶ rad for the last term: no differences of +// the computed function resolve the derivative of its highest terms. Its slope is checked on the +// sums of its first 1 to 7 terms (to 3⁶ times the first frequency; the phases' rounding shows +// from 3⁷ on, 3e-8 at 3⁷, 6e-8 at 3⁹, 4e-6 at 3¹⁰), with Richardson's extrapolation at steps of +// a hundredth of a radian of the highest term: within 1e-8 (relative, absolute below 1). Every +// term has the same formula, and the gradient is the slope of the whole sum at each gene +#[test] +fn the_weierstrass_gradient_matches_differences_of_its_partial_sums() { + let mut rng = StreamRng::seed_from_u64(7); + let mut largest: f64 = 0.0; + for terms in 1..=7 { + let frequency = 2.0 * PI * WEIERSTRASS_B.powi(terms - 1); + let h = exp2_floor(0.01 / frequency); + for _ in 0..200 { + let x = rng.unit_f64() - 0.5; + let slope = super::weierstrass_slope(x, terms); + let estimate = richardson(|x| weierstrass_sum(x, terms), x, h); + let error = (slope - estimate).abs() / slope.abs().max(estimate.abs()).max(1.0); + largest = largest.max(error); + } + } + assert!(largest <= 1e-8, "{largest}"); + // the whole function's gradient: each gene's slope of all 21 terms + let problem = Weierstrass::new(5); + let x = problem.representation().random_genome(&mut rng); + let mut gradient = vec![0.0; 5]; + let value = problem.evaluate_with(&x, &mut Extras::with_gradient(&mut gradient)); + assert_eq!(value.to_bits(), problem.evaluate(&x).to_bits()); + for (g, &xi) in gradient.iter().zip(x.iter()) { + assert_eq!(*g, super::weierstrass_slope(xi, WEIERSTRASS_TERMS)); + } + // 0 at the minimum, where every sin(2π 3ᵏ x) is 0 + problem.evaluate_with( + &Reals::from(vec![0.0; 5]), + &mut Extras::with_gradient(&mut gradient), + ); + assert_eq!(gradient, [0.0; 5]); +} + +// the noisy quartic's value jumps between any two genomes: it has no gradient +#[test] +fn the_noisy_quartic_has_no_gradient() { + assert!(Quartic::new(3).provides().gradient); + assert!(!Quartic::noisy(3).provides().gradient); + let x = Reals::from(vec![0.5, -0.25, 1.0]); + let mut gradient = [7.0; 3]; + let noisy = Quartic::noisy(3); + let value = noisy.evaluate_with(&x, &mut Extras::with_gradient(&mut gradient)); + assert_eq!(value.to_bits(), noisy.evaluate(&x).to_bits()); + assert_eq!(gradient, [7.0; 3]); + // without noise, 4 i xᵢ³ + Quartic::new(3).evaluate_with(&x, &mut Extras::with_gradient(&mut gradient)); + assert_eq!(gradient, [0.5, -0.125, 12.0]); +} + +// the functions whose derivative is undefined or 0 on sets of positive measure provide none +#[test] +fn the_flat_and_rugged_functions_have_no_gradient() { + assert!(!Step::default().provides().gradient); + assert!(!NonContinuousRastrigin::default().provides().gradient); + assert!(!Katsuura::default().provides().gradient); +} + +// at the kinks, each gradient's term with the kink is 0: the cones' and cusps' apexes +#[test] +fn gradients_at_the_kinks() { + let gradient_at = |problem: Box, x: &[f64]| { + let mut gradient = vec![f64::NAN; x.len()]; + problem.evaluate_with( + &Reals::from(x.to_vec()), + &mut Extras::with_gradient(&mut gradient), + ); + gradient + }; + // the apexes are the minima + assert_eq!( + gradient_at(problems::boxed(DifferentPowers::new(3)), &[0.0; 3]), + [0.0; 3] + ); + assert_eq!( + gradient_at(problems::boxed(BucheRastrigin::new(3)), &[0.0; 3]), + [0.0; 3] + ); + assert_eq!( + gradient_at(problems::boxed(HappyCat::new(3)), &[-1.0; 3]), + [0.0; 3] + ); + assert_eq!( + gradient_at(problems::boxed(HgBat::new(3)), &[-1.0; 3]), + [0.0; 3] + ); + assert_eq!( + gradient_at(problems::boxed(SchafferF7::new(3)), &[0.0; 3]), + [0.0; 3] + ); + // on HappyCat's sphere ‖x‖² = n elsewhere, only the slope (xᵢ + 1) / n is left + assert_eq!( + gradient_at(problems::boxed(HappyCat::new(2)), &[1.0, -1.0]), + [1.0, 0.0] + ); + // Schaffer F7's pair (0, 0) adds nothing, the pair (0, 3) its slope along the third gene + let schaffer = gradient_at(problems::boxed(SchafferF7::new(3)), &[0.0, 0.0, 3.0]); + assert_eq!(schaffer[0], 0.0); + assert!(schaffer[1] == 0.0 && schaffer[2] != 0.0); +} diff --git a/src/problems/transform.rs b/src/problems/transform.rs index 9c91a3af..3cff98a3 100644 --- a/src/problems/transform.rs +++ b/src/problems/transform.rs @@ -4,7 +4,7 @@ use super::{Constraints, Optimum, Problem}; use crate::Objective; use crate::StreamRng; -use crate::engine::FitnessFunction; +use crate::engine::{Extras, FitnessFunction, Provided}; use crate::genome::{Real, Reals, Representation}; // the streams of the seed's generator from which the shift and the rotation are drawn, so that a @@ -66,6 +66,13 @@ fn moved_optimum( /// in `x − o` can put a shifted solution a few ulps outside a constraint that's active there /// (6·10⁻¹⁴ at [`G06`](super::cec2006::G06)'s minimum with seed 3). /// +/// It [provides](FitnessFunction::provides) what the wrapped problem provides, and +/// [`evaluate_with`](FitnessFunction::evaluate_with) gives it at `x − o`: the gradient, the +/// constraints' values and their Jacobian, unchanged, since a shift moves the function without +/// turning or stretching it. A shifted smooth function keeps its analytic gradient, and a shifted +/// constrained problem gives its constraints' values to [`Bo`](crate::algorithm::Bo) (and, with +/// their Jacobian, to [`Mma`](crate::algorithm::Mma)). +/// /// ``` /// use genoxide::genome::Representation; /// use genoxide::problems::{Problem, Rastrigin, Shifted}; @@ -135,6 +142,21 @@ impl> FitnessFunction for Shifted

{ fn evaluate(&self, x: &Reals) -> P::Output { self.problem.evaluate(&self.unshifted(x)) } + + /// What the wrapped problem provides. + fn provides(&self) -> Provided { + self.problem.provides() + } + + /// The wrapped problem's [`evaluate_with`](FitnessFunction::evaluate_with) at `x − o`: its + /// fitness and extras, the same in `x` as in `x − o`. + /// + /// # Panics + /// + /// As the wrapped problem's. + fn evaluate_with(&self, x: &Reals, extras: &mut Extras<'_>) -> P::Output { + self.problem.evaluate_with(&self.unshifted(x), extras) + } } impl> Problem for Shifted

{ @@ -196,6 +218,13 @@ impl> Problem for Shifted

{ /// box can have lower values there, as for [`Shifted`]. The name and the reference are the wrapped /// problem's, and so is the constraint violation of a constrained problem, at `c + M (x − c)`. /// +/// It [provides](FitnessFunction::provides) what the wrapped problem provides, and +/// [`evaluate_with`](FitnessFunction::evaluate_with) gives it by the chain rule: the gradient +/// `Mᵀ ∇f(c + M (x − c))`, the constraints' values at `c + M (x − c)`, and their Jacobian `J M`, +/// each row (a constraint's gradient) turned as the gradient is. A rotated smooth function keeps +/// its analytic gradient, and a rotated constrained problem gives its constraints' values to +/// [`Bo`](crate::algorithm::Bo) (and, with their Jacobian, to [`Mma`](crate::algorithm::Mma)). +/// /// Rotating a shifted problem gives CEC 2005's shifted rotated functions, `f((x − o) M)`, e.g. its /// F10, the shifted rotated Rastrigin, whose rotation turns about the shifted optimum: /// @@ -261,6 +290,18 @@ impl> Rotated

{ &self.center } + // `Mᵀ v` into `out`: a gradient at c + M (x − c), or a row of the Jacobian there, as a + // gradient in x (∂/∂xⱼ = Σᵢ Mᵢⱼ ∂/∂yᵢ, y = c + M (x − c)) + fn turn_back(&self, v: &[f64], out: &mut [f64]) { + let n = self.center.len(); + out.fill(0.0); + for (row, vi) in self.matrix.chunks_exact(n.max(1)).zip(v) { + for (o, m) in out.iter_mut().zip(row) { + *o += m * vi; + } + } + } + // the point at which the wrapped problem is evaluated: c + M (x − c) fn rotated(&self, x: &Reals) -> Reals { let n = self.center.len(); @@ -311,6 +352,55 @@ impl> FitnessFunction for Rotated

{ fn evaluate(&self, x: &Reals) -> P::Output { self.problem.evaluate(&self.rotated(x)) } + + /// What the wrapped problem provides. + fn provides(&self) -> Provided { + self.problem.provides() + } + + /// The wrapped problem's [`evaluate_with`](FitnessFunction::evaluate_with) at + /// `c + M (x − c)`: its fitness and the constraints' values there, its gradient turned back + /// by `Mᵀ`, and the Jacobian times `M`. + /// + /// # Panics + /// + /// As the wrapped problem's, and if the gradient doesn't have a value per gene or the + /// Jacobian a whole number of rows of a value per gene. + fn evaluate_with(&self, x: &Reals, extras: &mut Extras<'_>) -> P::Output { + let n = self.center.len(); + let (gradient, inequalities, jacobian) = extras.buffers(); + // the wrapped problem's derivatives at c + M (x − c), in its own coordinates + let mut wrapped_gradient = gradient.as_ref().map(|gradient| { + assert_eq!(gradient.len(), n, "a gradient has a value per gene"); + vec![0.0; n] + }); + let mut wrapped_jacobian = jacobian.as_ref().map(|jacobian| { + assert_eq!( + jacobian.len() % n.max(1), + 0, + "a Jacobian has a row of a value per gene for each constraint" + ); + vec![0.0; jacobian.len()] + }); + let output = self.problem.evaluate_with( + &self.rotated(x), + &mut Extras::new( + wrapped_gradient.as_deref_mut(), + inequalities, + wrapped_jacobian.as_deref_mut(), + ), + ); + if let (Some(gradient), Some(wrapped)) = (gradient, &wrapped_gradient) { + self.turn_back(wrapped, gradient); + } + if let (Some(jacobian), Some(wrapped)) = (jacobian, &wrapped_jacobian) { + let rows = jacobian.chunks_exact_mut(n.max(1)); + for (row, wrapped) in rows.zip(wrapped.chunks_exact(n.max(1))) { + self.turn_back(wrapped, row); + } + } + output + } } impl> Problem for Rotated

{ @@ -363,8 +453,13 @@ impl> Problem for Rotated

{ #[cfg(test)] mod tests { use super::*; + use crate::gradient::Gradients; + use crate::prelude::{Engine, Lbfgsb, Stop, StopReason}; use crate::problems::cec2006::G06; - use crate::problems::{Branin, Rastrigin, Rosenbrock, Sphere}; + use crate::problems::{ + BentCigar, Branin, HighConditionedElliptic, Katsuura, Quartic, Rastrigin, Rosenbrock, + Sphere, Step, boxed, + }; fn assert_close(actual: f64, expected: f64, tolerance: f64) { assert!( @@ -507,4 +602,277 @@ mod tests { let center = Reals::from(rotated.center().to_vec()); assert_eq!(rotated.constraints(¢er), G06.constraints(¢er)); } + + // a smooth constrained problem with every extra: Σ (xᵢ − 1)² + x₀ x₁, subject to + // ‖x‖² − 4 <= 0 and x₀ + 2x₁ − x₂³ <= 0, with its gradient and its constraints' Jacobian + #[derive(Clone, Debug, PartialEq)] + struct Smooth; + + impl Smooth { + fn values(x: &[f64]) -> [f64; 2] { + let squares: f64 = x.iter().map(|xi| xi * xi).sum(); + [squares - 4.0, x[0] + 2.0 * x[1] - x[2] * x[2] * x[2]] + } + } + + impl FitnessFunction for Smooth { + type Output = (f64, f64); + + fn evaluate(&self, x: &Reals) -> (f64, f64) { + let value = x.iter().map(|xi| (xi - 1.0) * (xi - 1.0)).sum::() + x[0] * x[1]; + (value, self.constraints(x).violation(0.0)) + } + + fn provides(&self) -> Provided { + Provided::GRADIENT + .with_inequalities(2) + .with_constraint_jacobian() + } + + fn evaluate_with(&self, x: &Reals, extras: &mut Extras<'_>) -> (f64, f64) { + let (gradient, inequalities, jacobian) = extras.buffers(); + if let Some(gradient) = gradient { + for (g, xi) in gradient.iter_mut().zip(x.iter()) { + *g = 2.0 * (xi - 1.0); + } + gradient[0] += x[1]; + gradient[1] += x[0]; + } + if let Some(g) = inequalities { + g.copy_from_slice(&Self::values(x)); + } + if let Some(jacobian) = jacobian { + let (first, second) = jacobian.split_at_mut(3); + for (j, xi) in first.iter_mut().zip(x.iter()) { + *j = 2.0 * xi; + } + second.copy_from_slice(&[1.0, 2.0, -3.0 * x[2] * x[2]]); + } + self.evaluate(x) + } + } + + impl Problem for Smooth { + type Representation = Real; + + fn name(&self) -> &'static str { + "Smooth" + } + + fn representation(&self) -> Real { + Real::uniform(3, -2.0..=2.0).expect("valid bounds") + } + + fn optimum(&self) -> Option> { + None + } + + fn reference(&self) -> &'static str { + "a test problem" + } + + fn constraints(&self, x: &Reals) -> Constraints { + Constraints::new(Self::values(x).to_vec(), Vec::new()) + } + } + + // the gradient, the constraints' values and their Jacobian of `problem` at `x`, all wanted + fn extras_of

(problem: &P, x: &Reals) -> (P::Output, Vec, Vec, Vec) + where + P: Problem, + { + let (n, m) = (x.len(), problem.provides().inequalities); + let (mut gradient, mut g, mut jacobian) = (vec![0.0; n], vec![0.0; m], vec![0.0; m * n]); + let output = problem.evaluate_with( + x, + &mut Extras::new(Some(&mut gradient), Some(&mut g), Some(&mut jacobian)), + ); + (output, gradient, g, jacobian) + } + + // the slope of `f` along gene `gene` at `x`, by central differences with the step 2⁻¹⁸ + fn slope(f: impl Fn(&Reals) -> f64, x: &Reals, gene: usize) -> f64 { + let h = 2f64.powi(-18); + let (mut above, mut below) = (x.clone(), x.clone()); + above[gene] += h; + below[gene] -= h; + (f(&above) - f(&below)) / (above[gene] - below[gene]) + } + + // the wrappers' gradients and Jacobians against central differences of their values, and + // their values the same to the bit with and without extras + #[test] + fn the_wrappers_derivatives_match_central_differences() { + let mut rng = StreamRng::seed_from_u64(1); + for seed in 0..10 { + let shifted = Shifted::new(Smooth, seed); + let rotated = Rotated::new(Smooth, seed); + let both = Rotated::new(Shifted::new(Smooth, seed), seed); + for _ in 0..20 { + let x = Smooth.representation().random_genome(&mut rng); + check_derivatives(&shifted, &x); + check_derivatives(&rotated, &x); + check_derivatives(&both, &x); + } + } + } + + fn check_derivatives

(problem: &P, x: &Reals) + where + P: Problem, + { + let (output, gradient, g, jacobian) = extras_of(problem, x); + let plain = problem.evaluate(x); + assert_eq!(output.0.to_bits(), plain.0.to_bits()); + assert_eq!(output.1.to_bits(), plain.1.to_bits()); + assert_eq!(g, problem.constraints(x).inequalities()); + let close = |analytic: f64, estimate: f64| { + let error = (analytic - estimate).abs() / analytic.abs().max(1.0); + assert!(error < 1e-8, "{analytic} against {estimate}"); + }; + for gene in 0..x.len() { + close(gradient[gene], slope(|x| problem.evaluate(x).0, x, gene)); + for constraint in 0..2 { + let estimate = slope( + |x| problem.constraints(x).inequalities()[constraint], + x, + gene, + ); + close(jacobian[constraint * x.len() + gene], estimate); + } + } + } + + // a shift doesn't change the extras: the wrapped problem's at x − o, to the bit + #[test] + fn a_shift_passes_the_extras_through() { + let problem = Shifted::new(Smooth, 4); + let x = Reals::from(vec![0.5, -1.0, 1.5]); + let unshifted = problem.unshifted(&x); + assert_eq!(extras_of(&problem, &x), extras_of(&Smooth, &unshifted)); + // a test problem's analytic gradient, shifted + let rosenbrock = Shifted::new(Rosenbrock::new(4), 2); + let x = Reals::from(vec![0.3, -0.2, 1.1, 2.0]); + let mut gradient = [0.0; 4]; + let mut expected = [0.0; 4]; + let value = rosenbrock.evaluate_with(&x, &mut Extras::with_gradient(&mut gradient)); + let wrapped = Rosenbrock::new(4).evaluate_with( + &rosenbrock.unshifted(&x), + &mut Extras::with_gradient(&mut expected), + ); + assert_eq!((value, gradient), (wrapped, expected)); + } + + // a rotation turns the gradient back by Mᵀ and the Jacobian's rows likewise + #[test] + fn a_rotation_turns_the_derivatives_back() { + let problem = Rotated::new(Smooth, 6); + let x = Reals::from(vec![0.5, -1.0, 1.5]); + let y = problem.rotated(&x); + let (output, gradient, g, jacobian) = extras_of(&problem, &x); + let (wrapped, inner_gradient, inner_g, inner_jacobian) = extras_of(&Smooth, &y); + assert_eq!((output, g), (wrapped, inner_g)); + let m = problem.matrix(); + let turned = |v: &[f64], j: usize| (0..3).map(|i| m[i * 3 + j] * v[i]).sum::(); + for j in 0..3 { + assert!((gradient[j] - turned(&inner_gradient, j)).abs() < 1e-14); + for row in 0..2 { + let expected = turned(&inner_jacobian[row * 3..row * 3 + 3], j); + assert!((jacobian[row * 3 + j] - expected).abs() < 1e-14); + } + } + // the gradient of a sphere about its center is turned to the same one: M orthogonal + let sphere = Rotated::new(Sphere::new(5), 3); + let x = Reals::from(vec![1.0, -2.0, 0.5, 3.0, -1.5]); + let mut gradient = [0.0; 5]; + sphere.evaluate_with(&x, &mut Extras::with_gradient(&mut gradient)); + for (g, xi) in gradient.iter().zip(x.iter()) { + assert!((g - 2.0 * xi).abs() < 1e-13, "{g} {xi}"); + } + } + + // the wrappers provide what they wrap + #[test] + fn the_wrappers_provide_what_they_wrap() { + let gradient = Provided::GRADIENT; + assert_eq!(Shifted::new(Sphere::new(3), 1).provides(), gradient); + assert_eq!(Rotated::new(Rastrigin::new(3), 1).provides(), gradient); + let both = Rotated::new(Shifted::new(BentCigar::new(3), 1), 1); + assert_eq!(both.provides(), gradient); + assert_eq!( + Rotated::new(Shifted::new(Smooth, 1), 1).provides(), + Smooth.provides() + ); + let values = Provided::NOTHING.with_inequalities(2); + assert_eq!(Shifted::new(G06, 1).provides(), values); + assert_eq!(Rotated::new(G06, 1).provides(), values); + assert_eq!(Rotated::new(Shifted::new(G06, 1), 1).provides(), values); + let nothing = Provided::NOTHING; + assert_eq!(Shifted::new(Step::new(3), 1).provides(), nothing); + assert_eq!(Rotated::new(Katsuura::new(3), 1).provides(), nothing); + assert_eq!( + Rotated::new(Shifted::new(Quartic::noisy(3), 1), 1).provides(), + nothing + ); + assert_eq!(Shifted::new(Quartic::new(3), 1).provides(), gradient); + // and so do they behind DynProblem + assert_eq!(boxed(Shifted::new(G06, 1)).provides(), values); + assert_eq!(boxed(Rotated::new(Sphere::new(2), 1)).provides(), gradient); + } + + // a shifted or rotated G06 gives its constraints' values: those of G06 at the point it's + // evaluated at + #[test] + fn a_wrapped_g06_keeps_its_constraint_values() { + let x = Reals::from(vec![20.0, 5.0]); + let shifted = Shifted::new(G06, 2); + let rotated = Rotated::new(G06, 2); + for (problem, at) in [ + (&shifted as &dyn DynFitness, shifted.unshifted(&x)), + (&rotated as &dyn DynFitness, rotated.rotated(&x)), + ] { + let mut g = [0.0; 2]; + let fitness = problem.evaluate_with(&x, &mut Extras::new(None, Some(&mut g), None)); + assert_eq!(fitness, G06.evaluate(&at)); + assert_eq!(&g[..], G06.constraints(&at).inequalities()); + } + } + + // a fitness function of (score, violation) on Real genomes, as a trait object + trait DynFitness { + fn evaluate_with(&self, x: &Reals, extras: &mut Extras<'_>) -> (f64, f64); + } + + impl> DynFitness for F { + fn evaluate_with(&self, x: &Reals, extras: &mut Extras<'_>) -> (f64, f64) { + FitnessFunction::evaluate_with(self, x, extras) + } + } + + // L-BFGS-B takes the analytic gradient through both wrappers, and needs no differences + #[test] + fn lbfgsb_runs_on_the_wrapped_gradient() -> crate::Result<()> { + let problem = Rotated::new(Shifted::new(HighConditionedElliptic::new(10), 1), 1); + let optimum = problem.optimum().expect("known"); + let lbfgsb = Lbfgsb::builder(problem.representation()) + .gradients(Gradients::Supplied) + .gradient_tolerance(1e-9) + .function_tolerance(0.0) + .minimize() + .seed(1) + .build()?; + let mut engine = Engine::new(lbfgsb, problem).stop_when(Stop::evaluations(10_000)); + let outcome = engine.run()?; + assert_eq!(outcome.stop_reason(), StopReason::Converged); + assert!(outcome.best_fitness().score().expect("valid") < 1e-16); + for (x, o) in outcome + .best_genome() + .iter() + .zip(optimum.solutions()[0].iter()) + { + assert!((x - o).abs() < 1e-8, "{x} {o}"); + } + assert_eq!(engine.algorithm().stencil_evaluations(), 0); + Ok(()) + } }