sce

Contents

sce#

The Shuffled Complex Evolution (SCE) global optimizer

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

The following functions are provided:

sce(func, x0, lb[, ub, mask, args, kwargs, ...])

Shuffled Complex Evolution algorithm for finding the minimum of a multivariate function
History

  • Written in Fortran by Q Duan, Sep 2004

  • Ported to Python by Stijn Van Hoey, 2011 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

  • random_sample(1)[0] to assure scalar, Jul 2023, Matthias Cuntz

  • call_func method to assure scalar output, Jul 2023, Matthias Cuntz

  • Require keyword names after mask, Jul 2023, Matthias Cuntz

  • Set maxn to N*log(N) formula per default, Sep 2023, Matthias Cuntz

sce(func, x0, lb, ub=None, mask=None, *, args=(), kwargs={}, sampling='half-open', maxn=0, 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='')[source]#

Shuffled Complex Evolution algorithm for finding the minimum of a multivariate function

The SCE or SCE-UA method is a general purpose global optimization, which can be used with high-dimensional problems. It was used successfully, for example, to calibrate hydrologic models with more than 50 parameters [5]. 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. The implementation here 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=scalar and args=(scalar) are not valid but should be, for example, args=(scalar,).

  • 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 one 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: 6400 + 160 * nopt * log10(nopt)).

  • 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 optimizer (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 optimizer (default: restartfile1 + ‘.txt’).

Returns:

res – 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.

Return type:

OptimizeResult

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

[5]

Cuntz M, Mai J, Zink M, Thober S, Kumar R, Schäfer D, Schrön M, Craven J, Rakovec O, Spieler D, Prykhodko V, Dalmasso G, Musuuza J, Langenberg B, Attinger S, and Samaniego L, Computationally inexpensive identification of noninformative model parameters by sequential screening. Water Resources Research 51, 6417–6441, 2015, https://doi.org/10.1002/2015WR016907

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) pairs for parameter bounds as well as setting a number of keyword parameters for the SCE algorithm is:

>>> 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
Contents