_tnc.py 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438
  1. # TNC Python interface
  2. # @(#) $Jeannot: tnc.py,v 1.11 2005/01/28 18:27:31 js Exp $
  3. # Copyright (c) 2004-2005, Jean-Sebastien Roy (js@jeannot.org)
  4. # Permission is hereby granted, free of charge, to any person obtaining a
  5. # copy of this software and associated documentation files (the
  6. # "Software"), to deal in the Software without restriction, including
  7. # without limitation the rights to use, copy, modify, merge, publish,
  8. # distribute, sublicense, and/or sell copies of the Software, and to
  9. # permit persons to whom the Software is furnished to do so, subject to
  10. # the following conditions:
  11. # The above copyright notice and this permission notice shall be included
  12. # in all copies or substantial portions of the Software.
  13. # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  14. # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  15. # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  16. # IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
  17. # CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  18. # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  19. # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  20. """
  21. TNC: A Python interface to the TNC non-linear optimizer
  22. TNC is a non-linear optimizer. To use it, you must provide a function to
  23. minimize. The function must take one argument: the list of coordinates where to
  24. evaluate the function; and it must return either a tuple, whose first element is the
  25. value of the function, and whose second argument is the gradient of the function
  26. (as a list of values); or None, to abort the minimization.
  27. """
  28. from scipy.optimize import _moduleTNC as moduleTNC
  29. from ._optimize import (MemoizeJac, OptimizeResult, _check_unknown_options,
  30. _prepare_scalar_function)
  31. from ._constraints import old_bound_to_new
  32. from scipy._lib._array_api import array_namespace
  33. from scipy._lib import array_api_extra as xpx
  34. from numpy import inf, array, zeros
  35. __all__ = ['fmin_tnc']
  36. MSG_NONE = 0 # No messages
  37. MSG_ITER = 1 # One line per iteration
  38. MSG_INFO = 2 # Informational messages
  39. MSG_VERS = 4 # Version info
  40. MSG_EXIT = 8 # Exit reasons
  41. MSG_ALL = MSG_ITER + MSG_INFO + MSG_VERS + MSG_EXIT
  42. MSGS = {
  43. MSG_NONE: "No messages",
  44. MSG_ITER: "One line per iteration",
  45. MSG_INFO: "Informational messages",
  46. MSG_VERS: "Version info",
  47. MSG_EXIT: "Exit reasons",
  48. MSG_ALL: "All messages"
  49. }
  50. INFEASIBLE = -1 # Infeasible (lower bound > upper bound)
  51. LOCALMINIMUM = 0 # Local minimum reached (|pg| ~= 0)
  52. FCONVERGED = 1 # Converged (|f_n-f_(n-1)| ~= 0)
  53. XCONVERGED = 2 # Converged (|x_n-x_(n-1)| ~= 0)
  54. MAXFUN = 3 # Max. number of function evaluations reached
  55. LSFAIL = 4 # Linear search failed
  56. CONSTANT = 5 # All lower bounds are equal to the upper bounds
  57. NOPROGRESS = 6 # Unable to progress
  58. USERABORT = 7 # User requested end of minimization
  59. RCSTRINGS = {
  60. INFEASIBLE: "Infeasible (lower bound > upper bound)",
  61. LOCALMINIMUM: "Local minimum reached (|pg| ~= 0)",
  62. FCONVERGED: "Converged (|f_n-f_(n-1)| ~= 0)",
  63. XCONVERGED: "Converged (|x_n-x_(n-1)| ~= 0)",
  64. MAXFUN: "Max. number of function evaluations reached",
  65. LSFAIL: "Linear search failed",
  66. CONSTANT: "All lower bounds are equal to the upper bounds",
  67. NOPROGRESS: "Unable to progress",
  68. USERABORT: "User requested end of minimization"
  69. }
  70. # Changes to interface made by Travis Oliphant, Apr. 2004 for inclusion in
  71. # SciPy
  72. def fmin_tnc(func, x0, fprime=None, args=(), approx_grad=0,
  73. bounds=None, epsilon=1e-8, scale=None, offset=None,
  74. messages=MSG_ALL, maxCGit=-1, maxfun=None, eta=-1,
  75. stepmx=0, accuracy=0, fmin=0, ftol=-1, xtol=-1, pgtol=-1,
  76. rescale=-1, disp=None, callback=None):
  77. r"""
  78. Minimize a function with variables subject to bounds, using
  79. gradient information in a truncated Newton algorithm. This
  80. method wraps a C implementation of the algorithm.
  81. Parameters
  82. ----------
  83. func : callable ``func(x, *args)``
  84. Function to minimize. Must do one of:
  85. 1. Return f and g, where f is the value of the function and g its
  86. gradient (a list of floats).
  87. 2. Return the function value but supply gradient function
  88. separately as `fprime`.
  89. 3. Return the function value and set ``approx_grad=True``.
  90. If the function returns None, the minimization
  91. is aborted.
  92. x0 : array_like
  93. Initial estimate of minimum.
  94. fprime : callable ``fprime(x, *args)``, optional
  95. Gradient of `func`. If None, then either `func` must return the
  96. function value and the gradient (``f,g = func(x, *args)``)
  97. or `approx_grad` must be True.
  98. args : tuple, optional
  99. Arguments to pass to function.
  100. approx_grad : bool, optional
  101. If true, approximate the gradient numerically.
  102. bounds : list, optional
  103. (min, max) pairs for each element in x0, defining the
  104. bounds on that parameter. Use None or +/-inf for one of
  105. min or max when there is no bound in that direction.
  106. epsilon : float, optional
  107. Used if approx_grad is True. The stepsize in a finite
  108. difference approximation for fprime.
  109. scale : array_like, optional
  110. Scaling factors to apply to each variable. If None, the
  111. factors are up-low for interval bounded variables and
  112. 1+|x| for the others. Defaults to None.
  113. offset : array_like, optional
  114. Value to subtract from each variable. If None, the
  115. offsets are (up+low)/2 for interval bounded variables
  116. and x for the others.
  117. messages : int, optional
  118. Bit mask used to select messages display during
  119. minimization values defined in the MSGS dict. Defaults to
  120. MGS_ALL.
  121. disp : int, optional
  122. Integer interface to messages. 0 = no message, 5 = all messages
  123. maxCGit : int, optional
  124. Maximum number of hessian*vector evaluations per main
  125. iteration. If maxCGit == 0, the direction chosen is
  126. -gradient if maxCGit < 0, maxCGit is set to
  127. max(1,min(50,n/2)). Defaults to -1.
  128. maxfun : int, optional
  129. Maximum number of function evaluation. If None, maxfun is
  130. set to max(100, 10*len(x0)). Defaults to None. Note that this function
  131. may violate the limit because of evaluating gradients by numerical
  132. differentiation.
  133. eta : float, optional
  134. Severity of the line search. If < 0 or > 1, set to 0.25.
  135. Defaults to -1.
  136. stepmx : float, optional
  137. Maximum step for the line search. May be increased during
  138. call. If too small, it will be set to 10.0. Defaults to 0.
  139. accuracy : float, optional
  140. Relative precision for finite difference calculations. If
  141. <= machine_precision, set to sqrt(machine_precision).
  142. Defaults to 0.
  143. fmin : float, optional
  144. Minimum function value estimate. Defaults to 0.
  145. ftol : float, optional
  146. Precision goal for the value of f in the stopping criterion.
  147. If ftol < 0.0, ftol is set to 0.0 defaults to -1.
  148. xtol : float, optional
  149. Precision goal for the value of x in the stopping
  150. criterion (after applying x scaling factors). If xtol <
  151. 0.0, xtol is set to sqrt(machine_precision). Defaults to
  152. -1.
  153. pgtol : float, optional
  154. Precision goal for the value of the projected gradient in
  155. the stopping criterion (after applying x scaling factors).
  156. If pgtol < 0.0, pgtol is set to 1e-2 * sqrt(accuracy).
  157. Setting it to 0.0 is not recommended. Defaults to -1.
  158. rescale : float, optional
  159. Scaling factor (in log10) used to trigger f value
  160. rescaling. If 0, rescale at each iteration. If a large
  161. value, never rescale. If < 0, rescale is set to 1.3.
  162. callback : callable, optional
  163. Called after each iteration, as callback(xk), where xk is the
  164. current parameter vector.
  165. Returns
  166. -------
  167. x : ndarray
  168. The solution.
  169. nfeval : int
  170. The number of function evaluations.
  171. rc : int
  172. Return code, see below
  173. See also
  174. --------
  175. minimize: Interface to minimization algorithms for multivariate
  176. functions. See the 'TNC' `method` in particular.
  177. Notes
  178. -----
  179. The underlying algorithm is truncated Newton, also called
  180. Newton Conjugate-Gradient. This method differs from
  181. scipy.optimize.fmin_ncg in that
  182. 1. it wraps a C implementation of the algorithm
  183. 2. it allows each variable to be given an upper and lower bound.
  184. The algorithm incorporates the bound constraints by determining
  185. the descent direction as in an unconstrained truncated Newton,
  186. but never taking a step-size large enough to leave the space
  187. of feasible x's. The algorithm keeps track of a set of
  188. currently active constraints, and ignores them when computing
  189. the minimum allowable step size. (The x's associated with the
  190. active constraint are kept fixed.) If the maximum allowable
  191. step size is zero then a new constraint is added. At the end
  192. of each iteration one of the constraints may be deemed no
  193. longer active and removed. A constraint is considered
  194. no longer active is if it is currently active
  195. but the gradient for that variable points inward from the
  196. constraint. The specific constraint removed is the one
  197. associated with the variable of largest index whose
  198. constraint is no longer active.
  199. Return codes are defined as follows:
  200. - ``-1`` : Infeasible (lower bound > upper bound)
  201. - ``0`` : Local minimum reached (:math:`|pg| \approx 0`)
  202. - ``1`` : Converged (:math:`|f_n-f_(n-1)| \approx 0`)
  203. - ``2`` : Converged (:math:`|x_n-x_(n-1)| \approx 0`)
  204. - ``3`` : Max. number of function evaluations reached
  205. - ``4`` : Linear search failed
  206. - ``5`` : All lower bounds are equal to the upper bounds
  207. - ``6`` : Unable to progress
  208. - ``7`` : User requested end of minimization
  209. References
  210. ----------
  211. Wright S., Nocedal J. (2006), 'Numerical Optimization'
  212. Nash S.G. (1984), "Newton-Type Minimization Via the Lanczos Method",
  213. SIAM Journal of Numerical Analysis 21, pp. 770-778
  214. """
  215. # handle fprime/approx_grad
  216. if approx_grad:
  217. fun = func
  218. jac = None
  219. elif fprime is None:
  220. fun = MemoizeJac(func)
  221. jac = fun.derivative
  222. else:
  223. fun = func
  224. jac = fprime
  225. if disp is not None: # disp takes precedence over messages
  226. mesg_num = disp
  227. else:
  228. mesg_num = {0:MSG_NONE, 1:MSG_ITER, 2:MSG_INFO, 3:MSG_VERS,
  229. 4:MSG_EXIT, 5:MSG_ALL}.get(messages, MSG_ALL)
  230. # build options
  231. opts = {'eps': epsilon,
  232. 'scale': scale,
  233. 'offset': offset,
  234. 'mesg_num': mesg_num,
  235. 'maxCGit': maxCGit,
  236. 'maxfun': maxfun,
  237. 'eta': eta,
  238. 'stepmx': stepmx,
  239. 'accuracy': accuracy,
  240. 'minfev': fmin,
  241. 'ftol': ftol,
  242. 'xtol': xtol,
  243. 'gtol': pgtol,
  244. 'rescale': rescale,
  245. 'disp': False}
  246. res = _minimize_tnc(fun, x0, args, jac, bounds, callback=callback, **opts)
  247. return res['x'], res['nfev'], res['status']
  248. def _minimize_tnc(fun, x0, args=(), jac=None, bounds=None,
  249. eps=1e-8, scale=None, offset=None, mesg_num=None,
  250. maxCGit=-1, eta=-1, stepmx=0, accuracy=0,
  251. minfev=0, ftol=-1, xtol=-1, gtol=-1, rescale=-1, disp=False,
  252. callback=None, finite_diff_rel_step=None, maxfun=None,
  253. workers=None,
  254. **unknown_options):
  255. """
  256. Minimize a scalar function of one or more variables using a truncated
  257. Newton (TNC) algorithm.
  258. Options
  259. -------
  260. eps : float or ndarray
  261. If `jac is None` the absolute step size used for numerical
  262. approximation of the jacobian via forward differences.
  263. scale : list of floats
  264. Scaling factors to apply to each variable. If None, the
  265. factors are up-low for interval bounded variables and
  266. 1+|x] for the others. Defaults to None.
  267. offset : float
  268. Value to subtract from each variable. If None, the
  269. offsets are (up+low)/2 for interval bounded variables
  270. and x for the others.
  271. disp : bool
  272. Set to True to print convergence messages.
  273. maxCGit : int
  274. Maximum number of hessian*vector evaluations per main
  275. iteration. If maxCGit == 0, the direction chosen is
  276. -gradient if maxCGit < 0, maxCGit is set to
  277. max(1,min(50,n/2)). Defaults to -1.
  278. eta : float
  279. Severity of the line search. If < 0 or > 1, set to 0.25.
  280. Defaults to -1.
  281. stepmx : float
  282. Maximum step for the line search. May be increased during
  283. call. If too small, it will be set to 10.0. Defaults to 0.
  284. accuracy : float
  285. Relative precision for finite difference calculations. If
  286. <= machine_precision, set to sqrt(machine_precision).
  287. Defaults to 0.
  288. minfev : float
  289. Minimum function value estimate. Defaults to 0.
  290. ftol : float
  291. Precision goal for the value of f in the stopping criterion.
  292. If ftol < 0.0, ftol is set to 0.0 defaults to -1.
  293. xtol : float
  294. Precision goal for the value of x in the stopping
  295. criterion (after applying x scaling factors). If xtol <
  296. 0.0, xtol is set to sqrt(machine_precision). Defaults to
  297. -1.
  298. gtol : float
  299. Precision goal for the value of the projected gradient in
  300. the stopping criterion (after applying x scaling factors).
  301. If gtol < 0.0, gtol is set to 1e-2 * sqrt(accuracy).
  302. Setting it to 0.0 is not recommended. Defaults to -1.
  303. rescale : float
  304. Scaling factor (in log10) used to trigger f value
  305. rescaling. If 0, rescale at each iteration. If a large
  306. value, never rescale. If < 0, rescale is set to 1.3.
  307. finite_diff_rel_step : None or array_like, optional
  308. If ``jac in ['2-point', '3-point', 'cs']`` the relative step size to
  309. use for numerical approximation of the jacobian. The absolute step
  310. size is computed as ``h = rel_step * sign(x) * max(1, abs(x))``,
  311. possibly adjusted to fit into the bounds. For ``method='3-point'``
  312. the sign of `h` is ignored. If None (default) then step is selected
  313. automatically.
  314. maxfun : int
  315. Maximum number of function evaluations. If None, `maxfun` is
  316. set to max(100, 10*len(x0)). Defaults to None.
  317. workers : int, map-like callable, optional
  318. A map-like callable, such as `multiprocessing.Pool.map` for evaluating
  319. any numerical differentiation in parallel.
  320. This evaluation is carried out as ``workers(fun, iterable)``.
  321. .. versionadded:: 1.16.0
  322. """
  323. _check_unknown_options(unknown_options)
  324. fmin = minfev
  325. pgtol = gtol
  326. xp = array_namespace(x0)
  327. x0 = xpx.atleast_nd(xp.asarray(x0), ndim=1, xp=xp)
  328. dtype = xp.float64
  329. if xp.isdtype(x0.dtype, "real floating"):
  330. dtype = x0.dtype
  331. x0 = xp.reshape(xp.astype(x0, dtype), -1)
  332. n = len(x0)
  333. if bounds is None:
  334. bounds = [(None,None)] * n
  335. if len(bounds) != n:
  336. raise ValueError('length of x0 != length of bounds')
  337. new_bounds = old_bound_to_new(bounds)
  338. if mesg_num is not None:
  339. messages = {0:MSG_NONE, 1:MSG_ITER, 2:MSG_INFO, 3:MSG_VERS,
  340. 4:MSG_EXIT, 5:MSG_ALL}.get(mesg_num, MSG_ALL)
  341. elif disp:
  342. messages = MSG_ALL
  343. else:
  344. messages = MSG_NONE
  345. sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
  346. finite_diff_rel_step=finite_diff_rel_step,
  347. bounds=new_bounds, workers=workers)
  348. func_and_grad = sf.fun_and_grad
  349. """
  350. low, up : the bounds (lists of floats)
  351. if low is None, the lower bounds are removed.
  352. if up is None, the upper bounds are removed.
  353. low and up defaults to None
  354. """
  355. low = zeros(n)
  356. up = zeros(n)
  357. for i in range(n):
  358. if bounds[i] is None:
  359. l, u = -inf, inf
  360. else:
  361. l,u = bounds[i]
  362. if l is None:
  363. low[i] = -inf
  364. else:
  365. low[i] = l
  366. if u is None:
  367. up[i] = inf
  368. else:
  369. up[i] = u
  370. if scale is None:
  371. scale = array([])
  372. if offset is None:
  373. offset = array([])
  374. if maxfun is None:
  375. maxfun = max(100, 10*len(x0))
  376. rc, nf, nit, x, funv, jacv = moduleTNC.tnc_minimize(
  377. func_and_grad, x0, low, up, scale,
  378. offset, messages, maxCGit, maxfun,
  379. eta, stepmx, accuracy, fmin, ftol,
  380. xtol, pgtol, rescale, callback
  381. )
  382. # the TNC documentation states: "On output, x, f and g may be very
  383. # slightly out of sync because of scaling". Therefore re-evaluate
  384. # func_and_grad so they are synced.
  385. funv, jacv = func_and_grad(x)
  386. return OptimizeResult(x=x, fun=funv, jac=jacv, nfev=sf.nfev,
  387. nit=nit, status=rc, message=RCSTRINGS[rc],
  388. success=(-1 < rc < 3))