from __future__ import annotations from typing import ( # noqa: UP035 Any, Callable, Iterable, TYPE_CHECKING ) import numpy as np from scipy.optimize import OptimizeResult from ._constraints import old_bound_to_new, Bounds from ._direct import direct as _direct # type: ignore if TYPE_CHECKING: import numpy.typing as npt __all__ = ['direct'] ERROR_MESSAGES = ( "Number of function evaluations done is larger than maxfun={}", "Number of iterations is larger than maxiter={}", "u[i] < l[i] for some i", "maxfun is too large", "Initialization failed", "There was an error in the creation of the sample points", "An error occurred while the function was sampled", "Maximum number of levels has been reached.", "Forced stop", "Invalid arguments", "Out of memory", ) SUCCESS_MESSAGES = ( ("The best function value found is within a relative error={} " "of the (known) global optimum f_min"), ("The volume of the hyperrectangle containing the lowest function value " "found is below vol_tol={}"), ("The side length measure of the hyperrectangle containing the lowest " "function value found is below len_tol={}"), ) def direct( func: Callable[[npt.ArrayLike, tuple[Any]], float], bounds: Iterable | Bounds, *, args: tuple = (), eps: float = 1e-4, maxfun: int | None = None, maxiter: int = 1000, locally_biased: bool = True, f_min: float = -np.inf, f_min_rtol: float = 1e-4, vol_tol: float = 1e-16, len_tol: float = 1e-6, callback: Callable[[npt.ArrayLike], None] | None = None ) -> OptimizeResult: """ Finds the global minimum of a function using the DIRECT algorithm. Parameters ---------- func : callable The objective function to be minimized. ``func(x, *args) -> float`` where ``x`` is an 1-D array with shape (n,) and ``args`` is a tuple of the fixed parameters needed to completely specify the function. bounds : sequence or `Bounds` Bounds for variables. There are two ways to specify the bounds: 1. Instance of `Bounds` class. 2. ``(min, max)`` pairs for each element in ``x``. args : tuple, optional Any additional fixed parameters needed to completely specify the objective function. eps : float, optional Minimal required difference of the objective function values between the current best hyperrectangle and the next potentially optimal hyperrectangle to be divided. In consequence, `eps` serves as a tradeoff between local and global search: the smaller, the more local the search becomes. Default is 1e-4. maxfun : int or None, optional Approximate upper bound on objective function evaluations. If `None`, will be automatically set to ``1000 * N`` where ``N`` represents the number of dimensions. Will be capped if necessary to limit DIRECT's RAM usage to app. 1GiB. This will only occur for very high dimensional problems and excessive `max_fun`. Default is `None`. maxiter : int, optional Maximum number of iterations. Default is 1000. locally_biased : bool, optional If `True` (default), use the locally biased variant of the algorithm known as DIRECT_L. If `False`, use the original unbiased DIRECT algorithm. For hard problems with many local minima, `False` is recommended. f_min : float, optional Function value of the global optimum. Set this value only if the global optimum is known. Default is ``-np.inf``, so that this termination criterion is deactivated. f_min_rtol : float, optional Terminate the optimization once the relative error between the current best minimum `f` and the supplied global minimum `f_min` is smaller than `f_min_rtol`. This parameter is only used if `f_min` is also set. Must lie between 0 and 1. Default is 1e-4. vol_tol : float, optional Terminate the optimization once the volume of the hyperrectangle containing the lowest function value is smaller than `vol_tol` of the complete search space. Must lie between 0 and 1. Default is 1e-16. len_tol : float, optional If `locally_biased=True`, terminate the optimization once half of the normalized maximal side length of the hyperrectangle containing the lowest function value is smaller than `len_tol`. If `locally_biased=False`, terminate the optimization once half of the normalized diagonal of the hyperrectangle containing the lowest function value is smaller than `len_tol`. Must lie between 0 and 1. Default is 1e-6. callback : callable, optional A callback function with signature ``callback(xk)`` where ``xk`` represents the best function value found so far. Returns ------- res : OptimizeResult The optimization result represented as a ``OptimizeResult`` object. Important attributes are: ``x`` the solution array, ``success`` a Boolean flag indicating if the optimizer exited successfully and ``message`` which describes the cause of the termination. See `OptimizeResult` for a description of other attributes. Notes ----- DIviding RECTangles (DIRECT) is a deterministic global optimization algorithm capable of minimizing a black box function with its variables subject to lower and upper bound constraints by sampling potential solutions in the search space [1]_. The algorithm starts by normalising the search space to an n-dimensional unit hypercube. It samples the function at the center of this hypercube and at 2n (n is the number of variables) more points, 2 in each coordinate direction. Using these function values, DIRECT then divides the domain into hyperrectangles, each having exactly one of the sampling points as its center. In each iteration, DIRECT chooses, using the `eps` parameter which defaults to 1e-4, some of the existing hyperrectangles to be further divided. This division process continues until either the maximum number of iterations or maximum function evaluations allowed are exceeded, or the hyperrectangle containing the minimal value found so far becomes small enough. If `f_min` is specified, the optimization will stop once this function value is reached within a relative tolerance. The locally biased variant of DIRECT (originally called DIRECT_L) [2]_ is used by default. It makes the search more locally biased and more efficient for cases with only a few local minima. A note about termination criteria: `vol_tol` refers to the volume of the hyperrectangle containing the lowest function value found so far. This volume decreases exponentially with increasing dimensionality of the problem. Therefore `vol_tol` should be decreased to avoid premature termination of the algorithm for higher dimensions. This does not hold for `len_tol`: it refers either to half of the maximal side length (for ``locally_biased=True``) or half of the diagonal of the hyperrectangle (for ``locally_biased=False``). This code is based on the DIRECT 2.0.4 Fortran code by Gablonsky et al. at https://ctk.math.ncsu.edu/SOFTWARE/DIRECTv204.tar.gz . This original version was initially converted via f2c and then cleaned up and reorganized by Steven G. Johnson, August 2007, for the NLopt project. The `direct` function wraps the C implementation. .. versionadded:: 1.9.0 References ---------- .. [1] Jones, D.R., Perttunen, C.D. & Stuckman, B.E. Lipschitzian optimization without the Lipschitz constant. J Optim Theory Appl 79, 157-181 (1993). .. [2] Gablonsky, J., Kelley, C. A Locally-Biased form of the DIRECT Algorithm. Journal of Global Optimization 21, 27-37 (2001). Examples -------- The following example is a 2-D problem with four local minima: minimizing the Styblinski-Tang function (https://en.wikipedia.org/wiki/Test_functions_for_optimization). >>> from scipy.optimize import direct, Bounds >>> def styblinski_tang(pos): ... x, y = pos ... return 0.5 * (x**4 - 16*x**2 + 5*x + y**4 - 16*y**2 + 5*y) >>> bounds = Bounds([-4., -4.], [4., 4.]) >>> result = direct(styblinski_tang, bounds) >>> result.x, result.fun, result.nfev array([-2.90321597, -2.90321597]), -78.3323279095383, 2011 The correct global minimum was found but with a huge number of function evaluations (2011). Loosening the termination tolerances `vol_tol` and `len_tol` can be used to stop DIRECT earlier. >>> result = direct(styblinski_tang, bounds, len_tol=1e-3) >>> result.x, result.fun, result.nfev array([-2.9044353, -2.9044353]), -78.33230330754142, 207 """ # convert bounds to new Bounds class if necessary if not isinstance(bounds, Bounds): if isinstance(bounds, list) or isinstance(bounds, tuple): lb, ub = old_bound_to_new(bounds) bounds = Bounds(lb, ub) else: message = ("bounds must be a sequence or " "instance of Bounds class") raise ValueError(message) lb = np.ascontiguousarray(bounds.lb, dtype=np.float64) ub = np.ascontiguousarray(bounds.ub, dtype=np.float64) # validate bounds # check that lower bounds are smaller than upper bounds if not np.all(lb < ub): raise ValueError('Bounds are not consistent min < max') # check for infs if (np.any(np.isinf(lb)) or np.any(np.isinf(ub))): raise ValueError("Bounds must not be inf.") # validate tolerances if (vol_tol < 0 or vol_tol > 1): raise ValueError("vol_tol must be between 0 and 1.") if (len_tol < 0 or len_tol > 1): raise ValueError("len_tol must be between 0 and 1.") if (f_min_rtol < 0 or f_min_rtol > 1): raise ValueError("f_min_rtol must be between 0 and 1.") # validate maxfun and maxiter if maxfun is None: maxfun = 1000 * lb.shape[0] if not isinstance(maxfun, int): raise ValueError("maxfun must be of type int.") if maxfun < 0: raise ValueError("maxfun must be > 0.") if not isinstance(maxiter, int): raise ValueError("maxiter must be of type int.") if maxiter < 0: raise ValueError("maxiter must be > 0.") # validate boolean parameters if not isinstance(locally_biased, bool): raise ValueError("locally_biased must be True or False.") def _func_wrap(x, args=None): x = np.asarray(x) if args is None: f = func(x) else: f = func(x, *args) # always return a float return np.asarray(f).item() # TODO: fix disp argument x, fun, ret_code, nfev, nit = _direct( _func_wrap, np.asarray(lb), np.asarray(ub), args, False, eps, maxfun, maxiter, locally_biased, f_min, f_min_rtol, vol_tol, len_tol, callback ) format_val = (maxfun, maxiter, f_min_rtol, vol_tol, len_tol) if ret_code > 2: message = SUCCESS_MESSAGES[ret_code - 3].format( format_val[ret_code - 1]) elif 0 < ret_code <= 2: message = ERROR_MESSAGES[ret_code - 1].format(format_val[ret_code - 1]) elif 0 > ret_code > -100: message = ERROR_MESSAGES[abs(ret_code) + 1] else: message = ERROR_MESSAGES[ret_code + 99] return OptimizeResult(x=np.asarray(x), fun=fun, status=ret_code, success=ret_code > 2, message=message, nfev=nfev, nit=nit)