Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1""" 

2Method agnostic utility functions for linear progamming 

3""" 

4 

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, _remove_redundancy_sparse, _remove_redundancy_dense 

11 ) 

12from collections import namedtuple 

13 

14 

15_LPProblem = namedtuple('_LPProblem', 'c A_ub b_ub A_eq b_eq bounds x0') 

16_LPProblem.__new__.__defaults__ = (None,) * 6 # make c the only required arg 

17_LPProblem.__doc__ = \ 

18 """ Represents a linear-programming problem. 

19 

20 Attributes 

21 ---------- 

22 c : 1D array 

23 The coefficients of the linear objective function to be minimized. 

24 A_ub : 2D array, optional 

25 The inequality constraint matrix. Each row of ``A_ub`` specifies the 

26 coefficients of a linear inequality constraint on ``x``. 

27 b_ub : 1D array, optional 

28 The inequality constraint vector. Each element represents an 

29 upper bound on the corresponding value of ``A_ub @ x``. 

30 A_eq : 2D array, optional 

31 The equality constraint matrix. Each row of ``A_eq`` specifies the 

32 coefficients of a linear equality constraint on ``x``. 

33 b_eq : 1D array, optional 

34 The equality constraint vector. Each element of ``A_eq @ x`` must equal 

35 the corresponding element of ``b_eq``. 

36 bounds : various valid formats, optional 

37 The bounds of ``x``, as ``min`` and ``max`` pairs. 

38 If bounds are specified for all N variables separately, valid formats 

39 are: 

40 * a 2D array (N x 2); 

41 * a sequence of N sequences, each with 2 values. 

42 If all variables have the same bounds, the bounds can be specified as 

43 a 1-D or 2-D array or sequence with 2 scalar values. 

44 If all variables have a lower bound of 0 and no upper bound, the bounds 

45 parameter can be omitted (or given as None). 

46 Absent lower and/or upper bounds can be specified as -numpy.inf (no 

47 lower bound), numpy.inf (no upper bound) or None (both). 

48 x0 : 1D array, optional 

49 Guess values of the decision variables, which will be refined by 

50 the optimization algorithm. This argument is currently used only by the 

51 'revised simplex' method, and can only be used if `x0` represents a 

52 basic feasible solution. 

53 

54 Notes 

55 ----- 

56 This namedtuple supports 2 ways of initialization: 

57 >>> lp1 = _LPProblem(c=[-1, 4], A_ub=[[-3, 1], [1, 2]], b_ub=[6, 4]) 

58 >>> lp2 = _LPProblem([-1, 4], [[-3, 1], [1, 2]], [6, 4]) 

59 

60 Note that only ``c`` is a required argument here, whereas all other arguments 

61 ``A_ub``, ``b_ub``, ``A_eq``, ``b_eq``, ``bounds``, ``x0`` are optional with 

62 default values of None. 

63 For example, ``A_eq`` and ``b_eq`` can be set without ``A_ub`` or ``b_ub``: 

64 >>> lp3 = _LPProblem(c=[-1, 4], A_eq=[[2, 1]], b_eq=[10]) 

65 """ 

66 

67 

68def _check_sparse_inputs(options, A_ub, A_eq): 

69 """ 

70 Check the provided ``A_ub`` and ``A_eq`` matrices conform to the specified 

71 optional sparsity variables. 

72 

73 Parameters 

74 ---------- 

75 A_ub : 2-D array, optional 

76 2-D array such that ``A_ub @ x`` gives the values of the upper-bound 

77 inequality constraints at ``x``. 

78 A_eq : 2-D array, optional 

79 2-D array such that ``A_eq @ x`` gives the values of the equality 

80 constraints at ``x``. 

81 options : dict 

82 A dictionary of solver options. All methods accept the following 

83 generic options: 

84 

85 maxiter : int 

86 Maximum number of iterations to perform. 

87 disp : bool 

88 Set to True to print convergence messages. 

89 

90 For method-specific options, see :func:`show_options('linprog')`. 

91 

92 Returns 

93 ------- 

94 A_ub : 2-D array, optional 

95 2-D array such that ``A_ub @ x`` gives the values of the upper-bound 

96 inequality constraints at ``x``. 

97 A_eq : 2-D array, optional 

98 2-D array such that ``A_eq @ x`` gives the values of the equality 

99 constraints at ``x``. 

100 options : dict 

101 A dictionary of solver options. All methods accept the following 

102 generic options: 

103 

104 maxiter : int 

105 Maximum number of iterations to perform. 

106 disp : bool 

107 Set to True to print convergence messages. 

108 

109 For method-specific options, see :func:`show_options('linprog')`. 

110 """ 

111 # This is an undocumented option for unit testing sparse presolve 

112 _sparse_presolve = options.pop('_sparse_presolve', False) 

113 if _sparse_presolve and A_eq is not None: 

114 A_eq = sps.coo_matrix(A_eq) 

115 if _sparse_presolve and A_ub is not None: 

116 A_ub = sps.coo_matrix(A_ub) 

117 

118 sparse = options.get('sparse', False) 

119 if not sparse and (sps.issparse(A_eq) or sps.issparse(A_ub)): 

120 options['sparse'] = True 

121 warn("Sparse constraint matrix detected; setting 'sparse':True.", 

122 OptimizeWarning, stacklevel=4) 

123 return options, A_ub, A_eq 

124 

125 

126def _format_A_constraints(A, n_x, sparse_lhs=False): 

127 """Format the left hand side of the constraints to a 2-D array 

128 

129 Parameters 

130 ---------- 

131 A : 2-D array 

132 2-D array such that ``A @ x`` gives the values of the upper-bound 

133 (in)equality constraints at ``x``. 

134 n_x : int 

135 The number of variables in the linear programming problem. 

136 sparse_lhs : bool 

137 Whether either of `A_ub` or `A_eq` are sparse. If true return a 

138 coo_matrix instead of a numpy array. 

139 

140 Returns 

141 ------- 

142 np.ndarray or sparse.coo_matrix 

143 2-D array such that ``A @ x`` gives the values of the upper-bound 

144 (in)equality constraints at ``x``. 

145 

146 """ 

147 if sparse_lhs: 

148 return sps.coo_matrix( 

149 (0, n_x) if A is None else A, dtype=float, copy=True 

150 ) 

151 elif A is None: 

152 return np.zeros((0, n_x), dtype=float) 

153 else: 

154 return np.array(A, dtype=float, copy=True) 

155 

156 

157def _format_b_constraints(b): 

158 """Format the upper bounds of the constraints to a 1-D array 

159 

160 Parameters 

161 ---------- 

162 b : 1-D array 

163 1-D array of values representing the upper-bound of each (in)equality 

164 constraint (row) in ``A``. 

165 

166 Returns 

167 ------- 

168 1-D np.array 

169 1-D array of values representing the upper-bound of each (in)equality 

170 constraint (row) in ``A``. 

171 

172 """ 

173 if b is None: 

174 return np.array([], dtype=float) 

175 b = np.array(b, dtype=float, copy=True).squeeze() 

176 return b if b.size != 1 else b.reshape((-1)) 

177 

178 

179def _clean_inputs(lp): 

