Coverage for picos/solvers/solver_scip.py : 69.66%
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) 2018-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:`SCIPSolver`."""
21import cvxopt
23from ..apidoc import api_end, api_start
24from ..constraints import (AffineConstraint, ConvexQuadraticConstraint,
25 DummyConstraint, NonconvexQuadraticConstraint,
26 RSOCConstraint, SOCConstraint)
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_INFEASIBLE, SS_OPTIMAL,
31 SS_PREMATURE, SS_UNKNOWN)
32from .solver import OptionValueError, Solver
34_API_START = api_start(globals())
35# -------------------------------
38class SCIPSolver(Solver):
39 """Interface to the SCIP solver via the PySCIPOpt Python interface."""
41 SUPPORTED = Specification(
42 objectives=[
43 AffineExpression],
44 constraints=[
45 DummyConstraint,
46 AffineConstraint,
47 SOCConstraint,
48 RSOCConstraint,
49 ConvexQuadraticConstraint,
50 NonconvexQuadraticConstraint])
52 @classmethod
53 def supports(cls, footprint, explain=False):
54 """Implement :meth:`~.solver.Solver.supports`."""
55 result = Solver.supports(footprint, explain)
56 if not result or (explain and not result[0]):
57 return result
59 # HACK: SCIP has trouble with continuous problems, in particular it only
60 # produces dual values for LPs (wehich are sometimes incorrect).
61 # Nonconvex quadratics are supported as no other solver can deal
62 # with them at this point.
63 if footprint.continuous \
64 and not ("con", NonconvexQuadraticConstraint) in footprint:
65 if explain:
66 return (False,
67 "Continuous problems, except for nonconvex quadratic ones "
68 "(PICOS' choice).")
69 else:
70 return False
72 if footprint not in cls.SUPPORTED:
73 if explain:
74 return False, cls.SUPPORTED.mismatch_reason(footprint)
75 else:
76 return False
78 return (True, None) if explain else True
80 @classmethod
81 def default_penalty(cls):
82 """Implement :meth:`~.solver.Solver.default_penalty`."""
83 return 1.0 # Stable free/open source solver.
85 @classmethod
86 def test_availability(cls):
87 """Implement :meth:`~.solver.Solver.test_availability`."""
88 cls.check_import("pyscipopt")
90 @classmethod
91 def names(cls):
92 """Implement :meth:`~.solver.Solver.names`."""
93 return "scip", "SCIP", "SCIP Optimization Suite"
95 @classmethod
96 def is_free(cls):
97 """Implement :meth:`~.solver.Solver.is_free`."""
98 return False
100 def __init__(self, problem):
101 """Initialize a SCIP solver interface.
103 :param ~picos.Problem problem: The problem to be solved.
104 """
105 super(SCIPSolver, self).__init__(problem)
107 self._scipVar = dict()
108 """Maps PICOS variable indices to SCIP variables."""
110 self._scipCons = dict()
111 """
112 Maps PICOS constraints to lists of SCIP constraints.
114 For PICOS second order cone constraints, the first entry in the list is
115 a SCIP quadratic constraint and the second entry is a SCIP linear
116 auxiliary constraint.
117 """
119 self._scipQuadObjAuxVar = None
120 """A SCIP auxiliary variable to support a PICOS quadratic objective."""
122 self._scipQuadObjAuxCon = None
123 """A SCIP auxiliary constraint to support a PICOS quadr. objective."""
125 def reset_problem(self):
126 """Implement :meth:`~.solver.Solver.reset_problem`."""
127 self.int = None
129 self._scipVar.clear()
130 self._scipCons.clear()
132 self._scipQuadObjAuxVar = None
133 self._scipQuadObjAuxCon = None
135 def _make_scip_var_names(self, picosVar, localIndex=None):
136 """Make SCIP variable names.
138 Converts a PICOS variable to a list of SCIP variable names, each
139 corresponding to one scalar variable contained in the PICOS variable.
140 If localIndex is given, then only the name of the SCIP variable
141 representing the scalar variable with that offset is returned.
142 The name format is "picosName[localIndex]".
143 """
144 # TODO: This function appears in multiple solvers, move it to the Solver
145 # base class as "_make_scalar_var_names".
146 if localIndex is not None:
147 return "{}[{}]".format(picosVar.name, localIndex)
148 else:
149 return [
150 self._make_scip_var_names(picosVar, localIndex)
151 for localIndex in range(len(picosVar))]
153 def _import_variable(self, picosVar):
154 from math import ceil, floor
156 dim = picosVar.dim
157 inf = self.int.infinity()
159 # Retrieve type code.
160 if isinstance(picosVar, (BinaryVariable, IntegerVariable)):
161 scipVarType = "I"
162 else:
163 scipVarType = "C"
165 # Retrieve bounds.
166 lowerBounds = [-inf]*dim
167 upperBounds = [inf]*dim
168 lower, upper = picosVar.bound_dicts
169 for i, b in lower.items():
170 lowerBounds[i] = b
171 for i, b in upper.items():
172 upperBounds[i] = b
174 # Refine bounds.
175 if isinstance(picosVar, IntegerVariable):
176 lowerBounds = [int(ceil(b)) for b in lowerBounds]
177 upperBounds = [int(floor(b)) for b in upperBounds]
179 # Give names.
180 scipNames = self._make_scip_var_names(picosVar)
182 # Import scalar variables.
183 for i in range(dim):
184 self._scipVar[picosVar.id_at(i)] = self.int.addVar(
185 name=scipNames[i], vtype=scipVarType,
186 lb=lowerBounds[i], ub=upperBounds[i])
188 def _import_variable_values(self):
189 # TODO: Import variable values with SCIP.
190 raise NotImplementedError
192 def _reset_variable_values(self):
193 # TODO: Import variable values with SCIP.
194 raise NotImplementedError
196 def _affinexp_pic2scip(self, picosExpression):
197 """Transform an affine expression from PICOS to SCIP.
199 Multidimensional PICOS expressions are converted to a number of scalar
200 SCIP expressions.
202 :returns: A :class:`list` of :class:`SCIP expressions
203 <pyscipopt.scip.Expr>`.
204 """
205 from pyscipopt.scip import Expr
207 if picosExpression is None:
208 yield Expr()
209 return
211 for scipVars, coefficients, constant in picosExpression.sparse_rows(
212 None, indexFunction=lambda picosVar, localIndex:
213 self._scipVar[picosVar.id_at(localIndex)]):
214 scipExpression = Expr()
215 for scipVar, coefficient in zip(scipVars, coefficients):
216 scipExpression += coefficient * scipVar
217 scipExpression += constant
218 yield scipExpression
220 def _scalar_affinexp_pic2scip(self, picosExpression):
221 """Transform a PICOS scalar affine expression into a SCIP expression.
223 :returns: A :class:`SCIP expression <pyscipopt.scip.Expr>`.
224 """
225 assert len(picosExpression) == 1
226 return next(self._affinexp_pic2scip(picosExpression))
228 def _quadexp_pic2scip(self, picosExpression):
229 # Convert affine part.
230 scipExpression = self._scalar_affinexp_pic2scip(picosExpression.aff)
232 # Convert quadratic part.
233 for (pVar1, pVar2), pCoefficients \
234 in picosExpression._sparse_quads.items():
235 for sparseIndex in range(len(pCoefficients)):
236 localVar1Index = pCoefficients.I[sparseIndex]
237 localVar2Index = pCoefficients.J[sparseIndex]
238 localCoefficient = pCoefficients.V[sparseIndex]
239 scipVar1 = self._scipVar[pVar1.id_at(localVar1Index)]
240 scipVar2 = self._scipVar[pVar2.id_at(localVar2Index)]
241 scipExpression += localCoefficient * scipVar1 * scipVar2
243 return scipExpression
245 def _import_linear_constraint(self, picosConstraint):
246 assert isinstance(picosConstraint, AffineConstraint)
248 scipCons = []
249 picosLHS = picosConstraint.lhs - picosConstraint.rhs
250 for scipLHS in self._affinexp_pic2scip(picosLHS):
251 if picosConstraint.is_increasing():
252 scipCons.append(self.int.addCons(scipLHS <= 0.0))
253 elif picosConstraint.is_decreasing():
254 scipCons.append(self.int.addCons(scipLHS >= 0.0))
255 elif picosConstraint.is_equality():
256 scipCons.append(self.int.addCons(scipLHS == 0.0))
257 else:
258 assert False, "Unexpected constraint relation."
260 return scipCons
262 def _import_quad_constraint(self, picosConstraint):
263 assert isinstance(picosConstraint,
264 (ConvexQuadraticConstraint, NonconvexQuadraticConstraint))
266 scipLHS = self._quadexp_pic2scip(picosConstraint.le0)
267 scipCons = [self.int.addCons(scipLHS <= 0.0)]
269 return scipCons
271 def _import_socone_constraint(self, picosConstraint):
272 scipCons = []
273 scipQuadLHS = self._quadexp_pic2scip(
274 (picosConstraint.ne | picosConstraint.ne)
275 - (picosConstraint.ub * picosConstraint.ub))
276 scipCons.append(self.int.addCons(scipQuadLHS <= 0.0))
278 scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub)
279 if scipAuxLHS.degree() > 0:
280 scipCons.append(self.int.addCons(scipAuxLHS >= 0.0))
282 return scipCons
284 def _import_rscone_constraint(self, picosConstraint):
285 scipCons = []
286 scipLHS = self._quadexp_pic2scip(
287 (picosConstraint.ne | picosConstraint.ne)
288 - (picosConstraint.ub1 * picosConstraint.ub2))
289 scipCons.append(self.int.addCons(scipLHS <= 0.0))
291 # Make sure that either the RHS is constant, or one of the two
292 # expressions is non-negative.
293 scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub1)
294 if scipAuxLHS.degree() > 0:
295 scipCons.append(self.int.addCons(scipAuxLHS >= 0.0))
296 else:
297 scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub2)
298 if scipAuxLHS.degree() > 0:
299 scipCons.append(self.int.addCons(scipAuxLHS >= 0.0))
301 return scipCons
303 def _import_constraint(self, picosConstraint):
304 # Import constraint based on type.
305 if isinstance(picosConstraint, AffineConstraint):
306 scipCons = self._import_linear_constraint(picosConstraint)
307 elif isinstance(picosConstraint,
308 (ConvexQuadraticConstraint, NonconvexQuadraticConstraint)):
309 scipCons = self._import_quad_constraint(picosConstraint)
310 elif isinstance(picosConstraint, SOCConstraint):
311 scipCons = self._import_socone_constraint(picosConstraint)
312 elif isinstance(picosConstraint, RSOCConstraint):
313 scipCons = self._import_rscone_constraint(picosConstraint)
314 else:
315 assert isinstance(picosConstraint, DummyConstraint), \
316 "Unexpected constraint type: {}".format(
317 picosConstraint.__class__.__name__)
319 # Map PICOS constraints to lists of SCIP constraints.
320 self._scipCons[picosConstraint] = scipCons
322 def _import_objective(self):
323 picosSense, picosObjective = self.ext.no
325 # Retrieve objective sense.
326 if picosSense == "min":
327 scipSense = "minimize"
328 else:
329 assert picosSense == "max"
330 scipSense = "maximize"
332 # Import objective function.
333 scipObjective = self._scalar_affinexp_pic2scip(picosObjective)
335 self.int.setObjective(scipObjective, scipSense)
337 def _import_problem(self):
338 import pyscipopt as scip
340 # Create a problem instance.
341 self.int = scip.Model()
343 # Import variables.
344 for variable in self.ext.variables.values():
345 self._import_variable(variable)
347 # Import constraints.
348 for constraint in self.ext.constraints.values():
349 self._import_constraint(constraint)
351 # Set objective.
352 self._import_objective()
354 def _update_problem(self):
355 # TODO: Support all problem updates supported by SCIP.
356 raise NotImplementedError
358 def _solve(self):
359 import pyscipopt as scip
361 # TODO: Give Problem an interface for checks like this.
362 isLP = isinstance(self.ext.no.function, AffineExpression) \
363 and all(isinstance(constraint, AffineConstraint)
364 for constraint in self.ext.constraints.values())
366 # Reset options.
367 self.int.resetParams()
369 # verbosity
370 picosVerbosity = self.verbosity()
371 if picosVerbosity <= -1:
372 scipVerbosity = 0
373 elif picosVerbosity == 0:
374 scipVerbosity = 2
375 elif picosVerbosity == 1:
376 scipVerbosity = 3
377 elif picosVerbosity >= 2:
378 scipVerbosity = 5
379 self.int.setIntParam("display/verblevel", scipVerbosity)
381 # TODO: "numerics/lpfeastol" has been replaced in SCIP 7.0.0 (PySCIPOpt
382 # 3.0.0) by "numerics/lpfeastolfactor", which is multiplied with
383 # some unknown default feasibility (maybe "numerics/feastol"?).
384 # Looking at it now it is not clear to me which of SCIP's
385 # feasibility tolerances does what exactly and whether they are
386 # absolute or relative measures. The workaround is to maintain
387 # PICOS' old guess but ignore abs_prim_fsb_tol with the new SCIP.
388 if int(scip.__version__.split(".")[0]) < 3:
389 # abs_prim_fsb_tol
390 if self.ext.options.abs_prim_fsb_tol is not None:
391 self.int.setRealParam("numerics/lpfeastol",
392 self.ext.options.abs_prim_fsb_tol)
394 # abs_dual_fsb_tol
395 if self.ext.options.abs_dual_fsb_tol is not None:
396 self.int.setRealParam("numerics/dualfeastol",
397 self.ext.options.abs_dual_fsb_tol)
399 # rel_prim_fsb_tol, rel_dual_fsb_tol
400 relFsbTols = [tol for tol in (self.ext.options.rel_prim_fsb_tol,
401 self.ext.options.rel_dual_fsb_tol) if tol is not None]
402 if relFsbTols:
403 self.int.setRealParam("numerics/feastol", min(relFsbTols))
405 # abs_ipm_opt_tol, abs_bnb_opt_tol
406 absOptTols = [tol for tol in (self.ext.options.abs_ipm_opt_tol,
407 self.ext.options.abs_bnb_opt_tol) if tol is not None]
408 if absOptTols:
409 self.int.setRealParam("limits/absgap", min(absOptTols))
411 # rel_ipm_opt_tol, rel_bnb_opt_tol
412 relOptTols = [tol for tol in (self.ext.options.rel_ipm_opt_tol,
413 self.ext.options.rel_bnb_opt_tol) if tol is not None]
414 if relOptTols:
415 self.int.setRealParam("limits/gap", min(relOptTols))
417 # timelimit
418 if self.ext.options.timelimit is not None:
419 self.int.setRealParam("limits/time", self.ext.options.timelimit)
421 # treememory
422 if self.ext.options.treememory is not None:
423 self.int.setRealParam("limits/memory", self.ext.options.treememory)
425 # max_fsb_nodes
426 if self.ext.options.max_fsb_nodes is not None:
427 self.int.setRealParam(
428 "limits/solutions", float(self.ext.options.max_fsb_nodes))
430 # Handle SCIP-specific options.
431 for key, value in self.ext.options.scip_params.items():
432 try:
433 if isinstance(value, bool):
434 self.int.setBoolParam(key, value)
435 elif isinstance(value, str):
436 if len(value) == 1:
437 try:
438 self.int.setCharParam(key, ord(value))
439 except LookupError:
440 self.int.setStringParam(key, value)
441 else:
442 self.int.setStringParam(key, value)
443 elif isinstance(value, float):
444 self.int.setRealParam(key, value)
445 elif isinstance(value, int):
446 try:
447 self.int.setIntParam(key, value)
448 except LookupError:
449 try:
450 self.int.setLongintParam(key, value)
451 except LookupError:
452 self.int.setRealParam(key, float(value))
453 else:
454 raise TypeError(
455 "The value '{}' has an uncommon type.".format(value))
456 except KeyError as error:
457 self._handle_bad_solver_specific_option_key(key, error)
458 except ValueError as error:
459 self._handle_bad_solver_specific_option_value(key, value, error)
460 except (TypeError, LookupError) as error:
461 picos_option = "{}_params".format(self.name)
462 assert picos_option in self.ext.options, \
463 "The PICOS option '{}' does not exist.".format(picos_option)
465 raise OptionValueError(
466 "Failed to guess type of {} option '{}' set via '{}'."
467 .format(self.display_name, key, picos_option)) from error
469 # Handle unsupported options.
470 self._handle_unsupported_options(
471 "hotstart", "lp_node_method", "lp_root_method", "max_iterations")
473 # In the case of a pure LP, disable presolve to get duals.
474 if self.ext.options.duals is not False and isLP:
475 self._debug("Disabling SCIP's presolve, heuristics, and propagation"
476 " features to allow for LP duals.")
478 self.int.setPresolve(scip.SCIP_PARAMSETTING.OFF)
479 self.int.setHeuristics(scip.SCIP_PARAMSETTING.OFF)
481 # Note that this is a helper to set options, so they get reset at
482 # the beginning of the function instead of in the else-scope below.
483 self.int.disablePropagation()
484 else:
485 self.int.setPresolve(scip.SCIP_PARAMSETTING.DEFAULT)
486 self.int.setHeuristics(scip.SCIP_PARAMSETTING.DEFAULT)
488 # Attempt to solve the problem.
489 with self._header(), self._stopwatch():
490 self.int.optimize()
492 # Retrieve primals.
493 primals = {}
494 if self.ext.options.primals is not False:
495 for picosVar in self.ext.variables.values():
496 value = []
497 for localIndex in range(picosVar.dim):
498 scipVar = self._scipVar[picosVar.id_at(localIndex)]
499 value.append(self.int.getVal(scipVar))
501 primals[picosVar] = value
503 # Retrieve duals for LPs.
504 duals = {}
505 if self.ext.options.duals is not False and isLP:
506 for picosConstraint in self.ext.constraints.values():
507 if isinstance(picosConstraint, DummyConstraint):
508 duals[picosConstraint] = cvxopt.spmatrix(
509 [], [], [], picosConstraint.size)
510 continue
512 assert isinstance(picosConstraint, AffineConstraint)
514 # Retrieve dual value for constraint.
515 scipDuals = []
516 for scipCon in self._scipCons[picosConstraint]:
517 scipDuals.append(self.int.getDualsolLinear(scipCon))
518 picosDual = cvxopt.matrix(scipDuals, picosConstraint.size)
520 # Flip sign based on constraint relation.
521 if not picosConstraint.is_increasing():
522 picosDual = -picosDual
524 # Flip sign based on objective sense.
525 if picosDual and self.ext.no.direction == "min":
526 picosDual = -picosDual
528 duals[picosConstraint] = picosDual
530 # Retrieve objective value.
531 value = self.int.getObjVal()
533 # Retrieve solution status.
534 status = self.int.getStatus()
535 if status == "optimal":
536 primalStatus = SS_OPTIMAL
537 dualStatus = SS_OPTIMAL
538 problemStatus = PS_FEASIBLE
539 elif status == "timelimit":
540 primalStatus = SS_PREMATURE
541 dualStatus = SS_PREMATURE
542 problemStatus = PS_UNKNOWN
543 elif status == "infeasible":
544 primalStatus = SS_INFEASIBLE
545 dualStatus = SS_UNKNOWN
546 problemStatus = PS_INFEASIBLE
547 elif status == "unbounded":
548 primalStatus = SS_UNKNOWN
549 dualStatus = SS_INFEASIBLE
550 problemStatus = PS_UNBOUNDED
551 else:
552 primalStatus = SS_UNKNOWN
553 dualStatus = SS_UNKNOWN
554 problemStatus = PS_UNKNOWN
556 return self._make_solution(
557 value, primals, duals, primalStatus, dualStatus, problemStatus)
560# --------------------------------------
561__all__ = api_end(_API_START, globals())