Hide keyboard shortcuts

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# 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. 

5 

6""" 

7A Dual Annealing global optimization algorithm 

8""" 

9 

10import numpy as np 

11from scipy.optimize import OptimizeResult 

12from scipy.optimize import minimize 

13from scipy.special import gammaln 

14from scipy._lib._util import check_random_state 

15 

16 

17__all__ = ['dual_annealing'] 

18 

19 

20class VisitingDistribution(object): 

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. 

26 

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 (0, 3]. It's 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. 

44 """ 

45 TAIL_LIMIT = 1.e8 

46 MIN_VISIT_BOUND = 1.e-10 

47 

48 def __init__(self, lb, ub, visiting_param, rand_gen): 

49 # if you wish to make _visiting_param adjustable during the life of 

50 # the object then _factor2, _factor3, _factor5, _d1, _factor6 will 

51 # have to be dynamically calculated in `visit_fn`. They're factored 

52 # out here so they don't need to be recalculated all the time. 

53 self._visiting_param = visiting_param 

54 self.rand_gen = rand_gen 

55 self.lower = lb 

56 self.upper = ub 

57 self.bound_range = ub - lb 

58 

59 # these are invariant numbers unless visiting_param changes 

60 self._factor2 = np.exp((4.0 - self._visiting_param) * np.log( 

61 self._visiting_param - 1.0)) 

62 self._factor3 = np.exp((2.0 - self._visiting_param) * np.log(2.0) 

63 / (self._visiting_param - 1.0)) 

64 self._factor4_p = np.sqrt(np.pi) * self._factor2 / (self._factor3 * ( 

65 3.0 - self._visiting_param)) 

66 

67 self._factor5 = 1.0 / (self._visiting_param - 1.0) - 0.5 

68 self._d1 = 2.0 - self._factor5 

69 self._factor6 = np.pi * (1.0 - self._factor5) / np.sin( 

70 np.pi * (1.0 - self._factor5)) / np.exp(gammaln(self._d1)) 

71 

72 def visiting(self, x, step, temperature): 

73 """ Based on the step in the strategy chain, new coordinated are 

74 generated by changing all components is the same time or only 

75 one of them, the new values are computed with visit_fn method 

76 """ 

77 dim = x.size 

78 if step < dim: 

79 # Changing all coordinates with a new visiting value 

80 visits = self.visit_fn(temperature, dim) 

81 upper_sample, lower_sample = self.rand_gen.uniform(size=2) 

82 visits[visits > self.TAIL_LIMIT] = self.TAIL_LIMIT * upper_sample 

83 visits[visits < -self.TAIL_LIMIT] = -self.TAIL_LIMIT * lower_sample 

84 x_visit = visits + x 

85 a = x_visit - self.lower 

86 b = np.fmod(a, self.bound_range) + self.bound_range 

87 x_visit = np.fmod(b, self.bound_range) + self.lower 

88 x_visit[np.fabs( 

89 x_visit - self.lower) < self.MIN_VISIT_BOUND] += 1.e-10 

90 else: 

91 # Changing only one coordinate at a time based on strategy 

92 # chain step 

93 x_visit = np.copy(x) 

94 visit = self.visit_fn(temperature, 1) 

95 if visit > self.TAIL_LIMIT: 

96 visit = self.TAIL_LIMIT * self.rand_gen.uniform() 

97 elif visit < -self.TAIL_LIMIT: 

98 visit = -self.TAIL_LIMIT * self.rand_gen.uniform() 

99 index = step - dim 

100 x_visit[index] = visit + x[index] 

101 a = x_visit[index] - self.lower[index] 

102 b = np.fmod(a, self.bound_range[index]) + self.bound_range[index] 

103 x_visit[index] = np.fmod(b, self.bound_range[ 

104 index]) + self.lower[index] 

105 if np.fabs(x_visit[index] - self.lower[ 

106 index]) < self.MIN_VISIT_BOUND: 

107 x_visit[index] += self.MIN_VISIT_BOUND 

108 return x_visit 

109 

110 def visit_fn(self, temperature, dim): 

111 """ Formula Visita from p. 405 of reference [2] """ 

112 x, y = self.rand_gen.normal(size=(dim, 2)).T 

113 

114 factor1 = np.exp(np.log(temperature) / (self._visiting_param - 1.0)) 

115 factor4 = self._factor4_p * factor1 

116 

117 # sigmax 

118 x *= np.exp(-(self._visiting_param - 1.0) * np.log( 

119 self._factor6 / factor4) / (3.0 - self._visiting_param)) 

120 

121 den = np.exp((self._visiting_param - 1.0) * np.log(np.fabs(y)) / 

122 (3.0 - self._visiting_param)) 

123 

124 return x / den 

125 

126 

127class EnergyState(object): 

128 """ 