180 """ 

181 Given user inputs for a linear programming problem, return the 

182 objective vector, upper bound constraints, equality constraints, 

183 and simple bounds in a preferred format. 

184 

185 Parameters 

186 ---------- 

187 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 

188 

189 c : 1D array 

190 The coefficients of the linear objective function to be minimized. 

191 A_ub : 2D array, optional 

192 The inequality constraint matrix. Each row of ``A_ub`` specifies the 

193 coefficients of a linear inequality constraint on ``x``. 

194 b_ub : 1D array, optional 

195 The inequality constraint vector. Each element represents an 

196 upper bound on the corresponding value of ``A_ub @ x``. 

197 A_eq : 2D array, optional 

198 The equality constraint matrix. Each row of ``A_eq`` specifies the 

199 coefficients of a linear equality constraint on ``x``. 

200 b_eq : 1D array, optional 

201 The equality constraint vector. Each element of ``A_eq @ x`` must equal 

202 the corresponding element of ``b_eq``. 

203 bounds : various valid formats, optional 

204 The bounds of ``x``, as ``min`` and ``max`` pairs. 

205 If bounds are specified for all N variables separately, valid formats are: 

206 * a 2D array (2 x N or N x 2); 

207 * a sequence of N sequences, each with 2 values. 

208 If all variables have the same bounds, a single pair of values can 

209 be specified. Valid formats are: 

210 * a sequence with 2 scalar values; 

211 * a sequence with a single element containing 2 scalar values. 

212 If all variables have a lower bound of 0 and no upper bound, the bounds 

213 parameter can be omitted (or given as None). 

214 x0 : 1D array, optional 

215 Guess values of the decision variables, which will be refined by 

216 the optimization algorithm. This argument is currently used only by the 

217 'revised simplex' method, and can only be used if `x0` represents a 

218 basic feasible solution. 

219 

220 Returns 

221 ------- 

222 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 

223 

224 c : 1D array 

225 The coefficients of the linear objective function to be minimized. 

226 A_ub : 2D array, optional 

227 The inequality constraint matrix. Each row of ``A_ub`` specifies the 

228 coefficients of a linear inequality constraint on ``x``. 

229 b_ub : 1D array, optional 

230 The inequality constraint vector. Each element represents an 

231 upper bound on the corresponding value of ``A_ub @ x``. 

232 A_eq : 2D array, optional 

233 The equality constraint matrix. Each row of ``A_eq`` specifies the 

234 coefficients of a linear equality constraint on ``x``. 

235 b_eq : 1D array, optional 

236 The equality constraint vector. Each element of ``A_eq @ x`` must equal 

237 the corresponding element of ``b_eq``. 

238 bounds : 2D array 

239 The bounds of ``x``, as ``min`` and ``max`` pairs, one for each of the N 

240 elements of ``x``. The N x 2 array contains lower bounds in the first 

241 column and upper bounds in the 2nd. Unbounded variables have lower 

242 bound -np.inf and/or upper bound np.inf. 

243 x0 : 1D array, optional 

244 Guess values of the decision variables, which will be refined by 

245 the optimization algorithm. This argument is currently used only by the 

246 'revised simplex' method, and can only be used if `x0` represents a 

247 basic feasible solution. 

248 

249 """ 

250 c, A_ub, b_ub, A_eq, b_eq, bounds, x0 = lp 

251 

252 if c is None: 

253 raise TypeError 

254 

255 try: 

256 c = np.array(c, dtype=np.float, copy=True).squeeze() 

257 except ValueError: 

258 raise TypeError( 

259 "Invalid input for linprog: c must be a 1-D array of numerical " 

260 "coefficients") 

261 else: 

262 # If c is a single value, convert it to a 1-D array. 

263 if c.size == 1: 

264 c = c.reshape((-1)) 

265 

266 n_x = len(c) 

267 if n_x == 0 or len(c.shape) != 1: 

268 raise ValueError( 

269 "Invalid input for linprog: c must be a 1-D array and must " 

270 "not have more than one non-singleton dimension") 

271 if not(np.isfinite(c).all()): 

272 raise ValueError( 

273 "Invalid input for linprog: c must not contain values " 

274 "inf, nan, or None") 

275 

276 sparse_lhs = sps.issparse(A_eq) or sps.issparse(A_ub) 

277 try: 

278 A_ub = _format_A_constraints(A_ub, n_x, sparse_lhs=sparse_lhs) 

279 except ValueError: 

280 raise TypeError( 

281 "Invalid input for linprog: A_ub must be a 2-D array " 

282 "of numerical values") 

283 else: 

284 n_ub = A_ub.shape[0] 

285 if len(A_ub.shape) != 2 or A_ub.shape[1] != n_x: 

286 raise ValueError( 

287 "Invalid input for linprog: A_ub must have exactly two " 

288 "dimensions, and the number of columns in A_ub must be " 

289 "equal to the size of c") 

290 if (sps.issparse(A_ub) and not np.isfinite(A_ub.data).all() 

291 or not sps.issparse(A_ub) and not np.isfinite(A_ub).all()): 

292 raise ValueError( 

293 "Invalid input for linprog: A_ub must not contain values " 

294 "inf, nan, or None") 

295 

296 try: 

297 b_ub = _format_b_constraints(b_ub) 

298 except ValueError: 

299 raise TypeError( 

300 "Invalid input for linprog: b_ub must be a 1-D array of " 

301 "numerical values, each representing the upper bound of an " 

302 "inequality constraint (row) in A_ub") 

303 else: 

304 if b_ub.shape != (n_ub,): 

305 raise ValueError( 

306 "Invalid input for linprog: b_ub must be a 1-D array; b_ub " 

307 "must not have more than one non-singleton dimension and " 

308 "the number of rows in A_ub must equal the number of values " 

309 "in b_ub") 

310 if not(np.isfinite(b_ub).all()): 

311 raise ValueError( 

312 "Invalid input for linprog: b_ub must not contain values " 

313 "inf, nan, or None") 

314 

315 try: 

316 A_eq = _format_A_constraints(A_eq, n_x, sparse_lhs=sparse_lhs) 

317 except ValueError: 

318 raise TypeError( 

319 "Invalid input for linprog: A_eq must be a 2-D array " 

320 "of numerical values") 

321 else: 

322 n_eq = A_eq.shape[0] 

323 if len(A_eq.shape) != 2 or A_eq.shape[1] != n_x: 

324 raise ValueError( 

325 "Invalid input for linprog: A_eq must have exactly two " 

326 "dimensions, and the number of columns in A_eq must be " 

327 "equal to the size of c") 

328 

329 if (sps.issparse(A_eq) and not np.isfinite(A_eq.data).all() 

330 or not sps.issparse(A_eq) and not np.isfinite(A_eq).all()): 

331 raise ValueError( 

332 "Invalid input for linprog: A_eq must not contain values " 

333 "inf, nan, or None") 

334 

335 try: 

336 b_eq = _format_b_constraints(b_eq) 

337 except ValueError: 

338 raise TypeError( 

339 "Invalid input for linprog: b_eq must be a 1-D array of " 

340 "numerical values, each representing the upper bound of an " 

341 "inequality constraint (row) in A_eq") 

342 else: 

343 if b_eq.shape != (n_eq,): 

344 raise ValueError( 

345 "Invalid input for linprog: b_eq must be a 1-D array; b_eq " 

346 "must not have more than one non-singleton dimension and " 

347 "the number of rows in A_eq must equal the number of values " 

348 "in b_eq") 

349 if not(np.isfinite(b_eq).all()): 

350 raise ValueError( 

351 "Invalid input for linprog: b_eq must not contain values " 

352 "inf, nan, or None") 

353 

354 # x0 gives a (optional) starting solution to the solver. If x0 is None, 

355 # skip the checks. Initial solution will be generated automatically. 

356 if x0 is not None: 

357 try: 

358 x0 = np.array(x0, dtype=float, copy=True).squeeze() 

359 except ValueError: 

360 raise TypeError( 

361 "Invalid input for linprog: x0 must be a 1-D array of " 

362 "numerical coefficients") 

363 if x0.ndim == 0: 

364 x0 = x0.reshape((-1)) 

365 if len(x0) == 0 or x0.ndim != 1: 

366 raise ValueError( 

367 "Invalid input for linprog: x0 should be a 1-D array; it " 

368 "must not have more than one non-singleton dimension") 

369 if not x0.size == c.size: 

370 raise ValueError( 

371 "Invalid input for linprog: x0 and c should contain the " 

372 "same number of elements") 

373 if not np.isfinite(x0).all(): 

374 raise ValueError( 

375 "Invalid input for linprog: x0 must not contain values " 

376 "inf, nan, or None") 

377 

378 # Bounds can be one of these formats: 

379 # (1) a 2-D array or sequence, with shape N x 2 

380 # (2) a 1-D or 2-D sequence or array with 2 scalars 

381 # (3) None (or an empty sequence or array) 

382 # Unspecified bounds can be represented by None or (-)np.inf. 

383 # All formats are converted into a N x 2 np.array with (-)np.inf where 

384 # bounds are unspecified. 

385 

386 # Prepare clean bounds array 

387 bounds_clean = np.zeros((n_x, 2), dtype=float) 

388 

389 # Convert to a numpy array. 

390 # np.array(..,dtype=float) raises an error if dimensions are inconsistent 

