Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/optimize/_differentialevolution.py : 14%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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
12from scipy.optimize._constraints import (Bounds, new_bounds_to_old,
13 NonlinearConstraint, LinearConstraint)
14from scipy.sparse import issparse
17__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=()):
28 """Finds the global minimum of a multivariate function.
30 Differential Evolution is stochastic in nature (does not use gradient
31 methods) to find the minimum, and can search large areas of candidate
32 space, but often requires larger numbers of function evaluations than
33 conventional gradient-based techniques.
35 The algorithm is due to Storn and Price [1]_.
37 Parameters
38 ----------
39 func : callable
40 The objective function to be minimized. Must be in the form
41 ``f(x, *args)``, where ``x`` is the argument in the form of a 1-D array
42 and ``args`` is a tuple of any additional fixed parameters needed to
43 completely specify the function.
44 bounds : sequence or `Bounds`, optional
45 Bounds for variables. There are two ways to specify the bounds:
46 1. Instance of `Bounds` class.
47 2. ``(min, max)`` pairs for each element in ``x``, defining the finite
48 lower and upper bounds for the optimizing argument of `func`. It is
49 required to have ``len(bounds) == len(x)``. ``len(bounds)`` is used
50 to determine the number of parameters in ``x``.
51 args : tuple, optional
52 Any additional fixed parameters needed to
53 completely specify the objective function.
54 strategy : str, optional
55 The differential evolution strategy to use. Should be one of:
57 - 'best1bin'
58 - 'best1exp'
59 - 'rand1exp'
60 - 'randtobest1exp'
61 - 'currenttobest1exp'
62 - 'best2exp'
63 - 'rand2exp'
64 - 'randtobest1bin'
65 - 'currenttobest1bin'
66 - 'best2bin'
67 - 'rand2bin'
68 - 'rand1bin'
70 The default is 'best1bin'.
71 maxiter : int, optional
72 The maximum number of generations over which the entire population is
73 evolved. The maximum number of function evaluations (with no polishing)
74 is: ``(maxiter + 1) * popsize * len(x)``
75 popsize : int, optional
76 A multiplier for setting the total population size. The population has
77 ``popsize * len(x)`` individuals (unless the initial population is
78 supplied via the `init` keyword).
79 tol : float, optional
80 Relative tolerance for convergence, the solving stops when
81 ``np.std(pop) <= atol + tol * np.abs(np.mean(population_energies))``,
82 where and `atol` and `tol` are the absolute and relative tolerance
83 respectively.
84 mutation : float or tuple(float, float), optional
85 The mutation constant. In the literature this is also known as
86 differential weight, being denoted by F.
87 If specified as a float it should be in the range [0, 2].
88 If specified as a tuple ``(min, max)`` dithering is employed. Dithering
89 randomly changes the mutation constant on a generation by generation
90 basis. The mutation constant for that generation is taken from
91 ``U[min, max)``. Dithering can help speed convergence significantly.
92 Increasing the mutation constant increases the search radius, but will
93 slow down convergence.
94 recombination : float, optional
95 The recombination constant, should be in the range [0, 1]. In the
96 literature this is also known as the crossover probability, being
97 denoted by CR. Increasing this value allows a larger number of mutants
98 to progress into the next generation, but at the risk of population
99 stability.
100 seed : {int, `~np.random.RandomState`, `~np.random.Generator`}, optional
101 If `seed` is not specified the `~np.random.RandomState` singleton is
102 used.
103 If `seed` is an int, a new ``RandomState`` instance is used,
104 seeded with seed.
105 If `seed` is already a ``RandomState`` or a ``Generator`` instance,
106 then that object is used.
107 Specify `seed` for repeatable minimizations.
108 disp : bool, optional
109 Prints the evaluated `func` at every iteration.
110 callback : callable, `callback(xk, convergence=val)`, optional
111 A function to follow the progress of the minimization. ``xk`` is
112 the current value of ``x0``. ``val`` represents the fractional
113 value of the population convergence. When ``val`` is greater than one
114 the function halts. If callback returns `True`, then the minimization
115 is halted (any polishing is still carried out).
116 polish : bool, optional
117 If True (default), then `scipy.optimize.minimize` with the `L-BFGS-B`
118 method is used to polish the best population member at the end, which
119 can improve the minimization slightly. If a constrained problem is
120 being studied then the `trust-constr` method is used instead.
121 init : str or array-like, optional
122 Specify which type of population initialization is performed. Should be
123 one of:
125 - 'latinhypercube'
126 - 'random'
127 - array specifying the initial population. The array should have
128 shape ``(M, len(x))``, where M is the total population size and
129 len(x) is the number of parameters.
130 `init` is clipped to `bounds` before use.
132 The default is 'latinhypercube'. Latin Hypercube sampling tries to
133 maximize coverage of the available parameter space. 'random'
134 initializes the population randomly - this has the drawback that
135 clustering can occur, preventing the whole of parameter space being
136 covered. Use of an array to specify a population subset could be used,
137 for example, to create a tight bunch of initial guesses in an location
138 where the solution is known to exist, thereby reducing time for
139 convergence.
140 atol : float, optional
141 Absolute tolerance for convergence, the solving stops when
142 ``np.std(pop) <= atol + tol * np.abs(np.mean(population_energies))``,
143 where and `atol` and `tol` are the absolute and relative tolerance
144 respectively.
145 updating : {'immediate', 'deferred'}, optional
146 If ``'immediate'``, the best solution vector is continuously updated
147 within a single generation [4]_. This can lead to faster convergence as
148 trial vectors can take advantage of continuous improvements in the best
149 solution.
150 With ``'deferred'``, the best solution vector is updated once per
151 generation. Only ``'deferred'`` is compatible with parallelization, and
152 the `workers` keyword can over-ride this option.
154 .. versionadded:: 1.2.0
156 workers : int or map-like callable, optional
157 If `workers` is an int the population is subdivided into `workers`
158 sections and evaluated in parallel
159 (uses `multiprocessing.Pool <multiprocessing>`).
160 Supply -1 to use all available CPU cores.
161 Alternatively supply a map-like callable, such as
162 `multiprocessing.Pool.map` for evaluating the population in parallel.
163 This evaluation is carried out as ``workers(func, iterable)``.
164 This option will override the `updating` keyword to
165 ``updating='deferred'`` if ``workers != 1``.
166 Requires that `func` be pickleable.
168 .. versionadded:: 1.2.0
170 constraints : {NonLinearConstraint, LinearConstraint, Bounds}
171 Constraints on the solver, over and above those applied by the `bounds`
172 kwd. Uses the approach by Lampinen [5]_.
174 .. versionadded:: 1.4.0
176 Returns
177 -------
178 res : OptimizeResult
179 The optimization result represented as a `OptimizeResult` object.
180 Important attributes are: ``x`` the solution array, ``success`` a
181 Boolean flag indicating if the optimizer exited successfully and
182 ``message`` which describes the cause of the termination. See
183 `OptimizeResult` for a description of other attributes. If `polish`
184 was employed, and a lower minimum was obtained by the polishing, then
185 OptimizeResult also contains the ``jac`` attribute.
186 If the eventual solution does not satisfy the applied constraints
187 ``success`` will be `False`.
189 Notes
190 -----
191 Differential evolution is a stochastic population based method that is
192 useful for global optimization problems. At each pass through the population
193 the algorithm mutates each candidate solution by mixing with other candidate
194 solutions to create a trial candidate. There are several strategies [2]_ for
195 creating trial candidates, which suit some problems more than others. The
196 'best1bin' strategy is a good starting point for many systems. In this
197 strategy two members of the population are randomly chosen. Their difference
198 is used to mutate the best member (the 'best' in 'best1bin'), :math:`b_0`,
199 so far:
201 .. math::
203 b' = b_0 + mutation * (population[rand0] - population[rand1])
205 A trial vector is then constructed. Starting with a randomly chosen ith
206 parameter the trial is sequentially filled (in modulo) with parameters from
207 ``b'`` or the original candidate. The choice of whether to use ``b'`` or the
208 original candidate is made with a binomial distribution (the 'bin' in
209 'best1bin') - a random number in [0, 1) is generated. If this number is
210 less than the `recombination` constant then the parameter is loaded from
211 ``b'``, otherwise it is loaded from the original candidate. The final
212 parameter is always loaded from ``b'``. Once the trial candidate is built
213 its fitness is assessed. If the trial is better than the original candidate
214 then it takes its place. If it is also better than the best overall
215 candidate it also replaces that.
216 To improve your chances of finding a global minimum use higher `popsize`
217 values, with higher `mutation` and (dithering), but lower `recombination`
218 values. This has the effect of widening the search radius, but slowing
219 convergence.
220 By default the best solution vector is updated continuously within a single
221 iteration (``updating='immediate'``). This is a modification [4]_ of the
222 original differential evolution algorithm which can lead to faster
223 convergence as trial vectors can immediately benefit from improved
224 solutions. To use the original Storn and Price behaviour, updating the best
225 solution once per iteration, set ``updating='deferred'``.
227 .. versionadded:: 0.15.0
229 Examples
230 --------
231 Let us consider the problem of minimizing the Rosenbrock function. This
232 function is implemented in `rosen` in `scipy.optimize`.
234 >>> from scipy.optimize import rosen, differential_evolution
235 >>> bounds = [(0,2), (0, 2), (0, 2), (0, 2), (0, 2)]
236 >>> result = differential_evolution(rosen, bounds)
237 >>> result.x, result.fun
238 (array([1., 1., 1., 1., 1.]), 1.9216496320061384e-19)
240 Now repeat, but with parallelization.
242 >>> bounds = [(0,2), (0, 2), (0, 2), (0, 2), (0, 2)]
243 >>> result = differential_evolution(rosen, bounds, updating='deferred',
244 ... workers=2)
245 >>> result.x, result.fun
246 (array([1., 1., 1., 1., 1.]), 1.9216496320061384e-19)
248 Let's try and do a constrained minimization
250 >>> from scipy.optimize import NonlinearConstraint, Bounds
251 >>> def constr_f(x):
252 ... return np.array(x[0] + x[1])
253 >>>
254 >>> # the sum of x[0] and x[1] must be less than 1.9
255 >>> nlc = NonlinearConstraint(constr_f, -np.inf, 1.9)
256 >>> # specify limits using a `Bounds` object.
257 >>> bounds = Bounds([0., 0.], [2., 2.])
258 >>> result = differential_evolution(rosen, bounds, constraints=(nlc),
259 ... seed=1)
260 >>> result.x, result.fun
261 (array([0.96633867, 0.93363577]), 0.0011361355854792312)
263 Next find the minimum of the Ackley function
264 (https://en.wikipedia.org/wiki/Test_functions_for_optimization).
266 >>> from scipy.optimize import differential_evolution
267 >>> import numpy as np
268 >>> def ackley(x):
269 ... arg1 = -0.2 * np.sqrt(0.5 * (x[0] ** 2 + x[1] ** 2))
270 ... arg2 = 0.5 * (np.cos(2. * np.pi * x[0]) + np.cos(2. * np.pi * x[1]))
271 ... return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e
272 >>> bounds = [(-5, 5), (-5, 5)]
273 >>> result = differential_evolution(ackley, bounds)
274 >>> result.x, result.fun
275 (array([ 0., 0.]), 4.4408920985006262e-16)
277 References
278 ----------
279 .. [1] Storn, R and Price, K, Differential Evolution - a Simple and
280 Efficient Heuristic for Global Optimization over Continuous Spaces,
281 Journal of Global Optimization, 1997, 11, 341 - 359.
282 .. [2] http://www1.icsi.berkeley.edu/~storn/code.html
283 .. [3] http://en.wikipedia.org/wiki/Differential_evolution
284 .. [4] Wormington, M., Panaccione, C., Matney, K. M., Bowen, D. K., -
285 Characterization of structures from X-ray scattering data using
286 genetic algorithms, Phil. Trans. R. Soc. Lond. A, 1999, 357,
287 2827-2848
288 .. [5] Lampinen, J., A constraint handling approach for the differential
289 evolution algorithm. Proceedings of the 2002 Congress on
290 Evolutionary Computation. CEC'02 (Cat. No. 02TH8600). Vol. 2. IEEE,
291 2002.
292 """
294 # using a context manager means that any created Pool objects are
295 # cleared up.
296 with DifferentialEvolutionSolver(func, bounds, args=args,
297 strategy=strategy,
298 maxiter=maxiter,
299 popsize=popsize, tol=tol,
300 mutation=mutation,
301 recombination=recombination,
302 seed=seed, polish=polish,
303 callback=callback,
304 disp=disp, init=init, atol=atol,
305 updating=updating,
306 workers=workers,
307 constraints=constraints) as solver:
308 ret = solver.solve()
310 return ret
313class DifferentialEvolutionSolver(object):
315 """This class implements the differential evolution solver
317 Parameters
318 ----------
319 func : callable
320 The objective function to be minimized. Must be in the form
321 ``f(x, *args)``, where ``x`` is the argument in the form of a 1-D array
322 and ``args`` is a tuple of any additional fixed parameters needed to
323 completely specify the function.
324 bounds : sequence or `Bounds`, optional
325 Bounds for variables. There are two ways to specify the bounds:
326 1. Instance of `Bounds` class.
327 2. ``(min, max)`` pairs for each element in ``x``, defining the finite
328 lower and upper bounds for the optimizing argument of `func`. It is
329 required to have ``len(bounds) == len(x)``. ``len(bounds)`` is used
330 to determine the number of parameters in ``x``.
331 args : tuple, optional
332 Any additional fixed parameters needed to
333 completely specify the objective function.
334 strategy : str, optional
335 The differential evolution strategy to use. Should be one of:
337 - 'best1bin'
338 - 'best1exp'
339 - 'rand1exp'
340 - 'randtobest1exp'
341 - 'currenttobest1exp'
342 - 'best2exp'
343 - 'rand2exp'
344 - 'randtobest1bin'
345 - 'currenttobest1bin'
346 - 'best2bin'
347 - 'rand2bin'
348 - 'rand1bin'
350 The default is 'best1bin'
352 maxiter : int, optional
353 The maximum number of generations over which the entire population is
354 evolved. The maximum number of function evaluations (with no polishing)
355 is: ``(maxiter + 1) * popsize * len(x)``
356 popsize : int, optional
357 A multiplier for setting the total population size. The population has
358 ``popsize * len(x)`` individuals (unless the initial population is
359 supplied via the `init` keyword).
360 tol : float, optional
361 Relative tolerance for convergence, the solving stops when
362 ``np.std(pop) <= atol + tol * np.abs(np.mean(population_energies))``,
363 where and `atol` and `tol` are the absolute and relative tolerance
364 respectively.
365 mutation : float or tuple(float, float), optional
366 The mutation constant. In the literature this is also known as
367 differential weight, being denoted by F.
368 If specified as a float it should be in the range [0, 2].
369 If specified as a tuple ``(min, max)`` dithering is employed. Dithering
370 randomly changes the mutation constant on a generation by generation
371 basis. The mutation constant for that generation is taken from
372 U[min, max). Dithering can help speed convergence significantly.
373 Increasing the mutation constant increases the search radius, but will
374 slow down convergence.
375 recombination : float, optional
376 The recombination constant, should be in the range [0, 1]. In the
377 literature this is also known as the crossover probability, being
378 denoted by CR. Increasing this value allows a larger number of mutants
379 to progress into the next generation, but at the risk of population
380 stability.
381 seed : {int, `~np.random.RandomState`, `~np.random.Generator`}, optional
382 If `seed` is not specified the `~np.random.RandomState` singleton is
383 used.
384 If `seed` is an int, a new ``RandomState`` instance is used,
385 seeded with seed.
386 If `seed` is already a ``RandomState`` or a ``Generator`` instance,
387 then that object is used.
388 Specify `seed` for repeatable minimizations.
389 disp : bool, optional
390 Prints the evaluated `func` at every iteration.
391 callback : callable, `callback(xk, convergence=val)`, optional
392 A function to follow the progress of the minimization. ``xk`` is
393 the current value of ``x0``. ``val`` represents the fractional
394 value of the population convergence. When ``val`` is greater than one
395 the function halts. If callback returns `True`, then the minimization
396 is halted (any polishing is still carried out).
397 polish : bool, optional
398 If True (default), then `scipy.optimize.minimize` with the `L-BFGS-B`
399 method is used to polish the best population member at the end, which
400 can improve the minimization slightly. If a constrained problem is
401 being studied then the `trust-constr` method is used instead.
402 maxfun : int, optional
403 Set the maximum number of function evaluations. However, it probably
404 makes more sense to set `maxiter` instead.
405 init : str or array-like, optional
406 Specify which type of population initialization is performed. Should be
407 one of:
409 - 'latinhypercube'
410 - 'random'
411 - array specifying the initial population. The array should have
412 shape ``(M, len(x))``, where M is the total population size and
413 len(x) is the number of parameters.
414 `init` is clipped to `bounds` before use.
416 The default is 'latinhypercube'. Latin Hypercube sampling tries to
417 maximize coverage of the available parameter space. 'random'
418 initializes the population randomly - this has the drawback that
419 clustering can occur, preventing the whole of parameter space being
420 covered. Use of an array to specify a population could be used, for
421 example, to create a tight bunch of initial guesses in an location
422 where the solution is known to exist, thereby reducing time for
423 convergence.
424 atol : float, optional
425 Absolute tolerance for convergence, the solving stops when
426 ``np.std(pop) <= atol + tol * np.abs(np.mean(population_energies))``,
427 where and `atol` and `tol` are the absolute and relative tolerance
428 respectively.
429 updating : {'immediate', 'deferred'}, optional
430 If `immediate` the best solution vector is continuously updated within
431 a single generation. This can lead to faster convergence as trial
432 vectors can take advantage of continuous improvements in the best
433 solution.
434 With `deferred` the best solution vector is updated once per
435 generation. Only `deferred` is compatible with parallelization, and the
436 `workers` keyword can over-ride this option.
437 workers : int or map-like callable, optional
438 If `workers` is an int the population is subdivided into `workers`
439 sections and evaluated in parallel
440 (uses `multiprocessing.Pool <multiprocessing>`).
441 Supply `-1` to use all cores available to the Process.
442 Alternatively supply a map-like callable, such as
443 `multiprocessing.Pool.map` for evaluating the population in parallel.
444 This evaluation is carried out as ``workers(func, iterable)``.
445 This option will override the `updating` keyword to
446 `updating='deferred'` if `workers != 1`.
447 Requires that `func` be pickleable.
448 constraints : {NonLinearConstraint, LinearConstraint, Bounds}
449 Constraints on the solver, over and above those applied by the `bounds`
450 kwd. Uses the approach by Lampinen.
451 """
453 # Dispatch of mutation strategy method (binomial or exponential).
454 _binomial = {'best1bin': '_best1',
455 'randtobest1bin': '_randtobest1',
456 'currenttobest1bin': '_currenttobest1',
457 'best2bin': '_best2',
458 'rand2bin': '_rand2',
459 'rand1bin': '_rand1'}
460 _exponential = {'best1exp': '_best1',
461 'rand1exp': '_rand1',
462 'randtobest1exp': '_randtobest1',
463 'currenttobest1exp': '_currenttobest1',
464 'best2exp': '_best2',
465 'rand2exp': '_rand2'}
467 __init_error_msg = ("The population initialization method must be one of "
468 "'latinhypercube' or 'random', or an array of shape "
469 "(M, N) where N is the number of parameters and M>5")
471 def __init__(self, func, bounds, args=(),
472 strategy='best1bin', maxiter=1000, popsize=15,
473 tol=0.01, mutation=(0.5, 1), recombination=0.7, seed=None,
474 maxfun=np.inf, callback=None, disp=False, polish=True,
475 init='latinhypercube', atol=0, updating='immediate',
476 workers=1, constraints=()):
478 if strategy in self._binomial:
479 self.mutation_func = getattr(self, self._binomial[strategy])
480 elif strategy in self._exponential:
481 self.mutation_func = getattr(self, self._exponential[strategy])
482 else:
483 raise ValueError("Please select a valid mutation strategy")
484 self.strategy = strategy
486 self.callback = callback
487 self.polish = polish
489 # set the updating / parallelisation options
490 if updating in ['immediate', 'deferred']:
491 self._updating = updating
493 # want to use parallelisation, but updating is immediate
494 if workers != 1 and updating == 'immediate':
495 warnings.warn("differential_evolution: the 'workers' keyword has"
496 " overridden updating='immediate' to"
497 " updating='deferred'", UserWarning)
498 self._updating = 'deferred'
500 # an object with a map method.
501 self._mapwrapper = MapWrapper(workers)
503 # relative and absolute tolerances for convergence
504 self.tol, self.atol = tol, atol
506 # Mutation constant should be in [0, 2). If specified as a sequence
507 # then dithering is performed.
508 self.scale = mutation
509 if (not np.all(np.isfinite(mutation)) or
510 np.any(np.array(mutation) >= 2) or
511 np.any(np.array(mutation) < 0)):
512 raise ValueError('The mutation constant must be a float in '
513 'U[0, 2), or specified as a tuple(min, max)'
514 ' where min < max and min, max are in U[0, 2).')
516 self.dither = None
517 if hasattr(mutation, '__iter__') and len(mutation) > 1:
518 self.dither = [mutation[0], mutation[1]]
519 self.dither.sort()
521 self.cross_over_probability = recombination
523 # we create a wrapped function to allow the use of map (and Pool.map
524 # in the future)
525 self.func = _FunctionWrapper(func, args)
526 self.args = args
528 # convert tuple of lower and upper bounds to limits
529 # [(low_0, high_0), ..., (low_n, high_n]
530 # -> [[low_0, ..., low_n], [high_0, ..., high_n]]
531 if isinstance(bounds, Bounds):
532 self.limits = np.array(new_bounds_to_old(bounds.lb,
533 bounds.ub,
534 len(bounds.lb)),
535 dtype=float).T
536 else:
537 self.limits = np.array(bounds, dtype='float').T
539 if (np.size(self.limits, 0) != 2 or not
540 np.all(np.isfinite(self.limits))):
541 raise ValueError('bounds should be a sequence containing '
542 'real valued (min, max) pairs for each value'
543 ' in x')
545 if maxiter is None: # the default used to be None
546 maxiter = 1000
547 self.maxiter = maxiter
548 if maxfun is None: # the default used to be None
549 maxfun = np.inf
550 self.maxfun = maxfun
552 # population is scaled to between [0, 1].
553 # We have to scale between parameter <-> population
554 # save these arguments for _scale_parameter and
555 # _unscale_parameter. This is an optimization
556 self.__scale_arg1 = 0.5 * (self.limits[0] + self.limits[1])
557 self.__scale_arg2 = np.fabs(self.limits[0] - self.limits[1])
559 self.parameter_count = np.size(self.limits, 1)
561 self.random_number_generator = check_random_state(seed)
563 # default population initialization is a latin hypercube design, but
564 # there are other population initializations possible.
565 # the minimum is 5 because 'best2bin' requires a population that's at
566 # least 5 long
567 self.num_population_members = max(5, popsize * self.parameter_count)
569 self.population_shape = (self.num_population_members,
570 self.parameter_count)
572 self._nfev = 0
573 if isinstance(init, str):
574 if init == 'latinhypercube':
575 self.init_population_lhs()
576 elif init == 'random':
577 self.init_population_random()
578 else:
579 raise ValueError(self.__init_error_msg)
580 else:
581 self.init_population_array(init)
583 # infrastructure for constraints
584 # dummy parameter vector for preparing constraints, this is required so
585 # that the number of constraints is known.
586 x0 = self._scale_parameters(self.population[0])
588 self.constraints = constraints
589 self._wrapped_constraints = []
591 if hasattr(constraints, '__len__'):
592 # sequence of constraints, this will also deal with default
593 # keyword parameter
594 for c in constraints:
595 self._wrapped_constraints.append(_ConstraintWrapper(c, x0))
596 else:
597 self._wrapped_constraints = [_ConstraintWrapper(constraints, x0)]
599 self.constraint_violation = np.zeros((self.num_population_members, 1))
600 self.feasible = np.ones(self.num_population_members, bool)
602 self.disp = disp
604 def init_population_lhs(self):
605 """
606 Initializes the population with Latin Hypercube Sampling.
607 Latin Hypercube Sampling ensures that each parameter is uniformly
608 sampled over its range.
609 """
610 rng = self.random_number_generator
612 # Each parameter range needs to be sampled uniformly. The scaled
613 # parameter range ([0, 1)) needs to be split into
614 # `self.num_population_members` segments, each of which has the following
615 # size:
616 segsize = 1.0 / self.num_population_members
618 # Within each segment we sample from a uniform random distribution.
619 # We need to do this sampling for each parameter.
620 samples = (segsize * rng.uniform(size=self.population_shape)
622 # Offset each segment to cover the entire parameter range [0, 1)
623 + np.linspace(0., 1., self.num_population_members,
624 endpoint=False)[:, np.newaxis])
626 # Create an array for population of candidate solutions.
627 self.population = np.zeros_like(samples)
629 # Initialize population of candidate solutions by permutation of the
630 # random samples.
631 for j in range(self.parameter_count):
632 order = rng.permutation(range(self.num_population_members))
633 self.population[:, j] = samples[order, j]
635 # reset population energies
636 self.population_energies = np.full(self.num_population_members,
637 np.inf)
639 # reset number of function evaluations counter
640 self._nfev = 0
642 def init_population_random(self):
643 """
644 Initializes the population at random. This type of initialization
645 can possess clustering, Latin Hypercube sampling is generally better.
646 """
647 rng = self.random_number_generator
648 self.population = rng.uniform(size=self.population_shape)
650 # reset population energies
651 self.population_energies = np.full(self.num_population_members,
652 np.inf)
654 # reset number of function evaluations counter
655 self._nfev = 0
657 def init_population_array(self, init):
658 """
659 Initializes the population with a user specified population.
661 Parameters
662 ----------
663 init : np.ndarray
664 Array specifying subset of the initial population. The array should
665 have shape (M, len(x)), where len(x) is the number of parameters.
666 The population is clipped to the lower and upper bounds.
667 """
668 # make sure you're using a float array
669 popn = np.asfarray(init)
671 if (np.size(popn, 0) < 5 or
672 popn.shape[1] != self.parameter_count or
673 len(popn.shape) != 2):
674 raise ValueError("The population supplied needs to have shape"
675 " (M, len(x)), where M > 4.")
677 # scale values and clip to bounds, assigning to population
678 self.population = np.clip(self._unscale_parameters(popn), 0, 1)
680 self.num_population_members = np.size(self.population, 0)
682 self.population_shape = (self.num_population_members,
683 self.parameter_count)
685 # reset population energies
686 self.population_energies = np.full(self.num_population_members,
687 np.inf)
689 # reset number of function evaluations counter
690 self._nfev = 0
692 @property
693 def x(self):
694 """
695 The best solution from the solver
696 """
697 return self._scale_parameters(self.population[0])
699 @property
700 def convergence(self):
701 """
702 The standard deviation of the population energies divided by their
703 mean.
704 """
705 if np.any(np.isinf(self.population_energies)):
706 return np.inf
707 return (np.std(self.population_energies) /
708 np.abs(np.mean(self.population_energies) + _MACHEPS))
710 def converged(self):
711 """
712 Return True if the solver has converged.
713 """
714 if np.any(np.isinf(self.population_energies)):
715 return False
717 return (np.std(self.population_energies) <=
718 self.atol +
719 self.tol * np.abs(np.mean(self.population_energies)))
721 def solve(self):
722 """
723 Runs the DifferentialEvolutionSolver.
725 Returns
726 -------
727 res : OptimizeResult
728 The optimization result represented as a ``OptimizeResult`` object.
729 Important attributes are: ``x`` the solution array, ``success`` a
730 Boolean flag indicating if the optimizer exited successfully and
731 ``message`` which describes the cause of the termination. See
732 `OptimizeResult` for a description of other attributes. If `polish`
733 was employed, and a lower minimum was obtained by the polishing,
734 then OptimizeResult also contains the ``jac`` attribute.
735 """
736 nit, warning_flag = 0, False
737 status_message = _status_message['success']
739 # The population may have just been initialized (all entries are
740 # np.inf). If it has you have to calculate the initial energies.
741 # Although this is also done in the evolve generator it's possible
742 # that someone can set maxiter=0, at which point we still want the
743 # initial energies to be calculated (the following loop isn't run).
744 if np.all(np.isinf(self.population_energies)):
745 self.feasible, self.constraint_violation = (
746 self._calculate_population_feasibilities(self.population))
748 # only work out population energies for feasible solutions
749 self.population_energies[self.feasible] = (
750 self._calculate_population_energies(
751 self.population[self.feasible]))
753 self._promote_lowest_energy()
755 # do the optimization.
756 for nit in range(1, self.maxiter + 1):
757 # evolve the population by a generation
758 try:
759 next(self)
760 except StopIteration:
761 warning_flag = True
762 if self._nfev > self.maxfun:
763 status_message = _status_message['maxfev']
764 elif self._nfev == self.maxfun:
765 status_message = ('Maximum number of function evaluations'
766 ' has been reached.')
767 break
769 if self.disp:
770 print("differential_evolution step %d: f(x)= %g"
771 % (nit,
772 self.population_energies[0]))
774 if self.callback:
775 c = self.tol / (self.convergence + _MACHEPS)
776 warning_flag = bool(self.callback(self.x, convergence=c))
777 if warning_flag:
778 status_message = ('callback function requested stop early'
779 ' by returning True')
781 # should the solver terminate?
782 if warning_flag or self.converged():
783 break
785 else:
786 status_message = _status_message['maxiter']
787 warning_flag = True
789 DE_result = OptimizeResult(
790 x=self.x,
791 fun=self.population_energies[0],
792 nfev=self._nfev,
793 nit=nit,
794 message=status_message,
795 success=(warning_flag is not True))
797 if self.polish:
798 polish_method = 'L-BFGS-B'
800 if self._wrapped_constraints:
801 polish_method = 'trust-constr'
803 constr_violation = self._constraint_violation_fn(DE_result.x)
804 if np.any(constr_violation > 0.):
805 warnings.warn("differential evolution didn't find a"
806 " solution satisfying the constraints,"
807 " attempting to polish from the least"
808 " infeasible solution", UserWarning)
810 result = minimize(self.func,
811 np.copy(DE_result.x),
812 method=polish_method,
813 bounds=self.limits.T,
814 constraints=self.constraints)
816 self._nfev += result.nfev
817 DE_result.nfev = self._nfev
819 # Polishing solution is only accepted if there is an improvement in
820 # cost function, the polishing was successful and the solution lies
821 # within the bounds.
822 if (result.fun < DE_result.fun and
823 result.success and
824 np.all(result.x <= self.limits[1]) and
825 np.all(self.limits[0] <= result.x)):
826 DE_result.fun = result.fun
827 DE_result.x = result.x
828 DE_result.jac = result.jac
829 # to keep internal state consistent
830 self.population_energies[0] = result.fun
831 self.population[0] = self._unscale_parameters(result.x)
833 if self._wrapped_constraints:
834 DE_result.constr = [c.violation(DE_result.x) for
835 c in self._wrapped_constraints]
836 DE_result.constr_violation = np.max(
837 np.concatenate(DE_result.constr))
838 DE_result.maxcv = DE_result.constr_violation
839 if DE_result.maxcv > 0:
840 # if the result is infeasible then success must be False
841 DE_result.success = False
842 DE_result.message = ("The solution does not satisfy the"
843 " constraints, MAXCV = " % DE_result.maxcv)
845 return DE_result
847 def _calculate_population_energies(self, population):
848 """
849 Calculate the energies of a population.
851 Parameters
852 ----------
853 population : ndarray
854 An array of parameter vectors normalised to [0, 1] using lower
855 and upper limits. Has shape ``(np.size(population, 0), len(x))``.
857 Returns
858 -------
859 energies : ndarray
860 An array of energies corresponding to each population member. If
861 maxfun will be exceeded during this call, then the number of
862 function evaluations will be reduced and energies will be
863 right-padded with np.inf. Has shape ``(np.size(population, 0),)``
864 """
865 num_members = np.size(population, 0)
866 nfevs = min(num_members,
867 self.maxfun - num_members)
869 energies = np.full(num_members, np.inf)
871 parameters_pop = self._scale_parameters(population)
872 try:
873 calc_energies = list(self._mapwrapper(self.func,
874 parameters_pop[0:nfevs]))
875 energies[0:nfevs] = calc_energies
876 except (TypeError, ValueError):
877 # wrong number of arguments for _mapwrapper
878 # or wrong length returned from the mapper
879 raise RuntimeError("The map-like callable must be of the"
880 " form f(func, iterable), returning a sequence"
881 " of numbers the same length as 'iterable'")
883 self._nfev += nfevs
885 return energies
887 def _promote_lowest_energy(self):
888 # swaps 'best solution' into first population entry
890 idx = np.arange(self.num_population_members)
891 feasible_solutions = idx[self.feasible]
892 if feasible_solutions.size:
893 # find the best feasible solution
894 idx_t = np.argmin(self.population_energies[feasible_solutions])
895 l = feasible_solutions[idx_t]
896 else:
897 # no solution was feasible, use 'best' infeasible solution, which
898 # will violate constraints the least
899 l = np.argmin(np.sum(self.constraint_violation, axis=1))
901 self.population_energies[[0, l]] = self.population_energies[[l, 0]]
902 self.population[[0, l], :] = self.population[[l, 0], :]
903 self.feasible[[0, l]] = self.feasible[[l, 0]]
904 self.constraint_violation[[0, l], :] = (
905 self.constraint_violation[[l, 0], :])
907 def _constraint_violation_fn(self, x):
908 """
909 Calculates total constraint violation for all the constraints, for a given
910 solution.
912 Parameters
913 ----------
914 x : ndarray
915 Solution vector
917 Returns
918 -------
919 cv : ndarray
920 Total violation of constraints. Has shape ``(M,)``, where M is the
921 number of constraints (if each constraint function only returns one
922 value)
923 """
924 return np.concatenate([c.violation(x) for c in self._wrapped_constraints])
926 def _calculate_population_feasibilities(self, population):
927 """
928 Calculate the feasibilities of a population.
930 Parameters
931 ----------
932 population : ndarray
933 An array of parameter vectors normalised to [0, 1] using lower
934 and upper limits. Has shape ``(np.size(population, 0), len(x))``.
936 Returns
937 -------
938 feasible, constraint_violation : ndarray, ndarray
939 Boolean array of feasibility for each population member, and an
940 array of the constraint violation for each population member.
941 constraint_violation has shape ``(np.size(population, 0), M)``,
942 where M is the number of constraints.
943 """
944 num_members = np.size(population, 0)
945 if not self._wrapped_constraints:
946 # shortcut for no constraints
947 return np.ones(num_members, bool), np.zeros((num_members, 1))
949 parameters_pop = self._scale_parameters(population)
951 constraint_violation = np.array([self._constraint_violation_fn(x)
952 for x in parameters_pop])
953 feasible = ~(np.sum(constraint_violation, axis=1) > 0)
955 return feasible, constraint_violation
957 def __iter__(self):
958 return self
960 def __enter__(self):
961 return self
963 def __exit__(self, *args):
964 # to make sure resources are closed down
965 self._mapwrapper.close()
966 self._mapwrapper.terminate()
968 def __del__(self):
969 # to make sure resources are closed down
970 self._mapwrapper.close()
971 self._mapwrapper.terminate()
973 def _accept_trial(self, energy_trial, feasible_trial, cv_trial,
974 energy_orig, feasible_orig, cv_orig):
975 """
976 Trial is accepted if:
977 * it satisfies all constraints and provides a lower or equal objective
978 function value, while both the compared solutions are feasible
979 - or -
980 * it is feasible while the original solution is infeasible,
981 - or -
982 * it is infeasible, but provides a lower or equal constraint violation
983 for all constraint functions.
985 This test corresponds to section III of Lampinen [1]_.
987 Parameters
988 ----------
989 energy_trial : float
990 Energy of the trial solution
991 feasible_trial : float
992 Feasibility of trial solution
993 cv_trial : array-like
994 Excess constraint violation for the trial solution
995 energy_orig : float
996 Energy of the original solution
997 feasible_orig : float
998 Feasibility of original solution
999 cv_orig : array-like
1000 Excess constraint violation for the original solution
1002 Returns
1003 -------
1004 accepted : bool
1006 """
1007 if feasible_orig and feasible_trial:
1008 return energy_trial <= energy_orig
1009 elif feasible_trial and not feasible_orig:
1010 return True
1011 elif not feasible_trial and (cv_trial <= cv_orig).all():
1012 # cv_trial < cv_orig would imply that both trial and orig are not
1013 # feasible
1014 return True
1016 return False
1018 def __next__(self):
1019 """
1020 Evolve the population by a single generation
1022 Returns
1023 -------
1024 x : ndarray
1025 The best solution from the solver.
1026 fun : float
1027 Value of objective function obtained from the best solution.
1028 """
1029 # the population may have just been initialized (all entries are
1030 # np.inf). If it has you have to calculate the initial energies
1031 if np.all(np.isinf(self.population_energies)):
1032 self.feasible, self.constraint_violation = (
1033 self._calculate_population_feasibilities(self.population))
1035 # only need to work out population energies for those that are
1036 # feasible
1037 self.population_energies[self.feasible] = (
1038 self._calculate_population_energies(
1039 self.population[self.feasible]))
1041 self._promote_lowest_energy()
1043 if self.dither is not None:
1044 self.scale = self.random_number_generator.uniform(self.dither[0],
1045 self.dither[1])
1047 if self._updating == 'immediate':
1048 # update best solution immediately
1049 for candidate in range(self.num_population_members):
1050 if self._nfev > self.maxfun:
1051 raise StopIteration
1053 # create a trial solution
1054 trial = self._mutate(candidate)
1056 # ensuring that it's in the range [0, 1)
1057 self._ensure_constraint(trial)
1059 # scale from [0, 1) to the actual parameter value
1060 parameters = self._scale_parameters(trial)
1062 # determine the energy of the objective function
1063 if self._wrapped_constraints:
1064 cv = self._constraint_violation_fn(parameters)
1065 feasible = False
1066 energy = np.inf
1067 if not np.sum(cv) > 0:
1068 # solution is feasible
1069 feasible = True
1070 energy = self.func(parameters)
1071 self._nfev += 1
1072 else:
1073 feasible = True
1074 cv = np.atleast_2d([0.])
1075 energy = self.func(parameters)
1076 self._nfev += 1
1078 # compare trial and population member
1079 if self._accept_trial(energy, feasible, cv,
1080 self.population_energies[candidate],
1081 self.feasible[candidate],
1082 self.constraint_violation[candidate]):
1083 self.population[candidate] = trial
1084 self.population_energies[candidate] = energy
1085 self.feasible[candidate] = feasible
1086 self.constraint_violation[candidate] = cv
1088 # if the trial candidate is also better than the best
1089 # solution then promote it.
1090 if self._accept_trial(energy, feasible, cv,
1091 self.population_energies[0],
1092 self.feasible[0],
1093 self.constraint_violation[0]):
1094 self._promote_lowest_energy()
1096 elif self._updating == 'deferred':
1097 # update best solution once per generation
1098 if self._nfev >= self.maxfun:
1099 raise StopIteration
1101 # 'deferred' approach, vectorised form.
1102 # create trial solutions
1103 trial_pop = np.array(
1104 [self._mutate(i) for i in range(self.num_population_members)])
1106 # enforce bounds
1107 self._ensure_constraint(trial_pop)
1109 # determine the energies of the objective function, but only for
1110 # feasible trials
1111 feasible, cv = self._calculate_population_feasibilities(trial_pop)
1112 trial_energies = np.full(self.num_population_members, np.inf)
1114 # only calculate for feasible entries
1115 trial_energies[feasible] = self._calculate_population_energies(
1116 trial_pop[feasible])
1118 # which solutions are 'improved'?
1119 loc = [self._accept_trial(*val) for val in
1120 zip(trial_energies, feasible, cv, self.population_energies,
1121 self.feasible, self.constraint_violation)]
1122 loc = np.array(loc)
1123 self.population = np.where(loc[:, np.newaxis],
1124 trial_pop,
1125 self.population)
1126 self.population_energies = np.where(loc,
1127 trial_energies,
1128 self.population_energies)
1129 self.feasible = np.where(loc,
1130 feasible,
1131 self.feasible)
1132 self.constraint_violation = np.where(loc[:, np.newaxis],
1133 cv,
1134 self.constraint_violation)
1136 # make sure the best solution is updated if updating='deferred'.
1137 # put the lowest energy into the best solution position.
1138 self._promote_lowest_energy()
1140 return self.x, self.population_energies[0]
1142 next = __next__
1144 def _scale_parameters(self, trial):
1145 """Scale from a number between 0 and 1 to parameters."""
1146 return self.__scale_arg1 + (trial - 0.5) * self.__scale_arg2
1148 def _unscale_parameters(self, parameters):
1149 """Scale from parameters to a number between 0 and 1."""
1150 return (parameters - self.__scale_arg1) / self.__scale_arg2 + 0.5
1152 def _ensure_constraint(self, trial):
1153 """Make sure the parameters lie between the limits."""
1154 mask = np.where((trial > 1) | (trial < 0))
1155 trial[mask] = self.random_number_generator.uniform(size=mask[0].shape)
1157 def _mutate(self, candidate):
1158 """Create a trial vector based on a mutation strategy."""
1159 trial = np.copy(self.population[candidate])
1161 rng = self.random_number_generator
1163 fill_point = rng.choice(self.parameter_count)
1165 if self.strategy in ['currenttobest1exp', 'currenttobest1bin']:
1166 bprime = self.mutation_func(candidate,
1167 self._select_samples(candidate, 5))
1168 else:
1169 bprime = self.mutation_func(self._select_samples(candidate, 5))
1171 if self.strategy in self._binomial:
1172 crossovers = rng.uniform(size=self.parameter_count)
1173 crossovers = crossovers < self.cross_over_probability
1174 # the last one is always from the bprime vector for binomial
1175 # If you fill in modulo with a loop you have to set the last one to
1176 # true. If you don't use a loop then you can have any random entry
1177 # be True.
1178 crossovers[fill_point] = True
1179 trial = np.where(crossovers, bprime, trial)
1180 return trial
1182 elif self.strategy in self._exponential:
1183 i = 0
1184 crossovers = rng.uniform(size=self.parameter_count)
1185 crossovers = crossovers < self.cross_over_probability
1186 while (i < self.parameter_count and crossovers[i]):
1187 trial[fill_point] = bprime[fill_point]
1188 fill_point = (fill_point + 1) % self.parameter_count
1189 i += 1
1191 return trial
1193 def _best1(self, samples):
1194 """best1bin, best1exp"""
1195 r0, r1 = samples[:2]
1196 return (self.population[0] + self.scale *
1197 (self.population[r0] - self.population[r1]))
1199 def _rand1(self, samples):
1200 """rand1bin, rand1exp"""
1201 r0, r1, r2 = samples[:3]
1202 return (self.population[r0] + self.scale *
1203 (self.population[r1] - self.population[r2]))
1205 def _randtobest1(self, samples):
1206 """randtobest1bin, randtobest1exp"""
1207 r0, r1, r2 = samples[:3]
1208 bprime = np.copy(self.population[r0])
1209 bprime += self.scale * (self.population[0] - bprime)
1210 bprime += self.scale * (self.population[r1] -
1211 self.population[r2])
1212 return bprime
1214 def _currenttobest1(self, candidate, samples):
1215 """currenttobest1bin, currenttobest1exp"""
1216 r0, r1 = samples[:2]
1217 bprime = (self.population[candidate] + self.scale *
1218 (self.population[0] - self.population[candidate] +
1219 self.population[r0] - self.population[r1]))
1220 return bprime
1222 def _best2(self, samples):
1223 """best2bin, best2exp"""
1224 r0, r1, r2, r3 = samples[:4]
1225 bprime = (self.population[0] + self.scale *
1226 (self.population[r0] + self.population[r1] -
1227 self.population[r2] - self.population[r3]))
1229 return bprime
1231 def _rand2(self, samples):
1232 """rand2bin, rand2exp"""
1233 r0, r1, r2, r3, r4 = samples
1234 bprime = (self.population[r0] + self.scale *
1235 (self.population[r1] + self.population[r2] -
1236 self.population[r3] - self.population[r4]))
1238 return bprime
1240 def _select_samples(self, candidate, number_samples):
1241 """
1242 obtain random integers from range(self.num_population_members),
1243 without replacement. You can't have the original candidate either.
1244 """
1245 idxs = list(range(self.num_population_members))
1246 idxs.remove(candidate)
1247 self.random_number_generator.shuffle(idxs)
1248 idxs = idxs[:number_samples]
1249 return idxs
1252class _FunctionWrapper(object):
1253 """
1254 Object to wrap user cost function, allowing picklability
1255 """
1256 def __init__(self, f, args):
1257 self.f = f
1258 self.args = [] if args is None else args
1260 def __call__(self, x):
1261 return self.f(x, *self.args)
1264class _ConstraintWrapper(object):
1265 """Object to wrap/evaluate user defined constraints.
1267 Very similar in practice to `PreparedConstraint`, except that no evaluation
1268 of jac/hess is performed (explicit or implicit).
1270 If created successfully, it will contain the attributes listed below.
1272 Parameters
1273 ----------
1274 constraint : {`NonlinearConstraint`, `LinearConstraint`, `Bounds`}
1275 Constraint to check and prepare.
1276 x0 : array_like
1277 Initial vector of independent variables.
1279 Attributes
1280 ----------
1281 fun : callable
1282 Function defining the constraint wrapped by one of the convenience
1283 classes.
1284 bounds : 2-tuple
1285 Contains lower and upper bounds for the constraints --- lb and ub.
1286 These are converted to ndarray and have a size equal to the number of
1287 the constraints.
1288 """
1289 def __init__(self, constraint, x0):
1290 self.constraint = constraint
1292 if isinstance(constraint, NonlinearConstraint):
1293 def fun(x):
1294 return np.atleast_1d(constraint.fun(x))
1295 elif isinstance(constraint, LinearConstraint):
1296 def fun(x):
1297 if issparse(constraint.A):
1298 A = constraint.A
1299 else:
1300 A = np.atleast_2d(constraint.A)
1301 return A.dot(x)
1302 elif isinstance(constraint, Bounds):
1303 def fun(x):
1304 return x
1305 else:
1306 raise ValueError("`constraint` of an unknown type is passed.")
1308 self.fun = fun
1310 lb = np.asarray(constraint.lb, dtype=float)
1311 ub = np.asarray(constraint.ub, dtype=float)
1313 f0 = fun(x0)
1314 m = f0.size
1316 if lb.ndim == 0:
1317 lb = np.resize(lb, m)
1318 if ub.ndim == 0:
1319 ub = np.resize(ub, m)
1321 self.bounds = (lb, ub)
1323 def __call__(self, x):
1324 return np.atleast_1d(self.fun(x))
1326 def violation(self, x):
1327 """How much the constraint is exceeded by.
1329 Parameters
1330 ----------
1331 x : array-like
1332 Vector of independent variables
1334 Returns
1335 -------
1336 excess : array-like
1337 How much the constraint is exceeded by, for each of the
1338 constraints specified by `_ConstraintWrapper.fun`.
1339 """
1340 ev = self.fun(np.asarray(x))
1342 excess_lb = np.maximum(self.bounds[0] - ev, 0)
1343 excess_ub = np.maximum(ev - self.bounds[1], 0)
1345 return excess_lb + excess_ub