_cobyla_py.py 9.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293
  1. """
  2. Interface to Constrained Optimization By Linear Approximation
  3. Functions
  4. ---------
  5. .. autosummary::
  6. :toctree: generated/
  7. fmin_cobyla
  8. """
  9. import functools
  10. from threading import RLock
  11. import numpy as np
  12. from scipy.optimize import _cobyla as cobyla
  13. from ._optimize import OptimizeResult, _check_unknown_options
  14. try:
  15. from itertools import izip
  16. except ImportError:
  17. izip = zip
  18. __all__ = ['fmin_cobyla']
  19. # Workarund as _cobyla.minimize is not threadsafe
  20. # due to an unknown f2py bug and can segfault,
  21. # see gh-9658.
  22. _module_lock = RLock()
  23. def synchronized(func):
  24. @functools.wraps(func)
  25. def wrapper(*args, **kwargs):
  26. with _module_lock:
  27. return func(*args, **kwargs)
  28. return wrapper
  29. @synchronized
  30. def fmin_cobyla(func, x0, cons, args=(), consargs=None, rhobeg=1.0,
  31. rhoend=1e-4, maxfun=1000, disp=None, catol=2e-4,
  32. *, callback=None):
  33. """
  34. Minimize a function using the Constrained Optimization By Linear
  35. Approximation (COBYLA) method. This method wraps a FORTRAN
  36. implementation of the algorithm.
  37. Parameters
  38. ----------
  39. func : callable
  40. Function to minimize. In the form func(x, \\*args).
  41. x0 : ndarray
  42. Initial guess.
  43. cons : sequence
  44. Constraint functions; must all be ``>=0`` (a single function
  45. if only 1 constraint). Each function takes the parameters `x`
  46. as its first argument, and it can return either a single number or
  47. an array or list of numbers.
  48. args : tuple, optional
  49. Extra arguments to pass to function.
  50. consargs : tuple, optional
  51. Extra arguments to pass to constraint functions (default of None means
  52. use same extra arguments as those passed to func).
  53. Use ``()`` for no extra arguments.
  54. rhobeg : float, optional
  55. Reasonable initial changes to the variables.
  56. rhoend : float, optional
  57. Final accuracy in the optimization (not precisely guaranteed). This
  58. is a lower bound on the size of the trust region.
  59. disp : {0, 1, 2, 3}, optional
  60. Controls the frequency of output; 0 implies no output.
  61. maxfun : int, optional
  62. Maximum number of function evaluations.
  63. catol : float, optional
  64. Absolute tolerance for constraint violations.
  65. callback : callable, optional
  66. Called after each iteration, as ``callback(x)``, where ``x`` is the
  67. current parameter vector.
  68. Returns
  69. -------
  70. x : ndarray
  71. The argument that minimises `f`.
  72. See also
  73. --------
  74. minimize: Interface to minimization algorithms for multivariate
  75. functions. See the 'COBYLA' `method` in particular.
  76. Notes
  77. -----
  78. This algorithm is based on linear approximations to the objective
  79. function and each constraint. We briefly describe the algorithm.
  80. Suppose the function is being minimized over k variables. At the
  81. jth iteration the algorithm has k+1 points v_1, ..., v_(k+1),
  82. an approximate solution x_j, and a radius RHO_j.
  83. (i.e., linear plus a constant) approximations to the objective
  84. function and constraint functions such that their function values
  85. agree with the linear approximation on the k+1 points v_1,.., v_(k+1).
  86. This gives a linear program to solve (where the linear approximations
  87. of the constraint functions are constrained to be non-negative).
  88. However, the linear approximations are likely only good
  89. approximations near the current simplex, so the linear program is
  90. given the further requirement that the solution, which
  91. will become x_(j+1), must be within RHO_j from x_j. RHO_j only
  92. decreases, never increases. The initial RHO_j is rhobeg and the
  93. final RHO_j is rhoend. In this way COBYLA's iterations behave
  94. like a trust region algorithm.
  95. Additionally, the linear program may be inconsistent, or the
  96. approximation may give poor improvement. For details about
  97. how these issues are resolved, as well as how the points v_i are
  98. updated, refer to the source code or the references below.
  99. References
  100. ----------
  101. Powell M.J.D. (1994), "A direct search optimization method that models
  102. the objective and constraint functions by linear interpolation.", in
  103. Advances in Optimization and Numerical Analysis, eds. S. Gomez and
  104. J-P Hennart, Kluwer Academic (Dordrecht), pp. 51-67
  105. Powell M.J.D. (1998), "Direct search algorithms for optimization
  106. calculations", Acta Numerica 7, 287-336
  107. Powell M.J.D. (2007), "A view of algorithms for optimization without
  108. derivatives", Cambridge University Technical Report DAMTP 2007/NA03
  109. Examples
  110. --------
  111. Minimize the objective function f(x,y) = x*y subject
  112. to the constraints x**2 + y**2 < 1 and y > 0::
  113. >>> def objective(x):
  114. ... return x[0]*x[1]
  115. ...
  116. >>> def constr1(x):
  117. ... return 1 - (x[0]**2 + x[1]**2)
  118. ...
  119. >>> def constr2(x):
  120. ... return x[1]
  121. ...
  122. >>> from scipy.optimize import fmin_cobyla
  123. >>> fmin_cobyla(objective, [0.0, 0.1], [constr1, constr2], rhoend=1e-7)
  124. array([-0.70710685, 0.70710671])
  125. The exact solution is (-sqrt(2)/2, sqrt(2)/2).
  126. """
  127. err = "cons must be a sequence of callable functions or a single"\
  128. " callable function."
  129. try:
  130. len(cons)
  131. except TypeError as e:
  132. if callable(cons):
  133. cons = [cons]
  134. else:
  135. raise TypeError(err) from e
  136. else:
  137. for thisfunc in cons:
  138. if not callable(thisfunc):
  139. raise TypeError(err)
  140. if consargs is None:
  141. consargs = args
  142. # build constraints
  143. con = tuple({'type': 'ineq', 'fun': c, 'args': consargs} for c in cons)
  144. # options
  145. opts = {'rhobeg': rhobeg,
  146. 'tol': rhoend,
  147. 'disp': disp,
  148. 'maxiter': maxfun,
  149. 'catol': catol,
  150. 'callback': callback}
  151. sol = _minimize_cobyla(func, x0, args, constraints=con,
  152. **opts)
  153. if disp and not sol['success']:
  154. print("COBYLA failed to find a solution: %s" % (sol.message,))
  155. return sol['x']
  156. @synchronized
  157. def _minimize_cobyla(fun, x0, args=(), constraints=(),
  158. rhobeg=1.0, tol=1e-4, maxiter=1000,
  159. disp=False, catol=2e-4, callback=None,
  160. **unknown_options):
  161. """
  162. Minimize a scalar function of one or more variables using the
  163. Constrained Optimization BY Linear Approximation (COBYLA) algorithm.
  164. Options
  165. -------
  166. rhobeg : float
  167. Reasonable initial changes to the variables.
  168. tol : float
  169. Final accuracy in the optimization (not precisely guaranteed).
  170. This is a lower bound on the size of the trust region.
  171. disp : bool
  172. Set to True to print convergence messages. If False,
  173. `verbosity` is ignored as set to 0.
  174. maxiter : int
  175. Maximum number of function evaluations.
  176. catol : float
  177. Tolerance (absolute) for constraint violations
  178. """
  179. _check_unknown_options(unknown_options)
  180. maxfun = maxiter
  181. rhoend = tol
  182. iprint = int(bool(disp))
  183. # check constraints
  184. if isinstance(constraints, dict):
  185. constraints = (constraints, )
  186. for ic, con in enumerate(constraints):
  187. # check type
  188. try:
  189. ctype = con['type'].lower()
  190. except KeyError as e:
  191. raise KeyError('Constraint %d has no type defined.' % ic) from e
  192. except TypeError as e:
  193. raise TypeError('Constraints must be defined using a '
  194. 'dictionary.') from e
  195. except AttributeError as e:
  196. raise TypeError("Constraint's type must be a string.") from e
  197. else:
  198. if ctype != 'ineq':
  199. raise ValueError("Constraints of type '%s' not handled by "
  200. "COBYLA." % con['type'])
  201. # check function
  202. if 'fun' not in con:
  203. raise KeyError('Constraint %d has no function defined.' % ic)
  204. # check extra arguments
  205. if 'args' not in con:
  206. con['args'] = ()
  207. # m is the total number of constraint values
  208. # it takes into account that some constraints may be vector-valued
  209. cons_lengths = []
  210. for c in constraints:
  211. f = c['fun'](x0, *c['args'])
  212. try:
  213. cons_length = len(f)
  214. except TypeError:
  215. cons_length = 1
  216. cons_lengths.append(cons_length)
  217. m = sum(cons_lengths)
  218. def calcfc(x, con):
  219. f = fun(np.copy(x), *args)
  220. i = 0
  221. for size, c in izip(cons_lengths, constraints):
  222. con[i: i + size] = c['fun'](x, *c['args'])
  223. i += size
  224. return f
  225. def wrapped_callback(x):
  226. if callback is not None:
  227. callback(np.copy(x))
  228. info = np.zeros(4, np.float64)
  229. xopt, info = cobyla.minimize(calcfc, m=m, x=np.copy(x0), rhobeg=rhobeg,
  230. rhoend=rhoend, iprint=iprint, maxfun=maxfun,
  231. dinfo=info, callback=wrapped_callback)
  232. if info[3] > catol:
  233. # Check constraint violation
  234. info[0] = 4
  235. return OptimizeResult(x=xopt,
  236. status=int(info[0]),
  237. success=info[0] == 1,
  238. message={1: 'Optimization terminated successfully.',
  239. 2: 'Maximum number of function evaluations '
  240. 'has been exceeded.',
  241. 3: 'Rounding errors are becoming damaging '
  242. 'in COBYLA subroutine.',
  243. 4: 'Did not converge to a solution '
  244. 'satisfying the constraints. See '
  245. '`maxcv` for magnitude of violation.',
  246. 5: 'NaN result encountered.'
  247. }.get(info[0], 'Unknown exit status.'),
  248. nfev=int(info[1]),
  249. fun=info[2],
  250. maxcv=info[3])