391 # or if there are invalid data types in bounds. Just add a linprog prefix 

392 # to the error and re-raise. 

393 # Creating at least a 2-D array simplifies the cases to distinguish below. 

394 if bounds is None or np.array_equal(bounds, []) or np.array_equal(bounds, [[]]): 

395 bounds = (0, np.inf) 

396 try: 

397 bounds_conv = np.atleast_2d(np.array(bounds, dtype=float)) 

398 except ValueError as e: 

399 raise ValueError( 

400 "Invalid input for linprog: unable to interpret bounds, " 

401 "check values and dimensions: " + e.args[0]) 

402 except TypeError as e: 

403 raise TypeError( 

404 "Invalid input for linprog: unable to interpret bounds, " 

405 "check values and dimensions: " + e.args[0]) 

406 

407 # Check bounds options 

408 bsh = bounds_conv.shape 

409 if len(bsh) > 2: 

410 # Do not try to handle multidimensional bounds input 

411 raise ValueError( 

412 "Invalid input for linprog: provide a 2-D array for bounds, " 

413 "not a {:d}-D array.".format(len(bsh))) 

414 elif np.all(bsh == (n_x, 2)): 

415 # Regular N x 2 array 

416 bounds_clean = bounds_conv 

417 elif (np.all(bsh == (2, 1)) or np.all(bsh == (1, 2))): 

418 # 2 values: interpret as overall lower and upper bound 

419 bounds_flat = bounds_conv.flatten() 

420 bounds_clean[:, 0] = bounds_flat[0] 

421 bounds_clean[:, 1] = bounds_flat[1] 

422 elif np.all(bsh == (2, n_x)): 

423 # Reject a 2 x N array 

424 raise ValueError( 

425 "Invalid input for linprog: provide a {:d} x 2 array for bounds, " 

426 "not a 2 x {:d} array.".format(n_x, n_x)) 

427 else: 

428 raise ValueError( 

429 "Invalid input for linprog: unable to interpret bounds with this " 

430 "dimension tuple: {0}.".format(bsh)) 

431 

432 # The process above creates nan-s where the input specified None 

433 # Convert the nan-s in the 1st column to -np.inf and in the 2nd column 

434 # to np.inf 

435 i_none = np.isnan(bounds_clean[:, 0]) 

436 bounds_clean[i_none, 0] = -np.inf 

437 i_none = np.isnan(bounds_clean[:, 1]) 

438 bounds_clean[i_none, 1] = np.inf 

439 

440 return _LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds_clean, x0) 

441 

442 

443def _presolve(lp, rr, tol=1e-9): 

444 """ 

445 Given inputs for a linear programming problem in preferred format, 

446 presolve the problem: identify trivial infeasibilities, redundancies, 

447 and unboundedness, tighten bounds where possible, and eliminate fixed 

448 variables. 

449 

450 Parameters 

451 ---------- 

452 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 

453 

454 c : 1D array 

455 The coefficients of the linear objective function to be minimized. 

456 A_ub : 2D array, optional 

457 The inequality constraint matrix. Each row of ``A_ub`` specifies the 

458 coefficients of a linear inequality constraint on ``x``. 

459 b_ub : 1D array, optional 

460 The inequality constraint vector. Each element represents an 

461 upper bound on the corresponding value of ``A_ub @ x``. 

462 A_eq : 2D array, optional 

463 The equality constraint matrix. Each row of ``A_eq`` specifies the 

464 coefficients of a linear equality constraint on ``x``. 

465 b_eq : 1D array, optional 

466 The equality constraint vector. Each element of ``A_eq @ x`` must equal 

467 the corresponding element of ``b_eq``. 

468 bounds : 2D array 

469 The bounds of ``x``, as ``min`` and ``max`` pairs, one for each of the N 

470 elements of ``x``. The N x 2 array contains lower bounds in the first 

471 column and upper bounds in the 2nd. Unbounded variables have lower 

472 bound -np.inf and/or upper bound np.inf. 

473 x0 : 1D array, optional 

474 Guess values of the decision variables, which will be refined by 

475 the optimization algorithm. This argument is currently used only by the 

476 'revised simplex' method, and can only be used if `x0` represents a 

477 basic feasible solution. 

478 

479 rr : bool 

480 If ``True`` attempts to eliminate any redundant rows in ``A_eq``. 

481 Set False if ``A_eq`` is known to be of full row rank, or if you are 

482 looking for a potential speedup (at the expense of reliability). 

483 tol : float 

484 The tolerance which determines when a solution is "close enough" to 

485 zero in Phase 1 to be considered a basic feasible solution or close 

486 enough to positive to serve as an optimal solution. 

487 

488 Returns 

489 ------- 

490 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 

491 

492 c : 1D array 

493 The coefficients of the linear objective function to be minimized. 

494 A_ub : 2D array, optional 

495 The inequality constraint matrix. Each row of ``A_ub`` specifies the 

496 coefficients of a linear inequality constraint on ``x``. 

497 b_ub : 1D array, optional 

498 The inequality constraint vector. Each element represents an 

499 upper bound on the corresponding value of ``A_ub @ x``. 

500 A_eq : 2D array, optional 

501 The equality constraint matrix. Each row of ``A_eq`` specifies the 

502 coefficients of a linear equality constraint on ``x``. 

503 b_eq : 1D array, optional 

504 The equality constraint vector. Each element of ``A_eq @ x`` must equal 

505 the corresponding element of ``b_eq``. 

506 bounds : 2D array 

507 The bounds of ``x``, as ``min`` and ``max`` pairs, possibly tightened. 

508 x0 : 1D array, optional 

509 Guess values of the decision variables, which will be refined by 

510 the optimization algorithm. This argument is currently used only by the 

511 'revised simplex' method, and can only be used if `x0` represents a 

512 basic feasible solution. 

513 

514 c0 : 1D array 

515 Constant term in objective function due to fixed (and eliminated) 

516 variables. 

517 x : 1D array 

518 Solution vector (when the solution is trivial and can be determined 

519 in presolve) 

520 undo: list of tuples 

521 (index, value) pairs that record the original index and fixed value 

522 for each variable removed from the problem 

523 complete: bool 

524 Whether the solution is complete (solved or determined to be infeasible 

525 or unbounded in presolve) 

526 status : int 

527 An integer representing the exit status of the optimization:: 

528 

529 0 : Optimization terminated successfully 

530 1 : Iteration limit reached 

531 2 : Problem appears to be infeasible 

532 3 : Problem appears to be unbounded 

533 4 : Serious numerical difficulties encountered 

534 

535 message : str 

536 A string descriptor of the exit status of the optimization. 

537 

538 References 

539 ---------- 

540 .. [5] Andersen, Erling D. "Finding all linearly dependent rows in 

541 large-scale linear programming." Optimization Methods and Software 

542 6.3 (1995): 219-227. 

543 .. [8] Andersen, Erling D., and Knud D. Andersen. "Presolving in linear 

544 programming." Mathematical Programming 71.2 (1995): 221-245. 

545 

546 """ 

547 # ideas from Reference [5] by Andersen and Andersen 

548 # however, unlike the reference, this is performed before converting 

549 # problem to standard form 

550 # There are a few advantages: 

551 # * artificial variables have not been added, so matrices are smaller 

552 # * bounds have not been converted to constraints yet. (It is better to 

553 # do that after presolve because presolve may adjust the simple bounds.) 

554 # There are many improvements that can be made, namely: 

555 # * implement remaining checks from [5] 

556 # * loop presolve until no additional changes are made 

557 # * implement additional efficiency improvements in redundancy removal [2] 

558 

559 c, A_ub, b_ub, A_eq, b_eq, bounds, x0 = lp 

560 

561 undo = [] # record of variables eliminated from problem 

562 # constant term in cost function may be added if variables are eliminated 

563 c0 = 0 

564 complete = False # complete is True if detected infeasible/unbounded 

565 x = np.zeros(c.shape) # this is solution vector if completed in presolve 

566 

567 status = 0 # all OK unless determined otherwise 

568 message = "" 

569 

570 # Lower and upper bounds 

571 lb = bounds[:, 0] 

572 ub = bounds[:, 1] 

573 

574 m_eq, n = A_eq.shape 

575 m_ub, n = A_ub.shape 

576 

