Source code for pyjams.sce

#!/usr/bin/env python
"""
Shuffled-Complex-Evolution (SCE) algorithm for function minimization

This code is based on a Fortran program of Qingyun Duan (2004), ported to
Python by Stijn Van Hoey (2011). It was taken up, debugged, enhanced and is
maintained by Matthias Cuntz while at Department of Computational Hydrosystems
(CHS), Helmholtz Centre for Environmental Research - UFZ, Leipzig, Germany, and
continued while at Institut National de Recherche pour l'Agriculture,
l'Alimentation et l'Environnement (INRAE), Nancy, France.

:copyright: Copyright 2004-2023 Qingyun Duan [1]_, Stijn Van Hoey, Matthias Cuntz, see AUTHORS.rst for details.
:license: MIT License, see LICENSE for details.

References
----------
.. [1] Duan, Sorooshian and Gupta (1992) Effective and efficient global
       optimization for conceptual rainfall-runoff models, Water Resour Res 28,
       1015-1031, https://doi.org/10.1029/91WR02985

.. moduleauthor:: Matthias Cuntz

The following functions are provided:

.. autosummary::
   sce

History
    * Written in Fortran by Q Duan, Sep 2004
    * Ported to Python by Stijn Van Hoey, 2011
      https://github.com/stijnvanhoey/Optimization_SCE
    * Synchronised with enhanced Fortran version of CHS,
      Oct 2013, Matthias Cuntz
    * Added functionality to call external executable, Nov 2016, Matthias Cuntz
    * Treat NaN and Inf in function output, Nov 2016, Matthias Cuntz
    * Possibility to exclude (mask) parameters from optimisation,
      Nov 2016, Matthias Cuntz
    * Added restart possibility, Nov 2016, Matthias Cuntz
    * Return also function value of best parameter set if maxit==True,
      Nov 2016, Matthias Cuntz
    * Removed functionality to call external executable,
      Dec 2017, Matthias Cuntz
    * Print out number of function evaluations with printit=1,
      Mar 2018, Matthias Cuntz
    * Mask parameters with degenerated ranges, e.g. upper<lower bound,
      Mar 2018, Matthias Cuntz
    * Use only masked parameters in calculation of geometric range,
      Mar 2018, Matthias Cuntz
    * Removed bug that calculated the size of the complexes using all
      parameters and not only the masked parameters, Mar 2018, Matthias Cuntz
    * Fixed bug where masked parameters were always out of bounds,
      Mar 2018, Matthias Cuntz
    * Allow scalar bounds, which will be taken for all parameters,
      Mar 2018, Matthias Cuntz
    * Get number of parameters from lower boundary in SampleInputMatrix,
      Mar 2018, Matthias Cuntz
    * Define number of parameters by inquiring x0 in case of restart,
      May 2018, Matthias Cuntz
    * Removed multiplication with one hundred in criter_change,
      regarded a bug compared to Fortran code, May 2018, Matthias Cuntz
    * Removed exec command to make restart work with Python 3,
      May 2020, Matthias Cuntz
    * Use underscore before private routines, May 2020, Matthias Cuntz
    * Code refactoring, Sep 2021, Matthias Cuntz
    * Added keywords args and kwargs to pass to function,
      Apr 2022, Matthias Cuntz
    * Pass RandomState to _SampleInputMatrix, Jul 2022, Matthias Cuntz
    * Different sampling of input parameters, Jul 2022, Matthias Cuntz
    * Use helper class to pass arguments to objective function,
      Dec 2022, Matthias Cuntz
    * Copy strtobool from distutils.util because distutils is deprecated,
      Dec 2022, Matthias Cuntz
    * Include number of complexes ngs in restartfile, Dec 2022, Matthias Cuntz
    * No restart files written by default, Dec 2022, Matthias Cuntz
    * Rewrite into class SCESolver, Dec 2022, Matthias Cuntz
    * Output OptimizeResult class, Dec 2022, Matthias Cuntz
    * Polish results with L-BFGS-B, Dec 2022, Matthias Cuntz
    * Warn only if lb > ub, simply set mask if lb == ub,
      May 2023, Matthias Cuntz
    * Exit if initial population failed twice, May 2023, Matthias Cuntz

"""
import warnings
import numpy as np
from scipy.optimize import OptimizeResult, minimize
from scipy._lib._util import check_random_state
from scipy.optimize._constraints import Bounds

# ToDo:
# - write tmp/population files (of Fortran code)

__all__ = ['sce']


# extended scipy/_util.py version
class _FunctionWrapper:
    """
    Wrap user function with arguments and keywords, allowing picklability

    Parameters
    ----------
    func : callable
        Function in the form ``func(x, *args, **kwargs)``, where ``x`` are
        the parameters in the form of an iterable.
        ``args`` and ``kwargs`` are passed to the function via the usual
        unpacking operators.

    """
    def __init__(self, func, *args, **kwargs):
        self.func = func
        self.args = args
        self.kwargs = kwargs
        self.nargs = len(args)
        self.nkwargs = len(kwargs)
        if self.nkwargs == 0:
            if ( (len(self.args) == 2) and
                 isinstance(self.args[-1], dict) and
                 (len(self.args[-1]) == 0) ):
                # if kwargs={} then **kwargs={} and hence counted as args
                self.args = self.args[0]
                self.nargs = len(args)

    def __call__(self, x):
        if (self.nargs > 0) and (self.nkwargs > 0):
            return self.func(x, *self.args, **self.kwargs)
        elif (self.nargs > 0) and (self.nkwargs == 0):
            return self.func(x, *self.args)
        elif (self.nargs == 0) and (self.nkwargs > 0):
            return self.func(x, **self.kwargs)
        else:
            return self.func(x)


