Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_optimize.py: 6%
1413 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1#__docformat__ = "restructuredtext en"
2# ******NOTICE***************
3# optimize.py module by Travis E. Oliphant
4#
5# You may copy and use this module as you see fit with no
6# guarantee implied provided you keep this notice in all copies.
7# *****END NOTICE************
9# A collection of optimization algorithms. Version 0.5
10# CHANGES
11# Added fminbound (July 2001)
12# Added brute (Aug. 2002)
13# Finished line search satisfying strong Wolfe conditions (Mar. 2004)
14# Updated strong Wolfe conditions line search to use
15# cubic-interpolation (Mar. 2004)
18# Minimization routines
20__all__ = ['fmin', 'fmin_powell', 'fmin_bfgs', 'fmin_ncg', 'fmin_cg',
21 'fminbound', 'brent', 'golden', 'bracket', 'rosen', 'rosen_der',
22 'rosen_hess', 'rosen_hess_prod', 'brute', 'approx_fprime',
23 'line_search', 'check_grad', 'OptimizeResult', 'show_options',
24 'OptimizeWarning']
26__docformat__ = "restructuredtext en"
28import warnings
29import sys
30import inspect
31from numpy import (atleast_1d, eye, argmin, zeros, shape, squeeze,
32 asarray, sqrt, Inf)
33import numpy as np
34from scipy.sparse.linalg import LinearOperator
35from ._linesearch import (line_search_wolfe1, line_search_wolfe2,
36 line_search_wolfe2 as line_search,
37 LineSearchWarning)
38from ._numdiff import approx_derivative
39from ._hessian_update_strategy import HessianUpdateStrategy
40from scipy._lib._util import getfullargspec_no_self as _getfullargspec
41from scipy._lib._util import MapWrapper, check_random_state
42from scipy.optimize._differentiable_functions import ScalarFunction, FD_METHODS
45# standard status messages of optimizers
46_status_message = {'success': 'Optimization terminated successfully.',
47 'maxfev': 'Maximum number of function evaluations has '
48 'been exceeded.',
49 'maxiter': 'Maximum number of iterations has been '
50 'exceeded.',
51 'pr_loss': 'Desired error not necessarily achieved due '
52 'to precision loss.',
53 'nan': 'NaN result encountered.',
54 'out_of_bounds': 'The result is outside of the provided '
55 'bounds.'}
58class MemoizeJac:
59 """ Decorator that caches the return values of a function returning `(fun, grad)`
60 each time it is called. """
62 def __init__(self, fun):
63 self.fun = fun
64 self.jac = None
65 self._value = None
66 self.x = None
68 def _compute_if_needed(self, x, *args):
69 if not np.all(x == self.x) or self._value is None or self.jac is None:
70 self.x = np.asarray(x).copy()
71 fg = self.fun(x, *args)
72 self.jac = fg[1]
73 self._value = fg[0]
75 def __call__(self, x, *args):
76 """ returns the function value """
77 self._compute_if_needed(x, *args)
78 return self._value
80 def derivative(self, x, *args):
81 self._compute_if_needed(x, *args)
82 return self.jac
85def _indenter(s, n=0):
86 """
87 Ensures that lines after the first are indented by the specified amount
88 """
89 split = s.split("\n")
90 indent = " "*n
91 return ("\n" + indent).join(split)
94def _float_formatter_10(x):
95 """
96 Returns a string representation of a float with exactly ten characters
97 """
98 if np.isposinf(x):
99 return " inf"
100 elif np.isneginf(x):
101 return " -inf"
102 elif np.isnan(x):
103 return " nan"
104 return np.format_float_scientific(x, precision=3, pad_left=2, unique=False)
107def _dict_formatter(d, n=0, mplus=1, sorter=None):
108 """
109 Pretty printer for dictionaries
111 `n` keeps track of the starting indentation;
112 lines are indented by this much after a line break.
113 `mplus` is additional left padding applied to keys
114 """
115 if isinstance(d, dict):
116 m = max(map(len, list(d.keys()))) + mplus # width to print keys
117 s = '\n'.join([k.rjust(m) + ': ' + # right justified, width m
118 _indenter(_dict_formatter(v, m+n+2, 0, sorter), m+2)
119 for k, v in sorter(d)]) # +2 for ': '
120 else:
121 # By default, NumPy arrays print with linewidth=76. `n` is
122 # the indent at which a line begins printing, so it is subtracted
123 # from the default to avoid exceeding 76 characters total.
124 # `edgeitems` is the number of elements to include before and after
125 # ellipses when arrays are not shown in full.
126 # `threshold` is the maximum number of elements for which an
127 # array is shown in full.
128 # These values tend to work well for use with OptimizeResult.
129 with np.printoptions(linewidth=76-n, edgeitems=2, threshold=12,
130 formatter={'float_kind': _float_formatter_10}):
131 s = str(d)
132 return s
135def _wrap_callback(callback, method=None):
136 """Wrap a user-provided callback so that attributes can be attached."""
137 if callback is None or method in {'tnc', 'slsqp', 'cobyla'}:
138 return callback # don't wrap
140 sig = inspect.signature(callback)
142 if set(sig.parameters) == {'intermediate_result'}:
143 def wrapped_callback(res):
144 return callback(intermediate_result=res)
145 elif method == 'trust-constr':
146 def wrapped_callback(res):
147 return callback(np.copy(res.x), res)
148 else:
149 def wrapped_callback(res):
150 return callback(np.copy(res.x))
152 wrapped_callback.stop_iteration = False
153 return wrapped_callback
156def _call_callback_maybe_halt(callback, res):
157 """Call wrapped callback; return True if minimization should stop.
159 Parameters
160 ----------
161 callback : callable or None
162 A user-provided callback wrapped with `_wrap_callback`
163 res : OptimizeResult
164 Information about the current iterate
166 Returns
167 -------
168 halt : bool
169 True if minimization should stop
171 """
172 if callback is None:
173 return False
174 try:
175 callback(res)
176 return False
177 except StopIteration:
178 callback.stop_iteration = True # make `minimize` override status/msg
179 return True
182class OptimizeResult(dict):
183 """ Represents the optimization result.
185 Attributes
186 ----------
187 x : ndarray
188 The solution of the optimization.
189 success : bool
190 Whether or not the optimizer exited successfully.
191 status : int
192 Termination status of the optimizer. Its value depends on the
193 underlying solver. Refer to `message` for details.
194 message : str
195 Description of the cause of the termination.
196 fun, jac, hess: ndarray
197 Values of objective function, its Jacobian and its Hessian (if
198 available). The Hessians may be approximations, see the documentation
199 of the function in question.
200 hess_inv : object
201 Inverse of the objective function's Hessian; may be an approximation.
202 Not available for all solvers. The type of this attribute may be
203 either np.ndarray or scipy.sparse.linalg.LinearOperator.
204 nfev, njev, nhev : int
205 Number of evaluations of the objective functions and of its
206 Jacobian and Hessian.
207 nit : int
208 Number of iterations performed by the optimizer.
209 maxcv : float
210 The maximum constraint violation.
212 Notes
213 -----
214 Depending on the specific solver being used, `OptimizeResult` may
215 not have all attributes listed here, and they may have additional
216 attributes not listed here. Since this class is essentially a
217 subclass of dict with attribute accessors, one can see which
218 attributes are available using the `OptimizeResult.keys` method.
219 """
221 def __getattr__(self, name):
222 try:
223 return self[name]
224 except KeyError as e:
225 raise AttributeError(name) from e
227 __setattr__ = dict.__setitem__
228 __delattr__ = dict.__delitem__
230 def __repr__(self):
231 order_keys = ['message', 'success', 'status', 'fun', 'funl', 'x', 'xl',
232 'col_ind', 'nit', 'lower', 'upper', 'eqlin', 'ineqlin',
233 'converged', 'flag', 'function_calls', 'iterations',
234 'root']
235 # 'slack', 'con' are redundant with residuals
236 # 'crossover_nit' is probably not interesting to most users
237 omit_keys = {'slack', 'con', 'crossover_nit'}
239 def key(item):
240 try:
241 return order_keys.index(item[0].lower())
242 except ValueError: # item not in list
243 return np.inf
245 def omit_redundant(items):
246 for item in items:
247 if item[0] in omit_keys:
248 continue
249 yield item
251 def item_sorter(d):
252 return sorted(omit_redundant(d.items()), key=key)
254 if self.keys():
255 return _dict_formatter(self, sorter=item_sorter)
256 else:
257 return self.__class__.__name__ + "()"
259 def __dir__(self):
260 return list(self.keys())
263class OptimizeWarning(UserWarning):
264 pass
267def _check_unknown_options(unknown_options):
268 if unknown_options:
269 msg = ", ".join(map(str, unknown_options.keys()))
270 # Stack level 4: this is called from _minimize_*, which is
271 # called from another function in SciPy. Level 4 is the first
272 # level in user code.
273 warnings.warn("Unknown solver options: %s" % msg, OptimizeWarning, 4)
276def is_finite_scalar(x):
277 """Test whether `x` is either a finite scalar or a finite array scalar.
279 """
280 return np.size(x) == 1 and np.isfinite(x)
283_epsilon = sqrt(np.finfo(float).eps)
286def vecnorm(x, ord=2):
287 if ord == Inf:
288 return np.amax(np.abs(x))
289 elif ord == -Inf:
290 return np.amin(np.abs(x))
291 else:
292 return np.sum(np.abs(x)**ord, axis=0)**(1.0 / ord)
295def _prepare_scalar_function(fun, x0, jac=None, args=(), bounds=None,
296 epsilon=None, finite_diff_rel_step=None,
297 hess=None):
298 """
299 Creates a ScalarFunction object for use with scalar minimizers
300 (BFGS/LBFGSB/SLSQP/TNC/CG/etc).
302 Parameters
303 ----------
304 fun : callable
305 The objective function to be minimized.
307 ``fun(x, *args) -> float``
309 where ``x`` is an 1-D array with shape (n,) and ``args``
310 is a tuple of the fixed parameters needed to completely
311 specify the function.
312 x0 : ndarray, shape (n,)
313 Initial guess. Array of real elements of size (n,),
314 where 'n' is the number of independent variables.
315 jac : {callable, '2-point', '3-point', 'cs', None}, optional
316 Method for computing the gradient vector. If it is a callable, it
317 should be a function that returns the gradient vector:
319 ``jac(x, *args) -> array_like, shape (n,)``
321 If one of `{'2-point', '3-point', 'cs'}` is selected then the gradient
322 is calculated with a relative step for finite differences. If `None`,
323 then two-point finite differences with an absolute step is used.
324 args : tuple, optional
325 Extra arguments passed to the objective function and its
326 derivatives (`fun`, `jac` functions).
327 bounds : sequence, optional
328 Bounds on variables. 'new-style' bounds are required.
329 eps : float or ndarray
330 If `jac is None` the absolute step size used for numerical
331 approximation of the jacobian via forward differences.
332 finite_diff_rel_step : None or array_like, optional
333 If `jac in ['2-point', '3-point', 'cs']` the relative step size to
334 use for numerical approximation of the jacobian. The absolute step
335 size is computed as ``h = rel_step * sign(x0) * max(1, abs(x0))``,
336 possibly adjusted to fit into the bounds. For ``method='3-point'``
337 the sign of `h` is ignored. If None (default) then step is selected
338 automatically.
339 hess : {callable, '2-point', '3-point', 'cs', None}
340 Computes the Hessian matrix. If it is callable, it should return the
341 Hessian matrix:
343 ``hess(x, *args) -> {LinearOperator, spmatrix, array}, (n, n)``
345 Alternatively, the keywords {'2-point', '3-point', 'cs'} select a
346 finite difference scheme for numerical estimation.
347 Whenever the gradient is estimated via finite-differences, the Hessian
348 cannot be estimated with options {'2-point', '3-point', 'cs'} and needs
349 to be estimated using one of the quasi-Newton strategies.
351 Returns
352 -------
353 sf : ScalarFunction
354 """
355 if callable(jac):
356 grad = jac
357 elif jac in FD_METHODS:
358 # epsilon is set to None so that ScalarFunction is made to use
359 # rel_step
360 epsilon = None
361 grad = jac
362 else:
363 # default (jac is None) is to do 2-point finite differences with
364 # absolute step size. ScalarFunction has to be provided an
365 # epsilon value that is not None to use absolute steps. This is
366 # normally the case from most _minimize* methods.
367 grad = '2-point'
368 epsilon = epsilon
370 if hess is None:
371 # ScalarFunction requires something for hess, so we give a dummy
372 # implementation here if nothing is provided, return a value of None
373 # so that downstream minimisers halt. The results of `fun.hess`
374 # should not be used.
375 def hess(x, *args):
376 return None
378 if bounds is None:
379 bounds = (-np.inf, np.inf)
381 # ScalarFunction caches. Reuse of fun(x) during grad
382 # calculation reduces overall function evaluations.
383 sf = ScalarFunction(fun, x0, args, grad, hess,
384 finite_diff_rel_step, bounds, epsilon=epsilon)
386 return sf
389def _clip_x_for_func(func, bounds):
390 # ensures that x values sent to func are clipped to bounds
392 # this is used as a mitigation for gh11403, slsqp/tnc sometimes
393 # suggest a move that is outside the limits by 1 or 2 ULP. This
394 # unclean fix makes sure x is strictly within bounds.
395 def eval(x):
396 x = _check_clip_x(x, bounds)
397 return func(x)
399 return eval
402def _check_clip_x(x, bounds):
403 if (x < bounds[0]).any() or (x > bounds[1]).any():
404 warnings.warn("Values in x were outside bounds during a "
405 "minimize step, clipping to bounds", RuntimeWarning)
406 x = np.clip(x, bounds[0], bounds[1])
407 return x
409 return x
412def rosen(x):
413 """
414 The Rosenbrock function.
416 The function computed is::
418 sum(100.0*(x[1:] - x[:-1]**2.0)**2.0 + (1 - x[:-1])**2.0)
420 Parameters
421 ----------
422 x : array_like
423 1-D array of points at which the Rosenbrock function is to be computed.
425 Returns
426 -------
427 f : float
428 The value of the Rosenbrock function.
430 See Also
431 --------
432 rosen_der, rosen_hess, rosen_hess_prod
434 Examples
435 --------
436 >>> import numpy as np
437 >>> from scipy.optimize import rosen
438 >>> X = 0.1 * np.arange(10)
439 >>> rosen(X)
440 76.56
442 For higher-dimensional input ``rosen`` broadcasts.
443 In the following example, we use this to plot a 2D landscape.
444 Note that ``rosen_hess`` does not broadcast in this manner.
446 >>> import matplotlib.pyplot as plt
447 >>> from mpl_toolkits.mplot3d import Axes3D
448 >>> x = np.linspace(-1, 1, 50)
449 >>> X, Y = np.meshgrid(x, x)
450 >>> ax = plt.subplot(111, projection='3d')
451 >>> ax.plot_surface(X, Y, rosen([X, Y]))
452 >>> plt.show()
453 """
454 x = asarray(x)
455 r = np.sum(100.0 * (x[1:] - x[:-1]**2.0)**2.0 + (1 - x[:-1])**2.0,
456 axis=0)
457 return r
460def rosen_der(x):
461 """
462 The derivative (i.e. gradient) of the Rosenbrock function.
464 Parameters
465 ----------
466 x : array_like
467 1-D array of points at which the derivative is to be computed.
469 Returns
470 -------
471 rosen_der : (N,) ndarray
472 The gradient of the Rosenbrock function at `x`.
474 See Also
475 --------
476 rosen, rosen_hess, rosen_hess_prod
478 Examples
479 --------
480 >>> import numpy as np
481 >>> from scipy.optimize import rosen_der
482 >>> X = 0.1 * np.arange(9)
483 >>> rosen_der(X)
484 array([ -2. , 10.6, 15.6, 13.4, 6.4, -3. , -12.4, -19.4, 62. ])
486 """
487 x = asarray(x)
488 xm = x[1:-1]
489 xm_m1 = x[:-2]
490 xm_p1 = x[2:]
491 der = np.zeros_like(x)
492 der[1:-1] = (200 * (xm - xm_m1**2) -
493 400 * (xm_p1 - xm**2) * xm - 2 * (1 - xm))
494 der[0] = -400 * x[0] * (x[1] - x[0]**2) - 2 * (1 - x[0])
495 der[-1] = 200 * (x[-1] - x[-2]**2)
496 return der
499def rosen_hess(x):
500 """
501 The Hessian matrix of the Rosenbrock function.
503 Parameters
504 ----------
505 x : array_like
506 1-D array of points at which the Hessian matrix is to be computed.
508 Returns
509 -------
510 rosen_hess : ndarray
511 The Hessian matrix of the Rosenbrock function at `x`.
513 See Also
514 --------
515 rosen, rosen_der, rosen_hess_prod
517 Examples
518 --------
519 >>> import numpy as np
520 >>> from scipy.optimize import rosen_hess
521 >>> X = 0.1 * np.arange(4)
522 >>> rosen_hess(X)
523 array([[-38., 0., 0., 0.],
524 [ 0., 134., -40., 0.],
525 [ 0., -40., 130., -80.],
526 [ 0., 0., -80., 200.]])
528 """
529 x = atleast_1d(x)
530 H = np.diag(-400 * x[:-1], 1) - np.diag(400 * x[:-1], -1)
531 diagonal = np.zeros(len(x), dtype=x.dtype)
532 diagonal[0] = 1200 * x[0]**2 - 400 * x[1] + 2
533 diagonal[-1] = 200
534 diagonal[1:-1] = 202 + 1200 * x[1:-1]**2 - 400 * x[2:]
535 H = H + np.diag(diagonal)
536 return H
539def rosen_hess_prod(x, p):
540 """
541 Product of the Hessian matrix of the Rosenbrock function with a vector.
543 Parameters
544 ----------
545 x : array_like
546 1-D array of points at which the Hessian matrix is to be computed.
547 p : array_like
548 1-D array, the vector to be multiplied by the Hessian matrix.
550 Returns
551 -------
552 rosen_hess_prod : ndarray
553 The Hessian matrix of the Rosenbrock function at `x` multiplied
554 by the vector `p`.
556 See Also
557 --------
558 rosen, rosen_der, rosen_hess
560 Examples
561 --------
562 >>> import numpy as np
563 >>> from scipy.optimize import rosen_hess_prod
564 >>> X = 0.1 * np.arange(9)
565 >>> p = 0.5 * np.arange(9)
566 >>> rosen_hess_prod(X, p)
567 array([ -0., 27., -10., -95., -192., -265., -278., -195., -180.])
569 """
570 x = atleast_1d(x)
571 Hp = np.zeros(len(x), dtype=x.dtype)
572 Hp[0] = (1200 * x[0]**2 - 400 * x[1] + 2) * p[0] - 400 * x[0] * p[1]
573 Hp[1:-1] = (-400 * x[:-2] * p[:-2] +
574 (202 + 1200 * x[1:-1]**2 - 400 * x[2:]) * p[1:-1] -
575 400 * x[1:-1] * p[2:])
576 Hp[-1] = -400 * x[-2] * p[-2] + 200*p[-1]
577 return Hp
580def _wrap_scalar_function(function, args):
581 # wraps a minimizer function to count number of evaluations
582 # and to easily provide an args kwd.
583 ncalls = [0]
584 if function is None:
585 return ncalls, None
587 def function_wrapper(x, *wrapper_args):
588 ncalls[0] += 1
589 # A copy of x is sent to the user function (gh13740)
590 fx = function(np.copy(x), *(wrapper_args + args))
591 # Ideally, we'd like to a have a true scalar returned from f(x). For
592 # backwards-compatibility, also allow np.array([1.3]), np.array([[1.3]]) etc.
593 if not np.isscalar(fx):
594 try:
595 fx = np.asarray(fx).item()
596 except (TypeError, ValueError) as e:
597 raise ValueError("The user-provided objective function "
598 "must return a scalar value.") from e
599 return fx
601 return ncalls, function_wrapper
604class _MaxFuncCallError(RuntimeError):
605 pass
608def _wrap_scalar_function_maxfun_validation(function, args, maxfun):
609 # wraps a minimizer function to count number of evaluations
610 # and to easily provide an args kwd.
611 ncalls = [0]
612 if function is None:
613 return ncalls, None
615 def function_wrapper(x, *wrapper_args):
616 if ncalls[0] >= maxfun:
617 raise _MaxFuncCallError("Too many function calls")
618 ncalls[0] += 1
619 # A copy of x is sent to the user function (gh13740)
620 fx = function(np.copy(x), *(wrapper_args + args))
621 # Ideally, we'd like to a have a true scalar returned from f(x). For
622 # backwards-compatibility, also allow np.array([1.3]),
623 # np.array([[1.3]]) etc.
624 if not np.isscalar(fx):
625 try:
626 fx = np.asarray(fx).item()
627 except (TypeError, ValueError) as e:
628 raise ValueError("The user-provided objective function "
629 "must return a scalar value.") from e
630 return fx
632 return ncalls, function_wrapper
635def fmin(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, maxfun=None,
636 full_output=0, disp=1, retall=0, callback=None, initial_simplex=None):
637 """
638 Minimize a function using the downhill simplex algorithm.
640 This algorithm only uses function values, not derivatives or second
641 derivatives.
643 Parameters
644 ----------
645 func : callable func(x,*args)
646 The objective function to be minimized.
647 x0 : ndarray
648 Initial guess.
649 args : tuple, optional
650 Extra arguments passed to func, i.e., ``f(x,*args)``.
651 xtol : float, optional
652 Absolute error in xopt between iterations that is acceptable for
653 convergence.
654 ftol : number, optional
655 Absolute error in func(xopt) between iterations that is acceptable for
656 convergence.
657 maxiter : int, optional
658 Maximum number of iterations to perform.
659 maxfun : number, optional
660 Maximum number of function evaluations to make.
661 full_output : bool, optional
662 Set to True if fopt and warnflag outputs are desired.
663 disp : bool, optional
664 Set to True to print convergence messages.
665 retall : bool, optional
666 Set to True to return list of solutions at each iteration.
667 callback : callable, optional
668 Called after each iteration, as callback(xk), where xk is the
669 current parameter vector.
670 initial_simplex : array_like of shape (N + 1, N), optional
671 Initial simplex. If given, overrides `x0`.
672 ``initial_simplex[j,:]`` should contain the coordinates of
673 the jth vertex of the ``N+1`` vertices in the simplex, where
674 ``N`` is the dimension.
676 Returns
677 -------
678 xopt : ndarray
679 Parameter that minimizes function.
680 fopt : float
681 Value of function at minimum: ``fopt = func(xopt)``.
682 iter : int
683 Number of iterations performed.
684 funcalls : int
685 Number of function calls made.
686 warnflag : int
687 1 : Maximum number of function evaluations made.
688 2 : Maximum number of iterations reached.
689 allvecs : list
690 Solution at each iteration.
692 See also
693 --------
694 minimize: Interface to minimization algorithms for multivariate
695 functions. See the 'Nelder-Mead' `method` in particular.
697 Notes
698 -----
699 Uses a Nelder-Mead simplex algorithm to find the minimum of function of
700 one or more variables.
702 This algorithm has a long history of successful use in applications.
703 But it will usually be slower than an algorithm that uses first or
704 second derivative information. In practice, it can have poor
705 performance in high-dimensional problems and is not robust to
706 minimizing complicated functions. Additionally, there currently is no
707 complete theory describing when the algorithm will successfully
708 converge to the minimum, or how fast it will if it does. Both the ftol and
709 xtol criteria must be met for convergence.
711 Examples
712 --------
713 >>> def f(x):
714 ... return x**2
716 >>> from scipy import optimize
718 >>> minimum = optimize.fmin(f, 1)
719 Optimization terminated successfully.
720 Current function value: 0.000000
721 Iterations: 17
722 Function evaluations: 34
723 >>> minimum[0]
724 -8.8817841970012523e-16
726 References
727 ----------
728 .. [1] Nelder, J.A. and Mead, R. (1965), "A simplex method for function
729 minimization", The Computer Journal, 7, pp. 308-313
731 .. [2] Wright, M.H. (1996), "Direct Search Methods: Once Scorned, Now
732 Respectable", in Numerical Analysis 1995, Proceedings of the
733 1995 Dundee Biennial Conference in Numerical Analysis, D.F.
734 Griffiths and G.A. Watson (Eds.), Addison Wesley Longman,
735 Harlow, UK, pp. 191-208.
737 """
738 opts = {'xatol': xtol,
739 'fatol': ftol,
740 'maxiter': maxiter,
741 'maxfev': maxfun,
742 'disp': disp,
743 'return_all': retall,
744 'initial_simplex': initial_simplex}
746 callback = _wrap_callback(callback)
747 res = _minimize_neldermead(func, x0, args, callback=callback, **opts)
748 if full_output:
749 retlist = res['x'], res['fun'], res['nit'], res['nfev'], res['status']
750 if retall:
751 retlist += (res['allvecs'], )
752 return retlist
753 else:
754 if retall:
755 return res['x'], res['allvecs']
756 else:
757 return res['x']
760def _minimize_neldermead(func, x0, args=(), callback=None,
761 maxiter=None, maxfev=None, disp=False,
762 return_all=False, initial_simplex=None,
763 xatol=1e-4, fatol=1e-4, adaptive=False, bounds=None,
764 **unknown_options):
765 """
766 Minimization of scalar function of one or more variables using the
767 Nelder-Mead algorithm.
769 Options
770 -------
771 disp : bool
772 Set to True to print convergence messages.
773 maxiter, maxfev : int
774 Maximum allowed number of iterations and function evaluations.
775 Will default to ``N*200``, where ``N`` is the number of
776 variables, if neither `maxiter` or `maxfev` is set. If both
777 `maxiter` and `maxfev` are set, minimization will stop at the
778 first reached.
779 return_all : bool, optional
780 Set to True to return a list of the best solution at each of the
781 iterations.
782 initial_simplex : array_like of shape (N + 1, N)
783 Initial simplex. If given, overrides `x0`.
784 ``initial_simplex[j,:]`` should contain the coordinates of
785 the jth vertex of the ``N+1`` vertices in the simplex, where
786 ``N`` is the dimension.
787 xatol : float, optional
788 Absolute error in xopt between iterations that is acceptable for
789 convergence.
790 fatol : number, optional
791 Absolute error in func(xopt) between iterations that is acceptable for
792 convergence.
793 adaptive : bool, optional
794 Adapt algorithm parameters to dimensionality of problem. Useful for
795 high-dimensional minimization [1]_.
796 bounds : sequence or `Bounds`, optional
797 Bounds on variables. There are two ways to specify the bounds:
799 1. Instance of `Bounds` class.
800 2. Sequence of ``(min, max)`` pairs for each element in `x`. None
801 is used to specify no bound.
803 Note that this just clips all vertices in simplex based on
804 the bounds.
806 References
807 ----------
808 .. [1] Gao, F. and Han, L.
809 Implementing the Nelder-Mead simplex algorithm with adaptive
810 parameters. 2012. Computational Optimization and Applications.
811 51:1, pp. 259-277
813 """
814 _check_unknown_options(unknown_options)
815 maxfun = maxfev
816 retall = return_all
818 x0 = np.atleast_1d(x0).flatten()
819 x0 = np.asfarray(x0, x0.dtype)
821 if adaptive:
822 dim = float(len(x0))
823 rho = 1
824 chi = 1 + 2/dim
825 psi = 0.75 - 1/(2*dim)
826 sigma = 1 - 1/dim
827 else:
828 rho = 1
829 chi = 2
830 psi = 0.5
831 sigma = 0.5
833 nonzdelt = 0.05
834 zdelt = 0.00025
836 if bounds is not None:
837 lower_bound, upper_bound = bounds.lb, bounds.ub
838 # check bounds
839 if (lower_bound > upper_bound).any():
840 raise ValueError("Nelder Mead - one of the lower bounds is greater than an upper bound.")
841 if np.any(lower_bound > x0) or np.any(x0 > upper_bound):
842 warnings.warn("Initial guess is not within the specified bounds",
843 OptimizeWarning, 3)
845 if bounds is not None:
846 x0 = np.clip(x0, lower_bound, upper_bound)
848 if initial_simplex is None:
849 N = len(x0)
851 sim = np.empty((N + 1, N), dtype=x0.dtype)
852 sim[0] = x0
853 for k in range(N):
854 y = np.array(x0, copy=True)
855 if y[k] != 0:
856 y[k] = (1 + nonzdelt)*y[k]
857 else:
858 y[k] = zdelt
859 sim[k + 1] = y
860 else:
861 sim = np.atleast_2d(initial_simplex).copy()
862 sim = np.asfarray(sim, sim.dtype)
863 if sim.ndim != 2 or sim.shape[0] != sim.shape[1] + 1:
864 raise ValueError("`initial_simplex` should be an array of shape (N+1,N)")
865 if len(x0) != sim.shape[1]:
866 raise ValueError("Size of `initial_simplex` is not consistent with `x0`")
867 N = sim.shape[1]
869 if retall:
870 allvecs = [sim[0]]
872 # If neither are set, then set both to default
873 if maxiter is None and maxfun is None:
874 maxiter = N * 200
875 maxfun = N * 200
876 elif maxiter is None:
877 # Convert remaining Nones, to np.inf, unless the other is np.inf, in
878 # which case use the default to avoid unbounded iteration
879 if maxfun == np.inf:
880 maxiter = N * 200
881 else:
882 maxiter = np.inf
883 elif maxfun is None:
884 if maxiter == np.inf:
885 maxfun = N * 200
886 else:
887 maxfun = np.inf
889 if bounds is not None:
890 sim = np.clip(sim, lower_bound, upper_bound)
892 one2np1 = list(range(1, N + 1))
893 fsim = np.full((N + 1,), np.inf, dtype=float)
895 fcalls, func = _wrap_scalar_function_maxfun_validation(func, args, maxfun)
897 try:
898 for k in range(N + 1):
899 fsim[k] = func(sim[k])
900 except _MaxFuncCallError:
901 pass
902 finally:
903 ind = np.argsort(fsim)
904 sim = np.take(sim, ind, 0)
905 fsim = np.take(fsim, ind, 0)
907 ind = np.argsort(fsim)
908 fsim = np.take(fsim, ind, 0)
909 # sort so sim[0,:] has the lowest function value
910 sim = np.take(sim, ind, 0)
912 iterations = 1
914 while (fcalls[0] < maxfun and iterations < maxiter):
915 try:
916 if (np.max(np.ravel(np.abs(sim[1:] - sim[0]))) <= xatol and
917 np.max(np.abs(fsim[0] - fsim[1:])) <= fatol):
918 break
920 xbar = np.add.reduce(sim[:-1], 0) / N
921 xr = (1 + rho) * xbar - rho * sim[-1]
922 if bounds is not None:
923 xr = np.clip(xr, lower_bound, upper_bound)
924 fxr = func(xr)
925 doshrink = 0
927 if fxr < fsim[0]:
928 xe = (1 + rho * chi) * xbar - rho * chi * sim[-1]
929 if bounds is not None:
930 xe = np.clip(xe, lower_bound, upper_bound)
931 fxe = func(xe)
933 if fxe < fxr:
934 sim[-1] = xe
935 fsim[-1] = fxe
936 else:
937 sim[-1] = xr
938 fsim[-1] = fxr
939 else: # fsim[0] <= fxr
940 if fxr < fsim[-2]:
941 sim[-1] = xr
942 fsim[-1] = fxr
943 else: # fxr >= fsim[-2]
944 # Perform contraction
945 if fxr < fsim[-1]:
946 xc = (1 + psi * rho) * xbar - psi * rho * sim[-1]
947 if bounds is not None:
948 xc = np.clip(xc, lower_bound, upper_bound)
949 fxc = func(xc)
951 if fxc <= fxr:
952 sim[-1] = xc
953 fsim[-1] = fxc
954 else:
955 doshrink = 1
956 else:
957 # Perform an inside contraction
958 xcc = (1 - psi) * xbar + psi * sim[-1]
959 if bounds is not None:
960 xcc = np.clip(xcc, lower_bound, upper_bound)
961 fxcc = func(xcc)
963 if fxcc < fsim[-1]:
964 sim[-1] = xcc
965 fsim[-1] = fxcc
966 else:
967 doshrink = 1
969 if doshrink:
970 for j in one2np1:
971 sim[j] = sim[0] + sigma * (sim[j] - sim[0])
972 if bounds is not None:
973 sim[j] = np.clip(
974 sim[j], lower_bound, upper_bound)
975 fsim[j] = func(sim[j])
976 iterations += 1
977 except _MaxFuncCallError:
978 pass
979 finally:
980 ind = np.argsort(fsim)
981 sim = np.take(sim, ind, 0)
982 fsim = np.take(fsim, ind, 0)
983 if retall:
984 allvecs.append(sim[0])
985 intermediate_result = OptimizeResult(x=sim[0], fun=fsim[0])
986 if _call_callback_maybe_halt(callback, intermediate_result):
987 break
989 x = sim[0]
990 fval = np.min(fsim)
991 warnflag = 0
993 if fcalls[0] >= maxfun:
994 warnflag = 1
995 msg = _status_message['maxfev']
996 if disp:
997 warnings.warn(msg, RuntimeWarning, 3)
998 elif iterations >= maxiter:
999 warnflag = 2
1000 msg = _status_message['maxiter']
1001 if disp:
1002 warnings.warn(msg, RuntimeWarning, 3)
1003 else:
1004 msg = _status_message['success']
1005 if disp:
1006 print(msg)
1007 print(" Current function value: %f" % fval)
1008 print(" Iterations: %d" % iterations)
1009 print(" Function evaluations: %d" % fcalls[0])
1011 result = OptimizeResult(fun=fval, nit=iterations, nfev=fcalls[0],
1012 status=warnflag, success=(warnflag == 0),
1013 message=msg, x=x, final_simplex=(sim, fsim))
1014 if retall:
1015 result['allvecs'] = allvecs
1016 return result
1019def approx_fprime(xk, f, epsilon=_epsilon, *args):
1020 """Finite difference approximation of the derivatives of a
1021 scalar or vector-valued function.
1023 If a function maps from :math:`R^n` to :math:`R^m`, its derivatives form
1024 an m-by-n matrix
1025 called the Jacobian, where an element :math:`(i, j)` is a partial
1026 derivative of f[i] with respect to ``xk[j]``.
1028 Parameters
1029 ----------
1030 xk : array_like
1031 The coordinate vector at which to determine the gradient of `f`.
1032 f : callable
1033 Function of which to estimate the derivatives of. Has the signature
1034 ``f(xk, *args)`` where `xk` is the argument in the form of a 1-D array
1035 and `args` is a tuple of any additional fixed parameters needed to
1036 completely specify the function. The argument `xk` passed to this
1037 function is an ndarray of shape (n,) (never a scalar even if n=1).
1038 It must return a 1-D array_like of shape (m,) or a scalar.
1040 .. versionchanged:: 1.9.0
1041 `f` is now able to return a 1-D array-like, with the :math:`(m, n)`
1042 Jacobian being estimated.
1044 epsilon : {float, array_like}, optional
1045 Increment to `xk` to use for determining the function gradient.
1046 If a scalar, uses the same finite difference delta for all partial
1047 derivatives. If an array, should contain one value per element of
1048 `xk`. Defaults to ``sqrt(np.finfo(float).eps)``, which is approximately
1049 1.49e-08.
1050 \\*args : args, optional
1051 Any other arguments that are to be passed to `f`.
1053 Returns
1054 -------
1055 jac : ndarray
1056 The partial derivatives of `f` to `xk`.
1058 See Also
1059 --------
1060 check_grad : Check correctness of gradient function against approx_fprime.
1062 Notes
1063 -----
1064 The function gradient is determined by the forward finite difference
1065 formula::
1067 f(xk[i] + epsilon[i]) - f(xk[i])
1068 f'[i] = ---------------------------------
1069 epsilon[i]
1071 Examples
1072 --------
1073 >>> import numpy as np
1074 >>> from scipy import optimize
1075 >>> def func(x, c0, c1):
1076 ... "Coordinate vector `x` should be an array of size two."
1077 ... return c0 * x[0]**2 + c1*x[1]**2
1079 >>> x = np.ones(2)
1080 >>> c0, c1 = (1, 200)
1081 >>> eps = np.sqrt(np.finfo(float).eps)
1082 >>> optimize.approx_fprime(x, func, [eps, np.sqrt(200) * eps], c0, c1)
1083 array([ 2. , 400.00004198])
1085 """
1086 xk = np.asarray(xk, float)
1087 f0 = f(xk, *args)
1089 return approx_derivative(f, xk, method='2-point', abs_step=epsilon,
1090 args=args, f0=f0)
1093def check_grad(func, grad, x0, *args, epsilon=_epsilon,
1094 direction='all', seed=None):
1095 """Check the correctness of a gradient function by comparing it against a
1096 (forward) finite-difference approximation of the gradient.
1098 Parameters
1099 ----------
1100 func : callable ``func(x0, *args)``
1101 Function whose derivative is to be checked.
1102 grad : callable ``grad(x0, *args)``
1103 Jacobian of `func`.
1104 x0 : ndarray
1105 Points to check `grad` against forward difference approximation of grad
1106 using `func`.
1107 args : \\*args, optional
1108 Extra arguments passed to `func` and `grad`.
1109 epsilon : float, optional
1110 Step size used for the finite difference approximation. It defaults to
1111 ``sqrt(np.finfo(float).eps)``, which is approximately 1.49e-08.
1112 direction : str, optional
1113 If set to ``'random'``, then gradients along a random vector
1114 are used to check `grad` against forward difference approximation
1115 using `func`. By default it is ``'all'``, in which case, all
1116 the one hot direction vectors are considered to check `grad`.
1117 If `func` is a vector valued function then only ``'all'`` can be used.
1118 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
1119 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
1120 singleton is used.
1121 If `seed` is an int, a new ``RandomState`` instance is used,
1122 seeded with `seed`.
1123 If `seed` is already a ``Generator`` or ``RandomState`` instance then
1124 that instance is used.
1125 Specify `seed` for reproducing the return value from this function.
1126 The random numbers generated with this seed affect the random vector
1127 along which gradients are computed to check ``grad``. Note that `seed`
1128 is only used when `direction` argument is set to `'random'`.
1130 Returns
1131 -------
1132 err : float
1133 The square root of the sum of squares (i.e., the 2-norm) of the
1134 difference between ``grad(x0, *args)`` and the finite difference
1135 approximation of `grad` using func at the points `x0`.
1137 See Also
1138 --------
1139 approx_fprime
1141 Examples
1142 --------
1143 >>> import numpy as np
1144 >>> def func(x):
1145 ... return x[0]**2 - 0.5 * x[1]**3
1146 >>> def grad(x):
1147 ... return [2 * x[0], -1.5 * x[1]**2]
1148 >>> from scipy.optimize import check_grad
1149 >>> check_grad(func, grad, [1.5, -1.5])
1150 2.9802322387695312e-08 # may vary
1151 >>> rng = np.random.default_rng()
1152 >>> check_grad(func, grad, [1.5, -1.5],
1153 ... direction='random', seed=rng)
1154 2.9802322387695312e-08
1156 """
1157 step = epsilon
1158 x0 = np.asarray(x0)
1160 def g(w, func, x0, v, *args):
1161 return func(x0 + w*v, *args)
1163 if direction == 'random':
1164 _grad = np.asanyarray(grad(x0, *args))
1165 if _grad.ndim > 1:
1166 raise ValueError("'random' can only be used with scalar valued"
1167 " func")
1168 random_state = check_random_state(seed)
1169 v = random_state.normal(0, 1, size=(x0.shape))
1170 _args = (func, x0, v) + args
1171 _func = g
1172 vars = np.zeros((1,))
1173 analytical_grad = np.dot(_grad, v)
1174 elif direction == 'all':
1175 _args = args
1176 _func = func
1177 vars = x0
1178 analytical_grad = grad(x0, *args)
1179 else:
1180 raise ValueError("{} is not a valid string for "
1181 "``direction`` argument".format(direction))
1183 return np.sqrt(np.sum(np.abs(
1184 (analytical_grad - approx_fprime(vars, _func, step, *_args))**2
1185 )))
1188def approx_fhess_p(x0, p, fprime, epsilon, *args):
1189 # calculate fprime(x0) first, as this may be cached by ScalarFunction
1190 f1 = fprime(*((x0,) + args))
1191 f2 = fprime(*((x0 + epsilon*p,) + args))
1192 return (f2 - f1) / epsilon
1195class _LineSearchError(RuntimeError):
1196 pass
1199def _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval,
1200 **kwargs):
1201 """
1202 Same as line_search_wolfe1, but fall back to line_search_wolfe2 if
1203 suitable step length is not found, and raise an exception if a
1204 suitable step length is not found.
1206 Raises
1207 ------
1208 _LineSearchError
1209 If no suitable step size is found
1211 """
1213 extra_condition = kwargs.pop('extra_condition', None)
1215 ret = line_search_wolfe1(f, fprime, xk, pk, gfk,
1216 old_fval, old_old_fval,
1217 **kwargs)
1219 if ret[0] is not None and extra_condition is not None:
1220 xp1 = xk + ret[0] * pk
1221 if not extra_condition(ret[0], xp1, ret[3], ret[5]):
1222 # Reject step if extra_condition fails
1223 ret = (None,)
1225 if ret[0] is None:
1226 # line search failed: try different one.
1227 with warnings.catch_warnings():
1228 warnings.simplefilter('ignore', LineSearchWarning)
1229 kwargs2 = {}
1230 for key in ('c1', 'c2', 'amax'):
1231 if key in kwargs:
1232 kwargs2[key] = kwargs[key]
1233 ret = line_search_wolfe2(f, fprime, xk, pk, gfk,
1234 old_fval, old_old_fval,
1235 extra_condition=extra_condition,
1236 **kwargs2)
1238 if ret[0] is None:
1239 raise _LineSearchError()
1241 return ret
1244def fmin_bfgs(f, x0, fprime=None, args=(), gtol=1e-5, norm=Inf,
1245 epsilon=_epsilon, maxiter=None, full_output=0, disp=1,
1246 retall=0, callback=None, xrtol=0):
1247 """
1248 Minimize a function using the BFGS algorithm.
1250 Parameters
1251 ----------
1252 f : callable ``f(x,*args)``
1253 Objective function to be minimized.
1254 x0 : ndarray
1255 Initial guess.
1256 fprime : callable ``f'(x,*args)``, optional
1257 Gradient of f.
1258 args : tuple, optional
1259 Extra arguments passed to f and fprime.
1260 gtol : float, optional
1261 Terminate successfully if gradient norm is less than `gtol`
1262 norm : float, optional
1263 Order of norm (Inf is max, -Inf is min)
1264 epsilon : int or ndarray, optional
1265 If `fprime` is approximated, use this value for the step size.
1266 callback : callable, optional
1267 An optional user-supplied function to call after each
1268 iteration. Called as ``callback(xk)``, where ``xk`` is the
1269 current parameter vector.
1270 maxiter : int, optional
1271 Maximum number of iterations to perform.
1272 full_output : bool, optional
1273 If True, return ``fopt``, ``func_calls``, ``grad_calls``, and
1274 ``warnflag`` in addition to ``xopt``.
1275 disp : bool, optional
1276 Print convergence message if True.
1277 retall : bool, optional
1278 Return a list of results at each iteration if True.
1279 xrtol : float, default: 0
1280 Relative tolerance for `x`. Terminate successfully if step
1281 size is less than ``xk * xrtol`` where ``xk`` is the current
1282 parameter vector.
1284 Returns
1285 -------
1286 xopt : ndarray
1287 Parameters which minimize f, i.e., ``f(xopt) == fopt``.
1288 fopt : float
1289 Minimum value.
1290 gopt : ndarray
1291 Value of gradient at minimum, f'(xopt), which should be near 0.
1292 Bopt : ndarray
1293 Value of 1/f''(xopt), i.e., the inverse Hessian matrix.
1294 func_calls : int
1295 Number of function_calls made.
1296 grad_calls : int
1297 Number of gradient calls made.
1298 warnflag : integer
1299 1 : Maximum number of iterations exceeded.
1300 2 : Gradient and/or function calls not changing.
1301 3 : NaN result encountered.
1302 allvecs : list
1303 The value of `xopt` at each iteration. Only returned if `retall` is
1304 True.
1306 Notes
1307 -----
1308 Optimize the function, `f`, whose gradient is given by `fprime`
1309 using the quasi-Newton method of Broyden, Fletcher, Goldfarb,
1310 and Shanno (BFGS).
1312 See Also
1313 --------
1314 minimize: Interface to minimization algorithms for multivariate
1315 functions. See ``method='BFGS'`` in particular.
1317 References
1318 ----------
1319 Wright, and Nocedal 'Numerical Optimization', 1999, p. 198.
1321 Examples
1322 --------
1323 >>> import numpy as np
1324 >>> from scipy.optimize import fmin_bfgs
1325 >>> def quadratic_cost(x, Q):
1326 ... return x @ Q @ x
1327 ...
1328 >>> x0 = np.array([-3, -4])
1329 >>> cost_weight = np.diag([1., 10.])
1330 >>> # Note that a trailing comma is necessary for a tuple with single element
1331 >>> fmin_bfgs(quadratic_cost, x0, args=(cost_weight,))
1332 Optimization terminated successfully.
1333 Current function value: 0.000000
1334 Iterations: 7 # may vary
1335 Function evaluations: 24 # may vary
1336 Gradient evaluations: 8 # may vary
1337 array([ 2.85169950e-06, -4.61820139e-07])
1339 >>> def quadratic_cost_grad(x, Q):
1340 ... return 2 * Q @ x
1341 ...
1342 >>> fmin_bfgs(quadratic_cost, x0, quadratic_cost_grad, args=(cost_weight,))
1343 Optimization terminated successfully.
1344 Current function value: 0.000000
1345 Iterations: 7
1346 Function evaluations: 8
1347 Gradient evaluations: 8
1348 array([ 2.85916637e-06, -4.54371951e-07])
1350 """
1351 opts = {'gtol': gtol,
1352 'norm': norm,
1353 'eps': epsilon,
1354 'disp': disp,
1355 'maxiter': maxiter,
1356 'return_all': retall,
1357 'xrtol': xrtol}
1359 callback = _wrap_callback(callback)
1360 res = _minimize_bfgs(f, x0, args, fprime, callback=callback, **opts)
1362 if full_output:
1363 retlist = (res['x'], res['fun'], res['jac'], res['hess_inv'],
1364 res['nfev'], res['njev'], res['status'])
1365 if retall:
1366 retlist += (res['allvecs'], )
1367 return retlist
1368 else:
1369 if retall:
1370 return res['x'], res['allvecs']
1371 else:
1372 return res['x']
1375def _minimize_bfgs(fun, x0, args=(), jac=None, callback=None,
1376 gtol=1e-5, norm=Inf, eps=_epsilon, maxiter=None,
1377 disp=False, return_all=False, finite_diff_rel_step=None,
1378 xrtol=0, **unknown_options):
1379 """
1380 Minimization of scalar function of one or more variables using the
1381 BFGS algorithm.
1383 Options
1384 -------
1385 disp : bool
1386 Set to True to print convergence messages.
1387 maxiter : int
1388 Maximum number of iterations to perform.
1389 gtol : float
1390 Terminate successfully if gradient norm is less than `gtol`.
1391 norm : float
1392 Order of norm (Inf is max, -Inf is min).
1393 eps : float or ndarray
1394 If `jac is None` the absolute step size used for numerical
1395 approximation of the jacobian via forward differences.
1396 return_all : bool, optional
1397 Set to True to return a list of the best solution at each of the
1398 iterations.
1399 finite_diff_rel_step : None or array_like, optional
1400 If `jac in ['2-point', '3-point', 'cs']` the relative step size to
1401 use for numerical approximation of the jacobian. The absolute step
1402 size is computed as ``h = rel_step * sign(x) * max(1, abs(x))``,
1403 possibly adjusted to fit into the bounds. For ``method='3-point'``
1404 the sign of `h` is ignored. If None (default) then step is selected
1405 automatically.
1406 xrtol : float, default: 0
1407 Relative tolerance for `x`. Terminate successfully if step size is
1408 less than ``xk * xrtol`` where ``xk`` is the current parameter vector.
1409 """
1410 _check_unknown_options(unknown_options)
1411 retall = return_all
1413 x0 = asarray(x0).flatten()
1414 if x0.ndim == 0:
1415 x0.shape = (1,)
1416 if maxiter is None:
1417 maxiter = len(x0) * 200
1419 sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,
1420 finite_diff_rel_step=finite_diff_rel_step)
1422 f = sf.fun
1423 myfprime = sf.grad
1425 old_fval = f(x0)
1426 gfk = myfprime(x0)
1428 k = 0
1429 N = len(x0)
1430 I = np.eye(N, dtype=int)
1431 Hk = I
1433 # Sets the initial step guess to dx ~ 1
1434 old_old_fval = old_fval + np.linalg.norm(gfk) / 2
1436 xk = x0
1437 if retall:
1438 allvecs = [x0]
1439 warnflag = 0
1440 gnorm = vecnorm(gfk, ord=norm)
1441 while (gnorm > gtol) and (k < maxiter):
1442 pk = -np.dot(Hk, gfk)
1443 try:
1444 alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
1445 _line_search_wolfe12(f, myfprime, xk, pk, gfk,
1446 old_fval, old_old_fval, amin=1e-100, amax=1e100)
1447 except _LineSearchError:
1448 # Line search failed to find a better solution.
1449 warnflag = 2
1450 break
1452 sk = alpha_k * pk
1453 xkp1 = xk + sk
1455 if retall:
1456 allvecs.append(xkp1)
1457 xk = xkp1
1458 if gfkp1 is None:
1459 gfkp1 = myfprime(xkp1)
1461 yk = gfkp1 - gfk
1462 gfk = gfkp1
1463 k += 1
1464 intermediate_result = OptimizeResult(x=xk, fun=old_fval)
1465 if _call_callback_maybe_halt(callback, intermediate_result):
1466 break
1467 gnorm = vecnorm(gfk, ord=norm)
1468 if (gnorm <= gtol):
1469 break
1471 # See Chapter 5 in P.E. Frandsen, K. Jonasson, H.B. Nielsen,
1472 # O. Tingleff: "Unconstrained Optimization", IMM, DTU. 1999.
1473 # These notes are available here:
1474 # http://www2.imm.dtu.dk/documents/ftp/publlec.html
1475 if (alpha_k*vecnorm(pk) <= xrtol*(xrtol + vecnorm(xk))):
1476 break
1478 if not np.isfinite(old_fval):
1479 # We correctly found +-Inf as optimal value, or something went
1480 # wrong.
1481 warnflag = 2
1482 break
1484 rhok_inv = np.dot(yk, sk)
1485 # this was handled in numeric, let it remaines for more safety
1486 # Cryptic comment above is preserved for posterity. Future reader:
1487 # consider change to condition below proposed in gh-1261/gh-17345.
1488 if rhok_inv == 0.:
1489 rhok = 1000.0
1490 if disp:
1491 msg = "Divide-by-zero encountered: rhok assumed large"
1492 _print_success_message_or_warn(True, msg)
1493 else:
1494 rhok = 1. / rhok_inv
1496 A1 = I - sk[:, np.newaxis] * yk[np.newaxis, :] * rhok
1497 A2 = I - yk[:, np.newaxis] * sk[np.newaxis, :] * rhok
1498 Hk = np.dot(A1, np.dot(Hk, A2)) + (rhok * sk[:, np.newaxis] *
1499 sk[np.newaxis, :])
1501 fval = old_fval
1503 if warnflag == 2:
1504 msg = _status_message['pr_loss']
1505 elif k >= maxiter:
1506 warnflag = 1
1507 msg = _status_message['maxiter']
1508 elif np.isnan(gnorm) or np.isnan(fval) or np.isnan(xk).any():
1509 warnflag = 3
1510 msg = _status_message['nan']
1511 else:
1512 msg = _status_message['success']
1514 if disp:
1515 _print_success_message_or_warn(warnflag, msg)
1516 print(" Current function value: %f" % fval)
1517 print(" Iterations: %d" % k)
1518 print(" Function evaluations: %d" % sf.nfev)
1519 print(" Gradient evaluations: %d" % sf.ngev)
1521 result = OptimizeResult(fun=fval, jac=gfk, hess_inv=Hk, nfev=sf.nfev,
1522 njev=sf.ngev, status=warnflag,
1523 success=(warnflag == 0), message=msg, x=xk,
1524 nit=k)
1525 if retall:
1526 result['allvecs'] = allvecs
1527 return result
1530def _print_success_message_or_warn(warnflag, message, warntype=None):
1531 if not warnflag:
1532 print(message)
1533 else:
1534 warnings.warn(message, warntype or OptimizeWarning, stacklevel=3)
1537def fmin_cg(f, x0, fprime=None, args=(), gtol=1e-5, norm=Inf, epsilon=_epsilon,
1538 maxiter=None, full_output=0, disp=1, retall=0, callback=None):
1539 """
1540 Minimize a function using a nonlinear conjugate gradient algorithm.
1542 Parameters
1543 ----------
1544 f : callable, ``f(x, *args)``
1545 Objective function to be minimized. Here `x` must be a 1-D array of
1546 the variables that are to be changed in the search for a minimum, and
1547 `args` are the other (fixed) parameters of `f`.
1548 x0 : ndarray
1549 A user-supplied initial estimate of `xopt`, the optimal value of `x`.
1550 It must be a 1-D array of values.
1551 fprime : callable, ``fprime(x, *args)``, optional
1552 A function that returns the gradient of `f` at `x`. Here `x` and `args`
1553 are as described above for `f`. The returned value must be a 1-D array.
1554 Defaults to None, in which case the gradient is approximated
1555 numerically (see `epsilon`, below).
1556 args : tuple, optional
1557 Parameter values passed to `f` and `fprime`. Must be supplied whenever
1558 additional fixed parameters are needed to completely specify the
1559 functions `f` and `fprime`.
1560 gtol : float, optional
1561 Stop when the norm of the gradient is less than `gtol`.
1562 norm : float, optional
1563 Order to use for the norm of the gradient
1564 (``-np.inf`` is min, ``np.inf`` is max).
1565 epsilon : float or ndarray, optional
1566 Step size(s) to use when `fprime` is approximated numerically. Can be a
1567 scalar or a 1-D array. Defaults to ``sqrt(eps)``, with eps the
1568 floating point machine precision. Usually ``sqrt(eps)`` is about
1569 1.5e-8.
1570 maxiter : int, optional
1571 Maximum number of iterations to perform. Default is ``200 * len(x0)``.
1572 full_output : bool, optional
1573 If True, return `fopt`, `func_calls`, `grad_calls`, and `warnflag` in
1574 addition to `xopt`. See the Returns section below for additional
1575 information on optional return values.
1576 disp : bool, optional
1577 If True, return a convergence message, followed by `xopt`.
1578 retall : bool, optional
1579 If True, add to the returned values the results of each iteration.
1580 callback : callable, optional
1581 An optional user-supplied function, called after each iteration.
1582 Called as ``callback(xk)``, where ``xk`` is the current value of `x0`.
1584 Returns
1585 -------
1586 xopt : ndarray
1587 Parameters which minimize f, i.e., ``f(xopt) == fopt``.
1588 fopt : float, optional
1589 Minimum value found, f(xopt). Only returned if `full_output` is True.
1590 func_calls : int, optional
1591 The number of function_calls made. Only returned if `full_output`
1592 is True.
1593 grad_calls : int, optional
1594 The number of gradient calls made. Only returned if `full_output` is
1595 True.
1596 warnflag : int, optional
1597 Integer value with warning status, only returned if `full_output` is
1598 True.
1600 0 : Success.
1602 1 : The maximum number of iterations was exceeded.
1604 2 : Gradient and/or function calls were not changing. May indicate
1605 that precision was lost, i.e., the routine did not converge.
1607 3 : NaN result encountered.
1609 allvecs : list of ndarray, optional
1610 List of arrays, containing the results at each iteration.
1611 Only returned if `retall` is True.
1613 See Also
1614 --------
1615 minimize : common interface to all `scipy.optimize` algorithms for
1616 unconstrained and constrained minimization of multivariate
1617 functions. It provides an alternative way to call
1618 ``fmin_cg``, by specifying ``method='CG'``.
1620 Notes
1621 -----
1622 This conjugate gradient algorithm is based on that of Polak and Ribiere
1623 [1]_.
1625 Conjugate gradient methods tend to work better when:
1627 1. `f` has a unique global minimizing point, and no local minima or
1628 other stationary points,
1629 2. `f` is, at least locally, reasonably well approximated by a
1630 quadratic function of the variables,
1631 3. `f` is continuous and has a continuous gradient,
1632 4. `fprime` is not too large, e.g., has a norm less than 1000,
1633 5. The initial guess, `x0`, is reasonably close to `f` 's global
1634 minimizing point, `xopt`.
1636 References
1637 ----------
1638 .. [1] Wright & Nocedal, "Numerical Optimization", 1999, pp. 120-122.
1640 Examples
1641 --------
1642 Example 1: seek the minimum value of the expression
1643 ``a*u**2 + b*u*v + c*v**2 + d*u + e*v + f`` for given values
1644 of the parameters and an initial guess ``(u, v) = (0, 0)``.
1646 >>> import numpy as np
1647 >>> args = (2, 3, 7, 8, 9, 10) # parameter values
1648 >>> def f(x, *args):
1649 ... u, v = x
1650 ... a, b, c, d, e, f = args
1651 ... return a*u**2 + b*u*v + c*v**2 + d*u + e*v + f
1652 >>> def gradf(x, *args):
1653 ... u, v = x
1654 ... a, b, c, d, e, f = args
1655 ... gu = 2*a*u + b*v + d # u-component of the gradient
1656 ... gv = b*u + 2*c*v + e # v-component of the gradient
1657 ... return np.asarray((gu, gv))
1658 >>> x0 = np.asarray((0, 0)) # Initial guess.
1659 >>> from scipy import optimize
1660 >>> res1 = optimize.fmin_cg(f, x0, fprime=gradf, args=args)
1661 Optimization terminated successfully.
1662 Current function value: 1.617021
1663 Iterations: 4
1664 Function evaluations: 8
1665 Gradient evaluations: 8
1666 >>> res1
1667 array([-1.80851064, -0.25531915])
1669 Example 2: solve the same problem using the `minimize` function.
1670 (This `myopts` dictionary shows all of the available options,
1671 although in practice only non-default values would be needed.
1672 The returned value will be a dictionary.)
1674 >>> opts = {'maxiter' : None, # default value.
1675 ... 'disp' : True, # non-default value.
1676 ... 'gtol' : 1e-5, # default value.
1677 ... 'norm' : np.inf, # default value.
1678 ... 'eps' : 1.4901161193847656e-08} # default value.
1679 >>> res2 = optimize.minimize(f, x0, jac=gradf, args=args,
1680 ... method='CG', options=opts)
1681 Optimization terminated successfully.
1682 Current function value: 1.617021
1683 Iterations: 4
1684 Function evaluations: 8
1685 Gradient evaluations: 8
1686 >>> res2.x # minimum found
1687 array([-1.80851064, -0.25531915])
1689 """
1690 opts = {'gtol': gtol,
1691 'norm': norm,
1692 'eps': epsilon,
1693 'disp': disp,
1694 'maxiter': maxiter,
1695 'return_all': retall}
1697 callback = _wrap_callback(callback)
1698 res = _minimize_cg(f, x0, args, fprime, callback=callback, **opts)
1700 if full_output:
1701 retlist = res['x'], res['fun'], res['nfev'], res['njev'], res['status']
1702 if retall:
1703 retlist += (res['allvecs'], )
1704 return retlist
1705 else:
1706 if retall:
1707 return res['x'], res['allvecs']
1708 else:
1709 return res['x']
1712def _minimize_cg(fun, x0, args=(), jac=None, callback=None,
1713 gtol=1e-5, norm=Inf, eps=_epsilon, maxiter=None,
1714 disp=False, return_all=False, finite_diff_rel_step=None,
1715 **unknown_options):
1716 """
1717 Minimization of scalar function of one or more variables using the
1718 conjugate gradient algorithm.
1720 Options
1721 -------
1722 disp : bool
1723 Set to True to print convergence messages.
1724 maxiter : int
1725 Maximum number of iterations to perform.
1726 gtol : float
1727 Gradient norm must be less than `gtol` before successful
1728 termination.
1729 norm : float
1730 Order of norm (Inf is max, -Inf is min).
1731 eps : float or ndarray
1732 If `jac is None` the absolute step size used for numerical
1733 approximation of the jacobian via forward differences.
1734 return_all : bool, optional
1735 Set to True to return a list of the best solution at each of the
1736 iterations.
1737 finite_diff_rel_step : None or array_like, optional
1738 If `jac in ['2-point', '3-point', 'cs']` the relative step size to
1739 use for numerical approximation of the jacobian. The absolute step
1740 size is computed as ``h = rel_step * sign(x) * max(1, abs(x))``,
1741 possibly adjusted to fit into the bounds. For ``method='3-point'``
1742 the sign of `h` is ignored. If None (default) then step is selected
1743 automatically.
1744 """
1745 _check_unknown_options(unknown_options)
1747 retall = return_all
1749 x0 = asarray(x0).flatten()
1750 if maxiter is None:
1751 maxiter = len(x0) * 200
1753 sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
1754 finite_diff_rel_step=finite_diff_rel_step)
1756 f = sf.fun
1757 myfprime = sf.grad
1759 old_fval = f(x0)
1760 gfk = myfprime(x0)
1762 k = 0
1763 xk = x0
1764 # Sets the initial step guess to dx ~ 1
1765 old_old_fval = old_fval + np.linalg.norm(gfk) / 2
1767 if retall:
1768 allvecs = [xk]
1769 warnflag = 0
1770 pk = -gfk
1771 gnorm = vecnorm(gfk, ord=norm)
1773 sigma_3 = 0.01
1775 while (gnorm > gtol) and (k < maxiter):
1776 deltak = np.dot(gfk, gfk)
1778 cached_step = [None]
1780 def polak_ribiere_powell_step(alpha, gfkp1=None):
1781 xkp1 = xk + alpha * pk
1782 if gfkp1 is None:
1783 gfkp1 = myfprime(xkp1)
1784 yk = gfkp1 - gfk
1785 beta_k = max(0, np.dot(yk, gfkp1) / deltak)
1786 pkp1 = -gfkp1 + beta_k * pk
1787 gnorm = vecnorm(gfkp1, ord=norm)
1788 return (alpha, xkp1, pkp1, gfkp1, gnorm)
1790 def descent_condition(alpha, xkp1, fp1, gfkp1):
1791 # Polak-Ribiere+ needs an explicit check of a sufficient
1792 # descent condition, which is not guaranteed by strong Wolfe.
1793 #
1794 # See Gilbert & Nocedal, "Global convergence properties of
1795 # conjugate gradient methods for optimization",
1796 # SIAM J. Optimization 2, 21 (1992).
1797 cached_step[:] = polak_ribiere_powell_step(alpha, gfkp1)
1798 alpha, xk, pk, gfk, gnorm = cached_step
1800 # Accept step if it leads to convergence.
1801 if gnorm <= gtol:
1802 return True
1804 # Accept step if sufficient descent condition applies.
1805 return np.dot(pk, gfk) <= -sigma_3 * np.dot(gfk, gfk)
1807 try:
1808 alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
1809 _line_search_wolfe12(f, myfprime, xk, pk, gfk, old_fval,
1810 old_old_fval, c2=0.4, amin=1e-100, amax=1e100,
1811 extra_condition=descent_condition)
1812 except _LineSearchError:
1813 # Line search failed to find a better solution.
1814 warnflag = 2
1815 break
1817 # Reuse already computed results if possible
1818 if alpha_k == cached_step[0]:
1819 alpha_k, xk, pk, gfk, gnorm = cached_step
1820 else:
1821 alpha_k, xk, pk, gfk, gnorm = polak_ribiere_powell_step(alpha_k, gfkp1)
1823 if retall:
1824 allvecs.append(xk)
1825 k += 1
1826 intermediate_result = OptimizeResult(x=xk, fun=old_fval)
1827 if _call_callback_maybe_halt(callback, intermediate_result):
1828 break
1830 fval = old_fval
1831 if warnflag == 2:
1832 msg = _status_message['pr_loss']
1833 elif k >= maxiter:
1834 warnflag = 1
1835 msg = _status_message['maxiter']
1836 elif np.isnan(gnorm) or np.isnan(fval) or np.isnan(xk).any():
1837 warnflag = 3
1838 msg = _status_message['nan']
1839 else:
1840 msg = _status_message['success']
1842 if disp:
1843 _print_success_message_or_warn(warnflag, msg)
1844 print(" Current function value: %f" % fval)
1845 print(" Iterations: %d" % k)
1846 print(" Function evaluations: %d" % sf.nfev)
1847 print(" Gradient evaluations: %d" % sf.ngev)
1849 result = OptimizeResult(fun=fval, jac=gfk, nfev=sf.nfev,
1850 njev=sf.ngev, status=warnflag,
1851 success=(warnflag == 0), message=msg, x=xk,
1852 nit=k)
1853 if retall:
1854 result['allvecs'] = allvecs
1855 return result
1858def fmin_ncg(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5,
1859 epsilon=_epsilon, maxiter=None, full_output=0, disp=1, retall=0,
1860 callback=None):
1861 """
1862 Unconstrained minimization of a function using the Newton-CG method.
1864 Parameters
1865 ----------
1866 f : callable ``f(x, *args)``
1867 Objective function to be minimized.
1868 x0 : ndarray
1869 Initial guess.
1870 fprime : callable ``f'(x, *args)``
1871 Gradient of f.
1872 fhess_p : callable ``fhess_p(x, p, *args)``, optional
1873 Function which computes the Hessian of f times an
1874 arbitrary vector, p.
1875 fhess : callable ``fhess(x, *args)``, optional
1876 Function to compute the Hessian matrix of f.
1877 args : tuple, optional
1878 Extra arguments passed to f, fprime, fhess_p, and fhess
1879 (the same set of extra arguments is supplied to all of
1880 these functions).
1881 epsilon : float or ndarray, optional
1882 If fhess is approximated, use this value for the step size.
1883 callback : callable, optional
1884 An optional user-supplied function which is called after
1885 each iteration. Called as callback(xk), where xk is the
1886 current parameter vector.
1887 avextol : float, optional
1888 Convergence is assumed when the average relative error in
1889 the minimizer falls below this amount.
1890 maxiter : int, optional
1891 Maximum number of iterations to perform.
1892 full_output : bool, optional
1893 If True, return the optional outputs.
1894 disp : bool, optional
1895 If True, print convergence message.
1896 retall : bool, optional
1897 If True, return a list of results at each iteration.
1899 Returns
1900 -------
1901 xopt : ndarray
1902 Parameters which minimize f, i.e., ``f(xopt) == fopt``.
1903 fopt : float
1904 Value of the function at xopt, i.e., ``fopt = f(xopt)``.
1905 fcalls : int
1906 Number of function calls made.
1907 gcalls : int
1908 Number of gradient calls made.
1909 hcalls : int
1910 Number of Hessian calls made.
1911 warnflag : int
1912 Warnings generated by the algorithm.
1913 1 : Maximum number of iterations exceeded.
1914 2 : Line search failure (precision loss).
1915 3 : NaN result encountered.
1916 allvecs : list
1917 The result at each iteration, if retall is True (see below).
1919 See also
1920 --------
1921 minimize: Interface to minimization algorithms for multivariate
1922 functions. See the 'Newton-CG' `method` in particular.
1924 Notes
1925 -----
1926 Only one of `fhess_p` or `fhess` need to be given. If `fhess`
1927 is provided, then `fhess_p` will be ignored. If neither `fhess`
1928 nor `fhess_p` is provided, then the hessian product will be
1929 approximated using finite differences on `fprime`. `fhess_p`
1930 must compute the hessian times an arbitrary vector. If it is not
1931 given, finite-differences on `fprime` are used to compute
1932 it.
1934 Newton-CG methods are also called truncated Newton methods. This
1935 function differs from scipy.optimize.fmin_tnc because
1937 1. scipy.optimize.fmin_ncg is written purely in Python using NumPy
1938 and scipy while scipy.optimize.fmin_tnc calls a C function.
1939 2. scipy.optimize.fmin_ncg is only for unconstrained minimization
1940 while scipy.optimize.fmin_tnc is for unconstrained minimization
1941 or box constrained minimization. (Box constraints give
1942 lower and upper bounds for each variable separately.)
1944 References
1945 ----------
1946 Wright & Nocedal, 'Numerical Optimization', 1999, p. 140.
1948 """
1949 opts = {'xtol': avextol,
1950 'eps': epsilon,
1951 'maxiter': maxiter,
1952 'disp': disp,
1953 'return_all': retall}
1955 callback = _wrap_callback(callback)
1956 res = _minimize_newtoncg(f, x0, args, fprime, fhess, fhess_p,
1957 callback=callback, **opts)
1959 if full_output:
1960 retlist = (res['x'], res['fun'], res['nfev'], res['njev'],
1961 res['nhev'], res['status'])
1962 if retall:
1963 retlist += (res['allvecs'], )
1964 return retlist
1965 else:
1966 if retall:
1967 return res['x'], res['allvecs']
1968 else:
1969 return res['x']
1972def _minimize_newtoncg(fun, x0, args=(), jac=None, hess=None, hessp=None,
1973 callback=None, xtol=1e-5, eps=_epsilon, maxiter=None,
1974 disp=False, return_all=False,
1975 **unknown_options):
1976 """
1977 Minimization of scalar function of one or more variables using the
1978 Newton-CG algorithm.
1980 Note that the `jac` parameter (Jacobian) is required.
1982 Options
1983 -------
1984 disp : bool
1985 Set to True to print convergence messages.
1986 xtol : float
1987 Average relative error in solution `xopt` acceptable for
1988 convergence.
1989 maxiter : int
1990 Maximum number of iterations to perform.
1991 eps : float or ndarray
1992 If `hessp` is approximated, use this value for the step size.
1993 return_all : bool, optional
1994 Set to True to return a list of the best solution at each of the
1995 iterations.
1996 """
1997 _check_unknown_options(unknown_options)
1998 if jac is None:
1999 raise ValueError('Jacobian is required for Newton-CG method')
2000 fhess_p = hessp
2001 fhess = hess
2002 avextol = xtol
2003 epsilon = eps
2004 retall = return_all
2006 x0 = asarray(x0).flatten()
2007 # TODO: add hessp (callable or FD) to ScalarFunction?
2008 sf = _prepare_scalar_function(
2009 fun, x0, jac, args=args, epsilon=eps, hess=hess
2010 )
2011 f = sf.fun
2012 fprime = sf.grad
2013 _h = sf.hess(x0)
2015 # Logic for hess/hessp
2016 # - If a callable(hess) is provided, then use that
2017 # - If hess is a FD_METHOD, or the output fom hess(x) is a LinearOperator
2018 # then create a hessp function using those.
2019 # - If hess is None but you have callable(hessp) then use the hessp.
2020 # - If hess and hessp are None then approximate hessp using the grad/jac.
2022 if (hess in FD_METHODS or isinstance(_h, LinearOperator)):
2023 fhess = None
2025 def _hessp(x, p, *args):
2026 return sf.hess(x).dot(p)
2028 fhess_p = _hessp
2030 def terminate(warnflag, msg):
2031 if disp:
2032 _print_success_message_or_warn(warnflag, msg)
2033 print(" Current function value: %f" % old_fval)
2034 print(" Iterations: %d" % k)
2035 print(" Function evaluations: %d" % sf.nfev)
2036 print(" Gradient evaluations: %d" % sf.ngev)
2037 print(" Hessian evaluations: %d" % hcalls)
2038 fval = old_fval
2039 result = OptimizeResult(fun=fval, jac=gfk, nfev=sf.nfev,
2040 njev=sf.ngev, nhev=hcalls, status=warnflag,
2041 success=(warnflag == 0), message=msg, x=xk,
2042 nit=k)
2043 if retall:
2044 result['allvecs'] = allvecs
2045 return result
2047 hcalls = 0
2048 if maxiter is None:
2049 maxiter = len(x0)*200
2050 cg_maxiter = 20*len(x0)
2052 xtol = len(x0) * avextol
2053 update = [2 * xtol]
2054 xk = x0
2055 if retall:
2056 allvecs = [xk]
2057 k = 0
2058 gfk = None
2059 old_fval = f(x0)
2060 old_old_fval = None
2061 float64eps = np.finfo(np.float64).eps
2062 while np.add.reduce(np.abs(update)) > xtol:
2063 if k >= maxiter:
2064 msg = "Warning: " + _status_message['maxiter']
2065 return terminate(1, msg)
2066 # Compute a search direction pk by applying the CG method to
2067 # del2 f(xk) p = - grad f(xk) starting from 0.
2068 b = -fprime(xk)
2069 maggrad = np.add.reduce(np.abs(b))
2070 eta = np.min([0.5, np.sqrt(maggrad)])
2071 termcond = eta * maggrad
2072 xsupi = zeros(len(x0), dtype=x0.dtype)
2073 ri = -b
2074 psupi = -ri
2075 i = 0
2076 dri0 = np.dot(ri, ri)
2078 if fhess is not None: # you want to compute hessian once.
2079 A = sf.hess(xk)
2080 hcalls = hcalls + 1
2082 for k2 in range(cg_maxiter):
2083 if np.add.reduce(np.abs(ri)) <= termcond:
2084 break
2085 if fhess is None:
2086 if fhess_p is None:
2087 Ap = approx_fhess_p(xk, psupi, fprime, epsilon)
2088 else:
2089 Ap = fhess_p(xk, psupi, *args)
2090 hcalls = hcalls + 1
2091 else:
2092 if isinstance(A, HessianUpdateStrategy):
2093 # if hess was supplied as a HessianUpdateStrategy
2094 Ap = A.dot(psupi)
2095 else:
2096 Ap = np.dot(A, psupi)
2097 # check curvature
2098 Ap = asarray(Ap).squeeze() # get rid of matrices...
2099 curv = np.dot(psupi, Ap)
2100 if 0 <= curv <= 3 * float64eps:
2101 break
2102 elif curv < 0:
2103 if (i > 0):
2104 break
2105 else:
2106 # fall back to steepest descent direction
2107 xsupi = dri0 / (-curv) * b
2108 break
2109 alphai = dri0 / curv
2110 xsupi = xsupi + alphai * psupi
2111 ri = ri + alphai * Ap
2112 dri1 = np.dot(ri, ri)
2113 betai = dri1 / dri0
2114 psupi = -ri + betai * psupi
2115 i = i + 1
2116 dri0 = dri1 # update np.dot(ri,ri) for next time.
2117 else:
2118 # curvature keeps increasing, bail out
2119 msg = ("Warning: CG iterations didn't converge. The Hessian is not "
2120 "positive definite.")
2121 return terminate(3, msg)
2123 pk = xsupi # search direction is solution to system.
2124 gfk = -b # gradient at xk
2126 try:
2127 alphak, fc, gc, old_fval, old_old_fval, gfkp1 = \
2128 _line_search_wolfe12(f, fprime, xk, pk, gfk,
2129 old_fval, old_old_fval)
2130 except _LineSearchError:
2131 # Line search failed to find a better solution.
2132 msg = "Warning: " + _status_message['pr_loss']
2133 return terminate(2, msg)
2135 update = alphak * pk
2136 xk = xk + update # upcast if necessary
2137 if retall:
2138 allvecs.append(xk)
2139 k += 1
2140 intermediate_result = OptimizeResult(x=xk, fun=old_fval)
2141 if _call_callback_maybe_halt(callback, intermediate_result):
2142 return terminate(5, "")
2144 else:
2145 if np.isnan(old_fval) or np.isnan(update).any():
2146 return terminate(3, _status_message['nan'])
2148 msg = _status_message['success']
2149 return terminate(0, msg)
2152def fminbound(func, x1, x2, args=(), xtol=1e-5, maxfun=500,
2153 full_output=0, disp=1):
2154 """Bounded minimization for scalar functions.
2156 Parameters
2157 ----------
2158 func : callable f(x,*args)
2159 Objective function to be minimized (must accept and return scalars).
2160 x1, x2 : float or array scalar
2161 Finite optimization bounds.
2162 args : tuple, optional
2163 Extra arguments passed to function.
2164 xtol : float, optional
2165 The convergence tolerance.
2166 maxfun : int, optional
2167 Maximum number of function evaluations allowed.
2168 full_output : bool, optional
2169 If True, return optional outputs.
2170 disp : int, optional
2171 If non-zero, print messages.
2172 0 : no message printing.
2173 1 : non-convergence notification messages only.
2174 2 : print a message on convergence too.
2175 3 : print iteration results.
2178 Returns
2179 -------
2180 xopt : ndarray
2181 Parameters (over given interval) which minimize the
2182 objective function.
2183 fval : number
2184 (Optional output) The function value evaluated at the minimizer.
2185 ierr : int
2186 (Optional output) An error flag (0 if converged, 1 if maximum number of
2187 function calls reached).
2188 numfunc : int
2189 (Optional output) The number of function calls made.
2191 See also
2192 --------
2193 minimize_scalar: Interface to minimization algorithms for scalar
2194 univariate functions. See the 'Bounded' `method` in particular.
2196 Notes
2197 -----
2198 Finds a local minimizer of the scalar function `func` in the
2199 interval x1 < xopt < x2 using Brent's method. (See `brent`
2200 for auto-bracketing.)
2202 References
2203 ----------
2204 .. [1] Forsythe, G.E., M. A. Malcolm, and C. B. Moler. "Computer Methods
2205 for Mathematical Computations." Prentice-Hall Series in Automatic
2206 Computation 259 (1977).
2207 .. [2] Brent, Richard P. Algorithms for Minimization Without Derivatives.
2208 Courier Corporation, 2013.
2210 Examples
2211 --------
2212 `fminbound` finds the minimizer of the function in the given range.
2213 The following examples illustrate this.
2215 >>> from scipy import optimize
2216 >>> def f(x):
2217 ... return (x-1)**2
2218 >>> minimizer = optimize.fminbound(f, -4, 4)
2219 >>> minimizer
2220 1.0
2221 >>> minimum = f(minimizer)
2222 >>> minimum
2223 0.0
2224 >>> res = optimize.fminbound(f, 3, 4, full_output=True)
2225 >>> minimizer, fval, ierr, numfunc = res
2226 >>> minimizer
2227 3.000005960860986
2228 >>> minimum = f(minimizer)
2229 >>> minimum, fval
2230 (4.000023843479476, 4.000023843479476)
2231 """
2232 options = {'xatol': xtol,
2233 'maxiter': maxfun,
2234 'disp': disp}
2236 res = _minimize_scalar_bounded(func, (x1, x2), args, **options)
2237 if full_output:
2238 return res['x'], res['fun'], res['status'], res['nfev']
2239 else:
2240 return res['x']
2243def _minimize_scalar_bounded(func, bounds, args=(),
2244 xatol=1e-5, maxiter=500, disp=0,
2245 **unknown_options):
2246 """
2247 Options
2248 -------
2249 maxiter : int
2250 Maximum number of iterations to perform.
2251 disp: int, optional
2252 If non-zero, print messages.
2253 0 : no message printing.
2254 1 : non-convergence notification messages only.
2255 2 : print a message on convergence too.
2256 3 : print iteration results.
2257 xatol : float
2258 Absolute error in solution `xopt` acceptable for convergence.
2260 """
2261 _check_unknown_options(unknown_options)
2262 maxfun = maxiter
2263 # Test bounds are of correct form
2264 if len(bounds) != 2:
2265 raise ValueError('bounds must have two elements.')
2266 x1, x2 = bounds
2268 if not (is_finite_scalar(x1) and is_finite_scalar(x2)):
2269 raise ValueError("Optimization bounds must be finite scalars.")
2271 if x1 > x2:
2272 raise ValueError("The lower bound exceeds the upper bound.")
2274 flag = 0
2275 header = ' Func-count x f(x) Procedure'
2276 step = ' initial'
2278 sqrt_eps = sqrt(2.2e-16)
2279 golden_mean = 0.5 * (3.0 - sqrt(5.0))
2280 a, b = x1, x2
2281 fulc = a + golden_mean * (b - a)
2282 nfc, xf = fulc, fulc
2283 rat = e = 0.0
2284 x = xf
2285 fx = func(x, *args)
2286 num = 1
2287 fmin_data = (1, xf, fx)
2288 fu = np.inf
2290 ffulc = fnfc = fx
2291 xm = 0.5 * (a + b)
2292 tol1 = sqrt_eps * np.abs(xf) + xatol / 3.0
2293 tol2 = 2.0 * tol1
2295 if disp > 2:
2296 print(" ")
2297 print(header)
2298 print("%5.0f %12.6g %12.6g %s" % (fmin_data + (step,)))
2300 while (np.abs(xf - xm) > (tol2 - 0.5 * (b - a))):
2301 golden = 1
2302 # Check for parabolic fit
2303 if np.abs(e) > tol1:
2304 golden = 0
2305 r = (xf - nfc) * (fx - ffulc)
2306 q = (xf - fulc) * (fx - fnfc)
2307 p = (xf - fulc) * q - (xf - nfc) * r
2308 q = 2.0 * (q - r)
2309 if q > 0.0:
2310 p = -p
2311 q = np.abs(q)
2312 r = e
2313 e = rat
2315 # Check for acceptability of parabola
2316 if ((np.abs(p) < np.abs(0.5*q*r)) and (p > q*(a - xf)) and
2317 (p < q * (b - xf))):
2318 rat = (p + 0.0) / q
2319 x = xf + rat
2320 step = ' parabolic'
2322 if ((x - a) < tol2) or ((b - x) < tol2):
2323 si = np.sign(xm - xf) + ((xm - xf) == 0)
2324 rat = tol1 * si
2325 else: # do a golden-section step
2326 golden = 1
2328 if golden: # do a golden-section step
2329 if xf >= xm:
2330 e = a - xf
2331 else:
2332 e = b - xf
2333 rat = golden_mean*e
2334 step = ' golden'
2336 si = np.sign(rat) + (rat == 0)
2337 x = xf + si * np.maximum(np.abs(rat), tol1)
2338 fu = func(x, *args)
2339 num += 1
2340 fmin_data = (num, x, fu)
2341 if disp > 2:
2342 print("%5.0f %12.6g %12.6g %s" % (fmin_data + (step,)))
2344 if fu <= fx:
2345 if x >= xf:
2346 a = xf
2347 else:
2348 b = xf
2349 fulc, ffulc = nfc, fnfc
2350 nfc, fnfc = xf, fx
2351 xf, fx = x, fu
2352 else:
2353 if x < xf:
2354 a = x
2355 else:
2356 b = x
2357 if (fu <= fnfc) or (nfc == xf):
2358 fulc, ffulc = nfc, fnfc
2359 nfc, fnfc = x, fu
2360 elif (fu <= ffulc) or (fulc == xf) or (fulc == nfc):
2361 fulc, ffulc = x, fu
2363 xm = 0.5 * (a + b)
2364 tol1 = sqrt_eps * np.abs(xf) + xatol / 3.0
2365 tol2 = 2.0 * tol1
2367 if num >= maxfun:
2368 flag = 1
2369 break
2371 if np.isnan(xf) or np.isnan(fx) or np.isnan(fu):
2372 flag = 2
2374 fval = fx
2375 if disp > 0:
2376 _endprint(x, flag, fval, maxfun, xatol, disp)
2378 result = OptimizeResult(fun=fval, status=flag, success=(flag == 0),
2379 message={0: 'Solution found.',
2380 1: 'Maximum number of function calls '
2381 'reached.',
2382 2: _status_message['nan']}.get(flag, ''),
2383 x=xf, nfev=num, nit=num)
2385 return result
2388class Brent:
2389 #need to rethink design of __init__
2390 def __init__(self, func, args=(), tol=1.48e-8, maxiter=500,
2391 full_output=0, disp=0):
2392 self.func = func
2393 self.args = args
2394 self.tol = tol
2395 self.maxiter = maxiter
2396 self._mintol = 1.0e-11
2397 self._cg = 0.3819660
2398 self.xmin = None
2399 self.fval = None
2400 self.iter = 0
2401 self.funcalls = 0
2402 self.disp = disp
2404 # need to rethink design of set_bracket (new options, etc.)
2405 def set_bracket(self, brack=None):
2406 self.brack = brack
2408 def get_bracket_info(self):
2409 #set up
2410 func = self.func
2411 args = self.args
2412 brack = self.brack
2413 ### BEGIN core bracket_info code ###
2414 ### carefully DOCUMENT any CHANGES in core ##
2415 if brack is None:
2416 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, args=args)
2417 elif len(brack) == 2:
2418 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],
2419 xb=brack[1], args=args)
2420 elif len(brack) == 3:
2421 xa, xb, xc = brack
2422 if (xa > xc): # swap so xa < xc can be assumed
2423 xc, xa = xa, xc
2424 if not ((xa < xb) and (xb < xc)):
2425 raise ValueError(
2426 "Bracketing values (xa, xb, xc) do not"
2427 " fulfill this requirement: (xa < xb) and (xb < xc)"
2428 )
2429 fa = func(*((xa,) + args))
2430 fb = func(*((xb,) + args))
2431 fc = func(*((xc,) + args))
2432 if not ((fb < fa) and (fb < fc)):
2433 raise ValueError(
2434 "Bracketing values (xa, xb, xc) do not fulfill"
2435 " this requirement: (f(xb) < f(xa)) and (f(xb) < f(xc))"
2436 )
2438 funcalls = 3
2439 else:
2440 raise ValueError("Bracketing interval must be "
2441 "length 2 or 3 sequence.")
2442 ### END core bracket_info code ###
2444 return xa, xb, xc, fa, fb, fc, funcalls
2446 def optimize(self):
2447 # set up for optimization
2448 func = self.func
2449 xa, xb, xc, fa, fb, fc, funcalls = self.get_bracket_info()
2450 _mintol = self._mintol
2451 _cg = self._cg
2452 #################################
2453 #BEGIN CORE ALGORITHM
2454 #################################
2455 x = w = v = xb
2456 fw = fv = fx = fb
2457 if (xa < xc):
2458 a = xa
2459 b = xc
2460 else:
2461 a = xc
2462 b = xa
2463 deltax = 0.0
2464 iter = 0
2466 if self.disp > 2:
2467 print(" ")
2468 print(f"{'Func-count':^12} {'x':^12} {'f(x)': ^12}")
2469 print(f"{funcalls:^12g} {x:^12.6g} {fx:^12.6g}")
2471 while (iter < self.maxiter):
2472 tol1 = self.tol * np.abs(x) + _mintol
2473 tol2 = 2.0 * tol1
2474 xmid = 0.5 * (a + b)
2475 # check for convergence
2476 if np.abs(x - xmid) < (tol2 - 0.5 * (b - a)):
2477 break
2478 # XXX In the first iteration, rat is only bound in the true case
2479 # of this conditional. This used to cause an UnboundLocalError
2480 # (gh-4140). It should be set before the if (but to what?).
2481 if (np.abs(deltax) <= tol1):
2482 if (x >= xmid):
2483 deltax = a - x # do a golden section step
2484 else:
2485 deltax = b - x
2486 rat = _cg * deltax
2487 else: # do a parabolic step
2488 tmp1 = (x - w) * (fx - fv)
2489 tmp2 = (x - v) * (fx - fw)
2490 p = (x - v) * tmp2 - (x - w) * tmp1
2491 tmp2 = 2.0 * (tmp2 - tmp1)
2492 if (tmp2 > 0.0):
2493 p = -p
2494 tmp2 = np.abs(tmp2)
2495 dx_temp = deltax
2496 deltax = rat
2497 # check parabolic fit
2498 if ((p > tmp2 * (a - x)) and (p < tmp2 * (b - x)) and
2499 (np.abs(p) < np.abs(0.5 * tmp2 * dx_temp))):
2500 rat = p * 1.0 / tmp2 # if parabolic step is useful.
2501 u = x + rat
2502 if ((u - a) < tol2 or (b - u) < tol2):
2503 if xmid - x >= 0:
2504 rat = tol1
2505 else:
2506 rat = -tol1
2507 else:
2508 if (x >= xmid):
2509 deltax = a - x # if it's not do a golden section step
2510 else:
2511 deltax = b - x
2512 rat = _cg * deltax
2514 if (np.abs(rat) < tol1): # update by at least tol1
2515 if rat >= 0:
2516 u = x + tol1
2517 else:
2518 u = x - tol1
2519 else:
2520 u = x + rat
2521 fu = func(*((u,) + self.args)) # calculate new output value
2522 funcalls += 1
2524 if (fu > fx): # if it's bigger than current
2525 if (u < x):
2526 a = u
2527 else:
2528 b = u
2529 if (fu <= fw) or (w == x):
2530 v = w
2531 w = u
2532 fv = fw
2533 fw = fu
2534 elif (fu <= fv) or (v == x) or (v == w):
2535 v = u
2536 fv = fu
2537 else:
2538 if (u >= x):
2539 a = x
2540 else:
2541 b = x
2542 v = w
2543 w = x
2544 x = u
2545 fv = fw
2546 fw = fx
2547 fx = fu
2549 if self.disp > 2:
2550 print(f"{funcalls:^12g} {x:^12.6g} {fx:^12.6g}")
2552 iter += 1
2553 #################################
2554 #END CORE ALGORITHM
2555 #################################
2557 self.xmin = x
2558 self.fval = fx
2559 self.iter = iter
2560 self.funcalls = funcalls
2562 def get_result(self, full_output=False):
2563 if full_output:
2564 return self.xmin, self.fval, self.iter, self.funcalls
2565 else:
2566 return self.xmin
2569def brent(func, args=(), brack=None, tol=1.48e-8, full_output=0, maxiter=500):
2570 """
2571 Given a function of one variable and a possible bracket, return
2572 a local minimizer of the function isolated to a fractional precision
2573 of tol.
2575 Parameters
2576 ----------
2577 func : callable f(x,*args)
2578 Objective function.
2579 args : tuple, optional
2580 Additional arguments (if present).
2581 brack : tuple, optional
2582 Either a triple ``(xa, xb, xc)`` satisfying ``xa < xb < xc`` and
2583 ``func(xb) < func(xa) and func(xb) < func(xc)``, or a pair
2584 ``(xa, xb)`` to be used as initial points for a downhill bracket search
2585 (see `scipy.optimize.bracket`).
2586 The minimizer ``x`` will not necessarily satisfy ``xa <= x <= xb``.
2587 tol : float, optional
2588 Relative error in solution `xopt` acceptable for convergence.
2589 full_output : bool, optional
2590 If True, return all output args (xmin, fval, iter,
2591 funcalls).
2592 maxiter : int, optional
2593 Maximum number of iterations in solution.
2595 Returns
2596 -------
2597 xmin : ndarray
2598 Optimum point.
2599 fval : float
2600 (Optional output) Optimum function value.
2601 iter : int
2602 (Optional output) Number of iterations.
2603 funcalls : int
2604 (Optional output) Number of objective function evaluations made.
2606 See also
2607 --------
2608 minimize_scalar: Interface to minimization algorithms for scalar
2609 univariate functions. See the 'Brent' `method` in particular.
2611 Notes
2612 -----
2613 Uses inverse parabolic interpolation when possible to speed up
2614 convergence of golden section method.
2616 Does not ensure that the minimum lies in the range specified by
2617 `brack`. See `scipy.optimize.fminbound`.
2619 Examples
2620 --------
2621 We illustrate the behaviour of the function when `brack` is of
2622 size 2 and 3 respectively. In the case where `brack` is of the
2623 form ``(xa, xb)``, we can see for the given values, the output does
2624 not necessarily lie in the range ``(xa, xb)``.
2626 >>> def f(x):
2627 ... return (x-1)**2
2629 >>> from scipy import optimize
2631 >>> minimizer = optimize.brent(f, brack=(1, 2))
2632 >>> minimizer
2633 1
2634 >>> res = optimize.brent(f, brack=(-1, 0.5, 2), full_output=True)
2635 >>> xmin, fval, iter, funcalls = res
2636 >>> f(xmin), fval
2637 (0.0, 0.0)
2639 """
2640 options = {'xtol': tol,
2641 'maxiter': maxiter}
2642 res = _minimize_scalar_brent(func, brack, args, **options)
2643 if full_output:
2644 return res['x'], res['fun'], res['nit'], res['nfev']
2645 else:
2646 return res['x']
2649def _minimize_scalar_brent(func, brack=None, args=(), xtol=1.48e-8,
2650 maxiter=500, disp=0,
2651 **unknown_options):
2652 """
2653 Options
2654 -------
2655 maxiter : int
2656 Maximum number of iterations to perform.
2657 xtol : float
2658 Relative error in solution `xopt` acceptable for convergence.
2659 disp: int, optional
2660 If non-zero, print messages.
2661 0 : no message printing.
2662 1 : non-convergence notification messages only.
2663 2 : print a message on convergence too.
2664 3 : print iteration results.
2665 Notes
2666 -----
2667 Uses inverse parabolic interpolation when possible to speed up
2668 convergence of golden section method.
2670 """
2671 _check_unknown_options(unknown_options)
2672 tol = xtol
2673 if tol < 0:
2674 raise ValueError('tolerance should be >= 0, got %r' % tol)
2676 brent = Brent(func=func, args=args, tol=tol,
2677 full_output=True, maxiter=maxiter, disp=disp)
2678 brent.set_bracket(brack)
2679 brent.optimize()
2680 x, fval, nit, nfev = brent.get_result(full_output=True)
2682 success = nit < maxiter and not (np.isnan(x) or np.isnan(fval))
2684 if success:
2685 message = ("\nOptimization terminated successfully;\n"
2686 "The returned value satisfies the termination criteria\n"
2687 f"(using xtol = {xtol} )")
2688 else:
2689 if nit >= maxiter:
2690 message = "\nMaximum number of iterations exceeded"
2691 if np.isnan(x) or np.isnan(fval):
2692 message = f"{_status_message['nan']}"
2694 if disp:
2695 _print_success_message_or_warn(not success, message)
2697 return OptimizeResult(fun=fval, x=x, nit=nit, nfev=nfev,
2698 success=success, message=message)
2701def golden(func, args=(), brack=None, tol=_epsilon,
2702 full_output=0, maxiter=5000):
2703 """
2704 Return the minimizer of a function of one variable using the golden section
2705 method.
2707 Given a function of one variable and a possible bracketing interval,
2708 return a minimizer of the function isolated to a fractional precision of
2709 tol.
2711 Parameters
2712 ----------
2713 func : callable func(x,*args)
2714 Objective function to minimize.
2715 args : tuple, optional
2716 Additional arguments (if present), passed to func.
2717 brack : tuple, optional
2718 Either a triple ``(xa, xb, xc)`` where ``xa < xb < xc`` and
2719 ``func(xb) < func(xa) and func(xb) < func(xc)``, or a pair (xa, xb)
2720 to be used as initial points for a downhill bracket search (see
2721 `scipy.optimize.bracket`).
2722 The minimizer ``x`` will not necessarily satisfy ``xa <= x <= xb``.
2723 tol : float, optional
2724 x tolerance stop criterion
2725 full_output : bool, optional
2726 If True, return optional outputs.
2727 maxiter : int
2728 Maximum number of iterations to perform.
2730 Returns
2731 -------
2732 xmin : ndarray
2733 Optimum point.
2734 fval : float
2735 (Optional output) Optimum function value.
2736 funcalls : int
2737 (Optional output) Number of objective function evaluations made.
2739 See also
2740 --------
2741 minimize_scalar: Interface to minimization algorithms for scalar
2742 univariate functions. See the 'Golden' `method` in particular.
2744 Notes
2745 -----
2746 Uses analog of bisection method to decrease the bracketed
2747 interval.
2749 Examples
2750 --------
2751 We illustrate the behaviour of the function when `brack` is of
2752 size 2 and 3, respectively. In the case where `brack` is of the
2753 form (xa,xb), we can see for the given values, the output need
2754 not necessarily lie in the range ``(xa, xb)``.
2756 >>> def f(x):
2757 ... return (x-1)**2
2759 >>> from scipy import optimize
2761 >>> minimizer = optimize.golden(f, brack=(1, 2))
2762 >>> minimizer
2763 1
2764 >>> res = optimize.golden(f, brack=(-1, 0.5, 2), full_output=True)
2765 >>> xmin, fval, funcalls = res
2766 >>> f(xmin), fval
2767 (9.925165290385052e-18, 9.925165290385052e-18)
2769 """
2770 options = {'xtol': tol, 'maxiter': maxiter}
2771 res = _minimize_scalar_golden(func, brack, args, **options)
2772 if full_output:
2773 return res['x'], res['fun'], res['nfev']
2774 else:
2775 return res['x']
2778def _minimize_scalar_golden(func, brack=None, args=(),
2779 xtol=_epsilon, maxiter=5000, disp=0,
2780 **unknown_options):
2781 """
2782 Options
2783 -------
2784 xtol : float
2785 Relative error in solution `xopt` acceptable for convergence.
2786 maxiter : int
2787 Maximum number of iterations to perform.
2788 disp: int, optional
2789 If non-zero, print messages.
2790 0 : no message printing.
2791 1 : non-convergence notification messages only.
2792 2 : print a message on convergence too.
2793 3 : print iteration results.
2794 """
2795 _check_unknown_options(unknown_options)
2796 tol = xtol
2797 if brack is None:
2798 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, args=args)
2799 elif len(brack) == 2:
2800 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],
2801 xb=brack[1], args=args)
2802 elif len(brack) == 3:
2803 xa, xb, xc = brack
2804 if (xa > xc): # swap so xa < xc can be assumed
2805 xc, xa = xa, xc
2806 if not ((xa < xb) and (xb < xc)):
2807 raise ValueError(
2808 "Bracketing values (xa, xb, xc) do not"
2809 " fulfill this requirement: (xa < xb) and (xb < xc)"
2810 )
2811 fa = func(*((xa,) + args))
2812 fb = func(*((xb,) + args))
2813 fc = func(*((xc,) + args))
2814 if not ((fb < fa) and (fb < fc)):
2815 raise ValueError(
2816 "Bracketing values (xa, xb, xc) do not fulfill"
2817 " this requirement: (f(xb) < f(xa)) and (f(xb) < f(xc))"
2818 )
2819 funcalls = 3
2820 else:
2821 raise ValueError("Bracketing interval must be length 2 or 3 sequence.")
2823 _gR = 0.61803399 # golden ratio conjugate: 2.0/(1.0+sqrt(5.0))
2824 _gC = 1.0 - _gR
2825 x3 = xc
2826 x0 = xa
2827 if (np.abs(xc - xb) > np.abs(xb - xa)):
2828 x1 = xb
2829 x2 = xb + _gC * (xc - xb)
2830 else:
2831 x2 = xb
2832 x1 = xb - _gC * (xb - xa)
2833 f1 = func(*((x1,) + args))
2834 f2 = func(*((x2,) + args))
2835 funcalls += 2
2836 nit = 0
2838 if disp > 2:
2839 print(" ")
2840 print(f"{'Func-count':^12} {'x':^12} {'f(x)': ^12}")
2842 for i in range(maxiter):
2843 if np.abs(x3 - x0) <= tol * (np.abs(x1) + np.abs(x2)):
2844 break
2845 if (f2 < f1):
2846 x0 = x1
2847 x1 = x2
2848 x2 = _gR * x1 + _gC * x3
2849 f1 = f2
2850 f2 = func(*((x2,) + args))
2851 else:
2852 x3 = x2
2853 x2 = x1
2854 x1 = _gR * x2 + _gC * x0
2855 f2 = f1
2856 f1 = func(*((x1,) + args))
2857 funcalls += 1
2858 if disp > 2:
2859 if (f1 < f2):
2860 xmin, fval = x1, f1
2861 else:
2862 xmin, fval = x2, f2
2863 print(f"{funcalls:^12g} {xmin:^12.6g} {fval:^12.6g}")
2865 nit += 1
2866 # end of iteration loop
2868 if (f1 < f2):
2869 xmin = x1
2870 fval = f1
2871 else:
2872 xmin = x2
2873 fval = f2
2875 success = nit < maxiter and not (np.isnan(fval) or np.isnan(xmin))
2877 if success:
2878 message = ("\nOptimization terminated successfully;\n"
2879 "The returned value satisfies the termination criteria\n"
2880 f"(using xtol = {xtol} )")
2881 else:
2882 if nit >= maxiter:
2883 message = "\nMaximum number of iterations exceeded"
2884 if np.isnan(xmin) or np.isnan(fval):
2885 message = f"{_status_message['nan']}"
2887 if disp:
2888 _print_success_message_or_warn(not success, message)
2890 return OptimizeResult(fun=fval, nfev=funcalls, x=xmin, nit=nit,
2891 success=success, message=message)
2894def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
2895 """
2896 Bracket the minimum of a function.
2898 Given a function and distinct initial points, search in the
2899 downhill direction (as defined by the initial points) and return
2900 three points that bracket the minimum of the function.
2902 Parameters
2903 ----------
2904 func : callable f(x,*args)
2905 Objective function to minimize.
2906 xa, xb : float, optional
2907 Initial points. Defaults `xa` to 0.0, and `xb` to 1.0.
2908 A local minimum need not be contained within this interval.
2909 args : tuple, optional
2910 Additional arguments (if present), passed to `func`.
2911 grow_limit : float, optional
2912 Maximum grow limit. Defaults to 110.0
2913 maxiter : int, optional
2914 Maximum number of iterations to perform. Defaults to 1000.
2916 Returns
2917 -------
2918 xa, xb, xc : float
2919 Final points of the bracket.
2920 fa, fb, fc : float
2921 Objective function values at the bracket points.
2922 funcalls : int
2923 Number of function evaluations made.
2925 Raises
2926 ------
2927 BracketError
2928 If no valid bracket is found before the algorithm terminates.
2929 See notes for conditions of a valid bracket.
2931 Notes
2932 -----
2933 The algorithm attempts to find three strictly ordered points (i.e.
2934 :math:`x_a < x_b < x_c` or :math:`x_c < x_b < x_a`) satisfying
2935 :math:`f(x_b) ≤ f(x_a)` and :math:`f(x_b) ≤ f(x_c)`, where one of the
2936 inequalities must be satistfied strictly and all :math:`x_i` must be
2937 finite.
2939 Examples
2940 --------
2941 This function can find a downward convex region of a function:
2943 >>> import numpy as np
2944 >>> import matplotlib.pyplot as plt
2945 >>> from scipy.optimize import bracket
2946 >>> def f(x):
2947 ... return 10*x**2 + 3*x + 5
2948 >>> x = np.linspace(-2, 2)
2949 >>> y = f(x)
2950 >>> init_xa, init_xb = 0.1, 1
2951 >>> xa, xb, xc, fa, fb, fc, funcalls = bracket(f, xa=init_xa, xb=init_xb)
2952 >>> plt.axvline(x=init_xa, color="k", linestyle="--")
2953 >>> plt.axvline(x=init_xb, color="k", linestyle="--")
2954 >>> plt.plot(x, y, "-k")
2955 >>> plt.plot(xa, fa, "bx")
2956 >>> plt.plot(xb, fb, "rx")
2957 >>> plt.plot(xc, fc, "bx")
2958 >>> plt.show()
2960 Note that both initial points were to the right of the minimum, and the
2961 third point was found in the "downhill" direction: the direction
2962 in which the function appeared to be decreasing (to the left).
2963 The final points are strictly ordered, and the function value
2964 at the middle point is less than the function values at the endpoints;
2965 it follows that a minimum must lie within the bracket.
2967 """
2968 _gold = 1.618034 # golden ratio: (1.0+sqrt(5.0))/2.0
2969 _verysmall_num = 1e-21
2970 # convert to numpy floats if not already
2971 xa, xb = np.asarray([xa, xb])
2972 fa = func(*(xa,) + args)
2973 fb = func(*(xb,) + args)
2974 if (fa < fb): # Switch so fa > fb
2975 xa, xb = xb, xa
2976 fa, fb = fb, fa
2977 xc = xb + _gold * (xb - xa)
2978 fc = func(*((xc,) + args))
2979 funcalls = 3
2980 iter = 0
2981 while (fc < fb):
2982 tmp1 = (xb - xa) * (fb - fc)
2983 tmp2 = (xb - xc) * (fb - fa)
2984 val = tmp2 - tmp1
2985 if np.abs(val) < _verysmall_num:
2986 denom = 2.0 * _verysmall_num
2987 else:
2988 denom = 2.0 * val
2989 w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom
2990 wlim = xb + grow_limit * (xc - xb)
2991 msg = ("No valid bracket was found before the iteration limit was "
2992 "reached. Consider trying different initial points or "
2993 "increasing `maxiter`.")
2994 if iter > maxiter:
2995 raise RuntimeError(msg)
2996 iter += 1
2997 if (w - xc) * (xb - w) > 0.0:
2998 fw = func(*((w,) + args))
2999 funcalls += 1
3000 if (fw < fc):
3001 xa = xb
3002 xb = w
3003 fa = fb
3004 fb = fw
3005 break
3006 elif (fw > fb):
3007 xc = w
3008 fc = fw
3009 break
3010 w = xc + _gold * (xc - xb)
3011 fw = func(*((w,) + args))
3012 funcalls += 1
3013 elif (w - wlim)*(wlim - xc) >= 0.0:
3014 w = wlim
3015 fw = func(*((w,) + args))
3016 funcalls += 1
3017 elif (w - wlim)*(xc - w) > 0.0:
3018 fw = func(*((w,) + args))
3019 funcalls += 1
3020 if (fw < fc):
3021 xb = xc
3022 xc = w
3023 w = xc + _gold * (xc - xb)
3024 fb = fc
3025 fc = fw
3026 fw = func(*((w,) + args))
3027 funcalls += 1
3028 else:
3029 w = xc + _gold * (xc - xb)
3030 fw = func(*((w,) + args))
3031 funcalls += 1
3032 xa = xb
3033 xb = xc
3034 xc = w
3035 fa = fb
3036 fb = fc
3037 fc = fw
3039 # three conditions for a valid bracket
3040 cond1 = (fb < fc and fb <= fa) or (fb < fa and fb <= fc)
3041 cond2 = (xa < xb < xc or xc < xb < xa)
3042 cond3 = np.isfinite(xa) and np.isfinite(xb) and np.isfinite(xc)
3043 msg = ("The algorithm terminated without finding a valid bracket. "
3044 "Consider trying different initial points.")
3045 if not (cond1 and cond2 and cond3):
3046 e = BracketError(msg)
3047 e.data = (xa, xb, xc, fa, fb, fc, funcalls)
3048 raise e
3050 return xa, xb, xc, fa, fb, fc, funcalls
3053class BracketError(RuntimeError):
3054 pass
3057def _recover_from_bracket_error(solver, fun, bracket, args, **options):
3058 # `bracket` was originally written without checking whether the resulting
3059 # bracket is valid. `brent` and `golden` built on top of it without
3060 # checking the returned bracket for validity, and their output can be
3061 # incorrect without warning/error if the original bracket is invalid.
3062 # gh-14858 noticed the problem, and the following is the desired
3063 # behavior:
3064 # - `scipy.optimize.bracket`, `scipy.optimize.brent`, and
3065 # `scipy.optimize.golden` should raise an error if the bracket is
3066 # invalid, as opposed to silently returning garbage
3067 # - `scipy.optimize.minimize_scalar` should return with `success=False`
3068 # and other information
3069 # The changes that would be required to achieve this the traditional
3070 # way (`return`ing all the required information from bracket all the way
3071 # up to `minimizer_scalar`) are extensive and invasive. (See a6aa40d.)
3072 # We can achieve the same thing by raising the error in `bracket`, but
3073 # storing the information needed by `minimize_scalar` in the error object,
3074 # and intercepting it here.
3075 try:
3076 res = solver(fun, bracket, args, **options)
3077 except BracketError as e:
3078 msg = str(e)
3079 xa, xb, xc, fa, fb, fc, funcalls = e.data
3080 xs, fs = [xa, xb, xc], [fa, fb, fc]
3081 if np.any(np.isnan([xs, fs])):
3082 x, fun = np.nan, np.nan
3083 else:
3084 imin = np.argmin(fs)
3085 x, fun = xs[imin], fs[imin]
3086 return OptimizeResult(fun=fun, nfev=funcalls, x=x,
3087 nit=0, success=False, message=msg)
3088 return res
3091def _line_for_search(x0, alpha, lower_bound, upper_bound):
3092 """
3093 Given a parameter vector ``x0`` with length ``n`` and a direction
3094 vector ``alpha`` with length ``n``, and lower and upper bounds on
3095 each of the ``n`` parameters, what are the bounds on a scalar
3096 ``l`` such that ``lower_bound <= x0 + alpha * l <= upper_bound``.
3099 Parameters
3100 ----------
3101 x0 : np.array.
3102 The vector representing the current location.
3103 Note ``np.shape(x0) == (n,)``.
3104 alpha : np.array.
3105 The vector representing the direction.
3106 Note ``np.shape(alpha) == (n,)``.
3107 lower_bound : np.array.
3108 The lower bounds for each parameter in ``x0``. If the ``i``th
3109 parameter in ``x0`` is unbounded below, then ``lower_bound[i]``
3110 should be ``-np.inf``.
3111 Note ``np.shape(lower_bound) == (n,)``.
3112 upper_bound : np.array.
3113 The upper bounds for each parameter in ``x0``. If the ``i``th
3114 parameter in ``x0`` is unbounded above, then ``upper_bound[i]``
3115 should be ``np.inf``.
3116 Note ``np.shape(upper_bound) == (n,)``.
3118 Returns
3119 -------
3120 res : tuple ``(lmin, lmax)``
3121 The bounds for ``l`` such that
3122 ``lower_bound[i] <= x0[i] + alpha[i] * l <= upper_bound[i]``
3123 for all ``i``.
3125 """
3126 # get nonzero indices of alpha so we don't get any zero division errors.
3127 # alpha will not be all zero, since it is called from _linesearch_powell
3128 # where we have a check for this.
3129 nonzero, = alpha.nonzero()
3130 lower_bound, upper_bound = lower_bound[nonzero], upper_bound[nonzero]
3131 x0, alpha = x0[nonzero], alpha[nonzero]
3132 low = (lower_bound - x0) / alpha
3133 high = (upper_bound - x0) / alpha
3135 # positive and negative indices
3136 pos = alpha > 0
3138 lmin_pos = np.where(pos, low, 0)
3139 lmin_neg = np.where(pos, 0, high)
3140 lmax_pos = np.where(pos, high, 0)
3141 lmax_neg = np.where(pos, 0, low)
3143 lmin = np.max(lmin_pos + lmin_neg)
3144 lmax = np.min(lmax_pos + lmax_neg)
3146 # if x0 is outside the bounds, then it is possible that there is
3147 # no way to get back in the bounds for the parameters being updated
3148 # with the current direction alpha.
3149 # when this happens, lmax < lmin.
3150 # If this is the case, then we can just return (0, 0)
3151 return (lmin, lmax) if lmax >= lmin else (0, 0)
3154def _linesearch_powell(func, p, xi, tol=1e-3,
3155 lower_bound=None, upper_bound=None, fval=None):
3156 """Line-search algorithm using fminbound.
3158 Find the minimium of the function ``func(x0 + alpha*direc)``.
3160 lower_bound : np.array.
3161 The lower bounds for each parameter in ``x0``. If the ``i``th
3162 parameter in ``x0`` is unbounded below, then ``lower_bound[i]``
3163 should be ``-np.inf``.
3164 Note ``np.shape(lower_bound) == (n,)``.
3165 upper_bound : np.array.
3166 The upper bounds for each parameter in ``x0``. If the ``i``th
3167 parameter in ``x0`` is unbounded above, then ``upper_bound[i]``
3168 should be ``np.inf``.
3169 Note ``np.shape(upper_bound) == (n,)``.
3170 fval : number.
3171 ``fval`` is equal to ``func(p)``, the idea is just to avoid
3172 recomputing it so we can limit the ``fevals``.
3174 """
3175 def myfunc(alpha):
3176 return func(p + alpha*xi)
3178 # if xi is zero, then don't optimize
3179 if not np.any(xi):
3180 return ((fval, p, xi) if fval is not None else (func(p), p, xi))
3181 elif lower_bound is None and upper_bound is None:
3182 # non-bounded minimization
3183 res = _recover_from_bracket_error(_minimize_scalar_brent,
3184 myfunc, None, tuple(), xtol=tol)
3185 alpha_min, fret = res.x, res.fun
3186 xi = alpha_min * xi
3187 return squeeze(fret), p + xi, xi
3188 else:
3189 bound = _line_for_search(p, xi, lower_bound, upper_bound)
3190 if np.isneginf(bound[0]) and np.isposinf(bound[1]):
3191 # equivalent to unbounded
3192 return _linesearch_powell(func, p, xi, fval=fval, tol=tol)
3193 elif not np.isneginf(bound[0]) and not np.isposinf(bound[1]):
3194 # we can use a bounded scalar minimization
3195 res = _minimize_scalar_bounded(myfunc, bound, xatol=tol / 100)
3196 xi = res.x * xi
3197 return squeeze(res.fun), p + xi, xi
3198 else:
3199 # only bounded on one side. use the tangent function to convert
3200 # the infinity bound to a finite bound. The new bounded region
3201 # is a subregion of the region bounded by -np.pi/2 and np.pi/2.
3202 bound = np.arctan(bound[0]), np.arctan(bound[1])
3203 res = _minimize_scalar_bounded(
3204 lambda x: myfunc(np.tan(x)),
3205 bound,
3206 xatol=tol / 100)
3207 xi = np.tan(res.x) * xi
3208 return squeeze(res.fun), p + xi, xi
3211def fmin_powell(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None,
3212 maxfun=None, full_output=0, disp=1, retall=0, callback=None,
3213 direc=None):
3214 """
3215 Minimize a function using modified Powell's method.
3217 This method only uses function values, not derivatives.
3219 Parameters
3220 ----------
3221 func : callable f(x,*args)
3222 Objective function to be minimized.
3223 x0 : ndarray
3224 Initial guess.
3225 args : tuple, optional
3226 Extra arguments passed to func.
3227 xtol : float, optional
3228 Line-search error tolerance.
3229 ftol : float, optional
3230 Relative error in ``func(xopt)`` acceptable for convergence.
3231 maxiter : int, optional
3232 Maximum number of iterations to perform.
3233 maxfun : int, optional
3234 Maximum number of function evaluations to make.
3235 full_output : bool, optional
3236 If True, ``fopt``, ``xi``, ``direc``, ``iter``, ``funcalls``, and
3237 ``warnflag`` are returned.
3238 disp : bool, optional
3239 If True, print convergence messages.
3240 retall : bool, optional
3241 If True, return a list of the solution at each iteration.
3242 callback : callable, optional
3243 An optional user-supplied function, called after each
3244 iteration. Called as ``callback(xk)``, where ``xk`` is the
3245 current parameter vector.
3246 direc : ndarray, optional
3247 Initial fitting step and parameter order set as an (N, N) array, where N
3248 is the number of fitting parameters in `x0`. Defaults to step size 1.0
3249 fitting all parameters simultaneously (``np.eye((N, N))``). To
3250 prevent initial consideration of values in a step or to change initial
3251 step size, set to 0 or desired step size in the Jth position in the Mth
3252 block, where J is the position in `x0` and M is the desired evaluation
3253 step, with steps being evaluated in index order. Step size and ordering
3254 will change freely as minimization proceeds.
3256 Returns
3257 -------
3258 xopt : ndarray
3259 Parameter which minimizes `func`.
3260 fopt : number
3261 Value of function at minimum: ``fopt = func(xopt)``.
3262 direc : ndarray
3263 Current direction set.
3264 iter : int
3265 Number of iterations.
3266 funcalls : int
3267 Number of function calls made.
3268 warnflag : int
3269 Integer warning flag:
3270 1 : Maximum number of function evaluations.
3271 2 : Maximum number of iterations.
3272 3 : NaN result encountered.
3273 4 : The result is out of the provided bounds.
3274 allvecs : list
3275 List of solutions at each iteration.
3277 See also
3278 --------
3279 minimize: Interface to unconstrained minimization algorithms for
3280 multivariate functions. See the 'Powell' method in particular.
3282 Notes
3283 -----
3284 Uses a modification of Powell's method to find the minimum of
3285 a function of N variables. Powell's method is a conjugate
3286 direction method.
3288 The algorithm has two loops. The outer loop merely iterates over the inner
3289 loop. The inner loop minimizes over each current direction in the direction
3290 set. At the end of the inner loop, if certain conditions are met, the
3291 direction that gave the largest decrease is dropped and replaced with the
3292 difference between the current estimated x and the estimated x from the
3293 beginning of the inner-loop.
3295 The technical conditions for replacing the direction of greatest
3296 increase amount to checking that
3298 1. No further gain can be made along the direction of greatest increase
3299 from that iteration.
3300 2. The direction of greatest increase accounted for a large sufficient
3301 fraction of the decrease in the function value from that iteration of
3302 the inner loop.
3304 References
3305 ----------
3306 Powell M.J.D. (1964) An efficient method for finding the minimum of a
3307 function of several variables without calculating derivatives,
3308 Computer Journal, 7 (2):155-162.
3310 Press W., Teukolsky S.A., Vetterling W.T., and Flannery B.P.:
3311 Numerical Recipes (any edition), Cambridge University Press
3313 Examples
3314 --------
3315 >>> def f(x):
3316 ... return x**2
3318 >>> from scipy import optimize
3320 >>> minimum = optimize.fmin_powell(f, -1)
3321 Optimization terminated successfully.
3322 Current function value: 0.000000
3323 Iterations: 2
3324 Function evaluations: 16
3325 >>> minimum
3326 array(0.0)
3328 """
3329 opts = {'xtol': xtol,
3330 'ftol': ftol,
3331 'maxiter': maxiter,
3332 'maxfev': maxfun,
3333 'disp': disp,
3334 'direc': direc,
3335 'return_all': retall}
3337 callback = _wrap_callback(callback)
3338 res = _minimize_powell(func, x0, args, callback=callback, **opts)
3340 if full_output:
3341 retlist = (res['x'], res['fun'], res['direc'], res['nit'],
3342 res['nfev'], res['status'])
3343 if retall:
3344 retlist += (res['allvecs'], )
3345 return retlist
3346 else:
3347 if retall:
3348 return res['x'], res['allvecs']
3349 else:
3350 return res['x']
3353def _minimize_powell(func, x0, args=(), callback=None, bounds=None,
3354 xtol=1e-4, ftol=1e-4, maxiter=None, maxfev=None,
3355 disp=False, direc=None, return_all=False,
3356 **unknown_options):
3357 """
3358 Minimization of scalar function of one or more variables using the
3359 modified Powell algorithm.
3361 Parameters
3362 ----------
3363 fun : callable
3364 The objective function to be minimized.
3366 ``fun(x, *args) -> float``
3368 where ``x`` is a 1-D array with shape (n,) and ``args``
3369 is a tuple of the fixed parameters needed to completely
3370 specify the function.
3371 x0 : ndarray, shape (n,)
3372 Initial guess. Array of real elements of size (n,),
3373 where ``n`` is the number of independent variables.
3374 args : tuple, optional
3375 Extra arguments passed to the objective function and its
3376 derivatives (`fun`, `jac` and `hess` functions).
3377 method : str or callable, optional
3378 The present documentation is specific to ``method='powell'``, but other
3379 options are available. See documentation for `scipy.optimize.minimize`.
3380 bounds : sequence or `Bounds`, optional
3381 Bounds on decision variables. There are two ways to specify the bounds:
3383 1. Instance of `Bounds` class.
3384 2. Sequence of ``(min, max)`` pairs for each element in `x`. None
3385 is used to specify no bound.
3387 If bounds are not provided, then an unbounded line search will be used.
3388 If bounds are provided and the initial guess is within the bounds, then
3389 every function evaluation throughout the minimization procedure will be
3390 within the bounds. If bounds are provided, the initial guess is outside
3391 the bounds, and `direc` is full rank (or left to default), then some
3392 function evaluations during the first iteration may be outside the
3393 bounds, but every function evaluation after the first iteration will be
3394 within the bounds. If `direc` is not full rank, then some parameters
3395 may not be optimized and the solution is not guaranteed to be within
3396 the bounds.
3398 options : dict, optional
3399 A dictionary of solver options. All methods accept the following
3400 generic options:
3402 maxiter : int
3403 Maximum number of iterations to perform. Depending on the
3404 method each iteration may use several function evaluations.
3405 disp : bool
3406 Set to True to print convergence messages.
3408 See method-specific options for ``method='powell'`` below.
3409 callback : callable, optional
3410 Called after each iteration. The signature is:
3412 ``callback(xk)``
3414 where ``xk`` is the current parameter vector.
3416 Returns
3417 -------
3418 res : OptimizeResult
3419 The optimization result represented as a ``OptimizeResult`` object.
3420 Important attributes are: ``x`` the solution array, ``success`` a
3421 Boolean flag indicating if the optimizer exited successfully and
3422 ``message`` which describes the cause of the termination. See
3423 `OptimizeResult` for a description of other attributes.
3425 Options
3426 -------
3427 disp : bool
3428 Set to True to print convergence messages.
3429 xtol : float
3430 Relative error in solution `xopt` acceptable for convergence.
3431 ftol : float
3432 Relative error in ``fun(xopt)`` acceptable for convergence.
3433 maxiter, maxfev : int
3434 Maximum allowed number of iterations and function evaluations.
3435 Will default to ``N*1000``, where ``N`` is the number of
3436 variables, if neither `maxiter` or `maxfev` is set. If both
3437 `maxiter` and `maxfev` are set, minimization will stop at the
3438 first reached.
3439 direc : ndarray
3440 Initial set of direction vectors for the Powell method.
3441 return_all : bool, optional
3442 Set to True to return a list of the best solution at each of the
3443 iterations.
3444 """
3445 _check_unknown_options(unknown_options)
3446 maxfun = maxfev
3447 retall = return_all
3449 x = asarray(x0).flatten()
3450 if retall:
3451 allvecs = [x]
3452 N = len(x)
3453 # If neither are set, then set both to default
3454 if maxiter is None and maxfun is None:
3455 maxiter = N * 1000
3456 maxfun = N * 1000
3457 elif maxiter is None:
3458 # Convert remaining Nones, to np.inf, unless the other is np.inf, in
3459 # which case use the default to avoid unbounded iteration
3460 if maxfun == np.inf:
3461 maxiter = N * 1000
3462 else:
3463 maxiter = np.inf
3464 elif maxfun is None:
3465 if maxiter == np.inf:
3466 maxfun = N * 1000
3467 else:
3468 maxfun = np.inf
3470 # we need to use a mutable object here that we can update in the
3471 # wrapper function
3472 fcalls, func = _wrap_scalar_function_maxfun_validation(func, args, maxfun)
3474 if direc is None:
3475 direc = eye(N, dtype=float)
3476 else:
3477 direc = asarray(direc, dtype=float)
3478 if np.linalg.matrix_rank(direc) != direc.shape[0]:
3479 warnings.warn("direc input is not full rank, some parameters may "
3480 "not be optimized",
3481 OptimizeWarning, 3)
3483 if bounds is None:
3484 # don't make these arrays of all +/- inf. because
3485 # _linesearch_powell will do an unnecessary check of all the elements.
3486 # just keep them None, _linesearch_powell will not have to check
3487 # all the elements.
3488 lower_bound, upper_bound = None, None
3489 else:
3490 # bounds is standardized in _minimize.py.
3491 lower_bound, upper_bound = bounds.lb, bounds.ub
3492 if np.any(lower_bound > x0) or np.any(x0 > upper_bound):
3493 warnings.warn("Initial guess is not within the specified bounds",
3494 OptimizeWarning, 3)
3496 fval = squeeze(func(x))
3497 x1 = x.copy()
3498 iter = 0
3499 while True:
3500 try:
3501 fx = fval
3502 bigind = 0
3503 delta = 0.0
3504 for i in range(N):
3505 direc1 = direc[i]
3506 fx2 = fval
3507 fval, x, direc1 = _linesearch_powell(func, x, direc1,
3508 tol=xtol * 100,
3509 lower_bound=lower_bound,
3510 upper_bound=upper_bound,
3511 fval=fval)
3512 if (fx2 - fval) > delta:
3513 delta = fx2 - fval
3514 bigind = i
3515 iter += 1
3516 if retall:
3517 allvecs.append(x)
3518 intermediate_result = OptimizeResult(x=x, fun=fval)
3519 if _call_callback_maybe_halt(callback, intermediate_result):
3520 break
3521 bnd = ftol * (np.abs(fx) + np.abs(fval)) + 1e-20
3522 if 2.0 * (fx - fval) <= bnd:
3523 break
3524 if fcalls[0] >= maxfun:
3525 break
3526 if iter >= maxiter:
3527 break
3528 if np.isnan(fx) and np.isnan(fval):
3529 # Ended up in a nan-region: bail out
3530 break
3532 # Construct the extrapolated point
3533 direc1 = x - x1
3534 x1 = x.copy()
3535 # make sure that we don't go outside the bounds when extrapolating
3536 if lower_bound is None and upper_bound is None:
3537 lmax = 1
3538 else:
3539 _, lmax = _line_for_search(x, direc1, lower_bound, upper_bound)
3540 x2 = x + min(lmax, 1) * direc1
3541 fx2 = squeeze(func(x2))
3543 if (fx > fx2):
3544 t = 2.0*(fx + fx2 - 2.0*fval)
3545 temp = (fx - fval - delta)
3546 t *= temp*temp
3547 temp = fx - fx2
3548 t -= delta*temp*temp
3549 if t < 0.0:
3550 fval, x, direc1 = _linesearch_powell(
3551 func, x, direc1,
3552 tol=xtol * 100,
3553 lower_bound=lower_bound,
3554 upper_bound=upper_bound,
3555 fval=fval
3556 )
3557 if np.any(direc1):
3558 direc[bigind] = direc[-1]
3559 direc[-1] = direc1
3560 except _MaxFuncCallError:
3561 break
3563 warnflag = 0
3564 msg = _status_message['success']
3565 # out of bounds is more urgent than exceeding function evals or iters,
3566 # but I don't want to cause inconsistencies by changing the
3567 # established warning flags for maxfev and maxiter, so the out of bounds
3568 # warning flag becomes 3, but is checked for first.
3569 if bounds and (np.any(lower_bound > x) or np.any(x > upper_bound)):
3570 warnflag = 4
3571 msg = _status_message['out_of_bounds']
3572 elif fcalls[0] >= maxfun:
3573 warnflag = 1
3574 msg = _status_message['maxfev']
3575 elif iter >= maxiter:
3576 warnflag = 2
3577 msg = _status_message['maxiter']
3578 elif np.isnan(fval) or np.isnan(x).any():
3579 warnflag = 3
3580 msg = _status_message['nan']
3582 if disp:
3583 _print_success_message_or_warn(warnflag, msg, RuntimeWarning)
3584 print(" Current function value: %f" % fval)
3585 print(" Iterations: %d" % iter)
3586 print(" Function evaluations: %d" % fcalls[0])
3588 result = OptimizeResult(fun=fval, direc=direc, nit=iter, nfev=fcalls[0],
3589 status=warnflag, success=(warnflag == 0),
3590 message=msg, x=x)
3591 if retall:
3592 result['allvecs'] = allvecs
3593 return result
3596def _endprint(x, flag, fval, maxfun, xtol, disp):
3597 if flag == 0:
3598 if disp > 1:
3599 print("\nOptimization terminated successfully;\n"
3600 "The returned value satisfies the termination criteria\n"
3601 "(using xtol = ", xtol, ")")
3602 return
3604 if flag == 1:
3605 msg = ("\nMaximum number of function evaluations exceeded --- "
3606 "increase maxfun argument.\n")
3607 elif flag == 2:
3608 msg = "\n{}".format(_status_message['nan'])
3610 _print_success_message_or_warn(flag, msg)
3611 return
3614def brute(func, ranges, args=(), Ns=20, full_output=0, finish=fmin,
3615 disp=False, workers=1):
3616 """Minimize a function over a given range by brute force.
3618 Uses the "brute force" method, i.e., computes the function's value
3619 at each point of a multidimensional grid of points, to find the global
3620 minimum of the function.
3622 The function is evaluated everywhere in the range with the datatype of the
3623 first call to the function, as enforced by the ``vectorize`` NumPy
3624 function. The value and type of the function evaluation returned when
3625 ``full_output=True`` are affected in addition by the ``finish`` argument
3626 (see Notes).
3628 The brute force approach is inefficient because the number of grid points
3629 increases exponentially - the number of grid points to evaluate is
3630 ``Ns ** len(x)``. Consequently, even with coarse grid spacing, even
3631 moderately sized problems can take a long time to run, and/or run into
3632 memory limitations.
3634 Parameters
3635 ----------
3636 func : callable
3637 The objective function to be minimized. Must be in the
3638 form ``f(x, *args)``, where ``x`` is the argument in
3639 the form of a 1-D array and ``args`` is a tuple of any
3640 additional fixed parameters needed to completely specify
3641 the function.
3642 ranges : tuple
3643 Each component of the `ranges` tuple must be either a
3644 "slice object" or a range tuple of the form ``(low, high)``.
3645 The program uses these to create the grid of points on which
3646 the objective function will be computed. See `Note 2` for
3647 more detail.
3648 args : tuple, optional
3649 Any additional fixed parameters needed to completely specify
3650 the function.
3651 Ns : int, optional
3652 Number of grid points along the axes, if not otherwise
3653 specified. See `Note2`.
3654 full_output : bool, optional
3655 If True, return the evaluation grid and the objective function's
3656 values on it.
3657 finish : callable, optional
3658 An optimization function that is called with the result of brute force
3659 minimization as initial guess. `finish` should take `func` and
3660 the initial guess as positional arguments, and take `args` as
3661 keyword arguments. It may additionally take `full_output`
3662 and/or `disp` as keyword arguments. Use None if no "polishing"
3663 function is to be used. See Notes for more details.
3664 disp : bool, optional
3665 Set to True to print convergence messages from the `finish` callable.
3666 workers : int or map-like callable, optional
3667 If `workers` is an int the grid is subdivided into `workers`
3668 sections and evaluated in parallel (uses
3669 `multiprocessing.Pool <multiprocessing>`).
3670 Supply `-1` to use all cores available to the Process.
3671 Alternatively supply a map-like callable, such as
3672 `multiprocessing.Pool.map` for evaluating the grid in parallel.
3673 This evaluation is carried out as ``workers(func, iterable)``.
3674 Requires that `func` be pickleable.
3676 .. versionadded:: 1.3.0
3678 Returns
3679 -------
3680 x0 : ndarray
3681 A 1-D array containing the coordinates of a point at which the
3682 objective function had its minimum value. (See `Note 1` for
3683 which point is returned.)
3684 fval : float
3685 Function value at the point `x0`. (Returned when `full_output` is
3686 True.)
3687 grid : tuple
3688 Representation of the evaluation grid. It has the same
3689 length as `x0`. (Returned when `full_output` is True.)
3690 Jout : ndarray
3691 Function values at each point of the evaluation
3692 grid, i.e., ``Jout = func(*grid)``. (Returned
3693 when `full_output` is True.)
3695 See Also
3696 --------
3697 basinhopping, differential_evolution
3699 Notes
3700 -----
3701 *Note 1*: The program finds the gridpoint at which the lowest value
3702 of the objective function occurs. If `finish` is None, that is the
3703 point returned. When the global minimum occurs within (or not very far
3704 outside) the grid's boundaries, and the grid is fine enough, that
3705 point will be in the neighborhood of the global minimum.
3707 However, users often employ some other optimization program to
3708 "polish" the gridpoint values, i.e., to seek a more precise
3709 (local) minimum near `brute's` best gridpoint.
3710 The `brute` function's `finish` option provides a convenient way to do
3711 that. Any polishing program used must take `brute's` output as its
3712 initial guess as a positional argument, and take `brute's` input values
3713 for `args` as keyword arguments, otherwise an error will be raised.
3714 It may additionally take `full_output` and/or `disp` as keyword arguments.
3716 `brute` assumes that the `finish` function returns either an
3717 `OptimizeResult` object or a tuple in the form:
3718 ``(xmin, Jmin, ... , statuscode)``, where ``xmin`` is the minimizing
3719 value of the argument, ``Jmin`` is the minimum value of the objective
3720 function, "..." may be some other returned values (which are not used
3721 by `brute`), and ``statuscode`` is the status code of the `finish` program.
3723 Note that when `finish` is not None, the values returned are those
3724 of the `finish` program, *not* the gridpoint ones. Consequently,
3725 while `brute` confines its search to the input grid points,
3726 the `finish` program's results usually will not coincide with any
3727 gridpoint, and may fall outside the grid's boundary. Thus, if a
3728 minimum only needs to be found over the provided grid points, make
3729 sure to pass in `finish=None`.
3731 *Note 2*: The grid of points is a `numpy.mgrid` object.
3732 For `brute` the `ranges` and `Ns` inputs have the following effect.
3733 Each component of the `ranges` tuple can be either a slice object or a
3734 two-tuple giving a range of values, such as (0, 5). If the component is a
3735 slice object, `brute` uses it directly. If the component is a two-tuple
3736 range, `brute` internally converts it to a slice object that interpolates
3737 `Ns` points from its low-value to its high-value, inclusive.
3739 Examples
3740 --------
3741 We illustrate the use of `brute` to seek the global minimum of a function
3742 of two variables that is given as the sum of a positive-definite
3743 quadratic and two deep "Gaussian-shaped" craters. Specifically, define
3744 the objective function `f` as the sum of three other functions,
3745 ``f = f1 + f2 + f3``. We suppose each of these has a signature
3746 ``(z, *params)``, where ``z = (x, y)``, and ``params`` and the functions
3747 are as defined below.
3749 >>> import numpy as np
3750 >>> params = (2, 3, 7, 8, 9, 10, 44, -1, 2, 26, 1, -2, 0.5)
3751 >>> def f1(z, *params):
3752 ... x, y = z
3753 ... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
3754 ... return (a * x**2 + b * x * y + c * y**2 + d*x + e*y + f)
3756 >>> def f2(z, *params):
3757 ... x, y = z
3758 ... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
3759 ... return (-g*np.exp(-((x-h)**2 + (y-i)**2) / scale))
3761 >>> def f3(z, *params):
3762 ... x, y = z
3763 ... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
3764 ... return (-j*np.exp(-((x-k)**2 + (y-l)**2) / scale))
3766 >>> def f(z, *params):
3767 ... return f1(z, *params) + f2(z, *params) + f3(z, *params)
3769 Thus, the objective function may have local minima near the minimum
3770 of each of the three functions of which it is composed. To
3771 use `fmin` to polish its gridpoint result, we may then continue as
3772 follows:
3774 >>> rranges = (slice(-4, 4, 0.25), slice(-4, 4, 0.25))
3775 >>> from scipy import optimize
3776 >>> resbrute = optimize.brute(f, rranges, args=params, full_output=True,
3777 ... finish=optimize.fmin)
3778 >>> resbrute[0] # global minimum
3779 array([-1.05665192, 1.80834843])
3780 >>> resbrute[1] # function value at global minimum
3781 -3.4085818767
3783 Note that if `finish` had been set to None, we would have gotten the
3784 gridpoint [-1.0 1.75] where the rounded function value is -2.892.
3786 """
3787 N = len(ranges)
3788 if N > 40:
3789 raise ValueError("Brute Force not possible with more "
3790 "than 40 variables.")
3791 lrange = list(ranges)
3792 for k in range(N):
3793 if not isinstance(lrange[k], slice):
3794 if len(lrange[k]) < 3:
3795 lrange[k] = tuple(lrange[k]) + (complex(Ns),)
3796 lrange[k] = slice(*lrange[k])
3797 if (N == 1):
3798 lrange = lrange[0]
3800 grid = np.mgrid[lrange]
3802 # obtain an array of parameters that is iterable by a map-like callable
3803 inpt_shape = grid.shape
3804 if (N > 1):
3805 grid = np.reshape(grid, (inpt_shape[0], np.prod(inpt_shape[1:]))).T
3807 if not np.iterable(args):
3808 args = (args,)
3810 wrapped_func = _Brute_Wrapper(func, args)
3812 # iterate over input arrays, possibly in parallel
3813 with MapWrapper(pool=workers) as mapper:
3814 Jout = np.array(list(mapper(wrapped_func, grid)))
3815 if (N == 1):
3816 grid = (grid,)
3817 Jout = np.squeeze(Jout)
3818 elif (N > 1):
3819 Jout = np.reshape(Jout, inpt_shape[1:])
3820 grid = np.reshape(grid.T, inpt_shape)
3822 Nshape = shape(Jout)
3824 indx = argmin(Jout.ravel(), axis=-1)
3825 Nindx = np.empty(N, int)
3826 xmin = np.empty(N, float)
3827 for k in range(N - 1, -1, -1):
3828 thisN = Nshape[k]
3829 Nindx[k] = indx % Nshape[k]
3830 indx = indx // thisN
3831 for k in range(N):
3832 xmin[k] = grid[k][tuple(Nindx)]
3834 Jmin = Jout[tuple(Nindx)]
3835 if (N == 1):
3836 grid = grid[0]
3837 xmin = xmin[0]
3839 if callable(finish):
3840 # set up kwargs for `finish` function
3841 finish_args = _getfullargspec(finish).args
3842 finish_kwargs = dict()
3843 if 'full_output' in finish_args:
3844 finish_kwargs['full_output'] = 1
3845 if 'disp' in finish_args:
3846 finish_kwargs['disp'] = disp
3847 elif 'options' in finish_args:
3848 # pass 'disp' as `options`
3849 # (e.g., if `finish` is `minimize`)
3850 finish_kwargs['options'] = {'disp': disp}
3852 # run minimizer
3853 res = finish(func, xmin, args=args, **finish_kwargs)
3855 if isinstance(res, OptimizeResult):
3856 xmin = res.x
3857 Jmin = res.fun
3858 success = res.success
3859 else:
3860 xmin = res[0]
3861 Jmin = res[1]
3862 success = res[-1] == 0
3863 if not success:
3864 if disp:
3865 warnings.warn(
3866 "Either final optimization did not succeed "
3867 "or `finish` does not return `statuscode` as its last "
3868 "argument.", RuntimeWarning, 2)
3870 if full_output:
3871 return xmin, Jmin, grid, Jout
3872 else:
3873 return xmin
3876class _Brute_Wrapper:
3877 """
3878 Object to wrap user cost function for optimize.brute, allowing picklability
3879 """
3881 def __init__(self, f, args):
3882 self.f = f
3883 self.args = [] if args is None else args
3885 def __call__(self, x):
3886 # flatten needed for one dimensional case.
3887 return self.f(np.asarray(x).flatten(), *self.args)
3890def show_options(solver=None, method=None, disp=True):
3891 """
3892 Show documentation for additional options of optimization solvers.
3894 These are method-specific options that can be supplied through the
3895 ``options`` dict.
3897 Parameters
3898 ----------
3899 solver : str
3900 Type of optimization solver. One of 'minimize', 'minimize_scalar',
3901 'root', 'root_scalar', 'linprog', or 'quadratic_assignment'.
3902 method : str, optional
3903 If not given, shows all methods of the specified solver. Otherwise,
3904 show only the options for the specified method. Valid values
3905 corresponds to methods' names of respective solver (e.g., 'BFGS' for
3906 'minimize').
3907 disp : bool, optional
3908 Whether to print the result rather than returning it.
3910 Returns
3911 -------
3912 text
3913 Either None (for disp=True) or the text string (disp=False)
3915 Notes
3916 -----
3917 The solver-specific methods are:
3919 `scipy.optimize.minimize`
3921 - :ref:`Nelder-Mead <optimize.minimize-neldermead>`
3922 - :ref:`Powell <optimize.minimize-powell>`
3923 - :ref:`CG <optimize.minimize-cg>`
3924 - :ref:`BFGS <optimize.minimize-bfgs>`
3925 - :ref:`Newton-CG <optimize.minimize-newtoncg>`
3926 - :ref:`L-BFGS-B <optimize.minimize-lbfgsb>`
3927 - :ref:`TNC <optimize.minimize-tnc>`
3928 - :ref:`COBYLA <optimize.minimize-cobyla>`
3929 - :ref:`SLSQP <optimize.minimize-slsqp>`
3930 - :ref:`dogleg <optimize.minimize-dogleg>`
3931 - :ref:`trust-ncg <optimize.minimize-trustncg>`
3933 `scipy.optimize.root`
3935 - :ref:`hybr <optimize.root-hybr>`
3936 - :ref:`lm <optimize.root-lm>`
3937 - :ref:`broyden1 <optimize.root-broyden1>`
3938 - :ref:`broyden2 <optimize.root-broyden2>`
3939 - :ref:`anderson <optimize.root-anderson>`
3940 - :ref:`linearmixing <optimize.root-linearmixing>`
3941 - :ref:`diagbroyden <optimize.root-diagbroyden>`
3942 - :ref:`excitingmixing <optimize.root-excitingmixing>`
3943 - :ref:`krylov <optimize.root-krylov>`
3944 - :ref:`df-sane <optimize.root-dfsane>`
3946 `scipy.optimize.minimize_scalar`
3948 - :ref:`brent <optimize.minimize_scalar-brent>`
3949 - :ref:`golden <optimize.minimize_scalar-golden>`
3950 - :ref:`bounded <optimize.minimize_scalar-bounded>`
3952 `scipy.optimize.root_scalar`
3954 - :ref:`bisect <optimize.root_scalar-bisect>`
3955 - :ref:`brentq <optimize.root_scalar-brentq>`
3956 - :ref:`brenth <optimize.root_scalar-brenth>`
3957 - :ref:`ridder <optimize.root_scalar-ridder>`
3958 - :ref:`toms748 <optimize.root_scalar-toms748>`
3959 - :ref:`newton <optimize.root_scalar-newton>`
3960 - :ref:`secant <optimize.root_scalar-secant>`
3961 - :ref:`halley <optimize.root_scalar-halley>`
3963 `scipy.optimize.linprog`
3965 - :ref:`simplex <optimize.linprog-simplex>`
3966 - :ref:`interior-point <optimize.linprog-interior-point>`
3967 - :ref:`revised simplex <optimize.linprog-revised_simplex>`
3968 - :ref:`highs <optimize.linprog-highs>`
3969 - :ref:`highs-ds <optimize.linprog-highs-ds>`
3970 - :ref:`highs-ipm <optimize.linprog-highs-ipm>`
3972 `scipy.optimize.quadratic_assignment`
3974 - :ref:`faq <optimize.qap-faq>`
3975 - :ref:`2opt <optimize.qap-2opt>`
3977 Examples
3978 --------
3979 We can print documentations of a solver in stdout:
3981 >>> from scipy.optimize import show_options
3982 >>> show_options(solver="minimize")
3983 ...
3985 Specifying a method is possible:
3987 >>> show_options(solver="minimize", method="Nelder-Mead")
3988 ...
3990 We can also get the documentations as a string:
3992 >>> show_options(solver="minimize", method="Nelder-Mead", disp=False)
3993 Minimization of scalar function of one or more variables using the ...
3995 """
3996 import textwrap
3998 doc_routines = {
3999 'minimize': (
4000 ('bfgs', 'scipy.optimize._optimize._minimize_bfgs'),
4001 ('cg', 'scipy.optimize._optimize._minimize_cg'),
4002 ('cobyla', 'scipy.optimize._cobyla_py._minimize_cobyla'),
4003 ('dogleg', 'scipy.optimize._trustregion_dogleg._minimize_dogleg'),
4004 ('l-bfgs-b', 'scipy.optimize._lbfgsb_py._minimize_lbfgsb'),
4005 ('nelder-mead', 'scipy.optimize._optimize._minimize_neldermead'),
4006 ('newton-cg', 'scipy.optimize._optimize._minimize_newtoncg'),
4007 ('powell', 'scipy.optimize._optimize._minimize_powell'),
4008 ('slsqp', 'scipy.optimize._slsqp_py._minimize_slsqp'),
4009 ('tnc', 'scipy.optimize._tnc._minimize_tnc'),
4010 ('trust-ncg',
4011 'scipy.optimize._trustregion_ncg._minimize_trust_ncg'),
4012 ('trust-constr',
4013 'scipy.optimize._trustregion_constr.'
4014 '_minimize_trustregion_constr'),
4015 ('trust-exact',
4016 'scipy.optimize._trustregion_exact._minimize_trustregion_exact'),
4017 ('trust-krylov',
4018 'scipy.optimize._trustregion_krylov._minimize_trust_krylov'),
4019 ),
4020 'root': (
4021 ('hybr', 'scipy.optimize._minpack_py._root_hybr'),
4022 ('lm', 'scipy.optimize._root._root_leastsq'),
4023 ('broyden1', 'scipy.optimize._root._root_broyden1_doc'),
4024 ('broyden2', 'scipy.optimize._root._root_broyden2_doc'),
4025 ('anderson', 'scipy.optimize._root._root_anderson_doc'),
4026 ('diagbroyden', 'scipy.optimize._root._root_diagbroyden_doc'),
4027 ('excitingmixing', 'scipy.optimize._root._root_excitingmixing_doc'),
4028 ('linearmixing', 'scipy.optimize._root._root_linearmixing_doc'),
4029 ('krylov', 'scipy.optimize._root._root_krylov_doc'),
4030 ('df-sane', 'scipy.optimize._spectral._root_df_sane'),
4031 ),
4032 'root_scalar': (
4033 ('bisect', 'scipy.optimize._root_scalar._root_scalar_bisect_doc'),
4034 ('brentq', 'scipy.optimize._root_scalar._root_scalar_brentq_doc'),
4035 ('brenth', 'scipy.optimize._root_scalar._root_scalar_brenth_doc'),
4036 ('ridder', 'scipy.optimize._root_scalar._root_scalar_ridder_doc'),
4037 ('toms748', 'scipy.optimize._root_scalar._root_scalar_toms748_doc'),
4038 ('secant', 'scipy.optimize._root_scalar._root_scalar_secant_doc'),
4039 ('newton', 'scipy.optimize._root_scalar._root_scalar_newton_doc'),
4040 ('halley', 'scipy.optimize._root_scalar._root_scalar_halley_doc'),
4041 ),
4042 'linprog': (
4043 ('simplex', 'scipy.optimize._linprog._linprog_simplex_doc'),
4044 ('interior-point', 'scipy.optimize._linprog._linprog_ip_doc'),
4045 ('revised simplex', 'scipy.optimize._linprog._linprog_rs_doc'),
4046 ('highs-ipm', 'scipy.optimize._linprog._linprog_highs_ipm_doc'),
4047 ('highs-ds', 'scipy.optimize._linprog._linprog_highs_ds_doc'),
4048 ('highs', 'scipy.optimize._linprog._linprog_highs_doc'),
4049 ),
4050 'quadratic_assignment': (
4051 ('faq', 'scipy.optimize._qap._quadratic_assignment_faq'),
4052 ('2opt', 'scipy.optimize._qap._quadratic_assignment_2opt'),
4053 ),
4054 'minimize_scalar': (
4055 ('brent', 'scipy.optimize._optimize._minimize_scalar_brent'),
4056 ('bounded', 'scipy.optimize._optimize._minimize_scalar_bounded'),
4057 ('golden', 'scipy.optimize._optimize._minimize_scalar_golden'),
4058 ),
4059 }
4061 if solver is None:
4062 text = ["\n\n\n========\n", "minimize\n", "========\n"]
4063 text.append(show_options('minimize', disp=False))
4064 text.extend(["\n\n===============\n", "minimize_scalar\n",
4065 "===============\n"])
4066 text.append(show_options('minimize_scalar', disp=False))
4067 text.extend(["\n\n\n====\n", "root\n",
4068 "====\n"])
4069 text.append(show_options('root', disp=False))
4070 text.extend(['\n\n\n=======\n', 'linprog\n',
4071 '=======\n'])
4072 text.append(show_options('linprog', disp=False))
4073 text = "".join(text)
4074 else:
4075 solver = solver.lower()
4076 if solver not in doc_routines:
4077 raise ValueError(f'Unknown solver {solver!r}')
4079 if method is None:
4080 text = []
4081 for name, _ in doc_routines[solver]:
4082 text.extend(["\n\n" + name, "\n" + "="*len(name) + "\n\n"])
4083 text.append(show_options(solver, name, disp=False))
4084 text = "".join(text)
4085 else:
4086 method = method.lower()
4087 methods = dict(doc_routines[solver])
4088 if method not in methods:
4089 raise ValueError(f"Unknown method {method!r}")
4090 name = methods[method]
4092 # Import function object
4093 parts = name.split('.')
4094 mod_name = ".".join(parts[:-1])
4095 __import__(mod_name)
4096 obj = getattr(sys.modules[mod_name], parts[-1])
4098 # Get doc
4099 doc = obj.__doc__
4100 if doc is not None:
4101 text = textwrap.dedent(doc).strip()
4102 else:
4103 text = ""
4105 if disp:
4106 print(text)
4107 return
4108 else:
4109 return text