577 if sps.issparse(A_eq): 

578 A_eq = A_eq.tocsr() 

579 A_ub = A_ub.tocsr() 

580 

581 def where(A): 

582 return A.nonzero() 

583 

584 vstack = sps.vstack 

585 else: 

586 where = np.where 

587 vstack = np.vstack 

588 

589 # zero row in equality constraints 

590 zero_row = np.array(np.sum(A_eq != 0, axis=1) == 0).flatten() 

591 if np.any(zero_row): 

592 if np.any( 

593 np.logical_and( 

594 zero_row, 

595 np.abs(b_eq) > tol)): # test_zero_row_1 

596 # infeasible if RHS is not zero 

597 status = 2 

598 message = ("The problem is (trivially) infeasible due to a row " 

599 "of zeros in the equality constraint matrix with a " 

600 "nonzero corresponding constraint value.") 

601 complete = True 

602 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 

603 c0, x, undo, complete, status, message) 

604 else: # test_zero_row_2 

605 # if RHS is zero, we can eliminate this equation entirely 

606 A_eq = A_eq[np.logical_not(zero_row), :] 

607 b_eq = b_eq[np.logical_not(zero_row)] 

608 

609 # zero row in inequality constraints 

610 zero_row = np.array(np.sum(A_ub != 0, axis=1) == 0).flatten() 

611 if np.any(zero_row): 

612 if np.any(np.logical_and(zero_row, b_ub < -tol)): # test_zero_row_1 

613 # infeasible if RHS is less than zero (because LHS is zero) 

614 status = 2 

615 message = ("The problem is (trivially) infeasible due to a row " 

616 "of zeros in the equality constraint matrix with a " 

617 "nonzero corresponding constraint value.") 

618 complete = True 

619 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 

620 c0, x, undo, complete, status, message) 

621 else: # test_zero_row_2 

622 # if LHS is >= 0, we can eliminate this constraint entirely 

623 A_ub = A_ub[np.logical_not(zero_row), :] 

624 b_ub = b_ub[np.logical_not(zero_row)] 

625 

626 # zero column in (both) constraints 

627 # this indicates that a variable isn't constrained and can be removed 

628 A = vstack((A_eq, A_ub)) 

629 if A.shape[0] > 0: 

630 zero_col = np.array(np.sum(A != 0, axis=0) == 0).flatten() 

631 # variable will be at upper or lower bound, depending on objective 

632 x[np.logical_and(zero_col, c < 0)] = ub[ 

633 np.logical_and(zero_col, c < 0)] 

634 x[np.logical_and(zero_col, c > 0)] = lb[ 

635 np.logical_and(zero_col, c > 0)] 

636 if np.any(np.isinf(x)): # if an unconstrained variable has no bound 

637 status = 3 

638 message = ("If feasible, the problem is (trivially) unbounded " 

639 "due to a zero column in the constraint matrices. If " 

640 "you wish to check whether the problem is infeasible, " 

641 "turn presolve off.") 

642 complete = True 

643 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 

644 c0, x, undo, complete, status, message) 

645 # variables will equal upper/lower bounds will be removed later 

646 lb[np.logical_and(zero_col, c < 0)] = ub[ 

647 np.logical_and(zero_col, c < 0)] 

648 ub[np.logical_and(zero_col, c > 0)] = lb[ 

649 np.logical_and(zero_col, c > 0)] 

650 

651 # row singleton in equality constraints 

652 # this fixes a variable and removes the constraint 

653 singleton_row = np.array(np.sum(A_eq != 0, axis=1) == 1).flatten() 

654 rows = where(singleton_row)[0] 

655 cols = where(A_eq[rows, :])[1] 

656 if len(rows) > 0: 

657 for row, col in zip(rows, cols): 

658 val = b_eq[row] / A_eq[row, col] 

659 if not lb[col] - tol <= val <= ub[col] + tol: 

660 # infeasible if fixed value is not within bounds 

661 status = 2 

662 message = ("The problem is (trivially) infeasible because a " 

663 "singleton row in the equality constraints is " 

664 "inconsistent with the bounds.") 

665 complete = True 

666 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 

667 c0, x, undo, complete, status, message) 

668 else: 

669 # sets upper and lower bounds at that fixed value - variable 

670 # will be removed later 

671 lb[col] = val 

672 ub[col] = val 

673 A_eq = A_eq[np.logical_not(singleton_row), :] 

674 b_eq = b_eq[np.logical_not(singleton_row)] 

675 

676 # row singleton in inequality constraints 

677 # this indicates a simple bound and the constraint can be removed 

678 # simple bounds may be adjusted here 

679 # After all of the simple bound information is combined here, get_Abc will 

680 # turn the simple bounds into constraints 

681 singleton_row = np.array(np.sum(A_ub != 0, axis=1) == 1).flatten() 

682 cols = where(A_ub[singleton_row, :])[1] 

683 rows = where(singleton_row)[0] 

684 if len(rows) > 0: 

685 for row, col in zip(rows, cols): 

686 val = b_ub[row] / A_ub[row, col] 

687 if A_ub[row, col] > 0: # upper bound 

688 if val < lb[col] - tol: # infeasible 

689 complete = True 

690 elif val < ub[col]: # new upper bound 

691 ub[col] = val 

692 else: # lower bound 

693 if val > ub[col] + tol: # infeasible 

694 complete = True 

695 elif val > lb[col]: # new lower bound 

696 lb[col] = val 

697 if complete: 

698 status = 2 

699 message = ("The problem is (trivially) infeasible because a " 

700 "singleton row in the upper bound constraints is " 

701 "inconsistent with the bounds.") 

702 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 

703 c0, x, undo, complete, status, message) 

704 A_ub = A_ub[np.logical_not(singleton_row), :] 

705 b_ub = b_ub[np.logical_not(singleton_row)] 

706 

707 # identical bounds indicate that variable can be removed 

708 i_f = np.abs(lb - ub) < tol # indices of "fixed" variables 

709 i_nf = np.logical_not(i_f) # indices of "not fixed" variables 

710 

711 # test_bounds_equal_but_infeasible 

712 if np.all(i_f): # if bounds define solution, check for consistency 

713 residual = b_eq - A_eq.dot(lb) 

714 slack = b_ub - A_ub.dot(lb) 

715 if ((A_ub.size > 0 and np.any(slack < 0)) or 

716 (A_eq.size > 0 and not np.allclose(residual, 0))): 

717 status = 2 

718 message = ("The problem is (trivially) infeasible because the " 

719 "bounds fix all variables to values inconsistent with " 

720 "the constraints") 

721 complete = True 

722 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 

723 c0, x, undo, complete, status, message) 

724 

725 ub_mod = ub 

726 lb_mod = lb 

727 if np.any(i_f): 

728 c0 += c[i_f].dot(lb[i_f]) 

729 b_eq = b_eq - A_eq[:, i_f].dot(lb[i_f]) 

730 b_ub = b_ub - A_ub[:, i_f].dot(lb[i_f]) 

731 c = c[i_nf] 

732 x = x[i_nf] 

733 # user guess x0 stays separate from presolve solution x 

734 if x0 is not None: 

735 x0 = x0[i_nf] 

736 A_eq = A_eq[:, i_nf] 

737 A_ub = A_ub[:, i_nf] 

738 # record of variables to be added back in 

739 undo = [np.nonzero(i_f)[0], lb[i_f]] 

740 # don't remove these entries from bounds; they'll be used later. 

741 # but we _also_ need a version of the bounds with these removed 

742 lb_mod = lb[i_nf] 

743 ub_mod = ub[i_nf] 

744 

745 # no constraints indicates that problem is trivial 

746 if A_eq.size == 0 and A_ub.size == 0: 

747 b_eq = np.array([]) 

748 b_ub = np.array([]) 

749 # test_empty_constraint_1 

750 if c.size == 0: 

751 status = 0 

752 message = ("The solution was determined in presolve as there are " 

753 "no non-trivial constraints.") 

754 elif (np.any(np.logical_and(c < 0, ub_mod == np.inf)) or 

755 np.any(np.logical_and(c > 0, lb_mod == -np.inf))): 

756 # test_no_constraints() 

757 # test_unbounded_no_nontrivial_constraints_1 

758 # test_unbounded_no_nontrivial_constraints_2 