129 Class used to record the energy state. At any time, it knows what is the 

130 currently used coordinates and the most recent best location. 

131 

132 Parameters 

133 ---------- 

134 lower : array_like 

135 A 1-D NumPy ndarray containing lower bounds for generating an initial 

136 random components in the `reset` method. 

137 upper : array_like 

138 A 1-D NumPy ndarray containing upper bounds for generating an initial 

139 random components in the `reset` method 

140 components. Neither NaN or inf are allowed. 

141 callback : callable, ``callback(x, f, context)``, optional 

142 A callback function which will be called for all minima found. 

143 ``x`` and ``f`` are the coordinates and function value of the 

144 latest minimum found, and `context` has value in [0, 1, 2] 

145 """ 

146 # Maximimum number of trials for generating a valid starting point 

147 MAX_REINIT_COUNT = 1000 

148 

149 def __init__(self, lower, upper, callback=None): 

150 self.ebest = None 

151 self.current_energy = None 

152 self.current_location = None 

153 self.xbest = None 

154 self.lower = lower 

155 self.upper = upper 

156 self.callback = callback 

157 

158 def reset(self, func_wrapper, rand_gen, x0=None): 

159 """ 

160 Initialize current location is the search domain. If `x0` is not 

161 provided, a random location within the bounds is generated. 

162 """ 

163 if x0 is None: 

164 self.current_location = rand_gen.uniform(self.lower, self.upper, 

165 size=len(self.lower)) 

166 else: 

167 self.current_location = np.copy(x0) 

168 init_error = True 

169 reinit_counter = 0 

170 while init_error: 

171 self.current_energy = func_wrapper.fun(self.current_location) 

172 if self.current_energy is None: 

173 raise ValueError('Objective function is returning None') 

174 if (not np.isfinite(self.current_energy) or np.isnan( 

175 self.current_energy)): 

176 if reinit_counter >= EnergyState.MAX_REINIT_COUNT: 

177 init_error = False 

178 message = ( 

179 'Stopping algorithm because function ' 

180 'create NaN or (+/-) infinity values even with ' 

181 'trying new random parameters' 

182 ) 

183 raise ValueError(message) 

184 self.current_location = rand_gen.uniform(self.lower, 

185 self.upper, 

186 size=self.lower.size) 

187 reinit_counter += 1 

188 else: 

189 init_error = False 

190 # If first time reset, initialize ebest and xbest 

191 if self.ebest is None and self.xbest is None: 

192 self.ebest = self.current_energy 

193 self.xbest = np.copy(self.current_location) 

194 # Otherwise, we keep them in case of reannealing reset 

195 

196 def update_best(self, e, x, context): 

197 self.ebest = e 

198 self.xbest = np.copy(x) 

199 if self.callback is not None: 

200 val = self.callback(x, e, context) 

201 if val is not None: 

202 if val: 

203 return('Callback function requested to stop early by ' 

204 'returning True') 

205 

206 def update_current(self, e, x): 

207 self.current_energy = e 

208 self.current_location = np.copy(x) 

209 

210 

211class StrategyChain(object): 

212 """ 

213 Class that implements within a Markov chain the strategy for location 

214 acceptance and local search decision making. 

215 

216 Parameters 

217 ---------- 

218 acceptance_param : float 

219 Parameter for acceptance distribution. It is used to control the 

220 probability of acceptance. The lower the acceptance parameter, the 

221 smaller the probability of acceptance. Default value is -5.0 with 

222 a range (-1e4, -5]. 

223 visit_dist : VisitingDistribution 

224 Instance of `VisitingDistribution` class. 

225 func_wrapper : ObjectiveFunWrapper 

226 Instance of `ObjectiveFunWrapper` class. 

227 minimizer_wrapper: LocalSearchWrapper 

