_basinhopping.py 29 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741
  1. """
  2. basinhopping: The basinhopping global optimization algorithm
  3. """
  4. import numpy as np
  5. import math
  6. import scipy.optimize
  7. from scipy._lib._util import check_random_state
  8. __all__ = ['basinhopping']
  9. class Storage:
  10. """
  11. Class used to store the lowest energy structure
  12. """
  13. def __init__(self, minres):
  14. self._add(minres)
  15. def _add(self, minres):
  16. self.minres = minres
  17. self.minres.x = np.copy(minres.x)
  18. def update(self, minres):
  19. if minres.fun < self.minres.fun:
  20. self._add(minres)
  21. return True
  22. else:
  23. return False
  24. def get_lowest(self):
  25. return self.minres
  26. class BasinHoppingRunner:
  27. """This class implements the core of the basinhopping algorithm.
  28. x0 : ndarray
  29. The starting coordinates.
  30. minimizer : callable
  31. The local minimizer, with signature ``result = minimizer(x)``.
  32. The return value is an `optimize.OptimizeResult` object.
  33. step_taking : callable
  34. This function displaces the coordinates randomly. Signature should
  35. be ``x_new = step_taking(x)``. Note that `x` may be modified in-place.
  36. accept_tests : list of callables
  37. Each test is passed the kwargs `f_new`, `x_new`, `f_old` and
  38. `x_old`. These tests will be used to judge whether or not to accept
  39. the step. The acceptable return values are True, False, or ``"force
  40. accept"``. If any of the tests return False then the step is rejected.
  41. If ``"force accept"``, then this will override any other tests in
  42. order to accept the step. This can be used, for example, to forcefully
  43. escape from a local minimum that ``basinhopping`` is trapped in.
  44. disp : bool, optional
  45. Display status messages.
  46. """
  47. def __init__(self, x0, minimizer, step_taking, accept_tests, disp=False):
  48. self.x = np.copy(x0)
  49. self.minimizer = minimizer
  50. self.step_taking = step_taking
  51. self.accept_tests = accept_tests
  52. self.disp = disp
  53. self.nstep = 0
  54. # initialize return object
  55. self.res = scipy.optimize.OptimizeResult()
  56. self.res.minimization_failures = 0
  57. # do initial minimization
  58. minres = minimizer(self.x)
  59. if not minres.success:
  60. self.res.minimization_failures += 1
  61. if self.disp:
  62. print("warning: basinhopping: local minimization failure")
  63. self.x = np.copy(minres.x)
  64. self.energy = minres.fun
  65. if self.disp:
  66. print("basinhopping step %d: f %g" % (self.nstep, self.energy))
  67. # initialize storage class
  68. self.storage = Storage(minres)
  69. if hasattr(minres, "nfev"):
  70. self.res.nfev = minres.nfev
  71. if hasattr(minres, "njev"):
  72. self.res.njev = minres.njev
  73. if hasattr(minres, "nhev"):
  74. self.res.nhev = minres.nhev
  75. def _monte_carlo_step(self):
  76. """Do one Monte Carlo iteration
  77. Randomly displace the coordinates, minimize, and decide whether
  78. or not to accept the new coordinates.
  79. """
  80. # Take a random step. Make a copy of x because the step_taking
  81. # algorithm might change x in place
  82. x_after_step = np.copy(self.x)
  83. x_after_step = self.step_taking(x_after_step)
  84. # do a local minimization
  85. minres = self.minimizer(x_after_step)
  86. x_after_quench = minres.x
  87. energy_after_quench = minres.fun
  88. if not minres.success:
  89. self.res.minimization_failures += 1
  90. if self.disp:
  91. print("warning: basinhopping: local minimization failure")
  92. if hasattr(minres, "nfev"):
  93. self.res.nfev += minres.nfev
  94. if hasattr(minres, "njev"):
  95. self.res.njev += minres.njev
  96. if hasattr(minres, "nhev"):
  97. self.res.nhev += minres.nhev
  98. # accept the move based on self.accept_tests. If any test is False,
  99. # then reject the step. If any test returns the special string
  100. # 'force accept', then accept the step regardless. This can be used
  101. # to forcefully escape from a local minimum if normal basin hopping
  102. # steps are not sufficient.
  103. accept = True
  104. for test in self.accept_tests:
  105. testres = test(f_new=energy_after_quench, x_new=x_after_quench,
  106. f_old=self.energy, x_old=self.x)
  107. if testres == 'force accept':
  108. accept = True
  109. break
  110. elif testres is None:
  111. raise ValueError("accept_tests must return True, False, or "
  112. "'force accept'")
  113. elif not testres:
  114. accept = False
  115. # Report the result of the acceptance test to the take step class.
  116. # This is for adaptive step taking
  117. if hasattr(self.step_taking, "report"):
  118. self.step_taking.report(accept, f_new=energy_after_quench,
  119. x_new=x_after_quench, f_old=self.energy,
  120. x_old=self.x)
  121. return accept, minres
  122. def one_cycle(self):
  123. """Do one cycle of the basinhopping algorithm
  124. """
  125. self.nstep += 1
  126. new_global_min = False
  127. accept, minres = self._monte_carlo_step()
  128. if accept:
  129. self.energy = minres.fun
  130. self.x = np.copy(minres.x)
  131. new_global_min = self.storage.update(minres)
  132. # print some information
  133. if self.disp:
  134. self.print_report(minres.fun, accept)
  135. if new_global_min:
  136. print("found new global minimum on step %d with function"
  137. " value %g" % (self.nstep, self.energy))
  138. # save some variables as BasinHoppingRunner attributes
  139. self.xtrial = minres.x
  140. self.energy_trial = minres.fun
  141. self.accept = accept
  142. return new_global_min
  143. def print_report(self, energy_trial, accept):
  144. """print a status update"""
  145. minres = self.storage.get_lowest()
  146. print("basinhopping step %d: f %g trial_f %g accepted %d "
  147. " lowest_f %g" % (self.nstep, self.energy, energy_trial,
  148. accept, minres.fun))
  149. class AdaptiveStepsize:
  150. """
  151. Class to implement adaptive stepsize.
  152. This class wraps the step taking class and modifies the stepsize to
  153. ensure the true acceptance rate is as close as possible to the target.
  154. Parameters
  155. ----------
  156. takestep : callable
  157. The step taking routine. Must contain modifiable attribute
  158. takestep.stepsize
  159. accept_rate : float, optional
  160. The target step acceptance rate
  161. interval : int, optional
  162. Interval for how often to update the stepsize
  163. factor : float, optional
  164. The step size is multiplied or divided by this factor upon each
  165. update.
  166. verbose : bool, optional
  167. Print information about each update
  168. """
  169. def __init__(self, takestep, accept_rate=0.5, interval=50, factor=0.9,
  170. verbose=True):
  171. self.takestep = takestep
  172. self.target_accept_rate = accept_rate
  173. self.interval = interval
  174. self.factor = factor
  175. self.verbose = verbose
  176. self.nstep = 0
  177. self.nstep_tot = 0
  178. self.naccept = 0
  179. def __call__(self, x):
  180. return self.take_step(x)
  181. def _adjust_step_size(self):
  182. old_stepsize = self.takestep.stepsize
  183. accept_rate = float(self.naccept) / self.nstep
  184. if accept_rate > self.target_accept_rate:
  185. # We're accepting too many steps. This generally means we're
  186. # trapped in a basin. Take bigger steps.
  187. self.takestep.stepsize /= self.factor
  188. else:
  189. # We're not accepting enough steps. Take smaller steps.
  190. self.takestep.stepsize *= self.factor
  191. if self.verbose:
  192. print("adaptive stepsize: acceptance rate %f target %f new "
  193. "stepsize %g old stepsize %g" % (accept_rate,
  194. self.target_accept_rate, self.takestep.stepsize,
  195. old_stepsize))
  196. def take_step(self, x):
  197. self.nstep += 1
  198. self.nstep_tot += 1
  199. if self.nstep % self.interval == 0:
  200. self._adjust_step_size()
  201. return self.takestep(x)
  202. def report(self, accept, **kwargs):
  203. "called by basinhopping to report the result of the step"
  204. if accept:
  205. self.naccept += 1
  206. class RandomDisplacement:
  207. """Add a random displacement of maximum size `stepsize` to each coordinate.
  208. Calling this updates `x` in-place.
  209. Parameters
  210. ----------
  211. stepsize : float, optional
  212. Maximum stepsize in any dimension
  213. random_gen : {None, int, `numpy.random.Generator`,
  214. `numpy.random.RandomState`}, optional
  215. If `seed` is None (or `np.random`), the `numpy.random.RandomState`
  216. singleton is used.
  217. If `seed` is an int, a new ``RandomState`` instance is used,
  218. seeded with `seed`.
  219. If `seed` is already a ``Generator`` or ``RandomState`` instance then
  220. that instance is used.
  221. """
  222. def __init__(self, stepsize=0.5, random_gen=None):
  223. self.stepsize = stepsize
  224. self.random_gen = check_random_state(random_gen)
  225. def __call__(self, x):
  226. x += self.random_gen.uniform(-self.stepsize, self.stepsize,
  227. np.shape(x))
  228. return x
  229. class MinimizerWrapper:
  230. """
  231. wrap a minimizer function as a minimizer class
  232. """
  233. def __init__(self, minimizer, func=None, **kwargs):
  234. self.minimizer = minimizer
  235. self.func = func
  236. self.kwargs = kwargs
  237. def __call__(self, x0):
  238. if self.func is None:
  239. return self.minimizer(x0, **self.kwargs)
  240. else:
  241. return self.minimizer(self.func, x0, **self.kwargs)
  242. class Metropolis:
  243. """Metropolis acceptance criterion.
  244. Parameters
  245. ----------
  246. T : float
  247. The "temperature" parameter for the accept or reject criterion.
  248. random_gen : {None, int, `numpy.random.Generator`,
  249. `numpy.random.RandomState`}, optional
  250. If `seed` is None (or `np.random`), the `numpy.random.RandomState`
  251. singleton is used.
  252. If `seed` is an int, a new ``RandomState`` instance is used,
  253. seeded with `seed`.
  254. If `seed` is already a ``Generator`` or ``RandomState`` instance then
  255. that instance is used.
  256. Random number generator used for acceptance test.
  257. """
  258. def __init__(self, T, random_gen=None):
  259. # Avoid ZeroDivisionError since "MBH can be regarded as a special case
  260. # of the BH framework with the Metropolis criterion, where temperature
  261. # T = 0." (Reject all steps that increase energy.)
  262. self.beta = 1.0 / T if T != 0 else float('inf')
  263. self.random_gen = check_random_state(random_gen)
  264. def accept_reject(self, energy_new, energy_old):
  265. """
  266. If new energy is lower than old, it will always be accepted.
  267. If new is higher than old, there is a chance it will be accepted,
  268. less likely for larger differences.
  269. """
  270. with np.errstate(invalid='ignore'):
  271. # The energy values being fed to Metropolis are 1-length arrays, and if
  272. # they are equal, their difference is 0, which gets multiplied by beta,
  273. # which is inf, and array([0]) * float('inf') causes
  274. #
  275. # RuntimeWarning: invalid value encountered in multiply
  276. #
  277. # Ignore this warning so when the algorithm is on a flat plane, it always
  278. # accepts the step, to try to move off the plane.
  279. prod = -(energy_new - energy_old) * self.beta
  280. w = math.exp(min(0, prod))
  281. rand = self.random_gen.uniform()
  282. return w >= rand
  283. def __call__(self, **kwargs):
  284. """
  285. f_new and f_old are mandatory in kwargs
  286. """
  287. return bool(self.accept_reject(kwargs["f_new"],
  288. kwargs["f_old"]))
  289. def basinhopping(func, x0, niter=100, T=1.0, stepsize=0.5,
  290. minimizer_kwargs=None, take_step=None, accept_test=None,
  291. callback=None, interval=50, disp=False, niter_success=None,
  292. seed=None, *, target_accept_rate=0.5, stepwise_factor=0.9):
  293. """Find the global minimum of a function using the basin-hopping algorithm.
  294. Basin-hopping is a two-phase method that combines a global stepping
  295. algorithm with local minimization at each step. Designed to mimic
  296. the natural process of energy minimization of clusters of atoms, it works
  297. well for similar problems with "funnel-like, but rugged" energy landscapes
  298. [5]_.
  299. As the step-taking, step acceptance, and minimization methods are all
  300. customizable, this function can also be used to implement other two-phase
  301. methods.
  302. Parameters
  303. ----------
  304. func : callable ``f(x, *args)``
  305. Function to be optimized. ``args`` can be passed as an optional item
  306. in the dict `minimizer_kwargs`
  307. x0 : array_like
  308. Initial guess.
  309. niter : integer, optional
  310. The number of basin-hopping iterations. There will be a total of
  311. ``niter + 1`` runs of the local minimizer.
  312. T : float, optional
  313. The "temperature" parameter for the acceptance or rejection criterion.
  314. Higher "temperatures" mean that larger jumps in function value will be
  315. accepted. For best results `T` should be comparable to the
  316. separation (in function value) between local minima.
  317. stepsize : float, optional
  318. Maximum step size for use in the random displacement.
  319. minimizer_kwargs : dict, optional
  320. Extra keyword arguments to be passed to the local minimizer
  321. `scipy.optimize.minimize` Some important options could be:
  322. method : str
  323. The minimization method (e.g. ``"L-BFGS-B"``)
  324. args : tuple
  325. Extra arguments passed to the objective function (`func`) and
  326. its derivatives (Jacobian, Hessian).
  327. take_step : callable ``take_step(x)``, optional
  328. Replace the default step-taking routine with this routine. The default
  329. step-taking routine is a random displacement of the coordinates, but
  330. other step-taking algorithms may be better for some systems.
  331. `take_step` can optionally have the attribute ``take_step.stepsize``.
  332. If this attribute exists, then `basinhopping` will adjust
  333. ``take_step.stepsize`` in order to try to optimize the global minimum
  334. search.
  335. accept_test : callable, ``accept_test(f_new=f_new, x_new=x_new, f_old=fold, x_old=x_old)``, optional
  336. Define a test which will be used to judge whether to accept the
  337. step. This will be used in addition to the Metropolis test based on
  338. "temperature" `T`. The acceptable return values are True,
  339. False, or ``"force accept"``. If any of the tests return False
  340. then the step is rejected. If the latter, then this will override any
  341. other tests in order to accept the step. This can be used, for example,
  342. to forcefully escape from a local minimum that `basinhopping` is
  343. trapped in.
  344. callback : callable, ``callback(x, f, accept)``, optional
  345. A callback function which will be called for all minima found. ``x``
  346. and ``f`` are the coordinates and function value of the trial minimum,
  347. and ``accept`` is whether that minimum was accepted. This can
  348. be used, for example, to save the lowest N minima found. Also,
  349. `callback` can be used to specify a user defined stop criterion by
  350. optionally returning True to stop the `basinhopping` routine.
  351. interval : integer, optional
  352. interval for how often to update the `stepsize`
  353. disp : bool, optional
  354. Set to True to print status messages
  355. niter_success : integer, optional
  356. Stop the run if the global minimum candidate remains the same for this
  357. number of iterations.
  358. seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
  359. If `seed` is None (or `np.random`), the `numpy.random.RandomState`
  360. singleton is used.
  361. If `seed` is an int, a new ``RandomState`` instance is used,
  362. seeded with `seed`.
  363. If `seed` is already a ``Generator`` or ``RandomState`` instance then
  364. that instance is used.
  365. Specify `seed` for repeatable minimizations. The random numbers
  366. generated with this seed only affect the default Metropolis
  367. `accept_test` and the default `take_step`. If you supply your own
  368. `take_step` and `accept_test`, and these functions use random
  369. number generation, then those functions are responsible for the state
  370. of their random number generator.
  371. target_accept_rate : float, optional
  372. The target acceptance rate that is used to adjust the `stepsize`.
  373. If the current acceptance rate is greater than the target,
  374. then the `stepsize` is increased. Otherwise, it is decreased.
  375. Range is (0, 1). Default is 0.5.
  376. .. versionadded:: 1.8.0
  377. stepwise_factor : float, optional
  378. The `stepsize` is multiplied or divided by this stepwise factor upon
  379. each update. Range is (0, 1). Default is 0.9.
  380. .. versionadded:: 1.8.0
  381. Returns
  382. -------
  383. res : OptimizeResult
  384. The optimization result represented as a `OptimizeResult` object.
  385. Important attributes are: ``x`` the solution array, ``fun`` the value
  386. of the function at the solution, and ``message`` which describes the
  387. cause of the termination. The ``OptimizeResult`` object returned by the
  388. selected minimizer at the lowest minimum is also contained within this
  389. object and can be accessed through the ``lowest_optimization_result``
  390. attribute. See `OptimizeResult` for a description of other attributes.
  391. See Also
  392. --------
  393. minimize :
  394. The local minimization function called once for each basinhopping step.
  395. `minimizer_kwargs` is passed to this routine.
  396. Notes
  397. -----
  398. Basin-hopping is a stochastic algorithm which attempts to find the global
  399. minimum of a smooth scalar function of one or more variables [1]_ [2]_ [3]_
  400. [4]_. The algorithm in its current form was described by David Wales and
  401. Jonathan Doye [2]_ http://www-wales.ch.cam.ac.uk/.
  402. The algorithm is iterative with each cycle composed of the following
  403. features
  404. 1) random perturbation of the coordinates
  405. 2) local minimization
  406. 3) accept or reject the new coordinates based on the minimized function
  407. value
  408. The acceptance test used here is the Metropolis criterion of standard Monte
  409. Carlo algorithms, although there are many other possibilities [3]_.
  410. This global minimization method has been shown to be extremely efficient
  411. for a wide variety of problems in physics and chemistry. It is
  412. particularly useful when the function has many minima separated by large
  413. barriers. See the `Cambridge Cluster Database
  414. <https://www-wales.ch.cam.ac.uk/CCD.html>`_ for databases of molecular
  415. systems that have been optimized primarily using basin-hopping. This
  416. database includes minimization problems exceeding 300 degrees of freedom.
  417. See the free software program `GMIN <https://www-wales.ch.cam.ac.uk/GMIN>`_
  418. for a Fortran implementation of basin-hopping. This implementation has many
  419. variations of the procedure described above, including more
  420. advanced step taking algorithms and alternate acceptance criterion.
  421. For stochastic global optimization there is no way to determine if the true
  422. global minimum has actually been found. Instead, as a consistency check,
  423. the algorithm can be run from a number of different random starting points
  424. to ensure the lowest minimum found in each example has converged to the
  425. global minimum. For this reason, `basinhopping` will by default simply
  426. run for the number of iterations `niter` and return the lowest minimum
  427. found. It is left to the user to ensure that this is in fact the global
  428. minimum.
  429. Choosing `stepsize`: This is a crucial parameter in `basinhopping` and
  430. depends on the problem being solved. The step is chosen uniformly in the
  431. region from x0-stepsize to x0+stepsize, in each dimension. Ideally, it
  432. should be comparable to the typical separation (in argument values) between
  433. local minima of the function being optimized. `basinhopping` will, by
  434. default, adjust `stepsize` to find an optimal value, but this may take
  435. many iterations. You will get quicker results if you set a sensible
  436. initial value for ``stepsize``.
  437. Choosing `T`: The parameter `T` is the "temperature" used in the
  438. Metropolis criterion. Basinhopping steps are always accepted if
  439. ``func(xnew) < func(xold)``. Otherwise, they are accepted with
  440. probability::
  441. exp( -(func(xnew) - func(xold)) / T )
  442. So, for best results, `T` should to be comparable to the typical
  443. difference (in function values) between local minima. (The height of
  444. "walls" between local minima is irrelevant.)
  445. If `T` is 0, the algorithm becomes Monotonic Basin-Hopping, in which all
  446. steps that increase energy are rejected.
  447. .. versionadded:: 0.12.0
  448. References
  449. ----------
  450. .. [1] Wales, David J. 2003, Energy Landscapes, Cambridge University Press,
  451. Cambridge, UK.
  452. .. [2] Wales, D J, and Doye J P K, Global Optimization by Basin-Hopping and
  453. the Lowest Energy Structures of Lennard-Jones Clusters Containing up to
  454. 110 Atoms. Journal of Physical Chemistry A, 1997, 101, 5111.
  455. .. [3] Li, Z. and Scheraga, H. A., Monte Carlo-minimization approach to the
  456. multiple-minima problem in protein folding, Proc. Natl. Acad. Sci. USA,
  457. 1987, 84, 6611.
  458. .. [4] Wales, D. J. and Scheraga, H. A., Global optimization of clusters,
  459. crystals, and biomolecules, Science, 1999, 285, 1368.
  460. .. [5] Olson, B., Hashmi, I., Molloy, K., and Shehu1, A., Basin Hopping as
  461. a General and Versatile Optimization Framework for the Characterization
  462. of Biological Macromolecules, Advances in Artificial Intelligence,
  463. Volume 2012 (2012), Article ID 674832, :doi:`10.1155/2012/674832`
  464. Examples
  465. --------
  466. The following example is a 1-D minimization problem, with many
  467. local minima superimposed on a parabola.
  468. >>> import numpy as np
  469. >>> from scipy.optimize import basinhopping
  470. >>> func = lambda x: np.cos(14.5 * x - 0.3) + (x + 0.2) * x
  471. >>> x0 = [1.]
  472. Basinhopping, internally, uses a local minimization algorithm. We will use
  473. the parameter `minimizer_kwargs` to tell basinhopping which algorithm to
  474. use and how to set up that minimizer. This parameter will be passed to
  475. `scipy.optimize.minimize`.
  476. >>> minimizer_kwargs = {"method": "BFGS"}
  477. >>> ret = basinhopping(func, x0, minimizer_kwargs=minimizer_kwargs,
  478. ... niter=200)
  479. >>> print("global minimum: x = %.4f, f(x) = %.4f" % (ret.x, ret.fun))
  480. global minimum: x = -0.1951, f(x) = -1.0009
  481. Next consider a 2-D minimization problem. Also, this time, we
  482. will use gradient information to significantly speed up the search.
  483. >>> def func2d(x):
  484. ... f = np.cos(14.5 * x[0] - 0.3) + (x[1] + 0.2) * x[1] + (x[0] +
  485. ... 0.2) * x[0]
  486. ... df = np.zeros(2)
  487. ... df[0] = -14.5 * np.sin(14.5 * x[0] - 0.3) + 2. * x[0] + 0.2
  488. ... df[1] = 2. * x[1] + 0.2
  489. ... return f, df
  490. We'll also use a different local minimization algorithm. Also, we must tell
  491. the minimizer that our function returns both energy and gradient (Jacobian).
  492. >>> minimizer_kwargs = {"method":"L-BFGS-B", "jac":True}
  493. >>> x0 = [1.0, 1.0]
  494. >>> ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs,
  495. ... niter=200)
  496. >>> print("global minimum: x = [%.4f, %.4f], f(x) = %.4f" % (ret.x[0],
  497. ... ret.x[1],
  498. ... ret.fun))
  499. global minimum: x = [-0.1951, -0.1000], f(x) = -1.0109
  500. Here is an example using a custom step-taking routine. Imagine you want
  501. the first coordinate to take larger steps than the rest of the coordinates.
  502. This can be implemented like so:
  503. >>> class MyTakeStep:
  504. ... def __init__(self, stepsize=0.5):
  505. ... self.stepsize = stepsize
  506. ... self.rng = np.random.default_rng()
  507. ... def __call__(self, x):
  508. ... s = self.stepsize
  509. ... x[0] += self.rng.uniform(-2.*s, 2.*s)
  510. ... x[1:] += self.rng.uniform(-s, s, x[1:].shape)
  511. ... return x
  512. Since ``MyTakeStep.stepsize`` exists basinhopping will adjust the magnitude
  513. of `stepsize` to optimize the search. We'll use the same 2-D function as
  514. before
  515. >>> mytakestep = MyTakeStep()
  516. >>> ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs,
  517. ... niter=200, take_step=mytakestep)
  518. >>> print("global minimum: x = [%.4f, %.4f], f(x) = %.4f" % (ret.x[0],
  519. ... ret.x[1],
  520. ... ret.fun))
  521. global minimum: x = [-0.1951, -0.1000], f(x) = -1.0109
  522. Now, let's do an example using a custom callback function which prints the
  523. value of every minimum found
  524. >>> def print_fun(x, f, accepted):
  525. ... print("at minimum %.4f accepted %d" % (f, int(accepted)))
  526. We'll run it for only 10 basinhopping steps this time.
  527. >>> rng = np.random.default_rng()
  528. >>> ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs,
  529. ... niter=10, callback=print_fun, seed=rng)
  530. at minimum 0.4159 accepted 1
  531. at minimum -0.4317 accepted 1
  532. at minimum -1.0109 accepted 1
  533. at minimum -0.9073 accepted 1
  534. at minimum -0.4317 accepted 0
  535. at minimum -0.1021 accepted 1
  536. at minimum -0.7425 accepted 1
  537. at minimum -0.9073 accepted 1
  538. at minimum -0.4317 accepted 0
  539. at minimum -0.7425 accepted 1
  540. at minimum -0.9073 accepted 1
  541. The minimum at -1.0109 is actually the global minimum, found already on the
  542. 8th iteration.
  543. """
  544. if target_accept_rate <= 0. or target_accept_rate >= 1.:
  545. raise ValueError('target_accept_rate has to be in range (0, 1)')
  546. if stepwise_factor <= 0. or stepwise_factor >= 1.:
  547. raise ValueError('stepwise_factor has to be in range (0, 1)')
  548. x0 = np.array(x0)
  549. # set up the np.random generator
  550. rng = check_random_state(seed)
  551. # set up minimizer
  552. if minimizer_kwargs is None:
  553. minimizer_kwargs = dict()
  554. wrapped_minimizer = MinimizerWrapper(scipy.optimize.minimize, func,
  555. **minimizer_kwargs)
  556. # set up step-taking algorithm
  557. if take_step is not None:
  558. if not callable(take_step):
  559. raise TypeError("take_step must be callable")
  560. # if take_step.stepsize exists then use AdaptiveStepsize to control
  561. # take_step.stepsize
  562. if hasattr(take_step, "stepsize"):
  563. take_step_wrapped = AdaptiveStepsize(
  564. take_step, interval=interval,
  565. accept_rate=target_accept_rate,
  566. factor=stepwise_factor,
  567. verbose=disp)
  568. else:
  569. take_step_wrapped = take_step
  570. else:
  571. # use default
  572. displace = RandomDisplacement(stepsize=stepsize, random_gen=rng)
  573. take_step_wrapped = AdaptiveStepsize(displace, interval=interval,
  574. accept_rate=target_accept_rate,
  575. factor=stepwise_factor,
  576. verbose=disp)
  577. # set up accept tests
  578. accept_tests = []
  579. if accept_test is not None:
  580. if not callable(accept_test):
  581. raise TypeError("accept_test must be callable")
  582. accept_tests = [accept_test]
  583. # use default
  584. metropolis = Metropolis(T, random_gen=rng)
  585. accept_tests.append(metropolis)
  586. if niter_success is None:
  587. niter_success = niter + 2
  588. bh = BasinHoppingRunner(x0, wrapped_minimizer, take_step_wrapped,
  589. accept_tests, disp=disp)
  590. # The wrapped minimizer is called once during construction of
  591. # BasinHoppingRunner, so run the callback
  592. if callable(callback):
  593. callback(bh.storage.minres.x, bh.storage.minres.fun, True)
  594. # start main iteration loop
  595. count, i = 0, 0
  596. message = ["requested number of basinhopping iterations completed"
  597. " successfully"]
  598. for i in range(niter):
  599. new_global_min = bh.one_cycle()
  600. if callable(callback):
  601. # should we pass a copy of x?
  602. val = callback(bh.xtrial, bh.energy_trial, bh.accept)
  603. if val is not None:
  604. if val:
  605. message = ["callback function requested stop early by"
  606. "returning True"]
  607. break
  608. count += 1
  609. if new_global_min:
  610. count = 0
  611. elif count > niter_success:
  612. message = ["success condition satisfied"]
  613. break
  614. # prepare return object
  615. res = bh.res
  616. res.lowest_optimization_result = bh.storage.get_lowest()
  617. res.x = np.copy(res.lowest_optimization_result.x)
  618. res.fun = res.lowest_optimization_result.fun
  619. res.message = message
  620. res.nit = i + 1
  621. res.success = res.lowest_optimization_result.success
  622. return res