Source code for gemseo.algos.opt.lib_scipy_global

# Copyright 2021 IRT Saint Exupéry,
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# Lesser General Public License for more details.
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
# Contributors:
#    INITIAL AUTHORS - initial API and implementation and/or initial
#                           documentation
#        :author: Matthias De Lozzo
"""A wrapper for the global optimization algorithms of the SciPy library."""
from __future__ import annotations

import sys
from dataclasses import dataclass
from typing import Any
from typing import Mapping
from typing import Union

from numpy import float64
from numpy import inf as np_inf
from numpy import int32
from numpy import isfinite
from numpy import real
from numpy.typing import NDArray
from scipy import optimize
from scipy.optimize import NonlinearConstraint

from gemseo.algos.opt.opt_lib import OptimizationAlgorithmDescription
from gemseo.algos.opt.opt_lib import OptimizationLibrary
from gemseo.algos.opt_result import OptimizationResult
from gemseo.core.mdofunctions.mdo_function import WrappedFunctionType
from gemseo.core.mdofunctions.mdo_function import WrappedJacobianType

InputType = NDArray[Union[float64, int32]]

[docs]@dataclass class SciPyGlobalAlgorithmDescription(OptimizationAlgorithmDescription): """The description of a global optimization algorithm from the SciPy library.""" library_name: str = "SciPy" website: str = ( "" )
[docs]class ScipyGlobalOpt(OptimizationLibrary): """A wrapper for the global optimization algorithms of the SciPy library.""" LIB_COMPUTE_GRAD = True LIBRARY_NAME = "SciPy" def __init__(self) -> None: # noqa:D107 super().__init__() doc = "" self.descriptions = { "DUAL_ANNEALING": SciPyGlobalAlgorithmDescription( algorithm_name="Dual annealing", description="Dual annealing", handle_integer_variables=True, internal_algorithm_name="dual_annealing", website=f"{doc}scipy.optimize.dual_annealing.html", ), "SHGO": SciPyGlobalAlgorithmDescription( algorithm_name="SHGO", description="Simplicial homology global optimization", handle_equality_constraints=True, handle_inequality_constraints=True, handle_integer_variables=True, internal_algorithm_name="shgo", positive_constraints=True, website=f"{doc}scipy.optimize.shgo.html", ), "DIFFERENTIAL_EVOLUTION": SciPyGlobalAlgorithmDescription( algorithm_name="Differential evolution", description="Differential Evolution algorithm", handle_equality_constraints=True, handle_inequality_constraints=True, handle_integer_variables=True, internal_algorithm_name="differential_evolution", website=f"{doc}scipy.optimize.differential_evolution.html", ), } # maximum function calls option passed to the algorithm self.max_func_calls = sys.maxsize def _get_options( self, max_iter: int = 999, max_time: float = 0.0, ftol_rel: float = 1e-9, ftol_abs: float = 1e-9, xtol_rel: float = 1e-9, xtol_abs: float = 1e-9, stop_crit_n_x: int = 3, workers: int = 1, updating: str = "immediate", atol: float = 0.0, init: str = "latinhypercube", recombination: float = 0.7, tol: float = 0.01, popsize: int = 15, strategy: str = "best1bin", sampling_method: str = "simplicial", niters: int = 1, n: int = 100, seed: int = 1, polish: bool = True, iters: int = 1, eq_tolerance: float = 1e-6, ineq_tolerance: float = 1e-6, normalize_design_space: bool = True, local_options: Mapping[str, Any] = None, **kwargs: Any, ) -> dict[str, Any]: # pylint: disable=W0221 r"""Set the options default values. To get the best and up-to-date information about algorithms options, go to scipy.optimize documentation: Args: max_iter: The maximum number of iterations, i.e. unique calls to f(x). max_time: The maximum runtime in seconds. The value 0 disables the cap on the runtime. ftol_rel: A stop criteria, the relative tolerance on the objective function. If abs(f(xk)-f(xk+1))/abs(f(xk))<= ftol_rel: stop. ftol_abs: A stop criteria, the absolute tolerance on the objective function. If abs(f(xk)-f(xk+1))<= ftol_rel: stop. xtol_rel: A stop criteria, the relative tolerance on the design variables. If norm(xk-xk+1)/norm(xk)<= xtol_rel: stop. xtol_abs: A stop criteria, the absolute tolerance on the design variables. If norm(xk-xk+1)<= xtol_abs: stop. stop_crit_n_x: The number of design vectors to take into account in the stopping criteria. workers: The number of processes for parallel execution. updating: The strategy to update the solution vector. If ``"immediate"``, the best solution vector is continuously updated within a single generation. With 'deferred', the best solution vector is updated once per generation. Only 'deferred' is compatible with parallelization, and the ``workers`` keyword can over-ride this option. atol: The absolute tolerance for convergence. init: Either the type of population initialization to be used or an array specifying the initial population. recombination: The recombination constant, should be in the range [0, 1]. tol: The relative tolerance for convergence. popsize: A multiplier for setting the total population size. The population has popsize * len(x) individuals. strategy: The differential evolution strategy to use. sampling_method: The method to compute the initial points. Current built in sampling method options are ``halton``, ``sobol`` and ``simplicial``. n: The number of sampling points used in the construction of the simplicial complex. seed: The seed to be used for repeatable minimizations. If None, the ``numpy.random.RandomState`` singleton is used. polish: Whether to use the L-BFGS-B algorithm to polish the best population member at the end. mutation: The mutation constant. recombination: The recombination constant. initial_temp: The initial temperature. visit: The parameter for the visiting distribution. accept: The parameter for the acceptance distribution. niters: The number of iterations used in the construction of the simplicial complex. restart_temp_ratio: During the annealing process, temperature is decreasing, when it reaches ``initial_temp * restart_temp_ratio``, the reannealing process is triggered. normalize_design_space: If True, variables are scaled in [0, 1]. eq_tolerance: The tolerance on equality constraints. ineq_tolerance: The tolerance on inequality constraints. n: The number of sampling points used in the construction of the simplicial complex. iters: The number of iterations used in the construction of the simplicial complex. local_options: The options for the local optimization algorithm, only for shgo, see scipy.optimize doc. **kwargs: The other algorithms options. """ return self._process_options( max_iter=max_iter, max_time=max_time, ftol_rel=ftol_rel, ftol_abs=ftol_abs, xtol_rel=xtol_rel, xtol_abs=xtol_abs, stop_crit_n_x=stop_crit_n_x, workers=workers, updating=updating, init=init, tol=tol, atol=atol, recombination=recombination, seed=seed, popsize=popsize, strategy=strategy, sampling_method=sampling_method, iters=iters, n=n, polish=polish, eq_tolerance=eq_tolerance, ineq_tolerance=ineq_tolerance, normalize_design_space=normalize_design_space, local_options=local_options, **kwargs, )
[docs] def iter_callback(self, x_vect: InputType) -> None: """Call the objective and constraints functions. Args: x_vect: The input data with which to call the functions. """ if self.normalize_ds: x_vect = self.problem.design_space.normalize_vect(x_vect) self.problem.objective(x_vect) for constraint in self.problem.constraints: constraint(x_vect)
[docs] def real_part_obj_fun(self, x: InputType) -> int | float: """Wrap the function and return the real part. Args: x: The values to be given to the function. Returns: The real part of the evaluation of the objective function. """ return real(self.problem.objective.func(x))
def _run(self, **options: Any) -> OptimizationResult: """Run the algorithm, to be overloaded by subclasses. Args: **options: The options for the algorithm. Returns: The optimization result. """ # remove normalization from options for algo self.normalize_ds = options.pop(self.NORMALIZE_DESIGN_SPACE_OPTION, True) # Get the normalized bounds: _, l_b, u_b = self.get_x0_and_bounds_vects(self.normalize_ds) # Replace infinite values with None: l_b = [val if isfinite(val) else None for val in l_b] u_b = [val if isfinite(val) else None for val in u_b] bounds = list(zip(l_b, u_b)) # This is required because some algorithms do not # call the objective very often when the problem # is very constrained (Power2) and OptProblem may fail # to detect the optimum. if self.problem.has_constraints(): self.problem.add_callback(self.iter_callback) if self.internal_algo_name == "dual_annealing": opt_result = optimize.dual_annealing( func=self.real_part_obj_fun, bounds=bounds, maxiter=self.max_func_calls, initial_temp=5230.0, restart_temp_ratio=2e-05, maxfun=self.max_func_calls, seed=options["seed"], ) elif self.internal_algo_name == "shgo": opt_result = optimize.shgo( func=self.real_part_obj_fun, bounds=bounds, args=(), n=options["n"], iters=options["iters"], sampling_method=options["sampling_method"], constraints=self.__get_constraints_as_scipy_dictionary(), options=options.get("local_options"), ) elif self.internal_algo_name == "differential_evolution": opt_result = optimize.differential_evolution( func=self.real_part_obj_fun, bounds=bounds, args=(), strategy=options["strategy"], maxiter=self.max_func_calls, popsize=options["popsize"], tol=options["tol"], atol=options["atol"], mutation=options.get("mutation", (0.5, 1)), recombination=options["recombination"], seed=options["seed"], polish=options["polish"], init=options["init"], updating=options["updating"], workers=options["workers"], constraints=self.__get_non_linear_constraints(), ) else: # pragma: no cover raise ValueError(f"Unknown algorithm: {self.internal_algo_name}.") return self.get_optimum_from_database(opt_result.message, opt_result.success) def __get_non_linear_constraints(self) -> tuple[NonlinearConstraint]: """Return the SciPy nonlinear constraints. Returns: The SciPy nonlinear constraints. """ eq_tolerance = self.problem.eq_tolerance constraints = [ NonlinearConstraint(constr, -eq_tolerance, eq_tolerance, jac=constr.jac) for constr in self.problem.get_eq_constraints() ] ineq_tolerance = self.problem.ineq_tolerance constraints.extend( [ NonlinearConstraint(constr, -np_inf, ineq_tolerance, jac=constr.jac) for constr in self.problem.get_ineq_constraints() ] ) return tuple(constraints) def __get_constraints_as_scipy_dictionary( self, ) -> list[dict[str, str | WrappedFunctionType | WrappedJacobianType]]: """Create the constraints to be passed to a SciPy algorithm as dictionaries. Returns: The constraints. """ return [ {"type": constraint.f_type, "fun": constraint.func, "jac": constraint.jac} for constraint in self.get_right_sign_constraints() ]