diff --git a/CHANGELOG.md b/CHANGELOG.md index addf21992..1257d32b5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,10 @@ - Added enableDebugSol() and disableDebugSol() for controlling the debug solution mechanism if DEBUGSOL=true - Added getVarPseudocostScore() and getVarPseudocost() - Added getNBranchings() and getNBranchingsCurrentRun() +- Added isActive() which wraps SCIPvarIsActive() and test +- Added aggregateVars() and tests +- Added example shiftbound.py +- Added a tutorial in ./docs on the presolver plugin ### Fixed - Raised an error when an expression is used when a variable is required - Fixed some compile warnings diff --git a/docs/tutorials/presolver.rst b/docs/tutorials/presolver.rst new file mode 100644 index 000000000..aab9f319e --- /dev/null +++ b/docs/tutorials/presolver.rst @@ -0,0 +1,227 @@ +########### +Presolvers +########### + +For the following let us assume that a Model object is available, which is created as follows: + +.. code-block:: python + + from pyscipopt import Model, Presol, SCIP_RESULT, SCIP_PRESOLTIMING + + scip = Model() + +.. contents:: Contents +---------------------- + + +What is Presolving? +=================== + +Presolving simplifies a problem before the actual search starts. Typical +transformations include: + +- tightening bounds, +- removing redundant variables/constraints, +- aggregating variables, +- detecting infeasibility early. + +This can reduce numerical issues and simplify constraints and objective +expressions without changing the solution space. + + +The Presol Plugin Interface (Python) +==================================== + +A presolver in PySCIPOpt is a subclass of ``pyscipopt.Presol`` that implements the method: + +- ``presolexec(self, nrounds, presoltiming)`` + +and is registered on a ``pyscipopt.Model`` via +the class method ``pyscipopt.Model.includePresol``. + +Here is a high-level flow: + +1. Subclass ``MyPresolver`` and capture any parameters in ``__init__``. +2. Implement ``presolexec``: inspect variables, compute transformations, call SCIP aggregation APIs, and return a result code. +3. Register your presolver using ``includePresol`` with a priority, maximal rounds, and timing. +4. Solve the model, e.g. by calling ``presolve`` or ``optimize``. + + +A Minimal Skeleton +------------------ + +.. code-block:: python + + from pyscipopt import Presol, SCIP_RESULT + + class MyPresolver(Presol): + def __init__(self, someparam=123): + self.someparam = someparam + + def presolexec(self, nrounds, presoltiming): + scip = self.model + + # ... inspect model, change bounds, aggregate variables, etc. ... + + return {"result": SCIP_RESULT.SUCCESS} # or DIDNOTFIND, DIDNOTRUN, CUTOFF + + +Example: Writing a Custom Presolver +=================================== + +This tutorial shows how to write a presolver entirely in Python using +PySCIPOpt's ``Presol`` plugin interface. We will implement a small +presolver that shifts variable bounds from ``[a, b]`` to ``[0, b - a]`` +and optionally flips signs to reduce constant offsets. + +For educational purposes, we keep our example as close as possible to SCIP's implementation, which can be found `here `__. However, one may implement Boundshift differently as SCIP's logic does not translate perfectly to Python. To avoid any confusion with the already implemented version of Boundshift, we will call our custom presolver *Shiftbound*. + +A complete working example can be found in the directory: + +- ``examples/finished/shiftbound.py`` + + +Implementing Shiftbound +----------------------- + +Below we walk through the important parts to illustrate design decisions to translate the Boundshift presolver to PySCIPOpt. + +We want to provide parameters to control the presolver's behaviour: + +- ``maxshift``: maximum length of interval ``b - a`` we are willing to shift, +- ``flipping``: allow sign flips for better numerics, +- ``integer``: only shift integer-ranged variables if true. + +We will put these parameters into the ``__init__`` method to help us initialise the attributes of the presolver class. Then, in ``presolexec``, we implement the algorithm our custom presolver must follow. + +.. code-block:: python + + import math + from pyscipopt import SCIP_RESULT, Presol + + class ShiftboundPresolver(Presol): + def __init__(self, maxshift=float("inf"), flipping=True, integer=True): + self.maxshift = maxshift + self.flipping = flipping + self.integer = integer + + def presolexec(self, nrounds, presoltiming): + scip = self.model + + # Utility replacements for a few SCIP helpers which are not exposed to PySCIPOpt + # Emulate SCIP's absolute real value + def REALABS(x): return math.fabs(x) + + # Emulate SCIP's "is integral" using the model's epsilon value + def SCIPisIntegral(val): + return val - math.floor(val + scip.epsilon()) <= scip.epsilon() + + # Emulate adjusted bound rounding for integral variables + def SCIPadjustedVarBound(var, val): + if val < 0 and -val >= scip.infinity(): + return -scip.infinity() + if val > 0 and val >= scip.infinity(): + return scip.infinity() + if var.vtype() != "CONTINUOUS": + return scip.feasCeil(val) + if REALABS(val) <= scip.epsilon(): + return 0.0 + return val + + # Respect global presolve switches (here, if aggregation disabled) + if scip.getParam("presolving/donotaggr"): + return {"result": SCIP_RESULT.DIDNOTRUN} + + # We want to operate on non-binary active variables only + scipvars = scip.getVars() + nbin = scip.getNBinVars() + vars = scipvars[nbin:] # SCIP orders by type: binaries first + + result = SCIP_RESULT.DIDNOTFIND + + for var in reversed(vars): + if var.vtype() == "BINARY": + continue + if not var.isActive(): + continue + + lb = var.getLbGlobal() + ub = var.getUbGlobal() + + # For integral types: round to feasible integers to avoid noise + if var.vtype() != "CONTINUOUS": + assert SCIPisIntegral(lb) + assert SCIPisIntegral(ub) + lb = SCIPadjustedVarBound(var, lb) + ub = SCIPadjustedVarBound(var, ub) + + # Is the variable already fixed? + if scip.isEQ(lb, ub): + continue + + # If demanded by the parameters, restrict to integral-length intervals + if self.integer and not SCIPisIntegral(ub - lb): + continue + + # Only shift "reasonable" finite bounds + MAXABSBOUND = 1000.0 + shiftable = all(( + not scip.isEQ(lb, 0.0), + scip.isLT(ub, scip.infinity()), + scip.isGT(lb, -scip.infinity()), + scip.isLT(ub - lb, self.maxshift), + scip.isLE(REALABS(lb), MAXABSBOUND), + scip.isLE(REALABS(ub), MAXABSBOUND), + )) + if not shiftable: + continue + + # Create a new variable y with bounds [0, ub-lb], and same type + newvar = scip.addVar( + name=f"{var.name}_shift", + vtype=var.vtype(), + lb=0.0, + ub=(ub - lb), + obj=0.0, + ) + + # Aggregate old variable with new variable: + # x = y + lb (no flip), or + # x = -y + ub (flip), whichever yields smaller |offset| + if self.flipping and abs(ub) < abs(lb): + infeasible, redundant, aggregated = scip.aggregateVars(var, newvar, 1.0, 1.0, ub) + else: + infeasible, redundant, aggregated = scip.aggregateVars(var, newvar, 1.0, -1.0, lb) + + # Has the problem become infeasible? + if infeasible: + return {"result": SCIP_RESULT.CUTOFF} + + # Aggregation succeeded; SCIP marks x as redundant and keeps y for further search + assert redundant + assert aggregated + result = SCIP_RESULT.SUCCESS + + return {"result": result} + +Registering the Presolver +------------------------- + +After having initialised our ``model``, we instantiate an object based on our ``ShiftboundPresolver`` including the parameters we wish our presolver's behaviour to be set to. +Lastly, we register the custom presolver by including ``presolver``, followed by a name and a description, as well as specifying its priority, maximum rounds to be called (where ``-1`` specifies no limit), and timing mode. + +.. code-block:: python + + from pyscipopt import Model, SCIP_PRESOLTIMING, SCIP_PARAMSETTING + + model = Model() + + presolver = ShiftboundPresolver(maxshift=float("inf"), flipping=True, integer=True) + model.includePresol( + presolver, + "shiftbound", + "converts variables with domain [a,b] to variables with domain [0,b-a]", + priority=7900000, + maxrounds=-1, + timing=SCIP_PRESOLTIMING.FAST, + ) diff --git a/examples/finished/shiftbound.py b/examples/finished/shiftbound.py new file mode 100644 index 000000000..7a29e13fa --- /dev/null +++ b/examples/finished/shiftbound.py @@ -0,0 +1,292 @@ +""" +Example showing a custom presolver using PySCIPOpt's Presol plugin. + +This example reproduces the logic of boundshift.c from the SCIP source +as closely as possible. Some additional wrappers were required to expose +SCIP functionality to Python (e.g., scip.aggregateVars()). Other logic +(e.g., REALABS()) can be implemented in Python to mirror SCIP behaviour +without extra wrappers. This illustrates that users can implement +presolvers in Python by combining existing PySCIPOpt wrappers with +Python code and SCIP documentation. + +A simple knapsack problem was chosen to let the presolver plugin +operate on. +""" + +import math +from pyscipopt import ( + Model, + SCIP_PARAMSETTING, + SCIP_PRESOLTIMING, + Presol, + SCIP_RESULT +) +from typing import List, Optional + + +class ShiftboundPresolver(Presol): + """ + A presolver that converts variable domains from [a, b] to [0, b - a]. + + Attributes: + maxshift: float - Maximum absolute shift allowed. + flipping: bool - Whether to allow flipping (multiplying by -1) for + differentiation. + integer: bool - Whether to shift only integer ranges. + """ + + def __init__( + self, + maxshift: float = float("inf"), + flipping: bool = True, + integer: bool = True, + ): + self.maxshift = maxshift + self.flipping = flipping + self.integer = integer + + def presolexec(self, nrounds, presoltiming): + # the greatest absolute value by which bounds can be shifted to avoid + # large constant offsets + MAXABSBOUND = 1000.0 + + scip = self.model + + def REALABS(x): + return math.fabs(x) + + # scip.isIntegral() not implemented in wrapper. Work-around: + # compute integrality using epsilon from SCIP settings. + def SCIPisIntegral(val): + return val - math.floor(val + scip.epsilon()) <= scip.epsilon() + + # SCIPadjustedVarLb() not implemented in wrapper. Work-around: + # return the adjusted (i.e., rounded, if var is integral type) bound. + # Does not change the bounds of the variable. + def SCIPadjustedVarBound(var, val): + if val < 0 and -val >= scip.infinity(): + return -scip.infinity() + elif val > 0 and val >= scip.infinity(): + return scip.infinity() + elif var.vtype() != "CONTINUOUS": + return scip.feasCeil(val) + elif REALABS(val) <= scip.epsilon(): + return 0.0 + else: + return val + + # check whether aggregation of variables is not allowed + if scip.getParam("presolving/donotaggr"): + return {"result": SCIP_RESULT.DIDNOTRUN} + + scipvars = scip.getVars() + nbinvars = scip.getNBinVars() # number of binary variables + # infer number of non-binary variables + nvars = scip.getNVars() - nbinvars + + # if number of non-binary variables equals zero + if nvars == 0: + return {"result": SCIP_RESULT.DIDNOTRUN} + + # copy the non-binary variables into a separate list. + # this slice works because SCIP orders variables by type, starting with + # binary variables + vars = scipvars[nbinvars:] + + # loop over the non-binary variables + for var in reversed(vars): + # sanity check that variable is indeed not binary + assert var.vtype() != "BINARY" + + # do not shift non-active (fixed or (multi-)aggregated) variables + if not var.isActive(): + continue + + # get current variable's bounds + lb = var.getLbGlobal() + ub = var.getUbGlobal() + + # It can happen that integer variable bounds have not been + # propagated yet or contain small noise. This could result in an + # aggregation that might trigger assertions when updating bounds of + # aggregated variables (floating-point rounding errors). + # check if variable is integer + if var.vtype() != "CONTINUOUS": + # assert if bounds are integral + assert SCIPisIntegral(lb) + assert SCIPisIntegral(ub) + + # round the bound values for integral variables + lb = SCIPadjustedVarBound(var, lb) + ub = SCIPadjustedVarBound(var, ub) + + # sanity check lb < ub + assert scip.isLE(lb, ub) + # check if variable is already fixed + if scip.isEQ(lb, ub): + continue + # only operate on integer variables + if self.integer and not SCIPisIntegral(ub - lb): + continue + + # bounds are shiftable if all following conditions hold + cases = [ + not scip.isEQ(lb, 0.0), + scip.isLT(ub, scip.infinity()), + scip.isGT(lb, -scip.infinity()), + scip.isLT(ub - lb, self.maxshift), + scip.isLE(REALABS(lb), MAXABSBOUND), + scip.isLE(REALABS(ub), MAXABSBOUND), + ] + if all(cases): + # indicators for status of aggregation + infeasible = False + redundant = False + aggregated = False + + # create new variable with same properties as the current + # variable but with an added "_shift" suffix + orig_name = var.name + newvar = scip.addVar( + name=f"{orig_name}_shift", + vtype=f"{var.vtype()}", + lb=0.0, + ub=(ub - lb), + obj=0.0, + ) + + # aggregate old variable with new variable + # check if self.flipping is True + if self.flipping: + # check if |ub| < |lb| + if REALABS(ub) < REALABS(lb): + infeasible, redundant, aggregated = scip.aggregateVars( + var, newvar, 1.0, 1.0, ub + ) + else: + infeasible, redundant, aggregated = scip.aggregateVars( + var, newvar, 1.0, -1.0, lb + ) + else: + infeasible, redundant, aggregated = scip.aggregateVars( + var, newvar, 1.0, -1.0, lb + ) + + # problem has now become infeasible + if infeasible: + result = SCIP_RESULT.CUTOFF + else: + # sanity check flags + assert redundant + assert aggregated + + result = SCIP_RESULT.SUCCESS + + else: + result = SCIP_RESULT.DIDNOTFIND + + return {"result": result} + + +def knapsack( + instance_name: str, + sizes: List[int], + values: List[int], + upper_bound: List[int], + lower_bound: List[int], + capacity: int, + vtypes: Optional[List[str]], +) -> tuple[Model, dict]: + """ + Model an instance of the knapsack problem + + Parameters: + sizes: List[int] - the sizes of the items + values: List[int] - the values of the items + upper_bound: List[int] - upper bounds per variable + lower_bound: List[int] - lower bounds per variable + capacity: int - the knapsack capacity + vtypes: Optional[List[str]] - variable types ("B", "I", "C") + + Returns: + tuple[Model, dict] - the SCIP model and the variables dictionary + """ + + m = Model(f"Knapsack: {instance_name}") + x = {} + for i in range(len(sizes)): + assert isinstance(sizes[i], int) + assert isinstance(values[i], int) + assert isinstance(upper_bound[i], int) + assert isinstance(lower_bound[i], int) + + vt = "I" + if vtypes is not None: + assert len(vtypes) == len(sizes) + assert isinstance(vtypes[i], str) or (vtypes[i] is None) + vt = vtypes[i] + + x[i] = m.addVar( + vtype=vt, + obj=values[i], + lb=lower_bound[i], + ub=upper_bound[i], + name=f"x{i}", + ) + + assert isinstance(capacity, int) + m.addCons(sum(sizes[i] * x[i] for i in range(len(sizes))) <= capacity) + + m.setMaximize() + + return m, x + + +if __name__ == "__main__": + instance_name = "Knapsack" + sizes = [2, 1, 3] + values = [2, 3, 1] + upper_bounds = [1, 4, 1] + lower_bounds = [0, 2, 0] + capacity = 3 + + model, var_list = knapsack( + instance_name, sizes, values, upper_bounds, lower_bounds, capacity + ) + + + # isolate test: disable many automatic presolvers/propagators + model.setSeparating(SCIP_PARAMSETTING.OFF) + model.setHeuristics(SCIP_PARAMSETTING.OFF) + model.disablePropagation() + for key in ( + "presolving/boundshift/maxrounds", + "presolving/domcol/maxrounds", + "presolving/dualsparsify/maxrounds", + "presolving/implics/maxrounds", + "presolving/inttobinary/maxrounds", + "presolving/milp/maxrounds", + "presolving/sparsify/maxrounds", + "presolving/trivial/maxrounds", + "propagating/dualfix/maxprerounds", + "propagating/probing/maxprerounds", + "propagating/symmetry/maxprerounds", + "constraints/linear/maxprerounds", + ): + model.setParam(key, 0) + + # register and apply custom boundshift presolver + presolver = ShiftboundPresolver( + maxshift=float("inf"), flipping=True, integer=True + ) + model.includePresol( + presolver, + "shiftbound", + "converts variables with domain [a,b] to variables with domain [0,b-a]", + priority=7900000, + maxrounds=-1, + timing=SCIP_PRESOLTIMING.FAST, + ) + + # run presolve on instance + model.presolve() \ No newline at end of file diff --git a/src/pyscipopt/scip.pxd b/src/pyscipopt/scip.pxd index 444ea743f..659cc64de 100644 --- a/src/pyscipopt/scip.pxd +++ b/src/pyscipopt/scip.pxd @@ -816,6 +816,16 @@ cdef extern from "scip/scip.h": SCIP_Longint SCIPvarGetNBranchingsCurrentRun(SCIP_VAR* var, SCIP_BRANCHDIR dir) SCIP_Bool SCIPvarMayRoundUp(SCIP_VAR* var) SCIP_Bool SCIPvarMayRoundDown(SCIP_VAR* var) + SCIP_Bool SCIPvarIsActive(SCIP_VAR* var) + SCIP_RETCODE SCIPaggregateVars(SCIP* scip, + SCIP_VAR* varx, + SCIP_VAR* vary, + SCIP_Real scalarx, + SCIP_Real scalary, + SCIP_Real rhs, + SCIP_Bool* infeasible, + SCIP_Bool* redundant, + SCIP_Bool* aggregated) # LP Methods SCIP_RETCODE SCIPgetLPColsData(SCIP* scip, SCIP_COL*** cols, int* ncols) diff --git a/src/pyscipopt/scip.pxi b/src/pyscipopt/scip.pxi index c2603d2c3..39c089039 100644 --- a/src/pyscipopt/scip.pxi +++ b/src/pyscipopt/scip.pxi @@ -1727,6 +1727,16 @@ cdef class Variable(Expr): """ return SCIPvarIsDeletable(self.scip_var) + def isActive(self): + """ + Returns whether variable is an active (neither fixed nor aggregated) variable. + + Returns + ------- + boolean + """ + return SCIPvarIsActive(self.scip_var) + def getNLocksDown(self): """ Returns the number of locks for rounding down. @@ -4120,6 +4130,39 @@ cdef class Model: PY_SCIP_CALL(SCIPdelVar(self._scip, var.scip_var, &deleted)) return deleted + def aggregateVars(self, Variable varx, Variable vary, scalarx=1.0, scalary=-1.0, rhs=0.0): + """ + Aggregate varx and vary: calls SCIPaggregateVars and returns + (infeasible, redundant, aggregated) as Python bools. + + Parameters + ---------- + varx : Variable + vary : Variable + scalarx : float + scalary : float + rhs : float + + Returns + ------- + infeasible : bool + redundant : bool + aggregated : bool + """ + cdef SCIP_Bool infeasible + cdef SCIP_Bool redundant + cdef SCIP_Bool aggregated + PY_SCIP_CALL(SCIPaggregateVars(self._scip, + varx.scip_var, + vary.scip_var, + scalarx, + scalary, + rhs, + &infeasible, + &redundant, + &aggregated)) + return infeasible, redundant, aggregated + def tightenVarLb(self, Variable var, lb, force=False): """ Tighten the lower bound in preprocessing or current node, if the bound is tighter. diff --git a/tests/test_aggregate_vars.py b/tests/test_aggregate_vars.py new file mode 100644 index 000000000..b4069b27c --- /dev/null +++ b/tests/test_aggregate_vars.py @@ -0,0 +1,95 @@ +""" +Tests for Model.aggregateVars (wrapper around SCIPaggregateVars). +""" + +from pyscipopt import ( + Model, + Presol, + SCIP_PRESOLTIMING, + SCIP_RESULT, +) + + +class _AggPresol(Presol): + """ + Minimal presolver that aggregates two given variables and records the flags. + """ + + def __init__(self, varx, vary, scalarx, scalary, rhs): + self._args = (varx, vary, scalarx, scalary, rhs) + self.last = None # (infeasible, redundant, aggregated) + + def presolexec(self, nrounds, presoltiming): + x, y, ax, ay, rhs = self._args + infeas, redun, aggr = self.model.aggregateVars(x, y, ax, ay, rhs) + self.last = (bool(infeas), bool(redun), bool(aggr)) + # return SUCCESS to indicate presolver did work + return {"result": SCIP_RESULT.SUCCESS} + + +def _build_model_xy(vtype="C", lbx=0.0, ubx=10.0, lby=0.0, uby=10.0): + """ + Build a tiny model with two variables x and y. + """ + m = Model("agg-vars-test") + m.hideOutput() + x = m.addVar(name="x", vtype=vtype, lb=lbx, ub=ubx) + y = m.addVar(name="y", vtype=vtype, lb=lby, ub=uby) + # trivial objective to have a complete model + m.setMaximize() + return m, x, y + + +def test_aggregate_vars_success(): + """ + Aggregation succeeds for x - y = 0 on continuous variables with + compatible bounds, when called from a presolver. + """ + model, x, y = _build_model_xy( + vtype="C", lbx=0.0, ubx=10.0, lby=0.0, uby=10.0 + ) + + presol = _AggPresol(x, y, 1.0, -1.0, 0.0) + model.includePresol( + presol, + "agg-test", + "aggregate x and y", + priority=10**7, + maxrounds=1, + timing=SCIP_PRESOLTIMING.FAST, + ) + + model.presolve() + assert presol.last is not None + infeasible, redundant, aggregated = presol.last + + assert not infeasible + assert aggregated + # model should stay consistent + model.optimize() + + +def test_aggregate_vars_infeasible_binary_sum_exceeds_domain(): + """ + Aggregation detects infeasibility for x + y = 3 on binary variables, + since max(x + y) = 2 in {0,1} x {0,1}. + """ + model, x, y = _build_model_xy( + vtype="B", lbx=0.0, ubx=1.0, lby=0.0, uby=1.0 + ) + + presol = _AggPresol(x, y, 1.0, 1.0, 3.0) + model.includePresol( + presol, + "agg-infeas", + "aggregate x and y to infeasibility", + priority=10**7, + maxrounds=1, + timing=SCIP_PRESOLTIMING.FAST, + ) + + model.presolve() + assert presol.last is not None + infeasible, redundant, aggregated = presol.last + + assert infeasible \ No newline at end of file diff --git a/tests/test_vars.py b/tests/test_vars.py index 604c7870b..236426dd1 100644 --- a/tests/test_vars.py +++ b/tests/test_vars.py @@ -111,3 +111,14 @@ def test_getNBranchingsCurrentRun(): n_branchings += var.getNBranchingsCurrentRun(SCIP_BRANCHDIR.DOWNWARDS) assert n_branchings == m.getNNodes() - 1 + +def test_isActive(): + m = Model() + x = m.addVar(vtype='C', lb=0.0, ub=1.0) + # newly added variables should be active + assert x.isActive() + m.freeProb() + + # TODO lacks tests for cases when returned false due to + # - fixed (probably during probing) + # - aggregated