Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_constraints.py: 13%
205 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Constraints definition for minimize."""
2import numpy as np
3from ._hessian_update_strategy import BFGS
4from ._differentiable_functions import (
5 VectorFunction, LinearVectorFunction, IdentityVectorFunction)
6from ._optimize import OptimizeWarning
7from warnings import warn, catch_warnings, simplefilter
8from numpy.testing import suppress_warnings
9from scipy.sparse import issparse
12def _arr_to_scalar(x):
13 # If x is a numpy array, return x.item(). This will
14 # fail if the array has more than one element.
15 return x.item() if isinstance(x, np.ndarray) else x
18class NonlinearConstraint:
19 """Nonlinear constraint on the variables.
21 The constraint has the general inequality form::
23 lb <= fun(x) <= ub
25 Here the vector of independent variables x is passed as ndarray of shape
26 (n,) and ``fun`` returns a vector with m components.
28 It is possible to use equal bounds to represent an equality constraint or
29 infinite bounds to represent a one-sided constraint.
31 Parameters
32 ----------
33 fun : callable
34 The function defining the constraint.
35 The signature is ``fun(x) -> array_like, shape (m,)``.
36 lb, ub : array_like
37 Lower and upper bounds on the constraint. Each array must have the
38 shape (m,) or be a scalar, in the latter case a bound will be the same
39 for all components of the constraint. Use ``np.inf`` with an
40 appropriate sign to specify a one-sided constraint.
41 Set components of `lb` and `ub` equal to represent an equality
42 constraint. Note that you can mix constraints of different types:
43 interval, one-sided or equality, by setting different components of
44 `lb` and `ub` as necessary.
45 jac : {callable, '2-point', '3-point', 'cs'}, optional
46 Method of computing the Jacobian matrix (an m-by-n matrix,
47 where element (i, j) is the partial derivative of f[i] with
48 respect to x[j]). The keywords {'2-point', '3-point',
49 'cs'} select a finite difference scheme for the numerical estimation.
50 A callable must have the following signature:
51 ``jac(x) -> {ndarray, sparse matrix}, shape (m, n)``.
52 Default is '2-point'.
53 hess : {callable, '2-point', '3-point', 'cs', HessianUpdateStrategy, None}, optional
54 Method for computing the Hessian matrix. The keywords
55 {'2-point', '3-point', 'cs'} select a finite difference scheme for
56 numerical estimation. Alternatively, objects implementing
57 `HessianUpdateStrategy` interface can be used to approximate the
58 Hessian. Currently available implementations are:
60 - `BFGS` (default option)
61 - `SR1`
63 A callable must return the Hessian matrix of ``dot(fun, v)`` and
64 must have the following signature:
65 ``hess(x, v) -> {LinearOperator, sparse matrix, array_like}, shape (n, n)``.
66 Here ``v`` is ndarray with shape (m,) containing Lagrange multipliers.
67 keep_feasible : array_like of bool, optional
68 Whether to keep the constraint components feasible throughout
69 iterations. A single value set this property for all components.
70 Default is False. Has no effect for equality constraints.
71 finite_diff_rel_step: None or array_like, optional
72 Relative step size for the finite difference approximation. Default is
73 None, which will select a reasonable value automatically depending
74 on a finite difference scheme.
75 finite_diff_jac_sparsity: {None, array_like, sparse matrix}, optional
76 Defines the sparsity structure of the Jacobian matrix for finite
77 difference estimation, its shape must be (m, n). If the Jacobian has
78 only few non-zero elements in *each* row, providing the sparsity
79 structure will greatly speed up the computations. A zero entry means
80 that a corresponding element in the Jacobian is identically zero.
81 If provided, forces the use of 'lsmr' trust-region solver.
82 If None (default) then dense differencing will be used.
84 Notes
85 -----
86 Finite difference schemes {'2-point', '3-point', 'cs'} may be used for
87 approximating either the Jacobian or the Hessian. We, however, do not allow
88 its use for approximating both simultaneously. Hence whenever the Jacobian
89 is estimated via finite-differences, we require the Hessian to be estimated
90 using one of the quasi-Newton strategies.
92 The scheme 'cs' is potentially the most accurate, but requires the function
93 to correctly handles complex inputs and be analytically continuable to the
94 complex plane. The scheme '3-point' is more accurate than '2-point' but
95 requires twice as many operations.
97 Examples
98 --------
99 Constrain ``x[0] < sin(x[1]) + 1.9``
101 >>> from scipy.optimize import NonlinearConstraint
102 >>> import numpy as np
103 >>> con = lambda x: x[0] - np.sin(x[1])
104 >>> nlc = NonlinearConstraint(con, -np.inf, 1.9)
106 """
107 def __init__(self, fun, lb, ub, jac='2-point', hess=BFGS(),
108 keep_feasible=False, finite_diff_rel_step=None,
109 finite_diff_jac_sparsity=None):
110 self.fun = fun
111 self.lb = lb
112 self.ub = ub
113 self.finite_diff_rel_step = finite_diff_rel_step
114 self.finite_diff_jac_sparsity = finite_diff_jac_sparsity
115 self.jac = jac
116 self.hess = hess
117 self.keep_feasible = keep_feasible
120class LinearConstraint:
121 """Linear constraint on the variables.
123 The constraint has the general inequality form::
125 lb <= A.dot(x) <= ub
127 Here the vector of independent variables x is passed as ndarray of shape
128 (n,) and the matrix A has shape (m, n).
130 It is possible to use equal bounds to represent an equality constraint or
131 infinite bounds to represent a one-sided constraint.
133 Parameters
134 ----------
135 A : {array_like, sparse matrix}, shape (m, n)
136 Matrix defining the constraint.
137 lb, ub : array_like, optional
138 Lower and upper limits on the constraint. Each array must have the
139 shape (m,) or be a scalar, in the latter case a bound will be the same
140 for all components of the constraint. Use ``np.inf`` with an
141 appropriate sign to specify a one-sided constraint.
142 Set components of `lb` and `ub` equal to represent an equality
143 constraint. Note that you can mix constraints of different types:
144 interval, one-sided or equality, by setting different components of
145 `lb` and `ub` as necessary. Defaults to ``lb = -np.inf``
146 and ``ub = np.inf`` (no limits).
147 keep_feasible : array_like of bool, optional
148 Whether to keep the constraint components feasible throughout
149 iterations. A single value set this property for all components.
150 Default is False. Has no effect for equality constraints.
151 """
152 def _input_validation(self):
153 if self.A.ndim != 2:
154 message = "`A` must have exactly two dimensions."
155 raise ValueError(message)
157 try:
158 shape = self.A.shape[0:1]
159 self.lb = np.broadcast_to(self.lb, shape)
160 self.ub = np.broadcast_to(self.ub, shape)
161 self.keep_feasible = np.broadcast_to(self.keep_feasible, shape)
162 except ValueError:
163 message = ("`lb`, `ub`, and `keep_feasible` must be broadcastable "
164 "to shape `A.shape[0:1]`")
165 raise ValueError(message)
167 def __init__(self, A, lb=-np.inf, ub=np.inf, keep_feasible=False):
168 if not issparse(A):
169 # In some cases, if the constraint is not valid, this emits a
170 # VisibleDeprecationWarning about ragged nested sequences
171 # before eventually causing an error. `scipy.optimize.milp` would
172 # prefer that this just error out immediately so it can handle it
173 # rather than concerning the user.
174 with catch_warnings():
175 simplefilter("error")
176 self.A = np.atleast_2d(A).astype(np.float64)
177 else:
178 self.A = A
179 self.lb = np.atleast_1d(lb).astype(np.float64)
180 self.ub = np.atleast_1d(ub).astype(np.float64)
181 self.keep_feasible = np.atleast_1d(keep_feasible).astype(bool)
182 self._input_validation()
184 def residual(self, x):
185 """
186 Calculate the residual between the constraint function and the limits
188 For a linear constraint of the form::
190 lb <= A@x <= ub
192 the lower and upper residuals between ``A@x`` and the limits are values
193 ``sl`` and ``sb`` such that::
195 lb + sl == A@x == ub - sb
197 When all elements of ``sl`` and ``sb`` are positive, all elements of
198 the constraint are satisfied; a negative element in ``sl`` or ``sb``
199 indicates that the corresponding element of the constraint is not
200 satisfied.
202 Parameters
203 ----------
204 x: array_like
205 Vector of independent variables
207 Returns
208 -------
209 sl, sb : array-like
210 The lower and upper residuals
211 """
212 return self.A@x - self.lb, self.ub - self.A@x
215class Bounds:
216 """Bounds constraint on the variables.
218 The constraint has the general inequality form::
220 lb <= x <= ub
222 It is possible to use equal bounds to represent an equality constraint or
223 infinite bounds to represent a one-sided constraint.
225 Parameters
226 ----------
227 lb, ub : array_like, optional
228 Lower and upper bounds on independent variables. `lb`, `ub`, and
229 `keep_feasible` must be the same shape or broadcastable.
230 Set components of `lb` and `ub` equal
231 to fix a variable. Use ``np.inf`` with an appropriate sign to disable
232 bounds on all or some variables. Note that you can mix constraints of
233 different types: interval, one-sided or equality, by setting different
234 components of `lb` and `ub` as necessary. Defaults to ``lb = -np.inf``
235 and ``ub = np.inf`` (no bounds).
236 keep_feasible : array_like of bool, optional
237 Whether to keep the constraint components feasible throughout
238 iterations. Must be broadcastable with `lb` and `ub`.
239 Default is False. Has no effect for equality constraints.
240 """
241 def _input_validation(self):
242 try:
243 res = np.broadcast_arrays(self.lb, self.ub, self.keep_feasible)
244 self.lb, self.ub, self.keep_feasible = res
245 except ValueError:
246 message = "`lb`, `ub`, and `keep_feasible` must be broadcastable."
247 raise ValueError(message)
249 def __init__(self, lb=-np.inf, ub=np.inf, keep_feasible=False):
250 self.lb = np.atleast_1d(lb)
251 self.ub = np.atleast_1d(ub)
252 self.keep_feasible = np.atleast_1d(keep_feasible).astype(bool)
253 self._input_validation()
255 def __repr__(self):
256 start = f"{type(self).__name__}({self.lb!r}, {self.ub!r}"
257 if np.any(self.keep_feasible):
258 end = f", keep_feasible={self.keep_feasible!r})"
259 else:
260 end = ")"
261 return start + end
263 def residual(self, x):
264 """Calculate the residual (slack) between the input and the bounds
266 For a bound constraint of the form::
268 lb <= x <= ub
270 the lower and upper residuals between `x` and the bounds are values
271 ``sl`` and ``sb`` such that::
273 lb + sl == x == ub - sb
275 When all elements of ``sl`` and ``sb`` are positive, all elements of
276 ``x`` lie within the bounds; a negative element in ``sl`` or ``sb``
277 indicates that the corresponding element of ``x`` is out of bounds.
279 Parameters
280 ----------
281 x: array_like
282 Vector of independent variables
284 Returns
285 -------
286 sl, sb : array-like
287 The lower and upper residuals
288 """
289 return x - self.lb, self.ub - x
292class PreparedConstraint:
293 """Constraint prepared from a user defined constraint.
295 On creation it will check whether a constraint definition is valid and
296 the initial point is feasible. If created successfully, it will contain
297 the attributes listed below.
299 Parameters
300 ----------
301 constraint : {NonlinearConstraint, LinearConstraint`, Bounds}
302 Constraint to check and prepare.
303 x0 : array_like
304 Initial vector of independent variables.
305 sparse_jacobian : bool or None, optional
306 If bool, then the Jacobian of the constraint will be converted
307 to the corresponded format if necessary. If None (default), such
308 conversion is not made.
309 finite_diff_bounds : 2-tuple, optional
310 Lower and upper bounds on the independent variables for the finite
311 difference approximation, if applicable. Defaults to no bounds.
313 Attributes
314 ----------
315 fun : {VectorFunction, LinearVectorFunction, IdentityVectorFunction}
316 Function defining the constraint wrapped by one of the convenience
317 classes.
318 bounds : 2-tuple
319 Contains lower and upper bounds for the constraints --- lb and ub.
320 These are converted to ndarray and have a size equal to the number of
321 the constraints.
322 keep_feasible : ndarray
323 Array indicating which components must be kept feasible with a size
324 equal to the number of the constraints.
325 """
326 def __init__(self, constraint, x0, sparse_jacobian=None,
327 finite_diff_bounds=(-np.inf, np.inf)):
328 if isinstance(constraint, NonlinearConstraint):
329 fun = VectorFunction(constraint.fun, x0,
330 constraint.jac, constraint.hess,
331 constraint.finite_diff_rel_step,
332 constraint.finite_diff_jac_sparsity,
333 finite_diff_bounds, sparse_jacobian)
334 elif isinstance(constraint, LinearConstraint):
335 fun = LinearVectorFunction(constraint.A, x0, sparse_jacobian)
336 elif isinstance(constraint, Bounds):
337 fun = IdentityVectorFunction(x0, sparse_jacobian)
338 else:
339 raise ValueError("`constraint` of an unknown type is passed.")
341 m = fun.m
343 lb = np.asarray(constraint.lb, dtype=float)
344 ub = np.asarray(constraint.ub, dtype=float)
345 keep_feasible = np.asarray(constraint.keep_feasible, dtype=bool)
347 lb = np.broadcast_to(lb, m)
348 ub = np.broadcast_to(ub, m)
349 keep_feasible = np.broadcast_to(keep_feasible, m)
351 if keep_feasible.shape != (m,):
352 raise ValueError("`keep_feasible` has a wrong shape.")
354 mask = keep_feasible & (lb != ub)
355 f0 = fun.f
356 if np.any(f0[mask] < lb[mask]) or np.any(f0[mask] > ub[mask]):
357 raise ValueError("`x0` is infeasible with respect to some "
358 "inequality constraint with `keep_feasible` "
359 "set to True.")
361 self.fun = fun
362 self.bounds = (lb, ub)
363 self.keep_feasible = keep_feasible
365 def violation(self, x):
366 """How much the constraint is exceeded by.
368 Parameters
369 ----------
370 x : array-like
371 Vector of independent variables
373 Returns
374 -------
375 excess : array-like
376 How much the constraint is exceeded by, for each of the
377 constraints specified by `PreparedConstraint.fun`.
378 """
379 with suppress_warnings() as sup:
380 sup.filter(UserWarning)
381 ev = self.fun.fun(np.asarray(x))
383 excess_lb = np.maximum(self.bounds[0] - ev, 0)
384 excess_ub = np.maximum(ev - self.bounds[1], 0)
386 return excess_lb + excess_ub
389def new_bounds_to_old(lb, ub, n):
390 """Convert the new bounds representation to the old one.
392 The new representation is a tuple (lb, ub) and the old one is a list
393 containing n tuples, ith containing lower and upper bound on a ith
394 variable.
395 If any of the entries in lb/ub are -np.inf/np.inf they are replaced by
396 None.
397 """
398 lb = np.broadcast_to(lb, n)
399 ub = np.broadcast_to(ub, n)
401 lb = [float(x) if x > -np.inf else None for x in lb]
402 ub = [float(x) if x < np.inf else None for x in ub]
404 return list(zip(lb, ub))
407def old_bound_to_new(bounds):
408 """Convert the old bounds representation to the new one.
410 The new representation is a tuple (lb, ub) and the old one is a list
411 containing n tuples, ith containing lower and upper bound on a ith
412 variable.
413 If any of the entries in lb/ub are None they are replaced by
414 -np.inf/np.inf.
415 """
416 lb, ub = zip(*bounds)
418 # Convert occurrences of None to -inf or inf, and replace occurrences of
419 # any numpy array x with x.item(). Then wrap the results in numpy arrays.
420 lb = np.array([float(_arr_to_scalar(x)) if x is not None else -np.inf
421 for x in lb])
422 ub = np.array([float(_arr_to_scalar(x)) if x is not None else np.inf
423 for x in ub])
425 return lb, ub
428def strict_bounds(lb, ub, keep_feasible, n_vars):
429 """Remove bounds which are not asked to be kept feasible."""
430 strict_lb = np.resize(lb, n_vars).astype(float)
431 strict_ub = np.resize(ub, n_vars).astype(float)
432 keep_feasible = np.resize(keep_feasible, n_vars)
433 strict_lb[~keep_feasible] = -np.inf
434 strict_ub[~keep_feasible] = np.inf
435 return strict_lb, strict_ub
438def new_constraint_to_old(con, x0):
439 """
440 Converts new-style constraint objects to old-style constraint dictionaries.
441 """
442 if isinstance(con, NonlinearConstraint):
443 if (con.finite_diff_jac_sparsity is not None or
444 con.finite_diff_rel_step is not None or
445 not isinstance(con.hess, BFGS) or # misses user specified BFGS
446 con.keep_feasible):
447 warn("Constraint options `finite_diff_jac_sparsity`, "
448 "`finite_diff_rel_step`, `keep_feasible`, and `hess`"
449 "are ignored by this method.", OptimizeWarning)
451 fun = con.fun
452 if callable(con.jac):
453 jac = con.jac
454 else:
455 jac = None
457 else: # LinearConstraint
458 if np.any(con.keep_feasible):
459 warn("Constraint option `keep_feasible` is ignored by this "
460 "method.", OptimizeWarning)
462 A = con.A
463 if issparse(A):
464 A = A.toarray()
465 fun = lambda x: np.dot(A, x)
466 jac = lambda x: A
468 # FIXME: when bugs in VectorFunction/LinearVectorFunction are worked out,
469 # use pcon.fun.fun and pcon.fun.jac. Until then, get fun/jac above.
470 pcon = PreparedConstraint(con, x0)
471 lb, ub = pcon.bounds
473 i_eq = lb == ub
474 i_bound_below = np.logical_xor(lb != -np.inf, i_eq)
475 i_bound_above = np.logical_xor(ub != np.inf, i_eq)
476 i_unbounded = np.logical_and(lb == -np.inf, ub == np.inf)
478 if np.any(i_unbounded):
479 warn("At least one constraint is unbounded above and below. Such "
480 "constraints are ignored.", OptimizeWarning)
482 ceq = []
483 if np.any(i_eq):
484 def f_eq(x):
485 y = np.array(fun(x)).flatten()
486 return y[i_eq] - lb[i_eq]
487 ceq = [{"type": "eq", "fun": f_eq}]
489 if jac is not None:
490 def j_eq(x):
491 dy = jac(x)
492 if issparse(dy):
493 dy = dy.toarray()
494 dy = np.atleast_2d(dy)
495 return dy[i_eq, :]
496 ceq[0]["jac"] = j_eq
498 cineq = []
499 n_bound_below = np.sum(i_bound_below)
500 n_bound_above = np.sum(i_bound_above)
501 if n_bound_below + n_bound_above:
502 def f_ineq(x):
503 y = np.zeros(n_bound_below + n_bound_above)
504 y_all = np.array(fun(x)).flatten()
505 y[:n_bound_below] = y_all[i_bound_below] - lb[i_bound_below]
506 y[n_bound_below:] = -(y_all[i_bound_above] - ub[i_bound_above])
507 return y
508 cineq = [{"type": "ineq", "fun": f_ineq}]
510 if jac is not None:
511 def j_ineq(x):
512 dy = np.zeros((n_bound_below + n_bound_above, len(x0)))
513 dy_all = jac(x)
514 if issparse(dy_all):
515 dy_all = dy_all.toarray()
516 dy_all = np.atleast_2d(dy_all)
517 dy[:n_bound_below, :] = dy_all[i_bound_below]
518 dy[n_bound_below:, :] = -dy_all[i_bound_above]
519 return dy
520 cineq[0]["jac"] = j_ineq
522 old_constraints = ceq + cineq
524 if len(old_constraints) > 1:
525 warn("Equality and inequality constraints are specified in the same "
526 "element of the constraint list. For efficient use with this "
527 "method, equality and inequality constraints should be specified "
528 "in separate elements of the constraint list. ", OptimizeWarning)
529 return old_constraints
532def old_constraint_to_new(ic, con):
533 """
534 Converts old-style constraint dictionaries to new-style constraint objects.
535 """
536 # check type
537 try:
538 ctype = con['type'].lower()
539 except KeyError as e:
540 raise KeyError('Constraint %d has no type defined.' % ic) from e
541 except TypeError as e:
542 raise TypeError(
543 'Constraints must be a sequence of dictionaries.'
544 ) from e
545 except AttributeError as e:
546 raise TypeError("Constraint's type must be a string.") from e
547 else:
548 if ctype not in ['eq', 'ineq']:
549 raise ValueError("Unknown constraint type '%s'." % con['type'])
550 if 'fun' not in con:
551 raise ValueError('Constraint %d has no function defined.' % ic)
553 lb = 0
554 if ctype == 'eq':
555 ub = 0
556 else:
557 ub = np.inf
559 jac = '2-point'
560 if 'args' in con:
561 args = con['args']
562 fun = lambda x: con['fun'](x, *args)
563 if 'jac' in con:
564 jac = lambda x: con['jac'](x, *args)
565 else:
566 fun = con['fun']
567 if 'jac' in con:
568 jac = con['jac']
570 return NonlinearConstraint(fun, lb, ub, jac)