759 status = 3 

760 message = ("The problem is (trivially) unbounded " 

761 "because there are no non-trivial constraints and " 

762 "a) at least one decision variable is unbounded " 

763 "above and its corresponding cost is negative, or " 

764 "b) at least one decision variable is unbounded below " 

765 "and its corresponding cost is positive. ") 

766 else: # test_empty_constraint_2 

767 status = 0 

768 message = ("The solution was determined in presolve as there are " 

769 "no non-trivial constraints.") 

770 complete = True 

771 x[c < 0] = ub_mod[c < 0] 

772 x[c > 0] = lb_mod[c > 0] 

773 # where c is zero, set x to a finite bound or zero 

774 x_zero_c = ub_mod[c == 0] 

775 x_zero_c[np.isinf(x_zero_c)] = ub_mod[c == 0][np.isinf(x_zero_c)] 

776 x_zero_c[np.isinf(x_zero_c)] = 0 

777 x[c == 0] = x_zero_c 

778 # if this is not the last step of presolve, should convert bounds back 

779 # to array and return here 

780 

781 # Convert lb and ub back into Nx2 bounds 

782 bounds = np.hstack((lb[:, np.newaxis], ub[:, np.newaxis])) 

783 

784 # remove redundant (linearly dependent) rows from equality constraints 

785 n_rows_A = A_eq.shape[0] 

786 redundancy_warning = ("A_eq does not appear to be of full row rank. To " 

787 "improve performance, check the problem formulation " 

788 "for redundant equality constraints.") 

789 if (sps.issparse(A_eq)): 

790 if rr and A_eq.size > 0: # TODO: Fast sparse rank check? 

791 A_eq, b_eq, status, message = _remove_redundancy_sparse(A_eq, b_eq) 

792 if A_eq.shape[0] < n_rows_A: 

793 warn(redundancy_warning, OptimizeWarning, stacklevel=1) 

794 if status != 0: 

795 complete = True 

796 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 

797 c0, x, undo, complete, status, message) 

798 

799 # This is a wild guess for which redundancy removal algorithm will be 

800 # faster. More testing would be good. 

801 small_nullspace = 5 

802 if rr and A_eq.size > 0: 

803 try: # TODO: instead use results of first SVD in _remove_redundancy 

804 rank = np.linalg.matrix_rank(A_eq) 

805 except Exception: # oh well, we'll have to go with _remove_redundancy_dense 

806 rank = 0 

807 if rr and A_eq.size > 0 and rank < A_eq.shape[0]: 

808 warn(redundancy_warning, OptimizeWarning, stacklevel=3) 

809 dim_row_nullspace = A_eq.shape[0]-rank 

810 if dim_row_nullspace <= small_nullspace: 

811 A_eq, b_eq, status, message = _remove_redundancy(A_eq, b_eq) 

812 if dim_row_nullspace > small_nullspace or status == 4: 

813 A_eq, b_eq, status, message = _remove_redundancy_dense(A_eq, b_eq) 

814 if A_eq.shape[0] < rank: 

815 message = ("Due to numerical issues, redundant equality " 

816 "constraints could not be removed automatically. " 

817 "Try providing your constraint matrices as sparse " 

818 "matrices to activate sparse presolve, try turning " 

819 "off redundancy removal, or try turning off presolve " 

820 "altogether.") 

821 status = 4 

822 if status != 0: 

823 complete = True 

824 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 

825 c0, x, undo, complete, status, message) 

826 

827 

828def _parse_linprog(lp, options): 

829 """ 

830 Parse the provided linear programming problem 

831 

832 ``_parse_linprog`` employs two main steps ``_check_sparse_inputs`` and 

833 ``_clean_inputs``. ``_check_sparse_inputs`` checks for sparsity in the 

834 provided constraints (``A_ub`` and ``A_eq) and if these match the provided 

835 sparsity optional values. 

836 

837 ``_clean inputs`` checks of the provided inputs. If no violations are 

838 identified the objective vector, upper bound constraints, equality 

839 constraints, and simple bounds are returned in the expected format. 

840 

841 Parameters 

842 ---------- 

843 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 

844 

845 c : 1D array 

846 The coefficients of the linear objective function to be minimized. 

847 A_ub : 2D array, optional 

848 The inequality constraint matrix. Each row of ``A_ub`` specifies the 

849 coefficients of a linear inequality constraint on ``x``. 

850 b_ub : 1D array, optional 

851 The inequality constraint vector. Each element represents an 

852 upper bound on the corresponding value of ``A_ub @ x``. 

853 A_eq : 2D array, optional 

854 The equality constraint matrix. Each row of ``A_eq`` specifies the 

855 coefficients of a linear equality constraint on ``x``. 

856 b_eq : 1D array, optional 

857 The equality constraint vector. Each element of ``A_eq @ x`` must equal 

858 the corresponding element of ``b_eq``. 

859 bounds : various valid formats, optional 

860 The bounds of ``x``, as ``min`` and ``max`` pairs. 

861 If bounds are specified for all N variables separately, valid formats are: 

862 * a 2D array (2 x N or N x 2); 

863 * a sequence of N sequences, each with 2 values. 

864 If all variables have the same bounds, a single pair of values can 

865 be specified. Valid formats are: 

866 * a sequence with 2 scalar values; 

867 * a sequence with a single element containing 2 scalar values. 

868 If all variables have a lower bound of 0 and no upper bound, the bounds 

869 parameter can be omitted (or given as None). 

870 x0 : 1D array, optional 

871 Guess values of the decision variables, which will be refined by 

872 the optimization algorithm. This argument is currently used only by the 

873 'revised simplex' method, and can only be used if `x0` represents a 

874 basic feasible solution. 

875 

876 options : dict 

877 A dictionary of solver options. All methods accept the following 

878 generic options: 

879 

880 maxiter : int 

881 Maximum number of iterations to perform. 

882 disp : bool 

883 Set to True to print convergence messages. 

884 

885 For method-specific options, see :func:`show_options('linprog')`. 

886 

887 Returns 

888 ------- 

889 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 

890 

891 c : 1D array 

892 The coefficients of the linear objective function to be minimized. 

893 A_ub : 2D array, optional 

894 The inequality constraint matrix. Each row of ``A_ub`` specifies the 

895 coefficients of a linear inequality constraint on ``x``. 

896 b_ub : 1D array, optional 

897 The inequality constraint vector. Each element represents an 

898 upper bound on the corresponding value of ``A_ub @ x``. 

899 A_eq : 2D array, optional 

900 The equality constraint matrix. Each row of ``A_eq`` specifies the 

901 coefficients of a linear equality constraint on ``x``. 

902 b_eq : 1D array, optional 

903 The equality constraint vector. Each element of ``A_eq @ x`` must equal 

904 the corresponding element of ``b_eq``. 

905 bounds : 2D array 

906 The bounds of ``x``, as ``min`` and ``max`` pairs, one for each of the N 

907 elements of ``x``. The N x 2 array contains lower bounds in the first 

908 column and upper bounds in the 2nd. Unbounded variables have lower 

909 bound -np.inf and/or upper bound np.inf. 

910 x0 : 1D array, optional 

911 Guess values of the decision variables, which will be refined by 

912 the optimization algorithm. This argument is currently used only by the 

913 'revised simplex' method, and can only be used if `x0` represents a 

914 basic feasible solution. 

915 

916 options : dict, optional 

917 A dictionary of solver options. All methods accept the following 

918 generic options: 

919 

920 maxiter : int 

921 Maximum number of iterations to perform. 

922 disp : bool 

923 Set to True to print convergence messages. 

924 

925 For method-specific options, see :func:`show_options('linprog')`. 

926 

927 """ 

928 if options is None: 

929 options = {} 

930 

931 solver_options = {k: v for k, v in options.items()} 

932 solver_options, A_ub, A_eq = _check_sparse_inputs(solver_options, lp.A_ub, lp.A_eq) 

933 # Convert lists to numpy arrays, etc... 

934 lp = _clean_inputs(lp._replace(A_ub=A_ub, A_eq=A_eq)) 

935 return lp, solver_options 

936 

937 

938def _get_Abc(lp, c0, undo=[]): 

