# ------------------------------------------------------------------------------
# Copyright (C) 2018-2019 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 ..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 (wehich are sometimes incorrect).
# Nonconvex quadratics are supported as no other solver can deal
# with them at this point.
if 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"
[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._scipVar = dict()
"""Maps PICOS variable indices to SCIP variables."""
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._scipVar.clear()
self._scipCons.clear()
self._scipQuadObjAuxVar = None
self._scipQuadObjAuxCon = None
def _make_scip_var_names(self, picosVar, localIndex=None):
"""Make SCIP variable names.
Converts a PICOS variable to a list of SCIP variable names, each
corresponding to one scalar variable contained in the PICOS variable.
If localIndex is given, then only the name of the SCIP variable
representing the scalar variable with that offset is returned.
The name format is "picosName[localIndex]".
"""
# TODO: This function appears in multiple solvers, move it to the Solver
# base class as "_make_scalar_var_names".
if localIndex is not None:
return "{}[{}]".format(picosVar.name, localIndex)
else:
return [
self._make_scip_var_names(picosVar, localIndex)
for localIndex in range(len(picosVar))]
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]
# Give names.
scipNames = self._make_scip_var_names(picosVar)
# Import scalar variables.
for i in range(dim):
self._scipVar[picosVar.id_at(i)] = self.int.addVar(
name=scipNames[i], 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
for scipVars, coefficients, constant in picosExpression.sparse_rows(
None, indexFunction=lambda picosVar, localIndex:
self._scipVar[picosVar.id_at(localIndex)]):
scipExpression = Expr()
for scipVar, coefficient in zip(scipVars, coefficients):
scipExpression += coefficient * scipVar
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._scipVar[pVar1.id_at(localVar1Index)]
scipVar2 = self._scipVar[pVar2.id_at(localVar2Index)]
scipExpression += localCoefficient * scipVar1 * scipVar2
return scipExpression
def _import_linear_constraint(self, picosConstraint):
assert isinstance(picosConstraint, AffineConstraint)
scipCons = []
picosLHS = picosConstraint.lhs - picosConstraint.rhs
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.display_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():
value = []
for localIndex in range(picosVar.dim):
scipVar = self._scipVar[picosVar.id_at(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())