Coverage for picos/solvers/solver_glpk.py: 79.37%
378 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
1# ------------------------------------------------------------------------------
2# Copyright (C) 2017-2019 Maximilian Stahlberg
3#
4# This file is part of PICOS.
5#
6# PICOS is free software: you can redistribute it and/or modify it under the
7# terms of the GNU General Public License as published by the Free Software
8# Foundation, either version 3 of the License, or (at your option) any later
9# version.
10#
11# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
12# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License along with
16# this program. If not, see <http://www.gnu.org/licenses/>.
17# ------------------------------------------------------------------------------
19"""Implementation of :class:`GLPKSolver`."""
21from collections import namedtuple
23import cvxopt
25from ..apidoc import api_end, api_start
26from ..constraints import AffineConstraint, DummyConstraint
27from ..expressions import AffineExpression, BinaryVariable, IntegerVariable
28from ..modeling.footprint import Specification
29from ..modeling.solution import (PS_FEASIBLE, PS_INFEASIBLE, PS_UNBOUNDED,
30 PS_UNKNOWN, SS_EMPTY, SS_FEASIBLE,
31 SS_INFEASIBLE, SS_OPTIMAL, SS_UNKNOWN)
32from .solver import Solver
34_API_START = api_start(globals())
35# -------------------------------
38class GLPKSolver(Solver):
39 """Interface to the GLPK solver via swiglpk."""
41 SUPPORTED = Specification(
42 objectives=[
43 AffineExpression],
44 constraints=[
45 DummyConstraint,
46 AffineConstraint])
48 UNUSED_VAR = namedtuple("UnusedVar", ("dim",))(dim=1)
50 @classmethod
51 def supports(cls, footprint, explain=False):
52 """Implement :meth:`~.solver.Solver.supports`."""
53 result = Solver.supports(footprint, explain)
54 if not result or (explain and not result[0]):
55 return result
57 if footprint not in cls.SUPPORTED:
58 if explain:
59 return False, cls.SUPPORTED.mismatch_reason(footprint)
60 else:
61 return False
63 return (True, None) if explain else True
65 @classmethod
66 def default_penalty(cls):
67 """Implement :meth:`~.solver.Solver.default_penalty`."""
68 return 1.0 # Stable free/open source solver.
70 @classmethod
71 def test_availability(cls):
72 """Implement :meth:`~.solver.Solver.test_availability`."""
73 cls.check_import("swiglpk")
75 @classmethod
76 def names(cls):
77 """Implement :meth:`~.solver.Solver.names`."""
78 return "glpk", "GLPK", "GNU Linear Programming Kit", "swiglpk"
80 @classmethod
81 def is_free(cls):
82 """Implement :meth:`~.solver.Solver.is_free`."""
83 return True
85 def __init__(self, problem):
86 """Initialize a GLPK solver interface.
88 :param ~picos.Problem problem: The problem to be solved.
89 """
90 self._glpkVarOffset = {}
91 self._glpkConOffset = {}
93 super(GLPKSolver, self).__init__(problem)
95 def __del__(self):
96 try:
97 import swiglpk as glpk
98 except ImportError:
99 # Happens when python is shutting down. In this case, no garbage
100 # collection is necessary, anyway.
101 return
103 if self.int is not None:
104 glpk.glp_delete_prob(self.int)
106 def reset_problem(self):
107 """Implement :meth:`~.solver.Solver.reset_problem`."""
108 import swiglpk as glpk
110 self._glpkVarOffset.clear()
111 self._glpkConOffset.clear()
113 if self.int is not None:
114 glpk.glp_delete_prob(self.int)
115 self.int = None
117 def _import_variable(self, picosVar):
118 import swiglpk as glpk
120 p = self.int
121 dim = picosVar.dim
123 # Append columns to the constraint matrix.
124 offset = glpk.glp_add_cols(p, dim)
126 # Retrieve the type.
127 if isinstance(picosVar, IntegerVariable):
128 varType = glpk.GLP_IV
129 elif isinstance(picosVar, BinaryVariable):
130 varType = glpk.GLP_BV
131 else:
132 varType = glpk.GLP_CV
134 # Retrieve bounds.
135 if not isinstance(picosVar, BinaryVariable):
136 boundKeys = [glpk.GLP_FR]*dim
137 lowerBounds = [0]*dim
138 upperBounds = [0]*dim
140 lower, upper = picosVar.bound_dicts
142 for i, b in lower.items():
143 boundKeys[i] = glpk.GLP_LO
144 lowerBounds[i] = b
146 for i, b in upper.items():
147 if boundKeys[i] == glpk.GLP_FR:
148 boundKeys[i] = glpk.GLP_UP
149 else: # Also has a lower bound.
150 if lowerBounds[i] == b:
151 boundKeys[i] = glpk.GLP_FX
152 else:
153 boundKeys[i] = glpk.GLP_DB
155 upperBounds[i] = b
157 # Import scalar variables.
158 for i in range(picosVar.dim):
159 glpkIndex = offset + i
161 # Assign a name.
162 glpk.glp_set_col_name(
163 p, glpkIndex, "{}_{}".format(picosVar.name, i))
165 # Set the type.
166 glpk.glp_set_col_kind(p, glpkIndex, varType)
168 # Assign the bounds.
169 if not isinstance(picosVar, BinaryVariable):
170 glpk.glp_set_col_bnds(
171 p, glpkIndex, boundKeys[i], lowerBounds[i], upperBounds[i])
173 self._glpkVarOffset[picosVar] = offset
175 def _remove_variable(self, picosVar):
176 import swiglpk as glpk
178 dim = picosVar.dim
180 offset = self._glpkVarOffset.pop(picosVar)
182 glpkVars = glpk.intArray(dim + 1)
183 for i in range(dim):
184 glpkVars[i + 1] = offset + i # Index 0 is unused.
186 glpk.glp_del_cols(self.int, dim, glpkVars)
188 for otherVar, otherOffset in self._glpkVarOffset.items():
189 if otherOffset > offset:
190 self._glpkVarOffset[otherVar] -= dim
192 def _import_constraint(self, picosCon):
193 import swiglpk as glpk
195 p = self.int
197 # Append rows to the constraint matrix.
198 rowOffset = glpk.glp_add_rows(p, len(picosCon))
200 # Import scalar constraints.
201 lmr_rows = picosCon.lmr.sparse_rows(self._glpkVarOffset)
202 for localConIndex, (glpkVarIndices, coefs, c) in enumerate(lmr_rows):
203 rhs = -c
204 glpkConIndex = rowOffset + localConIndex
205 numColumns = len(glpkVarIndices)
207 # Assign a name.
208 glpk.glp_set_row_name(
209 p, glpkConIndex, "{}_{}".format(picosCon.id, localConIndex))
211 # Set the constant term.
212 if picosCon.is_equality():
213 glpk.glp_set_row_bnds(p, glpkConIndex, glpk.GLP_FX, rhs, rhs)
214 elif picosCon.is_increasing():
215 glpk.glp_set_row_bnds(p, glpkConIndex, glpk.GLP_UP, 0, rhs)
216 elif picosCon.is_decreasing():
217 glpk.glp_set_row_bnds(p, glpkConIndex, glpk.GLP_LO, rhs, 0)
218 else:
219 assert False, "Unexpected constraint relation."
221 # Set coefficients.
222 # NOTE: GLPK requires a glpk.intArray containing column indices and
223 # a glpk.doubleArray of same size containing the coefficients
224 # for the listed column index. The first element of both
225 # arrays (with index 0) is skipped by GLPK.
226 glpkVarIndicesArray = glpk.intArray(numColumns + 1)
227 for i in range(numColumns):
228 glpkVarIndicesArray[i + 1] = glpkVarIndices[i]
230 coefficientsArray = glpk.doubleArray(numColumns + 1)
231 for i in range(numColumns):
232 coefficientsArray[i + 1] = coefs[i]
234 glpk.glp_set_mat_row(p, glpkConIndex, numColumns,
235 glpkVarIndicesArray, coefficientsArray)
237 self._glpkConOffset[picosCon] = rowOffset
239 def _remove_constraint(self, picosCon):
240 import swiglpk as glpk
242 length = len(picosCon)
244 offset = self._glpkConOffset.pop(picosCon)
246 glpkCons = glpk.intArray(length + 1)
247 for i in range(length):
248 glpkCons[i + 1] = offset + i # Index 0 is unused.
250 glpk.glp_del_rows(self.int, length, glpkCons)
252 for otherCon, otherOffset in self._glpkConOffset.items():
253 if otherOffset > offset:
254 self._glpkConOffset[otherCon] -= length
256 def _import_objective(self):
257 import swiglpk as glpk
259 p = self.int
260 direction, objective = self.ext.no
262 # Set optimization direction.
263 if direction == "min":
264 glpk.glp_set_obj_dir(p, glpk.GLP_MIN)
265 else:
266 assert direction == "max"
267 glpk.glp_set_obj_dir(p, glpk.GLP_MAX)
269 # Set objective function shift (index 0).
270 glpk.glp_set_obj_coef(p, 0, objective._constant_coef[0])
272 # Set objective function coefficient of the scalar variable.
273 for picosVar, picosCoef in objective._linear_coefs.items():
274 for localIndex in range(picosVar.dim):
275 if picosCoef[localIndex]:
276 glpkIndex = self._glpkVarOffset[picosVar] + localIndex
277 glpk.glp_set_obj_coef(p, glpkIndex, picosCoef[localIndex])
279 def _reset_objective(self):
280 import swiglpk as glpk
282 p = self.int
284 # Zero the objective.
285 scalarVars = glpk.glp_get_num_cols(p)
286 for i in range(scalarVars + 1): # Index 0 refers to the constant term.
287 glpk.glp_set_obj_coef(p, i, 0)
289 def _import_problem(self):
290 import swiglpk as glpk
292 if self.verbosity() >= 1:
293 glpk.glp_term_out(glpk.GLP_ON)
294 else:
295 glpk.glp_term_out(glpk.GLP_OFF)
297 # Create a problem instance.
298 self.int = glpk.glp_create_prob()
300 # Add a dummy variable since index zero is unused.
301 # This is necessary for ComplexAffineExpression.sparse_row to work.
302 self._glpkVarOffset[self.UNUSED_VAR] = 0
304 # Import variables.
305 for variable in self.ext.variables.values():
306 self._import_variable(variable)
308 # Import constraints.
309 for constraint in self.ext.constraints.values():
310 if not isinstance(constraint, DummyConstraint):
311 self._import_constraint(constraint)
313 # Set objective.
314 self._import_objective()
316 def _update_problem(self):
317 import swiglpk as glpk
319 resetBasis = False
321 for oldConstraint in self._removed_constraints():
322 self._remove_constraint(oldConstraint)
323 resetBasis = True
325 for oldVariable in self._removed_variables():
326 self._remove_variable(oldVariable)
327 resetBasis = True
329 for newVariable in self._new_variables():
330 self._import_variable(newVariable)
332 for newConstraint in self._new_constraints():
333 self._import_constraint(newConstraint)
335 if self._objective_has_changed():
336 self._reset_objective()
337 self._import_objective()
339 if resetBasis:
340 # TODO: Repair the basis in _remove_constraint, _remove_variable.
341 glpk.glp_cpx_basis(self.int)
343 def _solve(self):
344 import swiglpk as glpk
346 p = self.int
348 continuous = self.ext.is_continuous()
349 minimizing = glpk.glp_get_obj_dir(p) == glpk.GLP_MIN
351 # Select LP solver (Simplex or Interior Point Method).
352 if continuous:
353 if self.ext.options.lp_root_method == "interior":
354 interior = True
355 else:
356 # Default to Simplex.
357 interior = False
358 simplex = not interior
359 else:
360 simplex = interior = False
362 # Select appropriate options container.
363 if simplex:
364 options = glpk.glp_smcp()
365 glpk.glp_init_smcp(options)
366 elif interior:
367 options = glpk.glp_iptcp()
368 glpk.glp_init_iptcp(options)
369 else:
370 options = glpk.glp_iocp()
371 glpk.glp_init_iocp(options)
373 # verbosity
374 verbosity = self.verbosity()
375 if verbosity < 0:
376 options.msg_lev = glpk.GLP_MSG_OFF
377 elif verbosity == 0:
378 options.msg_lev = glpk.GLP_MSG_ERR
379 elif verbosity == 1:
380 options.msg_lev = glpk.GLP_MSG_ON
381 elif verbosity >= 2:
382 options.msg_lev = glpk.GLP_MSG_ALL
384 # abs_prim_fsb_tol
385 if self.ext.options.abs_prim_fsb_tol is not None:
386 options.tol_bnd = self.ext.options.abs_prim_fsb_tol
388 # abs_dual_fsb_tol
389 if self.ext.options.abs_dual_fsb_tol is not None:
390 options.tol_dj = self.ext.options.abs_dual_fsb_tol
392 # rel_bnb_opt_tol
393 # Note that the option is silently ignored if passed alongside an LP;
394 # while the solver does not allow us to pass the option in that case, it
395 # is still technically a valid option as every LP is also a MIP.
396 if self.ext.options.rel_bnb_opt_tol is not None:
397 if not continuous:
398 options.mip_gap = self.ext.options.rel_bnb_opt_tol
400 # max_iterations
401 if not simplex:
402 self._handle_unsupported_option("max_iterations",
403 "GLPK supports the 'max_iterations' option only with Simplex.")
404 elif self.ext.options.max_iterations is not None:
405 options.it_lim = int(self.ext.options.max_iterations)
407 # lp_root_method
408 # Note that the PICOS option is explicitly also meant for the MIP
409 # preprocessing step but GLPK does not support it in that scenario.
410 if not continuous:
411 self._handle_unsupported_option("lp_root_method",
412 "GLPK supports the 'lp_root_method' option only for LPs.")
413 elif self.ext.options.lp_root_method is not None:
414 if self.ext.options.lp_root_method == "interior":
415 # Handled above.
416 pass
417 elif self.ext.options.lp_root_method == "psimplex":
418 assert simplex
419 options.meth = glpk.GLP_PRIMAL
420 elif self.ext.options.lp_root_method == "dsimplex":
421 assert simplex
422 options.meth = glpk.GLP_DUAL
423 else:
424 assert False, "Unexpected lp_root_method value."
426 # timelimit
427 if interior:
428 self._handle_unsupported_option("timelimit",
429 "GLPK does not support the 'timelimit' option with the "
430 "Interior Point Method.")
431 elif self.ext.options.timelimit is not None:
432 options.tm_lim = int(1000 * self.ext.options.timelimit)
434 # Handle unsupported options.
435 self._handle_unsupported_options(
436 "lp_node_method", "treememory", "max_fsb_nodes", "hotstart")
438 # TODO: Add GLPK-sepcific options. Candidates are:
439 # For both Simplex and MIPs:
440 # tol_*, out_*
441 # For Simplex:
442 # pricing, r_test, obj_*
443 # For the Interior Point Method:
444 # ord_alg
445 # For MIPs:
446 # *_tech, *_heur, ps_tm_lim, *_cuts, cb_size, binarize
448 # Attempt to solve the problem.
449 with self._header():
450 with self._stopwatch():
451 if simplex:
452 # TODO: Support glp_exact.
453 error = glpk.glp_simplex(p, options)
454 elif interior:
455 error = glpk.glp_interior(p, options)
456 else:
457 options.presolve = glpk.GLP_ON
458 error = glpk.glp_intopt(p, options)
460 # Conert error codes to text output.
461 # Note that by printing it above the footer, this output is made to
462 # look like it's coming from GLPK, which is technically wrong but
463 # semantically correct.
464 if error == glpk.GLP_EBADB:
465 self._warn("Unable to start the search, because the initial "
466 "basis specified in the problem object is invalid.")
467 elif error == glpk.GLP_ESING:
468 self._warn("Unable to start the search, because the basis "
469 "matrix corresponding to the initial basis is singular "
470 "within the working precision.")
471 elif error == glpk.GLP_ECOND:
472 self._warn("Unable to start the search, because the basis "
473 "matrix corresponding to the initial basis is "
474 "ill-conditioned.")
475 elif error == glpk.GLP_EBOUND:
476 self._warn("Unable to start the search, because some double-"
477 "bounded variables have incorrect bounds.")
478 elif error == glpk.GLP_EFAIL:
479 self._warn("The search was prematurely terminated due to a "
480 "solver failure.")
481 elif error == glpk.GLP_EOBJLL:
482 self._warn("The search was prematurely terminated, because the "
483 "objective function being maximized has reached its lower "
484 "limit and continues decreasing.")
485 elif error == glpk.GLP_EOBJUL:
486 self._warn("The search was prematurely terminated, because the "
487 "objective function being minimized has reached its upper "
488 "limit and continues increasing.")
489 elif error == glpk.GLP_EITLIM:
490 self._warn("The search was prematurely terminated, because the "
491 "simplex iteration limit has been exceeded.")
492 elif error == glpk.GLP_ETMLIM:
493 self._warn("The search was prematurely terminated, because the "
494 "time limit has been exceeded.")
495 elif error == glpk.GLP_ENOPFS:
496 self._verbose("The LP has no primal feasible solution.")
497 elif error == glpk.GLP_ENODFS:
498 self._verbose("The LP has no dual feasible solution.")
499 elif error != 0:
500 self._warn("GLPK error {:d}.".format(error))
502 # Retrieve primals.
503 primals = {}
504 if self.ext.options.primals is not False:
505 for variable in self.ext.variables.values():
506 value = []
508 for localIndex in range(variable.dim):
509 glpkIndex = self._glpkVarOffset[variable] + localIndex
511 if simplex:
512 localValue = glpk.glp_get_col_prim(p, glpkIndex)
513 elif interior:
514 localValue = glpk.glp_ipt_col_prim(p, glpkIndex)
515 else:
516 localValue = glpk.glp_mip_col_val(p, glpkIndex)
518 value.append(localValue)
520 primals[variable] = value
522 # Retrieve duals.
523 duals = {}
524 if self.ext.options.duals is not False and continuous:
525 for constraint in self.ext.constraints.values():
526 if isinstance(constraint, DummyConstraint):
527 duals[constraint] = cvxopt.spmatrix(
528 [], [], [], constraint.size)
529 continue
531 value = []
533 for localIndex in range(len(constraint)):
534 glpkIndex = self._glpkConOffset[constraint] + localIndex
536 if simplex:
537 localValue = glpk.glp_get_row_dual(p, glpkIndex)
538 elif interior:
539 localValue = glpk.glp_ipt_row_dual(p, glpkIndex)
540 else:
541 assert False
543 value.append(localValue)
545 dual = cvxopt.matrix(value, constraint.size)
546 if (not constraint.is_increasing()) ^ minimizing:
547 dual = -dual
549 duals[constraint] = dual
551 # Retrieve objective value.
552 if simplex:
553 value = glpk.glp_get_obj_val(p)
554 elif interior:
555 value = glpk.glp_ipt_obj_val(p)
556 else:
557 value = glpk.glp_mip_obj_val(p)
559 # Retrieve solution status.
560 if simplex:
561 probStatusCode = glpk.glp_get_status(p)
562 elif interior:
563 probStatusCode = glpk.glp_ipt_status(p)
564 else:
565 probStatusCode = glpk.glp_mip_status(p)
567 if probStatusCode == glpk.GLP_OPT:
568 # simplex, interior, mip
569 problemStatus = PS_FEASIBLE
570 elif probStatusCode == glpk.GLP_FEAS:
571 # simplex, mip
572 problemStatus = PS_FEASIBLE
573 elif probStatusCode == glpk.GLP_INFEAS:
574 # simplex, interior
575 problemStatus = PS_UNKNOWN
576 elif probStatusCode == glpk.GLP_NOFEAS:
577 # simplex, interior, mip
578 problemStatus = PS_INFEASIBLE
579 elif probStatusCode == glpk.GLP_UNBND:
580 # simplex
581 problemStatus = PS_UNBOUNDED
582 elif probStatusCode == glpk.GLP_UNDEF:
583 # simplex, interior, mip
584 problemStatus = PS_UNKNOWN
585 else:
586 problemStatus = PS_UNKNOWN
588 if simplex:
589 prmlStatusCode = glpk.glp_get_prim_stat(p)
590 dualStatusCode = glpk.glp_get_dual_stat(p)
592 if prmlStatusCode == glpk.GLP_FEAS:
593 if probStatusCode == glpk.GLP_OPT:
594 prmlStatus = SS_OPTIMAL
595 else:
596 prmlStatus = SS_FEASIBLE
597 elif prmlStatusCode == glpk.GLP_INFEAS:
598 prmlStatus = SS_INFEASIBLE
599 elif prmlStatusCode == glpk.GLP_NOFEAS:
600 prmlStatus = PS_INFEASIBLE
601 elif prmlStatusCode == glpk.GLP_UNDEF:
602 prmlStatus = SS_EMPTY
603 else:
604 prmlStatus = SS_UNKNOWN
606 if dualStatusCode == glpk.GLP_FEAS:
607 if probStatusCode == glpk.GLP_OPT:
608 dualStatus = SS_OPTIMAL
609 else:
610 dualStatus = SS_FEASIBLE
611 elif dualStatusCode == glpk.GLP_INFEAS:
612 dualStatus = SS_INFEASIBLE
613 elif dualStatusCode == glpk.GLP_NOFEAS:
614 dualStatus = PS_INFEASIBLE
615 elif dualStatusCode == glpk.GLP_UNDEF:
616 dualStatus = SS_EMPTY
617 else:
618 dualStatus = SS_UNKNOWN
619 elif interior:
620 if probStatusCode == glpk.GLP_UNDEF:
621 prmlStatus = SS_EMPTY
622 dualStatus = SS_EMPTY
623 elif probStatusCode == glpk.GLP_OPT:
624 prmlStatus = SS_OPTIMAL
625 dualStatus = SS_OPTIMAL
626 elif probStatusCode == glpk.GLP_FEAS:
627 prmlStatus = SS_FEASIBLE
628 dualStatus = SS_FEASIBLE
629 elif probStatusCode == glpk.GLP_NOFEAS:
630 prmlStatus = SS_INFEASIBLE
631 dualStatus = SS_INFEASIBLE
632 else:
633 prmlStatus = SS_UNKNOWN
634 dualStatus = SS_UNKNOWN
635 else: # MIP
636 if probStatusCode == glpk.GLP_UNDEF:
637 prmlStatus = SS_EMPTY
638 elif probStatusCode == glpk.GLP_OPT:
639 prmlStatus = SS_OPTIMAL
640 elif probStatusCode == glpk.GLP_FEAS:
641 prmlStatus = SS_FEASIBLE
642 elif probStatusCode == glpk.GLP_NOFEAS:
643 prmlStatus = SS_INFEASIBLE
644 else:
645 prmlStatus = SS_UNKNOWN
647 dualStatus = SS_EMPTY
649 return self._make_solution(
650 value, primals, duals, prmlStatus, dualStatus, problemStatus)
653# --------------------------------------
654__all__ = api_end(_API_START, globals())