939 """ 

940 Given a linear programming problem of the form: 

941 

942 Minimize:: 

943 

944 c @ x 

945 

946 Subject to:: 

947 

948 A_ub @ x <= b_ub 

949 A_eq @ x == b_eq 

950 lb <= x <= ub 

951 

952 where ``lb = 0`` and ``ub = None`` unless set in ``bounds``. 

953 

954 Return the problem in standard form: 

955 

956 Minimize:: 

957 

958 c @ x 

959 

960 Subject to:: 

961 

962 A @ x == b 

963 x >= 0 

964 

965 by adding slack variables and making variable substitutions as necessary. 

966 

967 Parameters 

968 ---------- 

969 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 

970 

971 c : 1D array 

972 The coefficients of the linear objective function to be minimized. 

973 A_ub : 2D array, optional 

974 The inequality constraint matrix. Each row of ``A_ub`` specifies the 

975 coefficients of a linear inequality constraint on ``x``. 

976 b_ub : 1D array, optional 

977 The inequality constraint vector. Each element represents an 

978 upper bound on the corresponding value of ``A_ub @ x``. 

979 A_eq : 2D array, optional 

980 The equality constraint matrix. Each row of ``A_eq`` specifies the 

981 coefficients of a linear equality constraint on ``x``. 

982 b_eq : 1D array, optional 

983 The equality constraint vector. Each element of ``A_eq @ x`` must equal 

984 the corresponding element of ``b_eq``. 

985 bounds : 2D array 

986 The bounds of ``x``, lower bounds in the 1st column, upper 

987 bounds in the 2nd column. The bounds are possibly tightened 

988 by the presolve procedure. 

989 x0 : 1D array, optional 

990 Guess values of the decision variables, which will be refined by 

991 the optimization algorithm. This argument is currently used only by the 

992 'revised simplex' method, and can only be used if `x0` represents a 

993 basic feasible solution. 

994 

995 c0 : float 

996 Constant term in objective function due to fixed (and eliminated) 

997 variables. 

998 

999 undo: list of tuples 

1000 (`index`, `value`) pairs that record the original index and fixed value 

1001 for each variable removed from the problem 

1002 

1003 Returns 

1004 ------- 

1005 A : 2-D array 

1006 2-D array such that ``A`` @ ``x``, gives the values of the equality 

1007 constraints at ``x``. 

1008 b : 1-D array 

1009 1-D array of values representing the RHS of each equality constraint 

1010 (row) in A (for standard form problem). 

1011 c : 1-D array 

1012 Coefficients of the linear objective function to be minimized (for 

1013 standard form problem). 

1014 c0 : float 

1015 Constant term in objective function due to fixed (and eliminated) 

1016 variables. 

1017 x0 : 1-D array 

1018 Starting values of the independent variables, which will be refined by 

1019 the optimization algorithm 

1020 

1021 References 

1022 ---------- 

1023 .. [9] Bertsimas, Dimitris, and J. Tsitsiklis. "Introduction to linear 

1024 programming." Athena Scientific 1 (1997): 997. 

1025 

1026 """ 

1027 c, A_ub, b_ub, A_eq, b_eq, bounds, x0 = lp 

1028 

1029 if sps.issparse(A_eq): 

1030 sparse = True 

1031 A_eq = sps.csr_matrix(A_eq) 

1032 A_ub = sps.csr_matrix(A_ub) 

1033 

1034 def hstack(blocks): 

1035 return sps.hstack(blocks, format="csr") 

1036 

1037 def vstack(blocks): 

1038 return sps.vstack(blocks, format="csr") 

1039 

1040 zeros = sps.csr_matrix 

1041 eye = sps.eye 

1042 else: 

1043 sparse = False 

1044 hstack = np.hstack 

1045 vstack = np.vstack 

1046 zeros = np.zeros 

1047 eye = np.eye 

1048 

1049 # bounds will be modified, create a copy 

1050 bounds = np.array(bounds, copy=True) 

1051 # undo[0] contains indices of variables removed from the problem 

1052 # however, their bounds are still part of the bounds list 

1053 # they are needed elsewhere, but not here 

1054 if undo is not None and undo != []: 

1055 bounds = np.delete(bounds, undo[0], 0) 

1056 

1057 # modify problem such that all variables have only non-negativity bounds 

1058 lbs = bounds[:, 0] 

1059 ubs = bounds[:, 1] 

1060 m_ub, n_ub = A_ub.shape 

1061 

1062 lb_none = np.equal(lbs, -np.inf) 

1063 ub_none = np.equal(ubs, np.inf) 

1064 lb_some = np.logical_not(lb_none) 

1065 ub_some = np.logical_not(ub_none) 

1066 

1067 # if preprocessing is on, lb == ub can't happen 

1068 # if preprocessing is off, then it would be best to convert that 

1069 # to an equality constraint, but it's tricky to make the other 

1070 # required modifications from inside here. 

1071 

1072 # unbounded below: substitute xi = -xi' (unbounded above) 

1073 # if -inf <= xi <= ub, then -ub <= -xi <= inf, so swap and invert bounds 

1074 l_nolb_someub = np.logical_and(lb_none, ub_some) 

1075 i_nolb = np.nonzero(l_nolb_someub)[0] 

1076 lbs[l_nolb_someub], ubs[l_nolb_someub] = ( 

1077 -ubs[l_nolb_someub], -lbs[l_nolb_someub]) 

1078 lb_none = np.equal(lbs, -np.inf) 

1079 ub_none = np.equal(ubs, np.inf) 

1080 lb_some = np.logical_not(lb_none) 

1081 ub_some = np.logical_not(ub_none) 

1082 c[i_nolb] *= -1 

1083 if x0 is not None: 

1084 x0[i_nolb] *= -1 

1085 if len(i_nolb) > 0: 

1086 if A_ub.shape[0] > 0: # sometimes needed for sparse arrays... weird 

1087 A_ub[:, i_nolb] *= -1 

1088 if A_eq.shape[0] > 0: 

1089 A_eq[:, i_nolb] *= -1 

1090 

1091 # upper bound: add inequality constraint 

1092 i_newub, = ub_some.nonzero() 

1093 ub_newub = ubs[ub_some] 

1094 n_bounds = len(i_newub) 

1095 if n_bounds > 0: 

1096 shape = (n_bounds, A_ub.shape[1]) 

1097 if sparse: 

1098 idxs = (np.arange(n_bounds), i_newub) 

1099 A_ub = vstack((A_ub, sps.csr_matrix((np.ones(n_bounds), idxs), 

1100 shape=shape))) 

1101 else: 

1102 A_ub = vstack((A_ub, np.zeros(shape))) 

1103 A_ub[np.arange(m_ub, A_ub.shape[0]), i_newub] = 1 

1104 b_ub = np.concatenate((b_ub, np.zeros(n_bounds))) 

1105 b_ub[m_ub:] = ub_newub 

1106 

1107 A1 = vstack((A_ub, A_eq)) 

1108 b = np.concatenate((b_ub, b_eq)) 

1109 c = np.concatenate((c, np.zeros((A_ub.shape[0],)))) 

1110 if x0 is not None: 

1111 x0 = np.concatenate((x0, np.zeros((A_ub.shape[0],)))) 

1112 # unbounded: substitute xi = xi+ + xi- 

1113 l_free = np.logical_and(lb_none, ub_none) 

1114 i_free = np.nonzero(l_free)[0] 

1115 n_free = len(i_free) 

1116 c = np.concatenate((c, np.zeros(n_free))) 

1117 if x0 is not None: 

1118 x0 = np.concatenate((x0, np.zeros(n_free))) 

1119 A1 = hstack((A1[:, :n_ub], -A1[:, i_free])) 

1120 c[n_ub:n_ub+n_free] = -c[i_free] 

1121 if x0 is not None: 

1122 i_free_neg = x0[i_free] < 0 

1123 x0[np.arange(n_ub, A1.shape[1])[i_free_neg]] = -x0[i_free[i_free_neg]] 

1124 x0[i_free[i_free_neg]] = 0 

1125 

1126 # add slack variables 

1127 A2 = vstack([eye(A_ub.shape[0]), zeros((A_eq.shape[0], A_ub.shape[0]))]) 

1128 

1129 A = hstack([A1, A2]) 

1130 

1131 # lower bound: substitute xi = xi' + lb 

1132 # now there is a constant term in objective 

