Source code for iddefix.solvers

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Created on Sat Dec  5 16:33:41 2020

@author: sjoly
@contributor: edelafue, babreufig
"""

from typing import Any, Callable

import numpy as np
from scipy.optimize import differential_evolution
from tqdm import tqdm

ParameterBounds = list[tuple[float, float]]
MinimizationFunction = Callable[[np.ndarray], float]


[docs]class ProgressBarCallback: def __init__(self, max_generations: int, desc: str = "Optimization") -> None: self.max_generations = max_generations self.current_generation = 0 self.pbar = tqdm(total=max_generations, desc=desc, unit="gen") def __call__(self, *args: Any) -> bool | None: # --- scipy differential_evolution: (xk, convergence) --- if len(args) == 2 and not hasattr(args[0], "evaluator"): xk, convergence = args self.current_generation += 1 self.pbar.update(1) self.pbar.set_postfix({"conv": f"{(100 * (convergence)):6.1f} %"}) # optional early stopping if convergence >= 1.0: self.pbar.close() return True # --- pymoo minimize: (algorithm) --- elif len(args) == 1: algorithm = args[0] self.current_generation += 1 self.pbar.update(1) # compute convergence metric F = np.atleast_1d(algorithm.pop.get("F")) gap = float(np.mean(F) - np.min(F)) self._gap0 = getattr(self, "_gap0", gap if gap > 0 else 1.0) convergence = float(np.clip(1.0 - gap / (self._gap0 + 1e-12), 0.0, 1.0)) self.pbar.set_postfix({"conv": f"{(100 * (convergence)):6.1f} %"})
[docs] def close(self) -> None: self.pbar.close()
[docs]def stop_criterion(solver: Any) -> float: """Stopping criterion used in SciPy's differential evolution. Based on the criterion used in ``scipy.optimize.differential_evolution``: it checks the ratio of the spread of the population fitness compared to its average. If most individuals converge to the same solution, an optimal solution is considered found. """ population_cost = np.vstack(solver)[:, 1] population_mean, population_std = ( np.mean(population_cost), np.std(population_cost), ) criterion = population_std / np.abs(population_mean) return criterion
[docs]class Solvers:
[docs] def run_scipy_solver( parameterBounds: ParameterBounds, minimization_function: MinimizationFunction, maxiter: int = 2000, popsize: int = 150, mutation: tuple[float, float] = (0.1, 0.5), crossover_rate: float = 0.8, tol: float = 0.01, **kwargs, ) -> tuple[np.ndarray, str]: """Run SciPy's ``differential_evolution`` to minimize a function. All arguments are detailed in the SciPy documentation: https://docs.scipy.org/doc/scipy/reference/generated/ scipy.optimize.differential_evolution.html. Setting ``workers=-1`` uses all available CPUs. Default parameters for the DE algorithm are taken from https://www.mdpi.com/2227-7390/9/4/427. Args: parameterBounds: List of ``(min, max)`` bounds for each parameter. minimization_function: Function to be minimized. maxiter: Maximum number of iterations to run the solver. popsize: Population size for the differential evolution algorithm. mutation: Tuple of two floats representing mutation factors. crossover_rate: Crossover rate for the differential evolution algorithm. tol: Tolerance for convergence. Returns: Tuple ``(solution, message)`` with the best solution and a status message. """ pbar = ProgressBarCallback(maxiter, desc="Differential Evolution") result = differential_evolution( minimization_function, parameterBounds, popsize=popsize, tol=tol, maxiter=maxiter, mutation=mutation, recombination=crossover_rate, polish=False, init="latinhypercube", strategy="rand1bin", callback=pbar, updating="deferred", workers=-1, **kwargs, ) pbar.close() solution, message = result.x, result.message return solution, message
[docs] def run_pyfde_solver( parameterBounds: ParameterBounds, minimization_function: MinimizationFunction, maxiter: int = 2000, popsize: int = 150, mutation: float = 0.45, crossover_rate: float = 0.8, tol: float = 0.01, **kwargs, ) -> tuple[np.ndarray, str]: """ Runs the pyfde ClassicDE solver to minimize a given function. Args: parameterBounds: A list of tuples representing the bounds for each parameter. minimization_function: The function to be minimized. maxiter: The maximum number of iterations to run the solver for. popsize: The population size for the differential evolution algorithm. mutation: A tuple of two floats representing the mutation factors. crossover_rate: The crossover rate for the differential evolution algorithm. tol: The tolerance for convergence. Returns: A tuple containing: - The solution found by the solver. - A message indicating the solver's status. """ try: from pyfde import ClassicDE except ImportError: raise ImportError( "Please install the pyfde package to use the pyfde solvers." ) solver = ClassicDE( minimization_function, n_dim=len(parameterBounds), n_pop=popsize * len(parameterBounds), limits=parameterBounds, minimize=True, ) solver.cr, solver.f = crossover_rate, np.mean(np.atleast_1d(mutation)) for i in tqdm(range(maxiter)): best, _ = solver.run(n_it=1) if stop_criterion(solver) < tol: break solution, message = ( best, "Convergence achieved" if i < maxiter else "Maximum iterations reached", ) return solution, message
[docs] def run_pyfde_jade_solver( parameterBounds: ParameterBounds, minimization_function: MinimizationFunction, maxiter: int = 2000, popsize: int = 150, tol: float = 0.01, **kwargs, ) -> tuple[np.ndarray, str]: """ Runs the pyfde JADE solver to minimize a given function. Args: parameterBounds: A list of tuples representing the bounds for each parameter. minimization_function: The function to be minimized. maxiter: The maximum number of iterations to run the solver for. popsize: The population size for the differential evolution algorithm. tol: The tolerance for convergence. Returns: A tuple containing: - The solution found by the solver. - A message indicating the solver's status. """ try: from pyfde import JADE except ImportError: raise ImportError( "Please install the pyfde package to use the pyfde solvers." ) solver = JADE( minimization_function, n_dim=len(parameterBounds), n_pop=popsize * len(parameterBounds), limits=parameterBounds, minimize=True, ) for i in tqdm(range(maxiter)): best, _ = solver.run(n_it=1) if stop_criterion(solver) < tol: break solution, message = ( best, "Convergence achieved" if i < maxiter else "Maximum iterations reached", ) return solution, message
[docs] def run_pymoo_cmaes_solver( parameterBounds: ParameterBounds, minimization_function: MinimizationFunction, sigma: float = 0.1, maxiter: int | None = None, # default: 100 + 150 * (N+3)**2 // popsize**0.5 popsize: int | None = None, # defaul: 4 + int(3 * np.log(len(parameterBounds))) verbose: bool = False, **kwargs, ) -> tuple[np.ndarray, str, Any]: """ Runs the pymoo CMAES solver to minimize a given function. Args: parameterBounds: A list of tuples representing the bounds for each parameter. minimization_function: The function to be minimized. sigma: The initial standard deviation for the CMA-ES algorithm. maxiter: The maximum number of iterations to run the solver for. popsize: The population size for the differential evolution algorithm. tol: The tolerance for convergence. Returns: A tuple containing: - The solution found by the solver. - A message indicating the solver's status. """ try: from pymoo.algorithms.soo.nonconvex.cmaes import CMAES from pymoo.core.problem import Problem from pymoo.optimize import minimize except ImportError: ImportError("""Please install the pymoo package to use the CMA-ES solver: >>> pip install pymoo """) class OptimizationProblem(Problem): def __init__( self, objective_function: MinimizationFunction, n_var: int, n_obj: int, xl: list[float], xu: list[float], ) -> None: super().__init__(n_var=n_var, n_obj=n_obj, xl=xl, xu=xu) self.objective_function = objective_function def _evaluate(self, x: np.ndarray, out: dict[str, Any]) -> None: out["F"] = [self.objective_function(xi) for xi in x] problem = OptimizationProblem( objective_function=minimization_function, n_var=len(parameterBounds), n_obj=1, xl=[bound[0] for bound in parameterBounds], xu=[bound[1] for bound in parameterBounds], ) # Calculate mean of parameter bounds as starting point x0 = np.mean(parameterBounds, axis=1) solver = CMAES( x0=x0, sigma=sigma, popsize=popsize, maxiter=maxiter, seed=42, restarts=3, restart_from_best=True, **kwargs, ) if not verbose: cb = ProgressBarCallback(maxiter, desc="CMA-ES evolution") res = minimize( problem, solver, seed=42, callback=cb, save_history=True, ) else: res = minimize( problem, solver, seed=42, verbose=True, save_history=True, ) if not verbose: cb.close() solution = res.X message = ( "Convergence achieved" if res.algorithm.n_gen < maxiter else "Maximum iterations reached" ) return solution, message, res