Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/optimize/optimize.py : 6%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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
30from numpy import (atleast_1d, eye, argmin, zeros, shape, squeeze,
31 asarray, sqrt, Inf, asfarray, isinf)
32import numpy as np
33from .linesearch import (line_search_wolfe1, line_search_wolfe2,
34 line_search_wolfe2 as line_search,
35 LineSearchWarning)
36from ._numdiff import approx_derivative
37from scipy._lib._util import getfullargspec_no_self as _getfullargspec
38from scipy._lib._util import MapWrapper
39from scipy.optimize._differentiable_functions import ScalarFunction, FD_METHODS
42# standard status messages of optimizers
43_status_message = {'success': 'Optimization terminated successfully.',
44 'maxfev': 'Maximum number of function evaluations has '
45 'been exceeded.',
46 'maxiter': 'Maximum number of iterations has been '
47 'exceeded.',
48 'pr_loss': 'Desired error not necessarily achieved due '
49 'to precision loss.',
50 'nan': 'NaN result encountered.',
51 'out_of_bounds': 'The result is outside of the provided '
52 'bounds.'}
55class MemoizeJac(object):
56 """ Decorator that caches the return values of a function returning `(fun, grad)`
57 each time it is called. """
59 def __init__(self, fun):
60 self.fun = fun
61 self.jac = None
62 self._value = None
63 self.x = None
65 def _compute_if_needed(self, x, *args):
66 if not np.all(x == self.x) or self._value is None or self.jac is None:
67 self.x = np.asarray(x).copy()
68 fg = self.fun(x, *args)
69 self.jac = fg[1]
70 self._value = fg[0]
72 def __call__(self, x, *args):
73 """ returns the the function value """
74 self._compute_if_needed(x, *args)
75 return self._value
77 def derivative(self, x, *args):
78 self._compute_if_needed(x, *args)
79 return self.jac
82class OptimizeResult(dict):
83 """ Represents the optimization result.
85 Attributes
86 ----------
87 x : ndarray
88 The solution of the optimization.
89 success : bool
90 Whether or not the optimizer exited successfully.
91 status : int
92 Termination status of the optimizer. Its value depends on the
93 underlying solver. Refer to `message` for details.
94 message : str
95 Description of the cause of the termination.
96 fun, jac, hess: ndarray
97 Values of objective function, its Jacobian and its Hessian (if
98 available). The Hessians may be approximations, see the documentation
99 of the function in question.
100 hess_inv : object
101 Inverse of the objective function's Hessian; may be an approximation.
102 Not available for all solvers. The type of this attribute may be
103 either np.ndarray or scipy.sparse.linalg.LinearOperator.
104 nfev, njev, nhev : int
105 Number of evaluations of the objective functions and of its
106 Jacobian and Hessian.
107 nit : int
108 Number of iterations performed by the optimizer.
109 maxcv : float
110 The maximum constraint violation.
112 Notes
113 -----
114 There may be additional attributes not listed above depending of the
115 specific solver. Since this class is essentially a subclass of dict
116 with attribute accessors, one can see which attributes are available
117 using the `keys()` method.
118 """
120 def __getattr__(self, name):
121 try:
122 return self[name]
123 except KeyError:
124 raise AttributeError(name)
126 __setattr__ = dict.__setitem__
127 __delattr__ = dict.__delitem__
129 def __repr__(self):
130 if self.keys():
131 m = max(map(len, list(self.keys()))) + 1
132 return '\n'.join([k.rjust(m) + ': ' + repr(v)
133 for k, v in sorted(self.items())])
134 else:
135 return self.__class__.__name__ + "()"
137 def __dir__(self):
138 return list(self.keys())
141class OptimizeWarning(UserWarning):
142 pass
145def _check_unknown_options(unknown_options):
146 if unknown_options:
147 msg = ", ".join(map(str, unknown_options.keys()))
148 # Stack level 4: this is called from _minimize_*, which is
149 # called from another function in SciPy. Level 4 is the first
150 # level in user code.
151 warnings.warn("Unknown solver options: %s" % msg, OptimizeWarning, 4)
154def is_array_scalar(x):
155 """Test whether `x` is either a scalar or an array scalar.
157 """
158 return np.size(x) == 1
161_epsilon = sqrt(np.finfo(float).eps)
164def vecnorm(x, ord=2):
165 if ord == Inf:
166 return np.amax(np.abs(x))
167 elif ord == -Inf:
168 return np.amin(np.abs(x))
169 else:
170 return np.sum(np.abs(x)**ord, axis=0)**(1.0 / ord)
173def _prepare_scalar_function(fun, x0, jac=None, args=(), bounds=None,
174 epsilon=None, finite_diff_rel_step=None,
175 hess=None):
176 """
177 Creates a ScalarFunction object for use with scalar minimizers
178 (BFGS/LBFGSB/SLSQP/TNC/CG/etc).
180 Parameters
181 ----------
182 fun : callable
183 The objective function to be minimized.
185 ``fun(x, *args) -> float``
187 where ``x`` is an 1-D array with shape (n,) and ``args``
188 is a tuple of the fixed parameters needed to completely
189 specify the function.
190 x0 : ndarray, shape (n,)
191 Initial guess. Array of real elements of size (n,),
192 where 'n' is the number of independent variables.
193 jac : {callable, '2-point', '3-point', 'cs', None}, optional
194 Method for computing the gradient vector. If it is a callable, it
195 should be a function that returns the gradient vector:
197 ``jac(x, *args) -> array_like, shape (n,)``
199 If one of `{'2-point', '3-point', 'cs'}` is selected then the gradient
200 is calculated with a relative step for finite differences. If `None`,
201 then two-point finite differences with an absolute step is used.
202 args : tuple, optional
203 Extra arguments passed to the objective function and its
204 derivatives (`fun`, `jac` functions).
205 bounds : sequence, optional
206 Bounds on variables. 'new-style' bounds are required.
207 eps : float or ndarray
208 If `jac is None` the absolute step size used for numerical
209 approximation of the jacobian via forward differences.
210 finite_diff_rel_step : None or array_like, optional
211 If `jac in ['2-point', '3-point', 'cs']` the relative step size to
212 use for numerical approximation of the jacobian. The absolute step
213 size is computed as ``h = rel_step * sign(x0) * max(1, abs(x0))``,
214 possibly adjusted to fit into the bounds. For ``method='3-point'``
215 the sign of `h` is ignored. If None (default) then step is selected
216 automatically.
217 hess : {callable, '2-point', '3-point', 'cs', None}
218 Computes the Hessian matrix. If it is callable, it should return the
219 Hessian matrix:
221 ``hess(x, *args) -> {LinearOperator, spmatrix, array}, (n, n)``
223 Alternatively, the keywords {'2-point', '3-point', 'cs'} select a
224 finite difference scheme for numerical estimation.
225 Whenever the gradient is estimated via finite-differences, the Hessian
226 cannot be estimated with options {'2-point', '3-point', 'cs'} and needs
227 to be estimated using one of the quasi-Newton strategies.
229 Returns
230 -------
231 sf : ScalarFunction
232 """
233 if callable(jac):
234 grad = jac
235 elif jac in FD_METHODS:
236 # epsilon is set to None so that ScalarFunction is made to use
237 # rel_step
238 epsilon = None
239 grad = jac
240 else:
241 # default (jac is None) is to do 2-point finite differences with
242 # absolute step size. ScalarFunction has to be provided an
243 # epsilon value that is not None to use absolute steps. This is
244 # normally the case from most _minimize* methods.
245 grad = '2-point'
246 epsilon = epsilon
248 if hess is None:
249 # ScalarFunction requires something for hess, so we give a dummy
250 # implementation here if nothing is provided, return a value of None
251 # so that downstream minimisers halt. The results of `fun.hess`
252 # should not be used.
253 def hess(x, *args):
254 return None
256 if bounds is None:
257 bounds = (-np.inf, np.inf)
259 # ScalarFunction caches. Reuse of fun(x) during grad
260 # calculation reduces overall function evaluations.
261 sf = ScalarFunction(fun, x0, args, grad, hess,
262 finite_diff_rel_step, bounds, epsilon=epsilon)
264 return sf
267def rosen(x):
268 """
269 The Rosenbrock function.
271 The function computed is::
273 sum(100.0*(x[1:] - x[:-1]**2.0)**2.0 + (1 - x[:-1])**2.0)
275 Parameters
276 ----------
277 x : array_like
278 1-D array of points at which the Rosenbrock function is to be computed.
280 Returns
281 -------
282 f : float
283 The value of the Rosenbrock function.
285 See Also
286 --------
287 rosen_der, rosen_hess, rosen_hess_prod
289 Examples
290 --------
291 >>> from scipy.optimize import rosen
292 >>> X = 0.1 * np.arange(10)
293 >>> rosen(X)
294 76.56
296 """
297 x = asarray(x)
298 r = np.sum(100.0 * (x[1:] - x[:-1]**2.0)**2.0 + (1 - x[:-1])**2.0,
299 axis=0)
300 return r
303def rosen_der(x):
304 """
305 The derivative (i.e. gradient) of the Rosenbrock function.
307 Parameters
308 ----------
309 x : array_like
310 1-D array of points at which the derivative is to be computed.
312 Returns
313 -------
314 rosen_der : (N,) ndarray
315 The gradient of the Rosenbrock function at `x`.
317 See Also
318 --------
319 rosen, rosen_hess, rosen_hess_prod
321 Examples
322 --------
323 >>> from scipy.optimize import rosen_der
324 >>> X = 0.1 * np.arange(9)
325 >>> rosen_der(X)
326 array([ -2. , 10.6, 15.6, 13.4, 6.4, -3. , -12.4, -19.4, 62. ])
328 """
329 x = asarray(x)
330 xm = x[1:-1]
331 xm_m1 = x[:-2]
332 xm_p1 = x[2:]
333 der = np.zeros_like(x)
334 der[1:-1] = (200 * (xm - xm_m1**2) -
335 400 * (xm_p1 - xm**2) * xm - 2 * (1 - xm))
336 der[0] = -400 * x[0] * (x[1] - x[0]**2) - 2 * (1 - x[0])
337 der[-1] = 200 * (x[-1] - x[-2]**2)
338 return der
341def rosen_hess(x):
342 """
343 The Hessian matrix of the Rosenbrock function.
345 Parameters
346 ----------
347 x : array_like
348 1-D array of points at which the Hessian matrix is to be computed.
350 Returns
351 -------
352 rosen_hess : ndarray
353 The Hessian matrix of the Rosenbrock function at `x`.
355 See Also
356 --------
357 rosen, rosen_der, rosen_hess_prod
359 Examples
360 --------
361 >>> from scipy.optimize import rosen_hess
362 >>> X = 0.1 * np.arange(4)
363 >>> rosen_hess(X)
364 array([[-38., 0., 0., 0.],
365 [ 0., 134., -40., 0.],
366 [ 0., -40., 130., -80.],
367 [ 0., 0., -80., 200.]])
369 """
370 x = atleast_1d(x)
371 H = np.diag(-400 * x[:-1], 1) - np.diag(400 * x[:-1], -1)
372 diagonal = np.zeros(len(x), dtype=x.dtype)
373 diagonal[0] = 1200 * x[0]**2 - 400 * x[1] + 2
374 diagonal[-1] = 200
375 diagonal[1:-1] = 202 + 1200 * x[1:-1]**2 - 400 * x[2:]
376 H = H + np.diag(diagonal)
377 return H
380def rosen_hess_prod(x, p):
381 """
382 Product of the Hessian matrix of the Rosenbrock function with a vector.
384 Parameters
385 ----------
386 x : array_like
387 1-D array of points at which the Hessian matrix is to be computed.
388 p : array_like
389 1-D array, the vector to be multiplied by the Hessian matrix.
391 Returns
392 -------
393 rosen_hess_prod : ndarray
394 The Hessian matrix of the Rosenbrock function at `x` multiplied
395 by the vector `p`.
397 See Also
398 --------
399 rosen, rosen_der, rosen_hess
401 Examples
402 --------
403 >>> from scipy.optimize import rosen_hess_prod
404 >>> X = 0.1 * np.arange(9)
405 >>> p = 0.5 * np.arange(9)
406 >>> rosen_hess_prod(X, p)
407 array([ -0., 27., -10., -95., -192., -265., -278., -195., -180.])
409 """
410 x = atleast_1d(x)
411 Hp = np.zeros(len(x), dtype=x.dtype)
412 Hp[0] = (1200 * x[0]**2 - 400 * x[1] + 2) * p[0] - 400 * x[0] * p[1]
413 Hp[1:-1] = (-400 * x[:-2] * p[:-2] +
414 (202 + 1200 * x[1:-1]**2 - 400 * x[2:]) * p[1:-1] -
415 400 * x[1:-1] * p[2:])
416 Hp[-1] = -400 * x[-2] * p[-2] + 200*p[-1]
417 return Hp
420def wrap_function(function, args):
421 ncalls = [0]
422 if function is None:
423 return ncalls, None
425 def function_wrapper(*wrapper_args):
426 ncalls[0] += 1
427 return function(*(wrapper_args + args))
429 return ncalls, function_wrapper
432def fmin(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, maxfun=None,
433 full_output=0, disp=1, retall=0, callback=None, initial_simplex=None):
434 """
435 Minimize a function using the downhill simplex algorithm.
437 This algorithm only uses function values, not derivatives or second
438 derivatives.
440 Parameters
441 ----------
442 func : callable func(x,*args)
443 The objective function to be minimized.
444 x0 : ndarray
445 Initial guess.
446 args : tuple, optional
447 Extra arguments passed to func, i.e., ``f(x,*args)``.
448 xtol : float, optional
449 Absolute error in xopt between iterations that is acceptable for
450 convergence.
451 ftol : number, optional
452 Absolute error in func(xopt) between iterations that is acceptable for
453 convergence.
454 maxiter : int, optional
455 Maximum number of iterations to perform.
456 maxfun : number, optional
457 Maximum number of function evaluations to make.
458 full_output : bool, optional
459 Set to True if fopt and warnflag outputs are desired.
460 disp : bool, optional
461 Set to True to print convergence messages.
462 retall : bool, optional
463 Set to True to return list of solutions at each iteration.
464 callback : callable, optional
465 Called after each iteration, as callback(xk), where xk is the
466 current parameter vector.
467 initial_simplex : array_like of shape (N + 1, N), optional
468 Initial simplex. If given, overrides `x0`.
469 ``initial_simplex[j,:]`` should contain the coordinates of
470 the jth vertex of the ``N+1`` vertices in the simplex, where
471 ``N`` is the dimension.
473 Returns
474 -------
475 xopt : ndarray
476 Parameter that minimizes function.
477 fopt : float
478 Value of function at minimum: ``fopt = func(xopt)``.
479 iter : int
480 Number of iterations performed.
481 funcalls : int
482 Number of function calls made.
483 warnflag : int
484 1 : Maximum number of function evaluations made.
485 2 : Maximum number of iterations reached.
486 allvecs : list
487 Solution at each iteration.
489 See also
490 --------
491 minimize: Interface to minimization algorithms for multivariate
492 functions. See the 'Nelder-Mead' `method` in particular.
494 Notes
495 -----
496 Uses a Nelder-Mead simplex algorithm to find the minimum of function of
497 one or more variables.
499 This algorithm has a long history of successful use in applications.
500 But it will usually be slower than an algorithm that uses first or
501 second derivative information. In practice, it can have poor
502 performance in high-dimensional problems and is not robust to
503 minimizing complicated functions. Additionally, there currently is no
504 complete theory describing when the algorithm will successfully
505 converge to the minimum, or how fast it will if it does. Both the ftol and
506 xtol criteria must be met for convergence.
508 Examples
509 --------
510 >>> def f(x):
511 ... return x**2
513 >>> from scipy import optimize
515 >>> minimum = optimize.fmin(f, 1)
516 Optimization terminated successfully.
517 Current function value: 0.000000
518 Iterations: 17
519 Function evaluations: 34
520 >>> minimum[0]
521 -8.8817841970012523e-16
523 References
524 ----------
525 .. [1] Nelder, J.A. and Mead, R. (1965), "A simplex method for function
526 minimization", The Computer Journal, 7, pp. 308-313
528 .. [2] Wright, M.H. (1996), "Direct Search Methods: Once Scorned, Now
529 Respectable", in Numerical Analysis 1995, Proceedings of the
530 1995 Dundee Biennial Conference in Numerical Analysis, D.F.
531 Griffiths and G.A. Watson (Eds.), Addison Wesley Longman,
532 Harlow, UK, pp. 191-208.
534 """
535 opts = {'xatol': xtol,
536 'fatol': ftol,
537 'maxiter': maxiter,
538 'maxfev': maxfun,
539 'disp': disp,
540 'return_all': retall,
541 'initial_simplex': initial_simplex}
543 res = _minimize_neldermead(func, x0, args, callback=callback, **opts)
544 if full_output:
545 retlist = res['x'], res['fun'], res['nit'], res['nfev'], res['status']
546 if retall:
547 retlist += (res['allvecs'], )
548 return retlist
549 else:
550 if retall:
551 return res['x'], res['allvecs']
552 else:
553 return res['x']
556def _minimize_neldermead(func, x0, args=(), callback=None,
557 maxiter=None, maxfev=None, disp=False,
558 return_all=False, initial_simplex=None,
559 xatol=1e-4, fatol=1e-4, adaptive=False,
560 **unknown_options):
561 """
562 Minimization of scalar function of one or more variables using the
563 Nelder-Mead algorithm.
565 Options
566 -------
567 disp : bool
568 Set to True to print convergence messages.
569 maxiter, maxfev : int
570 Maximum allowed number of iterations and function evaluations.
571 Will default to ``N*200``, where ``N`` is the number of
572 variables, if neither `maxiter` or `maxfev` is set. If both
573 `maxiter` and `maxfev` are set, minimization will stop at the
574 first reached.
575 return_all : bool, optional
576 Set to True to return a list of the best solution at each of the
577 iterations.
578 initial_simplex : array_like of shape (N + 1, N)
579 Initial simplex. If given, overrides `x0`.
580 ``initial_simplex[j,:]`` should contain the coordinates of
581 the jth vertex of the ``N+1`` vertices in the simplex, where
582 ``N`` is the dimension.
583 xatol : float, optional
584 Absolute error in xopt between iterations that is acceptable for
585 convergence.
586 fatol : number, optional
587 Absolute error in func(xopt) between iterations that is acceptable for
588 convergence.
589 adaptive : bool, optional
590 Adapt algorithm parameters to dimensionality of problem. Useful for
591 high-dimensional minimization [1]_.
593 References
594 ----------
595 .. [1] Gao, F. and Han, L.
596 Implementing the Nelder-Mead simplex algorithm with adaptive
597 parameters. 2012. Computational Optimization and Applications.
598 51:1, pp. 259-277
600 """
601 if 'ftol' in unknown_options:
602 warnings.warn("ftol is deprecated for Nelder-Mead,"
603 " use fatol instead. If you specified both, only"
604 " fatol is used.",
605 DeprecationWarning)
606 if (np.isclose(fatol, 1e-4) and
607 not np.isclose(unknown_options['ftol'], 1e-4)):
608 # only ftol was probably specified, use it.
609 fatol = unknown_options['ftol']
610 unknown_options.pop('ftol')
611 if 'xtol' in unknown_options:
612 warnings.warn("xtol is deprecated for Nelder-Mead,"
613 " use xatol instead. If you specified both, only"
614 " xatol is used.",
615 DeprecationWarning)
616 if (np.isclose(xatol, 1e-4) and
617 not np.isclose(unknown_options['xtol'], 1e-4)):
618 # only xtol was probably specified, use it.
619 xatol = unknown_options['xtol']
620 unknown_options.pop('xtol')
622 _check_unknown_options(unknown_options)
623 maxfun = maxfev
624 retall = return_all
626 fcalls, func = wrap_function(func, args)
628 if adaptive:
629 dim = float(len(x0))
630 rho = 1
631 chi = 1 + 2/dim
632 psi = 0.75 - 1/(2*dim)
633 sigma = 1 - 1/dim
634 else:
635 rho = 1
636 chi = 2
637 psi = 0.5
638 sigma = 0.5
640 nonzdelt = 0.05
641 zdelt = 0.00025
643 x0 = asfarray(x0).flatten()
645 if initial_simplex is None:
646 N = len(x0)
648 sim = np.zeros((N + 1, N), dtype=x0.dtype)
649 sim[0] = x0
650 for k in range(N):
651 y = np.array(x0, copy=True)
652 if y[k] != 0:
653 y[k] = (1 + nonzdelt)*y[k]
654 else:
655 y[k] = zdelt
656 sim[k + 1] = y
657 else:
658 sim = np.asfarray(initial_simplex).copy()
659 if sim.ndim != 2 or sim.shape[0] != sim.shape[1] + 1:
660 raise ValueError("`initial_simplex` should be an array of shape (N+1,N)")
661 if len(x0) != sim.shape[1]:
662 raise ValueError("Size of `initial_simplex` is not consistent with `x0`")
663 N = sim.shape[1]
665 if retall:
666 allvecs = [sim[0]]
668 # If neither are set, then set both to default
669 if maxiter is None and maxfun is None:
670 maxiter = N * 200
671 maxfun = N * 200
672 elif maxiter is None:
673 # Convert remaining Nones, to np.inf, unless the other is np.inf, in
674 # which case use the default to avoid unbounded iteration
675 if maxfun == np.inf:
676 maxiter = N * 200
677 else:
678 maxiter = np.inf
679 elif maxfun is None:
680 if maxiter == np.inf:
681 maxfun = N * 200
682 else:
683 maxfun = np.inf
685 one2np1 = list(range(1, N + 1))
686 fsim = np.zeros((N + 1,), float)
688 for k in range(N + 1):
689 fsim[k] = func(sim[k])
691 ind = np.argsort(fsim)
692 fsim = np.take(fsim, ind, 0)
693 # sort so sim[0,:] has the lowest function value
694 sim = np.take(sim, ind, 0)
696 iterations = 1
698 while (fcalls[0] < maxfun and iterations < maxiter):
699 if (np.max(np.ravel(np.abs(sim[1:] - sim[0]))) <= xatol and
700 np.max(np.abs(fsim[0] - fsim[1:])) <= fatol):
701 break
703 xbar = np.add.reduce(sim[:-1], 0) / N
704 xr = (1 + rho) * xbar - rho * sim[-1]
705 fxr = func(xr)
706 doshrink = 0
708 if fxr < fsim[0]:
709 xe = (1 + rho * chi) * xbar - rho * chi * sim[-1]
710 fxe = func(xe)
712 if fxe < fxr:
713 sim[-1] = xe
714 fsim[-1] = fxe
715 else:
716 sim[-1] = xr
717 fsim[-1] = fxr
718 else: # fsim[0] <= fxr
719 if fxr < fsim[-2]:
720 sim[-1] = xr
721 fsim[-1] = fxr
722 else: # fxr >= fsim[-2]
723 # Perform contraction
724 if fxr < fsim[-1]:
725 xc = (1 + psi * rho) * xbar - psi * rho * sim[-1]
726 fxc = func(xc)
728 if fxc <= fxr:
729 sim[-1] = xc
730 fsim[-1] = fxc
731 else:
732 doshrink = 1
733 else:
734 # Perform an inside contraction
735 xcc = (1 - psi) * xbar + psi * sim[-1]
736 fxcc = func(xcc)
738 if fxcc < fsim[-1]:
739 sim[-1] = xcc
740 fsim[-1] = fxcc
741 else:
742 doshrink = 1
744 if doshrink:
745 for j in one2np1:
746 sim[j] = sim[0] + sigma * (sim[j] - sim[0])
747 fsim[j] = func(sim[j])
749 ind = np.argsort(fsim)
750 sim = np.take(sim, ind, 0)
751 fsim = np.take(fsim, ind, 0)
752 if callback is not None:
753 callback(sim[0])
754 iterations += 1
755 if retall:
756 allvecs.append(sim[0])
758 x = sim[0]
759 fval = np.min(fsim)
760 warnflag = 0
762 if fcalls[0] >= maxfun:
763 warnflag = 1
764 msg = _status_message['maxfev']
765 if disp:
766 print('Warning: ' + msg)
767 elif iterations >= maxiter:
768 warnflag = 2
769 msg = _status_message['maxiter']
770 if disp:
771 print('Warning: ' + msg)
772 else:
773 msg = _status_message['success']
774 if disp:
775 print(msg)
776 print(" Current function value: %f" % fval)
777 print(" Iterations: %d" % iterations)
778 print(" Function evaluations: %d" % fcalls[0])
780 result = OptimizeResult(fun=fval, nit=iterations, nfev=fcalls[0],
781 status=warnflag, success=(warnflag == 0),
782 message=msg, x=x, final_simplex=(sim, fsim))
783 if retall:
784 result['allvecs'] = allvecs
785 return result
788def approx_fprime(xk, f, epsilon, *args):
789 """Finite-difference approximation of the gradient of a scalar function.
791 Parameters
792 ----------
793 xk : array_like
794 The coordinate vector at which to determine the gradient of `f`.
795 f : callable
796 The function of which to determine the gradient (partial derivatives).
797 Should take `xk` as first argument, other arguments to `f` can be
798 supplied in ``*args``. Should return a scalar, the value of the
799 function at `xk`.
800 epsilon : array_like
801 Increment to `xk` to use for determining the function gradient.
802 If a scalar, uses the same finite difference delta for all partial
803 derivatives. If an array, should contain one value per element of
804 `xk`.
805 \\*args : args, optional
806 Any other arguments that are to be passed to `f`.
808 Returns
809 -------
810 grad : ndarray
811 The partial derivatives of `f` to `xk`.
813 See Also
814 --------
815 check_grad : Check correctness of gradient function against approx_fprime.
817 Notes
818 -----
819 The function gradient is determined by the forward finite difference
820 formula::
822 f(xk[i] + epsilon[i]) - f(xk[i])
823 f'[i] = ---------------------------------
824 epsilon[i]
826 The main use of `approx_fprime` is in scalar function optimizers like
827 `fmin_bfgs`, to determine numerically the Jacobian of a function.
829 Examples
830 --------
831 >>> from scipy import optimize
832 >>> def func(x, c0, c1):
833 ... "Coordinate vector `x` should be an array of size two."
834 ... return c0 * x[0]**2 + c1*x[1]**2
836 >>> x = np.ones(2)
837 >>> c0, c1 = (1, 200)
838 >>> eps = np.sqrt(np.finfo(float).eps)
839 >>> optimize.approx_fprime(x, func, [eps, np.sqrt(200) * eps], c0, c1)
840 array([ 2. , 400.00004198])
842 """
843 xk = np.asarray(xk, float)
845 f0 = f(xk, *args)
846 if not np.isscalar(f0):
847 try:
848 f0 = f0.item()
849 except (ValueError, AttributeError):
850 raise ValueError("The user-provided "
851 "objective function must "
852 "return a scalar value.")
854 return approx_derivative(f, xk, method='2-point', abs_step=epsilon,
855 args=args, f0=f0)
858def check_grad(func, grad, x0, *args, **kwargs):
859 """Check the correctness of a gradient function by comparing it against a
860 (forward) finite-difference approximation of the gradient.
862 Parameters
863 ----------
864 func : callable ``func(x0, *args)``
865 Function whose derivative is to be checked.
866 grad : callable ``grad(x0, *args)``
867 Gradient of `func`.
868 x0 : ndarray
869 Points to check `grad` against forward difference approximation of grad
870 using `func`.
871 args : \\*args, optional
872 Extra arguments passed to `func` and `grad`.
873 epsilon : float, optional
874 Step size used for the finite difference approximation. It defaults to
875 ``sqrt(np.finfo(float).eps)``, which is approximately 1.49e-08.
877 Returns
878 -------
879 err : float
880 The square root of the sum of squares (i.e., the 2-norm) of the
881 difference between ``grad(x0, *args)`` and the finite difference
882 approximation of `grad` using func at the points `x0`.
884 See Also
885 --------
886 approx_fprime
888 Examples
889 --------
890 >>> def func(x):
891 ... return x[0]**2 - 0.5 * x[1]**3
892 >>> def grad(x):
893 ... return [2 * x[0], -1.5 * x[1]**2]
894 >>> from scipy.optimize import check_grad
895 >>> check_grad(func, grad, [1.5, -1.5])
896 2.9802322387695312e-08
898 """
899 step = kwargs.pop('epsilon', _epsilon)
900 if kwargs:
901 raise ValueError("Unknown keyword arguments: %r" %
902 (list(kwargs.keys()),))
903 return sqrt(sum((grad(x0, *args) -
904 approx_fprime(x0, func, step, *args))**2))
907def approx_fhess_p(x0, p, fprime, epsilon, *args):
908 # calculate fprime(x0) first, as this may be cached by ScalarFunction
909 f1 = fprime(*((x0,) + args))
910 f2 = fprime(*((x0 + epsilon*p,) + args))
911 return (f2 - f1) / epsilon
914class _LineSearchError(RuntimeError):
915 pass
918def _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval,
919 **kwargs):
920 """
921 Same as line_search_wolfe1, but fall back to line_search_wolfe2 if
922 suitable step length is not found, and raise an exception if a
923 suitable step length is not found.
925 Raises
926 ------
927 _LineSearchError
928 If no suitable step size is found
930 """
932 extra_condition = kwargs.pop('extra_condition', None)
934 ret = line_search_wolfe1(f, fprime, xk, pk, gfk,
935 old_fval, old_old_fval,
936 **kwargs)
938 if ret[0] is not None and extra_condition is not None:
939 xp1 = xk + ret[0] * pk
940 if not extra_condition(ret[0], xp1, ret[3], ret[5]):
941 # Reject step if extra_condition fails
942 ret = (None,)
944 if ret[0] is None:
945 # line search failed: try different one.
946 with warnings.catch_warnings():
947 warnings.simplefilter('ignore', LineSearchWarning)
948 kwargs2 = {}
949 for key in ('c1', 'c2', 'amax'):
950 if key in kwargs:
951 kwargs2[key] = kwargs[key]
952 ret = line_search_wolfe2(f, fprime, xk, pk, gfk,
953 old_fval, old_old_fval,
954 extra_condition=extra_condition,
955 **kwargs2)
957 if ret[0] is None:
958 raise _LineSearchError()
960 return ret
963def fmin_bfgs(f, x0, fprime=None, args=(), gtol=1e-5, norm=Inf,
964 epsilon=_epsilon, maxiter=None, full_output=0, disp=1,
965 retall=0, callback=None):
966 """
967 Minimize a function using the BFGS algorithm.
969 Parameters
970 ----------
971 f : callable f(x,*args)
972 Objective function to be minimized.
973 x0 : ndarray
974 Initial guess.
975 fprime : callable f'(x,*args), optional
976 Gradient of f.
977 args : tuple, optional
978 Extra arguments passed to f and fprime.
979 gtol : float, optional
980 Gradient norm must be less than gtol before successful termination.
981 norm : float, optional
982 Order of norm (Inf is max, -Inf is min)
983 epsilon : int or ndarray, optional
984 If fprime is approximated, use this value for the step size.
985 callback : callable, optional
986 An optional user-supplied function to call after each
987 iteration. Called as callback(xk), where xk is the
988 current parameter vector.
989 maxiter : int, optional
990 Maximum number of iterations to perform.
991 full_output : bool, optional
992 If True,return fopt, func_calls, grad_calls, and warnflag
993 in addition to xopt.
994 disp : bool, optional
995 Print convergence message if True.
996 retall : bool, optional
997 Return a list of results at each iteration if True.
999 Returns
1000 -------
1001 xopt : ndarray
1002 Parameters which minimize f, i.e., f(xopt) == fopt.
1003 fopt : float
1004 Minimum value.
1005 gopt : ndarray
1006 Value of gradient at minimum, f'(xopt), which should be near 0.
1007 Bopt : ndarray
1008 Value of 1/f''(xopt), i.e., the inverse Hessian matrix.
1009 func_calls : int
1010 Number of function_calls made.
1011 grad_calls : int
1012 Number of gradient calls made.
1013 warnflag : integer
1014 1 : Maximum number of iterations exceeded.
1015 2 : Gradient and/or function calls not changing.
1016 3 : NaN result encountered.
1017 allvecs : list
1018 The value of xopt at each iteration. Only returned if retall is True.
1020 See also
1021 --------
1022 minimize: Interface to minimization algorithms for multivariate
1023 functions. See the 'BFGS' `method` in particular.
1025 Notes
1026 -----
1027 Optimize the function, f, whose gradient is given by fprime
1028 using the quasi-Newton method of Broyden, Fletcher, Goldfarb,
1029 and Shanno (BFGS)
1031 References
1032 ----------
1033 Wright, and Nocedal 'Numerical Optimization', 1999, p. 198.
1035 """
1036 opts = {'gtol': gtol,
1037 'norm': norm,
1038 'eps': epsilon,
1039 'disp': disp,
1040 'maxiter': maxiter,
1041 'return_all': retall}
1043 res = _minimize_bfgs(f, x0, args, fprime, callback=callback, **opts)
1045 if full_output:
1046 retlist = (res['x'], res['fun'], res['jac'], res['hess_inv'],
1047 res['nfev'], res['njev'], res['status'])
1048 if retall:
1049 retlist += (res['allvecs'], )
1050 return retlist
1051 else:
1052 if retall:
1053 return res['x'], res['allvecs']
1054 else:
1055 return res['x']
1058def _minimize_bfgs(fun, x0, args=(), jac=None, callback=None,
1059 gtol=1e-5, norm=Inf, eps=_epsilon, maxiter=None,
1060 disp=False, return_all=False, finite_diff_rel_step=None,
1061 **unknown_options):
1062 """
1063 Minimization of scalar function of one or more variables using the
1064 BFGS algorithm.
1066 Options
1067 -------
1068 disp : bool
1069 Set to True to print convergence messages.
1070 maxiter : int
1071 Maximum number of iterations to perform.
1072 gtol : float
1073 Gradient norm must be less than `gtol` before successful
1074 termination.
1075 norm : float
1076 Order of norm (Inf is max, -Inf is min).
1077 eps : float or ndarray
1078 If `jac is None` the absolute step size used for numerical
1079 approximation of the jacobian via forward differences.
1080 return_all : bool, optional
1081 Set to True to return a list of the best solution at each of the
1082 iterations.
1083 finite_diff_rel_step : None or array_like, optional
1084 If `jac in ['2-point', '3-point', 'cs']` the relative step size to
1085 use for numerical approximation of the jacobian. The absolute step
1086 size is computed as ``h = rel_step * sign(x0) * max(1, abs(x0))``,
1087 possibly adjusted to fit into the bounds. For ``method='3-point'``
1088 the sign of `h` is ignored. If None (default) then step is selected
1089 automatically.
1091 """
1092 _check_unknown_options(unknown_options)
1093 retall = return_all
1095 x0 = asarray(x0).flatten()
1096 if x0.ndim == 0:
1097 x0.shape = (1,)
1098 if maxiter is None:
1099 maxiter = len(x0) * 200
1101 sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,
1102 finite_diff_rel_step=finite_diff_rel_step)
1104 f = sf.fun
1105 myfprime = sf.grad
1107 old_fval = f(x0)
1108 gfk = myfprime(x0)
1110 if not np.isscalar(old_fval):
1111 try:
1112 old_fval = old_fval.item()
1113 except (ValueError, AttributeError):
1114 raise ValueError("The user-provided "
1115 "objective function must "
1116 "return a scalar value.")
1118 k = 0
1119 N = len(x0)
1120 I = np.eye(N, dtype=int)
1121 Hk = I
1123 # Sets the initial step guess to dx ~ 1
1124 old_old_fval = old_fval + np.linalg.norm(gfk) / 2
1126 xk = x0
1127 if retall:
1128 allvecs = [x0]
1129 warnflag = 0
1130 gnorm = vecnorm(gfk, ord=norm)
1131 while (gnorm > gtol) and (k < maxiter):
1132 pk = -np.dot(Hk, gfk)
1133 try:
1134 alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
1135 _line_search_wolfe12(f, myfprime, xk, pk, gfk,
1136 old_fval, old_old_fval, amin=1e-100, amax=1e100)
1137 except _LineSearchError:
1138 # Line search failed to find a better solution.
1139 warnflag = 2
1140 break
1142 xkp1 = xk + alpha_k * pk
1143 if retall:
1144 allvecs.append(xkp1)
1145 sk = xkp1 - xk
1146 xk = xkp1
1147 if gfkp1 is None:
1148 gfkp1 = myfprime(xkp1)
1150 yk = gfkp1 - gfk
1151 gfk = gfkp1
1152 if callback is not None:
1153 callback(xk)
1154 k += 1
1155 gnorm = vecnorm(gfk, ord=norm)
1156 if (gnorm <= gtol):
1157 break
1159 if not np.isfinite(old_fval):
1160 # We correctly found +-Inf as optimal value, or something went
1161 # wrong.
1162 warnflag = 2
1163 break
1165 try: # this was handled in numeric, let it remaines for more safety
1166 rhok = 1.0 / (np.dot(yk, sk))
1167 except ZeroDivisionError:
1168 rhok = 1000.0
1169 if disp:
1170 print("Divide-by-zero encountered: rhok assumed large")
1171 if isinf(rhok): # this is patch for NumPy
1172 rhok = 1000.0
1173 if disp:
1174 print("Divide-by-zero encountered: rhok assumed large")
1175 A1 = I - sk[:, np.newaxis] * yk[np.newaxis, :] * rhok
1176 A2 = I - yk[:, np.newaxis] * sk[np.newaxis, :] * rhok
1177 Hk = np.dot(A1, np.dot(Hk, A2)) + (rhok * sk[:, np.newaxis] *
1178 sk[np.newaxis, :])
1180 fval = old_fval
1182 if warnflag == 2:
1183 msg = _status_message['pr_loss']
1184 elif k >= maxiter:
1185 warnflag = 1
1186 msg = _status_message['maxiter']
1187 elif np.isnan(gnorm) or np.isnan(fval) or np.isnan(xk).any():
1188 warnflag = 3
1189 msg = _status_message['nan']
1190 else:
1191 msg = _status_message['success']
1193 if disp:
1194 print("%s%s" % ("Warning: " if warnflag != 0 else "", msg))
1195 print(" Current function value: %f" % fval)
1196 print(" Iterations: %d" % k)
1197 print(" Function evaluations: %d" % sf.nfev)
1198 print(" Gradient evaluations: %d" % sf.ngev)
1200 result = OptimizeResult(fun=fval, jac=gfk, hess_inv=Hk, nfev=sf.nfev,
1201 njev=sf.ngev, status=warnflag,
1202 success=(warnflag == 0), message=msg, x=xk,
1203 nit=k)
1204 if retall:
1205 result['allvecs'] = allvecs
1206 return result
1209def fmin_cg(f, x0, fprime=None, args=(), gtol=1e-5, norm=Inf, epsilon=_epsilon,
1210 maxiter=None, full_output=0, disp=1, retall=0, callback=None):
1211 """
1212 Minimize a function using a nonlinear conjugate gradient algorithm.
1214 Parameters
1215 ----------
1216 f : callable, ``f(x, *args)``
1217 Objective function to be minimized. Here `x` must be a 1-D array of
1218 the variables that are to be changed in the search for a minimum, and
1219 `args` are the other (fixed) parameters of `f`.
1220 x0 : ndarray
1221 A user-supplied initial estimate of `xopt`, the optimal value of `x`.
1222 It must be a 1-D array of values.
1223 fprime : callable, ``fprime(x, *args)``, optional
1224 A function that returns the gradient of `f` at `x`. Here `x` and `args`
1225 are as described above for `f`. The returned value must be a 1-D array.
1226 Defaults to None, in which case the gradient is approximated
1227 numerically (see `epsilon`, below).
1228 args : tuple, optional
1229 Parameter values passed to `f` and `fprime`. Must be supplied whenever
1230 additional fixed parameters are needed to completely specify the
1231 functions `f` and `fprime`.
1232 gtol : float, optional
1233 Stop when the norm of the gradient is less than `gtol`.
1234 norm : float, optional
1235 Order to use for the norm of the gradient
1236 (``-np.Inf`` is min, ``np.Inf`` is max).
1237 epsilon : float or ndarray, optional
1238 Step size(s) to use when `fprime` is approximated numerically. Can be a
1239 scalar or a 1-D array. Defaults to ``sqrt(eps)``, with eps the
1240 floating point machine precision. Usually ``sqrt(eps)`` is about
1241 1.5e-8.
1242 maxiter : int, optional
1243 Maximum number of iterations to perform. Default is ``200 * len(x0)``.
1244 full_output : bool, optional
1245 If True, return `fopt`, `func_calls`, `grad_calls`, and `warnflag` in
1246 addition to `xopt`. See the Returns section below for additional
1247 information on optional return values.
1248 disp : bool, optional
1249 If True, return a convergence message, followed by `xopt`.
1250 retall : bool, optional
1251 If True, add to the returned values the results of each iteration.
1252 callback : callable, optional
1253 An optional user-supplied function, called after each iteration.
1254 Called as ``callback(xk)``, where ``xk`` is the current value of `x0`.
1256 Returns
1257 -------
1258 xopt : ndarray
1259 Parameters which minimize f, i.e., ``f(xopt) == fopt``.
1260 fopt : float, optional
1261 Minimum value found, f(xopt). Only returned if `full_output` is True.
1262 func_calls : int, optional
1263 The number of function_calls made. Only returned if `full_output`
1264 is True.
1265 grad_calls : int, optional
1266 The number of gradient calls made. Only returned if `full_output` is
1267 True.
1268 warnflag : int, optional
1269 Integer value with warning status, only returned if `full_output` is
1270 True.
1272 0 : Success.
1274 1 : The maximum number of iterations was exceeded.
1276 2 : Gradient and/or function calls were not changing. May indicate
1277 that precision was lost, i.e., the routine did not converge.
1279 3 : NaN result encountered.
1281 allvecs : list of ndarray, optional
1282 List of arrays, containing the results at each iteration.
1283 Only returned if `retall` is True.
1285 See Also
1286 --------
1287 minimize : common interface to all `scipy.optimize` algorithms for
1288 unconstrained and constrained minimization of multivariate
1289 functions. It provides an alternative way to call
1290 ``fmin_cg``, by specifying ``method='CG'``.
1292 Notes
1293 -----
1294 This conjugate gradient algorithm is based on that of Polak and Ribiere
1295 [1]_.
1297 Conjugate gradient methods tend to work better when:
1299 1. `f` has a unique global minimizing point, and no local minima or
1300 other stationary points,
1301 2. `f` is, at least locally, reasonably well approximated by a
1302 quadratic function of the variables,
1303 3. `f` is continuous and has a continuous gradient,
1304 4. `fprime` is not too large, e.g., has a norm less than 1000,
1305 5. The initial guess, `x0`, is reasonably close to `f` 's global
1306 minimizing point, `xopt`.
1308 References
1309 ----------
1310 .. [1] Wright & Nocedal, "Numerical Optimization", 1999, pp. 120-122.
1312 Examples
1313 --------
1314 Example 1: seek the minimum value of the expression
1315 ``a*u**2 + b*u*v + c*v**2 + d*u + e*v + f`` for given values
1316 of the parameters and an initial guess ``(u, v) = (0, 0)``.
1318 >>> args = (2, 3, 7, 8, 9, 10) # parameter values
1319 >>> def f(x, *args):
1320 ... u, v = x
1321 ... a, b, c, d, e, f = args
1322 ... return a*u**2 + b*u*v + c*v**2 + d*u + e*v + f
1323 >>> def gradf(x, *args):
1324 ... u, v = x
1325 ... a, b, c, d, e, f = args
1326 ... gu = 2*a*u + b*v + d # u-component of the gradient
1327 ... gv = b*u + 2*c*v + e # v-component of the gradient
1328 ... return np.asarray((gu, gv))
1329 >>> x0 = np.asarray((0, 0)) # Initial guess.
1330 >>> from scipy import optimize
1331 >>> res1 = optimize.fmin_cg(f, x0, fprime=gradf, args=args)
1332 Optimization terminated successfully.
1333 Current function value: 1.617021
1334 Iterations: 4
1335 Function evaluations: 8
1336 Gradient evaluations: 8
1337 >>> res1
1338 array([-1.80851064, -0.25531915])
1340 Example 2: solve the same problem using the `minimize` function.
1341 (This `myopts` dictionary shows all of the available options,
1342 although in practice only non-default values would be needed.
1343 The returned value will be a dictionary.)
1345 >>> opts = {'maxiter' : None, # default value.
1346 ... 'disp' : True, # non-default value.
1347 ... 'gtol' : 1e-5, # default value.
1348 ... 'norm' : np.inf, # default value.
1349 ... 'eps' : 1.4901161193847656e-08} # default value.
1350 >>> res2 = optimize.minimize(f, x0, jac=gradf, args=args,
1351 ... method='CG', options=opts)
1352 Optimization terminated successfully.
1353 Current function value: 1.617021
1354 Iterations: 4
1355 Function evaluations: 8
1356 Gradient evaluations: 8
1357 >>> res2.x # minimum found
1358 array([-1.80851064, -0.25531915])
1360 """
1361 opts = {'gtol': gtol,
1362 'norm': norm,
1363 'eps': epsilon,
1364 'disp': disp,
1365 'maxiter': maxiter,
1366 'return_all': retall}
1368 res = _minimize_cg(f, x0, args, fprime, callback=callback, **opts)
1370 if full_output:
1371 retlist = res['x'], res['fun'], res['nfev'], res['njev'], res['status']
1372 if retall:
1373 retlist += (res['allvecs'], )
1374 return retlist
1375 else:
1376 if retall:
1377 return res['x'], res['allvecs']
1378 else:
1379 return res['x']
1382def _minimize_cg(fun, x0, args=(), jac=None, callback=None,
1383 gtol=1e-5, norm=Inf, eps=_epsilon, maxiter=None,
1384 disp=False, return_all=False, finite_diff_rel_step=None,
1385 **unknown_options):
1386 """
1387 Minimization of scalar function of one or more variables using the
1388 conjugate gradient algorithm.
1390 Options
1391 -------
1392 disp : bool
1393 Set to True to print convergence messages.
1394 maxiter : int
1395 Maximum number of iterations to perform.
1396 gtol : float
1397 Gradient norm must be less than `gtol` before successful
1398 termination.
1399 norm : float
1400 Order of norm (Inf is max, -Inf is min).
1401 eps : float or ndarray
1402 If `jac is None` the absolute step size used for numerical
1403 approximation of the jacobian via forward differences.
1404 return_all : bool, optional
1405 Set to True to return a list of the best solution at each of the
1406 iterations.
1407 finite_diff_rel_step : None or array_like, optional
1408 If `jac in ['2-point', '3-point', 'cs']` the relative step size to
1409 use for numerical approximation of the jacobian. The absolute step
1410 size is computed as ``h = rel_step * sign(x0) * max(1, abs(x0))``,
1411 possibly adjusted to fit into the bounds. For ``method='3-point'``
1412 the sign of `h` is ignored. If None (default) then step is selected
1413 automatically.
1414 """
1415 _check_unknown_options(unknown_options)
1417 retall = return_all
1419 x0 = asarray(x0).flatten()
1420 if maxiter is None:
1421 maxiter = len(x0) * 200
1423 sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
1424 finite_diff_rel_step=finite_diff_rel_step)
1426 f = sf.fun
1427 myfprime = sf.grad
1429 old_fval = f(x0)
1430 gfk = myfprime(x0)
1432 if not np.isscalar(old_fval):
1433 try:
1434 old_fval = old_fval.item()
1435 except (ValueError, AttributeError):
1436 raise ValueError("The user-provided "
1437 "objective function must "
1438 "return a scalar value.")
1440 k = 0
1441 xk = x0
1442 # Sets the initial step guess to dx ~ 1
1443 old_old_fval = old_fval + np.linalg.norm(gfk) / 2
1445 if retall:
1446 allvecs = [xk]
1447 warnflag = 0
1448 pk = -gfk
1449 gnorm = vecnorm(gfk, ord=norm)
1451 sigma_3 = 0.01
1453 while (gnorm > gtol) and (k < maxiter):
1454 deltak = np.dot(gfk, gfk)
1456 cached_step = [None]
1458 def polak_ribiere_powell_step(alpha, gfkp1=None):
1459 xkp1 = xk + alpha * pk
1460 if gfkp1 is None:
1461 gfkp1 = myfprime(xkp1)
1462 yk = gfkp1 - gfk
1463 beta_k = max(0, np.dot(yk, gfkp1) / deltak)
1464 pkp1 = -gfkp1 + beta_k * pk
1465 gnorm = vecnorm(gfkp1, ord=norm)
1466 return (alpha, xkp1, pkp1, gfkp1, gnorm)
1468 def descent_condition(alpha, xkp1, fp1, gfkp1):
1469 # Polak-Ribiere+ needs an explicit check of a sufficient
1470 # descent condition, which is not guaranteed by strong Wolfe.
1471 #
1472 # See Gilbert & Nocedal, "Global convergence properties of
1473 # conjugate gradient methods for optimization",
1474 # SIAM J. Optimization 2, 21 (1992).
1475 cached_step[:] = polak_ribiere_powell_step(alpha, gfkp1)
1476 alpha, xk, pk, gfk, gnorm = cached_step
1478 # Accept step if it leads to convergence.
1479 if gnorm <= gtol:
1480 return True
1482 # Accept step if sufficient descent condition applies.
1483 return np.dot(pk, gfk) <= -sigma_3 * np.dot(gfk, gfk)
1485 try:
1486 alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
1487 _line_search_wolfe12(f, myfprime, xk, pk, gfk, old_fval,
1488 old_old_fval, c2=0.4, amin=1e-100, amax=1e100,
1489 extra_condition=descent_condition)
1490 except _LineSearchError:
1491 # Line search failed to find a better solution.
1492 warnflag = 2
1493 break
1495 # Reuse already computed results if possible
1496 if alpha_k == cached_step[0]:
1497 alpha_k, xk, pk, gfk, gnorm = cached_step
1498 else:
1499 alpha_k, xk, pk, gfk, gnorm = polak_ribiere_powell_step(alpha_k, gfkp1)
1501 if retall:
1502 allvecs.append(xk)
1503 if callback is not None:
1504 callback(xk)
1505 k += 1
1507 fval = old_fval
1508 if warnflag == 2:
1509 msg = _status_message['pr_loss']
1510 elif k >= maxiter:
1511 warnflag = 1
1512 msg = _status_message['maxiter']
1513 elif np.isnan(gnorm) or np.isnan(fval) or np.isnan(xk).any():
1514 warnflag = 3
1515 msg = _status_message['nan']
1516 else:
1517 msg = _status_message['success']
1519 if disp:
1520 print("%s%s" % ("Warning: " if warnflag != 0 else "", msg))
1521 print(" Current function value: %f" % fval)
1522 print(" Iterations: %d" % k)
1523 print(" Function evaluations: %d" % sf.nfev)
1524 print(" Gradient evaluations: %d" % sf.ngev)
1526 result = OptimizeResult(fun=fval, jac=gfk, nfev=sf.nfev,
1527 njev=sf.ngev, status=warnflag,
1528 success=(warnflag == 0), message=msg, x=xk,
1529 nit=k)
1530 if retall:
1531 result['allvecs'] = allvecs
1532 return result
1535def fmin_ncg(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5,
1536 epsilon=_epsilon, maxiter=None, full_output=0, disp=1, retall=0,
1537 callback=None):
1538 """
1539 Unconstrained minimization of a function using the Newton-CG method.
1541 Parameters
1542 ----------
1543 f : callable ``f(x, *args)``
1544 Objective function to be minimized.
1545 x0 : ndarray
1546 Initial guess.
1547 fprime : callable ``f'(x, *args)``
1548 Gradient of f.
1549 fhess_p : callable ``fhess_p(x, p, *args)``, optional
1550 Function which computes the Hessian of f times an
1551 arbitrary vector, p.
1552 fhess : callable ``fhess(x, *args)``, optional
1553 Function to compute the Hessian matrix of f.
1554 args : tuple, optional
1555 Extra arguments passed to f, fprime, fhess_p, and fhess
1556 (the same set of extra arguments is supplied to all of
1557 these functions).
1558 epsilon : float or ndarray, optional
1559 If fhess is approximated, use this value for the step size.
1560 callback : callable, optional
1561 An optional user-supplied function which is called after
1562 each iteration. Called as callback(xk), where xk is the
1563 current parameter vector.
1564 avextol : float, optional
1565 Convergence is assumed when the average relative error in
1566 the minimizer falls below this amount.
1567 maxiter : int, optional
1568 Maximum number of iterations to perform.
1569 full_output : bool, optional
1570 If True, return the optional outputs.
1571 disp : bool, optional
1572 If True, print convergence message.
1573 retall : bool, optional
1574 If True, return a list of results at each iteration.
1576 Returns
1577 -------
1578 xopt : ndarray
1579 Parameters which minimize f, i.e., ``f(xopt) == fopt``.
1580 fopt : float
1581 Value of the function at xopt, i.e., ``fopt = f(xopt)``.
1582 fcalls : int
1583 Number of function calls made.
1584 gcalls : int
1585 Number of gradient calls made.
1586 hcalls : int
1587 Number of Hessian calls made.
1588 warnflag : int
1589 Warnings generated by the algorithm.
1590 1 : Maximum number of iterations exceeded.
1591 2 : Line search failure (precision loss).
1592 3 : NaN result encountered.
1593 allvecs : list
1594 The result at each iteration, if retall is True (see below).
1596 See also
1597 --------
1598 minimize: Interface to minimization algorithms for multivariate
1599 functions. See the 'Newton-CG' `method` in particular.
1601 Notes
1602 -----
1603 Only one of `fhess_p` or `fhess` need to be given. If `fhess`
1604 is provided, then `fhess_p` will be ignored. If neither `fhess`
1605 nor `fhess_p` is provided, then the hessian product will be
1606 approximated using finite differences on `fprime`. `fhess_p`
1607 must compute the hessian times an arbitrary vector. If it is not
1608 given, finite-differences on `fprime` are used to compute
1609 it.
1611 Newton-CG methods are also called truncated Newton methods. This
1612 function differs from scipy.optimize.fmin_tnc because
1614 1. scipy.optimize.fmin_ncg is written purely in Python using NumPy
1615 and scipy while scipy.optimize.fmin_tnc calls a C function.
1616 2. scipy.optimize.fmin_ncg is only for unconstrained minimization
1617 while scipy.optimize.fmin_tnc is for unconstrained minimization
1618 or box constrained minimization. (Box constraints give
1619 lower and upper bounds for each variable separately.)
1621 References
1622 ----------
1623 Wright & Nocedal, 'Numerical Optimization', 1999, p. 140.
1625 """
1626 opts = {'xtol': avextol,
1627 'eps': epsilon,
1628 'maxiter': maxiter,
1629 'disp': disp,
1630 'return_all': retall}
1632 res = _minimize_newtoncg(f, x0, args, fprime, fhess, fhess_p,
1633 callback=callback, **opts)
1635 if full_output:
1636 retlist = (res['x'], res['fun'], res['nfev'], res['njev'],
1637 res['nhev'], res['status'])
1638 if retall:
1639 retlist += (res['allvecs'], )
1640 return retlist
1641 else:
1642 if retall:
1643 return res['x'], res['allvecs']
1644 else:
1645 return res['x']
1648def _minimize_newtoncg(fun, x0, args=(), jac=None, hess=None, hessp=None,
1649 callback=None, xtol=1e-5, eps=_epsilon, maxiter=None,
1650 disp=False, return_all=False,
1651 **unknown_options):
1652 """
1653 Minimization of scalar function of one or more variables using the
1654 Newton-CG algorithm.
1656 Note that the `jac` parameter (Jacobian) is required.
1658 Options
1659 -------
1660 disp : bool
1661 Set to True to print convergence messages.
1662 xtol : float
1663 Average relative error in solution `xopt` acceptable for
1664 convergence.
1665 maxiter : int
1666 Maximum number of iterations to perform.
1667 eps : float or ndarray
1668 If `hessp` is approximated, use this value for the step size.
1669 return_all : bool, optional
1670 Set to True to return a list of the best solution at each of the
1671 iterations.
1672 """
1673 _check_unknown_options(unknown_options)
1674 if jac is None:
1675 raise ValueError('Jacobian is required for Newton-CG method')
1676 fhess_p = hessp
1677 fhess = hess
1678 avextol = xtol
1679 epsilon = eps
1680 retall = return_all
1682 x0 = asarray(x0).flatten()
1683 # TODO: allow hess to be approximated by FD?
1684 # TODO: add hessp (callable or FD) to ScalarFunction?
1685 sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps, hess=fhess)
1686 f = sf.fun
1687 fprime = sf.grad
1689 def terminate(warnflag, msg):
1690 if disp:
1691 print(msg)
1692 print(" Current function value: %f" % old_fval)
1693 print(" Iterations: %d" % k)
1694 print(" Function evaluations: %d" % sf.nfev)
1695 print(" Gradient evaluations: %d" % sf.ngev)
1696 print(" Hessian evaluations: %d" % hcalls)
1697 fval = old_fval
1698 result = OptimizeResult(fun=fval, jac=gfk, nfev=sf.nfev,
1699 njev=sf.ngev, nhev=hcalls, status=warnflag,
1700 success=(warnflag == 0), message=msg, x=xk,
1701 nit=k)
1702 if retall:
1703 result['allvecs'] = allvecs
1704 return result
1706 hcalls = 0
1707 if maxiter is None:
1708 maxiter = len(x0)*200
1709 cg_maxiter = 20*len(x0)
1711 xtol = len(x0) * avextol
1712 update = [2 * xtol]
1713 xk = x0
1714 if retall:
1715 allvecs = [xk]
1716 k = 0
1717 gfk = None
1718 old_fval = f(x0)
1719 old_old_fval = None
1720 float64eps = np.finfo(np.float64).eps
1721 while np.add.reduce(np.abs(update)) > xtol:
1722 if k >= maxiter:
1723 msg = "Warning: " + _status_message['maxiter']
1724 return terminate(1, msg)
1725 # Compute a search direction pk by applying the CG method to
1726 # del2 f(xk) p = - grad f(xk) starting from 0.
1727 b = -fprime(xk)
1728 maggrad = np.add.reduce(np.abs(b))
1729 eta = np.min([0.5, np.sqrt(maggrad)])
1730 termcond = eta * maggrad
1731 xsupi = zeros(len(x0), dtype=x0.dtype)
1732 ri = -b
1733 psupi = -ri
1734 i = 0
1735 dri0 = np.dot(ri, ri)
1737 if fhess is not None: # you want to compute hessian once.
1738 A = sf.hess(xk)
1739 hcalls = hcalls + 1
1741 for k2 in range(cg_maxiter):
1742 if np.add.reduce(np.abs(ri)) <= termcond:
1743 break
1744 if fhess is None:
1745 if fhess_p is None:
1746 Ap = approx_fhess_p(xk, psupi, fprime, epsilon)
1747 else:
1748 Ap = fhess_p(xk, psupi, *args)
1749 hcalls = hcalls + 1
1750 else:
1751 Ap = np.dot(A, psupi)
1752 # check curvature
1753 Ap = asarray(Ap).squeeze() # get rid of matrices...
1754 curv = np.dot(psupi, Ap)
1755 if 0 <= curv <= 3 * float64eps:
1756 break
1757 elif curv < 0:
1758 if (i > 0):
1759 break
1760 else:
1761 # fall back to steepest descent direction
1762 xsupi = dri0 / (-curv) * b
1763 break
1764 alphai = dri0 / curv
1765 xsupi = xsupi + alphai * psupi
1766 ri = ri + alphai * Ap
1767 dri1 = np.dot(ri, ri)
1768 betai = dri1 / dri0
1769 psupi = -ri + betai * psupi
1770 i = i + 1
1771 dri0 = dri1 # update np.dot(ri,ri) for next time.
1772 else:
1773 # curvature keeps increasing, bail out
1774 msg = ("Warning: CG iterations didn't converge. The Hessian is not "
1775 "positive definite.")
1776 return terminate(3, msg)
1778 pk = xsupi # search direction is solution to system.
1779 gfk = -b # gradient at xk
1781 try:
1782 alphak, fc, gc, old_fval, old_old_fval, gfkp1 = \
1783 _line_search_wolfe12(f, fprime, xk, pk, gfk,
1784 old_fval, old_old_fval)
1785 except _LineSearchError:
1786 # Line search failed to find a better solution.
1787 msg = "Warning: " + _status_message['pr_loss']
1788 return terminate(2, msg)
1790 update = alphak * pk
1791 xk = xk + update # upcast if necessary
1792 if callback is not None:
1793 callback(xk)
1794 if retall:
1795 allvecs.append(xk)
1796 k += 1
1797 else:
1798 if np.isnan(old_fval) or np.isnan(update).any():
1799 return terminate(3, _status_message['nan'])
1801 msg = _status_message['success']
1802 return terminate(0, msg)
1805def fminbound(func, x1, x2, args=(), xtol=1e-5, maxfun=500,
1806 full_output=0, disp=1):
1807 """Bounded minimization for scalar functions.
1809 Parameters
1810 ----------
1811 func : callable f(x,*args)
1812 Objective function to be minimized (must accept and return scalars).
1813 x1, x2 : float or array scalar
1814 The optimization bounds.
1815 args : tuple, optional
1816 Extra arguments passed to function.
1817 xtol : float, optional
1818 The convergence tolerance.
1819 maxfun : int, optional
1820 Maximum number of function evaluations allowed.
1821 full_output : bool, optional
1822 If True, return optional outputs.
1823 disp : int, optional
1824 If non-zero, print messages.
1825 0 : no message printing.
1826 1 : non-convergence notification messages only.
1827 2 : print a message on convergence too.
1828 3 : print iteration results.
1831 Returns
1832 -------
1833 xopt : ndarray
1834 Parameters (over given interval) which minimize the
1835 objective function.
1836 fval : number
1837 The function value at the minimum point.
1838 ierr : int
1839 An error flag (0 if converged, 1 if maximum number of
1840 function calls reached).
1841 numfunc : int
1842 The number of function calls made.
1844 See also
1845 --------
1846 minimize_scalar: Interface to minimization algorithms for scalar
1847 univariate functions. See the 'Bounded' `method` in particular.
1849 Notes
1850 -----
1851 Finds a local minimizer of the scalar function `func` in the
1852 interval x1 < xopt < x2 using Brent's method. (See `brent`
1853 for auto-bracketing.)
1855 Examples
1856 --------
1857 `fminbound` finds the minimum of the function in the given range.
1858 The following examples illustrate the same
1860 >>> def f(x):
1861 ... return x**2
1863 >>> from scipy import optimize
1865 >>> minimum = optimize.fminbound(f, -1, 2)
1866 >>> minimum
1867 0.0
1868 >>> minimum = optimize.fminbound(f, 1, 2)
1869 >>> minimum
1870 1.0000059608609866
1871 """
1872 options = {'xatol': xtol,
1873 'maxiter': maxfun,
1874 'disp': disp}
1876 res = _minimize_scalar_bounded(func, (x1, x2), args, **options)
1877 if full_output:
1878 return res['x'], res['fun'], res['status'], res['nfev']
1879 else:
1880 return res['x']
1883def _minimize_scalar_bounded(func, bounds, args=(),
1884 xatol=1e-5, maxiter=500, disp=0,
1885 **unknown_options):
1886 """
1887 Options
1888 -------
1889 maxiter : int
1890 Maximum number of iterations to perform.
1891 disp: int, optional
1892 If non-zero, print messages.
1893 0 : no message printing.
1894 1 : non-convergence notification messages only.
1895 2 : print a message on convergence too.
1896 3 : print iteration results.
1897 xatol : float
1898 Absolute error in solution `xopt` acceptable for convergence.
1900 """
1901 _check_unknown_options(unknown_options)
1902 maxfun = maxiter
1903 # Test bounds are of correct form
1904 if len(bounds) != 2:
1905 raise ValueError('bounds must have two elements.')
1906 x1, x2 = bounds
1908 if not (is_array_scalar(x1) and is_array_scalar(x2)):
1909 raise ValueError("Optimization bounds must be scalars"
1910 " or array scalars.")
1911 if x1 > x2:
1912 raise ValueError("The lower bound exceeds the upper bound.")
1914 flag = 0
1915 header = ' Func-count x f(x) Procedure'
1916 step = ' initial'
1918 sqrt_eps = sqrt(2.2e-16)
1919 golden_mean = 0.5 * (3.0 - sqrt(5.0))
1920 a, b = x1, x2
1921 fulc = a + golden_mean * (b - a)
1922 nfc, xf = fulc, fulc
1923 rat = e = 0.0
1924 x = xf
1925 fx = func(x, *args)
1926 num = 1
1927 fmin_data = (1, xf, fx)
1928 fu = np.inf
1930 ffulc = fnfc = fx
1931 xm = 0.5 * (a + b)
1932 tol1 = sqrt_eps * np.abs(xf) + xatol / 3.0
1933 tol2 = 2.0 * tol1
1935 if disp > 2:
1936 print(" ")
1937 print(header)
1938 print("%5.0f %12.6g %12.6g %s" % (fmin_data + (step,)))
1940 while (np.abs(xf - xm) > (tol2 - 0.5 * (b - a))):
1941 golden = 1
1942 # Check for parabolic fit
1943 if np.abs(e) > tol1:
1944 golden = 0
1945 r = (xf - nfc) * (fx - ffulc)
1946 q = (xf - fulc) * (fx - fnfc)
1947 p = (xf - fulc) * q - (xf - nfc) * r
1948 q = 2.0 * (q - r)
1949 if q > 0.0:
1950 p = -p
1951 q = np.abs(q)
1952 r = e
1953 e = rat
1955 # Check for acceptability of parabola
1956 if ((np.abs(p) < np.abs(0.5*q*r)) and (p > q*(a - xf)) and
1957 (p < q * (b - xf))):
1958 rat = (p + 0.0) / q
1959 x = xf + rat
1960 step = ' parabolic'
1962 if ((x - a) < tol2) or ((b - x) < tol2):
1963 si = np.sign(xm - xf) + ((xm - xf) == 0)
1964 rat = tol1 * si
1965 else: # do a golden-section step
1966 golden = 1
1968 if golden: # do a golden-section step
1969 if xf >= xm:
1970 e = a - xf
1971 else:
1972 e = b - xf
1973 rat = golden_mean*e
1974 step = ' golden'
1976 si = np.sign(rat) + (rat == 0)
1977 x = xf + si * np.max([np.abs(rat), tol1])
1978 fu = func(x, *args)
1979 num += 1
1980 fmin_data = (num, x, fu)
1981 if disp > 2:
1982 print("%5.0f %12.6g %12.6g %s" % (fmin_data + (step,)))
1984 if fu <= fx:
1985 if x >= xf:
1986 a = xf
1987 else:
1988 b = xf
1989 fulc, ffulc = nfc, fnfc
1990 nfc, fnfc = xf, fx
1991 xf, fx = x, fu
1992 else:
1993 if x < xf:
1994 a = x
1995 else:
1996 b = x
1997 if (fu <= fnfc) or (nfc == xf):
1998 fulc, ffulc = nfc, fnfc
1999 nfc, fnfc = x, fu
2000 elif (fu <= ffulc) or (fulc == xf) or (fulc == nfc):
2001 fulc, ffulc = x, fu
2003 xm = 0.5 * (a + b)
2004 tol1 = sqrt_eps * np.abs(xf) + xatol / 3.0
2005 tol2 = 2.0 * tol1
2007 if num >= maxfun:
2008 flag = 1
2009 break
2011 if np.isnan(xf) or np.isnan(fx) or np.isnan(fu):
2012 flag = 2
2014 fval = fx
2015 if disp > 0:
2016 _endprint(x, flag, fval, maxfun, xatol, disp)
2018 result = OptimizeResult(fun=fval, status=flag, success=(flag == 0),
2019 message={0: 'Solution found.',
2020 1: 'Maximum number of function calls '
2021 'reached.',
2022 2: _status_message['nan']}.get(flag, ''),
2023 x=xf, nfev=num)
2025 return result
2028class Brent:
2029 #need to rethink design of __init__
2030 def __init__(self, func, args=(), tol=1.48e-8, maxiter=500,
2031 full_output=0):
2032 self.func = func
2033 self.args = args
2034 self.tol = tol
2035 self.maxiter = maxiter
2036 self._mintol = 1.0e-11
2037 self._cg = 0.3819660
2038 self.xmin = None
2039 self.fval = None
2040 self.iter = 0
2041 self.funcalls = 0
2043 # need to rethink design of set_bracket (new options, etc.)
2044 def set_bracket(self, brack=None):
2045 self.brack = brack
2047 def get_bracket_info(self):
2048 #set up
2049 func = self.func
2050 args = self.args
2051 brack = self.brack
2052 ### BEGIN core bracket_info code ###
2053 ### carefully DOCUMENT any CHANGES in core ##
2054 if brack is None:
2055 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, args=args)
2056 elif len(brack) == 2:
2057 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],
2058 xb=brack[1], args=args)
2059 elif len(brack) == 3:
2060 xa, xb, xc = brack
2061 if (xa > xc): # swap so xa < xc can be assumed
2062 xc, xa = xa, xc
2063 if not ((xa < xb) and (xb < xc)):
2064 raise ValueError("Not a bracketing interval.")
2065 fa = func(*((xa,) + args))
2066 fb = func(*((xb,) + args))
2067 fc = func(*((xc,) + args))
2068 if not ((fb < fa) and (fb < fc)):
2069 raise ValueError("Not a bracketing interval.")
2070 funcalls = 3
2071 else:
2072 raise ValueError("Bracketing interval must be "
2073 "length 2 or 3 sequence.")
2074 ### END core bracket_info code ###
2076 return xa, xb, xc, fa, fb, fc, funcalls
2078 def optimize(self):
2079 # set up for optimization
2080 func = self.func
2081 xa, xb, xc, fa, fb, fc, funcalls = self.get_bracket_info()
2082 _mintol = self._mintol
2083 _cg = self._cg
2084 #################################
2085 #BEGIN CORE ALGORITHM
2086 #################################
2087 x = w = v = xb
2088 fw = fv = fx = func(*((x,) + self.args))
2089 if (xa < xc):
2090 a = xa
2091 b = xc
2092 else:
2093 a = xc
2094 b = xa
2095 deltax = 0.0
2096 funcalls += 1
2097 iter = 0
2098 while (iter < self.maxiter):
2099 tol1 = self.tol * np.abs(x) + _mintol
2100 tol2 = 2.0 * tol1
2101 xmid = 0.5 * (a + b)
2102 # check for convergence
2103 if np.abs(x - xmid) < (tol2 - 0.5 * (b - a)):
2104 break
2105 # XXX In the first iteration, rat is only bound in the true case
2106 # of this conditional. This used to cause an UnboundLocalError
2107 # (gh-4140). It should be set before the if (but to what?).
2108 if (np.abs(deltax) <= tol1):
2109 if (x >= xmid):
2110 deltax = a - x # do a golden section step
2111 else:
2112 deltax = b - x
2113 rat = _cg * deltax
2114 else: # do a parabolic step
2115 tmp1 = (x - w) * (fx - fv)
2116 tmp2 = (x - v) * (fx - fw)
2117 p = (x - v) * tmp2 - (x - w) * tmp1
2118 tmp2 = 2.0 * (tmp2 - tmp1)
2119 if (tmp2 > 0.0):
2120 p = -p
2121 tmp2 = np.abs(tmp2)
2122 dx_temp = deltax
2123 deltax = rat
2124 # check parabolic fit
2125 if ((p > tmp2 * (a - x)) and (p < tmp2 * (b - x)) and
2126 (np.abs(p) < np.abs(0.5 * tmp2 * dx_temp))):
2127 rat = p * 1.0 / tmp2 # if parabolic step is useful.
2128 u = x + rat
2129 if ((u - a) < tol2 or (b - u) < tol2):
2130 if xmid - x >= 0:
2131 rat = tol1
2132 else:
2133 rat = -tol1
2134 else:
2135 if (x >= xmid):
2136 deltax = a - x # if it's not do a golden section step
2137 else:
2138 deltax = b - x
2139 rat = _cg * deltax
2141 if (np.abs(rat) < tol1): # update by at least tol1
2142 if rat >= 0:
2143 u = x + tol1
2144 else:
2145 u = x - tol1
2146 else:
2147 u = x + rat
2148 fu = func(*((u,) + self.args)) # calculate new output value
2149 funcalls += 1
2151 if (fu > fx): # if it's bigger than current
2152 if (u < x):
2153 a = u
2154 else:
2155 b = u
2156 if (fu <= fw) or (w == x):
2157 v = w
2158 w = u
2159 fv = fw
2160 fw = fu
2161 elif (fu <= fv) or (v == x) or (v == w):
2162 v = u
2163 fv = fu
2164 else:
2165 if (u >= x):
2166 a = x
2167 else:
2168 b = x
2169 v = w
2170 w = x
2171 x = u
2172 fv = fw
2173 fw = fx
2174 fx = fu
2176 iter += 1
2177 #################################
2178 #END CORE ALGORITHM
2179 #################################
2181 self.xmin = x
2182 self.fval = fx
2183 self.iter = iter
2184 self.funcalls = funcalls
2186 def get_result(self, full_output=False):
2187 if full_output:
2188 return self.xmin, self.fval, self.iter, self.funcalls
2189 else:
2190 return self.xmin
2193def brent(func, args=(), brack=None, tol=1.48e-8, full_output=0, maxiter=500):
2194 """
2195 Given a function of one variable and a possible bracket, return
2196 the local minimum of the function isolated to a fractional precision
2197 of tol.
2199 Parameters
2200 ----------
2201 func : callable f(x,*args)
2202 Objective function.
2203 args : tuple, optional
2204 Additional arguments (if present).
2205 brack : tuple, optional
2206 Either a triple (xa,xb,xc) where xa<xb<xc and func(xb) <
2207 func(xa), func(xc) or a pair (xa,xb) which are used as a
2208 starting interval for a downhill bracket search (see
2209 `bracket`). Providing the pair (xa,xb) does not always mean
2210 the obtained solution will satisfy xa<=x<=xb.
2211 tol : float, optional
2212 Stop if between iteration change is less than `tol`.
2213 full_output : bool, optional
2214 If True, return all output args (xmin, fval, iter,
2215 funcalls).
2216 maxiter : int, optional
2217 Maximum number of iterations in solution.
2219 Returns
2220 -------
2221 xmin : ndarray
2222 Optimum point.
2223 fval : float
2224 Optimum value.
2225 iter : int
2226 Number of iterations.
2227 funcalls : int
2228 Number of objective function evaluations made.
2230 See also
2231 --------
2232 minimize_scalar: Interface to minimization algorithms for scalar
2233 univariate functions. See the 'Brent' `method` in particular.
2235 Notes
2236 -----
2237 Uses inverse parabolic interpolation when possible to speed up
2238 convergence of golden section method.
2240 Does not ensure that the minimum lies in the range specified by
2241 `brack`. See `fminbound`.
2243 Examples
2244 --------
2245 We illustrate the behaviour of the function when `brack` is of
2246 size 2 and 3 respectively. In the case where `brack` is of the
2247 form (xa,xb), we can see for the given values, the output need
2248 not necessarily lie in the range (xa,xb).
2250 >>> def f(x):
2251 ... return x**2
2253 >>> from scipy import optimize
2255 >>> minimum = optimize.brent(f,brack=(1,2))
2256 >>> minimum
2257 0.0
2258 >>> minimum = optimize.brent(f,brack=(-1,0.5,2))
2259 >>> minimum
2260 -2.7755575615628914e-17
2262 """
2263 options = {'xtol': tol,
2264 'maxiter': maxiter}
2265 res = _minimize_scalar_brent(func, brack, args, **options)
2266 if full_output:
2267 return res['x'], res['fun'], res['nit'], res['nfev']
2268 else:
2269 return res['x']
2272def _minimize_scalar_brent(func, brack=None, args=(),
2273 xtol=1.48e-8, maxiter=500,
2274 **unknown_options):
2275 """
2276 Options
2277 -------
2278 maxiter : int
2279 Maximum number of iterations to perform.
2280 xtol : float
2281 Relative error in solution `xopt` acceptable for convergence.
2283 Notes
2284 -----
2285 Uses inverse parabolic interpolation when possible to speed up
2286 convergence of golden section method.
2288 """
2289 _check_unknown_options(unknown_options)
2290 tol = xtol
2291 if tol < 0:
2292 raise ValueError('tolerance should be >= 0, got %r' % tol)
2294 brent = Brent(func=func, args=args, tol=tol,
2295 full_output=True, maxiter=maxiter)
2296 brent.set_bracket(brack)
2297 brent.optimize()
2298 x, fval, nit, nfev = brent.get_result(full_output=True)
2300 success = nit < maxiter and not (np.isnan(x) or np.isnan(fval))
2302 return OptimizeResult(fun=fval, x=x, nit=nit, nfev=nfev,
2303 success=success)
2306def golden(func, args=(), brack=None, tol=_epsilon,
2307 full_output=0, maxiter=5000):
2308 """
2309 Return the minimum of a function of one variable using golden section
2310 method.
2312 Given a function of one variable and a possible bracketing interval,
2313 return the minimum of the function isolated to a fractional precision of
2314 tol.
2316 Parameters
2317 ----------
2318 func : callable func(x,*args)
2319 Objective function to minimize.
2320 args : tuple, optional
2321 Additional arguments (if present), passed to func.
2322 brack : tuple, optional
2323 Triple (a,b,c), where (a<b<c) and func(b) <
2324 func(a),func(c). If bracket consists of two numbers (a,
2325 c), then they are assumed to be a starting interval for a
2326 downhill bracket search (see `bracket`); it doesn't always
2327 mean that obtained solution will satisfy a<=x<=c.
2328 tol : float, optional
2329 x tolerance stop criterion
2330 full_output : bool, optional
2331 If True, return optional outputs.
2332 maxiter : int
2333 Maximum number of iterations to perform.
2335 See also
2336 --------
2337 minimize_scalar: Interface to minimization algorithms for scalar
2338 univariate functions. See the 'Golden' `method` in particular.
2340 Notes
2341 -----
2342 Uses analog of bisection method to decrease the bracketed
2343 interval.
2345 Examples
2346 --------
2347 We illustrate the behaviour of the function when `brack` is of
2348 size 2 and 3, respectively. In the case where `brack` is of the
2349 form (xa,xb), we can see for the given values, the output need
2350 not necessarily lie in the range ``(xa, xb)``.
2352 >>> def f(x):
2353 ... return x**2
2355 >>> from scipy import optimize
2357 >>> minimum = optimize.golden(f, brack=(1, 2))
2358 >>> minimum
2359 1.5717277788484873e-162
2360 >>> minimum = optimize.golden(f, brack=(-1, 0.5, 2))
2361 >>> minimum
2362 -1.5717277788484873e-162
2364 """
2365 options = {'xtol': tol, 'maxiter': maxiter}
2366 res = _minimize_scalar_golden(func, brack, args, **options)
2367 if full_output:
2368 return res['x'], res['fun'], res['nfev']
2369 else:
2370 return res['x']
2373def _minimize_scalar_golden(func, brack=None, args=(),
2374 xtol=_epsilon, maxiter=5000, **unknown_options):
2375 """
2376 Options
2377 -------
2378 maxiter : int
2379 Maximum number of iterations to perform.
2380 xtol : float
2381 Relative error in solution `xopt` acceptable for convergence.
2383 """
2384 _check_unknown_options(unknown_options)
2385 tol = xtol
2386 if brack is None:
2387 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, args=args)
2388 elif len(brack) == 2:
2389 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],
2390 xb=brack[1], args=args)
2391 elif len(brack) == 3:
2392 xa, xb, xc = brack
2393 if (xa > xc): # swap so xa < xc can be assumed
2394 xc, xa = xa, xc
2395 if not ((xa < xb) and (xb < xc)):
2396 raise ValueError("Not a bracketing interval.")
2397 fa = func(*((xa,) + args))
2398 fb = func(*((xb,) + args))
2399 fc = func(*((xc,) + args))
2400 if not ((fb < fa) and (fb < fc)):
2401 raise ValueError("Not a bracketing interval.")
2402 funcalls = 3
2403 else:
2404 raise ValueError("Bracketing interval must be length 2 or 3 sequence.")
2406 _gR = 0.61803399 # golden ratio conjugate: 2.0/(1.0+sqrt(5.0))
2407 _gC = 1.0 - _gR
2408 x3 = xc
2409 x0 = xa
2410 if (np.abs(xc - xb) > np.abs(xb - xa)):
2411 x1 = xb
2412 x2 = xb + _gC * (xc - xb)
2413 else:
2414 x2 = xb
2415 x1 = xb - _gC * (xb - xa)
2416 f1 = func(*((x1,) + args))
2417 f2 = func(*((x2,) + args))
2418 funcalls += 2
2419 nit = 0
2420 for i in range(maxiter):
2421 if np.abs(x3 - x0) <= tol * (np.abs(x1) + np.abs(x2)):
2422 break
2423 if (f2 < f1):
2424 x0 = x1
2425 x1 = x2
2426 x2 = _gR * x1 + _gC * x3
2427 f1 = f2
2428 f2 = func(*((x2,) + args))
2429 else:
2430 x3 = x2
2431 x2 = x1
2432 x1 = _gR * x2 + _gC * x0
2433 f2 = f1
2434 f1 = func(*((x1,) + args))
2435 funcalls += 1
2436 nit += 1
2437 if (f1 < f2):
2438 xmin = x1
2439 fval = f1
2440 else:
2441 xmin = x2
2442 fval = f2
2444 success = nit < maxiter and not (np.isnan(fval) or np.isnan(xmin))
2446 return OptimizeResult(fun=fval, nfev=funcalls, x=xmin, nit=nit,
2447 success=success)
2450def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
2451 """
2452 Bracket the minimum of the function.
2454 Given a function and distinct initial points, search in the
2455 downhill direction (as defined by the initial points) and return
2456 new points xa, xb, xc that bracket the minimum of the function
2457 f(xa) > f(xb) < f(xc). It doesn't always mean that obtained
2458 solution will satisfy xa<=x<=xb.
2460 Parameters
2461 ----------
2462 func : callable f(x,*args)
2463 Objective function to minimize.
2464 xa, xb : float, optional
2465 Bracketing interval. Defaults `xa` to 0.0, and `xb` to 1.0.
2466 args : tuple, optional
2467 Additional arguments (if present), passed to `func`.
2468 grow_limit : float, optional
2469 Maximum grow limit. Defaults to 110.0
2470 maxiter : int, optional
2471 Maximum number of iterations to perform. Defaults to 1000.
2473 Returns
2474 -------
2475 xa, xb, xc : float
2476 Bracket.
2477 fa, fb, fc : float
2478 Objective function values in bracket.
2479 funcalls : int
2480 Number of function evaluations made.
2482 Examples
2483 --------
2484 This function can find a downward convex region of a function:
2486 >>> import matplotlib.pyplot as plt
2487 >>> from scipy.optimize import bracket
2488 >>> def f(x):
2489 ... return 10*x**2 + 3*x + 5
2490 >>> x = np.linspace(-2, 2)
2491 >>> y = f(x)
2492 >>> init_xa, init_xb = 0, 1
2493 >>> xa, xb, xc, fa, fb, fc, funcalls = bracket(f, xa=init_xa, xb=init_xb)
2494 >>> plt.axvline(x=init_xa, color="k", linestyle="--")
2495 >>> plt.axvline(x=init_xb, color="k", linestyle="--")
2496 >>> plt.plot(x, y, "-k")
2497 >>> plt.plot(xa, fa, "bx")
2498 >>> plt.plot(xb, fb, "rx")
2499 >>> plt.plot(xc, fc, "bx")
2500 >>> plt.show()
2502 """
2503 _gold = 1.618034 # golden ratio: (1.0+sqrt(5.0))/2.0
2504 _verysmall_num = 1e-21
2505 fa = func(*(xa,) + args)
2506 fb = func(*(xb,) + args)
2507 if (fa < fb): # Switch so fa > fb
2508 xa, xb = xb, xa
2509 fa, fb = fb, fa
2510 xc = xb + _gold * (xb - xa)
2511 fc = func(*((xc,) + args))
2512 funcalls = 3
2513 iter = 0
2514 while (fc < fb):
2515 tmp1 = (xb - xa) * (fb - fc)
2516 tmp2 = (xb - xc) * (fb - fa)
2517 val = tmp2 - tmp1
2518 if np.abs(val) < _verysmall_num:
2519 denom = 2.0 * _verysmall_num
2520 else:
2521 denom = 2.0 * val
2522 w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom
2523 wlim = xb + grow_limit * (xc - xb)
2524 if iter > maxiter:
2525 raise RuntimeError("Too many iterations.")
2526 iter += 1
2527 if (w - xc) * (xb - w) > 0.0:
2528 fw = func(*((w,) + args))
2529 funcalls += 1
2530 if (fw < fc):
2531 xa = xb
2532 xb = w
2533 fa = fb
2534 fb = fw
2535 return xa, xb, xc, fa, fb, fc, funcalls
2536 elif (fw > fb):
2537 xc = w
2538 fc = fw
2539 return xa, xb, xc, fa, fb, fc, funcalls
2540 w = xc + _gold * (xc - xb)
2541 fw = func(*((w,) + args))
2542 funcalls += 1
2543 elif (w - wlim)*(wlim - xc) >= 0.0:
2544 w = wlim
2545 fw = func(*((w,) + args))
2546 funcalls += 1
2547 elif (w - wlim)*(xc - w) > 0.0:
2548 fw = func(*((w,) + args))
2549 funcalls += 1
2550 if (fw < fc):
2551 xb = xc
2552 xc = w
2553 w = xc + _gold * (xc - xb)
2554 fb = fc
2555 fc = fw
2556 fw = func(*((w,) + args))
2557 funcalls += 1
2558 else:
2559 w = xc + _gold * (xc - xb)
2560 fw = func(*((w,) + args))
2561 funcalls += 1
2562 xa = xb
2563 xb = xc
2564 xc = w
2565 fa = fb
2566 fb = fc
2567 fc = fw
2568 return xa, xb, xc, fa, fb, fc, funcalls
2571def _line_for_search(x0, alpha, lower_bound, upper_bound):
2572 """
2573 Given a parameter vector ``x0`` with length ``n`` and a direction
2574 vector ``alpha`` with length ``n``, and lower and upper bounds on
2575 each of the ``n`` parameters, what are the bounds on a scalar
2576 ``l`` such that ``lower_bound <= x0 + alpha * l <= upper_bound``.
2579 Parameters
2580 ----------
2581 x0 : np.array.
2582 The vector representing the current location.
2583 Note ``np.shape(x0) == (n,)``.
2584 alpha : np.array.
2585 The vector representing the direction.
2586 Note ``np.shape(alpha) == (n,)``.
2587 lower_bound : np.array.
2588 The lower bounds for each parameter in ``x0``. If the ``i``th
2589 parameter in ``x0`` is unbounded below, then ``lower_bound[i]``
2590 should be ``-np.inf``.
2591 Note ``np.shape(lower_bound) == (n,)``.
2592 upper_bound : np.array.
2593 The upper bounds for each parameter in ``x0``. If the ``i``th
2594 parameter in ``x0`` is unbounded above, then ``upper_bound[i]``
2595 should be ``np.inf``.
2596 Note ``np.shape(upper_bound) == (n,)``.
2598 Returns
2599 -------
2600 res : tuple ``(lmin, lmax)``
2601 The bounds for ``l`` such that
2602 ``lower_bound[i] <= x0[i] + alpha[i] * l <= upper_bound[i]``
2603 for all ``i``.
2605 """
2606 # get nonzero indices of alpha so we don't get any zero division errors.
2607 # alpha will not be all zero, since it is called from _linesearch_powell
2608 # where we have a check for this.
2609 nonzero, = alpha.nonzero()
2610 lower_bound, upper_bound = lower_bound[nonzero], upper_bound[nonzero]
2611 x0, alpha = x0[nonzero], alpha[nonzero]
2612 low = (lower_bound - x0) / alpha
2613 high = (upper_bound - x0) / alpha
2615 # positive and negative indices
2616 pos = alpha > 0
2618 lmin_pos = np.where(pos, low, 0)
2619 lmin_neg = np.where(pos, 0, high)
2620 lmax_pos = np.where(pos, high, 0)
2621 lmax_neg = np.where(pos, 0, low)
2623 lmin = np.max(lmin_pos + lmin_neg)
2624 lmax = np.min(lmax_pos + lmax_neg)
2626 # if x0 is outside the bounds, then it is possible that there is
2627 # no way to get back in the bounds for the parameters being updated
2628 # with the current direction alpha.
2629 # when this happens, lmax < lmin.
2630 # If this is the case, then we can just return (0, 0)
2631 return (lmin, lmax) if lmax >= lmin else (0, 0)
2634def _linesearch_powell(func, p, xi, tol=1e-3,
2635 lower_bound=None, upper_bound=None, fval=None):
2636 """Line-search algorithm using fminbound.
2638 Find the minimium of the function ``func(x0 + alpha*direc)``.
2640 lower_bound : np.array.
2641 The lower bounds for each parameter in ``x0``. If the ``i``th
2642 parameter in ``x0`` is unbounded below, then ``lower_bound[i]``
2643 should be ``-np.inf``.
2644 Note ``np.shape(lower_bound) == (n,)``.
2645 upper_bound : np.array.
2646 The upper bounds for each parameter in ``x0``. If the ``i``th
2647 parameter in ``x0`` is unbounded above, then ``upper_bound[i]``
2648 should be ``np.inf``.
2649 Note ``np.shape(upper_bound) == (n,)``.
2650 fval : number.
2651 ``fval`` is equal to ``func(p)``, the idea is just to avoid
2652 recomputing it so we can limit the ``fevals``.
2654 """
2655 def myfunc(alpha):
2656 return func(p + alpha*xi)
2658 # if xi is zero, then don't optimize
2659 if not np.any(xi):
2660 return ((fval, p, xi) if fval is not None else (func(p), p, xi))
2661 elif lower_bound is None and upper_bound is None:
2662 # non-bounded minimization
2663 alpha_min, fret, _, _ = brent(myfunc, full_output=1, tol=tol)
2664 xi = alpha_min * xi
2665 return squeeze(fret), p + xi, xi
2666 else:
2667 bound = _line_for_search(p, xi, lower_bound, upper_bound)
2668 if np.isneginf(bound[0]) and np.isposinf(bound[1]):
2669 # equivalent to unbounded
2670 return _linesearch_powell(func, p, xi, fval=fval, tol=tol)
2671 elif not np.isneginf(bound[0]) and not np.isposinf(bound[1]):
2672 # we can use a bounded scalar minimization
2673 res = _minimize_scalar_bounded(myfunc, bound, xatol=tol / 100)
2674 xi = res.x * xi
2675 return squeeze(res.fun), p + xi, xi
2676 else:
2677 # only bounded on one side. use the tangent function to convert
2678 # the infinity bound to a finite bound. The new bounded region
2679 # is a subregion of the region bounded by -np.pi/2 and np.pi/2.
2680 bound = np.arctan(bound[0]), np.arctan(bound[1])
2681 res = _minimize_scalar_bounded(
2682 lambda x: myfunc(np.tan(x)),
2683 bound,
2684 xatol=tol / 100)
2685 xi = np.tan(res.x) * xi
2686 return squeeze(res.fun), p + xi, xi
2689def fmin_powell(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None,
2690 maxfun=None, full_output=0, disp=1, retall=0, callback=None,
2691 direc=None):
2692 """
2693 Minimize a function using modified Powell's method.
2695 This method only uses function values, not derivatives.
2697 Parameters
2698 ----------
2699 func : callable f(x,*args)
2700 Objective function to be minimized.
2701 x0 : ndarray
2702 Initial guess.
2703 args : tuple, optional
2704 Extra arguments passed to func.
2705 xtol : float, optional
2706 Line-search error tolerance.
2707 ftol : float, optional
2708 Relative error in ``func(xopt)`` acceptable for convergence.
2709 maxiter : int, optional
2710 Maximum number of iterations to perform.
2711 maxfun : int, optional
2712 Maximum number of function evaluations to make.
2713 full_output : bool, optional
2714 If True, ``fopt``, ``xi``, ``direc``, ``iter``, ``funcalls``, and
2715 ``warnflag`` are returned.
2716 disp : bool, optional
2717 If True, print convergence messages.
2718 retall : bool, optional
2719 If True, return a list of the solution at each iteration.
2720 callback : callable, optional
2721 An optional user-supplied function, called after each
2722 iteration. Called as ``callback(xk)``, where ``xk`` is the
2723 current parameter vector.
2724 direc : ndarray, optional
2725 Initial fitting step and parameter order set as an (N, N) array, where N
2726 is the number of fitting parameters in `x0`. Defaults to step size 1.0
2727 fitting all parameters simultaneously (``np.ones((N, N))``). To
2728 prevent initial consideration of values in a step or to change initial
2729 step size, set to 0 or desired step size in the Jth position in the Mth
2730 block, where J is the position in `x0` and M is the desired evaluation
2731 step, with steps being evaluated in index order. Step size and ordering
2732 will change freely as minimization proceeds.
2734 Returns
2735 -------
2736 xopt : ndarray
2737 Parameter which minimizes `func`.
2738 fopt : number
2739 Value of function at minimum: ``fopt = func(xopt)``.
2740 direc : ndarray
2741 Current direction set.
2742 iter : int
2743 Number of iterations.
2744 funcalls : int
2745 Number of function calls made.
2746 warnflag : int
2747 Integer warning flag:
2748 1 : Maximum number of function evaluations.
2749 2 : Maximum number of iterations.
2750 3 : NaN result encountered.
2751 4 : The result is out of the provided bounds.
2752 allvecs : list
2753 List of solutions at each iteration.
2755 See also
2756 --------
2757 minimize: Interface to unconstrained minimization algorithms for
2758 multivariate functions. See the 'Powell' method in particular.
2760 Notes
2761 -----
2762 Uses a modification of Powell's method to find the minimum of
2763 a function of N variables. Powell's method is a conjugate
2764 direction method.
2766 The algorithm has two loops. The outer loop merely iterates over the inner
2767 loop. The inner loop minimizes over each current direction in the direction
2768 set. At the end of the inner loop, if certain conditions are met, the
2769 direction that gave the largest decrease is dropped and replaced with the
2770 difference between the current estimated x and the estimated x from the
2771 beginning of the inner-loop.
2773 The technical conditions for replacing the direction of greatest
2774 increase amount to checking that
2776 1. No further gain can be made along the direction of greatest increase
2777 from that iteration.
2778 2. The direction of greatest increase accounted for a large sufficient
2779 fraction of the decrease in the function value from that iteration of
2780 the inner loop.
2782 References
2783 ----------
2784 Powell M.J.D. (1964) An efficient method for finding the minimum of a
2785 function of several variables without calculating derivatives,
2786 Computer Journal, 7 (2):155-162.
2788 Press W., Teukolsky S.A., Vetterling W.T., and Flannery B.P.:
2789 Numerical Recipes (any edition), Cambridge University Press
2791 Examples
2792 --------
2793 >>> def f(x):
2794 ... return x**2
2796 >>> from scipy import optimize
2798 >>> minimum = optimize.fmin_powell(f, -1)
2799 Optimization terminated successfully.
2800 Current function value: 0.000000
2801 Iterations: 2
2802 Function evaluations: 18
2803 >>> minimum
2804 array(0.0)
2806 """
2807 opts = {'xtol': xtol,
2808 'ftol': ftol,
2809 'maxiter': maxiter,
2810 'maxfev': maxfun,
2811 'disp': disp,
2812 'direc': direc,
2813 'return_all': retall}
2815 res = _minimize_powell(func, x0, args, callback=callback, **opts)
2817 if full_output:
2818 retlist = (res['x'], res['fun'], res['direc'], res['nit'],
2819 res['nfev'], res['status'])
2820 if retall:
2821 retlist += (res['allvecs'], )
2822 return retlist
2823 else:
2824 if retall:
2825 return res['x'], res['allvecs']
2826 else:
2827 return res['x']
2830def _minimize_powell(func, x0, args=(), callback=None, bounds=None,
2831 xtol=1e-4, ftol=1e-4, maxiter=None, maxfev=None,
2832 disp=False, direc=None, return_all=False,
2833 **unknown_options):
2834 """
2835 Minimization of scalar function of one or more variables using the
2836 modified Powell algorithm.
2838 Options
2839 -------
2840 disp : bool
2841 Set to True to print convergence messages.
2842 xtol : float
2843 Relative error in solution `xopt` acceptable for convergence.
2844 ftol : float
2845 Relative error in ``fun(xopt)`` acceptable for convergence.
2846 maxiter, maxfev : int
2847 Maximum allowed number of iterations and function evaluations.
2848 Will default to ``N*1000``, where ``N`` is the number of
2849 variables, if neither `maxiter` or `maxfev` is set. If both
2850 `maxiter` and `maxfev` are set, minimization will stop at the
2851 first reached.
2852 direc : ndarray
2853 Initial set of direction vectors for the Powell method.
2854 return_all : bool, optional
2855 Set to True to return a list of the best solution at each of the
2856 iterations.
2857 bounds : `Bounds`
2858 If bounds are not provided, then an unbounded line search will be used.
2859 If bounds are provided and the initial guess is within the bounds, then
2860 every function evaluation throughout the minimization procedure will be
2861 within the bounds. If bounds are provided, the initial guess is outside
2862 the bounds, and `direc` is full rank (or left to default), then some
2863 function evaluations during the first iteration may be outside the
2864 bounds, but every function evaluation after the first iteration will be
2865 within the bounds. If `direc` is not full rank, then some parameters may
2866 not be optimized and the solution is not guaranteed to be within the
2867 bounds.
2868 return_all : bool, optional
2869 Set to True to return a list of the best solution at each of the
2870 iterations.
2871 """
2872 _check_unknown_options(unknown_options)
2873 maxfun = maxfev
2874 retall = return_all
2875 # we need to use a mutable object here that we can update in the
2876 # wrapper function
2877 fcalls, func = wrap_function(func, args)
2878 x = asarray(x0).flatten()
2879 if retall:
2880 allvecs = [x]
2881 N = len(x)
2882 # If neither are set, then set both to default
2883 if maxiter is None and maxfun is None:
2884 maxiter = N * 1000
2885 maxfun = N * 1000
2886 elif maxiter is None:
2887 # Convert remaining Nones, to np.inf, unless the other is np.inf, in
2888 # which case use the default to avoid unbounded iteration
2889 if maxfun == np.inf:
2890 maxiter = N * 1000
2891 else:
2892 maxiter = np.inf
2893 elif maxfun is None:
2894 if maxiter == np.inf:
2895 maxfun = N * 1000
2896 else:
2897 maxfun = np.inf
2899 if direc is None:
2900 direc = eye(N, dtype=float)
2901 else:
2902 direc = asarray(direc, dtype=float)
2903 if np.linalg.matrix_rank(direc) != direc.shape[0]:
2904 warnings.warn("direc input is not full rank, some parameters may "
2905 "not be optimized",
2906 OptimizeWarning, 3)
2908 if bounds is None:
2909 # don't make these arrays of all +/- inf. because
2910 # _linesearch_powell will do an unnecessary check of all the elements.
2911 # just keep them None, _linesearch_powell will not have to check
2912 # all the elements.
2913 lower_bound, upper_bound = None, None
2914 else:
2915 # bounds is standardized in _minimize.py.
2916 lower_bound, upper_bound = bounds.lb, bounds.ub
2917 if np.any(lower_bound > x0) or np.any(x0 > upper_bound):
2918 warnings.warn("Initial guess is not within the specified bounds",
2919 OptimizeWarning, 3)
2921 fval = squeeze(func(x))
2922 x1 = x.copy()
2923 iter = 0
2924 ilist = list(range(N))
2925 while True:
2926 fx = fval
2927 bigind = 0
2928 delta = 0.0
2929 for i in ilist:
2930 direc1 = direc[i]
2931 fx2 = fval
2932 fval, x, direc1 = _linesearch_powell(func, x, direc1,
2933 tol=xtol * 100,
2934 lower_bound=lower_bound,
2935 upper_bound=upper_bound,
2936 fval=fval)
2937 if (fx2 - fval) > delta:
2938 delta = fx2 - fval
2939 bigind = i
2940 iter += 1
2941 if callback is not None:
2942 callback(x)
2943 if retall:
2944 allvecs.append(x)
2945 bnd = ftol * (np.abs(fx) + np.abs(fval)) + 1e-20
2946 if 2.0 * (fx - fval) <= bnd:
2947 break
2948 if fcalls[0] >= maxfun:
2949 break
2950 if iter >= maxiter:
2951 break
2952 if np.isnan(fx) and np.isnan(fval):
2953 # Ended up in a nan-region: bail out
2954 break
2956 # Construct the extrapolated point
2957 direc1 = x - x1
2958 x2 = 2*x - x1
2959 x1 = x.copy()
2960 fx2 = squeeze(func(x2))
2962 if (fx > fx2):
2963 t = 2.0*(fx + fx2 - 2.0*fval)
2964 temp = (fx - fval - delta)
2965 t *= temp*temp
2966 temp = fx - fx2
2967 t -= delta*temp*temp
2968 if t < 0.0:
2969 fval, x, direc1 = _linesearch_powell(func, x, direc1,
2970 tol=xtol * 100,
2971 lower_bound=lower_bound,
2972 upper_bound=upper_bound,
2973 fval=fval)
2974 if np.any(direc1):
2975 direc[bigind] = direc[-1]
2976 direc[-1] = direc1
2978 warnflag = 0
2979 # out of bounds is more urgent than exceeding function evals or iters,
2980 # but I don't want to cause inconsistencies by changing the
2981 # established warning flags for maxfev and maxiter, so the out of bounds
2982 # warning flag becomes 3, but is checked for first.
2983 if bounds and (np.any(lower_bound > x) or np.any(x > upper_bound)):
2984 warnflag = 4
2985 msg = _status_message['out_of_bounds']
2986 elif fcalls[0] >= maxfun:
2987 warnflag = 1
2988 msg = _status_message['maxfev']
2989 if disp:
2990 print("Warning: " + msg)
2991 elif iter >= maxiter:
2992 warnflag = 2
2993 msg = _status_message['maxiter']
2994 if disp:
2995 print("Warning: " + msg)
2996 elif np.isnan(fval) or np.isnan(x).any():
2997 warnflag = 3
2998 msg = _status_message['nan']
2999 if disp:
3000 print("Warning: " + msg)
3001 else:
3002 msg = _status_message['success']
3003 if disp:
3004 print(msg)
3005 print(" Current function value: %f" % fval)
3006 print(" Iterations: %d" % iter)
3007 print(" Function evaluations: %d" % fcalls[0])
3009 result = OptimizeResult(fun=fval, direc=direc, nit=iter, nfev=fcalls[0],
3010 status=warnflag, success=(warnflag == 0),
3011 message=msg, x=x)
3012 if retall:
3013 result['allvecs'] = allvecs
3014 return result
3017def _endprint(x, flag, fval, maxfun, xtol, disp):
3018 if flag == 0:
3019 if disp > 1:
3020 print("\nOptimization terminated successfully;\n"
3021 "The returned value satisfies the termination criteria\n"
3022 "(using xtol = ", xtol, ")")
3023 if flag == 1:
3024 if disp:
3025 print("\nMaximum number of function evaluations exceeded --- "
3026 "increase maxfun argument.\n")
3027 if flag == 2:
3028 if disp:
3029 print("\n{}".format(_status_message['nan']))
3030 return
3033def brute(func, ranges, args=(), Ns=20, full_output=0, finish=fmin,
3034 disp=False, workers=1):
3035 """Minimize a function over a given range by brute force.
3037 Uses the "brute force" method, i.e., computes the function's value
3038 at each point of a multidimensional grid of points, to find the global
3039 minimum of the function.
3041 The function is evaluated everywhere in the range with the datatype of the
3042 first call to the function, as enforced by the ``vectorize`` NumPy
3043 function. The value and type of the function evaluation returned when
3044 ``full_output=True`` are affected in addition by the ``finish`` argument
3045 (see Notes).
3047 The brute force approach is inefficient because the number of grid points
3048 increases exponentially - the number of grid points to evaluate is
3049 ``Ns ** len(x)``. Consequently, even with coarse grid spacing, even
3050 moderately sized problems can take a long time to run, and/or run into
3051 memory limitations.
3053 Parameters
3054 ----------
3055 func : callable
3056 The objective function to be minimized. Must be in the
3057 form ``f(x, *args)``, where ``x`` is the argument in
3058 the form of a 1-D array and ``args`` is a tuple of any
3059 additional fixed parameters needed to completely specify
3060 the function.
3061 ranges : tuple
3062 Each component of the `ranges` tuple must be either a
3063 "slice object" or a range tuple of the form ``(low, high)``.
3064 The program uses these to create the grid of points on which
3065 the objective function will be computed. See `Note 2` for
3066 more detail.
3067 args : tuple, optional
3068 Any additional fixed parameters needed to completely specify
3069 the function.
3070 Ns : int, optional
3071 Number of grid points along the axes, if not otherwise
3072 specified. See `Note2`.
3073 full_output : bool, optional
3074 If True, return the evaluation grid and the objective function's
3075 values on it.
3076 finish : callable, optional
3077 An optimization function that is called with the result of brute force
3078 minimization as initial guess. `finish` should take `func` and
3079 the initial guess as positional arguments, and take `args` as
3080 keyword arguments. It may additionally take `full_output`
3081 and/or `disp` as keyword arguments. Use None if no "polishing"
3082 function is to be used. See Notes for more details.
3083 disp : bool, optional
3084 Set to True to print convergence messages from the `finish` callable.
3085 workers : int or map-like callable, optional
3086 If `workers` is an int the grid is subdivided into `workers`
3087 sections and evaluated in parallel (uses
3088 `multiprocessing.Pool <multiprocessing>`).
3089 Supply `-1` to use all cores available to the Process.
3090 Alternatively supply a map-like callable, such as
3091 `multiprocessing.Pool.map` for evaluating the grid in parallel.
3092 This evaluation is carried out as ``workers(func, iterable)``.
3093 Requires that `func` be pickleable.
3095 .. versionadded:: 1.3.0
3097 Returns
3098 -------
3099 x0 : ndarray
3100 A 1-D array containing the coordinates of a point at which the
3101 objective function had its minimum value. (See `Note 1` for
3102 which point is returned.)
3103 fval : float
3104 Function value at the point `x0`. (Returned when `full_output` is
3105 True.)
3106 grid : tuple
3107 Representation of the evaluation grid. It has the same
3108 length as `x0`. (Returned when `full_output` is True.)
3109 Jout : ndarray
3110 Function values at each point of the evaluation
3111 grid, i.e., ``Jout = func(*grid)``. (Returned
3112 when `full_output` is True.)
3114 See Also
3115 --------
3116 basinhopping, differential_evolution
3118 Notes
3119 -----
3120 *Note 1*: The program finds the gridpoint at which the lowest value
3121 of the objective function occurs. If `finish` is None, that is the
3122 point returned. When the global minimum occurs within (or not very far
3123 outside) the grid's boundaries, and the grid is fine enough, that
3124 point will be in the neighborhood of the global minimum.
3126 However, users often employ some other optimization program to
3127 "polish" the gridpoint values, i.e., to seek a more precise
3128 (local) minimum near `brute's` best gridpoint.
3129 The `brute` function's `finish` option provides a convenient way to do
3130 that. Any polishing program used must take `brute's` output as its
3131 initial guess as a positional argument, and take `brute's` input values
3132 for `args` as keyword arguments, otherwise an error will be raised.
3133 It may additionally take `full_output` and/or `disp` as keyword arguments.
3135 `brute` assumes that the `finish` function returns either an
3136 `OptimizeResult` object or a tuple in the form:
3137 ``(xmin, Jmin, ... , statuscode)``, where ``xmin`` is the minimizing
3138 value of the argument, ``Jmin`` is the minimum value of the objective
3139 function, "..." may be some other returned values (which are not used
3140 by `brute`), and ``statuscode`` is the status code of the `finish` program.
3142 Note that when `finish` is not None, the values returned are those
3143 of the `finish` program, *not* the gridpoint ones. Consequently,
3144 while `brute` confines its search to the input grid points,
3145 the `finish` program's results usually will not coincide with any
3146 gridpoint, and may fall outside the grid's boundary. Thus, if a
3147 minimum only needs to be found over the provided grid points, make
3148 sure to pass in `finish=None`.
3150 *Note 2*: The grid of points is a `numpy.mgrid` object.
3151 For `brute` the `ranges` and `Ns` inputs have the following effect.
3152 Each component of the `ranges` tuple can be either a slice object or a
3153 two-tuple giving a range of values, such as (0, 5). If the component is a
3154 slice object, `brute` uses it directly. If the component is a two-tuple
3155 range, `brute` internally converts it to a slice object that interpolates
3156 `Ns` points from its low-value to its high-value, inclusive.
3158 Examples
3159 --------
3160 We illustrate the use of `brute` to seek the global minimum of a function
3161 of two variables that is given as the sum of a positive-definite
3162 quadratic and two deep "Gaussian-shaped" craters. Specifically, define
3163 the objective function `f` as the sum of three other functions,
3164 ``f = f1 + f2 + f3``. We suppose each of these has a signature
3165 ``(z, *params)``, where ``z = (x, y)``, and ``params`` and the functions
3166 are as defined below.
3168 >>> params = (2, 3, 7, 8, 9, 10, 44, -1, 2, 26, 1, -2, 0.5)
3169 >>> def f1(z, *params):
3170 ... x, y = z
3171 ... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
3172 ... return (a * x**2 + b * x * y + c * y**2 + d*x + e*y + f)
3174 >>> def f2(z, *params):
3175 ... x, y = z
3176 ... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
3177 ... return (-g*np.exp(-((x-h)**2 + (y-i)**2) / scale))
3179 >>> def f3(z, *params):
3180 ... x, y = z
3181 ... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
3182 ... return (-j*np.exp(-((x-k)**2 + (y-l)**2) / scale))
3184 >>> def f(z, *params):
3185 ... return f1(z, *params) + f2(z, *params) + f3(z, *params)
3187 Thus, the objective function may have local minima near the minimum
3188 of each of the three functions of which it is composed. To
3189 use `fmin` to polish its gridpoint result, we may then continue as
3190 follows:
3192 >>> rranges = (slice(-4, 4, 0.25), slice(-4, 4, 0.25))
3193 >>> from scipy import optimize
3194 >>> resbrute = optimize.brute(f, rranges, args=params, full_output=True,
3195 ... finish=optimize.fmin)
3196 >>> resbrute[0] # global minimum
3197 array([-1.05665192, 1.80834843])
3198 >>> resbrute[1] # function value at global minimum
3199 -3.4085818767
3201 Note that if `finish` had been set to None, we would have gotten the
3202 gridpoint [-1.0 1.75] where the rounded function value is -2.892.
3204 """
3205 N = len(ranges)
3206 if N > 40:
3207 raise ValueError("Brute Force not possible with more "
3208 "than 40 variables.")
3209 lrange = list(ranges)
3210 for k in range(N):
3211 if type(lrange[k]) is not type(slice(None)):
3212 if len(lrange[k]) < 3:
3213 lrange[k] = tuple(lrange[k]) + (complex(Ns),)
3214 lrange[k] = slice(*lrange[k])
3215 if (N == 1):
3216 lrange = lrange[0]
3218 grid = np.mgrid[lrange]
3220 # obtain an array of parameters that is iterable by a map-like callable
3221 inpt_shape = grid.shape
3222 if (N > 1):
3223 grid = np.reshape(grid, (inpt_shape[0], np.prod(inpt_shape[1:]))).T
3225 wrapped_func = _Brute_Wrapper(func, args)
3227 # iterate over input arrays, possibly in parallel
3228 with MapWrapper(pool=workers) as mapper:
3229 Jout = np.array(list(mapper(wrapped_func, grid)))
3230 if (N == 1):
3231 grid = (grid,)
3232 Jout = np.squeeze(Jout)
3233 elif (N > 1):
3234 Jout = np.reshape(Jout, inpt_shape[1:])
3235 grid = np.reshape(grid.T, inpt_shape)
3237 Nshape = shape(Jout)
3239 indx = argmin(Jout.ravel(), axis=-1)
3240 Nindx = zeros(N, int)
3241 xmin = zeros(N, float)
3242 for k in range(N - 1, -1, -1):
3243 thisN = Nshape[k]
3244 Nindx[k] = indx % Nshape[k]
3245 indx = indx // thisN
3246 for k in range(N):
3247 xmin[k] = grid[k][tuple(Nindx)]
3249 Jmin = Jout[tuple(Nindx)]
3250 if (N == 1):
3251 grid = grid[0]
3252 xmin = xmin[0]
3254 if callable(finish):
3255 # set up kwargs for `finish` function
3256 finish_args = _getfullargspec(finish).args
3257 finish_kwargs = dict()
3258 if 'full_output' in finish_args:
3259 finish_kwargs['full_output'] = 1
3260 if 'disp' in finish_args:
3261 finish_kwargs['disp'] = disp
3262 elif 'options' in finish_args:
3263 # pass 'disp' as `options`
3264 # (e.g., if `finish` is `minimize`)
3265 finish_kwargs['options'] = {'disp': disp}
3267 # run minimizer
3268 res = finish(func, xmin, args=args, **finish_kwargs)
3270 if isinstance(res, OptimizeResult):
3271 xmin = res.x
3272 Jmin = res.fun
3273 success = res.success
3274 else:
3275 xmin = res[0]
3276 Jmin = res[1]
3277 success = res[-1] == 0
3278 if not success:
3279 if disp:
3280 print("Warning: Either final optimization did not succeed "
3281 "or `finish` does not return `statuscode` as its last "
3282 "argument.")
3284 if full_output:
3285 return xmin, Jmin, grid, Jout
3286 else:
3287 return xmin
3290class _Brute_Wrapper(object):
3291 """
3292 Object to wrap user cost function for optimize.brute, allowing picklability
3293 """
3295 def __init__(self, f, args):
3296 self.f = f
3297 self.args = [] if args is None else args
3299 def __call__(self, x):
3300 # flatten needed for one dimensional case.
3301 return self.f(np.asarray(x).flatten(), *self.args)
3304def show_options(solver=None, method=None, disp=True):
3305 """
3306 Show documentation for additional options of optimization solvers.
3308 These are method-specific options that can be supplied through the
3309 ``options`` dict.
3311 Parameters
3312 ----------
3313 solver : str
3314 Type of optimization solver. One of 'minimize', 'minimize_scalar',
3315 'root', or 'linprog'.
3316 method : str, optional
3317 If not given, shows all methods of the specified solver. Otherwise,
3318 show only the options for the specified method. Valid values
3319 corresponds to methods' names of respective solver (e.g., 'BFGS' for
3320 'minimize').
3321 disp : bool, optional
3322 Whether to print the result rather than returning it.
3324 Returns
3325 -------
3326 text
3327 Either None (for disp=True) or the text string (disp=False)
3329 Notes
3330 -----
3331 The solver-specific methods are:
3333 `scipy.optimize.minimize`
3335 - :ref:`Nelder-Mead <optimize.minimize-neldermead>`
3336 - :ref:`Powell <optimize.minimize-powell>`
3337 - :ref:`CG <optimize.minimize-cg>`
3338 - :ref:`BFGS <optimize.minimize-bfgs>`
3339 - :ref:`Newton-CG <optimize.minimize-newtoncg>`
3340 - :ref:`L-BFGS-B <optimize.minimize-lbfgsb>`
3341 - :ref:`TNC <optimize.minimize-tnc>`
3342 - :ref:`COBYLA <optimize.minimize-cobyla>`
3343 - :ref:`SLSQP <optimize.minimize-slsqp>`
3344 - :ref:`dogleg <optimize.minimize-dogleg>`
3345 - :ref:`trust-ncg <optimize.minimize-trustncg>`
3347 `scipy.optimize.root`
3349 - :ref:`hybr <optimize.root-hybr>`
3350 - :ref:`lm <optimize.root-lm>`
3351 - :ref:`broyden1 <optimize.root-broyden1>`
3352 - :ref:`broyden2 <optimize.root-broyden2>`
3353 - :ref:`anderson <optimize.root-anderson>`
3354 - :ref:`linearmixing <optimize.root-linearmixing>`
3355 - :ref:`diagbroyden <optimize.root-diagbroyden>`
3356 - :ref:`excitingmixing <optimize.root-excitingmixing>`
3357 - :ref:`krylov <optimize.root-krylov>`
3358 - :ref:`df-sane <optimize.root-dfsane>`
3360 `scipy.optimize.minimize_scalar`
3362 - :ref:`brent <optimize.minimize_scalar-brent>`
3363 - :ref:`golden <optimize.minimize_scalar-golden>`
3364 - :ref:`bounded <optimize.minimize_scalar-bounded>`
3366 `scipy.optimize.linprog`
3368 - :ref:`simplex <optimize.linprog-simplex>`
3369 - :ref:`interior-point <optimize.linprog-interior-point>`
3371 Examples
3372 --------
3373 We can print documentations of a solver in stdout:
3375 >>> from scipy.optimize import show_options
3376 >>> show_options(solver="minimize")
3377 ...
3379 Specifying a method is possible:
3381 >>> show_options(solver="minimize", method="Nelder-Mead")
3382 ...
3384 We can also get the documentations as a string:
3386 >>> show_options(solver="minimize", method="Nelder-Mead", disp=False)
3387 Minimization of scalar function of one or more variables using the ...
3389 """
3390 import textwrap
3392 doc_routines = {
3393 'minimize': (
3394 ('bfgs', 'scipy.optimize.optimize._minimize_bfgs'),
3395 ('cg', 'scipy.optimize.optimize._minimize_cg'),
3396 ('cobyla', 'scipy.optimize.cobyla._minimize_cobyla'),
3397 ('dogleg', 'scipy.optimize._trustregion_dogleg._minimize_dogleg'),
3398 ('l-bfgs-b', 'scipy.optimize.lbfgsb._minimize_lbfgsb'),
3399 ('nelder-mead', 'scipy.optimize.optimize._minimize_neldermead'),
3400 ('newton-cg', 'scipy.optimize.optimize._minimize_newtoncg'),
3401 ('powell', 'scipy.optimize.optimize._minimize_powell'),
3402 ('slsqp', 'scipy.optimize.slsqp._minimize_slsqp'),
3403 ('tnc', 'scipy.optimize.tnc._minimize_tnc'),
3404 ('trust-ncg', 'scipy.optimize._trustregion_ncg._minimize_trust_ncg'),
3405 ),
3406 'root': (
3407 ('hybr', 'scipy.optimize.minpack._root_hybr'),
3408 ('lm', 'scipy.optimize._root._root_leastsq'),
3409 ('broyden1', 'scipy.optimize._root._root_broyden1_doc'),
3410 ('broyden2', 'scipy.optimize._root._root_broyden2_doc'),
3411 ('anderson', 'scipy.optimize._root._root_anderson_doc'),
3412 ('diagbroyden', 'scipy.optimize._root._root_diagbroyden_doc'),
3413 ('excitingmixing', 'scipy.optimize._root._root_excitingmixing_doc'),
3414 ('linearmixing', 'scipy.optimize._root._root_linearmixing_doc'),
3415 ('krylov', 'scipy.optimize._root._root_krylov_doc'),
3416 ('df-sane', 'scipy.optimize._spectral._root_df_sane'),
3417 ),
3418 'root_scalar': (
3419 ('bisect', 'scipy.optimize._root_scalar._root_scalar_bisect_doc'),
3420 ('brentq', 'scipy.optimize._root_scalar._root_scalar_brentq_doc'),
3421 ('brenth', 'scipy.optimize._root_scalar._root_scalar_brenth_doc'),
3422 ('ridder', 'scipy.optimize._root_scalar._root_scalar_ridder_doc'),
3423 ('toms748', 'scipy.optimize._root_scalar._root_scalar_toms748_doc'),
3424 ('secant', 'scipy.optimize._root_scalar._root_scalar_secant_doc'),
3425 ('newton', 'scipy.optimize._root_scalar._root_scalar_newton_doc'),
3426 ('halley', 'scipy.optimize._root_scalar._root_scalar_halley_doc'),
3427 ),
3428 'linprog': (
3429 ('simplex', 'scipy.optimize._linprog._linprog_simplex'),
3430 ('interior-point', 'scipy.optimize._linprog._linprog_ip'),
3431 ),
3432 'minimize_scalar': (
3433 ('brent', 'scipy.optimize.optimize._minimize_scalar_brent'),
3434 ('bounded', 'scipy.optimize.optimize._minimize_scalar_bounded'),
3435 ('golden', 'scipy.optimize.optimize._minimize_scalar_golden'),
3436 ),
3437 }
3439 if solver is None:
3440 text = ["\n\n\n========\n", "minimize\n", "========\n"]
3441 text.append(show_options('minimize', disp=False))
3442 text.extend(["\n\n===============\n", "minimize_scalar\n",
3443 "===============\n"])
3444 text.append(show_options('minimize_scalar', disp=False))
3445 text.extend(["\n\n\n====\n", "root\n",
3446 "====\n"])
3447 text.append(show_options('root', disp=False))
3448 text.extend(['\n\n\n=======\n', 'linprog\n',
3449 '=======\n'])
3450 text.append(show_options('linprog', disp=False))
3451 text = "".join(text)
3452 else:
3453 solver = solver.lower()
3454 if solver not in doc_routines:
3455 raise ValueError('Unknown solver %r' % (solver,))
3457 if method is None:
3458 text = []
3459 for name, _ in doc_routines[solver]:
3460 text.extend(["\n\n" + name, "\n" + "="*len(name) + "\n\n"])
3461 text.append(show_options(solver, name, disp=False))
3462 text = "".join(text)
3463 else:
3464 method = method.lower()
3465 methods = dict(doc_routines[solver])
3466 if method not in methods:
3467 raise ValueError("Unknown method %r" % (method,))
3468 name = methods[method]
3470 # Import function object
3471 parts = name.split('.')
3472 mod_name = ".".join(parts[:-1])
3473 __import__(mod_name)
3474 obj = getattr(sys.modules[mod_name], parts[-1])
3476 # Get doc
3477 doc = obj.__doc__
3478 if doc is not None:
3479 text = textwrap.dedent(doc).strip()
3480 else:
3481 text = ""
3483 if disp:
3484 print(text)
3485 return
3486 else:
3487 return text
3490def main():
3491 import time
3493 times = []
3494 algor = []
3495 x0 = [0.8, 1.2, 0.7]
3496 print("Nelder-Mead Simplex")
3497 print("===================")
3498 start = time.time()
3499 x = fmin(rosen, x0)
3500 print(x)
3501 times.append(time.time() - start)
3502 algor.append('Nelder-Mead Simplex\t')
3504 print()
3505 print("Powell Direction Set Method")
3506 print("===========================")
3507 start = time.time()
3508 x = fmin_powell(rosen, x0)
3509 print(x)
3510 times.append(time.time() - start)
3511 algor.append('Powell Direction Set Method.')
3513 print()
3514 print("Nonlinear CG")
3515 print("============")
3516 start = time.time()
3517 x = fmin_cg(rosen, x0, fprime=rosen_der, maxiter=200)
3518 print(x)
3519 times.append(time.time() - start)
3520 algor.append('Nonlinear CG \t')
3522 print()
3523 print("BFGS Quasi-Newton")
3524 print("=================")
3525 start = time.time()
3526 x = fmin_bfgs(rosen, x0, fprime=rosen_der, maxiter=80)
3527 print(x)
3528 times.append(time.time() - start)
3529 algor.append('BFGS Quasi-Newton\t')
3531 print()
3532 print("BFGS approximate gradient")
3533 print("=========================")
3534 start = time.time()
3535 x = fmin_bfgs(rosen, x0, gtol=1e-4, maxiter=100)
3536 print(x)
3537 times.append(time.time() - start)
3538 algor.append('BFGS without gradient\t')
3540 print()
3541 print("Newton-CG with Hessian product")
3542 print("==============================")
3543 start = time.time()
3544 x = fmin_ncg(rosen, x0, rosen_der, fhess_p=rosen_hess_prod, maxiter=80)
3545 print(x)
3546 times.append(time.time() - start)
3547 algor.append('Newton-CG with hessian product')
3549 print()
3550 print("Newton-CG with full Hessian")
3551 print("===========================")
3552 start = time.time()
3553 x = fmin_ncg(rosen, x0, rosen_der, fhess=rosen_hess, maxiter=80)
3554 print(x)
3555 times.append(time.time() - start)
3556 algor.append('Newton-CG with full Hessian')
3558 print()
3559 print("\nMinimizing the Rosenbrock function of order 3\n")
3560 print(" Algorithm \t\t\t Seconds")
3561 print("===========\t\t\t =========")
3562 for k in range(len(algor)):
3563 print(algor[k], "\t -- ", times[k])
3566if __name__ == "__main__":
3567 main()