import numpy as np from scipy.optimize._optimize import OptimizeResult from scipy.optimize._zeros_py import (_scalar_optimization_initialize, # noqa: F401 _chandrupatla_iv, _scalar_optimization_loop, _ECONVERGED, _ESIGNERR, _ECONVERR, _EVALUEERR, _ECALLBACK, _EINPROGRESS) def _chandrupatla_minimize(func, x1, x2, x3, *, args=(), xatol=None, xrtol=None, fatol=None, frtol=None, maxiter=100, callback=None): """Find the minimizer of an elementwise function. For each element of the output of `func`, `_chandrupatla_minimize` seeks the scalar minimizer that minimizes the element. This function allows for `x1`, `x2`, `x3`, and the elements of `args` to be arrays of any broadcastable shapes. Parameters ---------- func : callable The function whose minimizer is desired. The signature must be:: func(x: ndarray, *args) -> ndarray where each element of ``x`` is a finite real and ``args`` is a tuple, which may contain an arbitrary number of arrays that are broadcastable with `x`. ``func`` must be an elementwise function: each element ``func(x)[i]`` must equal ``func(x[i])`` for all indices ``i``. `_chandrupatla` seeks an array ``x`` such that ``func(x)`` is an array of minima. x1, x2, x3 : array_like The abscissae of a standard scalar minimization bracket. A bracket is valid if ``x1 < x2 < x3`` and ``func(x1) > func(x2) <= func(x3)``. Must be broadcastable with one another and `args`. args : tuple, optional Additional positional arguments to be passed to `func`. Must be arrays broadcastable with `x1`, `x2`, and `x3`. If the callable to be differentiated requires arguments that are not broadcastable with `x`, wrap that callable with `func` such that `func` accepts only `x` and broadcastable arrays. xatol, xrtol, fatol, frtol : float, optional Absolute and relative tolerances on the minimizer and function value. See Notes for details. maxiter : int, optional The maximum number of iterations of the algorithm to perform. callback : callable, optional An optional user-supplied function to be called before the first iteration and after each iteration. Called as ``callback(res)``, where ``res`` is an ``OptimizeResult`` similar to that returned by `_chandrupatla_minimize` (but containing the current iterate's values of all variables). If `callback` raises a ``StopIteration``, the algorithm will terminate immediately and `_chandrupatla_minimize` will return a result. Returns ------- res : OptimizeResult An instance of `scipy.optimize.OptimizeResult` with the following attributes. (The descriptions are written as though the values will be scalars; however, if `func` returns an array, the outputs will be arrays of the same shape.) success : bool ``True`` when the algorithm terminated successfully (status ``0``). status : int An integer representing the exit status of the algorithm. ``0`` : The algorithm converged to the specified tolerances. ``-1`` : The algorithm encountered an invalid bracket. ``-2`` : The maximum number of iterations was reached. ``-3`` : A non-finite value was encountered. ``-4`` : Iteration was terminated by `callback`. ``1`` : The algorithm is proceeding normally (in `callback` only). x : float The minimizer of the function, if the algorithm terminated successfully. fun : float The value of `func` evaluated at `x`. nfev : int The number of points at which `func` was evaluated. nit : int The number of iterations of the algorithm that were performed. xl, xm, xr : float The final three-point bracket. fl, fm, fr : float The function value at the bracket points. Notes ----- Implemented based on Chandrupatla's original paper [1]_. If ``x1 < x2 < x3`` are the points of the bracket and ``f1 > f2 <= f3`` are the values of ``func`` at those points, then the algorithm is considered to have converged when ``x3 - x1 <= abs(x2)*xrtol + xatol`` or ``(f1 - 2*f2 + f3)/2 <= abs(f2)*frtol + fatol``. Note that first of these differs from the termination conditions described in [1]_. The default values of `xrtol` is the square root of the precision of the appropriate dtype, and ``xatol=fatol = frtol`` is the smallest normal number of the appropriate dtype. References ---------- .. [1] Chandrupatla, Tirupathi R. (1998). "An efficient quadratic fit-sectioning algorithm for minimization without derivatives". Computer Methods in Applied Mechanics and Engineering, 152 (1-2), 211-217. https://doi.org/10.1016/S0045-7825(97)00190-4 See Also -------- golden, brent, bounded Examples -------- >>> from scipy.optimize._chandrupatla import _chandrupatla_minimize >>> def f(x, args=1): ... return (x - args)**2 >>> res = _chandrupatla_minimize(f, -5, 0, 5) >>> res.x 1.0 >>> c = [1, 1.5, 2] >>> res = _chandrupatla_minimize(f, -5, 0, 5, args=(c,)) >>> res.x array([1. , 1.5, 2. ]) """ res = _chandrupatla_iv(func, args, xatol, xrtol, fatol, frtol, maxiter, callback) func, args, xatol, xrtol, fatol, frtol, maxiter, callback = res # Initialization xs = (x1, x2, x3) temp = _scalar_optimization_initialize(func, xs, args) xs, fs, args, shape, dtype = temp # line split for PEP8 x1, x2, x3 = xs f1, f2, f3 = fs phi = dtype.type(0.5 + 0.5*5**0.5) # golden ratio status = np.full_like(x1, _EINPROGRESS, dtype=int) # in progress nit, nfev = 0, 3 # three function evaluations performed above fatol = np.finfo(dtype).tiny if fatol is None else fatol frtol = np.finfo(dtype).tiny if frtol is None else frtol xatol = np.finfo(dtype).tiny if xatol is None else xatol xrtol = np.sqrt(np.finfo(dtype).eps) if xrtol is None else xrtol # Ensure that x1 < x2 < x3 initially. xs, fs = np.vstack((x1, x2, x3)), np.vstack((f1, f2, f3)) i = np.argsort(xs, axis=0) x1, x2, x3 = np.take_along_axis(xs, i, axis=0) f1, f2, f3 = np.take_along_axis(fs, i, axis=0) q0 = x3.copy() # "At the start, q0 is set at x3..." ([1] after (7)) work = OptimizeResult(x1=x1, f1=f1, x2=x2, f2=f2, x3=x3, f3=f3, phi=phi, xatol=xatol, xrtol=xrtol, fatol=fatol, frtol=frtol, nit=nit, nfev=nfev, status=status, q0=q0, args=args) res_work_pairs = [('status', 'status'), ('x', 'x2'), ('fun', 'f2'), ('nit', 'nit'), ('nfev', 'nfev'), ('xl', 'x1'), ('xm', 'x2'), ('xr', 'x3'), ('fl', 'f1'), ('fm', 'f2'), ('fr', 'f3')] def pre_func_eval(work): # `_check_termination` is called first -> `x3 - x2 > x2 - x1` # But let's calculate a few terms that we'll reuse x21 = work.x2 - work.x1 x32 = work.x3 - work.x2 # [1] Section 3. "The quadratic minimum point Q1 is calculated using # the relations developed in the previous section." [1] Section 2 (5/6) A = x21 * (work.f3 - work.f2) B = x32 * (work.f1 - work.f2) C = A / (A + B) # q1 = C * (work.x1 + work.x2) / 2 + (1 - C) * (work.x2 + work.x3) / 2 q1 = 0.5 * (C*(work.x1 - work.x3) + work.x2 + work.x3) # much faster # this is an array, so multiplying by 0.5 does not change dtype # "If Q1 and Q0 are sufficiently close... Q1 is accepted if it is # sufficiently away from the inside point x2" i = abs(q1 - work.q0) < 0.5 * abs(x21) # [1] (7) xi = q1[i] # Later, after (9), "If the point Q1 is in a +/- xtol neighborhood of # x2, the new point is chosen in the larger interval at a distance # tol away from x2." # See also QBASIC code after "Accept Ql adjust if close to X2". j = abs(q1[i] - work.x2[i]) <= work.xtol[i] xi[j] = work.x2[i][j] + np.sign(x32[i][j]) * work.xtol[i][j] # "If condition (7) is not satisfied, golden sectioning of the larger # interval is carried out to introduce the new point." # (For simplicity, we go ahead and calculate it for all points, but we # change the elements for which the condition was satisfied.) x = work.x2 + (2 - work.phi) * x32 x[i] = xi # "We define Q0 as the value of Q1 at the previous iteration." work.q0 = q1 return x def post_func_eval(x, f, work): # Standard logic for updating a three-point bracket based on a new # point. In QBASIC code, see "IF SGN(X-X2) = SGN(X3-X2) THEN...". # There is an awful lot of data copying going on here; this would # probably benefit from code optimization or implementation in Pythran. i = np.sign(x - work.x2) == np.sign(work.x3 - work.x2) xi, x1i, x2i, x3i = x[i], work.x1[i], work.x2[i], work.x3[i], fi, f1i, f2i, f3i = f[i], work.f1[i], work.f2[i], work.f3[i] j = fi > f2i x3i[j], f3i[j] = xi[j], fi[j] j = ~j x1i[j], f1i[j], x2i[j], f2i[j] = x2i[j], f2i[j], xi[j], fi[j] ni = ~i xni, x1ni, x2ni, x3ni = x[ni], work.x1[ni], work.x2[ni], work.x3[ni], fni, f1ni, f2ni, f3ni = f[ni], work.f1[ni], work.f2[ni], work.f3[ni] j = fni > f2ni x1ni[j], f1ni[j] = xni[j], fni[j] j = ~j x3ni[j], f3ni[j], x2ni[j], f2ni[j] = x2ni[j], f2ni[j], xni[j], fni[j] work.x1[i], work.x2[i], work.x3[i] = x1i, x2i, x3i work.f1[i], work.f2[i], work.f3[i] = f1i, f2i, f3i work.x1[ni], work.x2[ni], work.x3[ni] = x1ni, x2ni, x3ni, work.f1[ni], work.f2[ni], work.f3[ni] = f1ni, f2ni, f3ni def check_termination(work): # Check for all terminal conditions and record statuses. stop = np.zeros_like(work.x1, dtype=bool) # termination condition met # Bracket is invalid; stop and don't return minimizer/minimum i = ((work.f2 > work.f1) | (work.f2 > work.f3)) work.x2[i], work.f2[i] = np.nan, np.nan stop[i], work.status[i] = True, _ESIGNERR # Non-finite values; stop and don't return minimizer/minimum finite = np.isfinite(work.x1+work.x2+work.x3+work.f1+work.f2+work.f3) i = ~(finite | stop) work.x2[i], work.f2[i] = np.nan, np.nan stop[i], work.status[i] = True, _EVALUEERR # [1] Section 3 "Points 1 and 3 are interchanged if necessary to make # the (x2, x3) the larger interval." # Note: I had used np.choose; this is much faster. This would be a good # place to save e.g. `work.x3 - work.x2` for reuse, but I tried and # didn't notice a speed boost, so let's keep it simple. i = abs(work.x3 - work.x2) < abs(work.x2 - work.x1) temp = work.x1[i] work.x1[i] = work.x3[i] work.x3[i] = temp temp = work.f1[i] work.f1[i] = work.f3[i] work.f3[i] = temp # [1] Section 3 (bottom of page 212) # "We set a tolerance value xtol..." work.xtol = abs(work.x2) * work.xrtol + work.xatol # [1] (8) # "The convergence based on interval is achieved when..." # Note: Equality allowed in case of `xtol=0` i = abs(work.x3 - work.x2) <= 2 * work.xtol # [1] (9) # "We define ftol using..." ftol = abs(work.f2) * work.frtol + work.fatol # [1] (10) # "The convergence based on function values is achieved when..." # Note 1: modify in place to incorporate tolerance on function value. # Note 2: factor of 2 is not in the text; see QBASIC start of DO loop i |= (work.f1 - 2 * work.f2 + work.f3) <= 2*ftol # [1] (11) i &= ~stop stop[i], work.status[i] = True, _ECONVERGED return stop def post_termination_check(work): pass def customize_result(res, shape): xl, xr, fl, fr = res['xl'], res['xr'], res['fl'], res['fr'] i = res['xl'] < res['xr'] res['xl'] = np.choose(i, (xr, xl)) res['xr'] = np.choose(i, (xl, xr)) res['fl'] = np.choose(i, (fr, fl)) res['fr'] = np.choose(i, (fl, fr)) return shape return _scalar_optimization_loop(work, callback, shape, maxiter, func, args, dtype, pre_func_eval, post_func_eval, check_termination, post_termination_check, customize_result, res_work_pairs)