1133 i_shift = np.nonzero(lb_some)[0] 

1134 lb_shift = lbs[lb_some].astype(float) 

1135 c0 += np.sum(lb_shift * c[i_shift]) 

1136 if sparse: 

1137 b = b.reshape(-1, 1) 

1138 A = A.tocsc() 

1139 b -= (A[:, i_shift] * sps.diags(lb_shift)).sum(axis=1) 

1140 b = b.ravel() 

1141 else: 

1142 b -= (A[:, i_shift] * lb_shift).sum(axis=1) 

1143 if x0 is not None: 

1144 x0[i_shift] -= lb_shift 

1145 

1146 return A, b, c, c0, x0 

1147 

1148 

1149def _round_to_power_of_two(x): 

1150 """ 

1151 Round elements of the array to the nearest power of two. 

1152 """ 

1153 return 2**np.around(np.log2(x)) 

1154 

1155 

1156def _autoscale(A, b, c, x0): 

1157 """ 

1158 Scales the problem according to equilibration from [12]. 

1159 Also normalizes the right hand side vector by its maximum element. 

1160 """ 

1161 m, n = A.shape 

1162 

1163 C = 1 

1164 R = 1 

1165 

1166 if A.size > 0: 

1167 

1168 R = np.max(np.abs(A), axis=1) 

1169 if sps.issparse(A): 

1170 R = R.toarray().flatten() 

1171 R[R == 0] = 1 

1172 R = 1/_round_to_power_of_two(R) 

1173 A = sps.diags(R)*A if sps.issparse(A) else A*R.reshape(m, 1) 

1174 b = b*R 

1175 

1176 C = np.max(np.abs(A), axis=0) 

1177 if sps.issparse(A): 

1178 C = C.toarray().flatten() 

1179 C[C == 0] = 1 

1180 C = 1/_round_to_power_of_two(C) 

1181 A = A*sps.diags(C) if sps.issparse(A) else A*C 

1182 c = c*C 

1183 

1184 b_scale = np.max(np.abs(b)) if b.size > 0 else 1 

1185 if b_scale == 0: 

1186 b_scale = 1. 

1187 b = b/b_scale 

1188 

1189 if x0 is not None: 

1190 x0 = x0/b_scale*(1/C) 

1191 return A, b, c, x0, C, b_scale 

1192 

1193 

1194def _unscale(x, C, b_scale): 

1195 """ 

1196 Converts solution to _autoscale problem -> solution to original problem. 

1197 """ 

1198 

1199 try: 

1200 n = len(C) 

1201 # fails if sparse or scalar; that's OK. 

1202 # this is only needed for original simplex (never sparse) 

1203 except TypeError: 

1204 n = len(x) 

1205 

1206 return x[:n]*b_scale*C 

1207 

1208 

1209def _display_summary(message, status, fun, iteration): 

1210 """ 

1211 Print the termination summary of the linear program 

1212 

1213 Parameters 

1214 ---------- 

1215 message : str 

1216 A string descriptor of the exit status of the optimization. 

1217 status : int 

1218 An integer representing the exit status of the optimization:: 

1219 

1220 0 : Optimization terminated successfully 

1221 1 : Iteration limit reached 

1222 2 : Problem appears to be infeasible 

1223 3 : Problem appears to be unbounded 

1224 4 : Serious numerical difficulties encountered 

1225 

1226 fun : float 

1227 Value of the objective function. 

1228 iteration : iteration 

1229 The number of iterations performed. 

1230 """ 

1231 print(message) 

1232 if status in (0, 1): 

1233 print(" Current function value: {0: <12.6f}".format(fun)) 

1234 print(" Iterations: {0:d}".format(iteration)) 

1235 

1236 

1237def _postsolve(x, postsolve_args, complete=False, tol=1e-8, copy=False): 

1238 """ 

1239 Given solution x to presolved, standard form linear program x, add 

1240 fixed variables back into the problem and undo the variable substitutions 

1241 to get solution to original linear program. Also, calculate the objective 

1242 function value, slack in original upper bound constraints, and residuals 

1243 in original equality constraints. 

1244 

1245 Parameters 

1246 ---------- 

1247 x : 1-D array 

1248 Solution vector to the standard-form problem. 

1249 postsolve_args : tuple 

1250 Data needed by _postsolve to convert the solution to the standard-form 

1251 problem into the solution to the original problem, including: 

1252 

1253 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 

1254 

1255 c : 1D array 

1256 The coefficients of the linear objective function to be minimized. 

1257 A_ub : 2D array, optional 

1258 The inequality constraint matrix. Each row of ``A_ub`` specifies the 

1259 coefficients of a linear inequality constraint on ``x``. 

1260 b_ub : 1D array, optional 

1261 The inequality constraint vector. Each element represents an 

1262 upper bound on the corresponding value of ``A_ub @ x``. 

1263 A_eq : 2D array, optional 

1264 The equality constraint matrix. Each row of ``A_eq`` specifies the 

1265 coefficients of a linear equality constraint on ``x``. 

1266 b_eq : 1D array, optional 

1267 The equality constraint vector. Each element of ``A_eq @ x`` must equal 

1268 the corresponding element of ``b_eq``. 

1269 bounds : 2D array 

1270 The bounds of ``x``, lower bounds in the 1st column, upper 

1271 bounds in the 2nd column. The bounds are possibly tightened 

1272 by the presolve procedure. 

1273 x0 : 1D array, optional 

1274 Guess values of the decision variables, which will be refined by 

1275 the optimization algorithm. This argument is currently used only by the 

1276 'revised simplex' method, and can only be used if `x0` represents a 

1277 basic feasible solution. 

1278 

1279 undo: list of tuples 

1280 (`index`, `value`) pairs that record the original index and fixed value 

1281 for each variable removed from the problem 

1282 complete : bool 

1283 Whether the solution is was determined in presolve (``True`` if so) 

1284 tol : float 

1285 Termination tolerance; see [1]_ Section 4.5. 

1286 

1287 Returns 

1288 ------- 

1289 x : 1-D array 

1290 Solution vector to original linear programming problem 

1291 fun: float 

1292 optimal objective value for original problem 

1293 slack : 1-D array 

1294 The (non-negative) slack in the upper bound constraints, that is, 

1295 ``b_ub - A_ub @ x`` 

1296 con : 1-D array 

1297 The (nominally zero) residuals of the equality constraints, that is, 

1298 ``b - A_eq @ x`` 

1299 bounds : 2D array 

1300 The bounds on the original variables ``x`` 

1301 """ 

1302 # note that all the inputs are the ORIGINAL, unmodified versions 

1303 # no rows, columns have been removed 

1304 # the only exception is bounds; it has been modified 

1305 # we need these modified values to undo the variable substitutions 

1306 # in retrospect, perhaps this could have been simplified if the "undo" 

1307 # variable also contained information for undoing variable substitutions 

1308 

1309 (c, A_ub, b_ub, A_eq, b_eq, bounds, x0), undo, C, b_scale = postsolve_args 

1310 x = _unscale(x, C, b_scale) 

1311 

1312 n_x = len(c) 

1313 

1314 # we don't have to undo variable substitutions for fixed variables that 

1315 # were removed from the problem 

1316 no_adjust = set() 

1317 

1318 # if there were variables removed from the problem, add them back into the 

1319 # solution vector 

1320 if len(undo) > 0: 

1321 no_adjust = set(undo[0]) 

1322 x = x.tolist() 

1323 for i, val in zip(undo[0], undo[1]): 

1324 x.insert(i, val) 

1325 copy = True 

1326 if copy: 

1327 x = np.array(x, copy=True) 

1328 

1329 # now undo variable substitutions 

1330 # if "complete", problem was solved in presolve; don't do anything here 

1331 if not complete and bounds is not None: # bounds are never none, probably 

1332 n_unbounded = 0 

1333 for i, bi in enumerate(bounds): 

1334 if i in no_adjust: 

1335 continue 

1336 lbi = bi[0] 

1337 ubi = bi[1] 

1338 if lbi == -np.inf and ubi == np.inf: 

1339 n_unbounded += 1 

1340 x[i] = x[i] - x[n_x + n_unbounded - 1] 

1341 else: 

1342 if lbi == -np.inf: 

1343 x[i] = ubi - x[i] 

