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