# coding: utf-8
# ------------------------------------------------------------------------------
# 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 ..compat import assert_import_exists
from ..constraints import (AffineConstraint, ConvexQuadraticConstraint,
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 Solver
_API_START = api_start(globals())
# -------------------------------
[docs]class SCIPSolver(Solver):
"""Interface to the SCIP solver via the PySCIPOpt Python interface."""
SUPPORTED = Specification.compile(
directions=Specification.ALL,
objectives=[
AffineExpression],
variables=Specification.ALL,
constraints=[
AffineConstraint,
SOCConstraint,
RSOCConstraint,
ConvexQuadraticConstraint,
NonconvexQuadraticConstraint])
[docs] @classmethod
def supports(cls, footprint):
"""Implement :meth:`~.solver.Solver.supports`."""
if not Solver.supports(footprint):
return False
# HACK: SCIP does not produce dual values (except for LPs, where they
# are known to be sometimes incorrect). For now pretend SCIP does
# not support continuous problems at all (unless there are
# nonconvex quadratics where no other solver works).
# TODO: Support all continuous problems when ``duals=False``.
if footprint.continuous \
and not ("con", NonconvexQuadraticConstraint) in footprint:
return False
return footprint << cls.SUPPORTED
[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`."""
assert_import_exists("pyscipopt")
[docs] @classmethod
def names(cls):
"""Implement :meth:`~.solver.Solver.names`."""
return "scip", "SCIP", "SCIP Optimization Suite"
[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.quadratic_forms.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 False, "Constraint type belongs to unsupported problem type."
# 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)
# HACK: Remove a constant term from the objective as this is not allowed
# by SCIP under Python 2. (Not sure if the term is stored or
# dropped under Python 3, but no exception is thrown there.)
from pyscipopt.scip import Term
constantTerm = Term()
if constantTerm in scipObjective.terms:
scipObjective.terms.pop(constantTerm)
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 option, value in self.ext.options.scip_params.items():
try:
if isinstance(value, bool):
self.int.setBoolParam(option, value)
elif isinstance(value, str):
if len(value) == 1:
try:
self.int.setCharParam(option, ord(value))
except LookupError:
self.int.setStringParam(option, value)
else:
self.int.setStringParam(option, value)
elif isinstance(value, float):
self.int.setRealParam(option, value)
elif isinstance(value, int):
try:
self.int.setIntParam(option, value)
except LookupError:
try:
self.int.setLongintParam(option, value)
except LookupError:
self.int.setRealParam(option, float(value))
except KeyError:
self._handle_bad_option_value("scip_params",
"SCIP option '{}' does not exist.".format(option))
except ValueError:
self._handle_bad_option_value("scip_params",
"Invalid value '{}' for SCIP option '{}'."
.format(value, option))
except LookupError:
self._handle_bad_option_value("scip_params",
"Failed to guess type of SCIP option '{}'.".format(option))
# 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():
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 picosConstraint.is_decreasing():
picosDual = -picosDual
# Flip sign based on objective sense.
if picosDual and self.ext.no.direction == "min":
picosDual = -picosDual
duals[picosConstraint] = picosDual
# # Retrieve objective value.
# objectiveValue = self.int.getObjVal()
#
# # HACK: Add back the constant objective term; see above.
# picosObjective = self.ext.no.function
# picosObjectiveAffinePart = picosObjective.aff \
# if isinstance(picosObjective, QuadExp) else picosObjective
# if picosObjectiveAffinePart.constant:
# objectiveValue += picosObjectiveAffinePart.constant[0]
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(
primals, duals, primalStatus, dualStatus, problemStatus)
# --------------------------------------
__all__ = api_end(_API_START, globals())