1344 else: 

1345 x[i] += lbi 

1346 

1347 n_x = len(c) 

1348 x = x[:n_x] # all the rest of the variables were artificial 

1349 fun = x.dot(c) 

1350 slack = b_ub - A_ub.dot(x) # report slack for ORIGINAL UB constraints 

1351 # report residuals of ORIGINAL EQ constraints 

1352 con = b_eq - A_eq.dot(x) 

1353 

1354 return x, fun, slack, con, bounds 

1355 

1356 

1357def _check_result(x, fun, status, slack, con, bounds, tol, message): 

1358 """ 

1359 Check the validity of the provided solution. 

1360 

1361 A valid (optimal) solution satisfies all bounds, all slack variables are 

1362 negative and all equality constraint residuals are strictly non-zero. 

1363 Further, the lower-bounds, upper-bounds, slack and residuals contain 

1364 no nan values. 

1365 

1366 Parameters 

1367 ---------- 

1368 x : 1-D array 

1369 Solution vector to original linear programming problem 

1370 fun: float 

1371 optimal objective value for original problem 

1372 status : int 

1373 An integer representing the exit status of the optimization:: 

1374 

1375 0 : Optimization terminated successfully 

1376 1 : Iteration limit reached 

1377 2 : Problem appears to be infeasible 

1378 3 : Problem appears to be unbounded 

1379 4 : Serious numerical difficulties encountered 

1380 

1381 slack : 1-D array 

1382 The (non-negative) slack in the upper bound constraints, that is, 

1383 ``b_ub - A_ub @ x`` 

1384 con : 1-D array 

1385 The (nominally zero) residuals of the equality constraints, that is, 

1386 ``b - A_eq @ x`` 

1387 bounds : 2D array 

1388 The bounds on the original variables ``x`` 

1389 message : str 

1390 A string descriptor of the exit status of the optimization. 

1391 tol : float 

1392 Termination tolerance; see [1]_ Section 4.5. 

1393 

1394 Returns 

1395 ------- 

1396 status : int 

1397 An integer representing the exit status of the optimization:: 

1398 

1399 0 : Optimization terminated successfully 

1400 1 : Iteration limit reached 

1401 2 : Problem appears to be infeasible 

1402 3 : Problem appears to be unbounded 

1403 4 : Serious numerical difficulties encountered 

1404 

1405 message : str 

1406 A string descriptor of the exit status of the optimization. 

1407 """ 

1408 # Somewhat arbitrary, but status 5 is very unusual 

1409 tol = np.sqrt(tol) * 10 

1410 

1411 contains_nans = ( 

1412 np.isnan(x).any() 

1413 or np.isnan(fun) 

1414 or np.isnan(slack).any() 

1415 or np.isnan(con).any() 

1416 ) 

1417 

1418 if contains_nans: 

1419 is_feasible = False 

1420 else: 

1421 invalid_bounds = (x < bounds[:, 0] - tol).any() or (x > bounds[:, 1] + tol).any() 

1422 invalid_slack = status != 3 and (slack < -tol).any() 

1423 invalid_con = status != 3 and (np.abs(con) > tol).any() 

1424 is_feasible = not (invalid_bounds or invalid_slack or invalid_con) 

1425 

1426 if status == 0 and not is_feasible: 

1427 status = 4 

1428 message = ("The solution does not satisfy the constraints within the " 

1429 "required tolerance of " + "{:.2E}".format(tol) + ", yet " 

1430 "no errors were raised and there is no certificate of " 

1431 "infeasibility or unboundedness. This is known to occur " 

1432 "if the `presolve` option is False and the problem is " 

1433 "infeasible. This can also occur due to the limited " 

1434 "accuracy of the `interior-point` method. Check whether " 

1435 "the slack and constraint residuals are acceptable; " 

1436 "if not, consider enabling presolve, reducing option " 

1437 "`tol`, and/or using method `revised simplex`. " 

1438 "If you encounter this message under different " 

1439 "circumstances, please submit a bug report.") 

1440 elif status == 0 and contains_nans: 

1441 status = 4 

1442 message = ("Numerical difficulties were encountered but no errors " 

1443 "were raised. This is known to occur if the 'presolve' " 

1444 "option is False, 'sparse' is True, and A_eq includes " 

1445 "redundant rows. If you encounter this under different " 

1446 "circumstances, please submit a bug report. Otherwise, " 

1447 "remove linearly dependent equations from your equality " 

1448 "constraints or enable presolve.") 

1449 elif status == 2 and is_feasible: 

1450 # Occurs if the simplex method exits after phase one with a very 

1451 # nearly basic feasible solution. Postsolving can make the solution 

1452 # basic, however, this solution is NOT optimal 

1453 raise ValueError(message) 

1454 

1455 return status, message 

1456 

1457 

1458def _postprocess(x, postsolve_args, complete=False, status=0, message="", 

1459 tol=1e-8, iteration=None, disp=False): 

1460 """ 

1461 Given solution x to presolved, standard form linear program x, add 

1462 fixed variables back into the problem and undo the variable substitutions 

1463 to get solution to original linear program. Also, calculate the objective 

1464 function value, slack in original upper bound constraints, and residuals 

1465 in original equality constraints. 

1466 

1467 Parameters 

1468 ---------- 

1469 x : 1-D array 

1470 Solution vector to the standard-form problem. 

1471 c : 1-D array 

1472 Original coefficients of the linear objective function to be minimized. 

1473 A_ub : 2-D array, optional 

1474 2-D array such that ``A_ub @ x`` gives the values of the upper-bound 

1475 inequality constraints at ``x``. 

1476 b_ub : 1-D array, optional 

1477 1-D array of values representing the upper-bound of each inequality 

1478 constraint (row) in ``A_ub``. 

1479 A_eq : 2-D array, optional 

1480 2-D array such that ``A_eq @ x`` gives the values of the equality 

1481 constraints at ``x``. 

1482 b_eq : 1-D array, optional 

1483 1-D array of values representing the RHS of each equality constraint 

1484 (row) in ``A_eq``. 

1485 bounds : 2D array 

1486 The bounds of ``x``, lower bounds in the 1st column, upper 

1487 bounds in the 2nd column. The bounds are possibly tightened 

1488 by the presolve procedure. 

1489 complete : bool 

1490 Whether the solution is was determined in presolve (``True`` if so) 

1491 undo: list of tuples 

1492 (`index`, `value`) pairs that record the original index and fixed value 

1493 for each variable removed from the problem 

1494 status : int 

1495 An integer representing the exit status of the optimization:: 

1496 

1497 0 : Optimization terminated successfully 

1498 1 : Iteration limit reached 

1499 2 : Problem appears to be infeasible 

1500 3 : Problem appears to be unbounded 

1501 4 : Serious numerical difficulties encountered 

1502 

1503 message : str 

1504 A string descriptor of the exit status of the optimization. 

1505 tol : float 

1506 Termination tolerance; see [1]_ Section 4.5. 

1507 

1508 Returns 

1509 ------- 

1510 x : 1-D array 

1511 Solution vector to original linear programming problem 

1512 fun: float 

1513 optimal objective value for original problem 

1514 slack : 1-D array 

1515 The (non-negative) slack in the upper bound constraints, that is, 

1516 ``b_ub - A_ub @ x`` 

1517 con : 1-D array 

1518 The (nominally zero) residuals of the equality constraints, that is, 

1519 ``b - A_eq @ x`` 

1520 status : int 

1521 An integer representing the exit status of the optimization:: 

1522 

1523 0 : Optimization terminated successfully 

1524 1 : Iteration limit reached 

1525 2 : Problem appears to be infeasible 

1526 3 : Problem appears to be unbounded 

1527 4 : Serious numerical difficulties encountered 

1528 

1529 message : str 

1530 A string descriptor of the exit status of the optimization. 

1531 

1532 """ 

1533 

1534 x, fun, slack, con, bounds = _postsolve( 

1535 x, postsolve_args, complete, tol 

1536 ) 

1537 

1538 status, message = _check_result( 

1539 x, fun, status, slack, con, 

1540 bounds, tol, message 

1541 ) 

1542 

1543 if disp: 

1544 _display_summary(message, status, fun, iteration) 

1545 

1546 return x, fun, slack, con, status, message