Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_linprog_util.py: 5%
487 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""
2Method agnostic utility functions for linear progamming
3"""
5import numpy as np
6import scipy.sparse as sps
7from warnings import warn
8from ._optimize import OptimizeWarning
9from scipy.optimize._remove_redundancy import (
10 _remove_redundancy_svd, _remove_redundancy_pivot_sparse,
11 _remove_redundancy_pivot_dense, _remove_redundancy_id
12 )
13from collections import namedtuple
15_LPProblem = namedtuple('_LPProblem',
16 'c A_ub b_ub A_eq b_eq bounds x0 integrality')
17_LPProblem.__new__.__defaults__ = (None,) * 7 # make c the only required arg
18_LPProblem.__doc__ = \
19 """ Represents a linear-programming problem.
21 Attributes
22 ----------
23 c : 1D array
24 The coefficients of the linear objective function to be minimized.
25 A_ub : 2D array, optional
26 The inequality constraint matrix. Each row of ``A_ub`` specifies the
27 coefficients of a linear inequality constraint on ``x``.
28 b_ub : 1D array, optional
29 The inequality constraint vector. Each element represents an
30 upper bound on the corresponding value of ``A_ub @ x``.
31 A_eq : 2D array, optional
32 The equality constraint matrix. Each row of ``A_eq`` specifies the
33 coefficients of a linear equality constraint on ``x``.
34 b_eq : 1D array, optional
35 The equality constraint vector. Each element of ``A_eq @ x`` must equal
36 the corresponding element of ``b_eq``.
37 bounds : various valid formats, optional
38 The bounds of ``x``, as ``min`` and ``max`` pairs.
39 If bounds are specified for all N variables separately, valid formats
40 are:
41 * a 2D array (N x 2);
42 * a sequence of N sequences, each with 2 values.
43 If all variables have the same bounds, the bounds can be specified as
44 a 1-D or 2-D array or sequence with 2 scalar values.
45 If all variables have a lower bound of 0 and no upper bound, the bounds
46 parameter can be omitted (or given as None).
47 Absent lower and/or upper bounds can be specified as -numpy.inf (no
48 lower bound), numpy.inf (no upper bound) or None (both).
49 x0 : 1D array, optional
50 Guess values of the decision variables, which will be refined by
51 the optimization algorithm. This argument is currently used only by the
52 'revised simplex' method, and can only be used if `x0` represents a
53 basic feasible solution.
54 integrality : 1-D array or int, optional
55 Indicates the type of integrality constraint on each decision variable.
57 ``0`` : Continuous variable; no integrality constraint.
59 ``1`` : Integer variable; decision variable must be an integer
60 within `bounds`.
62 ``2`` : Semi-continuous variable; decision variable must be within
63 `bounds` or take value ``0``.
65 ``3`` : Semi-integer variable; decision variable must be an integer
66 within `bounds` or take value ``0``.
68 By default, all variables are continuous.
70 For mixed integrality constraints, supply an array of shape `c.shape`.
71 To infer a constraint on each decision variable from shorter inputs,
72 the argument will be broadcasted to `c.shape` using `np.broadcast_to`.
74 This argument is currently used only by the ``'highs'`` method and
75 ignored otherwise.
77 Notes
78 -----
79 This namedtuple supports 2 ways of initialization:
80 >>> lp1 = _LPProblem(c=[-1, 4], A_ub=[[-3, 1], [1, 2]], b_ub=[6, 4])
81 >>> lp2 = _LPProblem([-1, 4], [[-3, 1], [1, 2]], [6, 4])
83 Note that only ``c`` is a required argument here, whereas all other arguments
84 ``A_ub``, ``b_ub``, ``A_eq``, ``b_eq``, ``bounds``, ``x0`` are optional with
85 default values of None.
86 For example, ``A_eq`` and ``b_eq`` can be set without ``A_ub`` or ``b_ub``:
87 >>> lp3 = _LPProblem(c=[-1, 4], A_eq=[[2, 1]], b_eq=[10])
88 """
91def _check_sparse_inputs(options, meth, A_ub, A_eq):
92 """
93 Check the provided ``A_ub`` and ``A_eq`` matrices conform to the specified
94 optional sparsity variables.
96 Parameters
97 ----------
98 A_ub : 2-D array, optional
99 2-D array such that ``A_ub @ x`` gives the values of the upper-bound
100 inequality constraints at ``x``.
101 A_eq : 2-D array, optional
102 2-D array such that ``A_eq @ x`` gives the values of the equality
103 constraints at ``x``.
104 options : dict
105 A dictionary of solver options. All methods accept the following
106 generic options:
108 maxiter : int
109 Maximum number of iterations to perform.
110 disp : bool
111 Set to True to print convergence messages.
113 For method-specific options, see :func:`show_options('linprog')`.
114 method : str, optional
115 The algorithm used to solve the standard form problem.
117 Returns
118 -------
119 A_ub : 2-D array, optional
120 2-D array such that ``A_ub @ x`` gives the values of the upper-bound
121 inequality constraints at ``x``.
122 A_eq : 2-D array, optional
123 2-D array such that ``A_eq @ x`` gives the values of the equality
124 constraints at ``x``.
125 options : dict
126 A dictionary of solver options. All methods accept the following
127 generic options:
129 maxiter : int
130 Maximum number of iterations to perform.
131 disp : bool
132 Set to True to print convergence messages.
134 For method-specific options, see :func:`show_options('linprog')`.
135 """
136 # This is an undocumented option for unit testing sparse presolve
137 _sparse_presolve = options.pop('_sparse_presolve', False)
138 if _sparse_presolve and A_eq is not None:
139 A_eq = sps.coo_matrix(A_eq)
140 if _sparse_presolve and A_ub is not None:
141 A_ub = sps.coo_matrix(A_ub)
143 sparse_constraint = sps.issparse(A_eq) or sps.issparse(A_ub)
145 preferred_methods = {"highs", "highs-ds", "highs-ipm"}
146 dense_methods = {"simplex", "revised simplex"}
147 if meth in dense_methods and sparse_constraint:
148 raise ValueError(f"Method '{meth}' does not support sparse "
149 "constraint matrices. Please consider using one of "
150 f"{preferred_methods}.")
152 sparse = options.get('sparse', False)
153 if not sparse and sparse_constraint and meth == 'interior-point':
154 options['sparse'] = True
155 warn("Sparse constraint matrix detected; setting 'sparse':True.",
156 OptimizeWarning, stacklevel=4)
157 return options, A_ub, A_eq
160def _format_A_constraints(A, n_x, sparse_lhs=False):
161 """Format the left hand side of the constraints to a 2-D array
163 Parameters
164 ----------
165 A : 2-D array
166 2-D array such that ``A @ x`` gives the values of the upper-bound
167 (in)equality constraints at ``x``.
168 n_x : int
169 The number of variables in the linear programming problem.
170 sparse_lhs : bool
171 Whether either of `A_ub` or `A_eq` are sparse. If true return a
172 coo_matrix instead of a numpy array.
174 Returns
175 -------
176 np.ndarray or sparse.coo_matrix
177 2-D array such that ``A @ x`` gives the values of the upper-bound
178 (in)equality constraints at ``x``.
180 """
181 if sparse_lhs:
182 return sps.coo_matrix(
183 (0, n_x) if A is None else A, dtype=float, copy=True
184 )
185 elif A is None:
186 return np.zeros((0, n_x), dtype=float)
187 else:
188 return np.array(A, dtype=float, copy=True)
191def _format_b_constraints(b):
192 """Format the upper bounds of the constraints to a 1-D array
194 Parameters
195 ----------
196 b : 1-D array
197 1-D array of values representing the upper-bound of each (in)equality
198 constraint (row) in ``A``.
200 Returns
201 -------
202 1-D np.array
203 1-D array of values representing the upper-bound of each (in)equality
204 constraint (row) in ``A``.
206 """
207 if b is None:
208 return np.array([], dtype=float)
209 b = np.array(b, dtype=float, copy=True).squeeze()
210 return b if b.size != 1 else b.reshape((-1))
213def _clean_inputs(lp):
214 """
215 Given user inputs for a linear programming problem, return the
216 objective vector, upper bound constraints, equality constraints,
217 and simple bounds in a preferred format.
219 Parameters
220 ----------
221 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields:
223 c : 1D array
224 The coefficients of the linear objective function to be minimized.
225 A_ub : 2D array, optional
226 The inequality constraint matrix. Each row of ``A_ub`` specifies the
227 coefficients of a linear inequality constraint on ``x``.
228 b_ub : 1D array, optional
229 The inequality constraint vector. Each element represents an
230 upper bound on the corresponding value of ``A_ub @ x``.
231 A_eq : 2D array, optional
232 The equality constraint matrix. Each row of ``A_eq`` specifies the
233 coefficients of a linear equality constraint on ``x``.
234 b_eq : 1D array, optional
235 The equality constraint vector. Each element of ``A_eq @ x`` must equal
236 the corresponding element of ``b_eq``.
237 bounds : various valid formats, optional
238 The bounds of ``x``, as ``min`` and ``max`` pairs.
239 If bounds are specified for all N variables separately, valid formats are:
240 * a 2D array (2 x N or N x 2);
241 * a sequence of N sequences, each with 2 values.
242 If all variables have the same bounds, a single pair of values can
243 be specified. Valid formats are:
244 * a sequence with 2 scalar values;
245 * a sequence with a single element containing 2 scalar values.
246 If all variables have a lower bound of 0 and no upper bound, the bounds
247 parameter can be omitted (or given as None).
248 x0 : 1D array, optional
249 Guess values of the decision variables, which will be refined by
250 the optimization algorithm. This argument is currently used only by the
251 'revised simplex' method, and can only be used if `x0` represents a
252 basic feasible solution.
254 Returns
255 -------
256 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields:
258 c : 1D array
259 The coefficients of the linear objective function to be minimized.
260 A_ub : 2D array, optional
261 The inequality constraint matrix. Each row of ``A_ub`` specifies the
262 coefficients of a linear inequality constraint on ``x``.
263 b_ub : 1D array, optional
264 The inequality constraint vector. Each element represents an
265 upper bound on the corresponding value of ``A_ub @ x``.
266 A_eq : 2D array, optional
267 The equality constraint matrix. Each row of ``A_eq`` specifies the
268 coefficients of a linear equality constraint on ``x``.
269 b_eq : 1D array, optional
270 The equality constraint vector. Each element of ``A_eq @ x`` must equal
271 the corresponding element of ``b_eq``.
272 bounds : 2D array
273 The bounds of ``x``, as ``min`` and ``max`` pairs, one for each of the N
274 elements of ``x``. The N x 2 array contains lower bounds in the first
275 column and upper bounds in the 2nd. Unbounded variables have lower
276 bound -np.inf and/or upper bound np.inf.
277 x0 : 1D array, optional
278 Guess values of the decision variables, which will be refined by
279 the optimization algorithm. This argument is currently used only by the
280 'revised simplex' method, and can only be used if `x0` represents a
281 basic feasible solution.
283 """
284 c, A_ub, b_ub, A_eq, b_eq, bounds, x0, integrality = lp
286 if c is None:
287 raise TypeError
289 try:
290 c = np.array(c, dtype=np.float64, copy=True).squeeze()
291 except ValueError as e:
292 raise TypeError(
293 "Invalid input for linprog: c must be a 1-D array of numerical "
294 "coefficients") from e
295 else:
296 # If c is a single value, convert it to a 1-D array.
297 if c.size == 1:
298 c = c.reshape((-1))
300 n_x = len(c)
301 if n_x == 0 or len(c.shape) != 1:
302 raise ValueError(
303 "Invalid input for linprog: c must be a 1-D array and must "
304 "not have more than one non-singleton dimension")
305 if not np.isfinite(c).all():
306 raise ValueError(
307 "Invalid input for linprog: c must not contain values "
308 "inf, nan, or None")
310 sparse_lhs = sps.issparse(A_eq) or sps.issparse(A_ub)
311 try:
312 A_ub = _format_A_constraints(A_ub, n_x, sparse_lhs=sparse_lhs)
313 except ValueError as e:
314 raise TypeError(
315 "Invalid input for linprog: A_ub must be a 2-D array "
316 "of numerical values") from e
317 else:
318 n_ub = A_ub.shape[0]
319 if len(A_ub.shape) != 2 or A_ub.shape[1] != n_x:
320 raise ValueError(
321 "Invalid input for linprog: A_ub must have exactly two "
322 "dimensions, and the number of columns in A_ub must be "
323 "equal to the size of c")
324 if (sps.issparse(A_ub) and not np.isfinite(A_ub.data).all()
325 or not sps.issparse(A_ub) and not np.isfinite(A_ub).all()):
326 raise ValueError(
327 "Invalid input for linprog: A_ub must not contain values "
328 "inf, nan, or None")
330 try:
331 b_ub = _format_b_constraints(b_ub)
332 except ValueError as e:
333 raise TypeError(
334 "Invalid input for linprog: b_ub must be a 1-D array of "
335 "numerical values, each representing the upper bound of an "
336 "inequality constraint (row) in A_ub") from e
337 else:
338 if b_ub.shape != (n_ub,):
339 raise ValueError(
340 "Invalid input for linprog: b_ub must be a 1-D array; b_ub "
341 "must not have more than one non-singleton dimension and "
342 "the number of rows in A_ub must equal the number of values "
343 "in b_ub")
344 if not np.isfinite(b_ub).all():
345 raise ValueError(
346 "Invalid input for linprog: b_ub must not contain values "
347 "inf, nan, or None")
349 try:
350 A_eq = _format_A_constraints(A_eq, n_x, sparse_lhs=sparse_lhs)
351 except ValueError as e:
352 raise TypeError(
353 "Invalid input for linprog: A_eq must be a 2-D array "
354 "of numerical values") from e
355 else:
356 n_eq = A_eq.shape[0]
357 if len(A_eq.shape) != 2 or A_eq.shape[1] != n_x:
358 raise ValueError(
359 "Invalid input for linprog: A_eq must have exactly two "
360 "dimensions, and the number of columns in A_eq must be "
361 "equal to the size of c")
363 if (sps.issparse(A_eq) and not np.isfinite(A_eq.data).all()
364 or not sps.issparse(A_eq) and not np.isfinite(A_eq).all()):
365 raise ValueError(
366 "Invalid input for linprog: A_eq must not contain values "
367 "inf, nan, or None")
369 try:
370 b_eq = _format_b_constraints(b_eq)
371 except ValueError as e:
372 raise TypeError(
373 "Invalid input for linprog: b_eq must be a dense, 1-D array of "
374 "numerical values, each representing the right hand side of an "
375 "equality constraint (row) in A_eq") from e
376 else:
377 if b_eq.shape != (n_eq,):
378 raise ValueError(
379 "Invalid input for linprog: b_eq must be a 1-D array; b_eq "
380 "must not have more than one non-singleton dimension and "
381 "the number of rows in A_eq must equal the number of values "
382 "in b_eq")
383 if not np.isfinite(b_eq).all():
384 raise ValueError(
385 "Invalid input for linprog: b_eq must not contain values "
386 "inf, nan, or None")
388 # x0 gives a (optional) starting solution to the solver. If x0 is None,
389 # skip the checks. Initial solution will be generated automatically.
390 if x0 is not None:
391 try:
392 x0 = np.array(x0, dtype=float, copy=True).squeeze()
393 except ValueError as e:
394 raise TypeError(
395 "Invalid input for linprog: x0 must be a 1-D array of "
396 "numerical coefficients") from e
397 if x0.ndim == 0:
398 x0 = x0.reshape((-1))
399 if len(x0) == 0 or x0.ndim != 1:
400 raise ValueError(
401 "Invalid input for linprog: x0 should be a 1-D array; it "
402 "must not have more than one non-singleton dimension")
403 if not x0.size == c.size:
404 raise ValueError(
405 "Invalid input for linprog: x0 and c should contain the "
406 "same number of elements")
407 if not np.isfinite(x0).all():
408 raise ValueError(
409 "Invalid input for linprog: x0 must not contain values "
410 "inf, nan, or None")
412 # Bounds can be one of these formats:
413 # (1) a 2-D array or sequence, with shape N x 2
414 # (2) a 1-D or 2-D sequence or array with 2 scalars
415 # (3) None (or an empty sequence or array)
416 # Unspecified bounds can be represented by None or (-)np.inf.
417 # All formats are converted into a N x 2 np.array with (-)np.inf where
418 # bounds are unspecified.
420 # Prepare clean bounds array
421 bounds_clean = np.zeros((n_x, 2), dtype=float)
423 # Convert to a numpy array.
424 # np.array(..,dtype=float) raises an error if dimensions are inconsistent
425 # or if there are invalid data types in bounds. Just add a linprog prefix
426 # to the error and re-raise.
427 # Creating at least a 2-D array simplifies the cases to distinguish below.
428 if bounds is None or np.array_equal(bounds, []) or np.array_equal(bounds, [[]]):
429 bounds = (0, np.inf)
430 try:
431 bounds_conv = np.atleast_2d(np.array(bounds, dtype=float))
432 except ValueError as e:
433 raise ValueError(
434 "Invalid input for linprog: unable to interpret bounds, "
435 "check values and dimensions: " + e.args[0]) from e
436 except TypeError as e:
437 raise TypeError(
438 "Invalid input for linprog: unable to interpret bounds, "
439 "check values and dimensions: " + e.args[0]) from e
441 # Check bounds options
442 bsh = bounds_conv.shape
443 if len(bsh) > 2:
444 # Do not try to handle multidimensional bounds input
445 raise ValueError(
446 "Invalid input for linprog: provide a 2-D array for bounds, "
447 "not a {:d}-D array.".format(len(bsh)))
448 elif np.all(bsh == (n_x, 2)):
449 # Regular N x 2 array
450 bounds_clean = bounds_conv
451 elif (np.all(bsh == (2, 1)) or np.all(bsh == (1, 2))):
452 # 2 values: interpret as overall lower and upper bound
453 bounds_flat = bounds_conv.flatten()
454 bounds_clean[:, 0] = bounds_flat[0]
455 bounds_clean[:, 1] = bounds_flat[1]
456 elif np.all(bsh == (2, n_x)):
457 # Reject a 2 x N array
458 raise ValueError(
459 "Invalid input for linprog: provide a {:d} x 2 array for bounds, "
460 "not a 2 x {:d} array.".format(n_x, n_x))
461 else:
462 raise ValueError(
463 "Invalid input for linprog: unable to interpret bounds with this "
464 "dimension tuple: {0}.".format(bsh))
466 # The process above creates nan-s where the input specified None
467 # Convert the nan-s in the 1st column to -np.inf and in the 2nd column
468 # to np.inf
469 i_none = np.isnan(bounds_clean[:, 0])
470 bounds_clean[i_none, 0] = -np.inf
471 i_none = np.isnan(bounds_clean[:, 1])
472 bounds_clean[i_none, 1] = np.inf
474 return _LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds_clean, x0, integrality)
477def _presolve(lp, rr, rr_method, tol=1e-9):
478 """
479 Given inputs for a linear programming problem in preferred format,
480 presolve the problem: identify trivial infeasibilities, redundancies,
481 and unboundedness, tighten bounds where possible, and eliminate fixed
482 variables.
484 Parameters
485 ----------
486 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields:
488 c : 1D array
489 The coefficients of the linear objective function to be minimized.
490 A_ub : 2D array, optional
491 The inequality constraint matrix. Each row of ``A_ub`` specifies the
492 coefficients of a linear inequality constraint on ``x``.
493 b_ub : 1D array, optional
494 The inequality constraint vector. Each element represents an
495 upper bound on the corresponding value of ``A_ub @ x``.
496 A_eq : 2D array, optional
497 The equality constraint matrix. Each row of ``A_eq`` specifies the
498 coefficients of a linear equality constraint on ``x``.
499 b_eq : 1D array, optional
500 The equality constraint vector. Each element of ``A_eq @ x`` must equal
501 the corresponding element of ``b_eq``.
502 bounds : 2D array
503 The bounds of ``x``, as ``min`` and ``max`` pairs, one for each of the N
504 elements of ``x``. The N x 2 array contains lower bounds in the first
505 column and upper bounds in the 2nd. Unbounded variables have lower
506 bound -np.inf and/or upper bound np.inf.
507 x0 : 1D array, optional
508 Guess values of the decision variables, which will be refined by
509 the optimization algorithm. This argument is currently used only by the
510 'revised simplex' method, and can only be used if `x0` represents a
511 basic feasible solution.
513 rr : bool
514 If ``True`` attempts to eliminate any redundant rows in ``A_eq``.
515 Set False if ``A_eq`` is known to be of full row rank, or if you are
516 looking for a potential speedup (at the expense of reliability).
517 rr_method : string
518 Method used to identify and remove redundant rows from the
519 equality constraint matrix after presolve.
520 tol : float
521 The tolerance which determines when a solution is "close enough" to
522 zero in Phase 1 to be considered a basic feasible solution or close
523 enough to positive to serve as an optimal solution.
525 Returns
526 -------
527 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields:
529 c : 1D array
530 The coefficients of the linear objective function to be minimized.
531 A_ub : 2D array, optional
532 The inequality constraint matrix. Each row of ``A_ub`` specifies the
533 coefficients of a linear inequality constraint on ``x``.
534 b_ub : 1D array, optional
535 The inequality constraint vector. Each element represents an
536 upper bound on the corresponding value of ``A_ub @ x``.
537 A_eq : 2D array, optional
538 The equality constraint matrix. Each row of ``A_eq`` specifies the
539 coefficients of a linear equality constraint on ``x``.
540 b_eq : 1D array, optional
541 The equality constraint vector. Each element of ``A_eq @ x`` must equal
542 the corresponding element of ``b_eq``.
543 bounds : 2D array
544 The bounds of ``x``, as ``min`` and ``max`` pairs, possibly tightened.
545 x0 : 1D array, optional
546 Guess values of the decision variables, which will be refined by
547 the optimization algorithm. This argument is currently used only by the
548 'revised simplex' method, and can only be used if `x0` represents a
549 basic feasible solution.
551 c0 : 1D array
552 Constant term in objective function due to fixed (and eliminated)
553 variables.
554 x : 1D array
555 Solution vector (when the solution is trivial and can be determined
556 in presolve)
557 revstack: list of functions
558 the functions in the list reverse the operations of _presolve()
559 the function signature is x_org = f(x_mod), where x_mod is the result
560 of a presolve step and x_org the value at the start of the step
561 (currently, the revstack contains only one function)
562 complete: bool
563 Whether the solution is complete (solved or determined to be infeasible
564 or unbounded in presolve)
565 status : int
566 An integer representing the exit status of the optimization::
568 0 : Optimization terminated successfully
569 1 : Iteration limit reached
570 2 : Problem appears to be infeasible
571 3 : Problem appears to be unbounded
572 4 : Serious numerical difficulties encountered
574 message : str
575 A string descriptor of the exit status of the optimization.
577 References
578 ----------
579 .. [5] Andersen, Erling D. "Finding all linearly dependent rows in
580 large-scale linear programming." Optimization Methods and Software
581 6.3 (1995): 219-227.
582 .. [8] Andersen, Erling D., and Knud D. Andersen. "Presolving in linear
583 programming." Mathematical Programming 71.2 (1995): 221-245.
585 """
586 # ideas from Reference [5] by Andersen and Andersen
587 # however, unlike the reference, this is performed before converting
588 # problem to standard form
589 # There are a few advantages:
590 # * artificial variables have not been added, so matrices are smaller
591 # * bounds have not been converted to constraints yet. (It is better to
592 # do that after presolve because presolve may adjust the simple bounds.)
593 # There are many improvements that can be made, namely:
594 # * implement remaining checks from [5]
595 # * loop presolve until no additional changes are made
596 # * implement additional efficiency improvements in redundancy removal [2]
598 c, A_ub, b_ub, A_eq, b_eq, bounds, x0, _ = lp
600 revstack = [] # record of variables eliminated from problem
601 # constant term in cost function may be added if variables are eliminated
602 c0 = 0
603 complete = False # complete is True if detected infeasible/unbounded
604 x = np.zeros(c.shape) # this is solution vector if completed in presolve
606 status = 0 # all OK unless determined otherwise
607 message = ""
609 # Lower and upper bounds. Copy to prevent feedback.
610 lb = bounds[:, 0].copy()
611 ub = bounds[:, 1].copy()
613 m_eq, n = A_eq.shape
614 m_ub, n = A_ub.shape
616 if (rr_method is not None
617 and rr_method.lower() not in {"svd", "pivot", "id"}):
618 message = ("'" + str(rr_method) + "' is not a valid option "
619 "for redundancy removal. Valid options are 'SVD', "
620 "'pivot', and 'ID'.")
621 raise ValueError(message)
623 if sps.issparse(A_eq):
624 A_eq = A_eq.tocsr()
625 A_ub = A_ub.tocsr()
627 def where(A):
628 return A.nonzero()
630 vstack = sps.vstack
631 else:
632 where = np.where
633 vstack = np.vstack
635 # upper bounds > lower bounds
636 if np.any(ub < lb) or np.any(lb == np.inf) or np.any(ub == -np.inf):
637 status = 2
638 message = ("The problem is (trivially) infeasible since one "
639 "or more upper bounds are smaller than the corresponding "
640 "lower bounds, a lower bound is np.inf or an upper bound "
641 "is -np.inf.")
642 complete = True
643 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0),
644 c0, x, revstack, complete, status, message)
646 # zero row in equality constraints
647 zero_row = np.array(np.sum(A_eq != 0, axis=1) == 0).flatten()
648 if np.any(zero_row):
649 if np.any(
650 np.logical_and(
651 zero_row,
652 np.abs(b_eq) > tol)): # test_zero_row_1
653 # infeasible if RHS is not zero
654 status = 2
655 message = ("The problem is (trivially) infeasible due to a row "
656 "of zeros in the equality constraint matrix with a "
657 "nonzero corresponding constraint value.")
658 complete = True
659 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0),
660 c0, x, revstack, complete, status, message)
661 else: # test_zero_row_2
662 # if RHS is zero, we can eliminate this equation entirely
663 A_eq = A_eq[np.logical_not(zero_row), :]
664 b_eq = b_eq[np.logical_not(zero_row)]
666 # zero row in inequality constraints
667 zero_row = np.array(np.sum(A_ub != 0, axis=1) == 0).flatten()
668 if np.any(zero_row):
669 if np.any(np.logical_and(zero_row, b_ub < -tol)): # test_zero_row_1
670 # infeasible if RHS is less than zero (because LHS is zero)
671 status = 2
672 message = ("The problem is (trivially) infeasible due to a row "
673 "of zeros in the equality constraint matrix with a "
674 "nonzero corresponding constraint value.")
675 complete = True
676 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0),
677 c0, x, revstack, complete, status, message)
678 else: # test_zero_row_2
679 # if LHS is >= 0, we can eliminate this constraint entirely
680 A_ub = A_ub[np.logical_not(zero_row), :]
681 b_ub = b_ub[np.logical_not(zero_row)]
683 # zero column in (both) constraints
684 # this indicates that a variable isn't constrained and can be removed
685 A = vstack((A_eq, A_ub))
686 if A.shape[0] > 0:
687 zero_col = np.array(np.sum(A != 0, axis=0) == 0).flatten()
688 # variable will be at upper or lower bound, depending on objective
689 x[np.logical_and(zero_col, c < 0)] = ub[
690 np.logical_and(zero_col, c < 0)]
691 x[np.logical_and(zero_col, c > 0)] = lb[
692 np.logical_and(zero_col, c > 0)]
693 if np.any(np.isinf(x)): # if an unconstrained variable has no bound
694 status = 3
695 message = ("If feasible, the problem is (trivially) unbounded "
696 "due to a zero column in the constraint matrices. If "
697 "you wish to check whether the problem is infeasible, "
698 "turn presolve off.")
699 complete = True
700 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0),
701 c0, x, revstack, complete, status, message)
702 # variables will equal upper/lower bounds will be removed later
703 lb[np.logical_and(zero_col, c < 0)] = ub[
704 np.logical_and(zero_col, c < 0)]
705 ub[np.logical_and(zero_col, c > 0)] = lb[
706 np.logical_and(zero_col, c > 0)]
708 # row singleton in equality constraints
709 # this fixes a variable and removes the constraint
710 singleton_row = np.array(np.sum(A_eq != 0, axis=1) == 1).flatten()
711 rows = where(singleton_row)[0]
712 cols = where(A_eq[rows, :])[1]
713 if len(rows) > 0:
714 for row, col in zip(rows, cols):
715 val = b_eq[row] / A_eq[row, col]
716 if not lb[col] - tol <= val <= ub[col] + tol:
717 # infeasible if fixed value is not within bounds
718 status = 2
719 message = ("The problem is (trivially) infeasible because a "
720 "singleton row in the equality constraints is "
721 "inconsistent with the bounds.")
722 complete = True
723 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0),
724 c0, x, revstack, complete, status, message)
725 else:
726 # sets upper and lower bounds at that fixed value - variable
727 # will be removed later
728 lb[col] = val
729 ub[col] = val
730 A_eq = A_eq[np.logical_not(singleton_row), :]
731 b_eq = b_eq[np.logical_not(singleton_row)]
733 # row singleton in inequality constraints
734 # this indicates a simple bound and the constraint can be removed
735 # simple bounds may be adjusted here
736 # After all of the simple bound information is combined here, get_Abc will
737 # turn the simple bounds into constraints
738 singleton_row = np.array(np.sum(A_ub != 0, axis=1) == 1).flatten()
739 cols = where(A_ub[singleton_row, :])[1]
740 rows = where(singleton_row)[0]
741 if len(rows) > 0:
742 for row, col in zip(rows, cols):
743 val = b_ub[row] / A_ub[row, col]
744 if A_ub[row, col] > 0: # upper bound
745 if val < lb[col] - tol: # infeasible
746 complete = True
747 elif val < ub[col]: # new upper bound
748 ub[col] = val
749 else: # lower bound
750 if val > ub[col] + tol: # infeasible
751 complete = True
752 elif val > lb[col]: # new lower bound
753 lb[col] = val
754 if complete:
755 status = 2
756 message = ("The problem is (trivially) infeasible because a "
757 "singleton row in the upper bound constraints is "
758 "inconsistent with the bounds.")
759 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0),
760 c0, x, revstack, complete, status, message)
761 A_ub = A_ub[np.logical_not(singleton_row), :]
762 b_ub = b_ub[np.logical_not(singleton_row)]
764 # identical bounds indicate that variable can be removed
765 i_f = np.abs(lb - ub) < tol # indices of "fixed" variables
766 i_nf = np.logical_not(i_f) # indices of "not fixed" variables
768 # test_bounds_equal_but_infeasible
769 if np.all(i_f): # if bounds define solution, check for consistency
770 residual = b_eq - A_eq.dot(lb)
771 slack = b_ub - A_ub.dot(lb)
772 if ((A_ub.size > 0 and np.any(slack < 0)) or
773 (A_eq.size > 0 and not np.allclose(residual, 0))):
774 status = 2
775 message = ("The problem is (trivially) infeasible because the "
776 "bounds fix all variables to values inconsistent with "
777 "the constraints")
778 complete = True
779 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0),
780 c0, x, revstack, complete, status, message)
782 ub_mod = ub
783 lb_mod = lb
784 if np.any(i_f):
785 c0 += c[i_f].dot(lb[i_f])
786 b_eq = b_eq - A_eq[:, i_f].dot(lb[i_f])
787 b_ub = b_ub - A_ub[:, i_f].dot(lb[i_f])
788 c = c[i_nf]
789 x_undo = lb[i_f] # not x[i_f], x is just zeroes
790 x = x[i_nf]
791 # user guess x0 stays separate from presolve solution x
792 if x0 is not None:
793 x0 = x0[i_nf]
794 A_eq = A_eq[:, i_nf]
795 A_ub = A_ub[:, i_nf]
796 # modify bounds
797 lb_mod = lb[i_nf]
798 ub_mod = ub[i_nf]
800 def rev(x_mod):
801 # Function to restore x: insert x_undo into x_mod.
802 # When elements have been removed at positions k1, k2, k3, ...
803 # then these must be replaced at (after) positions k1-1, k2-2,
804 # k3-3, ... in the modified array to recreate the original
805 i = np.flatnonzero(i_f)
806 # Number of variables to restore
807 N = len(i)
808 index_offset = np.arange(N)
809 # Create insert indices
810 insert_indices = i - index_offset
811 x_rev = np.insert(x_mod.astype(float), insert_indices, x_undo)
812 return x_rev
814 # Use revstack as a list of functions, currently just this one.
815 revstack.append(rev)
817 # no constraints indicates that problem is trivial
818 if A_eq.size == 0 and A_ub.size == 0:
819 b_eq = np.array([])
820 b_ub = np.array([])
821 # test_empty_constraint_1
822 if c.size == 0:
823 status = 0
824 message = ("The solution was determined in presolve as there are "
825 "no non-trivial constraints.")
826 elif (np.any(np.logical_and(c < 0, ub_mod == np.inf)) or
827 np.any(np.logical_and(c > 0, lb_mod == -np.inf))):
828 # test_no_constraints()
829 # test_unbounded_no_nontrivial_constraints_1
830 # test_unbounded_no_nontrivial_constraints_2
831 status = 3
832 message = ("The problem is (trivially) unbounded "
833 "because there are no non-trivial constraints and "
834 "a) at least one decision variable is unbounded "
835 "above and its corresponding cost is negative, or "
836 "b) at least one decision variable is unbounded below "
837 "and its corresponding cost is positive. ")
838 else: # test_empty_constraint_2
839 status = 0
840 message = ("The solution was determined in presolve as there are "
841 "no non-trivial constraints.")
842 complete = True
843 x[c < 0] = ub_mod[c < 0]
844 x[c > 0] = lb_mod[c > 0]
845 # where c is zero, set x to a finite bound or zero
846 x_zero_c = ub_mod[c == 0]
847 x_zero_c[np.isinf(x_zero_c)] = ub_mod[c == 0][np.isinf(x_zero_c)]
848 x_zero_c[np.isinf(x_zero_c)] = 0
849 x[c == 0] = x_zero_c
850 # if this is not the last step of presolve, should convert bounds back
851 # to array and return here
853 # Convert modified lb and ub back into N x 2 bounds
854 bounds = np.hstack((lb_mod[:, np.newaxis], ub_mod[:, np.newaxis]))
856 # remove redundant (linearly dependent) rows from equality constraints
857 n_rows_A = A_eq.shape[0]
858 redundancy_warning = ("A_eq does not appear to be of full row rank. To "
859 "improve performance, check the problem formulation "
860 "for redundant equality constraints.")
861 if (sps.issparse(A_eq)):
862 if rr and A_eq.size > 0: # TODO: Fast sparse rank check?
863 rr_res = _remove_redundancy_pivot_sparse(A_eq, b_eq)
864 A_eq, b_eq, status, message = rr_res
865 if A_eq.shape[0] < n_rows_A:
866 warn(redundancy_warning, OptimizeWarning, stacklevel=1)
867 if status != 0:
868 complete = True
869 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0),
870 c0, x, revstack, complete, status, message)
872 # This is a wild guess for which redundancy removal algorithm will be
873 # faster. More testing would be good.
874 small_nullspace = 5
875 if rr and A_eq.size > 0:
876 try: # TODO: use results of first SVD in _remove_redundancy_svd
877 rank = np.linalg.matrix_rank(A_eq)
878 # oh well, we'll have to go with _remove_redundancy_pivot_dense
879 except Exception:
880 rank = 0
881 if rr and A_eq.size > 0 and rank < A_eq.shape[0]:
882 warn(redundancy_warning, OptimizeWarning, stacklevel=3)
883 dim_row_nullspace = A_eq.shape[0]-rank
884 if rr_method is None:
885 if dim_row_nullspace <= small_nullspace:
886 rr_res = _remove_redundancy_svd(A_eq, b_eq)
887 A_eq, b_eq, status, message = rr_res
888 if dim_row_nullspace > small_nullspace or status == 4:
889 rr_res = _remove_redundancy_pivot_dense(A_eq, b_eq)
890 A_eq, b_eq, status, message = rr_res
892 else:
893 rr_method = rr_method.lower()
894 if rr_method == "svd":
895 rr_res = _remove_redundancy_svd(A_eq, b_eq)
896 A_eq, b_eq, status, message = rr_res
897 elif rr_method == "pivot":
898 rr_res = _remove_redundancy_pivot_dense(A_eq, b_eq)
899 A_eq, b_eq, status, message = rr_res
900 elif rr_method == "id":
901 rr_res = _remove_redundancy_id(A_eq, b_eq, rank)
902 A_eq, b_eq, status, message = rr_res
903 else: # shouldn't get here; option validity checked above
904 pass
905 if A_eq.shape[0] < rank:
906 message = ("Due to numerical issues, redundant equality "
907 "constraints could not be removed automatically. "
908 "Try providing your constraint matrices as sparse "
909 "matrices to activate sparse presolve, try turning "
910 "off redundancy removal, or try turning off presolve "
911 "altogether.")
912 status = 4
913 if status != 0:
914 complete = True
915 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0),
916 c0, x, revstack, complete, status, message)
919def _parse_linprog(lp, options, meth):
920 """
921 Parse the provided linear programming problem
923 ``_parse_linprog`` employs two main steps ``_check_sparse_inputs`` and
924 ``_clean_inputs``. ``_check_sparse_inputs`` checks for sparsity in the
925 provided constraints (``A_ub`` and ``A_eq) and if these match the provided
926 sparsity optional values.
928 ``_clean inputs`` checks of the provided inputs. If no violations are
929 identified the objective vector, upper bound constraints, equality
930 constraints, and simple bounds are returned in the expected format.
932 Parameters
933 ----------
934 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields:
936 c : 1D array
937 The coefficients of the linear objective function to be minimized.
938 A_ub : 2D array, optional
939 The inequality constraint matrix. Each row of ``A_ub`` specifies the
940 coefficients of a linear inequality constraint on ``x``.
941 b_ub : 1D array, optional
942 The inequality constraint vector. Each element represents an
943 upper bound on the corresponding value of ``A_ub @ x``.
944 A_eq : 2D array, optional
945 The equality constraint matrix. Each row of ``A_eq`` specifies the
946 coefficients of a linear equality constraint on ``x``.
947 b_eq : 1D array, optional
948 The equality constraint vector. Each element of ``A_eq @ x`` must equal
949 the corresponding element of ``b_eq``.
950 bounds : various valid formats, optional
951 The bounds of ``x``, as ``min`` and ``max`` pairs.
952 If bounds are specified for all N variables separately, valid formats are:
953 * a 2D array (2 x N or N x 2);
954 * a sequence of N sequences, each with 2 values.
955 If all variables have the same bounds, a single pair of values can
956 be specified. Valid formats are:
957 * a sequence with 2 scalar values;
958 * a sequence with a single element containing 2 scalar values.
959 If all variables have a lower bound of 0 and no upper bound, the bounds
960 parameter can be omitted (or given as None).
961 x0 : 1D array, optional
962 Guess values of the decision variables, which will be refined by
963 the optimization algorithm. This argument is currently used only by the
964 'revised simplex' method, and can only be used if `x0` represents a
965 basic feasible solution.
967 options : dict
968 A dictionary of solver options. All methods accept the following
969 generic options:
971 maxiter : int
972 Maximum number of iterations to perform.
973 disp : bool
974 Set to True to print convergence messages.
976 For method-specific options, see :func:`show_options('linprog')`.
978 Returns
979 -------
980 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields:
982 c : 1D array
983 The coefficients of the linear objective function to be minimized.
984 A_ub : 2D array, optional
985 The inequality constraint matrix. Each row of ``A_ub`` specifies the
986 coefficients of a linear inequality constraint on ``x``.
987 b_ub : 1D array, optional
988 The inequality constraint vector. Each element represents an
989 upper bound on the corresponding value of ``A_ub @ x``.
990 A_eq : 2D array, optional
991 The equality constraint matrix. Each row of ``A_eq`` specifies the
992 coefficients of a linear equality constraint on ``x``.
993 b_eq : 1D array, optional
994 The equality constraint vector. Each element of ``A_eq @ x`` must equal
995 the corresponding element of ``b_eq``.
996 bounds : 2D array
997 The bounds of ``x``, as ``min`` and ``max`` pairs, one for each of the N
998 elements of ``x``. The N x 2 array contains lower bounds in the first
999 column and upper bounds in the 2nd. Unbounded variables have lower
1000 bound -np.inf and/or upper bound np.inf.
1001 x0 : 1D array, optional
1002 Guess values of the decision variables, which will be refined by
1003 the optimization algorithm. This argument is currently used only by the
1004 'revised simplex' method, and can only be used if `x0` represents a
1005 basic feasible solution.
1007 options : dict, optional
1008 A dictionary of solver options. All methods accept the following
1009 generic options:
1011 maxiter : int
1012 Maximum number of iterations to perform.
1013 disp : bool
1014 Set to True to print convergence messages.
1016 For method-specific options, see :func:`show_options('linprog')`.
1018 """
1019 if options is None:
1020 options = {}
1022 solver_options = {k: v for k, v in options.items()}
1023 solver_options, A_ub, A_eq = _check_sparse_inputs(solver_options, meth,
1024 lp.A_ub, lp.A_eq)
1025 # Convert lists to numpy arrays, etc...
1026 lp = _clean_inputs(lp._replace(A_ub=A_ub, A_eq=A_eq))
1027 return lp, solver_options
1030def _get_Abc(lp, c0):
1031 """
1032 Given a linear programming problem of the form:
1034 Minimize::
1036 c @ x
1038 Subject to::
1040 A_ub @ x <= b_ub
1041 A_eq @ x == b_eq
1042 lb <= x <= ub
1044 where ``lb = 0`` and ``ub = None`` unless set in ``bounds``.
1046 Return the problem in standard form:
1048 Minimize::
1050 c @ x
1052 Subject to::
1054 A @ x == b
1055 x >= 0
1057 by adding slack variables and making variable substitutions as necessary.
1059 Parameters
1060 ----------
1061 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields:
1063 c : 1D array
1064 The coefficients of the linear objective function to be minimized.
1065 A_ub : 2D array, optional
1066 The inequality constraint matrix. Each row of ``A_ub`` specifies the
1067 coefficients of a linear inequality constraint on ``x``.
1068 b_ub : 1D array, optional
1069 The inequality constraint vector. Each element represents an
1070 upper bound on the corresponding value of ``A_ub @ x``.
1071 A_eq : 2D array, optional
1072 The equality constraint matrix. Each row of ``A_eq`` specifies the
1073 coefficients of a linear equality constraint on ``x``.
1074 b_eq : 1D array, optional
1075 The equality constraint vector. Each element of ``A_eq @ x`` must equal
1076 the corresponding element of ``b_eq``.
1077 bounds : 2D array
1078 The bounds of ``x``, lower bounds in the 1st column, upper
1079 bounds in the 2nd column. The bounds are possibly tightened
1080 by the presolve procedure.
1081 x0 : 1D array, optional
1082 Guess values of the decision variables, which will be refined by
1083 the optimization algorithm. This argument is currently used only by the
1084 'revised simplex' method, and can only be used if `x0` represents a
1085 basic feasible solution.
1087 c0 : float
1088 Constant term in objective function due to fixed (and eliminated)
1089 variables.
1091 Returns
1092 -------
1093 A : 2-D array
1094 2-D array such that ``A`` @ ``x``, gives the values of the equality
1095 constraints at ``x``.
1096 b : 1-D array
1097 1-D array of values representing the RHS of each equality constraint
1098 (row) in A (for standard form problem).
1099 c : 1-D array
1100 Coefficients of the linear objective function to be minimized (for
1101 standard form problem).
1102 c0 : float
1103 Constant term in objective function due to fixed (and eliminated)
1104 variables.
1105 x0 : 1-D array
1106 Starting values of the independent variables, which will be refined by
1107 the optimization algorithm
1109 References
1110 ----------
1111 .. [9] Bertsimas, Dimitris, and J. Tsitsiklis. "Introduction to linear
1112 programming." Athena Scientific 1 (1997): 997.
1114 """
1115 c, A_ub, b_ub, A_eq, b_eq, bounds, x0, integrality = lp
1117 if sps.issparse(A_eq):
1118 sparse = True
1119 A_eq = sps.csr_matrix(A_eq)
1120 A_ub = sps.csr_matrix(A_ub)
1122 def hstack(blocks):
1123 return sps.hstack(blocks, format="csr")
1125 def vstack(blocks):
1126 return sps.vstack(blocks, format="csr")
1128 zeros = sps.csr_matrix
1129 eye = sps.eye
1130 else:
1131 sparse = False
1132 hstack = np.hstack
1133 vstack = np.vstack
1134 zeros = np.zeros
1135 eye = np.eye
1137 # Variables lbs and ubs (see below) may be changed, which feeds back into
1138 # bounds, so copy.
1139 bounds = np.array(bounds, copy=True)
1141 # modify problem such that all variables have only non-negativity bounds
1142 lbs = bounds[:, 0]
1143 ubs = bounds[:, 1]
1144 m_ub, n_ub = A_ub.shape
1146 lb_none = np.equal(lbs, -np.inf)
1147 ub_none = np.equal(ubs, np.inf)
1148 lb_some = np.logical_not(lb_none)
1149 ub_some = np.logical_not(ub_none)
1151 # unbounded below: substitute xi = -xi' (unbounded above)
1152 # if -inf <= xi <= ub, then -ub <= -xi <= inf, so swap and invert bounds
1153 l_nolb_someub = np.logical_and(lb_none, ub_some)
1154 i_nolb = np.nonzero(l_nolb_someub)[0]
1155 lbs[l_nolb_someub], ubs[l_nolb_someub] = (
1156 -ubs[l_nolb_someub], -lbs[l_nolb_someub])
1157 lb_none = np.equal(lbs, -np.inf)
1158 ub_none = np.equal(ubs, np.inf)
1159 lb_some = np.logical_not(lb_none)
1160 ub_some = np.logical_not(ub_none)
1161 c[i_nolb] *= -1
1162 if x0 is not None:
1163 x0[i_nolb] *= -1
1164 if len(i_nolb) > 0:
1165 if A_ub.shape[0] > 0: # sometimes needed for sparse arrays... weird
1166 A_ub[:, i_nolb] *= -1
1167 if A_eq.shape[0] > 0:
1168 A_eq[:, i_nolb] *= -1
1170 # upper bound: add inequality constraint
1171 i_newub, = ub_some.nonzero()
1172 ub_newub = ubs[ub_some]
1173 n_bounds = len(i_newub)
1174 if n_bounds > 0:
1175 shape = (n_bounds, A_ub.shape[1])
1176 if sparse:
1177 idxs = (np.arange(n_bounds), i_newub)
1178 A_ub = vstack((A_ub, sps.csr_matrix((np.ones(n_bounds), idxs),
1179 shape=shape)))
1180 else:
1181 A_ub = vstack((A_ub, np.zeros(shape)))
1182 A_ub[np.arange(m_ub, A_ub.shape[0]), i_newub] = 1
1183 b_ub = np.concatenate((b_ub, np.zeros(n_bounds)))
1184 b_ub[m_ub:] = ub_newub
1186 A1 = vstack((A_ub, A_eq))
1187 b = np.concatenate((b_ub, b_eq))
1188 c = np.concatenate((c, np.zeros((A_ub.shape[0],))))
1189 if x0 is not None:
1190 x0 = np.concatenate((x0, np.zeros((A_ub.shape[0],))))
1191 # unbounded: substitute xi = xi+ + xi-
1192 l_free = np.logical_and(lb_none, ub_none)
1193 i_free = np.nonzero(l_free)[0]
1194 n_free = len(i_free)
1195 c = np.concatenate((c, np.zeros(n_free)))
1196 if x0 is not None:
1197 x0 = np.concatenate((x0, np.zeros(n_free)))
1198 A1 = hstack((A1[:, :n_ub], -A1[:, i_free]))
1199 c[n_ub:n_ub+n_free] = -c[i_free]
1200 if x0 is not None:
1201 i_free_neg = x0[i_free] < 0
1202 x0[np.arange(n_ub, A1.shape[1])[i_free_neg]] = -x0[i_free[i_free_neg]]
1203 x0[i_free[i_free_neg]] = 0
1205 # add slack variables
1206 A2 = vstack([eye(A_ub.shape[0]), zeros((A_eq.shape[0], A_ub.shape[0]))])
1208 A = hstack([A1, A2])
1210 # lower bound: substitute xi = xi' + lb
1211 # now there is a constant term in objective
1212 i_shift = np.nonzero(lb_some)[0]
1213 lb_shift = lbs[lb_some].astype(float)
1214 c0 += np.sum(lb_shift * c[i_shift])
1215 if sparse:
1216 b = b.reshape(-1, 1)
1217 A = A.tocsc()
1218 b -= (A[:, i_shift] * sps.diags(lb_shift)).sum(axis=1)
1219 b = b.ravel()
1220 else:
1221 b -= (A[:, i_shift] * lb_shift).sum(axis=1)
1222 if x0 is not None:
1223 x0[i_shift] -= lb_shift
1225 return A, b, c, c0, x0
1228def _round_to_power_of_two(x):
1229 """
1230 Round elements of the array to the nearest power of two.
1231 """
1232 return 2**np.around(np.log2(x))
1235def _autoscale(A, b, c, x0):
1236 """
1237 Scales the problem according to equilibration from [12].
1238 Also normalizes the right hand side vector by its maximum element.
1239 """
1240 m, n = A.shape
1242 C = 1
1243 R = 1
1245 if A.size > 0:
1247 R = np.max(np.abs(A), axis=1)
1248 if sps.issparse(A):
1249 R = R.toarray().flatten()
1250 R[R == 0] = 1
1251 R = 1/_round_to_power_of_two(R)
1252 A = sps.diags(R)*A if sps.issparse(A) else A*R.reshape(m, 1)
1253 b = b*R
1255 C = np.max(np.abs(A), axis=0)
1256 if sps.issparse(A):
1257 C = C.toarray().flatten()
1258 C[C == 0] = 1
1259 C = 1/_round_to_power_of_two(C)
1260 A = A*sps.diags(C) if sps.issparse(A) else A*C
1261 c = c*C
1263 b_scale = np.max(np.abs(b)) if b.size > 0 else 1
1264 if b_scale == 0:
1265 b_scale = 1.
1266 b = b/b_scale
1268 if x0 is not None:
1269 x0 = x0/b_scale*(1/C)
1270 return A, b, c, x0, C, b_scale
1273def _unscale(x, C, b_scale):
1274 """
1275 Converts solution to _autoscale problem -> solution to original problem.
1276 """
1278 try:
1279 n = len(C)
1280 # fails if sparse or scalar; that's OK.
1281 # this is only needed for original simplex (never sparse)
1282 except TypeError:
1283 n = len(x)
1285 return x[:n]*b_scale*C
1288def _display_summary(message, status, fun, iteration):
1289 """
1290 Print the termination summary of the linear program
1292 Parameters
1293 ----------
1294 message : str
1295 A string descriptor of the exit status of the optimization.
1296 status : int
1297 An integer representing the exit status of the optimization::
1299 0 : Optimization terminated successfully
1300 1 : Iteration limit reached
1301 2 : Problem appears to be infeasible
1302 3 : Problem appears to be unbounded
1303 4 : Serious numerical difficulties encountered
1305 fun : float
1306 Value of the objective function.
1307 iteration : iteration
1308 The number of iterations performed.
1309 """
1310 print(message)
1311 if status in (0, 1):
1312 print(" Current function value: {0: <12.6f}".format(fun))
1313 print(" Iterations: {0:d}".format(iteration))
1316def _postsolve(x, postsolve_args, complete=False):
1317 """
1318 Given solution x to presolved, standard form linear program x, add
1319 fixed variables back into the problem and undo the variable substitutions
1320 to get solution to original linear program. Also, calculate the objective
1321 function value, slack in original upper bound constraints, and residuals
1322 in original equality constraints.
1324 Parameters
1325 ----------
1326 x : 1-D array
1327 Solution vector to the standard-form problem.
1328 postsolve_args : tuple
1329 Data needed by _postsolve to convert the solution to the standard-form
1330 problem into the solution to the original problem, including:
1332 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields:
1334 c : 1D array
1335 The coefficients of the linear objective function to be minimized.
1336 A_ub : 2D array, optional
1337 The inequality constraint matrix. Each row of ``A_ub`` specifies the
1338 coefficients of a linear inequality constraint on ``x``.
1339 b_ub : 1D array, optional
1340 The inequality constraint vector. Each element represents an
1341 upper bound on the corresponding value of ``A_ub @ x``.
1342 A_eq : 2D array, optional
1343 The equality constraint matrix. Each row of ``A_eq`` specifies the
1344 coefficients of a linear equality constraint on ``x``.
1345 b_eq : 1D array, optional
1346 The equality constraint vector. Each element of ``A_eq @ x`` must equal
1347 the corresponding element of ``b_eq``.
1348 bounds : 2D array
1349 The bounds of ``x``, lower bounds in the 1st column, upper
1350 bounds in the 2nd column. The bounds are possibly tightened
1351 by the presolve procedure.
1352 x0 : 1D array, optional
1353 Guess values of the decision variables, which will be refined by
1354 the optimization algorithm. This argument is currently used only by the
1355 'revised simplex' method, and can only be used if `x0` represents a
1356 basic feasible solution.
1358 revstack: list of functions
1359 the functions in the list reverse the operations of _presolve()
1360 the function signature is x_org = f(x_mod), where x_mod is the result
1361 of a presolve step and x_org the value at the start of the step
1362 complete : bool
1363 Whether the solution is was determined in presolve (``True`` if so)
1365 Returns
1366 -------
1367 x : 1-D array
1368 Solution vector to original linear programming problem
1369 fun: float
1370 optimal objective value for original problem
1371 slack : 1-D array
1372 The (non-negative) slack in the upper bound constraints, that is,
1373 ``b_ub - A_ub @ x``
1374 con : 1-D array
1375 The (nominally zero) residuals of the equality constraints, that is,
1376 ``b - A_eq @ x``
1377 """
1378 # note that all the inputs are the ORIGINAL, unmodified versions
1379 # no rows, columns have been removed
1381 c, A_ub, b_ub, A_eq, b_eq, bounds, x0, integrality = postsolve_args[0]
1382 revstack, C, b_scale = postsolve_args[1:]
1384 x = _unscale(x, C, b_scale)
1386 # Undo variable substitutions of _get_Abc()
1387 # if "complete", problem was solved in presolve; don't do anything here
1388 n_x = bounds.shape[0]
1389 if not complete and bounds is not None: # bounds are never none, probably
1390 n_unbounded = 0
1391 for i, bi in enumerate(bounds):
1392 lbi = bi[0]
1393 ubi = bi[1]
1394 if lbi == -np.inf and ubi == np.inf:
1395 n_unbounded += 1
1396 x[i] = x[i] - x[n_x + n_unbounded - 1]
1397 else:
1398 if lbi == -np.inf:
1399 x[i] = ubi - x[i]
1400 else:
1401 x[i] += lbi
1402 # all the rest of the variables were artificial
1403 x = x[:n_x]
1405 # If there were variables removed from the problem, add them back into the
1406 # solution vector
1407 # Apply the functions in revstack (reverse direction)
1408 for rev in reversed(revstack):
1409 x = rev(x)
1411 fun = x.dot(c)
1412 slack = b_ub - A_ub.dot(x) # report slack for ORIGINAL UB constraints
1413 # report residuals of ORIGINAL EQ constraints
1414 con = b_eq - A_eq.dot(x)
1416 return x, fun, slack, con
1419def _check_result(x, fun, status, slack, con, bounds, tol, message):
1420 """
1421 Check the validity of the provided solution.
1423 A valid (optimal) solution satisfies all bounds, all slack variables are
1424 negative and all equality constraint residuals are strictly non-zero.
1425 Further, the lower-bounds, upper-bounds, slack and residuals contain
1426 no nan values.
1428 Parameters
1429 ----------
1430 x : 1-D array
1431 Solution vector to original linear programming problem
1432 fun: float
1433 optimal objective value for original problem
1434 status : int
1435 An integer representing the exit status of the optimization::
1437 0 : Optimization terminated successfully
1438 1 : Iteration limit reached
1439 2 : Problem appears to be infeasible
1440 3 : Problem appears to be unbounded
1441 4 : Serious numerical difficulties encountered
1443 slack : 1-D array
1444 The (non-negative) slack in the upper bound constraints, that is,
1445 ``b_ub - A_ub @ x``
1446 con : 1-D array
1447 The (nominally zero) residuals of the equality constraints, that is,
1448 ``b - A_eq @ x``
1449 bounds : 2D array
1450 The bounds on the original variables ``x``
1451 message : str
1452 A string descriptor of the exit status of the optimization.
1453 tol : float
1454 Termination tolerance; see [1]_ Section 4.5.
1456 Returns
1457 -------
1458 status : int
1459 An integer representing the exit status of the optimization::
1461 0 : Optimization terminated successfully
1462 1 : Iteration limit reached
1463 2 : Problem appears to be infeasible
1464 3 : Problem appears to be unbounded
1465 4 : Serious numerical difficulties encountered
1467 message : str
1468 A string descriptor of the exit status of the optimization.
1469 """
1470 # Somewhat arbitrary
1471 tol = np.sqrt(tol) * 10
1473 if x is None:
1474 # HiGHS does not provide x if infeasible/unbounded
1475 if status == 0: # Observed with HiGHS Simplex Primal
1476 status = 4
1477 message = ("The solver did not provide a solution nor did it "
1478 "report a failure. Please submit a bug report.")
1479 return status, message
1481 contains_nans = (
1482 np.isnan(x).any()
1483 or np.isnan(fun)
1484 or np.isnan(slack).any()
1485 or np.isnan(con).any()
1486 )
1488 if contains_nans:
1489 is_feasible = False
1490 else:
1491 invalid_bounds = (x < bounds[:, 0] - tol).any() or (x > bounds[:, 1] + tol).any()
1492 invalid_slack = status != 3 and (slack < -tol).any()
1493 invalid_con = status != 3 and (np.abs(con) > tol).any()
1494 is_feasible = not (invalid_bounds or invalid_slack or invalid_con)
1496 if status == 0 and not is_feasible:
1497 status = 4
1498 message = ("The solution does not satisfy the constraints within the "
1499 "required tolerance of " + "{:.2E}".format(tol) + ", yet "
1500 "no errors were raised and there is no certificate of "
1501 "infeasibility or unboundedness. Check whether "
1502 "the slack and constraint residuals are acceptable; "
1503 "if not, consider enabling presolve, adjusting the "
1504 "tolerance option(s), and/or using a different method. "
1505 "Please consider submitting a bug report.")
1506 elif status == 2 and is_feasible:
1507 # Occurs if the simplex method exits after phase one with a very
1508 # nearly basic feasible solution. Postsolving can make the solution
1509 # basic, however, this solution is NOT optimal
1510 status = 4
1511 message = ("The solution is feasible, but the solver did not report "
1512 "that the solution was optimal. Please try a different "
1513 "method.")
1515 return status, message