228 Instance of `LocalSearchWrapper` class. 

229 rand_gen : {`~numpy.random.RandomState`, `~numpy.random.Generator`} 

230 A `~numpy.random.RandomState` or `~numpy.random.Generator` 

231 object for using the current state of the created random generator 

232 container. 

233 energy_state: EnergyState 

234 Instance of `EnergyState` class. 

235 """ 

236 def __init__(self, acceptance_param, visit_dist, func_wrapper, 

237 minimizer_wrapper, rand_gen, energy_state): 

238 # Local strategy chain minimum energy and location 

239 self.emin = energy_state.current_energy 

240 self.xmin = np.array(energy_state.current_location) 

241 # Global optimizer state 

242 self.energy_state = energy_state 

243 # Acceptance parameter 

244 self.acceptance_param = acceptance_param 

245 # Visiting distribution instance 

246 self.visit_dist = visit_dist 

247 # Wrapper to objective function 

248 self.func_wrapper = func_wrapper 

249 # Wrapper to the local minimizer 

250 self.minimizer_wrapper = minimizer_wrapper 

251 self.not_improved_idx = 0 

252 self.not_improved_max_idx = 1000 

253 self._rand_gen = rand_gen 

254 self.temperature_step = 0 

255 self.K = 100 * len(energy_state.current_location) 

256 

257 def accept_reject(self, j, e, x_visit): 

258 r = self._rand_gen.uniform() 

259 pqv_temp = (self.acceptance_param - 1.0) * ( 

260 e - self.energy_state.current_energy) / ( 

261 self.temperature_step + 1.) 

262 if pqv_temp <= 0.: 

263 pqv = 0. 

264 else: 

265 pqv = np.exp(np.log(pqv_temp) / ( 

266 1. - self.acceptance_param)) 

267 if r <= pqv: 

268 # We accept the new location and update state 

269 self.energy_state.update_current(e, x_visit) 

270 self.xmin = np.copy(self.energy_state.current_location) 

271 

272 # No improvement for a long time 

273 if self.not_improved_idx >= self.not_improved_max_idx: 

274 if j == 0 or self.energy_state.current_energy < self.emin: 

275 self.emin = self.energy_state.current_energy 

276 self.xmin = np.copy(self.energy_state.current_location) 

277 

278 def run(self, step, temperature): 

279 self.temperature_step = temperature / float(step + 1) 

280 self.not_improved_idx += 1 

281 for j in range(self.energy_state.current_location.size * 2): 

282 if j == 0: 

283 if step == 0: 

284 self.energy_state_improved = True 

285 else: 

286 self.energy_state_improved = False 

287 x_visit = self.visit_dist.visiting( 

288 self.energy_state.current_location, j, temperature) 

289 # Calling the objective function 

290 e = self.func_wrapper.fun(x_visit) 

291 if e < self.energy_state.current_energy: 

292 # We have got a better energy value 

293 self.energy_state.update_current(e, x_visit) 

294 if e < self.energy_state.ebest: 

295 val = self.energy_state.update_best(e, x_visit, 0) 

296 if val is not None: 

297 if val: 

298 return val 

299 self.energy_state_improved = True 

300 self.not_improved_idx = 0 

301 else: 

302 # We have not improved but do we accept the new location? 

303 self.accept_reject(j, e, x_visit) 

304 if self.func_wrapper.nfev >= self.func_wrapper.maxfun: 

305 return ('Maximum number of function call reached ' 

306 'during annealing') 

307 # End of StrategyChain loop 

308 

309 def local_search(self): 

310 # Decision making for performing a local search 

311 # based on strategy chain results 

312 # If energy has been improved or no improvement since too long, 

313 # performing a local search with the best strategy chain location 

314 if self.energy_state_improved: 

315 # Global energy has improved, let's see if LS improves further 

316 e, x = self.minimizer_wrapper.local_search(self.energy_state.xbest, 

317 self.energy_state.ebest) 

318 if e < self.energy_state.ebest: 

319 self.not_improved_idx = 0 

320 val = self.energy_state.update_best(e, x, 1) 

321 if val is not None: 

322 if val: 

323 return val 

324 self.energy_state.update_current(e, x) 

325 if self.func_wrapper.nfev >= self.func_wrapper.maxfun: 

326 return ('Maximum number of function call reached ' 

327 'during local search') 

328 # Check probability of a need to perform a LS even if no improvement 

329 do_ls = False 

330 if self.K < 90 * len(self.energy_state.current_location): 

331 pls = np.exp(self.K * ( 

332 self.energy_state.ebest - self.energy_state.current_energy) / 

333 self.temperature_step) 

334 if pls >= self._rand_gen.uniform(): 

335 do_ls = True 

336 # Global energy not improved, let's see what LS gives 

337 # on the best strategy chain location 

338 if self.not_improved_idx >= self.not_improved_max_idx: 

339 do_ls = True 

340 if do_ls: 

341 e, x = self.minimizer_wrapper.local_search(self.xmin, self.emin) 

342 self.xmin = np.copy(x) 

343 self.emin = e 

344 self.not_improved_idx = 0 

345 self.not_improved_max_idx = self.energy_state.current_location.size 

346 if e < self.energy_state.ebest: 

347 val = self.energy_state.update_best( 

348 self.emin, self.xmin, 2) 

349 if val is not None: 

350 if val: 

351 return val 

352 self.energy_state.update_current(e, x) 

353 if self.func_wrapper.nfev >= self.func_wrapper.maxfun: 

354 return ('Maximum number of function call reached ' 

355 'during dual annealing') 

356 

357 

358class ObjectiveFunWrapper(object): 

359 

360 def __init__(self, func, maxfun=1e7, *args): 

361 self.func = func 

362 self.args = args 

363 # Number of objective function evaluations 

364 self.nfev = 0 

365 # Number of gradient function evaluation if used 

366 self.ngev = 0 

367 # Number of hessian of the objective function if used 

368 self.nhev = 0 

369 self.maxfun = maxfun 

370 

371 def fun(self, x): 

372 self.nfev += 1 

373 return self.func(x, *self.args) 

374 

375 

376class LocalSearchWrapper(object): 

377 """ 

