Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_constraints.py: 12%
217 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
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 : dense 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 : dense 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 if issparse(lb) or issparse(ub):
180 raise ValueError("Constraint limits must be dense arrays.")
181 self.lb = np.atleast_1d(lb).astype(np.float64)
182 self.ub = np.atleast_1d(ub).astype(np.float64)
184 if issparse(keep_feasible):
185 raise ValueError("`keep_feasible` must be a dense array.")
186 self.keep_feasible = np.atleast_1d(keep_feasible).astype(bool)
187 self._input_validation()
189 def residual(self, x):
190 """
191 Calculate the residual between the constraint function and the limits
193 For a linear constraint of the form::
195 lb <= A@x <= ub
197 the lower and upper residuals between ``A@x`` and the limits are values
198 ``sl`` and ``sb`` such that::
200 lb + sl == A@x == ub - sb
202 When all elements of ``sl`` and ``sb`` are positive, all elements of
203 the constraint are satisfied; a negative element in ``sl`` or ``sb``
204 indicates that the corresponding element of the constraint is not
205 satisfied.
207 Parameters
208 ----------
209 x: array_like
210 Vector of independent variables
212 Returns
213 -------
214 sl, sb : array-like
215 The lower and upper residuals
216 """
217 return self.A@x - self.lb, self.ub - self.A@x
220class Bounds:
221 """Bounds constraint on the variables.
223 The constraint has the general inequality form::
225 lb <= x <= ub
227 It is possible to use equal bounds to represent an equality constraint or
228 infinite bounds to represent a one-sided constraint.
230 Parameters
231 ----------
232 lb, ub : dense array_like, optional
233 Lower and upper bounds on independent variables. `lb`, `ub`, and
234 `keep_feasible` must be the same shape or broadcastable.
235 Set components of `lb` and `ub` equal
236 to fix a variable. Use ``np.inf`` with an appropriate sign to disable
237 bounds on all or some variables. Note that you can mix constraints of
238 different types: interval, one-sided or equality, by setting different
239 components of `lb` and `ub` as necessary. Defaults to ``lb = -np.inf``
240 and ``ub = np.inf`` (no bounds).
241 keep_feasible : dense array_like of bool, optional
242 Whether to keep the constraint components feasible throughout
243 iterations. Must be broadcastable with `lb` and `ub`.
244 Default is False. Has no effect for equality constraints.
245 """
246 def _input_validation(self):
247 try:
248 res = np.broadcast_arrays(self.lb, self.ub, self.keep_feasible)
249 self.lb, self.ub, self.keep_feasible = res
250 except ValueError:
251 message = "`lb`, `ub`, and `keep_feasible` must be broadcastable."
252 raise ValueError(message)
254 def __init__(self, lb=-np.inf, ub=np.inf, keep_feasible=False):
255 if issparse(lb) or issparse(ub):
256 raise ValueError("Lower and upper bounds must be dense arrays.")
257 self.lb = np.atleast_1d(lb)
258 self.ub = np.atleast_1d(ub)
260 if issparse(keep_feasible):
261 raise ValueError("`keep_feasible` must be a dense array.")
262 self.keep_feasible = np.atleast_1d(keep_feasible).astype(bool)
263 self._input_validation()
265 def __repr__(self):
266 start = f"{type(self).__name__}({self.lb!r}, {self.ub!r}"
267 if np.any(self.keep_feasible):
268 end = f", keep_feasible={self.keep_feasible!r})"
269 else:
270 end = ")"
271 return start + end
273 def residual(self, x):
274 """Calculate the residual (slack) between the input and the bounds
276 For a bound constraint of the form::
278 lb <= x <= ub
280 the lower and upper residuals between `x` and the bounds are values
281 ``sl`` and ``sb`` such that::
283 lb + sl == x == ub - sb
285 When all elements of ``sl`` and ``sb`` are positive, all elements of
286 ``x`` lie within the bounds; a negative element in ``sl`` or ``sb``
287 indicates that the corresponding element of ``x`` is out of bounds.
289 Parameters
290 ----------
291 x: array_like
292 Vector of independent variables
294 Returns
295 -------
296 sl, sb : array-like
297 The lower and upper residuals
298 """
299 return x - self.lb, self.ub - x
302class PreparedConstraint:
303 """Constraint prepared from a user defined constraint.
305 On creation it will check whether a constraint definition is valid and
306 the initial point is feasible. If created successfully, it will contain
307 the attributes listed below.
309 Parameters
310 ----------
311 constraint : {NonlinearConstraint, LinearConstraint`, Bounds}
312 Constraint to check and prepare.
313 x0 : array_like
314 Initial vector of independent variables.
315 sparse_jacobian : bool or None, optional
316 If bool, then the Jacobian of the constraint will be converted
317 to the corresponded format if necessary. If None (default), such
318 conversion is not made.
319 finite_diff_bounds : 2-tuple, optional
320 Lower and upper bounds on the independent variables for the finite
321 difference approximation, if applicable. Defaults to no bounds.
323 Attributes
324 ----------
325 fun : {VectorFunction, LinearVectorFunction, IdentityVectorFunction}
326 Function defining the constraint wrapped by one of the convenience
327 classes.
328 bounds : 2-tuple
329 Contains lower and upper bounds for the constraints --- lb and ub.
330 These are converted to ndarray and have a size equal to the number of
331 the constraints.
332 keep_feasible : ndarray
333 Array indicating which components must be kept feasible with a size
334 equal to the number of the constraints.
335 """
336 def __init__(self, constraint, x0, sparse_jacobian=None,
337 finite_diff_bounds=(-np.inf, np.inf)):
338 if isinstance(constraint, NonlinearConstraint):
339 fun = VectorFunction(constraint.fun, x0,
340 constraint.jac, constraint.hess,
341 constraint.finite_diff_rel_step,
342 constraint.finite_diff_jac_sparsity,
343 finite_diff_bounds, sparse_jacobian)
344 elif isinstance(constraint, LinearConstraint):
345 fun = LinearVectorFunction(constraint.A, x0, sparse_jacobian)
346 elif isinstance(constraint, Bounds):
347 fun = IdentityVectorFunction(x0, sparse_jacobian)
348 else:
349 raise ValueError("`constraint` of an unknown type is passed.")
351 m = fun.m
353 lb = np.asarray(constraint.lb, dtype=float)
354 ub = np.asarray(constraint.ub, dtype=float)
355 keep_feasible = np.asarray(constraint.keep_feasible, dtype=bool)
357 lb = np.broadcast_to(lb, m)
358 ub = np.broadcast_to(ub, m)
359 keep_feasible = np.broadcast_to(keep_feasible, m)
361 if keep_feasible.shape != (m,):
362 raise ValueError("`keep_feasible` has a wrong shape.")
364 mask = keep_feasible & (lb != ub)
365 f0 = fun.f
366 if np.any(f0[mask] < lb[mask]) or np.any(f0[mask] > ub[mask]):
367 raise ValueError("`x0` is infeasible with respect to some "
368 "inequality constraint with `keep_feasible` "
369 "set to True.")
371 self.fun = fun
372 self.bounds = (lb, ub)
373 self.keep_feasible = keep_feasible
375 def violation(self, x):
376 """How much the constraint is exceeded by.
378 Parameters
379 ----------
380 x : array-like
381 Vector of independent variables
383 Returns
384 -------
385 excess : array-like
386 How much the constraint is exceeded by, for each of the
387 constraints specified by `PreparedConstraint.fun`.
388 """
389 with suppress_warnings() as sup:
390 sup.filter(UserWarning)
391 ev = self.fun.fun(np.asarray(x))
393 excess_lb = np.maximum(self.bounds[0] - ev, 0)
394 excess_ub = np.maximum(ev - self.bounds[1], 0)
396 return excess_lb + excess_ub
399def new_bounds_to_old(lb, ub, n):
400 """Convert the new bounds representation to the old one.
402 The new representation is a tuple (lb, ub) and the old one is a list
403 containing n tuples, ith containing lower and upper bound on a ith
404 variable.
405 If any of the entries in lb/ub are -np.inf/np.inf they are replaced by
406 None.
407 """
408 lb = np.broadcast_to(lb, n)
409 ub = np.broadcast_to(ub, n)
411 lb = [float(x) if x > -np.inf else None for x in lb]
412 ub = [float(x) if x < np.inf else None for x in ub]
414 return list(zip(lb, ub))
417def old_bound_to_new(bounds):
418 """Convert the old bounds representation to the new one.
420 The new representation is a tuple (lb, ub) and the old one is a list
421 containing n tuples, ith containing lower and upper bound on a ith
422 variable.
423 If any of the entries in lb/ub are None they are replaced by
424 -np.inf/np.inf.
425 """
426 lb, ub = zip(*bounds)
428 # Convert occurrences of None to -inf or inf, and replace occurrences of
429 # any numpy array x with x.item(). Then wrap the results in numpy arrays.
430 lb = np.array([float(_arr_to_scalar(x)) if x is not None else -np.inf
431 for x in lb])
432 ub = np.array([float(_arr_to_scalar(x)) if x is not None else np.inf
433 for x in ub])
435 return lb, ub
438def strict_bounds(lb, ub, keep_feasible, n_vars):
439 """Remove bounds which are not asked to be kept feasible."""
440 strict_lb = np.resize(lb, n_vars).astype(float)
441 strict_ub = np.resize(ub, n_vars).astype(float)
442 keep_feasible = np.resize(keep_feasible, n_vars)
443 strict_lb[~keep_feasible] = -np.inf
444 strict_ub[~keep_feasible] = np.inf
445 return strict_lb, strict_ub
448def new_constraint_to_old(con, x0):
449 """
450 Converts new-style constraint objects to old-style constraint dictionaries.
451 """
452 if isinstance(con, NonlinearConstraint):
453 if (con.finite_diff_jac_sparsity is not None or
454 con.finite_diff_rel_step is not None or
455 not isinstance(con.hess, BFGS) or # misses user specified BFGS
456 con.keep_feasible):
457 warn("Constraint options `finite_diff_jac_sparsity`, "
458 "`finite_diff_rel_step`, `keep_feasible`, and `hess`"
459 "are ignored by this method.", OptimizeWarning)
461 fun = con.fun
462 if callable(con.jac):
463 jac = con.jac
464 else:
465 jac = None
467 else: # LinearConstraint
468 if np.any(con.keep_feasible):
469 warn("Constraint option `keep_feasible` is ignored by this "
470 "method.", OptimizeWarning)
472 A = con.A
473 if issparse(A):
474 A = A.toarray()
475 def fun(x):
476 return np.dot(A, x)
477 def jac(x):
478 return A
480 # FIXME: when bugs in VectorFunction/LinearVectorFunction are worked out,
481 # use pcon.fun.fun and pcon.fun.jac. Until then, get fun/jac above.
482 pcon = PreparedConstraint(con, x0)
483 lb, ub = pcon.bounds
485 i_eq = lb == ub
486 i_bound_below = np.logical_xor(lb != -np.inf, i_eq)
487 i_bound_above = np.logical_xor(ub != np.inf, i_eq)
488 i_unbounded = np.logical_and(lb == -np.inf, ub == np.inf)
490 if np.any(i_unbounded):
491 warn("At least one constraint is unbounded above and below. Such "
492 "constraints are ignored.", OptimizeWarning)
494 ceq = []
495 if np.any(i_eq):
496 def f_eq(x):
497 y = np.array(fun(x)).flatten()
498 return y[i_eq] - lb[i_eq]
499 ceq = [{"type": "eq", "fun": f_eq}]
501 if jac is not None:
502 def j_eq(x):
503 dy = jac(x)
504 if issparse(dy):
505 dy = dy.toarray()
506 dy = np.atleast_2d(dy)
507 return dy[i_eq, :]
508 ceq[0]["jac"] = j_eq
510 cineq = []
511 n_bound_below = np.sum(i_bound_below)
512 n_bound_above = np.sum(i_bound_above)
513 if n_bound_below + n_bound_above:
514 def f_ineq(x):
515 y = np.zeros(n_bound_below + n_bound_above)
516 y_all = np.array(fun(x)).flatten()
517 y[:n_bound_below] = y_all[i_bound_below] - lb[i_bound_below]
518 y[n_bound_below:] = -(y_all[i_bound_above] - ub[i_bound_above])
519 return y
520 cineq = [{"type": "ineq", "fun": f_ineq}]
522 if jac is not None:
523 def j_ineq(x):
524 dy = np.zeros((n_bound_below + n_bound_above, len(x0)))
525 dy_all = jac(x)
526 if issparse(dy_all):
527 dy_all = dy_all.toarray()
528 dy_all = np.atleast_2d(dy_all)
529 dy[:n_bound_below, :] = dy_all[i_bound_below]
530 dy[n_bound_below:, :] = -dy_all[i_bound_above]
531 return dy
532 cineq[0]["jac"] = j_ineq
534 old_constraints = ceq + cineq
536 if len(old_constraints) > 1:
537 warn("Equality and inequality constraints are specified in the same "
538 "element of the constraint list. For efficient use with this "
539 "method, equality and inequality constraints should be specified "
540 "in separate elements of the constraint list. ", OptimizeWarning)
541 return old_constraints
544def old_constraint_to_new(ic, con):
545 """
546 Converts old-style constraint dictionaries to new-style constraint objects.
547 """
548 # check type
549 try:
550 ctype = con['type'].lower()
551 except KeyError as e:
552 raise KeyError('Constraint %d has no type defined.' % ic) from e
553 except TypeError as e:
554 raise TypeError(
555 'Constraints must be a sequence of dictionaries.'
556 ) from e
557 except AttributeError as e:
558 raise TypeError("Constraint's type must be a string.") from e
559 else:
560 if ctype not in ['eq', 'ineq']:
561 raise ValueError("Unknown constraint type '%s'." % con['type'])
562 if 'fun' not in con:
563 raise ValueError('Constraint %d has no function defined.' % ic)
565 lb = 0
566 if ctype == 'eq':
567 ub = 0
568 else:
569 ub = np.inf
571 jac = '2-point'
572 if 'args' in con:
573 args = con['args']
574 def fun(x):
575 return con["fun"](x, *args)
576 if 'jac' in con:
577 def jac(x):
578 return con["jac"](x, *args)
579 else:
580 fun = con['fun']
581 if 'jac' in con:
582 jac = con['jac']
584 return NonlinearConstraint(fun, lb, ub, jac)