Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_slsqp_py.py: 8%
168 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""
2This module implements the Sequential Least Squares Programming optimization
3algorithm (SLSQP), originally developed by Dieter Kraft.
4See http://www.netlib.org/toms/733
6Functions
7---------
8.. autosummary::
9 :toctree: generated/
11 approx_jacobian
12 fmin_slsqp
14"""
16__all__ = ['approx_jacobian', 'fmin_slsqp']
18import numpy as np
19from scipy.optimize._slsqp import slsqp
20from numpy import (zeros, array, linalg, append, asfarray, concatenate, finfo,
21 sqrt, vstack, isfinite, atleast_1d)
22from ._optimize import (OptimizeResult, _check_unknown_options,
23 _prepare_scalar_function, _clip_x_for_func,
24 _check_clip_x)
25from ._numdiff import approx_derivative
26from ._constraints import old_bound_to_new, _arr_to_scalar
29__docformat__ = "restructuredtext en"
31_epsilon = sqrt(finfo(float).eps)
34def approx_jacobian(x, func, epsilon, *args):
35 """
36 Approximate the Jacobian matrix of a callable function.
38 Parameters
39 ----------
40 x : array_like
41 The state vector at which to compute the Jacobian matrix.
42 func : callable f(x,*args)
43 The vector-valued function.
44 epsilon : float
45 The perturbation used to determine the partial derivatives.
46 args : sequence
47 Additional arguments passed to func.
49 Returns
50 -------
51 An array of dimensions ``(lenf, lenx)`` where ``lenf`` is the length
52 of the outputs of `func`, and ``lenx`` is the number of elements in
53 `x`.
55 Notes
56 -----
57 The approximation is done using forward differences.
59 """
60 # approx_derivative returns (m, n) == (lenf, lenx)
61 jac = approx_derivative(func, x, method='2-point', abs_step=epsilon,
62 args=args)
63 # if func returns a scalar jac.shape will be (lenx,). Make sure
64 # it's at least a 2D array.
65 return np.atleast_2d(jac)
68def fmin_slsqp(func, x0, eqcons=(), f_eqcons=None, ieqcons=(), f_ieqcons=None,
69 bounds=(), fprime=None, fprime_eqcons=None,
70 fprime_ieqcons=None, args=(), iter=100, acc=1.0E-6,
71 iprint=1, disp=None, full_output=0, epsilon=_epsilon,
72 callback=None):
73 """
74 Minimize a function using Sequential Least Squares Programming
76 Python interface function for the SLSQP Optimization subroutine
77 originally implemented by Dieter Kraft.
79 Parameters
80 ----------
81 func : callable f(x,*args)
82 Objective function. Must return a scalar.
83 x0 : 1-D ndarray of float
84 Initial guess for the independent variable(s).
85 eqcons : list, optional
86 A list of functions of length n such that
87 eqcons[j](x,*args) == 0.0 in a successfully optimized
88 problem.
89 f_eqcons : callable f(x,*args), optional
90 Returns a 1-D array in which each element must equal 0.0 in a
91 successfully optimized problem. If f_eqcons is specified,
92 eqcons is ignored.
93 ieqcons : list, optional
94 A list of functions of length n such that
95 ieqcons[j](x,*args) >= 0.0 in a successfully optimized
96 problem.
97 f_ieqcons : callable f(x,*args), optional
98 Returns a 1-D ndarray in which each element must be greater or
99 equal to 0.0 in a successfully optimized problem. If
100 f_ieqcons is specified, ieqcons is ignored.
101 bounds : list, optional
102 A list of tuples specifying the lower and upper bound
103 for each independent variable [(xl0, xu0),(xl1, xu1),...]
104 Infinite values will be interpreted as large floating values.
105 fprime : callable `f(x,*args)`, optional
106 A function that evaluates the partial derivatives of func.
107 fprime_eqcons : callable `f(x,*args)`, optional
108 A function of the form `f(x, *args)` that returns the m by n
109 array of equality constraint normals. If not provided,
110 the normals will be approximated. The array returned by
111 fprime_eqcons should be sized as ( len(eqcons), len(x0) ).
112 fprime_ieqcons : callable `f(x,*args)`, optional
113 A function of the form `f(x, *args)` that returns the m by n
114 array of inequality constraint normals. If not provided,
115 the normals will be approximated. The array returned by
116 fprime_ieqcons should be sized as ( len(ieqcons), len(x0) ).
117 args : sequence, optional
118 Additional arguments passed to func and fprime.
119 iter : int, optional
120 The maximum number of iterations.
121 acc : float, optional
122 Requested accuracy.
123 iprint : int, optional
124 The verbosity of fmin_slsqp :
126 * iprint <= 0 : Silent operation
127 * iprint == 1 : Print summary upon completion (default)
128 * iprint >= 2 : Print status of each iterate and summary
129 disp : int, optional
130 Overrides the iprint interface (preferred).
131 full_output : bool, optional
132 If False, return only the minimizer of func (default).
133 Otherwise, output final objective function and summary
134 information.
135 epsilon : float, optional
136 The step size for finite-difference derivative estimates.
137 callback : callable, optional
138 Called after each iteration, as ``callback(x)``, where ``x`` is the
139 current parameter vector.
141 Returns
142 -------
143 out : ndarray of float
144 The final minimizer of func.
145 fx : ndarray of float, if full_output is true
146 The final value of the objective function.
147 its : int, if full_output is true
148 The number of iterations.
149 imode : int, if full_output is true
150 The exit mode from the optimizer (see below).
151 smode : string, if full_output is true
152 Message describing the exit mode from the optimizer.
154 See also
155 --------
156 minimize: Interface to minimization algorithms for multivariate
157 functions. See the 'SLSQP' `method` in particular.
159 Notes
160 -----
161 Exit modes are defined as follows ::
163 -1 : Gradient evaluation required (g & a)
164 0 : Optimization terminated successfully
165 1 : Function evaluation required (f & c)
166 2 : More equality constraints than independent variables
167 3 : More than 3*n iterations in LSQ subproblem
168 4 : Inequality constraints incompatible
169 5 : Singular matrix E in LSQ subproblem
170 6 : Singular matrix C in LSQ subproblem
171 7 : Rank-deficient equality constraint subproblem HFTI
172 8 : Positive directional derivative for linesearch
173 9 : Iteration limit reached
175 Examples
176 --------
177 Examples are given :ref:`in the tutorial <tutorial-sqlsp>`.
179 """
180 if disp is not None:
181 iprint = disp
183 opts = {'maxiter': iter,
184 'ftol': acc,
185 'iprint': iprint,
186 'disp': iprint != 0,
187 'eps': epsilon,
188 'callback': callback}
190 # Build the constraints as a tuple of dictionaries
191 cons = ()
192 # 1. constraints of the 1st kind (eqcons, ieqcons); no Jacobian; take
193 # the same extra arguments as the objective function.
194 cons += tuple({'type': 'eq', 'fun': c, 'args': args} for c in eqcons)
195 cons += tuple({'type': 'ineq', 'fun': c, 'args': args} for c in ieqcons)
196 # 2. constraints of the 2nd kind (f_eqcons, f_ieqcons) and their Jacobian
197 # (fprime_eqcons, fprime_ieqcons); also take the same extra arguments
198 # as the objective function.
199 if f_eqcons:
200 cons += ({'type': 'eq', 'fun': f_eqcons, 'jac': fprime_eqcons,
201 'args': args}, )
202 if f_ieqcons:
203 cons += ({'type': 'ineq', 'fun': f_ieqcons, 'jac': fprime_ieqcons,
204 'args': args}, )
206 res = _minimize_slsqp(func, x0, args, jac=fprime, bounds=bounds,
207 constraints=cons, **opts)
208 if full_output:
209 return res['x'], res['fun'], res['nit'], res['status'], res['message']
210 else:
211 return res['x']
214def _minimize_slsqp(func, x0, args=(), jac=None, bounds=None,
215 constraints=(),
216 maxiter=100, ftol=1.0E-6, iprint=1, disp=False,
217 eps=_epsilon, callback=None, finite_diff_rel_step=None,
218 **unknown_options):
219 """
220 Minimize a scalar function of one or more variables using Sequential
221 Least Squares Programming (SLSQP).
223 Options
224 -------
225 ftol : float
226 Precision goal for the value of f in the stopping criterion.
227 eps : float
228 Step size used for numerical approximation of the Jacobian.
229 disp : bool
230 Set to True to print convergence messages. If False,
231 `verbosity` is ignored and set to 0.
232 maxiter : int
233 Maximum number of iterations.
234 finite_diff_rel_step : None or array_like, optional
235 If `jac in ['2-point', '3-point', 'cs']` the relative step size to
236 use for numerical approximation of `jac`. The absolute step
237 size is computed as ``h = rel_step * sign(x) * max(1, abs(x))``,
238 possibly adjusted to fit into the bounds. For ``method='3-point'``
239 the sign of `h` is ignored. If None (default) then step is selected
240 automatically.
241 """
242 _check_unknown_options(unknown_options)
243 iter = maxiter - 1
244 acc = ftol
245 epsilon = eps
247 if not disp:
248 iprint = 0
250 # Transform x0 into an array.
251 x = asfarray(x0).flatten()
253 # SLSQP is sent 'old-style' bounds, 'new-style' bounds are required by
254 # ScalarFunction
255 if bounds is None or len(bounds) == 0:
256 new_bounds = (-np.inf, np.inf)
257 else:
258 new_bounds = old_bound_to_new(bounds)
260 # clip the initial guess to bounds, otherwise ScalarFunction doesn't work
261 x = np.clip(x, new_bounds[0], new_bounds[1])
263 # Constraints are triaged per type into a dictionary of tuples
264 if isinstance(constraints, dict):
265 constraints = (constraints, )
267 cons = {'eq': (), 'ineq': ()}
268 for ic, con in enumerate(constraints):
269 # check type
270 try:
271 ctype = con['type'].lower()
272 except KeyError as e:
273 raise KeyError('Constraint %d has no type defined.' % ic) from e
274 except TypeError as e:
275 raise TypeError('Constraints must be defined using a '
276 'dictionary.') from e
277 except AttributeError as e:
278 raise TypeError("Constraint's type must be a string.") from e
279 else:
280 if ctype not in ['eq', 'ineq']:
281 raise ValueError("Unknown constraint type '%s'." % con['type'])
283 # check function
284 if 'fun' not in con:
285 raise ValueError('Constraint %d has no function defined.' % ic)
287 # check Jacobian
288 cjac = con.get('jac')
289 if cjac is None:
290 # approximate Jacobian function. The factory function is needed
291 # to keep a reference to `fun`, see gh-4240.
292 def cjac_factory(fun):
293 def cjac(x, *args):
294 x = _check_clip_x(x, new_bounds)
296 if jac in ['2-point', '3-point', 'cs']:
297 return approx_derivative(fun, x, method=jac, args=args,
298 rel_step=finite_diff_rel_step,
299 bounds=new_bounds)
300 else:
301 return approx_derivative(fun, x, method='2-point',
302 abs_step=epsilon, args=args,
303 bounds=new_bounds)
305 return cjac
306 cjac = cjac_factory(con['fun'])
308 # update constraints' dictionary
309 cons[ctype] += ({'fun': con['fun'],
310 'jac': cjac,
311 'args': con.get('args', ())}, )
313 exit_modes = {-1: "Gradient evaluation required (g & a)",
314 0: "Optimization terminated successfully",
315 1: "Function evaluation required (f & c)",
316 2: "More equality constraints than independent variables",
317 3: "More than 3*n iterations in LSQ subproblem",
318 4: "Inequality constraints incompatible",
319 5: "Singular matrix E in LSQ subproblem",
320 6: "Singular matrix C in LSQ subproblem",
321 7: "Rank-deficient equality constraint subproblem HFTI",
322 8: "Positive directional derivative for linesearch",
323 9: "Iteration limit reached"}
325 # Set the parameters that SLSQP will need
326 # meq, mieq: number of equality and inequality constraints
327 meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
328 for c in cons['eq']]))
329 mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
330 for c in cons['ineq']]))
331 # m = The total number of constraints
332 m = meq + mieq
333 # la = The number of constraints, or 1 if there are no constraints
334 la = array([1, m]).max()
335 # n = The number of independent variables
336 n = len(x)
338 # Define the workspaces for SLSQP
339 n1 = n + 1
340 mineq = m - meq + n1 + n1
341 len_w = (3*n1+m)*(n1+1)+(n1-meq+1)*(mineq+2) + 2*mineq+(n1+mineq)*(n1-meq) \
342 + 2*meq + n1 + ((n+1)*n)//2 + 2*m + 3*n + 3*n1 + 1
343 len_jw = mineq
344 w = zeros(len_w)
345 jw = zeros(len_jw)
347 # Decompose bounds into xl and xu
348 if bounds is None or len(bounds) == 0:
349 xl = np.empty(n, dtype=float)
350 xu = np.empty(n, dtype=float)
351 xl.fill(np.nan)
352 xu.fill(np.nan)
353 else:
354 bnds = array([(_arr_to_scalar(l), _arr_to_scalar(u))
355 for (l, u) in bounds], float)
356 if bnds.shape[0] != n:
357 raise IndexError('SLSQP Error: the length of bounds is not '
358 'compatible with that of x0.')
360 with np.errstate(invalid='ignore'):
361 bnderr = bnds[:, 0] > bnds[:, 1]
363 if bnderr.any():
364 raise ValueError('SLSQP Error: lb > ub in bounds %s.' %
365 ', '.join(str(b) for b in bnderr))
366 xl, xu = bnds[:, 0], bnds[:, 1]
368 # Mark infinite bounds with nans; the Fortran code understands this
369 infbnd = ~isfinite(bnds)
370 xl[infbnd[:, 0]] = np.nan
371 xu[infbnd[:, 1]] = np.nan
373 # ScalarFunction provides function and gradient evaluation
374 sf = _prepare_scalar_function(func, x, jac=jac, args=args, epsilon=eps,
375 finite_diff_rel_step=finite_diff_rel_step,
376 bounds=new_bounds)
377 # gh11403 SLSQP sometimes exceeds bounds by 1 or 2 ULP, make sure this
378 # doesn't get sent to the func/grad evaluator.
379 wrapped_fun = _clip_x_for_func(sf.fun, new_bounds)
380 wrapped_grad = _clip_x_for_func(sf.grad, new_bounds)
382 # Initialize the iteration counter and the mode value
383 mode = array(0, int)
384 acc = array(acc, float)
385 majiter = array(iter, int)
386 majiter_prev = 0
388 # Initialize internal SLSQP state variables
389 alpha = array(0, float)
390 f0 = array(0, float)
391 gs = array(0, float)
392 h1 = array(0, float)
393 h2 = array(0, float)
394 h3 = array(0, float)
395 h4 = array(0, float)
396 t = array(0, float)
397 t0 = array(0, float)
398 tol = array(0, float)
399 iexact = array(0, int)
400 incons = array(0, int)
401 ireset = array(0, int)
402 itermx = array(0, int)
403 line = array(0, int)
404 n1 = array(0, int)
405 n2 = array(0, int)
406 n3 = array(0, int)
408 # Print the header if iprint >= 2
409 if iprint >= 2:
410 print("%5s %5s %16s %16s" % ("NIT", "FC", "OBJFUN", "GNORM"))
412 # mode is zero on entry, so call objective, constraints and gradients
413 # there should be no func evaluations here because it's cached from
414 # ScalarFunction
415 fx = wrapped_fun(x)
416 g = append(wrapped_grad(x), 0.0)
417 c = _eval_constraint(x, cons)
418 a = _eval_con_normals(x, cons, la, n, m, meq, mieq)
420 while 1:
421 # Call SLSQP
422 slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw,
423 alpha, f0, gs, h1, h2, h3, h4, t, t0, tol,
424 iexact, incons, ireset, itermx, line,
425 n1, n2, n3)
427 if mode == 1: # objective and constraint evaluation required
428 fx = wrapped_fun(x)
429 c = _eval_constraint(x, cons)
431 if mode == -1: # gradient evaluation required
432 g = append(wrapped_grad(x), 0.0)
433 a = _eval_con_normals(x, cons, la, n, m, meq, mieq)
435 if majiter > majiter_prev:
436 # call callback if major iteration has incremented
437 if callback is not None:
438 callback(np.copy(x))
440 # Print the status of the current iterate if iprint > 2
441 if iprint >= 2:
442 print("%5i %5i % 16.6E % 16.6E" % (majiter, sf.nfev,
443 fx, linalg.norm(g)))
445 # If exit mode is not -1 or 1, slsqp has completed
446 if abs(mode) != 1:
447 break
449 majiter_prev = int(majiter)
451 # Optimization loop complete. Print status if requested
452 if iprint >= 1:
453 print(exit_modes[int(mode)] + " (Exit mode " + str(mode) + ')')
454 print(" Current function value:", fx)
455 print(" Iterations:", majiter)
456 print(" Function evaluations:", sf.nfev)
457 print(" Gradient evaluations:", sf.ngev)
459 return OptimizeResult(x=x, fun=fx, jac=g[:-1], nit=int(majiter),
460 nfev=sf.nfev, njev=sf.ngev, status=int(mode),
461 message=exit_modes[int(mode)], success=(mode == 0))
464def _eval_constraint(x, cons):
465 # Compute constraints
466 if cons['eq']:
467 c_eq = concatenate([atleast_1d(con['fun'](x, *con['args']))
468 for con in cons['eq']])
469 else:
470 c_eq = zeros(0)
472 if cons['ineq']:
473 c_ieq = concatenate([atleast_1d(con['fun'](x, *con['args']))
474 for con in cons['ineq']])
475 else:
476 c_ieq = zeros(0)
478 # Now combine c_eq and c_ieq into a single matrix
479 c = concatenate((c_eq, c_ieq))
480 return c
483def _eval_con_normals(x, cons, la, n, m, meq, mieq):
484 # Compute the normals of the constraints
485 if cons['eq']:
486 a_eq = vstack([con['jac'](x, *con['args'])
487 for con in cons['eq']])
488 else: # no equality constraint
489 a_eq = zeros((meq, n))
491 if cons['ineq']:
492 a_ieq = vstack([con['jac'](x, *con['args'])
493 for con in cons['ineq']])
494 else: # no inequality constraint
495 a_ieq = zeros((mieq, n))
497 # Now combine a_eq and a_ieq into a single a matrix
498 if m == 0: # no constraints
499 a = zeros((la, n))
500 else:
501 a = vstack((a_eq, a_ieq))
502 a = concatenate((a, zeros([la, 1])), 1)
504 return a