378 Class used to wrap around the minimizer used for local search 

379 Default local minimizer is SciPy minimizer L-BFGS-B 

380 """ 

381 

382 LS_MAXITER_RATIO = 6 

383 LS_MAXITER_MIN = 100 

384 LS_MAXITER_MAX = 1000 

385 

386 def __init__(self, bounds, func_wrapper, **kwargs): 

387 self.func_wrapper = func_wrapper 

388 self.kwargs = kwargs 

389 self.minimizer = minimize 

390 bounds_list = list(zip(*bounds)) 

391 self.lower = np.array(bounds_list[0]) 

392 self.upper = np.array(bounds_list[1]) 

393 

394 # If no minimizer specified, use SciPy minimize with 'L-BFGS-B' method 

395 if not self.kwargs: 

396 n = len(self.lower) 

397 ls_max_iter = min(max(n * self.LS_MAXITER_RATIO, 

398 self.LS_MAXITER_MIN), 

399 self.LS_MAXITER_MAX) 

400 self.kwargs['method'] = 'L-BFGS-B' 

401 self.kwargs['options'] = { 

402 'maxiter': ls_max_iter, 

403 } 

404 self.kwargs['bounds'] = list(zip(self.lower, self.upper)) 

405 

406 def local_search(self, x, e): 

407 # Run local search from the given x location where energy value is e 

408 x_tmp = np.copy(x) 

409 mres = self.minimizer(self.func_wrapper.fun, x, **self.kwargs) 

410 if 'njev' in mres.keys(): 

411 self.func_wrapper.ngev += mres.njev 

412 if 'nhev' in mres.keys(): 

413 self.func_wrapper.nhev += mres.nhev 

414 # Check if is valid value 

415 is_finite = np.all(np.isfinite(mres.x)) and np.isfinite(mres.fun) 

416 in_bounds = np.all(mres.x >= self.lower) and np.all( 

417 mres.x <= self.upper) 

418 is_valid = is_finite and in_bounds 

419 

420 # Use the new point only if it is valid and return a better results 

421 if is_valid and mres.fun < e: 

422 return mres.fun, mres.x 

423 else: 

424 return e, x_tmp 

425 

426 

427def dual_annealing(func, bounds, args=(), maxiter=1000, 

428 local_search_options={}, initial_temp=5230., 

429 restart_temp_ratio=2.e-5, visit=2.62, accept=-5.0, 

430 maxfun=1e7, seed=None, no_local_search=False, 

431 callback=None, x0=None): 

432 """ 

