Simulated annealing, genetic algorithms, and differential evolution

A population of candidates that searches a two-dimensional landscape, one generation per step.

• Name: evolutionary_optimisers
• Field: computer_science/evolutionary_computation
• Formalism: agent-based model
• References: Raman 2021; Storn and Price 1997; Kirkpatrick, Gelatt, and Vecchi 1983

Factory

@register_model(
    field='computer_science/evolutionary_computation',
    tags=('optimisation', 'evolutionary algorithm', 'annealing', 'fitness landscape'),
    references=('Raman 2021', 'Storn and Price 1997', 'Kirkpatrick, Gelatt, and Vecchi 1983'),
    defaults={'until': 100},
    observables=_search_observables,
)
def evolutionary_optimisers(
    n_candidates: int = 40,
    *,
    algorithm: str = 'differential-evolution',
    landscape: str = 'rastrigin',
    temperature: float = 5.0,
    width: float = 0.5,
    generations: float = 100.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return a population of candidates that searches a two-dimensional landscape, one generation per step.

    A population optimiser moves the candidates through `OptimizationInteraction`, and reports `search.best`. Simulated
    annealing samples every coordinate through `SamplingInteraction` at a temperature that falls geometrically to 1e-3
    of its start over `generations` steps, and `tracker/record.best` holds the lowest energy it has reached, from the
    observable `search/lowest` that the catalog lists, so a run passes it in `observables` (Raman 2021 §6.3.3.3,
    §6.3.4).

    Args:
        n_candidates:
            The number of candidates.

        algorithm:
            `'differential-evolution'`, `'particle-swarm'`, `'evolution-strategy'`, `'covariance-adaptation'`, or
            `'annealing'`.

        landscape:
            `'rastrigin'` or `'rosenbrock'`.

        temperature:
            The initial temperature of annealing.

        width:
            The standard deviation of a move of annealing.

        generations:
            The number of steps over which annealing cools.

        random_state:
            The seed or generator of the initial candidates.

    Returns:
        system:
            The model, with the group `candidates`.

    Raises:
        ValidationError:
            The algorithm or the landscape is unknown.

    Examples:
        >>> from skcomplex.models.computer_science.evolutionary_computation import evolutionary_optimisers
        >>> evolutionary_optimisers(10).interactions[0].label
        'search'
    """
    if (algorithm not in ALGORITHMS and algorithm != ANNEALING) or landscape not in LANDSCAPES:
        names = ', '.join((*ALGORITHMS, ANNEALING))
        error_msg = (
            f'evolutionary_optimisers has algorithm={algorithm!r} and landscape={landscape!r}. Use one of {names}, '
            f'and one of {", ".join(LANDSCAPES)}.'
        )
        raise ValidationError(error_msg)
    energy, objective, bound = LANDSCAPES[landscape]
    start = np.random.default_rng(random_state).uniform(-bound, bound, (n_candidates, N_DIM))
    if algorithm == ANNEALING:
        return _annealing(start, energy, temperature, width, generations)
    candidates = GroupElement('candidates', element=SearchCandidateElement, size=n_candidates, values={'x': start})
    search = OptimizationInteraction(
        'search',
        target='candidates',
        state='x',
        objective=objective,
        bounds=(-bound, bound),
        algorithm=ALGORITHMS[algorithm],
    )
    return System(elements=[candidates], interactions=[search], time=DiscreteTime())

Components

class SearchCandidateElement(Element):
    """Candidate solution that holds a vector.

    Args:
        label:
            The label of the candidate or the group.

        x:
            The vector of the candidate.

    Examples:
        >>> from skcomplex.models.computer_science.evolutionary_computation import SearchCandidateElement
        >>> SearchCandidateElement('candidate', x=(0.0, 1.0)).x
        (0.0, 1.0)
    """

    def __init__(self, label: str, x: Annotated[npt.ArrayLike, StateSpec(vector=True)] = (0.0, 0.0)) -> None:
        super().__init__(label)
        self.x = x

class OptimizationInteraction(Interaction):
    """Move a group of candidates on an objective, one generation of a pymoo algorithm per step of discrete time.

    Every candidate holds a vector state in a box, from which the first generation starts. Every step asks the algorithm
    for candidates, evaluates the objective on their array, tells the algorithm the values, and writes the candidates of
    the generation into the members. The interaction reports `best`, the lowest value so far, and `mean`, the mean value
    of the generation, every step.

    Args:
        label:
            The label of the interaction.

        target:
            The path of the group of candidates.

        state:
            The name of the vector state of the candidates.

        objective:
            The function of the array of candidates, one row per candidate, that returns their values.

        bounds:
            The lower and the upper bounds of the box, each a number or one number per component.

        algorithm:
            The settings of the algorithm, which select it.

        phase:
            The name of the phase of the step the rule runs in, or `None` for the default phase.

    Raises:
        ValidationError:
            The algorithm is not an algorithm spec, or the objective is not callable, raised by `check`.

    Examples:
        >>> from skcomplex.core.interactions import OptimizationInteraction
        >>> objective = lambda values: (values**2).sum(axis=1)
        >>> search = OptimizationInteraction('search', target='swarm', state='x', objective=objective, bounds=(-1, 1))
        >>> search.check().label
        'search'
    """

    reports: ClassVar[tuple[str, ...]] = ('best', 'mean')

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        state: str,
        objective: Callable[[Array], npt.ArrayLike],
        bounds: tuple[npt.ArrayLike, npt.ArrayLike],
        algorithm: AlgorithmSpec = DifferentialEvolutionSpec(),
        phase: str | None = None,
    ) -> None:
        super().__init__(label)
        self.target = target
        self.state = state
        self.objective = objective
        self.bounds = bounds
        self.algorithm = algorithm
        self.phase = phase

    def _check(self) -> None:
        """Raise an error for an algorithm that is not a spec, or an objective that is not callable.

        Raises:
            ValidationError:
                The algorithm is not an algorithm spec, or the objective is not callable.
        """
        super()._check()
        if not isinstance(self.algorithm, ALGORITHMS):
            names = ', '.join(item.__name__ for item in ALGORITHMS)
            error_msg = f'The algorithm {self.algorithm!r} of {self.label!r} is no algorithm spec. Use one of {names}.'
            raise ValidationError(error_msg)
        if not callable(self.objective):
            error_msg = f'The objective of {self.label!r} must be a function of the array of candidates.'
            raise ValidationError(error_msg)

    def _differential_evolution(self, spec: DifferentialEvolutionSpec, candidates: Array) -> Algorithm:
        """Return differential evolution from the candidates."""
        return DE(
            pop_size=len(candidates),
            sampling=candidates,
            variant=spec.variant,
            F=spec.weight,
            CR=spec.crossover,
        )

    def _particle_swarm(self, spec: ParticleSwarmSpec, candidates: Array) -> Algorithm:
        """Return particle swarm optimisation from the candidates."""
        return PSO(
            pop_size=len(candidates),
            sampling=candidates,
            w=spec.inertia,
            c1=spec.cognitive,
            c2=spec.social,
            adaptive=False,
        )

    def _evolution_strategy(self, spec: EvolutionStrategySpec, candidates: Array) -> Algorithm:
        """Return an evolution strategy from the candidates."""
        return ES(pop_size=len(candidates), n_offsprings=spec.offspring, sampling=candidates)

    def _covariance_adaptation(self, spec: CovarianceAdaptationSpec, candidates: Array) -> Algorithm:
        """Return covariance matrix adaptation from the mean of the candidates."""
        return CMAES(x0=candidates.mean(axis=0), sigma=spec.sigma, pop_size=len(candidates))

    def _get_algorithm(self, candidates: Array, seed: int) -> Algorithm:
        """Return the algorithm of the spec, set up on the box of the candidates with a seed."""
        makers: dict[type, Callable[..., Algorithm]] = {
            DifferentialEvolutionSpec: self._differential_evolution,
            ParticleSwarmSpec: self._particle_swarm,
            EvolutionStrategySpec: self._evolution_strategy,
            CovarianceAdaptationSpec: self._covariance_adaptation,
        }
        size = candidates.shape[1]
        lower, upper = (np.broadcast_to(np.asarray(item, dtype=float), (size,)) for item in self.bounds)
        problem = Problem(n_var=size, n_obj=1, xl=lower, xu=upper)
        algorithm = makers[type(self.algorithm)](self.algorithm, candidates)
        algorithm.setup(problem, termination=NoTermination(), seed=seed, verbose=False)
        return algorithm

    @staticmethod
    def _members_of(generation: Array, size: int) -> Array:
        """Return one row per member from the candidates of a generation, repeating them when they are fewer."""
        return generation[np.arange(size) % len(generation)]

    def _write(self, agents: Sequence[object], candidates: Array) -> None:
        """Write the candidates into the members."""
        for agent, row in zip(agents, self._members_of(candidates, len(agents)), strict=True):
            setattr(agent, self.state, row)

    def _advance(self, model: ModelViewProtocol) -> None:
        """Advance the population by one generation, unless the algorithm has converged, and report its values."""
        agents = model.agents(self._get_role('target'))
        if 'algorithm' not in model.storage:
            candidates = np.array([getattr(agent, self.state) for agent in agents], dtype=float)
            seed = int(model.rng.integers(SEEDS))
            model.storage.update(algorithm=self._get_algorithm(candidates, seed), best=np.inf)
        algorithm = cast(Algorithm, model.storage['algorithm'])
        if not algorithm.has_next():
            model.report('best', cast(float, model.storage['best']))
            model.report('mean', cast(float, model.storage['mean']))
            return
        population = algorithm.ask()
        generation = population.get('X')
        values = np.asarray(self.objective(generation), dtype=float).ravel()
        population.set('F', values[:, None])
        algorithm.tell(infills=population)
        best = min(cast(float, model.storage['best']), float(values.min()))
        model.storage.update(best=best, mean=float(values.mean()))
        self._write(agents, generation)
        model.report('best', best)
        model.report('mean', float(values.mean()))

    def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
        """Return the rule of one generation.

        Args:
            namespace:
                The namespace of the symbols the laws can use, which the rule does not read.

        Returns:
            laws:
                The rule that advances the population.
        """
        return [RuleLaw(self._advance, phase=self.phase)]