# ------------------------------------------------------------------------------
# Copyright (C) 2018-2022 Maximilian Stahlberg
#
# This file is part of PICOS.
#
# PICOS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# this program. If not, see <http://www.gnu.org/licenses/>.
# ------------------------------------------------------------------------------
"""Implementation of :class:`SCIPSolver`."""
import cvxopt
from .. import settings
from ..apidoc import api_end, api_start
from ..constraints import (AffineConstraint, ConvexQuadraticConstraint,
DummyConstraint, NonconvexQuadraticConstraint,
RSOCConstraint, SOCConstraint)
from ..expressions import AffineExpression, BinaryVariable, IntegerVariable
from ..modeling.footprint import Specification
from ..modeling.solution import (PS_FEASIBLE, PS_INFEASIBLE, PS_UNBOUNDED,
PS_UNKNOWN, SS_INFEASIBLE, SS_OPTIMAL,
SS_PREMATURE, SS_UNKNOWN)
from .solver import OptionValueError, Solver
_API_START = api_start(globals())
# -------------------------------
[docs]class SCIPSolver(Solver):
"""Interface to the SCIP solver via the PySCIPOpt Python interface."""
SUPPORTED = Specification(
objectives=[
AffineExpression],
constraints=[
DummyConstraint,
AffineConstraint,
SOCConstraint,
RSOCConstraint,
ConvexQuadraticConstraint,
NonconvexQuadraticConstraint])
[docs] @classmethod
def supports(cls, footprint, explain=False):
"""Implement :meth:`~.solver.Solver.supports`."""
result = Solver.supports(footprint, explain)
if not result or (explain and not result[0]):
return result
# HACK: SCIP has trouble with continuous problems, in particular it only
# produces dual values for LPs (which are sometimes incorrect).
# Nonconvex quadratics are supported as no other solver can deal
# with them at this point.
if not settings.UNRELIABLE_STRATEGIES \
and footprint.continuous \
and not ("con", NonconvexQuadraticConstraint) in footprint:
if explain:
return (False,
"Continuous problems, except for nonconvex quadratic ones "
"(PICOS' choice).")
else:
return False
if footprint not in cls.SUPPORTED:
if explain:
return False, cls.SUPPORTED.mismatch_reason(footprint)
else:
return False
return (True, None) if explain else True
[docs] @classmethod
def default_penalty(cls):
"""Implement :meth:`~.solver.Solver.default_penalty`."""
return 1.0 # Stable free/open source solver.
[docs] @classmethod
def test_availability(cls):
"""Implement :meth:`~.solver.Solver.test_availability`."""
cls.check_import("pyscipopt")
[docs] @classmethod
def names(cls):
"""Implement :meth:`~.solver.Solver.names`."""
return "scip", "SCIP", "SCIP Optimization Suite", "PySCIPOpt"
[docs] @classmethod
def is_free(cls):
"""Implement :meth:`~.solver.Solver.is_free`."""
return False
[docs] def __init__(self, problem):
"""Initialize a SCIP solver interface.
:param ~picos.Problem problem: The problem to be solved.
"""
super(SCIPSolver, self).__init__(problem)
self._scipVars = []
"""A list of all SCIP variables added."""
self._scipVarOffset = dict()
"""Maps PICOS variables to SCIP variable offsets."""
self._scipCons = dict()
"""
Maps PICOS constraints to lists of SCIP constraints.
For PICOS second order cone constraints, the first entry in the list is
a SCIP quadratic constraint and the second entry is a SCIP linear
auxiliary constraint.
"""
self._scipQuadObjAuxVar = None
"""A SCIP auxiliary variable to support a PICOS quadratic objective."""
self._scipQuadObjAuxCon = None
"""A SCIP auxiliary constraint to support a PICOS quadr. objective."""
[docs] def reset_problem(self):
"""Implement :meth:`~.solver.Solver.reset_problem`."""
self.int = None
self._scipVars.clear()
self._scipVarOffset.clear()
self._scipCons.clear()
self._scipQuadObjAuxVar = None
self._scipQuadObjAuxCon = None
def _import_variable(self, picosVar):
from math import ceil, floor
dim = picosVar.dim
inf = self.int.infinity()
# Retrieve type code.
if isinstance(picosVar, (BinaryVariable, IntegerVariable)):
scipVarType = "I"
else:
scipVarType = "C"
# Retrieve bounds.
lowerBounds = [-inf]*dim
upperBounds = [inf]*dim
lower, upper = picosVar.bound_dicts
for i, b in lower.items():
lowerBounds[i] = b
for i, b in upper.items():
upperBounds[i] = b
# Refine bounds.
if isinstance(picosVar, IntegerVariable):
lowerBounds = [int(ceil(b)) for b in lowerBounds]
upperBounds = [int(floor(b)) for b in upperBounds]
# Import scalar variables.
self._scipVarOffset[picosVar] = len(self._scipVars)
for i in range(dim):
self._scipVars.append(self.int.addVar(
vtype=scipVarType, lb=lowerBounds[i], ub=upperBounds[i]))
def _import_variable_values(self):
# TODO: Import variable values with SCIP.
raise NotImplementedError
def _reset_variable_values(self):
# TODO: Import variable values with SCIP.
raise NotImplementedError
def _affinexp_pic2scip(self, picosExpression):
"""Transform an affine expression from PICOS to SCIP.
Multidimensional PICOS expressions are converted to a number of scalar
SCIP expressions.
:returns: A :class:`list` of :class:`SCIP expressions
<pyscipopt.scip.Expr>`.
"""
from pyscipopt.scip import Expr
if picosExpression is None:
yield Expr()
return
rows = picosExpression.sparse_rows(self._scipVarOffset)
for scipVars, coefficients, constant in rows:
scipExpression = Expr()
for scipVarIndex, coefficient in zip(scipVars, coefficients):
scipExpression += coefficient * self._scipVars[scipVarIndex]
scipExpression += constant
yield scipExpression
def _scalar_affinexp_pic2scip(self, picosExpression):
"""Transform a PICOS scalar affine expression into a SCIP expression.
:returns: A :class:`SCIP expression <pyscipopt.scip.Expr>`.
"""
assert len(picosExpression) == 1
return next(self._affinexp_pic2scip(picosExpression))
def _quadexp_pic2scip(self, picosExpression):
# Convert affine part.
scipExpression = self._scalar_affinexp_pic2scip(picosExpression.aff)
# Convert quadratic part.
for (pVar1, pVar2), pCoefficients \
in picosExpression._sparse_quads.items():
for sparseIndex in range(len(pCoefficients)):
localVar1Index = pCoefficients.I[sparseIndex]
localVar2Index = pCoefficients.J[sparseIndex]
localCoefficient = pCoefficients.V[sparseIndex]
scipVar1 = self._scipVars[
self._scipVarOffset[pVar1] + localVar1Index]
scipVar2 = self._scipVars[
self._scipVarOffset[pVar2] + localVar2Index]
scipExpression += localCoefficient * scipVar1 * scipVar2
return scipExpression
def _import_linear_constraint(self, picosConstraint):
assert isinstance(picosConstraint, AffineConstraint)
scipCons = []
picosLHS = picosConstraint.lmr
for scipLHS in self._affinexp_pic2scip(picosLHS):
if picosConstraint.is_increasing():
scipCons.append(self.int.addCons(scipLHS <= 0.0))
elif picosConstraint.is_decreasing():
scipCons.append(self.int.addCons(scipLHS >= 0.0))
elif picosConstraint.is_equality():
scipCons.append(self.int.addCons(scipLHS == 0.0))
else:
assert False, "Unexpected constraint relation."
return scipCons
def _import_quad_constraint(self, picosConstraint):
assert isinstance(picosConstraint,
(ConvexQuadraticConstraint, NonconvexQuadraticConstraint))
scipLHS = self._quadexp_pic2scip(picosConstraint.le0)
scipCons = [self.int.addCons(scipLHS <= 0.0)]
return scipCons
def _import_socone_constraint(self, picosConstraint):
scipCons = []
scipQuadLHS = self._quadexp_pic2scip(
(picosConstraint.ne | picosConstraint.ne)
- (picosConstraint.ub * picosConstraint.ub))
scipCons.append(self.int.addCons(scipQuadLHS <= 0.0))
scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub)
if scipAuxLHS.degree() > 0:
scipCons.append(self.int.addCons(scipAuxLHS >= 0.0))
return scipCons
def _import_rscone_constraint(self, picosConstraint):
scipCons = []
scipLHS = self._quadexp_pic2scip(
(picosConstraint.ne | picosConstraint.ne)
- (picosConstraint.ub1 * picosConstraint.ub2))
scipCons.append(self.int.addCons(scipLHS <= 0.0))
# Make sure that either the RHS is constant, or one of the two
# expressions is non-negative.
scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub1)
if scipAuxLHS.degree() > 0:
scipCons.append(self.int.addCons(scipAuxLHS >= 0.0))
else:
scipAuxLHS = self._scalar_affinexp_pic2scip(picosConstraint.ub2)
if scipAuxLHS.degree() > 0:
scipCons.append(self.int.addCons(scipAuxLHS >= 0.0))
return scipCons
def _import_constraint(self, picosConstraint):
# Import constraint based on type.
if isinstance(picosConstraint, AffineConstraint):
scipCons = self._import_linear_constraint(picosConstraint)
elif isinstance(picosConstraint,
(ConvexQuadraticConstraint, NonconvexQuadraticConstraint)):
scipCons = self._import_quad_constraint(picosConstraint)
elif isinstance(picosConstraint, SOCConstraint):
scipCons = self._import_socone_constraint(picosConstraint)
elif isinstance(picosConstraint, RSOCConstraint):
scipCons = self._import_rscone_constraint(picosConstraint)
else:
assert isinstance(picosConstraint, DummyConstraint), \
"Unexpected constraint type: {}".format(
picosConstraint.__class__.__name__)
# Map PICOS constraints to lists of SCIP constraints.
self._scipCons[picosConstraint] = scipCons
def _import_objective(self):
picosSense, picosObjective = self.ext.no
# Retrieve objective sense.
if picosSense == "min":
scipSense = "minimize"
else:
assert picosSense == "max"
scipSense = "maximize"
# Import objective function.
scipObjective = self._scalar_affinexp_pic2scip(picosObjective)
self.int.setObjective(scipObjective, scipSense)
def _import_problem(self):
import pyscipopt as scip
# Create a problem instance.
self.int = scip.Model()
# Import variables.
for variable in self.ext.variables.values():
self._import_variable(variable)
# Import constraints.
for constraint in self.ext.constraints.values():
self._import_constraint(constraint)
# Set objective.
self._import_objective()
def _update_problem(self):
# TODO: Support all problem updates supported by SCIP.
raise NotImplementedError
def _solve(self):
import pyscipopt as scip
# TODO: Give Problem an interface for checks like this.
isLP = isinstance(self.ext.no.function, AffineExpression) \
and all(isinstance(constraint, AffineConstraint)
for constraint in self.ext.constraints.values())
# Reset options.
self.int.resetParams()
# verbosity
picosVerbosity = self.verbosity()
if picosVerbosity <= -1:
scipVerbosity = 0
elif picosVerbosity == 0:
scipVerbosity = 2
elif picosVerbosity == 1:
scipVerbosity = 3
elif picosVerbosity >= 2:
scipVerbosity = 5
self.int.setIntParam("display/verblevel", scipVerbosity)
# TODO: "numerics/lpfeastol" has been replaced in SCIP 7.0.0 (PySCIPOpt
# 3.0.0) by "numerics/lpfeastolfactor", which is multiplied with
# some unknown default feasibility (maybe "numerics/feastol"?).
# Looking at it now it is not clear to me which of SCIP's
# feasibility tolerances does what exactly and whether they are
# absolute or relative measures. The workaround is to maintain
# PICOS' old guess but ignore abs_prim_fsb_tol with the new SCIP.
if int(scip.__version__.split(".")[0]) < 3:
# abs_prim_fsb_tol
if self.ext.options.abs_prim_fsb_tol is not None:
self.int.setRealParam("numerics/lpfeastol",
self.ext.options.abs_prim_fsb_tol)
# abs_dual_fsb_tol
if self.ext.options.abs_dual_fsb_tol is not None:
self.int.setRealParam("numerics/dualfeastol",
self.ext.options.abs_dual_fsb_tol)
# rel_prim_fsb_tol, rel_dual_fsb_tol
relFsbTols = [tol for tol in (self.ext.options.rel_prim_fsb_tol,
self.ext.options.rel_dual_fsb_tol) if tol is not None]
if relFsbTols:
self.int.setRealParam("numerics/feastol", min(relFsbTols))
# abs_ipm_opt_tol, abs_bnb_opt_tol
absOptTols = [tol for tol in (self.ext.options.abs_ipm_opt_tol,
self.ext.options.abs_bnb_opt_tol) if tol is not None]
if absOptTols:
self.int.setRealParam("limits/absgap", min(absOptTols))
# rel_ipm_opt_tol, rel_bnb_opt_tol
relOptTols = [tol for tol in (self.ext.options.rel_ipm_opt_tol,
self.ext.options.rel_bnb_opt_tol) if tol is not None]
if relOptTols:
self.int.setRealParam("limits/gap", min(relOptTols))
# timelimit
if self.ext.options.timelimit is not None:
self.int.setRealParam("limits/time", self.ext.options.timelimit)
# treememory
if self.ext.options.treememory is not None:
self.int.setRealParam("limits/memory", self.ext.options.treememory)
# max_fsb_nodes
if self.ext.options.max_fsb_nodes is not None:
self.int.setRealParam(
"limits/solutions", float(self.ext.options.max_fsb_nodes))
# Handle SCIP-specific options.
for key, value in self.ext.options.scip_params.items():
try:
if isinstance(value, bool):
self.int.setBoolParam(key, value)
elif isinstance(value, str):
if len(value) == 1:
try:
self.int.setCharParam(key, ord(value))
except LookupError:
self.int.setStringParam(key, value)
else:
self.int.setStringParam(key, value)
elif isinstance(value, float):
self.int.setRealParam(key, value)
elif isinstance(value, int):
try:
self.int.setIntParam(key, value)
except LookupError:
try:
self.int.setLongintParam(key, value)
except LookupError:
self.int.setRealParam(key, float(value))
else:
raise TypeError(
"The value '{}' has an uncommon type.".format(value))
except KeyError as error:
self._handle_bad_solver_specific_option_key(key, error)
except ValueError as error:
self._handle_bad_solver_specific_option_value(key, value, error)
except (TypeError, LookupError) as error:
picos_option = "{}_params".format(self.name)
assert picos_option in self.ext.options, \
"The PICOS option '{}' does not exist.".format(picos_option)
raise OptionValueError(
"Failed to guess type of {} option '{}' set via '{}'."
.format(self.short_name, key, picos_option)) from error
# Handle unsupported options.
self._handle_unsupported_options(
"hotstart", "lp_node_method", "lp_root_method", "max_iterations")
# In the case of a pure LP, disable presolve to get duals.
if self.ext.options.duals is not False and isLP:
self._debug("Disabling SCIP's presolve, heuristics, and propagation"
" features to allow for LP duals.")
self.int.setPresolve(scip.SCIP_PARAMSETTING.OFF)
self.int.setHeuristics(scip.SCIP_PARAMSETTING.OFF)
# Note that this is a helper to set options, so they get reset at
# the beginning of the function instead of in the else-scope below.
self.int.disablePropagation()
else:
self.int.setPresolve(scip.SCIP_PARAMSETTING.DEFAULT)
self.int.setHeuristics(scip.SCIP_PARAMSETTING.DEFAULT)
# Attempt to solve the problem.
with self._header(), self._stopwatch():
self.int.optimize()
# Retrieve primals.
primals = {}
if self.ext.options.primals is not False:
for picosVar in self.ext.variables.values():
scipVarOffset = self._scipVarOffset[picosVar]
value = []
for localIndex in range(picosVar.dim):
scipVar = self._scipVars[scipVarOffset + localIndex]
value.append(self.int.getVal(scipVar))
primals[picosVar] = value
# Retrieve duals for LPs.
duals = {}
if self.ext.options.duals is not False and isLP:
for picosConstraint in self.ext.constraints.values():
if isinstance(picosConstraint, DummyConstraint):
duals[picosConstraint] = cvxopt.spmatrix(
[], [], [], picosConstraint.size)
continue
assert isinstance(picosConstraint, AffineConstraint)
# Retrieve dual value for constraint.
scipDuals = []
for scipCon in self._scipCons[picosConstraint]:
scipDuals.append(self.int.getDualsolLinear(scipCon))
picosDual = cvxopt.matrix(scipDuals, picosConstraint.size)
# Flip sign based on constraint relation.
if not picosConstraint.is_increasing():
picosDual = -picosDual
# Flip sign based on objective sense.
if picosDual and self.ext.no.direction == "min":
picosDual = -picosDual
duals[picosConstraint] = picosDual
# Retrieve objective value.
value = self.int.getObjVal()
# Retrieve solution status.
status = self.int.getStatus()
if status == "optimal":
primalStatus = SS_OPTIMAL
dualStatus = SS_OPTIMAL
problemStatus = PS_FEASIBLE
elif status == "timelimit":
primalStatus = SS_PREMATURE
dualStatus = SS_PREMATURE
problemStatus = PS_UNKNOWN
elif status == "infeasible":
primalStatus = SS_INFEASIBLE
dualStatus = SS_UNKNOWN
problemStatus = PS_INFEASIBLE
elif status == "unbounded":
primalStatus = SS_UNKNOWN
dualStatus = SS_INFEASIBLE
problemStatus = PS_UNBOUNDED
else:
primalStatus = SS_UNKNOWN
dualStatus = SS_UNKNOWN
problemStatus = PS_UNKNOWN
return self._make_solution(
value, primals, duals, primalStatus, dualStatus, problemStatus)
# --------------------------------------
__all__ = api_end(_API_START, globals())