# from distutils.util (Python 3.11.1)
def _strtobool(val):
    """
    Convert a string representation of truth to true (1) or false (0).
    True values are 'y', 'yes', 't', 'true', 'on', and '1'; false values
    are 'n', 'no', 'f', 'false', 'off', and '0'.  Raises ValueError if
    'val' is anything else.

    """
    val = val.lower()
    if val in ('y', 'yes', 't', 'true', 'on', '1'):
        return 1
    elif val in ('n', 'no', 'f', 'false', 'off', '0'):
        return 0
    else:
        raise ValueError("invalid truth value {!r}".format(val))


[docs]def sce(func, x0, lb, ub=None, mask=None, args=(), kwargs={}, sampling='half-open', maxn=1000, kstop=10, pcento=0.0001, peps=0.001, ngs=2, npg=0, nps=0, nspl=0, mings=0, seed=None, iniflg=True, alpha=0.8, beta=0.45, maxit=False, printit=2, polish=True, restart=False, restartfile1='', restartfile2=''): """ Shuffled-Complex-Evolution algorithm for function minimization The SCE or SCE-UA method is a general purpose global optimization. The algorithm has been described in detail in Duan et al. [1]_ and [2]_. Another paper of Duan et al. [3]_ discusses how to use the method effectively. This implementation also includes the recommendations of Behrangi et al. [4]_. Parameters ---------- func : callable Function in the form ``func(x, *args, **kwargs)``, where ``x`` are the parameters in the form of an iterable. ``args`` and ``kwargs`` are passed to the function via the usual unpacking operators. x0 : array_like Parameter values with `mask==1` used in initial complex if `iniflg==1`. lb : array_like, sequence or `Bounds` Lower bounds of parameters if ``ub`` given. If `ub` is not given, then `lb` are the bounds of the parameters, either 1) as an instance of the `Bounds` class, or 2) as ``(min, max)`` pairs for each element in ``x``, defining the finite lower and upper bounds for the parameters of `func`. ub : array_like, optional Upper bounds of parameters. mask : array_like, optional Include (1, True) or exclude (0, False) parameters in minimization (default: include all parameters). The number of parameters ``nopt`` is ``sum(mask)``. args : tuple, optional Extra arguments passed to the function *func*. Note that ``args`` must be iterable. `args=int`and `args=(int)` are not valid (with int being any scalar variable) but should be `args=(int,)`. kwargs : dict, optional Extra keyword arguments passed to the function `func`. sampling : string or array_like of strings, optional Options for sampling random numbers. Options can be on of: - 'half-open': same as 'right-half-open' [lb, ub) - 'left-half-open': sample random floats in half-open interval (lb, ub] - 'right-half-open': sample random floats in half-open interval [lb, ub) - 'open': sample random floats in open interval (lb, ub) - 'log': sample half-open interval [log(lb), log(ub)), which samples better intervals spanning orders or magnitude such as `lb=1e-9` and `ub=1e-4` The default is 'half-open'. maxn : int, optional Maximum number of function evaluations allowed during minimization (without polishing) (default: 1000). kstop : int, optional Maximum number of evolution loops before convergence (default: 10). pcento : float, optional Percentage change allowed in kstop loops before convergence (default: 0.0001). peps : float, optional Value of normalised geometric range needed for convergence (default: 0.001). ngs : int, optional Number of complexes (default: 2). npg : int, optional Number of points in each complex (default: `2*nopt+1`). nps : int, optional Number of points in each sub-complex (default: `nopt+1`). mings : int, optional Minimum number of complexes required if the number of complexes is allowed to reduce as the optimization proceeds (default: `ngs`). nspl : int, optional Number of evolution steps allowed for each complex before complex shuffling (default: `2*nopt+1`). seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional If `seed` is None (or `numpy.random`), the `numpy.random.RandomState` singleton is used. If `seed` is an int, a new `RandomState` instance is used, seeded with `seed`. If `seed` is already a `Generator` or `RandomState` instance then that instance is used. Specify `seed` for repeatable results. iniflg : bool, optional If True: include initial parameters ``x0`` in initial population (default: True). alpha : float, optional Parameter for reflection of points in complex (default: 0.8). beta : float, optional Parameter for contraction of points in complex (default: 0.45). maxit : bool, optional If True: maximize instead of minimize `func` (default: False). printit : int, optional Controlling print-out (default: 2): - 0: print information for the best point of the population - 1: print information for each function evaluation - 2: no printing. The default is 2. polish : bool, optional If True (default), then `scipy.optimize.minimize` is used with the `L-BFGS-B` method to polish the result at the end, which can improve the minimization slightly. For large problems, polishing can take a long time due to the computation of the Jacobian. restart : bool, optional If True, continue from saved state in `restartfile1` and `restartfile2` (default: False). restartfile1 : str, optional Filename for saving state of array variables of `SCE` (default: ''). If `restart==True` and `restartfile1==''` then `restartfile1='sce.restart.npz'` will be taken. restartfile2 : int, optional Filename for saving state of non-array variables of `SCE` (default: `restartfile1 + '.txt'`). 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 (0) and ``message`` which describes the cause of the termination. See `OptimizeResult` for a description of other attributes. If `polish` was employed, and a lower minimum was obtained by the polishing, then OptimizeResult also contains the ``jac`` attribute. References ---------- .. [1] Duan QY, Sorooshian S, and Gupta VK, Effective and efficient global optimization for conceptual rainfall-runoff models, Water Resources Research 28(4), 1015-1031, 1992, https://doi.org/10.1029/91WR02985 .. [2] Duan QY, Gupta VK, and Sorooshian S, Shuffled Complex Evolution approach for effective and efficient global minimization, Journal of Optimization Theory and its Applications 76(3), 501-521, 1993, https://doi.org/10.1007/BF00939380 .. [3] Duan QY, Gupta VK, and Sorooshian S, Optimal use of the SCE-UA global optimization method for calibrating watershed models, Journal of Hydrology 158, 265-284, 1994, https://doi.org/10.1016/0022-1694(94)90057-4 .. [4] Behrangi A, Khakbaz B, Vrugt JA, Duan Q, and Sorooshian S Comment on "Dynamically dimensioned search algorithm for computationally efficient watershed model calibration" by Bryan A. Tolson and Christine A. Shoemaker, Water Resources Research 44, W12603, 2008, http://doi.org/10.1029/2007WR006429 Examples -------- Search the minimum of the Rosenbrock function, implemented in `rosen` in `scipy.optimize`. >>> import numpy as np >>> from scipy.optimize import rosen One has to provide the function `rosen`, an initial guess of the parameters, and the lower and upper limits of the parameters. The 2D version is: >>> lb = np.array([-5., -2.]) >>> ub = np.array([5., 8.]) >>> x0 = np.array([-2., 7.]) >>> res = sce(rosen, x0, lb, ub, seed=1, maxn=1000, printit=2) >>> print(res.nfev) 298 >>> print('{:.3f}'.format(res.fun)) 0.001 A 10-dimensional version using `(min. max)` pars for the bounds and setting a number of parameters could be: >>> nopt = 10 >>> lb = np.full(10, -5.) >>> ub = np.full(10, 5.) >>> x0 = np.full(10, 0.5) >>> res = sce(rosen, x0, zip(lb, ub), maxn=30000, kstop=10, pcento=0.0001, ... seed=12358, ngs=5, npg=5*nopt+1, nps=nopt+1, nspl=5*nopt+1, ... mings=2, iniflg=True, printit=2, alpha=0.8, beta=0.45) >>> print(res.nfev) 30228 >>> print('{:.3g}'.format(res.fun)) 3.38e-12 """ # using a context manager means that any created Pool objects are # cleared up. with SCESolver(func, x0, lb, ub=ub, mask=mask, args=args, kwargs=kwargs, sampling=sampling, maxn=maxn, kstop=kstop, pcento=pcento, ngs=ngs, npg=npg, nps=nps, nspl=nspl, mings=mings, peps=peps, seed=seed, iniflg=iniflg, alpha=alpha, beta=beta, maxit=maxit, printit=printit, polish=polish, restart=restart, restartfile1=restartfile1, restartfile2=restartfile2) as solver: ret = solver.solve() return ret
class SCESolver: """ This class implements the Shuffled-Complex-Evolution algorithm Parameters ---------- func : callable Function in the form ``func(x, *args, **kwargs)``, where ``x`` are the parameters in the form of an iterable. ``args`` and ``kwargs`` are passed to the function via the usual unpacking operators. x0 : array_like Parameter values with `mask==1` used in initial complex if `iniflg==1`. lb : array_like, sequence or `Bounds` Lower bounds of parameters if ``ub`` given. If `ub` is not given, then `lb` are the bounds of the parameters, either 1) as an instance of the `Bounds` class, or 2) as ``(min, max)`` pairs for each element in ``x``, defining the finite lower and upper bounds for the parameters of `func`. ub : array_like Upper bounds of parameters. mask : array_like, optional Include (1, True) or exclude (0, False) parameters in minimization (default: include all parameters). The number of parameters ``nopt`` is ``sum(mask)``. args : tuple, optional Extra arguments passed to the function *func*. kwargs : dict, optional Extra keyword arguments passed to the function `func`. sampling : string or array_like of strings, optional Options for sampling random numbers. Options can be on of: - 'half-open': same as 'right-half-open' [lb, ub) - 'left-half-open': sample random floats in half-open interval (lb, ub] - 'right-half-open': sample random floats in half-open interval [lb, ub) - 'open': sample random floats in open interval (lb, ub) - 'log': sample half-open interval [log(lb), log(ub)), which samples better intervals spanning orders or magnitude such as `lb=1e-9` and `ub=1e-4` The default is 'half-open'. maxn : int, optional Maximum number of function evaluations allowed during minimization (without polishing) (default: 1000). kstop : int, optional Maximum number of evolution loops before convergence (default: 10). pcento : float, optional Percentage change allowed in kstop loops before convergence (default: 0.0001). peps : float, optional Value of normalised geometric range needed for convergence (default: 0.001). ngs : int, optional Number of complexes (default: 2). npg : int, optional Number of points in each complex (default: `2*nopt+1`). nps : int, optional Number of points in each sub-complex (default: `nopt+1`). mings : int, optional Minimum number of complexes required if the number of complexes is allowed to reduce as the optimization proceeds (default: `ngs`). nspl : int, optional Number of evolution steps allowed for each complex before complex shuffling (default: `2*nopt+1`). seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional If `seed` is None (or `numpy.random`), the `numpy.random.RandomState` singleton is used. If `seed` is an int, a new `RandomState` instance is used, seeded with `seed`. If `seed` is already a `Generator` or `RandomState` instance then that instance is used. Specify `seed` for repeatable results. iniflg : bool, optional If True: include initial parameters ``x0`` in initial population (default: True). alpha : float, optional Parameter for reflection of points in complex (default: 0.8). beta : float, optional Parameter for contraction of points in complex (default: 0.45). maxit : bool, optional If True: maximize instead of minimize `func` (default: False). printit : int, optional Controlling print-out (default: 2): - 0: print information for the best point of the population - 1: print information for each function evaluation - 2: no printing. The default is 2. polish : bool, optional If True (default), then `scipy.optimize.minimize` is used with the `L-BFGS-B` method to polish the result at the end, which can improve the minimization slightly. For large problems, polishing can take a long time due to the computation of the Jacobian. restart : bool, optional If True, continue from saved state in `restartfile1` and `restartfile2` (default: False). restartfile1 : str, optional Filename for saving state of array variables of `SCE` (default: ''). If `restart==True` and `restartfile1==''` then `restartfile1='sce.restart.npz'` will be taken. restartfile2 : int, optional Filename for saving state of non-array variables of `SCE` (default: `restartfile1 + '.txt'`). """ def __init__(self, func, x0, lb, ub=None, mask=None, args=(), kwargs={}, sampling='half-open', maxn=1000, kstop=10, pcento=0.0001, ngs=2, npg=0, nps=0, nspl=0, mings=0, peps=0.001, seed=None, iniflg=True, alpha=0.8, beta=0.45, maxit=False, printit=2, polish=True, restart=False, restartfile1='', restartfile2=''): # function to minimize self.func = _FunctionWrapper(func, *args, **kwargs) # seed random number generator - will be overwritten by restart # self.rnd = np.random.RandomState(seed=seed) self.rnd = check_random_state(seed) # parameters for initial run and for restart self.sampling = sampling self.maxn = maxn self.kstop = kstop self.pcento = pcento self.peps = peps self.alpha = alpha self.beta = beta self.maxit = maxit self.printit = printit self.polish = polish self.restartfile1 = restartfile1 self.restartfile2 = restartfile2 if restart and (not self.restartfile1): self.restartfile1 = 'sce.restart.npz' if self.restartfile1 and (not self.restartfile2): self.restartfile2 = self.restartfile1 + '.txt' if not restart: # initialize SCE parameters self.nn = len(x0) self.mask = np.ones(self.nn, dtype=bool) if mask is None else mask self.nopt = np.sum(self.mask) self.ngs = ngs if ngs > 0 else 2 self.npg = npg if npg > 0 else 2 * self.nopt + 1 self.nps = nps if nps > 0 else self.nopt + 1 self.nspl = nspl if nspl > 0 else 2 * self.nopt + 1 self.mings = mings if mings > 0 else self.ngs self.npt = self.npg * self.ngs # assure lb and ub are numpy arrays if ub is None: if isinstance(lb, Bounds): self.lb = lb.lb self.ub = lb.ub else: self.lb, self.ub = zip(*lb) else: self.lb = lb self.ub = ub self.lb = np.array(self.lb) self.ub = np.array(self.ub) # same bounds for all parameters try: if len(self.lb) == 1: self.lb = np.full(self.nn, self.lb[0]) except TypeError: # could be size 0 array if lb was a scalar self.lb = np.full(self.nn, self.lb) try: if len(self.ub) == 1: self.ub = np.full(self.nn, self.ub[0]) except TypeError: self.ub = np.full(self.nn, self.ub) bound = self.ub - self.lb # degenerated bounds if np.any(bound == 0.): ii = np.where(bound == 0.)[0] self.mask[ii] = False if np.any(bound < 0.): ii = np.where(bound < 0.)[0] warnings.warn( f'sce: found lower bound lb > upper bound ub for' f' parameter(s) {ii} with lb={self.lb[ii]} and' f' ub={self.ub[ii]} => masking the parameter(s).', UserWarning, stacklevel=2) self.mask[ii] = False self.ub[ii] = self.lb[ii] # set 'large' for function runs that return NaN self.large = 0.5 * np.finfo(float).max # create an initial population to fill array x(npt, nparams) self.x = self.sample_input_matrix(self.npt) for i in range(self.npt): self.x[i, :] = np.where(self.mask, self.x[i, :], x0) if iniflg == 1: self.x[0, :] = x0 self.icall = 0 self.xf = np.zeros(self.npt) for i in range(self.npt): fuc = self.func(self.x[i, :]) self.xf[i] = -fuc if self.maxit else fuc self.icall += 1 if self.printit == 1: print(' i, f, X: ', self.icall, self.xf[i], self.x[i, :]) # redo an initial population if all runs failed if not np.any(np.isfinite(self.xf)): if self.printit < 2: print('Redo initial population because all failed') self.x = self.sample_input_matrix(self.npt) for i in range(self.npt): self.x[i, :] = np.where(self.mask, self.x[i, :], x0) for i in range(self.npt): fuc = self.func(self.x[i, :]) self.xf[i] = -fuc if self.maxit else fuc self.icall += 1 if self.printit == 1: print(' i, f, X: ', self.icall, self.xf[i], self.x[i, :]) if not np.any(np.isfinite(self.xf)): raise ValueError( 'Did not succeed to produce initial population:' ' all function evaluations failed. Give an initial' ' value x0 that works and set iniflg=True (default).') # remember large for treating of NaNs self.large = self.xf[np.isfinite(self.xf)].max() self.large = (1.1 * self.large if self.large > 0. else 0.9 * self.large) # sort the population in order of increasing function values # and report best point self.xf = np.where(np.isfinite(self.xf), self.xf, self.large) idx = np.argsort(self.xf) self.xf = self.xf[idx] self.x = self.x[idx, :] self.bestx = self.x[0, :] self.bestf = self.xf[0] # compute the normalized geometric range of the parameters self.gnrng = self.calc_gnrng() # initialise evolution loops self.nloop = 0 self.criter = [] self.criter_change = 1e+5 # save restart self.write_restartfiles() else: # if not restart self.nn = len(x0) self.read_restartfiles() def __enter__(self): return self def __exit__(self, *args): return self # def __iter__(self): # return self def __next__(self): # loop on complexes (sub-populations) for igs in range(self.ngs): k1 = np.array(range(self.npg)) k2 = k1 * self.ngs + igs # partition the population into complexes (sub-populations) cx = np.zeros((self.npg, self.nn)) cf = np.zeros((self.npg)) k1 = np.array(range(self.npg)) k2 = k1 * self.ngs + igs cx[k1, :] = self.x[k2, :] cf[k1] = self.xf[k2] # evolve sub-population igs for nspl steps: for loop in range(self.nspl): # select simplex by sampling the complex according to a # linear probability distribution lcs = np.zeros(self.nps, dtype=int) lcs[0] = 1 for k3 in range(1, self.nps): for i in range(1000): lpos = int(np.floor( self.npg + 0.5 - np.sqrt((self.npg + 0.5)**2 - self.npg * (self.npg + 1) * self.rnd.random_sample(1)) )) # check if element was already chosen idx = (lcs[0:k3] == lpos).nonzero() if idx[0].size == 0: break lcs[k3] = lpos lcs.sort() # construct the simplex s = np.zeros((self.nps, self.nn)) s = cx[lcs, :] sf = cf[lcs] # remember large for treating of NaNs self.large = cf[np.isfinite(cf)].max() self.large = (1.1 * self.large if self.large > 0. else 0.9 * self.large) snew, fnew, icall = self.cce(s, sf) # replace the worst point in simplex with the new point s[-1, :] = snew sf[-1] = fnew # self.icall += icall self.icall += icall # reinsert the simplex into the complex cx[lcs, :] = s cf[lcs] = sf # sort the complex cf = np.where(np.isfinite(cf), cf, self.large) idx = np.argsort(cf) cf = cf[idx] cx = cx[idx, :] # end of inner loop for competitive evolution of simplexes: # for loop in range(self.nspl): # i.e. end of evolve sub-population igs for nspl steps self.x[k2, :] = cx[k1, :] self.xf[k2] = cf[k1] return def calc_criter_change(self): """ Computes the percentage change in function output """ if self.nloop >= self.kstop: criter_change = np.abs(self.criter[self.nloop - 1] - self.criter[self.nloop - self.kstop]) criter_change = ( criter_change / np.maximum( 1e-15, np.mean( np.abs( self.criter[self.nloop - self.kstop:self.nloop]) ) ) ) else: criter_change = self.criter_change return criter_change def calc_gnrng(self): """ Computes the normalized geometric range of the parameters """ bound = self.ub - self.lb rrange = (np.ma.array(self.x.max(axis=0) - self.x.min(axis=0), mask=~self.mask) / np.ma.array(bound, mask=~self.mask)) gnrng = np.ma.exp(np.ma.mean(np.ma.log(rrange))) return gnrng def cce(self, s, sf): """ Generate a new point in a simplex Parameters ---------- s : 2D-array The sorted simplex in order of increasing function values sf : 1D-array Function values in increasing order Results ------- new parameter set, function value of new set, number of function evaluations performed """ # best point sb = s[0, :] # fb = sf[0] # worst point and function value sw = s[-1, :] fw = sf[-1] # centroid of the simplex excluding worst point ce = np.mean(s[:-1, :], axis=0) # attempt a reflection point snew = ce + self.alpha * (ce - sw) # sb should have initial params at mask==False snew = np.where(self.mask, snew, sb) # check if new point is outside bounds ibound = 0 if np.ma.any(np.ma.array(snew - self.lb, mask=~self.mask) < 0.): ibound = 1 if np.ma.any(np.ma.array(self.ub - snew, mask=~self.mask) < 0.): ibound = 2 if ibound >= 1: snew = self.sample_input_matrix(1)[0, :] snew = np.where(self.mask, snew, sb) icall = 0 # calc function for reflection point fuc = self.func(snew) fnew = -fuc if self.maxit else fuc icall += 1 if self.printit == 1: print(' i, f, X: ', self.icall + icall, fnew, snew) # reflection failed: attempt a contraction point if fnew > fw: snew = sw + self.beta * (ce - sw) snew = np.where(self.mask, snew, sb) fuc = self.func(snew) fnew = -fuc if self.maxit else fuc icall += 1 if self.printit == 1: print(' i, f, X: ', self.icall + icall, fnew, snew) # both reflection and contraction have failed, attempt a random point if fnew > fw: snew = self.sample_input_matrix(1)[0, :] snew = np.where(self.mask, snew, sb) fuc = self.func(snew) fnew = -fuc if self.maxit else fuc icall += 1 if self.printit == 1: print(' i, f, X: ', self.icall + icall, fnew, snew) # end of cce return snew, fnew, icall def check_number_calls(self): """ Check for maximum number of function calls reached """ if self.icall < self.maxn: return True else: if self.printit < 2: print(f'Optimisation terminated because trial number' f' {self.maxn} reached maximum number of trials' f' {self.icall}.') return False def check_criter_change(self): """ Check if percentage change in function output is below pcento """ if self.criter_change < self.pcento: if self.printit < 2: print(f'The best point has improved by less then {self.pcento}' f' in the last {self.kstop} loops.') return True else: return False def check_geometric_range(self): """ Check if normalized geometric range of the parameters is below peps """ self.gnrng = self.calc_gnrng() if self.gnrng < self.peps: if self.printit < 2: print(f'The population has converged to a small parameter' f' space {self.gnrng} (<{self.peps}).') return True else: return False def read_restartfiles(self): """ Read from restart files """ with open(self.restartfile1, 'rb') as p1: pp = np.load(p1) self.lb = pp['lb'] self.ub = pp['ub'] self.mask = pp['mask'] self.criter = pp['criter'] self.x = pp['x'] self.xf = pp['xf'] rs2 = pp['rs2'] with open(self.restartfile2, 'r') as p2: (self.ngs, self.nopt, self.npg, self.nps, self.nspl, self.mings, self.npt, self.nloop, self.icall, rs3, rs4) = [ int(inn) for inn in p2.readline().rstrip().split(',') ] (self.gnrng, self.criter_change, self.bestf, rs5) = [ float(inn) for inn in p2.readline().rstrip().split(',') ] self.maxit = bool(_strtobool(p2.readline().rstrip())) rs1 = p2.readline().rstrip() self.rnd.set_state((rs1, rs2, rs3, rs4, rs5)) def write_restartfiles(self): """ Write restart files if demanded """ if self.restartfile1: # compute the normalized geometric range of the parameters self.gnrng = self.calc_gnrng() # compute percentage change in function output self.criter_change = self.calc_criter_change() # get the current state of the random number generator rs1, rs2, rs3, rs4, rs5 = self.rnd.get_state() # only arrays into savez_compressed - restartfile1 np.savez_compressed(self.restartfile1, lb=self.lb, ub=self.ub, mask=self.mask, criter=self.criter, x=self.x, xf=self.xf, bestx=self.bestx, rs2=rs2) # scalars into text file - restartfile2 with open(self.restartfile2, 'w') as p: iout = [self.ngs, self.nopt, self.npg, self.nps, self.nspl, self.mings, self.npt, self.nloop, self.icall, rs3, rs4] iform = '{:d},' * (len(iout) - 1) + '{:d}' print(iform.format(*iout), file=p) fout = [self.gnrng, self.criter_change, self.bestf, rs5] fform = '{:.15g},' * (len(fout) - 1) + '{:.15g}' print(fform.format(*fout), file=p) print(self.maxit, file=p) print(rs1, file=p) return def sample_input_matrix(self, npt=0): """ Create input parameter matrix (npt, npars) for npt simulations and npars parameters with bounds lb and ub Returns ------- parameter matrix (npt, npars) with new parameter samples """ if npt < 1: npt = self.npt npars = len(self.lb) if isinstance(self.sampling, str): isampling = [self.sampling] * npars else: isampling = self.sampling assert len(isampling) == npars, ( f'sampling must be string or list of strings' f' with {npars} entries') x = np.zeros((npt, npars)) bound = self.ub - self.lb for i in range(npt): irnd = self.rnd.random_sample(npars) for j in range(npars): opt = isampling[j].lower() if (opt == 'half-open') or (opt == 'right-half-open'): x[i, j] = self.lb[j] + irnd[j] * bound[j] elif opt == 'left-half-open': irnd[j] = 1. - irnd[j] x[i, j] = self.lb[j] + irnd[j] * bound[j] elif opt == 'open': iirnd = irnd[j] while not (iirnd > 0.): iirnd = self.rnd.random_sample(1) x[i, j] = self.lb[j] + iirnd * bound[j] elif opt == 'log': # x must be > 0. for ln(x) xshift = 0. if (self.lb[j] * self.ub[j]) < 0.: # lb < 0 and ub > 0 -> shift both > 0 xshift = 2. * np.maximum(np.abs(self.lb[j]), np.abs(self.ub[j])) elif (self.lb[j] * self.ub[j]) == 0.: if self.lb[j] == 0.: # lb == 0 and ub > 0 -> shift to [ub, 2*ub) xshift = self.ub[j] if self.ub[j] == 0.: # lb < 0 and ub == 0 -> shift to [-lb, -2*lb) > 0. xshift = -2. * self.lb[j] elif (self.lb[j] < 0.) and (self.ub[j] < 0.): # lb < 0 and ub < 0 -> shift both > 0 xshift = 2. * np.maximum(np.abs(self.lb[j]), np.abs(self.ub[j])) lnlb = np.log(self.lb[j] + xshift) lnub = np.log(self.ub[j] + xshift) x[i, j] = np.exp(lnlb + irnd[j] * (lnub - lnlb)) - xshift else: raise ValueError(f'unknown sampling option' f' {isampling[j]}.\n' f'Known samplings are: half-open,' f' left-half-open, right-half-open, open,' f' log') return x def solve(self): while (self.check_number_calls() and (not self.check_geometric_range()) and (not self.check_criter_change())): self.nloop += 1 # loop on complexes (sub-populations) next(self) # shuffle the complexes idx = np.argsort(self.xf) self.xf = self.xf[idx] self.x = self.x[idx, :] # record best and worst points self.bestx = self.x[0, :] self.bestf = self.xf[0] if self.printit < 2: print(f'\nEvolution loop {self.nloop}, trials {self.icall}.' f' Best f: {self.bestf:.2f}, worst f:' f' {self.xf[-1]:.2f}\n' f' best x:\n', self.bestx, '\n') self.criter = np.append(self.criter, self.bestf) self.criter_change = self.calc_criter_change() # save restart self.write_restartfiles() # finish up if self.printit < 2: print('Search stopped at trial number {0:d} with normalized' ' geometric range {1:f}. '.format(self.icall, self.gnrng)) print('The best point has improved by {:f} in the last' ' {:d} loops.'.format(self.criter_change, self.kstop)) if not self.check_number_calls(): success = False status = 2 message = f'Reached maximum number of trials {self.maxn}' if self.check_criter_change(): success = True status = 1 message = (f'Best point improved less than {self.pcento}' f' in last {self.kstop} loops') if self.check_geometric_range(): success = True status = 0 message = (f'Normalized geometric range of parameters' f' {self.gnrng} < {self.peps}') if self.maxit: self.bestf *= -1. sce_result = OptimizeResult( x=self.bestx, success=success, status=status, message=message, fun=self.bestf, nfev=self.icall, nit=self.nloop) # polish results, only works if no mask if self.polish and (len(self.lb) == self.nopt): polish_method = 'L-BFGS-B' result = minimize(self.func, np.copy(sce_result.x), method=polish_method, bounds=zip(self.lb, self.ub)) self.icall += result.nfev sce_result.nfev = self.icall # Polishing solution is only accepted if there is an improvement in # the cost function, the polishing was successful and the solution # lies within the bounds. if ( (result.fun < sce_result.fun) and result.success and np.all(result.x <= self.ub) and np.all(self.lb <= result.x) ): sce_result.fun = result.fun sce_result.x = result.x sce_result.jac = result.jac return sce_result if __name__ == '__main__': import doctest doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE) # from pyjams.functions import ackley, griewank, goldstein_price # from pyjams.functions import rastrigin, rosenbrock, six_hump_camelback # from scipy.optimize import Bounds # maxn = 10000 # seed = 12345 # """ # This is the Ackley Function # Global Optimum (n>=2): 0.0 at origin # """ # npara = 10 # lb = -10 * np.ones(npara) # ub = 10 * np.ones(npara) # x0 = np.random.rand(npara) * 10. # res = sce(ackley, x0, lb, ub, maxn=maxn, # restartfile1='', seed=seed, polish=False) # print('Ackley unpolished ', res.x, res.fun) # # print(' ', res.message) # res = sce(ackley, x0, lb, ub, maxn=maxn, # restartfile1='', seed=seed) # print('Ackley ', res.x, res.fun) # # print(' ', res.message) # """ # This is the Griewank Function (2-D or 10-D) # Bound: X(i)=[-600,600], for i=1,2,...,10 !for visualization only 2! # Global Optimum: 0, at origin # """ # npara = 10 # lb = -600 * np.ones(npara) # ub = 600 * np.ones(npara) # x0 = np.random.rand(npara) * 600. # res = sce(griewank, x0, lb, ub, maxn=maxn, # restartfile1='', seed=seed) # print('Griewank ', res.x, res.fun) # # print(' ', res.message) # """ # This is the Goldstein-Price Function # Bound X1=[-2,2], X2=[-2,2] # Global Optimum: 3.0,(0.0,-1.0) # """ # npara = 2 # lb = -2 * np.ones(npara) # ub = 2 * np.ones(npara) # x0 = np.random.rand(npara) * 2. # res = sce(goldstein_price, x0, zip(lb, ub), maxn=maxn, # restartfile1='', seed=seed) # print('Goldstein ', res.x, res.fun) # # print(' ', res.message) # """ # This is the Rastrigin Function # Bound: X1=[-1,1], X2=[-1,1] # Global Optimum: -2, (0,0) # """ # npara = 2 # lb = -1 * np.ones(npara) # ub = 1 * np.ones(npara) # x0 = np.random.rand(npara) * 1. # bounds = Bounds(lb, ub) # res = sce(rastrigin, x0, bounds, maxn=maxn, # restartfile1='', seed=seed) # print('Rastrigin ', res.x, res.fun) # # print(' ', res.message) # """ # This is the Rosenbrock Function # Bound: X1=[-5,5], X2=[-2,8]; Global Optimum: 0,(1,1) # lb=[-5 -5]; ub=[5 5]; x0=[1 1]; # """ # npara = 2 # lb = -2 * np.ones(npara) # ub = 5 * np.ones(npara) # x0 = np.random.rand(npara) * 2. # res = sce(rosenbrock, x0, lb, ub, maxn=maxn, # seed=seed) # print('Rosenbrock ', res.x, res.fun) # # print(' ', res.message) # """ # This is the Six-hump Camelback Function. # Bound: X1=[-5,5], X2=[-5,5] # True Optima: -1.031628453489877, (-0.08983,0.7126), (0.08983,-0.7126) # """ # npara = 2 # lb = -5 * np.ones(npara) # ub = 5 * np.ones(npara) # x0 = np.random.rand(npara) * 5. # res = sce(six_hump_camelback, x0, lb, ub, maxn=maxn, # seed=seed) # print('Six_hump_camelback ', res.x, res.fun) # # print(' ', res.message) # # Restart # maxn = 500 # npara = 2 # lb = -2 * np.ones(npara) # ub = 5 * np.ones(npara) # x0 = np.zeros(npara) # res = sce(rosenbrock, x0, lb, ub, maxn=maxn, # restartfile1='', # seed=seed, restart=False) # print('Rosenbrock Reference - ', res.x, res.fun) # # print(' ', res.message) # res = sce(rosenbrock, x0, lb, ub, maxn=maxn // 2, # seed=seed, restart=False, # restartfile1='sce.restart.npz') # print('Rosenbrock Restart 1 - ', res.x, res.fun) # res = sce(rosenbrock, x0, lb, ub, maxn=maxn, # seed=seed, restart=True, # restartfile1='sce.restart.npz') # print('Rosenbrock Restart 2 - ', res.x, res.fun) # # print(' ', res.message) # # Clean restart # import os # os.remove('sce.restart.npz') # os.remove('sce.restart.npz.txt') # # different sampling # npara = 2 # lb = -2 * np.ones(npara) # ub = 5 * np.ones(npara) # x0 = np.random.rand(npara) * 2. # # sampling is str # for sampling in ['half-open', 'left-half-open', # 'right-half-open', 'open', 'log']: # res = sce(rosenbrock, x0, lb, ub, sampling=sampling, # maxn=maxn, seed=seed) # print('Rosenbrock0 ', sampling, res.x, res.fun) # # print(' ', res.message) # # sampling is list # sampling = ['half-open', 'log'] # res = sce(rosenbrock, x0, lb, ub, sampling=sampling, # maxn=maxn, seed=seed) # print('Rosenbrock1 ', sampling, res.x, res.fun) # # print(' ', res.message) # # bounds include 0 # lb = [0., -2.] # ub = [5., 8.] # x0 = np.random.rand(npara) * 2. # for sampling in ['half-open', 'left-half-open', # 'right-half-open', 'open', 'log']: # res = sce(rosenbrock, x0, lb, ub, sampling=sampling, # maxn=maxn, seed=seed) # print('Rosenbrock2 ', sampling, res.x, res.fun) # # print(' ', res.message) # # bounds all > 0 # lb = [0.1, 0.1] # ub = [5., 8.] # x0 = np.random.rand(npara) * 2. # for sampling in ['half-open', 'left-half-open', # 'right-half-open', 'Open', 'log']: # res = sce(rosenbrock, x0, lb, ub, sampling=sampling, # maxn=maxn, restartfile1='', seed=seed) # print('Rosenbrock3 ', sampling, res.x, res.fun) # # print(' ', res.message)