Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_differentialevolution.py: 11%
443 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""
2differential_evolution: The differential evolution global optimization algorithm
3Added by Andrew Nelson 2014
4"""
5import warnings
7import numpy as np
8from scipy.optimize import OptimizeResult, minimize
9from scipy.optimize._optimize import _status_message
10from scipy._lib._util import check_random_state, MapWrapper, _FunctionWrapper
12from scipy.optimize._constraints import (Bounds, new_bounds_to_old,
13 NonlinearConstraint, LinearConstraint)
14from scipy.sparse import issparse
16__all__ = ['differential_evolution']
19_MACHEPS = np.finfo(np.float64).eps
22def differential_evolution(func, bounds, args=(), strategy='best1bin',
23 maxiter=1000, popsize=15, tol=0.01,
24 mutation=(0.5, 1), recombination=0.7, seed=None,
25 callback=None, disp=False, polish=True,
26 init='latinhypercube', atol=0, updating='immediate',
27 workers=1, constraints=(), x0=None, *,
28 integrality=None, vectorized=False):
29 """Finds the global minimum of a multivariate function.
31 The differential evolution method [1]_ is stochastic in nature. It does
32 not use gradient methods to find the minimum, and can search large areas
33 of candidate space, but often requires larger numbers of function
34 evaluations than conventional gradient-based techniques.
36 The algorithm is due to Storn and Price [2]_.
38 Parameters
39 ----------
40 func : callable
41 The objective function to be minimized. Must be in the form
42 ``f(x, *args)``, where ``x`` is the argument in the form of a 1-D array
43 and ``args`` is a tuple of any additional fixed parameters needed to
44 completely specify the function. The number of parameters, N, is equal
45 to ``len(x)``.
46 bounds : sequence or `Bounds`
47 Bounds for variables. There are two ways to specify the bounds:
48 1. Instance of `Bounds` class.
49 2. ``(min, max)`` pairs for each element in ``x``, defining the finite
50 lower and upper bounds for the optimizing argument of `func`.
51 The total number of bounds is used to determine the number of
52 parameters, N.
53 args : tuple, optional
54 Any additional fixed parameters needed to
55 completely specify the objective function.
56 strategy : str, optional
57 The differential evolution strategy to use. Should be one of:
59 - 'best1bin'
60 - 'best1exp'
61 - 'rand1exp'
62 - 'randtobest1exp'
63 - 'currenttobest1exp'
64 - 'best2exp'
65 - 'rand2exp'
66 - 'randtobest1bin'
67 - 'currenttobest1bin'
68 - 'best2bin'
69 - 'rand2bin'
70 - 'rand1bin'
72 The default is 'best1bin'.
73 maxiter : int, optional
74 The maximum number of generations over which the entire population is
75 evolved. The maximum number of function evaluations (with no polishing)
76 is: ``(maxiter + 1) * popsize * N``
77 popsize : int, optional
78 A multiplier for setting the total population size. The population has
79 ``popsize * N`` individuals. This keyword is overridden if an
80 initial population is supplied via the `init` keyword. When using
81 ``init='sobol'`` the population size is calculated as the next power
82 of 2 after ``popsize * N``.
83 tol : float, optional
84 Relative tolerance for convergence, the solving stops when
85 ``np.std(pop) <= atol + tol * np.abs(np.mean(population_energies))``,
86 where and `atol` and `tol` are the absolute and relative tolerance
87 respectively.
88 mutation : float or tuple(float, float), optional
89 The mutation constant. In the literature this is also known as
90 differential weight, being denoted by F.
91 If specified as a float it should be in the range [0, 2].
92 If specified as a tuple ``(min, max)`` dithering is employed. Dithering
93 randomly changes the mutation constant on a generation by generation
94 basis. The mutation constant for that generation is taken from
95 ``U[min, max)``. Dithering can help speed convergence significantly.
96 Increasing the mutation constant increases the search radius, but will
97 slow down convergence.
98 recombination : float, optional
99 The recombination constant, should be in the range [0, 1]. In the
100 literature this is also known as the crossover probability, being
101 denoted by CR. Increasing this value allows a larger number of mutants
102 to progress into the next generation, but at the risk of population
103 stability.
104 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
105 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
106 singleton is used.
107 If `seed` is an int, a new ``RandomState`` instance is used,
108 seeded with `seed`.
109 If `seed` is already a ``Generator`` or ``RandomState`` instance then
110 that instance is used.
111 Specify `seed` for repeatable minimizations.
112 disp : bool, optional
113 Prints the evaluated `func` at every iteration.
114 callback : callable, `callback(xk, convergence=val)`, optional
115 A function to follow the progress of the minimization. ``xk`` is
116 the best solution found so far. ``val`` represents the fractional
117 value of the population convergence. When ``val`` is greater than one
118 the function halts. If callback returns `True`, then the minimization
119 is halted (any polishing is still carried out).
120 polish : bool, optional
121 If True (default), then `scipy.optimize.minimize` with the `L-BFGS-B`
122 method is used to polish the best population member at the end, which
123 can improve the minimization slightly. If a constrained problem is
124 being studied then the `trust-constr` method is used instead. For large
125 problems with many constraints, polishing can take a long time due to
126 the Jacobian computations.
127 init : str or array-like, optional
128 Specify which type of population initialization is performed. Should be
129 one of:
131 - 'latinhypercube'
132 - 'sobol'
133 - 'halton'
134 - 'random'
135 - array specifying the initial population. The array should have
136 shape ``(S, N)``, where S is the total population size and N is
137 the number of parameters.
138 `init` is clipped to `bounds` before use.
140 The default is 'latinhypercube'. Latin Hypercube sampling tries to
141 maximize coverage of the available parameter space.
143 'sobol' and 'halton' are superior alternatives and maximize even more
144 the parameter space. 'sobol' will enforce an initial population
145 size which is calculated as the next power of 2 after
146 ``popsize * N``. 'halton' has no requirements but is a bit less
147 efficient. See `scipy.stats.qmc` for more details.
149 'random' initializes the population randomly - this has the drawback
150 that clustering can occur, preventing the whole of parameter space
151 being covered. Use of an array to specify a population could be used,
152 for example, to create a tight bunch of initial guesses in an location
153 where the solution is known to exist, thereby reducing time for
154 convergence.
155 atol : float, optional
156 Absolute tolerance for convergence, the solving stops when
157 ``np.std(pop) <= atol + tol * np.abs(np.mean(population_energies))``,
158 where and `atol` and `tol` are the absolute and relative tolerance
159 respectively.
160 updating : {'immediate', 'deferred'}, optional
161 If ``'immediate'``, the best solution vector is continuously updated
162 within a single generation [4]_. This can lead to faster convergence as
163 trial vectors can take advantage of continuous improvements in the best
164 solution.
165 With ``'deferred'``, the best solution vector is updated once per
166 generation. Only ``'deferred'`` is compatible with parallelization or
167 vectorization, and the `workers` and `vectorized` keywords can
168 over-ride this option.
170 .. versionadded:: 1.2.0
172 workers : int or map-like callable, optional
173 If `workers` is an int the population is subdivided into `workers`
174 sections and evaluated in parallel
175 (uses `multiprocessing.Pool <multiprocessing>`).
176 Supply -1 to use all available CPU cores.
177 Alternatively supply a map-like callable, such as
178 `multiprocessing.Pool.map` for evaluating the population in parallel.
179 This evaluation is carried out as ``workers(func, iterable)``.
180 This option will override the `updating` keyword to
181 ``updating='deferred'`` if ``workers != 1``.
182 This option overrides the `vectorized` keyword if ``workers != 1``.
183 Requires that `func` be pickleable.
185 .. versionadded:: 1.2.0
187 constraints : {NonLinearConstraint, LinearConstraint, Bounds}
188 Constraints on the solver, over and above those applied by the `bounds`
189 kwd. Uses the approach by Lampinen [5]_.
191 .. versionadded:: 1.4.0
193 x0 : None or array-like, optional
194 Provides an initial guess to the minimization. Once the population has
195 been initialized this vector replaces the first (best) member. This
196 replacement is done even if `init` is given an initial population.
197 ``x0.shape == (N,)``.
199 .. versionadded:: 1.7.0
201 integrality : 1-D array, optional
202 For each decision variable, a boolean value indicating whether the
203 decision variable is constrained to integer values. The array is
204 broadcast to ``(N,)``.
205 If any decision variables are constrained to be integral, they will not
206 be changed during polishing.
207 Only integer values lying between the lower and upper bounds are used.
208 If there are no integer values lying between the bounds then a
209 `ValueError` is raised.
211 .. versionadded:: 1.9.0
213 vectorized : bool, optional
214 If ``vectorized is True``, `func` is sent an `x` array with
215 ``x.shape == (N, S)``, and is expected to return an array of shape
216 ``(S,)``, where `S` is the number of solution vectors to be calculated.
217 If constraints are applied, each of the functions used to construct
218 a `Constraint` object should accept an `x` array with
219 ``x.shape == (N, S)``, and return an array of shape ``(M, S)``, where
220 `M` is the number of constraint components.
221 This option is an alternative to the parallelization offered by
222 `workers`, and may help in optimization speed by reducing interpreter
223 overhead from multiple function calls. This keyword is ignored if
224 ``workers != 1``.
225 This option will override the `updating` keyword to
226 ``updating='deferred'``.
227 See the notes section for further discussion on when to use
228 ``'vectorized'``, and when to use ``'workers'``.
230 .. versionadded:: 1.9.0
232 Returns
233 -------
234 res : OptimizeResult
235 The optimization result represented as a `OptimizeResult` object.
236 Important attributes are: ``x`` the solution array, ``success`` a
237 Boolean flag indicating if the optimizer exited successfully and
238 ``message`` which describes the cause of the termination. See
239 `OptimizeResult` for a description of other attributes. If `polish`
240 was employed, and a lower minimum was obtained by the polishing, then
241 OptimizeResult also contains the ``jac`` attribute.
242 If the eventual solution does not satisfy the applied constraints
243 ``success`` will be `False`.
245 Notes
246 -----
247 Differential evolution is a stochastic population based method that is
248 useful for global optimization problems. At each pass through the
249 population the algorithm mutates each candidate solution by mixing with
250 other candidate solutions to create a trial candidate. There are several
251 strategies [3]_ for creating trial candidates, which suit some problems
252 more than others. The 'best1bin' strategy is a good starting point for
253 many systems. In this strategy two members of the population are randomly
254 chosen. Their difference is used to mutate the best member (the 'best' in
255 'best1bin'), :math:`b_0`, so far:
257 .. math::
259 b' = b_0 + mutation * (population[rand0] - population[rand1])
261 A trial vector is then constructed. Starting with a randomly chosen ith
262 parameter the trial is sequentially filled (in modulo) with parameters
263 from ``b'`` or the original candidate. The choice of whether to use ``b'``
264 or the original candidate is made with a binomial distribution (the 'bin'
265 in 'best1bin') - a random number in [0, 1) is generated. If this number is
266 less than the `recombination` constant then the parameter is loaded from
267 ``b'``, otherwise it is loaded from the original candidate. The final
268 parameter is always loaded from ``b'``. Once the trial candidate is built
269 its fitness is assessed. If the trial is better than the original candidate
270 then it takes its place. If it is also better than the best overall
271 candidate it also replaces that.
272 To improve your chances of finding a global minimum use higher `popsize`
273 values, with higher `mutation` and (dithering), but lower `recombination`
274 values. This has the effect of widening the search radius, but slowing
275 convergence.
276 By default the best solution vector is updated continuously within a single
277 iteration (``updating='immediate'``). This is a modification [4]_ of the
278 original differential evolution algorithm which can lead to faster
279 convergence as trial vectors can immediately benefit from improved
280 solutions. To use the original Storn and Price behaviour, updating the best
281 solution once per iteration, set ``updating='deferred'``.
282 The ``'deferred'`` approach is compatible with both parallelization and
283 vectorization (``'workers'`` and ``'vectorized'`` keywords). These may
284 improve minimization speed by using computer resources more efficiently.
285 The ``'workers'`` distribute calculations over multiple processors. By
286 default the Python `multiprocessing` module is used, but other approaches
287 are also possible, such as the Message Passing Interface (MPI) used on
288 clusters [6]_ [7]_. The overhead from these approaches (creating new
289 Processes, etc) may be significant, meaning that computational speed
290 doesn't necessarily scale with the number of processors used.
291 Parallelization is best suited to computationally expensive objective
292 functions. If the objective function is less expensive, then
293 ``'vectorized'`` may aid by only calling the objective function once per
294 iteration, rather than multiple times for all the population members; the
295 interpreter overhead is reduced.
297 .. versionadded:: 0.15.0
299 References
300 ----------
301 .. [1] Differential evolution, Wikipedia,
302 http://en.wikipedia.org/wiki/Differential_evolution
303 .. [2] Storn, R and Price, K, Differential Evolution - a Simple and
304 Efficient Heuristic for Global Optimization over Continuous Spaces,
305 Journal of Global Optimization, 1997, 11, 341 - 359.
306 .. [3] http://www1.icsi.berkeley.edu/~storn/code.html
307 .. [4] Wormington, M., Panaccione, C., Matney, K. M., Bowen, D. K., -
308 Characterization of structures from X-ray scattering data using
309 genetic algorithms, Phil. Trans. R. Soc. Lond. A, 1999, 357,
310 2827-2848
311 .. [5] Lampinen, J., A constraint handling approach for the differential
312 evolution algorithm. Proceedings of the 2002 Congress on
313 Evolutionary Computation. CEC'02 (Cat. No. 02TH8600). Vol. 2. IEEE,
314 2002.
315 .. [6] https://mpi4py.readthedocs.io/en/stable/
316 .. [7] https://schwimmbad.readthedocs.io/en/latest/
318 Examples
319 --------
320 Let us consider the problem of minimizing the Rosenbrock function. This
321 function is implemented in `rosen` in `scipy.optimize`.
323 >>> import numpy as np
324 >>> from scipy.optimize import rosen, differential_evolution
325 >>> bounds = [(0,2), (0, 2), (0, 2), (0, 2), (0, 2)]
326 >>> result = differential_evolution(rosen, bounds)
327 >>> result.x, result.fun
328 (array([1., 1., 1., 1., 1.]), 1.9216496320061384e-19)
330 Now repeat, but with parallelization.
332 >>> result = differential_evolution(rosen, bounds, updating='deferred',
333 ... workers=2)
334 >>> result.x, result.fun
335 (array([1., 1., 1., 1., 1.]), 1.9216496320061384e-19)
337 Let's do a constrained minimization.
339 >>> from scipy.optimize import LinearConstraint, Bounds
341 We add the constraint that the sum of ``x[0]`` and ``x[1]`` must be less
342 than or equal to 1.9. This is a linear constraint, which may be written
343 ``A @ x <= 1.9``, where ``A = array([[1, 1]])``. This can be encoded as
344 a `LinearConstraint` instance:
346 >>> lc = LinearConstraint([[1, 1]], -np.inf, 1.9)
348 Specify limits using a `Bounds` object.
350 >>> bounds = Bounds([0., 0.], [2., 2.])
351 >>> result = differential_evolution(rosen, bounds, constraints=lc,
352 ... seed=1)
353 >>> result.x, result.fun
354 (array([0.96632622, 0.93367155]), 0.0011352416852625719)
356 Next find the minimum of the Ackley function
357 (https://en.wikipedia.org/wiki/Test_functions_for_optimization).
359 >>> def ackley(x):
360 ... arg1 = -0.2 * np.sqrt(0.5 * (x[0] ** 2 + x[1] ** 2))
361 ... arg2 = 0.5 * (np.cos(2. * np.pi * x[0]) + np.cos(2. * np.pi * x[1]))
362 ... return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e
363 >>> bounds = [(-5, 5), (-5, 5)]
364 >>> result = differential_evolution(ackley, bounds, seed=1)
365 >>> result.x, result.fun
366 (array([0., 0.]), 4.440892098500626e-16)
368 The Ackley function is written in a vectorized manner, so the
369 ``'vectorized'`` keyword can be employed. Note the reduced number of
370 function evaluations.
372 >>> result = differential_evolution(
373 ... ackley, bounds, vectorized=True, updating='deferred', seed=1
374 ... )
375 >>> result.x, result.fun
376 (array([0., 0.]), 4.440892098500626e-16)
378 """
380 # using a context manager means that any created Pool objects are
381 # cleared up.
382 with DifferentialEvolutionSolver(func, bounds, args=args,
383 strategy=strategy,
384 maxiter=maxiter,
385 popsize=popsize, tol=tol,
386 mutation=mutation,
387 recombination=recombination,
388 seed=seed, polish=polish,
389 callback=callback,
390 disp=disp, init=init, atol=atol,
391 updating=updating,
392 workers=workers,
393 constraints=constraints,
394 x0=x0,
395 integrality=integrality,
396 vectorized=vectorized) as solver:
397 ret = solver.solve()
399 return ret
402class DifferentialEvolutionSolver:
404 """This class implements the differential evolution solver
406 Parameters
407 ----------
408 func : callable
409 The objective function to be minimized. Must be in the form
410 ``f(x, *args)``, where ``x`` is the argument in the form of a 1-D array
411 and ``args`` is a tuple of any additional fixed parameters needed to
412 completely specify the function. The number of parameters, N, is equal
413 to ``len(x)``.
414 bounds : sequence or `Bounds`
415 Bounds for variables. There are two ways to specify the bounds:
416 1. Instance of `Bounds` class.
417 2. ``(min, max)`` pairs for each element in ``x``, defining the finite
418 lower and upper bounds for the optimizing argument of `func`.
419 The total number of bounds is used to determine the number of
420 parameters, N.
421 args : tuple, optional
422 Any additional fixed parameters needed to
423 completely specify the objective function.
424 strategy : str, optional
425 The differential evolution strategy to use. Should be one of:
427 - 'best1bin'
428 - 'best1exp'
429 - 'rand1exp'
430 - 'randtobest1exp'
431 - 'currenttobest1exp'
432 - 'best2exp'
433 - 'rand2exp'
434 - 'randtobest1bin'
435 - 'currenttobest1bin'
436 - 'best2bin'
437 - 'rand2bin'
438 - 'rand1bin'
440 The default is 'best1bin'
442 maxiter : int, optional
443 The maximum number of generations over which the entire population is
444 evolved. The maximum number of function evaluations (with no polishing)
445 is: ``(maxiter + 1) * popsize * N``
446 popsize : int, optional
447 A multiplier for setting the total population size. The population has
448 ``popsize * N`` individuals. This keyword is overridden if an
449 initial population is supplied via the `init` keyword. When using
450 ``init='sobol'`` the population size is calculated as the next power
451 of 2 after ``popsize * N``.
452 tol : float, optional
453 Relative tolerance for convergence, the solving stops when
454 ``np.std(pop) <= atol + tol * np.abs(np.mean(population_energies))``,
455 where and `atol` and `tol` are the absolute and relative tolerance
456 respectively.
457 mutation : float or tuple(float, float), optional
458 The mutation constant. In the literature this is also known as
459 differential weight, being denoted by F.
460 If specified as a float it should be in the range [0, 2].
461 If specified as a tuple ``(min, max)`` dithering is employed. Dithering
462 randomly changes the mutation constant on a generation by generation
463 basis. The mutation constant for that generation is taken from
464 U[min, max). Dithering can help speed convergence significantly.
465 Increasing the mutation constant increases the search radius, but will
466 slow down convergence.
467 recombination : float, optional
468 The recombination constant, should be in the range [0, 1]. In the
469 literature this is also known as the crossover probability, being
470 denoted by CR. Increasing this value allows a larger number of mutants
471 to progress into the next generation, but at the risk of population
472 stability.
473 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
474 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
475 singleton is used.
476 If `seed` is an int, a new ``RandomState`` instance is used,
477 seeded with `seed`.
478 If `seed` is already a ``Generator`` or ``RandomState`` instance then
479 that instance is used.
480 Specify `seed` for repeatable minimizations.
481 disp : bool, optional
482 Prints the evaluated `func` at every iteration.
483 callback : callable, `callback(xk, convergence=val)`, optional
484 A function to follow the progress of the minimization. ``xk`` is
485 the current value of ``x0``. ``val`` represents the fractional
486 value of the population convergence. When ``val`` is greater than one
487 the function halts. If callback returns `True`, then the minimization
488 is halted (any polishing is still carried out).
489 polish : bool, optional
490 If True (default), then `scipy.optimize.minimize` with the `L-BFGS-B`
491 method is used to polish the best population member at the end, which
492 can improve the minimization slightly. If a constrained problem is
493 being studied then the `trust-constr` method is used instead. For large
494 problems with many constraints, polishing can take a long time due to
495 the Jacobian computations.
496 maxfun : int, optional
497 Set the maximum number of function evaluations. However, it probably
498 makes more sense to set `maxiter` instead.
499 init : str or array-like, optional
500 Specify which type of population initialization is performed. Should be
501 one of:
503 - 'latinhypercube'
504 - 'sobol'
505 - 'halton'
506 - 'random'
507 - array specifying the initial population. The array should have
508 shape ``(S, N)``, where S is the total population size and
509 N is the number of parameters.
510 `init` is clipped to `bounds` before use.
512 The default is 'latinhypercube'. Latin Hypercube sampling tries to
513 maximize coverage of the available parameter space.
515 'sobol' and 'halton' are superior alternatives and maximize even more
516 the parameter space. 'sobol' will enforce an initial population
517 size which is calculated as the next power of 2 after
518 ``popsize * N``. 'halton' has no requirements but is a bit less
519 efficient. See `scipy.stats.qmc` for more details.
521 'random' initializes the population randomly - this has the drawback
522 that clustering can occur, preventing the whole of parameter space
523 being covered. Use of an array to specify a population could be used,
524 for example, to create a tight bunch of initial guesses in an location
525 where the solution is known to exist, thereby reducing time for
526 convergence.
527 atol : float, optional
528 Absolute tolerance for convergence, the solving stops when
529 ``np.std(pop) <= atol + tol * np.abs(np.mean(population_energies))``,
530 where and `atol` and `tol` are the absolute and relative tolerance
531 respectively.
532 updating : {'immediate', 'deferred'}, optional
533 If ``'immediate'``, the best solution vector is continuously updated
534 within a single generation [4]_. This can lead to faster convergence as
535 trial vectors can take advantage of continuous improvements in the best
536 solution.
537 With ``'deferred'``, the best solution vector is updated once per
538 generation. Only ``'deferred'`` is compatible with parallelization or
539 vectorization, and the `workers` and `vectorized` keywords can
540 over-ride this option.
541 workers : int or map-like callable, optional
542 If `workers` is an int the population is subdivided into `workers`
543 sections and evaluated in parallel
544 (uses `multiprocessing.Pool <multiprocessing>`).
545 Supply `-1` to use all cores available to the Process.
546 Alternatively supply a map-like callable, such as
547 `multiprocessing.Pool.map` for evaluating the population in parallel.
548 This evaluation is carried out as ``workers(func, iterable)``.
549 This option will override the `updating` keyword to
550 `updating='deferred'` if `workers != 1`.
551 Requires that `func` be pickleable.
552 constraints : {NonLinearConstraint, LinearConstraint, Bounds}
553 Constraints on the solver, over and above those applied by the `bounds`
554 kwd. Uses the approach by Lampinen.
555 x0 : None or array-like, optional
556 Provides an initial guess to the minimization. Once the population has
557 been initialized this vector replaces the first (best) member. This
558 replacement is done even if `init` is given an initial population.
559 ``x0.shape == (N,)``.
560 integrality : 1-D array, optional
561 For each decision variable, a boolean value indicating whether the
562 decision variable is constrained to integer values. The array is
563 broadcast to ``(N,)``.
564 If any decision variables are constrained to be integral, they will not
565 be changed during polishing.
566 Only integer values lying between the lower and upper bounds are used.
567 If there are no integer values lying between the bounds then a
568 `ValueError` is raised.
569 vectorized : bool, optional
570 If ``vectorized is True``, `func` is sent an `x` array with
571 ``x.shape == (N, S)``, and is expected to return an array of shape
572 ``(S,)``, where `S` is the number of solution vectors to be calculated.
573 If constraints are applied, each of the functions used to construct
574 a `Constraint` object should accept an `x` array with
575 ``x.shape == (N, S)``, and return an array of shape ``(M, S)``, where
576 `M` is the number of constraint components.
577 This option is an alternative to the parallelization offered by
578 `workers`, and may help in optimization speed. This keyword is
579 ignored if ``workers != 1``.
580 This option will override the `updating` keyword to
581 ``updating='deferred'``.
582 """
584 # Dispatch of mutation strategy method (binomial or exponential).
585 _binomial = {'best1bin': '_best1',
586 'randtobest1bin': '_randtobest1',
587 'currenttobest1bin': '_currenttobest1',
588 'best2bin': '_best2',
589 'rand2bin': '_rand2',
590 'rand1bin': '_rand1'}
591 _exponential = {'best1exp': '_best1',
592 'rand1exp': '_rand1',
593 'randtobest1exp': '_randtobest1',
594 'currenttobest1exp': '_currenttobest1',
595 'best2exp': '_best2',
596 'rand2exp': '_rand2'}
598 __init_error_msg = ("The population initialization method must be one of "
599 "'latinhypercube' or 'random', or an array of shape "
600 "(S, N) where N is the number of parameters and S>5")
602 def __init__(self, func, bounds, args=(),
603 strategy='best1bin', maxiter=1000, popsize=15,
604 tol=0.01, mutation=(0.5, 1), recombination=0.7, seed=None,
605 maxfun=np.inf, callback=None, disp=False, polish=True,
606 init='latinhypercube', atol=0, updating='immediate',
607 workers=1, constraints=(), x0=None, *, integrality=None,
608 vectorized=False):
610 if strategy in self._binomial:
611 self.mutation_func = getattr(self, self._binomial[strategy])
612 elif strategy in self._exponential:
613 self.mutation_func = getattr(self, self._exponential[strategy])
614 else:
615 raise ValueError("Please select a valid mutation strategy")
616 self.strategy = strategy
618 self.callback = callback
619 self.polish = polish
621 # set the updating / parallelisation options
622 if updating in ['immediate', 'deferred']:
623 self._updating = updating
625 self.vectorized = vectorized
627 # want to use parallelisation, but updating is immediate
628 if workers != 1 and updating == 'immediate':
629 warnings.warn("differential_evolution: the 'workers' keyword has"
630 " overridden updating='immediate' to"
631 " updating='deferred'", UserWarning, stacklevel=2)
632 self._updating = 'deferred'
634 if vectorized and workers != 1:
635 warnings.warn("differential_evolution: the 'workers' keyword"
636 " overrides the 'vectorized' keyword", stacklevel=2)
637 self.vectorized = vectorized = False
639 if vectorized and updating == 'immediate':
640 warnings.warn("differential_evolution: the 'vectorized' keyword"
641 " has overridden updating='immediate' to updating"
642 "='deferred'", UserWarning, stacklevel=2)
643 self._updating = 'deferred'
645 # an object with a map method.
646 if vectorized:
647 def maplike_for_vectorized_func(func, x):
648 # send an array (N, S) to the user func,
649 # expect to receive (S,). Transposition is required because
650 # internally the population is held as (S, N)
651 return np.atleast_1d(func(x.T))
652 workers = maplike_for_vectorized_func
654 self._mapwrapper = MapWrapper(workers)
656 # relative and absolute tolerances for convergence
657 self.tol, self.atol = tol, atol
659 # Mutation constant should be in [0, 2). If specified as a sequence
660 # then dithering is performed.
661 self.scale = mutation
662 if (not np.all(np.isfinite(mutation)) or
663 np.any(np.array(mutation) >= 2) or
664 np.any(np.array(mutation) < 0)):
665 raise ValueError('The mutation constant must be a float in '
666 'U[0, 2), or specified as a tuple(min, max)'
667 ' where min < max and min, max are in U[0, 2).')
669 self.dither = None
670 if hasattr(mutation, '__iter__') and len(mutation) > 1:
671 self.dither = [mutation[0], mutation[1]]
672 self.dither.sort()
674 self.cross_over_probability = recombination
676 # we create a wrapped function to allow the use of map (and Pool.map
677 # in the future)
678 self.func = _FunctionWrapper(func, args)
679 self.args = args
681 # convert tuple of lower and upper bounds to limits
682 # [(low_0, high_0), ..., (low_n, high_n]
683 # -> [[low_0, ..., low_n], [high_0, ..., high_n]]
684 if isinstance(bounds, Bounds):
685 self.limits = np.array(new_bounds_to_old(bounds.lb,
686 bounds.ub,
687 len(bounds.lb)),
688 dtype=float).T
689 else:
690 self.limits = np.array(bounds, dtype='float').T
692 if (np.size(self.limits, 0) != 2 or not
693 np.all(np.isfinite(self.limits))):
694 raise ValueError('bounds should be a sequence containing '
695 'real valued (min, max) pairs for each value'
696 ' in x')
698 if maxiter is None: # the default used to be None
699 maxiter = 1000
700 self.maxiter = maxiter
701 if maxfun is None: # the default used to be None
702 maxfun = np.inf
703 self.maxfun = maxfun
705 # population is scaled to between [0, 1].
706 # We have to scale between parameter <-> population
707 # save these arguments for _scale_parameter and
708 # _unscale_parameter. This is an optimization
709 self.__scale_arg1 = 0.5 * (self.limits[0] + self.limits[1])
710 self.__scale_arg2 = np.fabs(self.limits[0] - self.limits[1])
712 self.parameter_count = np.size(self.limits, 1)
714 self.random_number_generator = check_random_state(seed)
716 # Which parameters are going to be integers?
717 if np.any(integrality):
718 # # user has provided a truth value for integer constraints
719 integrality = np.broadcast_to(
720 integrality,
721 self.parameter_count
722 )
723 integrality = np.asarray(integrality, bool)
724 # For integrality parameters change the limits to only allow
725 # integer values lying between the limits.
726 lb, ub = np.copy(self.limits)
728 lb = np.ceil(lb)
729 ub = np.floor(ub)
730 if not (lb[integrality] <= ub[integrality]).all():
731 # there's a parameter that doesn't have an integer value
732 # lying between the limits
733 raise ValueError("One of the integrality constraints does not"
734 " have any possible integer values between"
735 " the lower/upper bounds.")
736 nlb = np.nextafter(lb[integrality] - 0.5, np.inf)
737 nub = np.nextafter(ub[integrality] + 0.5, -np.inf)
739 self.integrality = integrality
740 self.limits[0, self.integrality] = nlb
741 self.limits[1, self.integrality] = nub
742 else:
743 self.integrality = False
745 # default population initialization is a latin hypercube design, but
746 # there are other population initializations possible.
747 # the minimum is 5 because 'best2bin' requires a population that's at
748 # least 5 long
749 self.num_population_members = max(5, popsize * self.parameter_count)
750 self.population_shape = (self.num_population_members,
751 self.parameter_count)
753 self._nfev = 0
754 # check first str otherwise will fail to compare str with array
755 if isinstance(init, str):
756 if init == 'latinhypercube':
757 self.init_population_lhs()
758 elif init == 'sobol':
759 # must be Ns = 2**m for Sobol'
760 n_s = int(2 ** np.ceil(np.log2(self.num_population_members)))
761 self.num_population_members = n_s
762 self.population_shape = (self.num_population_members,
763 self.parameter_count)
764 self.init_population_qmc(qmc_engine='sobol')
765 elif init == 'halton':
766 self.init_population_qmc(qmc_engine='halton')
767 elif init == 'random':
768 self.init_population_random()
769 else:
770 raise ValueError(self.__init_error_msg)
771 else:
772 self.init_population_array(init)
774 if x0 is not None:
775 # scale to within unit interval and
776 # ensure parameters are within bounds.
777 x0_scaled = self._unscale_parameters(np.asarray(x0))
778 if ((x0_scaled > 1.0) | (x0_scaled < 0.0)).any():
779 raise ValueError(
780 "Some entries in x0 lay outside the specified bounds"
781 )
782 self.population[0] = x0_scaled
784 # infrastructure for constraints
785 self.constraints = constraints
786 self._wrapped_constraints = []
788 if hasattr(constraints, '__len__'):
789 # sequence of constraints, this will also deal with default
790 # keyword parameter
791 for c in constraints:
792 self._wrapped_constraints.append(
793 _ConstraintWrapper(c, self.x)
794 )
795 else:
796 self._wrapped_constraints = [
797 _ConstraintWrapper(constraints, self.x)
798 ]
799 self.total_constraints = np.sum(
800 [c.num_constr for c in self._wrapped_constraints]
801 )
802 self.constraint_violation = np.zeros((self.num_population_members, 1))
803 self.feasible = np.ones(self.num_population_members, bool)
805 self.disp = disp
807 def init_population_lhs(self):
808 """
809 Initializes the population with Latin Hypercube Sampling.
810 Latin Hypercube Sampling ensures that each parameter is uniformly
811 sampled over its range.
812 """
813 rng = self.random_number_generator
815 # Each parameter range needs to be sampled uniformly. The scaled
816 # parameter range ([0, 1)) needs to be split into
817 # `self.num_population_members` segments, each of which has the following
818 # size:
819 segsize = 1.0 / self.num_population_members
821 # Within each segment we sample from a uniform random distribution.
822 # We need to do this sampling for each parameter.
823 samples = (segsize * rng.uniform(size=self.population_shape)
825 # Offset each segment to cover the entire parameter range [0, 1)
826 + np.linspace(0., 1., self.num_population_members,
827 endpoint=False)[:, np.newaxis])
829 # Create an array for population of candidate solutions.
830 self.population = np.zeros_like(samples)
832 # Initialize population of candidate solutions by permutation of the
833 # random samples.
834 for j in range(self.parameter_count):
835 order = rng.permutation(range(self.num_population_members))
836 self.population[:, j] = samples[order, j]
838 # reset population energies
839 self.population_energies = np.full(self.num_population_members,
840 np.inf)
842 # reset number of function evaluations counter
843 self._nfev = 0
845 def init_population_qmc(self, qmc_engine):
846 """Initializes the population with a QMC method.
848 QMC methods ensures that each parameter is uniformly
849 sampled over its range.
851 Parameters
852 ----------
853 qmc_engine : str
854 The QMC method to use for initialization. Can be one of
855 ``latinhypercube``, ``sobol`` or ``halton``.
857 """
858 from scipy.stats import qmc
860 rng = self.random_number_generator
862 # Create an array for population of candidate solutions.
863 if qmc_engine == 'latinhypercube':
864 sampler = qmc.LatinHypercube(d=self.parameter_count, seed=rng)
865 elif qmc_engine == 'sobol':
866 sampler = qmc.Sobol(d=self.parameter_count, seed=rng)
867 elif qmc_engine == 'halton':
868 sampler = qmc.Halton(d=self.parameter_count, seed=rng)
869 else:
870 raise ValueError(self.__init_error_msg)
872 self.population = sampler.random(n=self.num_population_members)
874 # reset population energies
875 self.population_energies = np.full(self.num_population_members,
876 np.inf)
878 # reset number of function evaluations counter
879 self._nfev = 0
881 def init_population_random(self):
882 """
883 Initializes the population at random. This type of initialization
884 can possess clustering, Latin Hypercube sampling is generally better.
885 """
886 rng = self.random_number_generator
887 self.population = rng.uniform(size=self.population_shape)
889 # reset population energies
890 self.population_energies = np.full(self.num_population_members,
891 np.inf)
893 # reset number of function evaluations counter
894 self._nfev = 0
896 def init_population_array(self, init):
897 """
898 Initializes the population with a user specified population.
900 Parameters
901 ----------
902 init : np.ndarray
903 Array specifying subset of the initial population. The array should
904 have shape (S, N), where N is the number of parameters.
905 The population is clipped to the lower and upper bounds.
906 """
907 # make sure you're using a float array
908 popn = np.asfarray(init)
910 if (np.size(popn, 0) < 5 or
911 popn.shape[1] != self.parameter_count or
912 len(popn.shape) != 2):
913 raise ValueError("The population supplied needs to have shape"
914 " (S, len(x)), where S > 4.")
916 # scale values and clip to bounds, assigning to population
917 self.population = np.clip(self._unscale_parameters(popn), 0, 1)
919 self.num_population_members = np.size(self.population, 0)
921 self.population_shape = (self.num_population_members,
922 self.parameter_count)
924 # reset population energies
925 self.population_energies = np.full(self.num_population_members,
926 np.inf)
928 # reset number of function evaluations counter
929 self._nfev = 0
931 @property
932 def x(self):
933 """
934 The best solution from the solver
935 """
936 return self._scale_parameters(self.population[0])
938 @property
939 def convergence(self):
940 """
941 The standard deviation of the population energies divided by their
942 mean.
943 """
944 if np.any(np.isinf(self.population_energies)):
945 return np.inf
946 return (np.std(self.population_energies) /
947 np.abs(np.mean(self.population_energies) + _MACHEPS))
949 def converged(self):
950 """
951 Return True if the solver has converged.
952 """
953 if np.any(np.isinf(self.population_energies)):
954 return False
956 return (np.std(self.population_energies) <=
957 self.atol +
958 self.tol * np.abs(np.mean(self.population_energies)))
960 def solve(self):
961 """
962 Runs the DifferentialEvolutionSolver.
964 Returns
965 -------
966 res : OptimizeResult
967 The optimization result represented as a ``OptimizeResult`` object.
968 Important attributes are: ``x`` the solution array, ``success`` a
969 Boolean flag indicating if the optimizer exited successfully and
970 ``message`` which describes the cause of the termination. See
971 `OptimizeResult` for a description of other attributes. If `polish`
972 was employed, and a lower minimum was obtained by the polishing,
973 then OptimizeResult also contains the ``jac`` attribute.
974 """
975 nit, warning_flag = 0, False
976 status_message = _status_message['success']
978 # The population may have just been initialized (all entries are
979 # np.inf). If it has you have to calculate the initial energies.
980 # Although this is also done in the evolve generator it's possible
981 # that someone can set maxiter=0, at which point we still want the
982 # initial energies to be calculated (the following loop isn't run).
983 if np.all(np.isinf(self.population_energies)):
984 self.feasible, self.constraint_violation = (
985 self._calculate_population_feasibilities(self.population))
987 # only work out population energies for feasible solutions
988 self.population_energies[self.feasible] = (
989 self._calculate_population_energies(
990 self.population[self.feasible]))
992 self._promote_lowest_energy()
994 # do the optimization.
995 for nit in range(1, self.maxiter + 1):
996 # evolve the population by a generation
997 try:
998 next(self)
999 except StopIteration:
1000 warning_flag = True
1001 if self._nfev > self.maxfun:
1002 status_message = _status_message['maxfev']
1003 elif self._nfev == self.maxfun:
1004 status_message = ('Maximum number of function evaluations'
1005 ' has been reached.')
1006 break
1008 if self.disp:
1009 print("differential_evolution step %d: f(x)= %g"
1010 % (nit,
1011 self.population_energies[0]))
1013 if self.callback:
1014 c = self.tol / (self.convergence + _MACHEPS)
1015 warning_flag = bool(self.callback(self.x, convergence=c))
1016 if warning_flag:
1017 status_message = ('callback function requested stop early'
1018 ' by returning True')
1020 # should the solver terminate?
1021 if warning_flag or self.converged():
1022 break
1024 else:
1025 status_message = _status_message['maxiter']
1026 warning_flag = True
1028 DE_result = OptimizeResult(
1029 x=self.x,
1030 fun=self.population_energies[0],
1031 nfev=self._nfev,
1032 nit=nit,
1033 message=status_message,
1034 success=(warning_flag is not True))
1036 if self.polish and not np.all(self.integrality):
1037 # can't polish if all the parameters are integers
1038 if np.any(self.integrality):
1039 # set the lower/upper bounds equal so that any integrality
1040 # constraints work.
1041 limits, integrality = self.limits, self.integrality
1042 limits[0, integrality] = DE_result.x[integrality]
1043 limits[1, integrality] = DE_result.x[integrality]
1045 polish_method = 'L-BFGS-B'
1047 if self._wrapped_constraints:
1048 polish_method = 'trust-constr'
1050 constr_violation = self._constraint_violation_fn(DE_result.x)
1051 if np.any(constr_violation > 0.):
1052 warnings.warn("differential evolution didn't find a"
1053 " solution satisfying the constraints,"
1054 " attempting to polish from the least"
1055 " infeasible solution", UserWarning)
1056 if self.disp:
1057 print(f"Polishing solution with '{polish_method}'")
1058 result = minimize(self.func,
1059 np.copy(DE_result.x),
1060 method=polish_method,
1061 bounds=self.limits.T,
1062 constraints=self.constraints)
1064 self._nfev += result.nfev
1065 DE_result.nfev = self._nfev
1067 # Polishing solution is only accepted if there is an improvement in
1068 # cost function, the polishing was successful and the solution lies
1069 # within the bounds.
1070 if (result.fun < DE_result.fun and
1071 result.success and
1072 np.all(result.x <= self.limits[1]) and
1073 np.all(self.limits[0] <= result.x)):
1074 DE_result.fun = result.fun
1075 DE_result.x = result.x
1076 DE_result.jac = result.jac
1077 # to keep internal state consistent
1078 self.population_energies[0] = result.fun
1079 self.population[0] = self._unscale_parameters(result.x)
1081 if self._wrapped_constraints:
1082 DE_result.constr = [c.violation(DE_result.x) for
1083 c in self._wrapped_constraints]
1084 DE_result.constr_violation = np.max(
1085 np.concatenate(DE_result.constr))
1086 DE_result.maxcv = DE_result.constr_violation
1087 if DE_result.maxcv > 0:
1088 # if the result is infeasible then success must be False
1089 DE_result.success = False
1090 DE_result.message = ("The solution does not satisfy the "
1091 f"constraints, MAXCV = {DE_result.maxcv}")
1093 return DE_result
1095 def _calculate_population_energies(self, population):
1096 """
1097 Calculate the energies of a population.
1099 Parameters
1100 ----------
1101 population : ndarray
1102 An array of parameter vectors normalised to [0, 1] using lower
1103 and upper limits. Has shape ``(np.size(population, 0), N)``.
1105 Returns
1106 -------
1107 energies : ndarray
1108 An array of energies corresponding to each population member. If
1109 maxfun will be exceeded during this call, then the number of
1110 function evaluations will be reduced and energies will be
1111 right-padded with np.inf. Has shape ``(np.size(population, 0),)``
1112 """
1113 num_members = np.size(population, 0)
1114 # S is the number of function evals left to stay under the
1115 # maxfun budget
1116 S = min(num_members, self.maxfun - self._nfev)
1118 energies = np.full(num_members, np.inf)
1120 parameters_pop = self._scale_parameters(population)
1121 try:
1122 calc_energies = list(
1123 self._mapwrapper(self.func, parameters_pop[0:S])
1124 )
1125 calc_energies = np.squeeze(calc_energies)
1126 except (TypeError, ValueError) as e:
1127 # wrong number of arguments for _mapwrapper
1128 # or wrong length returned from the mapper
1129 raise RuntimeError(
1130 "The map-like callable must be of the form f(func, iterable), "
1131 "returning a sequence of numbers the same length as 'iterable'"
1132 ) from e
1134 if calc_energies.size != S:
1135 if self.vectorized:
1136 raise RuntimeError("The vectorized function must return an"
1137 " array of shape (S,) when given an array"
1138 " of shape (len(x), S)")
1139 raise RuntimeError("func(x, *args) must return a scalar value")
1141 energies[0:S] = calc_energies
1143 if self.vectorized:
1144 self._nfev += 1
1145 else:
1146 self._nfev += S
1148 return energies
1150 def _promote_lowest_energy(self):
1151 # swaps 'best solution' into first population entry
1153 idx = np.arange(self.num_population_members)
1154 feasible_solutions = idx[self.feasible]
1155 if feasible_solutions.size:
1156 # find the best feasible solution
1157 idx_t = np.argmin(self.population_energies[feasible_solutions])
1158 l = feasible_solutions[idx_t]
1159 else:
1160 # no solution was feasible, use 'best' infeasible solution, which
1161 # will violate constraints the least
1162 l = np.argmin(np.sum(self.constraint_violation, axis=1))
1164 self.population_energies[[0, l]] = self.population_energies[[l, 0]]
1165 self.population[[0, l], :] = self.population[[l, 0], :]
1166 self.feasible[[0, l]] = self.feasible[[l, 0]]
1167 self.constraint_violation[[0, l], :] = (
1168 self.constraint_violation[[l, 0], :])
1170 def _constraint_violation_fn(self, x):
1171 """
1172 Calculates total constraint violation for all the constraints, for a
1173 set of solutions.
1175 Parameters
1176 ----------
1177 x : ndarray
1178 Solution vector(s). Has shape (S, N), or (N,), where S is the
1179 number of solutions to investigate and N is the number of
1180 parameters.
1182 Returns
1183 -------
1184 cv : ndarray
1185 Total violation of constraints. Has shape ``(S, M)``, where M is
1186 the total number of constraint components (which is not necessarily
1187 equal to len(self._wrapped_constraints)).
1188 """
1189 # how many solution vectors you're calculating constraint violations
1190 # for
1191 S = np.size(x) // self.parameter_count
1192 _out = np.zeros((S, self.total_constraints))
1193 offset = 0
1194 for con in self._wrapped_constraints:
1195 # the input/output of the (vectorized) constraint function is
1196 # {(N, S), (N,)} --> (M, S)
1197 # The input to _constraint_violation_fn is (S, N) or (N,), so
1198 # transpose to pass it to the constraint. The output is transposed
1199 # from (M, S) to (S, M) for further use.
1200 c = con.violation(x.T).T
1202 # The shape of c should be (M,), (1, M), or (S, M). Check for
1203 # those shapes, as an incorrect shape indicates that the
1204 # user constraint function didn't return the right thing, and
1205 # the reshape operation will fail. Intercept the wrong shape
1206 # to give a reasonable error message. I'm not sure what failure
1207 # modes an inventive user will come up with.
1208 if c.shape[-1] != con.num_constr or (S > 1 and c.shape[0] != S):
1209 raise RuntimeError("An array returned from a Constraint has"
1210 " the wrong shape. If `vectorized is False`"
1211 " the Constraint should return an array of"
1212 " shape (M,). If `vectorized is True` then"
1213 " the Constraint must return an array of"
1214 " shape (M, S), where S is the number of"
1215 " solution vectors and M is the number of"
1216 " constraint components in a given"
1217 " Constraint object.")
1219 # the violation function may return a 1D array, but is it a
1220 # sequence of constraints for one solution (S=1, M>=1), or the
1221 # value of a single constraint for a sequence of solutions
1222 # (S>=1, M=1)
1223 c = np.reshape(c, (S, con.num_constr))
1224 _out[:, offset:offset + con.num_constr] = c
1225 offset += con.num_constr
1227 return _out
1229 def _calculate_population_feasibilities(self, population):
1230 """
1231 Calculate the feasibilities of a population.
1233 Parameters
1234 ----------
1235 population : ndarray
1236 An array of parameter vectors normalised to [0, 1] using lower
1237 and upper limits. Has shape ``(np.size(population, 0), N)``.
1239 Returns
1240 -------
1241 feasible, constraint_violation : ndarray, ndarray
1242 Boolean array of feasibility for each population member, and an
1243 array of the constraint violation for each population member.
1244 constraint_violation has shape ``(np.size(population, 0), M)``,
1245 where M is the number of constraints.
1246 """
1247 num_members = np.size(population, 0)
1248 if not self._wrapped_constraints:
1249 # shortcut for no constraints
1250 return np.ones(num_members, bool), np.zeros((num_members, 1))
1252 # (S, N)
1253 parameters_pop = self._scale_parameters(population)
1255 if self.vectorized:
1256 # (S, M)
1257 constraint_violation = np.array(
1258 self._constraint_violation_fn(parameters_pop)
1259 )
1260 else:
1261 # (S, 1, M)
1262 constraint_violation = np.array([self._constraint_violation_fn(x)
1263 for x in parameters_pop])
1264 # if you use the list comprehension in the line above it will
1265 # create an array of shape (S, 1, M), because each iteration
1266 # generates an array of (1, M). In comparison the vectorized
1267 # version returns (S, M). It's therefore necessary to remove axis 1
1268 constraint_violation = constraint_violation[:, 0]
1270 feasible = ~(np.sum(constraint_violation, axis=1) > 0)
1272 return feasible, constraint_violation
1274 def __iter__(self):
1275 return self
1277 def __enter__(self):
1278 return self
1280 def __exit__(self, *args):
1281 return self._mapwrapper.__exit__(*args)
1283 def _accept_trial(self, energy_trial, feasible_trial, cv_trial,
1284 energy_orig, feasible_orig, cv_orig):
1285 """
1286 Trial is accepted if:
1287 * it satisfies all constraints and provides a lower or equal objective
1288 function value, while both the compared solutions are feasible
1289 - or -
1290 * it is feasible while the original solution is infeasible,
1291 - or -
1292 * it is infeasible, but provides a lower or equal constraint violation
1293 for all constraint functions.
1295 This test corresponds to section III of Lampinen [1]_.
1297 Parameters
1298 ----------
1299 energy_trial : float
1300 Energy of the trial solution
1301 feasible_trial : float
1302 Feasibility of trial solution
1303 cv_trial : array-like
1304 Excess constraint violation for the trial solution
1305 energy_orig : float
1306 Energy of the original solution
1307 feasible_orig : float
1308 Feasibility of original solution
1309 cv_orig : array-like
1310 Excess constraint violation for the original solution
1312 Returns
1313 -------
1314 accepted : bool
1316 """
1317 if feasible_orig and feasible_trial:
1318 return energy_trial <= energy_orig
1319 elif feasible_trial and not feasible_orig:
1320 return True
1321 elif not feasible_trial and (cv_trial <= cv_orig).all():
1322 # cv_trial < cv_orig would imply that both trial and orig are not
1323 # feasible
1324 return True
1326 return False
1328 def __next__(self):
1329 """
1330 Evolve the population by a single generation
1332 Returns
1333 -------
1334 x : ndarray
1335 The best solution from the solver.
1336 fun : float
1337 Value of objective function obtained from the best solution.
1338 """
1339 # the population may have just been initialized (all entries are
1340 # np.inf). If it has you have to calculate the initial energies
1341 if np.all(np.isinf(self.population_energies)):
1342 self.feasible, self.constraint_violation = (
1343 self._calculate_population_feasibilities(self.population))
1345 # only need to work out population energies for those that are
1346 # feasible
1347 self.population_energies[self.feasible] = (
1348 self._calculate_population_energies(
1349 self.population[self.feasible]))
1351 self._promote_lowest_energy()
1353 if self.dither is not None:
1354 self.scale = self.random_number_generator.uniform(self.dither[0],
1355 self.dither[1])
1357 if self._updating == 'immediate':
1358 # update best solution immediately
1359 for candidate in range(self.num_population_members):
1360 if self._nfev > self.maxfun:
1361 raise StopIteration
1363 # create a trial solution
1364 trial = self._mutate(candidate)
1366 # ensuring that it's in the range [0, 1)
1367 self._ensure_constraint(trial)
1369 # scale from [0, 1) to the actual parameter value
1370 parameters = self._scale_parameters(trial)
1372 # determine the energy of the objective function
1373 if self._wrapped_constraints:
1374 cv = self._constraint_violation_fn(parameters)
1375 feasible = False
1376 energy = np.inf
1377 if not np.sum(cv) > 0:
1378 # solution is feasible
1379 feasible = True
1380 energy = self.func(parameters)
1381 self._nfev += 1
1382 else:
1383 feasible = True
1384 cv = np.atleast_2d([0.])
1385 energy = self.func(parameters)
1386 self._nfev += 1
1388 # compare trial and population member
1389 if self._accept_trial(energy, feasible, cv,
1390 self.population_energies[candidate],
1391 self.feasible[candidate],
1392 self.constraint_violation[candidate]):
1393 self.population[candidate] = trial
1394 self.population_energies[candidate] = np.squeeze(energy)
1395 self.feasible[candidate] = feasible
1396 self.constraint_violation[candidate] = cv
1398 # if the trial candidate is also better than the best
1399 # solution then promote it.
1400 if self._accept_trial(energy, feasible, cv,
1401 self.population_energies[0],
1402 self.feasible[0],
1403 self.constraint_violation[0]):
1404 self._promote_lowest_energy()
1406 elif self._updating == 'deferred':
1407 # update best solution once per generation
1408 if self._nfev >= self.maxfun:
1409 raise StopIteration
1411 # 'deferred' approach, vectorised form.
1412 # create trial solutions
1413 trial_pop = np.array(
1414 [self._mutate(i) for i in range(self.num_population_members)])
1416 # enforce bounds
1417 self._ensure_constraint(trial_pop)
1419 # determine the energies of the objective function, but only for
1420 # feasible trials
1421 feasible, cv = self._calculate_population_feasibilities(trial_pop)
1422 trial_energies = np.full(self.num_population_members, np.inf)
1424 # only calculate for feasible entries
1425 trial_energies[feasible] = self._calculate_population_energies(
1426 trial_pop[feasible])
1428 # which solutions are 'improved'?
1429 loc = [self._accept_trial(*val) for val in
1430 zip(trial_energies, feasible, cv, self.population_energies,
1431 self.feasible, self.constraint_violation)]
1432 loc = np.array(loc)
1433 self.population = np.where(loc[:, np.newaxis],
1434 trial_pop,
1435 self.population)
1436 self.population_energies = np.where(loc,
1437 trial_energies,
1438 self.population_energies)
1439 self.feasible = np.where(loc,
1440 feasible,
1441 self.feasible)
1442 self.constraint_violation = np.where(loc[:, np.newaxis],
1443 cv,
1444 self.constraint_violation)
1446 # make sure the best solution is updated if updating='deferred'.
1447 # put the lowest energy into the best solution position.
1448 self._promote_lowest_energy()
1450 return self.x, self.population_energies[0]
1452 def _scale_parameters(self, trial):
1453 """Scale from a number between 0 and 1 to parameters."""
1454 # trial either has shape (N, ) or (L, N), where L is the number of
1455 # solutions being scaled
1456 scaled = self.__scale_arg1 + (trial - 0.5) * self.__scale_arg2
1457 if np.any(self.integrality):
1458 i = np.broadcast_to(self.integrality, scaled.shape)
1459 scaled[i] = np.round(scaled[i])
1460 return scaled
1462 def _unscale_parameters(self, parameters):
1463 """Scale from parameters to a number between 0 and 1."""
1464 return (parameters - self.__scale_arg1) / self.__scale_arg2 + 0.5
1466 def _ensure_constraint(self, trial):
1467 """Make sure the parameters lie between the limits."""
1468 mask = np.where((trial > 1) | (trial < 0))
1469 trial[mask] = self.random_number_generator.uniform(size=mask[0].shape)
1471 def _mutate(self, candidate):
1472 """Create a trial vector based on a mutation strategy."""
1473 trial = np.copy(self.population[candidate])
1475 rng = self.random_number_generator
1477 fill_point = rng.choice(self.parameter_count)
1479 if self.strategy in ['currenttobest1exp', 'currenttobest1bin']:
1480 bprime = self.mutation_func(candidate,
1481 self._select_samples(candidate, 5))
1482 else:
1483 bprime = self.mutation_func(self._select_samples(candidate, 5))
1485 if self.strategy in self._binomial:
1486 crossovers = rng.uniform(size=self.parameter_count)
1487 crossovers = crossovers < self.cross_over_probability
1488 # the last one is always from the bprime vector for binomial
1489 # If you fill in modulo with a loop you have to set the last one to
1490 # true. If you don't use a loop then you can have any random entry
1491 # be True.
1492 crossovers[fill_point] = True
1493 trial = np.where(crossovers, bprime, trial)
1494 return trial
1496 elif self.strategy in self._exponential:
1497 i = 0
1498 crossovers = rng.uniform(size=self.parameter_count)
1499 crossovers = crossovers < self.cross_over_probability
1500 crossovers[0] = True
1501 while (i < self.parameter_count and crossovers[i]):
1502 trial[fill_point] = bprime[fill_point]
1503 fill_point = (fill_point + 1) % self.parameter_count
1504 i += 1
1506 return trial
1508 def _best1(self, samples):
1509 """best1bin, best1exp"""
1510 r0, r1 = samples[:2]
1511 return (self.population[0] + self.scale *
1512 (self.population[r0] - self.population[r1]))
1514 def _rand1(self, samples):
1515 """rand1bin, rand1exp"""
1516 r0, r1, r2 = samples[:3]
1517 return (self.population[r0] + self.scale *
1518 (self.population[r1] - self.population[r2]))
1520 def _randtobest1(self, samples):
1521 """randtobest1bin, randtobest1exp"""
1522 r0, r1, r2 = samples[:3]
1523 bprime = np.copy(self.population[r0])
1524 bprime += self.scale * (self.population[0] - bprime)
1525 bprime += self.scale * (self.population[r1] -
1526 self.population[r2])
1527 return bprime
1529 def _currenttobest1(self, candidate, samples):
1530 """currenttobest1bin, currenttobest1exp"""
1531 r0, r1 = samples[:2]
1532 bprime = (self.population[candidate] + self.scale *
1533 (self.population[0] - self.population[candidate] +
1534 self.population[r0] - self.population[r1]))
1535 return bprime
1537 def _best2(self, samples):
1538 """best2bin, best2exp"""
1539 r0, r1, r2, r3 = samples[:4]
1540 bprime = (self.population[0] + self.scale *
1541 (self.population[r0] + self.population[r1] -
1542 self.population[r2] - self.population[r3]))
1544 return bprime
1546 def _rand2(self, samples):
1547 """rand2bin, rand2exp"""
1548 r0, r1, r2, r3, r4 = samples
1549 bprime = (self.population[r0] + self.scale *
1550 (self.population[r1] + self.population[r2] -
1551 self.population[r3] - self.population[r4]))
1553 return bprime
1555 def _select_samples(self, candidate, number_samples):
1556 """
1557 obtain random integers from range(self.num_population_members),
1558 without replacement. You can't have the original candidate either.
1559 """
1560 idxs = list(range(self.num_population_members))
1561 idxs.remove(candidate)
1562 self.random_number_generator.shuffle(idxs)
1563 idxs = idxs[:number_samples]
1564 return idxs
1567class _ConstraintWrapper:
1568 """Object to wrap/evaluate user defined constraints.
1570 Very similar in practice to `PreparedConstraint`, except that no evaluation
1571 of jac/hess is performed (explicit or implicit).
1573 If created successfully, it will contain the attributes listed below.
1575 Parameters
1576 ----------
1577 constraint : {`NonlinearConstraint`, `LinearConstraint`, `Bounds`}
1578 Constraint to check and prepare.
1579 x0 : array_like
1580 Initial vector of independent variables, shape (N,)
1582 Attributes
1583 ----------
1584 fun : callable
1585 Function defining the constraint wrapped by one of the convenience
1586 classes.
1587 bounds : 2-tuple
1588 Contains lower and upper bounds for the constraints --- lb and ub.
1589 These are converted to ndarray and have a size equal to the number of
1590 the constraints.
1591 """
1592 def __init__(self, constraint, x0):
1593 self.constraint = constraint
1595 if isinstance(constraint, NonlinearConstraint):
1596 def fun(x):
1597 x = np.asarray(x)
1598 return np.atleast_1d(constraint.fun(x))
1599 elif isinstance(constraint, LinearConstraint):
1600 def fun(x):
1601 if issparse(constraint.A):
1602 A = constraint.A
1603 else:
1604 A = np.atleast_2d(constraint.A)
1605 return A.dot(x)
1606 elif isinstance(constraint, Bounds):
1607 def fun(x):
1608 return np.asarray(x)
1609 else:
1610 raise ValueError("`constraint` of an unknown type is passed.")
1612 self.fun = fun
1614 lb = np.asarray(constraint.lb, dtype=float)
1615 ub = np.asarray(constraint.ub, dtype=float)
1617 x0 = np.asarray(x0)
1619 # find out the number of constraints
1620 f0 = fun(x0)
1621 self.num_constr = m = f0.size
1622 self.parameter_count = x0.size
1624 if lb.ndim == 0:
1625 lb = np.resize(lb, m)
1626 if ub.ndim == 0:
1627 ub = np.resize(ub, m)
1629 self.bounds = (lb, ub)
1631 def __call__(self, x):
1632 return np.atleast_1d(self.fun(x))
1634 def violation(self, x):
1635 """How much the constraint is exceeded by.
1637 Parameters
1638 ----------
1639 x : array-like
1640 Vector of independent variables, (N, S), where N is number of
1641 parameters and S is the number of solutions to be investigated.
1643 Returns
1644 -------
1645 excess : array-like
1646 How much the constraint is exceeded by, for each of the
1647 constraints specified by `_ConstraintWrapper.fun`.
1648 Has shape (M, S) where M is the number of constraint components.
1649 """
1650 # expect ev to have shape (num_constr, S) or (num_constr,)
1651 ev = self.fun(np.asarray(x))
1653 try:
1654 excess_lb = np.maximum(self.bounds[0] - ev.T, 0)
1655 excess_ub = np.maximum(ev.T - self.bounds[1], 0)
1656 except ValueError as e:
1657 raise RuntimeError("An array returned from a Constraint has"
1658 " the wrong shape. If `vectorized is False`"
1659 " the Constraint should return an array of"
1660 " shape (M,). If `vectorized is True` then"
1661 " the Constraint must return an array of"
1662 " shape (M, S), where S is the number of"
1663 " solution vectors and M is the number of"
1664 " constraint components in a given"
1665 " Constraint object.") from e
1667 v = (excess_lb + excess_ub).T
1668 return v