433 Find the global minimum of a function using Dual Annealing. 

434 

435 Parameters 

436 ---------- 

437 func : callable 

438 The objective function to be minimized. Must be in the form 

439 ``f(x, *args)``, where ``x`` is the argument in the form of a 1-D array 

440 and ``args`` is a tuple of any additional fixed parameters needed to 

441 completely specify the function. 

442 bounds : sequence, shape (n, 2) 

443 Bounds for variables. ``(min, max)`` pairs for each element in ``x``, 

444 defining bounds for the objective function parameter. 

445 args : tuple, optional 

446 Any additional fixed parameters needed to completely specify the 

447 objective function. 

448 maxiter : int, optional 

449 The maximum number of global search iterations. Default value is 1000. 

450 local_search_options : dict, optional 

451 Extra keyword arguments to be passed to the local minimizer 

452 (`minimize`). Some important options could be: 

453 ``method`` for the minimizer method to use and ``args`` for 

454 objective function additional arguments. 

455 initial_temp : float, optional 

456 The initial temperature, use higher values to facilitates a wider 

457 search of the energy landscape, allowing dual_annealing to escape 

458 local minima that it is trapped in. Default value is 5230. Range is 

459 (0.01, 5.e4]. 

460 restart_temp_ratio : float, optional 

461 During the annealing process, temperature is decreasing, when it 

462 reaches ``initial_temp * restart_temp_ratio``, the reannealing process 

463 is triggered. Default value of the ratio is 2e-5. Range is (0, 1). 

464 visit : float, optional 

465 Parameter for visiting distribution. Default value is 2.62. Higher 

466 values give the visiting distribution a heavier tail, this makes 

467 the algorithm jump to a more distant region. The value range is (0, 3]. 

468 accept : float, optional 

469 Parameter for acceptance distribution. It is used to control the 

470 probability of acceptance. The lower the acceptance parameter, the 

471 smaller the probability of acceptance. Default value is -5.0 with 

472 a range (-1e4, -5]. 

473 maxfun : int, optional 

474 Soft limit for the number of objective function calls. If the 

475 algorithm is in the middle of a local search, this number will be 

476 exceeded, the algorithm will stop just after the local search is 

477 done. Default value is 1e7. 

478 seed : {int, `~numpy.random.RandomState`, `~numpy.random.Generator`}, optional 

479 If `seed` is not specified the `~numpy.random.RandomState` singleton is 

480 used. 

481 If `seed` is an int, a new ``RandomState`` instance is used, seeded 

482 with `seed`. 

483 If `seed` is already a ``RandomState`` or ``Generator`` instance, then 

484 that instance is used. 

485 Specify `seed` for repeatable minimizations. The random numbers 

486 generated with this seed only affect the visiting distribution function 

487 and new coordinates generation. 

488 no_local_search : bool, optional 

489 If `no_local_search` is set to True, a traditional Generalized 

490 Simulated Annealing will be performed with no local search 

491 strategy applied. 

492 callback : callable, optional 

493 A callback function with signature ``callback(x, f, context)``, 

494 which will be called for all minima found. 

495 ``x`` and ``f`` are the coordinates and function value of the 

496 latest minimum found, and ``context`` has value in [0, 1, 2], with the 

497 following meaning: 

498 

499 - 0: minimum detected in the annealing process. 

500 - 1: detection occurred in the local search process. 

501 - 2: detection done in the dual annealing process. 

502 

503 If the callback implementation returns True, the algorithm will stop. 

504 x0 : ndarray, shape(n,), optional 

505 Coordinates of a single N-D starting point. 

506 

507 Returns 

508 ------- 

509 res : OptimizeResult 

510 The optimization result represented as a `OptimizeResult` object. 

511 Important attributes are: ``x`` the solution array, ``fun`` the value 

512 of the function at the solution, and ``message`` which describes the 

513 cause of the termination. 

514 See `OptimizeResult` for a description of other attributes. 

515 

516 Notes 

517 ----- 

518 This function implements the Dual Annealing optimization. This stochastic 

