_basinhopping.py 29 KB

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