Source code for gemseo.algos.opt.lib_scipy_global

# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
#
# 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
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# 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
#    OTHER AUTHORS   - MACROSCOPIC CHANGES
"""A wrapper for the global optimization algorithms of the SciPy library."""

from __future__ import annotations

import sys
from dataclasses import dataclass
from typing import TYPE_CHECKING
from typing import Any
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.optimization_library import OptimizationAlgorithmDescription
from gemseo.algos.opt.optimization_library import OptimizationLibrary
from gemseo.utils.seeder import SEED

if TYPE_CHECKING:
    from collections.abc import Mapping

    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 = ( "https://docs.scipy.org/doc/scipy/reference/optimize.html#global-optimization" )
[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 = "https://docs.scipy.org/doc/scipy/reference/generated/" 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 = SEED, 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 = 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: https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html 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 msg = f"Unknown algorithm: {self.internal_algo_name}." raise ValueError(msg) 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() ]