519 approach derived from [3]_ combines the generalization of CSA (Classical 

520 Simulated Annealing) and FSA (Fast Simulated Annealing) [1]_ [2]_ coupled 

521 to a strategy for applying a local search on accepted locations [4]_. 

522 An alternative implementation of this same algorithm is described in [5]_ 

523 and benchmarks are presented in [6]_. This approach introduces an advanced 

524 method to refine the solution found by the generalized annealing 

525 process. This algorithm uses a distorted Cauchy-Lorentz visiting 

526 distribution, with its shape controlled by the parameter :math:`q_{v}` 

527 

528 .. math:: 

529 

530 g_{q_{v}}(\\Delta x(t)) \\propto \\frac{ \\ 

531 \\left[T_{q_{v}}(t) \\right]^{-\\frac{D}{3-q_{v}}}}{ \\ 

532 \\left[{1+(q_{v}-1)\\frac{(\\Delta x(t))^{2}} { \\ 

533 \\left[T_{q_{v}}(t)\\right]^{\\frac{2}{3-q_{v}}}}}\\right]^{ \\ 

534 \\frac{1}{q_{v}-1}+\\frac{D-1}{2}}} 

535 

536 Where :math:`t` is the artificial time. This visiting distribution is used 

537 to generate a trial jump distance :math:`\\Delta x(t)` of variable 

538 :math:`x(t)` under artificial temperature :math:`T_{q_{v}}(t)`. 

539 

540 From the starting point, after calling the visiting distribution 

541 function, the acceptance probability is computed as follows: 

542 

543 .. math:: 

544 

545 p_{q_{a}} = \\min{\\{1,\\left[1-(1-q_{a}) \\beta \\Delta E \\right]^{ \\ 

546 \\frac{1}{1-q_{a}}}\\}} 

547 

548 Where :math:`q_{a}` is a acceptance parameter. For :math:`q_{a}<1`, zero 

549 acceptance probability is assigned to the cases where 

550 

551 .. math:: 

552 

553 [1-(1-q_{a}) \\beta \\Delta E] < 0 

554 

555 The artificial temperature :math:`T_{q_{v}}(t)` is decreased according to 

556 

557 .. math:: 

558 

559 T_{q_{v}}(t) = T_{q_{v}}(1) \\frac{2^{q_{v}-1}-1}{\\left( \\ 

560 1 + t\\right)^{q_{v}-1}-1} 

561 

562 Where :math:`q_{v}` is the visiting parameter. 

563 

564 .. versionadded:: 1.2.0 

565 

566 References 

567 ---------- 

568 .. [1] Tsallis C. Possible generalization of Boltzmann-Gibbs 

569 statistics. Journal of Statistical Physics, 52, 479-487 (1998). 

570 .. [2] Tsallis C, Stariolo DA. Generalized Simulated Annealing. 

571 Physica A, 233, 395-406 (1996). 

572 .. [3] Xiang Y, Sun DY, Fan W, Gong XG. Generalized Simulated 

573 Annealing Algorithm and Its Application to the Thomson Model. 

574 Physics Letters A, 233, 216-220 (1997). 

575 .. [4] Xiang Y, Gong XG. Efficiency of Generalized Simulated 

576 Annealing. Physical Review E, 62, 4473 (2000). 

577 .. [5] Xiang Y, Gubian S, Suomela B, Hoeng J. Generalized 

578 Simulated Annealing for Efficient Global Optimization: the GenSA 

579 Package for R. The R Journal, Volume 5/1 (2013). 

580 .. [6] Mullen, K. Continuous Global Optimization in R. Journal of 

581 Statistical Software, 60(6), 1 - 45, (2014). DOI:10.18637/jss.v060.i06 

582 

583 Examples 

584 -------- 

585 The following example is a 10-D problem, with many local minima. 

586 The function involved is called Rastrigin 

