Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_dual_annealing.py: 12%
291 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1# Dual Annealing implementation.
2# Copyright (c) 2018 Sylvain Gubian <sylvain.gubian@pmi.com>,
3# Yang Xiang <yang.xiang@pmi.com>
4# Author: Sylvain Gubian, Yang Xiang, PMP S.A.
6"""
7A Dual Annealing global optimization algorithm
8"""
10import numpy as np
11from scipy.optimize import OptimizeResult
12from scipy.optimize import minimize, Bounds
13from scipy.special import gammaln
14from scipy._lib._util import check_random_state
15from scipy.optimize._constraints import new_bounds_to_old
17__all__ = ['dual_annealing']
20class VisitingDistribution:
21 """
22 Class used to generate new coordinates based on the distorted
23 Cauchy-Lorentz distribution. Depending on the steps within the strategy
24 chain, the class implements the strategy for generating new location
25 changes.
27 Parameters
28 ----------
29 lb : array_like
30 A 1-D NumPy ndarray containing lower bounds of the generated
31 components. Neither NaN or inf are allowed.
32 ub : array_like
33 A 1-D NumPy ndarray containing upper bounds for the generated
34 components. Neither NaN or inf are allowed.
35 visiting_param : float
36 Parameter for visiting distribution. Default value is 2.62.
37 Higher values give the visiting distribution a heavier tail, this
38 makes the algorithm jump to a more distant region.
39 The value range is (1, 3]. Its value is fixed for the life of the
40 object.
41 rand_gen : {`~numpy.random.RandomState`, `~numpy.random.Generator`}
42 A `~numpy.random.RandomState`, `~numpy.random.Generator` object
43 for using the current state of the created random generator container.
45 """
46 TAIL_LIMIT = 1.e8
47 MIN_VISIT_BOUND = 1.e-10
49 def __init__(self, lb, ub, visiting_param, rand_gen):
50 # if you wish to make _visiting_param adjustable during the life of
51 # the object then _factor2, _factor3, _factor5, _d1, _factor6 will
52 # have to be dynamically calculated in `visit_fn`. They're factored
53 # out here so they don't need to be recalculated all the time.
54 self._visiting_param = visiting_param
55 self.rand_gen = rand_gen
56 self.lower = lb
57 self.upper = ub
58 self.bound_range = ub - lb
60 # these are invariant numbers unless visiting_param changes
61 self._factor2 = np.exp((4.0 - self._visiting_param) * np.log(
62 self._visiting_param - 1.0))
63 self._factor3 = np.exp((2.0 - self._visiting_param) * np.log(2.0)
64 / (self._visiting_param - 1.0))
65 self._factor4_p = np.sqrt(np.pi) * self._factor2 / (self._factor3 * (
66 3.0 - self._visiting_param))
68 self._factor5 = 1.0 / (self._visiting_param - 1.0) - 0.5
69 self._d1 = 2.0 - self._factor5
70 self._factor6 = np.pi * (1.0 - self._factor5) / np.sin(
71 np.pi * (1.0 - self._factor5)) / np.exp(gammaln(self._d1))
73 def visiting(self, x, step, temperature):
74 """ Based on the step in the strategy chain, new coordinates are
75 generated by changing all components is the same time or only
76 one of them, the new values are computed with visit_fn method
77 """
78 dim = x.size
79 if step < dim:
80 # Changing all coordinates with a new visiting value
81 visits = self.visit_fn(temperature, dim)
82 upper_sample, lower_sample = self.rand_gen.uniform(size=2)
83 visits[visits > self.TAIL_LIMIT] = self.TAIL_LIMIT * upper_sample
84 visits[visits < -self.TAIL_LIMIT] = -self.TAIL_LIMIT * lower_sample
85 x_visit = visits + x
86 a = x_visit - self.lower
87 b = np.fmod(a, self.bound_range) + self.bound_range
88 x_visit = np.fmod(b, self.bound_range) + self.lower
89 x_visit[np.fabs(
90 x_visit - self.lower) < self.MIN_VISIT_BOUND] += 1.e-10
91 else:
92 # Changing only one coordinate at a time based on strategy
93 # chain step
94 x_visit = np.copy(x)
95 visit = self.visit_fn(temperature, 1)[0]
96 if visit > self.TAIL_LIMIT:
97 visit = self.TAIL_LIMIT * self.rand_gen.uniform()
98 elif visit < -self.TAIL_LIMIT:
99 visit = -self.TAIL_LIMIT * self.rand_gen.uniform()
100 index = step - dim
101 x_visit[index] = visit + x[index]
102 a = x_visit[index] - self.lower[index]
103 b = np.fmod(a, self.bound_range[index]) + self.bound_range[index]
104 x_visit[index] = np.fmod(b, self.bound_range[
105 index]) + self.lower[index]
106 if np.fabs(x_visit[index] - self.lower[
107 index]) < self.MIN_VISIT_BOUND:
108 x_visit[index] += self.MIN_VISIT_BOUND
109 return x_visit
111 def visit_fn(self, temperature, dim):
112 """ Formula Visita from p. 405 of reference [2] """
113 x, y = self.rand_gen.normal(size=(dim, 2)).T
115 factor1 = np.exp(np.log(temperature) / (self._visiting_param - 1.0))
116 factor4 = self._factor4_p * factor1
118 # sigmax
119 x *= np.exp(-(self._visiting_param - 1.0) * np.log(
120 self._factor6 / factor4) / (3.0 - self._visiting_param))
122 den = np.exp((self._visiting_param - 1.0) * np.log(np.fabs(y)) /
123 (3.0 - self._visiting_param))
125 return x / den
128class EnergyState:
129 """
130 Class used to record the energy state. At any time, it knows what is the
131 currently used coordinates and the most recent best location.
133 Parameters
134 ----------
135 lower : array_like
136 A 1-D NumPy ndarray containing lower bounds for generating an initial
137 random components in the `reset` method.
138 upper : array_like
139 A 1-D NumPy ndarray containing upper bounds for generating an initial
140 random components in the `reset` method
141 components. Neither NaN or inf are allowed.
142 callback : callable, ``callback(x, f, context)``, optional
143 A callback function which will be called for all minima found.
144 ``x`` and ``f`` are the coordinates and function value of the
145 latest minimum found, and `context` has value in [0, 1, 2]
146 """
147 # Maximum number of trials for generating a valid starting point
148 MAX_REINIT_COUNT = 1000
150 def __init__(self, lower, upper, callback=None):
151 self.ebest = None
152 self.current_energy = None
153 self.current_location = None
154 self.xbest = None
155 self.lower = lower
156 self.upper = upper
157 self.callback = callback
159 def reset(self, func_wrapper, rand_gen, x0=None):
160 """
161 Initialize current location is the search domain. If `x0` is not
162 provided, a random location within the bounds is generated.
163 """
164 if x0 is None:
165 self.current_location = rand_gen.uniform(self.lower, self.upper,
166 size=len(self.lower))
167 else:
168 self.current_location = np.copy(x0)
169 init_error = True
170 reinit_counter = 0
171 while init_error:
172 self.current_energy = func_wrapper.fun(self.current_location)
173 if self.current_energy is None:
174 raise ValueError('Objective function is returning None')
175 if (not np.isfinite(self.current_energy) or np.isnan(
176 self.current_energy)):
177 if reinit_counter >= EnergyState.MAX_REINIT_COUNT:
178 init_error = False
179 message = (
180 'Stopping algorithm because function '
181 'create NaN or (+/-) infinity values even with '
182 'trying new random parameters'
183 )
184 raise ValueError(message)
185 self.current_location = rand_gen.uniform(self.lower,
186 self.upper,
187 size=self.lower.size)
188 reinit_counter += 1
189 else:
190 init_error = False
191 # If first time reset, initialize ebest and xbest
192 if self.ebest is None and self.xbest is None:
193 self.ebest = self.current_energy
194 self.xbest = np.copy(self.current_location)
195 # Otherwise, we keep them in case of reannealing reset
197 def update_best(self, e, x, context):
198 self.ebest = e
199 self.xbest = np.copy(x)
200 if self.callback is not None:
201 val = self.callback(x, e, context)
202 if val is not None:
203 if val:
204 return ('Callback function requested to stop early by '
205 'returning True')
207 def update_current(self, e, x):
208 self.current_energy = e
209 self.current_location = np.copy(x)
212class StrategyChain:
213 """
214 Class that implements within a Markov chain the strategy for location
215 acceptance and local search decision making.
217 Parameters
218 ----------
219 acceptance_param : float
220 Parameter for acceptance distribution. It is used to control the
221 probability of acceptance. The lower the acceptance parameter, the
222 smaller the probability of acceptance. Default value is -5.0 with
223 a range (-1e4, -5].
224 visit_dist : VisitingDistribution
225 Instance of `VisitingDistribution` class.
226 func_wrapper : ObjectiveFunWrapper
227 Instance of `ObjectiveFunWrapper` class.
228 minimizer_wrapper: LocalSearchWrapper
229 Instance of `LocalSearchWrapper` class.
230 rand_gen : {None, int, `numpy.random.Generator`,
231 `numpy.random.RandomState`}, optional
233 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
234 singleton is used.
235 If `seed` is an int, a new ``RandomState`` instance is used,
236 seeded with `seed`.
237 If `seed` is already a ``Generator`` or ``RandomState`` instance then
238 that instance is used.
239 energy_state: EnergyState
240 Instance of `EnergyState` class.
242 """
244 def __init__(self, acceptance_param, visit_dist, func_wrapper,
245 minimizer_wrapper, rand_gen, energy_state):
246 # Local strategy chain minimum energy and location
247 self.emin = energy_state.current_energy
248 self.xmin = np.array(energy_state.current_location)
249 # Global optimizer state
250 self.energy_state = energy_state
251 # Acceptance parameter
252 self.acceptance_param = acceptance_param
253 # Visiting distribution instance
254 self.visit_dist = visit_dist
255 # Wrapper to objective function
256 self.func_wrapper = func_wrapper
257 # Wrapper to the local minimizer
258 self.minimizer_wrapper = minimizer_wrapper
259 self.not_improved_idx = 0
260 self.not_improved_max_idx = 1000
261 self._rand_gen = rand_gen
262 self.temperature_step = 0
263 self.K = 100 * len(energy_state.current_location)
265 def accept_reject(self, j, e, x_visit):
266 r = self._rand_gen.uniform()
267 pqv_temp = 1.0 - ((1.0 - self.acceptance_param) *
268 (e - self.energy_state.current_energy) / self.temperature_step)
269 if pqv_temp <= 0.:
270 pqv = 0.
271 else:
272 pqv = np.exp(np.log(pqv_temp) / (
273 1. - self.acceptance_param))
275 if r <= pqv:
276 # We accept the new location and update state
277 self.energy_state.update_current(e, x_visit)
278 self.xmin = np.copy(self.energy_state.current_location)
280 # No improvement for a long time
281 if self.not_improved_idx >= self.not_improved_max_idx:
282 if j == 0 or self.energy_state.current_energy < self.emin:
283 self.emin = self.energy_state.current_energy
284 self.xmin = np.copy(self.energy_state.current_location)
286 def run(self, step, temperature):
287 self.temperature_step = temperature / float(step + 1)
288 self.not_improved_idx += 1
289 for j in range(self.energy_state.current_location.size * 2):
290 if j == 0:
291 if step == 0:
292 self.energy_state_improved = True
293 else:
294 self.energy_state_improved = False
295 x_visit = self.visit_dist.visiting(
296 self.energy_state.current_location, j, temperature)
297 # Calling the objective function
298 e = self.func_wrapper.fun(x_visit)
299 if e < self.energy_state.current_energy:
300 # We have got a better energy value
301 self.energy_state.update_current(e, x_visit)
302 if e < self.energy_state.ebest:
303 val = self.energy_state.update_best(e, x_visit, 0)
304 if val is not None:
305 if val:
306 return val
307 self.energy_state_improved = True
308 self.not_improved_idx = 0
309 else:
310 # We have not improved but do we accept the new location?
311 self.accept_reject(j, e, x_visit)
312 if self.func_wrapper.nfev >= self.func_wrapper.maxfun:
313 return ('Maximum number of function call reached '
314 'during annealing')
315 # End of StrategyChain loop
317 def local_search(self):
318 # Decision making for performing a local search
319 # based on strategy chain results
320 # If energy has been improved or no improvement since too long,
321 # performing a local search with the best strategy chain location
322 if self.energy_state_improved:
323 # Global energy has improved, let's see if LS improves further
324 e, x = self.minimizer_wrapper.local_search(self.energy_state.xbest,
325 self.energy_state.ebest)
326 if e < self.energy_state.ebest:
327 self.not_improved_idx = 0
328 val = self.energy_state.update_best(e, x, 1)
329 if val is not None:
330 if val:
331 return val
332 self.energy_state.update_current(e, x)
333 if self.func_wrapper.nfev >= self.func_wrapper.maxfun:
334 return ('Maximum number of function call reached '
335 'during local search')
336 # Check probability of a need to perform a LS even if no improvement
337 do_ls = False
338 if self.K < 90 * len(self.energy_state.current_location):
339 pls = np.exp(self.K * (
340 self.energy_state.ebest - self.energy_state.current_energy) /
341 self.temperature_step)
342 if pls >= self._rand_gen.uniform():
343 do_ls = True
344 # Global energy not improved, let's see what LS gives
345 # on the best strategy chain location
346 if self.not_improved_idx >= self.not_improved_max_idx:
347 do_ls = True
348 if do_ls:
349 e, x = self.minimizer_wrapper.local_search(self.xmin, self.emin)
350 self.xmin = np.copy(x)
351 self.emin = e
352 self.not_improved_idx = 0
353 self.not_improved_max_idx = self.energy_state.current_location.size
354 if e < self.energy_state.ebest:
355 val = self.energy_state.update_best(
356 self.emin, self.xmin, 2)
357 if val is not None:
358 if val:
359 return val
360 self.energy_state.update_current(e, x)
361 if self.func_wrapper.nfev >= self.func_wrapper.maxfun:
362 return ('Maximum number of function call reached '
363 'during dual annealing')
366class ObjectiveFunWrapper:
368 def __init__(self, func, maxfun=1e7, *args):
369 self.func = func
370 self.args = args
371 # Number of objective function evaluations
372 self.nfev = 0
373 # Number of gradient function evaluation if used
374 self.ngev = 0
375 # Number of hessian of the objective function if used
376 self.nhev = 0
377 self.maxfun = maxfun
379 def fun(self, x):
380 self.nfev += 1
381 return self.func(x, *self.args)
384class LocalSearchWrapper:
385 """
386 Class used to wrap around the minimizer used for local search
387 Default local minimizer is SciPy minimizer L-BFGS-B
388 """
390 LS_MAXITER_RATIO = 6
391 LS_MAXITER_MIN = 100
392 LS_MAXITER_MAX = 1000
394 def __init__(self, search_bounds, func_wrapper, *args, **kwargs):
395 self.func_wrapper = func_wrapper
396 self.kwargs = kwargs
397 self.jac = self.kwargs.get('jac', None)
398 self.minimizer = minimize
399 bounds_list = list(zip(*search_bounds))
400 self.lower = np.array(bounds_list[0])
401 self.upper = np.array(bounds_list[1])
403 # If no minimizer specified, use SciPy minimize with 'L-BFGS-B' method
404 if not self.kwargs:
405 n = len(self.lower)
406 ls_max_iter = min(max(n * self.LS_MAXITER_RATIO,
407 self.LS_MAXITER_MIN),
408 self.LS_MAXITER_MAX)
409 self.kwargs['method'] = 'L-BFGS-B'
410 self.kwargs['options'] = {
411 'maxiter': ls_max_iter,
412 }
413 self.kwargs['bounds'] = list(zip(self.lower, self.upper))
414 elif callable(self.jac):
415 def wrapped_jac(x):
416 return self.jac(x, *args)
417 self.kwargs['jac'] = wrapped_jac
419 def local_search(self, x, e):
420 # Run local search from the given x location where energy value is e
421 x_tmp = np.copy(x)
422 mres = self.minimizer(self.func_wrapper.fun, x, **self.kwargs)
423 if 'njev' in mres:
424 self.func_wrapper.ngev += mres.njev
425 if 'nhev' in mres:
426 self.func_wrapper.nhev += mres.nhev
427 # Check if is valid value
428 is_finite = np.all(np.isfinite(mres.x)) and np.isfinite(mres.fun)
429 in_bounds = np.all(mres.x >= self.lower) and np.all(
430 mres.x <= self.upper)
431 is_valid = is_finite and in_bounds
433 # Use the new point only if it is valid and return a better results
434 if is_valid and mres.fun < e:
435 return mres.fun, mres.x
436 else:
437 return e, x_tmp
440def dual_annealing(func, bounds, args=(), maxiter=1000,
441 minimizer_kwargs=None, initial_temp=5230.,
442 restart_temp_ratio=2.e-5, visit=2.62, accept=-5.0,
443 maxfun=1e7, seed=None, no_local_search=False,
444 callback=None, x0=None):
445 """
446 Find the global minimum of a function using Dual Annealing.
448 Parameters
449 ----------
450 func : callable
451 The objective function to be minimized. Must be in the form
452 ``f(x, *args)``, where ``x`` is the argument in the form of a 1-D array
453 and ``args`` is a tuple of any additional fixed parameters needed to
454 completely specify the function.
455 bounds : sequence or `Bounds`
456 Bounds for variables. There are two ways to specify the bounds:
458 1. Instance of `Bounds` class.
459 2. Sequence of ``(min, max)`` pairs for each element in `x`.
461 args : tuple, optional
462 Any additional fixed parameters needed to completely specify the
463 objective function.
464 maxiter : int, optional
465 The maximum number of global search iterations. Default value is 1000.
466 minimizer_kwargs : dict, optional
467 Extra keyword arguments to be passed to the local minimizer
468 (`minimize`). Some important options could be:
469 ``method`` for the minimizer method to use and ``args`` for
470 objective function additional arguments.
471 initial_temp : float, optional
472 The initial temperature, use higher values to facilitates a wider
473 search of the energy landscape, allowing dual_annealing to escape
474 local minima that it is trapped in. Default value is 5230. Range is
475 (0.01, 5.e4].
476 restart_temp_ratio : float, optional
477 During the annealing process, temperature is decreasing, when it
478 reaches ``initial_temp * restart_temp_ratio``, the reannealing process
479 is triggered. Default value of the ratio is 2e-5. Range is (0, 1).
480 visit : float, optional
481 Parameter for visiting distribution. Default value is 2.62. Higher
482 values give the visiting distribution a heavier tail, this makes
483 the algorithm jump to a more distant region. The value range is (1, 3].
484 accept : float, optional
485 Parameter for acceptance distribution. It is used to control the
486 probability of acceptance. The lower the acceptance parameter, the
487 smaller the probability of acceptance. Default value is -5.0 with
488 a range (-1e4, -5].
489 maxfun : int, optional
490 Soft limit for the number of objective function calls. If the
491 algorithm is in the middle of a local search, this number will be
492 exceeded, the algorithm will stop just after the local search is
493 done. Default value is 1e7.
494 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
495 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
496 singleton is used.
497 If `seed` is an int, a new ``RandomState`` instance is used,
498 seeded with `seed`.
499 If `seed` is already a ``Generator`` or ``RandomState`` instance then
500 that instance is used.
501 Specify `seed` for repeatable minimizations. The random numbers
502 generated with this seed only affect the visiting distribution function
503 and new coordinates generation.
504 no_local_search : bool, optional
505 If `no_local_search` is set to True, a traditional Generalized
506 Simulated Annealing will be performed with no local search
507 strategy applied.
508 callback : callable, optional
509 A callback function with signature ``callback(x, f, context)``,
510 which will be called for all minima found.
511 ``x`` and ``f`` are the coordinates and function value of the
512 latest minimum found, and ``context`` has value in [0, 1, 2], with the
513 following meaning:
515 - 0: minimum detected in the annealing process.
516 - 1: detection occurred in the local search process.
517 - 2: detection done in the dual annealing process.
519 If the callback implementation returns True, the algorithm will stop.
520 x0 : ndarray, shape(n,), optional
521 Coordinates of a single N-D starting point.
523 Returns
524 -------
525 res : OptimizeResult
526 The optimization result represented as a `OptimizeResult` object.
527 Important attributes are: ``x`` the solution array, ``fun`` the value
528 of the function at the solution, and ``message`` which describes the
529 cause of the termination.
530 See `OptimizeResult` for a description of other attributes.
532 Notes
533 -----
534 This function implements the Dual Annealing optimization. This stochastic
535 approach derived from [3]_ combines the generalization of CSA (Classical
536 Simulated Annealing) and FSA (Fast Simulated Annealing) [1]_ [2]_ coupled
537 to a strategy for applying a local search on accepted locations [4]_.
538 An alternative implementation of this same algorithm is described in [5]_
539 and benchmarks are presented in [6]_. This approach introduces an advanced
540 method to refine the solution found by the generalized annealing
541 process. This algorithm uses a distorted Cauchy-Lorentz visiting
542 distribution, with its shape controlled by the parameter :math:`q_{v}`
544 .. math::
546 g_{q_{v}}(\\Delta x(t)) \\propto \\frac{ \\
547 \\left[T_{q_{v}}(t) \\right]^{-\\frac{D}{3-q_{v}}}}{ \\
548 \\left[{1+(q_{v}-1)\\frac{(\\Delta x(t))^{2}} { \\
549 \\left[T_{q_{v}}(t)\\right]^{\\frac{2}{3-q_{v}}}}}\\right]^{ \\
550 \\frac{1}{q_{v}-1}+\\frac{D-1}{2}}}
552 Where :math:`t` is the artificial time. This visiting distribution is used
553 to generate a trial jump distance :math:`\\Delta x(t)` of variable
554 :math:`x(t)` under artificial temperature :math:`T_{q_{v}}(t)`.
556 From the starting point, after calling the visiting distribution
557 function, the acceptance probability is computed as follows:
559 .. math::
561 p_{q_{a}} = \\min{\\{1,\\left[1-(1-q_{a}) \\beta \\Delta E \\right]^{ \\
562 \\frac{1}{1-q_{a}}}\\}}
564 Where :math:`q_{a}` is a acceptance parameter. For :math:`q_{a}<1`, zero
565 acceptance probability is assigned to the cases where
567 .. math::
569 [1-(1-q_{a}) \\beta \\Delta E] < 0
571 The artificial temperature :math:`T_{q_{v}}(t)` is decreased according to
573 .. math::
575 T_{q_{v}}(t) = T_{q_{v}}(1) \\frac{2^{q_{v}-1}-1}{\\left( \\
576 1 + t\\right)^{q_{v}-1}-1}
578 Where :math:`q_{v}` is the visiting parameter.
580 .. versionadded:: 1.2.0
582 References
583 ----------
584 .. [1] Tsallis C. Possible generalization of Boltzmann-Gibbs
585 statistics. Journal of Statistical Physics, 52, 479-487 (1998).
586 .. [2] Tsallis C, Stariolo DA. Generalized Simulated Annealing.
587 Physica A, 233, 395-406 (1996).
588 .. [3] Xiang Y, Sun DY, Fan W, Gong XG. Generalized Simulated
589 Annealing Algorithm and Its Application to the Thomson Model.
590 Physics Letters A, 233, 216-220 (1997).
591 .. [4] Xiang Y, Gong XG. Efficiency of Generalized Simulated
592 Annealing. Physical Review E, 62, 4473 (2000).
593 .. [5] Xiang Y, Gubian S, Suomela B, Hoeng J. Generalized
594 Simulated Annealing for Efficient Global Optimization: the GenSA
595 Package for R. The R Journal, Volume 5/1 (2013).
596 .. [6] Mullen, K. Continuous Global Optimization in R. Journal of
597 Statistical Software, 60(6), 1 - 45, (2014).
598 :doi:`10.18637/jss.v060.i06`
600 Examples
601 --------
602 The following example is a 10-D problem, with many local minima.
603 The function involved is called Rastrigin
604 (https://en.wikipedia.org/wiki/Rastrigin_function)
606 >>> import numpy as np
607 >>> from scipy.optimize import dual_annealing
608 >>> func = lambda x: np.sum(x*x - 10*np.cos(2*np.pi*x)) + 10*np.size(x)
609 >>> lw = [-5.12] * 10
610 >>> up = [5.12] * 10
611 >>> ret = dual_annealing(func, bounds=list(zip(lw, up)))
612 >>> ret.x
613 array([-4.26437714e-09, -3.91699361e-09, -1.86149218e-09, -3.97165720e-09,
614 -6.29151648e-09, -6.53145322e-09, -3.93616815e-09, -6.55623025e-09,
615 -6.05775280e-09, -5.00668935e-09]) # random
616 >>> ret.fun
617 0.000000
619 """
621 if isinstance(bounds, Bounds):
622 bounds = new_bounds_to_old(bounds.lb, bounds.ub, len(bounds.lb))
624 # noqa: E501
625 if x0 is not None and not len(x0) == len(bounds):
626 raise ValueError('Bounds size does not match x0')
628 lu = list(zip(*bounds))
629 lower = np.array(lu[0])
630 upper = np.array(lu[1])
631 # Check that restart temperature ratio is correct
632 if restart_temp_ratio <= 0. or restart_temp_ratio >= 1.:
633 raise ValueError('Restart temperature ratio has to be in range (0, 1)')
634 # Checking bounds are valid
635 if (np.any(np.isinf(lower)) or np.any(np.isinf(upper)) or np.any(
636 np.isnan(lower)) or np.any(np.isnan(upper))):
637 raise ValueError('Some bounds values are inf values or nan values')
638 # Checking that bounds are consistent
639 if not np.all(lower < upper):
640 raise ValueError('Bounds are not consistent min < max')
641 # Checking that bounds are the same length
642 if not len(lower) == len(upper):
643 raise ValueError('Bounds do not have the same dimensions')
645 # Wrapper for the objective function
646 func_wrapper = ObjectiveFunWrapper(func, maxfun, *args)
648 # minimizer_kwargs has to be a dict, not None
649 minimizer_kwargs = minimizer_kwargs or {}
651 minimizer_wrapper = LocalSearchWrapper(
652 bounds, func_wrapper, *args, **minimizer_kwargs)
654 # Initialization of random Generator for reproducible runs if seed provided
655 rand_state = check_random_state(seed)
656 # Initialization of the energy state
657 energy_state = EnergyState(lower, upper, callback)
658 energy_state.reset(func_wrapper, rand_state, x0)
659 # Minimum value of annealing temperature reached to perform
660 # re-annealing
661 temperature_restart = initial_temp * restart_temp_ratio
662 # VisitingDistribution instance
663 visit_dist = VisitingDistribution(lower, upper, visit, rand_state)
664 # Strategy chain instance
665 strategy_chain = StrategyChain(accept, visit_dist, func_wrapper,
666 minimizer_wrapper, rand_state, energy_state)
667 need_to_stop = False
668 iteration = 0
669 message = []
670 # OptimizeResult object to be returned
671 optimize_res = OptimizeResult()
672 optimize_res.success = True
673 optimize_res.status = 0
675 t1 = np.exp((visit - 1) * np.log(2.0)) - 1.0
676 # Run the search loop
677 while not need_to_stop:
678 for i in range(maxiter):
679 # Compute temperature for this step
680 s = float(i) + 2.0
681 t2 = np.exp((visit - 1) * np.log(s)) - 1.0
682 temperature = initial_temp * t1 / t2
683 if iteration >= maxiter:
684 message.append("Maximum number of iteration reached")
685 need_to_stop = True
686 break
687 # Need a re-annealing process?
688 if temperature < temperature_restart:
689 energy_state.reset(func_wrapper, rand_state)
690 break
691 # starting strategy chain
692 val = strategy_chain.run(i, temperature)
693 if val is not None:
694 message.append(val)
695 need_to_stop = True
696 optimize_res.success = False
697 break
698 # Possible local search at the end of the strategy chain
699 if not no_local_search:
700 val = strategy_chain.local_search()
701 if val is not None:
702 message.append(val)
703 need_to_stop = True
704 optimize_res.success = False
705 break
706 iteration += 1
708 # Setting the OptimizeResult values
709 optimize_res.x = energy_state.xbest
710 optimize_res.fun = energy_state.ebest
711 optimize_res.nit = iteration
712 optimize_res.nfev = func_wrapper.nfev
713 optimize_res.njev = func_wrapper.ngev
714 optimize_res.nhev = func_wrapper.nhev
715 optimize_res.message = message
716 return optimize_res