Coverage for picos/solvers/solver_glpk.py : 79.14%
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# ------------------------------------------------------------------------------
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`."""
21import cvxopt
23from ..apidoc import api_end, api_start
24from ..constraints import AffineConstraint, DummyConstraint
25from ..expressions import AffineExpression, BinaryVariable, IntegerVariable
26from ..modeling.footprint import Specification
27from ..modeling.solution import (PS_FEASIBLE, PS_INFEASIBLE, PS_UNBOUNDED,
28 PS_UNKNOWN, SS_EMPTY, SS_FEASIBLE,
29 SS_INFEASIBLE, SS_OPTIMAL, SS_UNKNOWN)
30from .solver import Solver
32_API_START = api_start(globals())
33# -------------------------------
36class GLPKSolver(Solver):
37 """Interface to the GLPK solver via swiglpk."""
39 SUPPORTED = Specification(
40 objectives=[
41 AffineExpression],
42 constraints=[
43 DummyConstraint,
44 AffineConstraint])
46 @classmethod
47 def supports(cls, footprint, explain=False):
48 """Implement :meth:`~.solver.Solver.supports`."""
49 result = Solver.supports(footprint, explain)
50 if not result or (explain and not result[0]):
51 return result
53 if footprint not in cls.SUPPORTED:
54 if explain:
55 return False, cls.SUPPORTED.mismatch_reason(footprint)
56 else:
57 return False
59 return (True, None) if explain else True
61 @classmethod
62 def default_penalty(cls):
63 """Implement :meth:`~.solver.Solver.default_penalty`."""
64 return 1.0 # Stable free/open source solver.
66 @classmethod
67 def test_availability(cls):
68 """Implement :meth:`~.solver.Solver.test_availability`."""
69 cls.check_import("swiglpk")
71 @classmethod
72 def names(cls):
73 """Implement :meth:`~.solver.Solver.names`."""
74 return "glpk", "GLPK", "GNU Linear Programming Kit"
76 @classmethod
77 def is_free(cls):
78 """Implement :meth:`~.solver.Solver.is_free`."""
79 return True
81 def __init__(self, problem):
82 """Initialize a GLPK solver interface.
84 :param ~picos.Problem problem: The problem to be solved.
85 """
86 self._glpkVarOffset = {}
87 self._glpkConOffset = {}
89 super(GLPKSolver, self).__init__(problem)
91 def __del__(self):
92 try:
93 import swiglpk as glpk
94 except ImportError:
95 # Happens when python is shutting down. In this case, no garbage
96 # collection is necessary, anyway.
97 return
99 if self.int is not None:
100 glpk.glp_delete_prob(self.int)
102 def reset_problem(self):
103 """Implement :meth:`~.solver.Solver.reset_problem`."""
104 import swiglpk as glpk
106 self._glpkVarOffset.clear()
107 self._glpkConOffset.clear()
109 if self.int is not None:
110 glpk.glp_delete_prob(self.int)
111 self.int = None
113 def _import_variable(self, picosVar):
114 import swiglpk as glpk
116 p = self.int
117 dim = picosVar.dim
119 # Append columns to the constraint matrix.
120 offset = glpk.glp_add_cols(p, dim)
122 # Retrieve the type.
123 if isinstance(picosVar, IntegerVariable):
124 varType = glpk.GLP_IV
125 elif isinstance(picosVar, BinaryVariable):
126 varType = glpk.GLP_BV
127 else:
128 varType = glpk.GLP_CV
130 # Retrieve bounds.
131 if not isinstance(picosVar, BinaryVariable):
132 boundKeys = [glpk.GLP_FR]*dim
133 lowerBounds = [0]*dim
134 upperBounds = [0]*dim
136 lower, upper = picosVar.bound_dicts
138 for i, b in lower.items():
139 boundKeys[i] = glpk.GLP_LO
140 lowerBounds[i] = b
142 for i, b in upper.items():
143 if boundKeys[i] == glpk.GLP_FR:
144 boundKeys[i] = glpk.GLP_UP
145 else: # Also has a lower bound.
146 if lowerBounds[i] == b:
147 boundKeys[i] = glpk.GLP_FX
148 else:
149 boundKeys[i] = glpk.GLP_DB
151 upperBounds[i] = b
153 # Import scalar variables.
154 for i in range(picosVar.dim):
155 glpkIndex = offset + i
157 # Assign a name.
158 glpk.glp_set_col_name(
159 p, glpkIndex, "{}_{}".format(picosVar.name, i))
161 # Set the type.
162 glpk.glp_set_col_kind(p, glpkIndex, varType)
164 # Assign the bounds.
165 if not isinstance(picosVar, BinaryVariable):
166 glpk.glp_set_col_bnds(
167 p, glpkIndex, boundKeys[i], lowerBounds[i], upperBounds[i])
169 self._glpkVarOffset[picosVar] = offset
171 def _remove_variable(self, picosVar):
172 import swiglpk as glpk
174 dim = picosVar.dim
176 offset = self._glpkVarOffset.pop(picosVar)
178 glpkVars = glpk.intArray(dim + 1)
179 for i in range(dim):
180 glpkVars[i + 1] = offset + i # Index 0 is unused.
182 glpk.glp_del_cols(self.int, dim, glpkVars)
184 for otherVar, otherOffset in self._glpkVarOffset.items():
185 if otherOffset > offset:
186 self._glpkVarOffset[otherVar] -= dim
188 def _import_constraint(self, picosCon):
189 import swiglpk as glpk
191 p = self.int
193 # Append rows to the constraint matrix.
194 rowOffset = glpk.glp_add_rows(p, len(picosCon))
196 # Import scalar constraints.
197 rows = picosCon.sparse_Ab_rows(self._glpkVarOffset)
198 for localConIndex, (glpkVarIndices, coefs, rhs) in enumerate(rows):
199 glpkConIndex = rowOffset + localConIndex
200 numColumns = len(glpkVarIndices)
202 # Assign a name.
203 glpk.glp_set_row_name(
204 p, glpkConIndex, "{}_{}".format(picosCon.id, localConIndex))
206 # Set the constant term.
207 if picosCon.is_equality():
208 glpk.glp_set_row_bnds(p, glpkConIndex, glpk.GLP_FX, rhs, rhs)
209 elif picosCon.is_increasing():
210 glpk.glp_set_row_bnds(p, glpkConIndex, glpk.GLP_UP, 0, rhs)
211 elif picosCon.is_decreasing():
212 glpk.glp_set_row_bnds(p, glpkConIndex, glpk.GLP_LO, rhs, 0)
213 else:
214 assert False, "Unexpected constraint relation."
216 # Set coefficients.
217 # NOTE: GLPK requires a glpk.intArray containing column indices and
218 # a glpk.doubleArray of same size containing the coefficients
219 # for the listed column index. The first element of both
220 # arrays (with index 0) is skipped by GLPK.
221 glpkVarIndicesArray = glpk.intArray(numColumns + 1)
222 for i in range(numColumns):
223 glpkVarIndicesArray[i + 1] = glpkVarIndices[i]
225 coefficientsArray = glpk.doubleArray(numColumns + 1)
226 for i in range(numColumns):
227 coefficientsArray[i + 1] = coefs[i]
229 glpk.glp_set_mat_row(p, glpkConIndex, numColumns,
230 glpkVarIndicesArray, coefficientsArray)
232 self._glpkConOffset[picosCon] = rowOffset
234 def _remove_constraint(self, picosCon):
235 import swiglpk as glpk
237 length = len(picosCon)
239 offset = self._glpkConOffset.pop(picosCon)
241 glpkCons = glpk.intArray(length + 1)
242 for i in range(length):
243 glpkCons[i + 1] = offset + i # Index 0 is unused.
245 glpk.glp_del_rows(self.int, length, glpkCons)
247 for otherCon, otherOffset in self._glpkConOffset.items():
248 if otherOffset > offset:
249 self._glpkConOffset[otherCon] -= length
251 def _import_objective(self):
252 import swiglpk as glpk
254 p = self.int
255 direction, objective = self.ext.no
257 # Set optimization direction.
258 if direction == "min":
259 glpk.glp_set_obj_dir(p, glpk.GLP_MIN)
260 else:
261 assert direction == "max"
262 glpk.glp_set_obj_dir(p, glpk.GLP_MAX)
264 # Set objective function shift (index 0).
265 glpk.glp_set_obj_coef(p, 0, objective._constant_coef[0])
267 # Set objective function coefficient of the scalar variable.
268 for picosVar, picosCoef in objective._linear_coefs.items():
269 for localIndex in range(picosVar.dim):
270 if picosCoef[localIndex]:
271 glpkIndex = self._glpkVarOffset[picosVar] + localIndex
272 glpk.glp_set_obj_coef(p, glpkIndex, picosCoef[localIndex])
274 def _reset_objective(self):
275 import swiglpk as glpk
277 p = self.int
279 # Zero the objective.
280 scalarVars = glpk.glp_get_num_cols(p)
281 for i in range(scalarVars + 1): # Index 0 refers to the constant term.
282 glpk.glp_set_obj_coef(p, i, 0)
284 def _import_problem(self):
285 import swiglpk as glpk
287 if self.verbosity() >= 1:
288 glpk.glp_term_out(glpk.GLP_ON)
289 else:
290 glpk.glp_term_out(glpk.GLP_OFF)
292 # Create a problem instance.
293 self.int = glpk.glp_create_prob()
295 # Import variables.
296 for variable in self.ext.variables.values():
297 self._import_variable(variable)
299 # Import constraints.
300 for constraint in self.ext.constraints.values():
301 if not isinstance(constraint, DummyConstraint):
302 self._import_constraint(constraint)
304 # Set objective.
305 self._import_objective()
307 def _update_problem(self):
308 import swiglpk as glpk
310 resetBasis = False
312 for oldConstraint in self._removed_constraints():
313 self._remove_constraint(oldConstraint)
314 resetBasis = True
316 for oldVariable in self._removed_variables():
317 self._remove_variable(oldVariable)
318 resetBasis = True
320 for newVariable in self._new_variables():
321 self._import_variable(newVariable)
323 for newConstraint in self._new_constraints():
324 self._import_constraint(newConstraint)
326 if self._objective_has_changed():
327 self._reset_objective()
328 self._import_objective()
330 if resetBasis:
331 # TODO: Repair the basis in _remove_constraint, _remove_variable.
332 glpk.glp_cpx_basis(self.int)
334 def _solve(self):
335 import swiglpk as glpk
337 p = self.int
339 continuous = self.ext.is_continuous()
340 minimizing = glpk.glp_get_obj_dir(p) == glpk.GLP_MIN
342 # Select LP solver (Simplex or Interior Point Method).
343 if continuous:
344 if self.ext.options.lp_root_method == "interior":
345 interior = True
346 else:
347 # Default to Simplex.
348 interior = False
349 simplex = not interior
350 else:
351 simplex = interior = False
353 # Select appropriate options container.
354 if simplex:
355 options = glpk.glp_smcp()
356 glpk.glp_init_smcp(options)
357 elif interior:
358 options = glpk.glp_iptcp()
359 glpk.glp_init_iptcp(options)
360 else:
361 options = glpk.glp_iocp()
362 glpk.glp_init_iocp(options)
364 # verbosity
365 verbosity = self.verbosity()
366 if verbosity < 0:
367 options.msg_lev = glpk.GLP_MSG_OFF
368 elif verbosity == 0:
369 options.msg_lev = glpk.GLP_MSG_ERR
370 elif verbosity == 1:
371 options.msg_lev = glpk.GLP_MSG_ON
372 elif verbosity >= 2:
373 options.msg_lev = glpk.GLP_MSG_ALL
375 # abs_prim_fsb_tol
376 if self.ext.options.abs_prim_fsb_tol is not None:
377 options.tol_bnd = self.ext.options.abs_prim_fsb_tol
379 # abs_dual_fsb_tol
380 if self.ext.options.abs_dual_fsb_tol is not None:
381 options.tol_dj = self.ext.options.abs_dual_fsb_tol
383 # rel_bnb_opt_tol
384 # Note that the option is silently ignored if passed alongside an LP;
385 # while the solver does not allow us to pass the option in that case, it
386 # is still technically a valid option as every LP is also a MIP.
387 if self.ext.options.rel_bnb_opt_tol is not None:
388 if not continuous:
389 options.mip_gap = self.ext.options.rel_bnb_opt_tol
391 # max_iterations
392 if not simplex:
393 self._handle_unsupported_option("max_iterations",
394 "GLPK supports the 'max_iterations' option only with Simplex.")
395 elif self.ext.options.max_iterations is not None:
396 options.it_lim = int(self.ext.options.max_iterations)
398 # lp_root_method
399 # Note that the PICOS option is explicitly also meant for the MIP
400 # preprocessing step but GLPK does not support it in that scenario.
401 if not continuous:
402 self._handle_unsupported_option("lp_root_method",
403 "GLPK supports the 'lp_root_method' option only for LPs.")
404 elif self.ext.options.lp_root_method is not None:
405 if self.ext.options.lp_root_method == "interior":
406 # Handled above.
407 pass
408 elif self.ext.options.lp_root_method == "psimplex":
409 assert simplex
410 options.meth = glpk.GLP_PRIMAL
411 elif self.ext.options.lp_root_method == "dsimplex":
412 assert simplex
413 options.meth = glpk.GLP_DUAL
414 else:
415 assert False, "Unexpected lp_root_method value."
417 # timelimit
418 if interior:
419 self._handle_unsupported_option("timelimit",
420 "GLPK does not support the 'timelimit' option with the "
421 "Interior Point Method.")
422 elif self.ext.options.timelimit is not None:
423 options.tm_lim = int(1000 * self.ext.options.timelimit)
425 # Handle unsupported options.
426 self._handle_unsupported_options(
427 "lp_node_method", "treememory", "max_fsb_nodes", "hotstart")
429 # TODO: Add GLPK-sepcific options. Candidates are:
430 # For both Simplex and MIPs:
431 # tol_*, out_*
432 # For Simplex:
433 # pricing, r_test, obj_*
434 # For the Interior Point Method:
435 # ord_alg
436 # For MIPs:
437 # *_tech, *_heur, ps_tm_lim, *_cuts, cb_size, binarize
439 # Attempt to solve the problem.
440 with self._header():
441 with self._stopwatch():
442 if simplex:
443 # TODO: Support glp_exact.
444 error = glpk.glp_simplex(p, options)
445 elif interior:
446 error = glpk.glp_interior(p, options)
447 else:
448 options.presolve = glpk.GLP_ON
449 error = glpk.glp_intopt(p, options)
451 # Conert error codes to text output.
452 # Note that by printing it above the footer, this output is made to
453 # look like it's coming from GLPK, which is technically wrong but
454 # semantically correct.
455 if error == glpk.GLP_EBADB:
456 self._warn("Unable to start the search, because the initial "
457 "basis specified in the problem object is invalid.")
458 elif error == glpk.GLP_ESING:
459 self._warn("Unable to start the search, because the basis "
460 "matrix corresponding to the initial basis is singular "
461 "within the working precision.")
462 elif error == glpk.GLP_ECOND:
463 self._warn("Unable to start the search, because the basis "
464 "matrix corresponding to the initial basis is "
465 "ill-conditioned.")
466 elif error == glpk.GLP_EBOUND:
467 self._warn("Unable to start the search, because some double-"
468 "bounded variables have incorrect bounds.")
469 elif error == glpk.GLP_EFAIL:
470 self._warn("The search was prematurely terminated due to a "
471 "solver failure.")
472 elif error == glpk.GLP_EOBJLL:
473 self._warn("The search was prematurely terminated, because the "
474 "objective function being maximized has reached its lower "
475 "limit and continues decreasing.")
476 elif error == glpk.GLP_EOBJUL:
477 self._warn("The search was prematurely terminated, because the "
478 "objective function being minimized has reached its upper "
479 "limit and continues increasing.")
480 elif error == glpk.GLP_EITLIM:
481 self._warn("The search was prematurely terminated, because the "
482 "simplex iteration limit has been exceeded.")
483 elif error == glpk.GLP_ETMLIM:
484 self._warn("The search was prematurely terminated, because the "
485 "time limit has been exceeded.")
486 elif error == glpk.GLP_ENOPFS:
487 self._verbose("The LP has no primal feasible solution.")
488 elif error == glpk.GLP_ENODFS:
489 self._verbose("The LP has no dual feasible solution.")
490 elif error != 0:
491 self._warn("GLPK error {:d}.".format(error))
493 # Retrieve primals.
494 primals = {}
495 if self.ext.options.primals is not False:
496 for variable in self.ext.variables.values():
497 value = []
499 for localIndex in range(variable.dim):
500 glpkIndex = self._glpkVarOffset[variable] + localIndex
502 if simplex:
503 localValue = glpk.glp_get_col_prim(p, glpkIndex)
504 elif interior:
505 localValue = glpk.glp_ipt_col_prim(p, glpkIndex)
506 else:
507 localValue = glpk.glp_mip_col_val(p, glpkIndex)
509 value.append(localValue)
511 primals[variable] = value
513 # Retrieve duals.
514 duals = {}
515 if self.ext.options.duals is not False and continuous:
516 for constraint in self.ext.constraints.values():
517 if isinstance(constraint, DummyConstraint):
518 duals[constraint] = cvxopt.spmatrix(
519 [], [], [], constraint.size)
520 continue
522 value = []
524 for localIndex in range(len(constraint)):
525 glpkIndex = self._glpkConOffset[constraint] + localIndex
527 if simplex:
528 localValue = glpk.glp_get_row_dual(p, glpkIndex)
529 elif interior:
530 localValue = glpk.glp_ipt_row_dual(p, glpkIndex)
531 else:
532 assert False
534 value.append(localValue)
536 dual = cvxopt.matrix(value, constraint.size)
537 if (not constraint.is_increasing()) ^ minimizing:
538 dual = -dual
540 duals[constraint] = dual
542 # Retrieve objective value.
543 if simplex:
544 value = glpk.glp_get_obj_val(p)
545 elif interior:
546 value = glpk.glp_ipt_obj_val(p)
547 else:
548 value = glpk.glp_mip_obj_val(p)
550 # Retrieve solution status.
551 if simplex:
552 probStatusCode = glpk.glp_get_status(p)
553 elif interior:
554 probStatusCode = glpk.glp_ipt_status(p)
555 else:
556 probStatusCode = glpk.glp_mip_status(p)
558 if probStatusCode == glpk.GLP_OPT:
559 # simplex, interior, mip
560 problemStatus = PS_FEASIBLE
561 elif probStatusCode == glpk.GLP_FEAS:
562 # simplex, mip
563 problemStatus = PS_FEASIBLE
564 elif probStatusCode == glpk.GLP_INFEAS:
565 # simplex, interior
566 problemStatus = PS_UNKNOWN
567 elif probStatusCode == glpk.GLP_NOFEAS:
568 # simplex, interior, mip
569 problemStatus = PS_INFEASIBLE
570 elif probStatusCode == glpk.GLP_UNBND:
571 # simplex
572 problemStatus = PS_UNBOUNDED
573 elif probStatusCode == glpk.GLP_UNDEF:
574 # simplex, interior, mip
575 problemStatus = PS_UNKNOWN
576 else:
577 problemStatus = PS_UNKNOWN
579 if simplex:
580 prmlStatusCode = glpk.glp_get_prim_stat(p)
581 dualStatusCode = glpk.glp_get_dual_stat(p)
583 if prmlStatusCode == glpk.GLP_FEAS:
584 if probStatusCode == glpk.GLP_OPT:
585 prmlStatus = SS_OPTIMAL
586 else:
587 prmlStatus = SS_FEASIBLE
588 elif prmlStatusCode == glpk.GLP_INFEAS:
589 prmlStatus = SS_INFEASIBLE
590 elif prmlStatusCode == glpk.GLP_NOFEAS:
591 prmlStatus = PS_INFEASIBLE
592 elif prmlStatusCode == glpk.GLP_UNDEF:
593 prmlStatus = SS_EMPTY
594 else:
595 prmlStatus = SS_UNKNOWN
597 if dualStatusCode == glpk.GLP_FEAS:
598 if probStatusCode == glpk.GLP_OPT:
599 dualStatus = SS_OPTIMAL
600 else:
601 dualStatus = SS_FEASIBLE
602 elif dualStatusCode == glpk.GLP_INFEAS:
603 dualStatus = SS_INFEASIBLE
604 elif dualStatusCode == glpk.GLP_NOFEAS:
605 dualStatus = PS_INFEASIBLE
606 elif dualStatusCode == glpk.GLP_UNDEF:
607 dualStatus = SS_EMPTY
608 else:
609 dualStatus = SS_UNKNOWN
610 elif interior:
611 if probStatusCode == glpk.GLP_UNDEF:
612 prmlStatus = SS_EMPTY
613 dualStatus = SS_EMPTY
614 elif probStatusCode == glpk.GLP_OPT:
615 prmlStatus = SS_OPTIMAL
616 dualStatus = SS_OPTIMAL
617 elif probStatusCode == glpk.GLP_FEAS:
618 prmlStatus = SS_FEASIBLE
619 dualStatus = SS_FEASIBLE
620 elif probStatusCode == glpk.GLP_NOFEAS:
621 prmlStatus = SS_INFEASIBLE
622 dualStatus = SS_INFEASIBLE
623 else:
624 prmlStatus = SS_UNKNOWN
625 dualStatus = SS_UNKNOWN
626 else: # MIP
627 if probStatusCode == glpk.GLP_UNDEF:
628 prmlStatus = SS_EMPTY
629 elif probStatusCode == glpk.GLP_OPT:
630 prmlStatus = SS_OPTIMAL
631 elif probStatusCode == glpk.GLP_FEAS:
632 prmlStatus = SS_FEASIBLE
633 elif probStatusCode == glpk.GLP_NOFEAS:
634 prmlStatus = SS_INFEASIBLE
635 else:
636 prmlStatus = SS_UNKNOWN
638 dualStatus = SS_EMPTY
640 return self._make_solution(
641 value, primals, duals, prmlStatus, dualStatus, problemStatus)
644# --------------------------------------
645__all__ = api_end(_API_START, globals())