587 (https://en.wikipedia.org/wiki/Rastrigin_function) 

588 

589 >>> from scipy.optimize import dual_annealing 

590 >>> func = lambda x: np.sum(x*x - 10*np.cos(2*np.pi*x)) + 10*np.size(x) 

591 >>> lw = [-5.12] * 10 

592 >>> up = [5.12] * 10 

593 >>> ret = dual_annealing(func, bounds=list(zip(lw, up)), seed=1234) 

594 >>> ret.x 

595 array([-4.26437714e-09, -3.91699361e-09, -1.86149218e-09, -3.97165720e-09, 

596 -6.29151648e-09, -6.53145322e-09, -3.93616815e-09, -6.55623025e-09, 

597 -6.05775280e-09, -5.00668935e-09]) # may vary 

598 >>> ret.fun 

599 0.000000 

600 

601 """ # noqa: E501 

602 if x0 is not None and not len(x0) == len(bounds): 

603 raise ValueError('Bounds size does not match x0') 

604 

605 lu = list(zip(*bounds)) 

606 lower = np.array(lu[0]) 

607 upper = np.array(lu[1]) 

608 # Check that restart temperature ratio is correct 

609 if restart_temp_ratio <= 0. or restart_temp_ratio >= 1.: 

610 raise ValueError('Restart temperature ratio has to be in range (0, 1)') 

611 # Checking bounds are valid 

612 if (np.any(np.isinf(lower)) or np.any(np.isinf(upper)) or np.any( 

613 np.isnan(lower)) or np.any(np.isnan(upper))): 

614 raise ValueError('Some bounds values are inf values or nan values') 

615 # Checking that bounds are consistent 

616 if not np.all(lower < upper): 

617 raise ValueError('Bounds are not consistent min < max') 

618 # Checking that bounds are the same length 

619 if not len(lower) == len(upper): 

620 raise ValueError('Bounds do not have the same dimensions') 

621 

622 # Wrapper for the objective function 

623 func_wrapper = ObjectiveFunWrapper(func, maxfun, *args) 

624 # Wrapper fot the minimizer 

625 minimizer_wrapper = LocalSearchWrapper( 

626 bounds, func_wrapper, **local_search_options) 

627 # Initialization of RandomState for reproducible runs if seed provided 

628 rand_state = check_random_state(seed) 

629 # Initialization of the energy state 

630 energy_state = EnergyState(lower, upper, callback) 

631 energy_state.reset(func_wrapper, rand_state, x0) 

632 # Minimum value of annealing temperature reached to perform 

633 # re-annealing 

634 temperature_restart = initial_temp * restart_temp_ratio 

635 # VisitingDistribution instance 

636 visit_dist = VisitingDistribution(lower, upper, visit, rand_state) 

637 # Strategy chain instance 

638 strategy_chain = StrategyChain(accept, visit_dist, func_wrapper, 

639 minimizer_wrapper, rand_state, energy_state) 

640 need_to_stop = False 

641 iteration = 0 

642 message = [] 

643 # OptimizeResult object to be returned 

644 optimize_res = OptimizeResult() 

645 optimize_res.success = True 

646 optimize_res.status = 0 

647 

648 t1 = np.exp((visit - 1) * np.log(2.0)) - 1.0 

649 # Run the search loop 

650 while(not need_to_stop): 

651 for i in range(maxiter): 

652 # Compute temperature for this step 

653 s = float(i) + 2.0 

654 t2 = np.exp((visit - 1) * np.log(s)) - 1.0 

655 temperature = initial_temp * t1 / t2 

656 if iteration >= maxiter: 

657 message.append("Maximum number of iteration reached") 

658 need_to_stop = True 

659 break 

660 # Need a re-annealing process? 

661 if temperature < temperature_restart: 

662 energy_state.reset(func_wrapper, rand_state) 

663 break 

664 # starting strategy chain 

665 val = strategy_chain.run(i, temperature) 

666 if val is not None: 

667 message.append(val) 

668 need_to_stop = True 

669 optimize_res.success = False 

670 break 

671 # Possible local search at the end of the strategy chain 

672 if not no_local_search: 

673 val = strategy_chain.local_search() 

674 if val is not None: 

675 message.append(val) 

676 need_to_stop = True 

677 optimize_res.success = False 

678 break 

679 iteration += 1 

680 

681 # Setting the OptimizeResult values 

682 optimize_res.x = energy_state.xbest 

683 optimize_res.fun = energy_state.ebest 

684 optimize_res.nit = iteration 

685 optimize_res.nfev = func_wrapper.nfev 

686 optimize_res.njev = func_wrapper.ngev 

687 optimize_res.nhev = func_wrapper.nhev 

688 optimize_res.message = message 

689 return optimize_res