feature/mus scip

cpmpy/solvers/scip.py

@@ -24,6 +24,7 @@
 """
 import warnings
 from typing import Optional
+import cpmpy as cp
 
 from .solver_interface import SolverInterface, SolverStatus, ExitStatus
 from ..exceptions import NotSupportedError
@@ -95,6 +96,8 @@ def __init__(self, cpm_model=None, subsolver=None):
 
         self.scip_model = scip.Model()
         self.scip_model.setParam("display/verblevel", 0)  # remove solver logs from output
+        self.scip_model.setParam("iis/silent", True)  # suppress native IIS statistics/logs
+        self.scip_model.hideOutput()
         self.objective_value_ = None
         super().__init__(name="scip", cpm_model=cpm_model)
 
@@ -294,8 +297,14 @@ def add(self, cpm_expr):
 
     __add__ = add
 
-    def _add_transformed_constraint(self, cpm_expr):
-        """Add already transformed CPMpy constraints to the solver. Some constraints are further transformed in this file, such as reified linear equality constraints `b -> ... == k` into `b -> ... >= k and b -> ... <= k`. In this case, we recursively call this function instead of `self.add`, which avoids both the full transformation pipeline overhead and also does not pollute `user_vars` with `b`."""
+    def _add_transformed_constraint(self, cpm_expr, name=None):
+        """Add already transformed CPMpy constraints to the solver and return the native SCIP constraint(s).
+
+        Some constraints are further transformed in this file, such as reified linear equality constraints
+        `b -> ... == k` into `b -> ... >= k and b -> ... <= k`. In this case, we recursively call this
+        function instead of `self.add`, which avoids both the full transformation pipeline overhead and also
+        does not pollute `user_vars` with `b`.
+        """
         if isinstance(cpm_expr, Comparison):
             lhs, rhs = cpm_expr.args
             lhs_is_operator = isinstance(lhs, Operator)
@@ -304,18 +313,18 @@ def _add_transformed_constraint(self, cpm_expr):
             if cpm_expr.name == '<=':
                 if (lhs_is_operator and lhs.name == "sum" and all(a.is_bool() and not isinstance(a, NegBoolView) for a in lhs.args)):
                     if rhs == 1:
-                        self.scip_model.addConsSOS1(self.solver_vars(lhs.args))
+                        return self.scip_model.addConsSOS1(self.solver_vars(lhs.args), name=name or "")
                     else:
-                        self.scip_model.addConsCardinality(self.solver_vars(lhs.args), int(rhs))
+                        return self.scip_model.addConsCardinality(self.solver_vars(lhs.args), int(rhs), name=name or "")
                 else:
                     sciplhs = self._make_numexpr(lhs)
-                    self.scip_model.addCons(sciplhs <= sciprhs)
+                    return self.scip_model.addCons(sciplhs <= sciprhs, name=name or "")
             elif cpm_expr.name == '>=':
                 sciplhs = self._make_numexpr(lhs)
-                self.scip_model.addCons(sciplhs >= sciprhs)
+                return self.scip_model.addCons(sciplhs >= sciprhs, name=name or "")
             elif cpm_expr.name == '==':
                 sciplhs = self._make_numexpr(lhs)
-                self.scip_model.addCons(sciplhs == sciprhs)
+                return self.scip_model.addCons(sciplhs == sciprhs, name=name or "")
             else:
                 raise NotImplementedError(
                     "Not a known supported scip comparison '{}' {}".format(cpm_expr.name, cpm_expr))
@@ -336,12 +345,14 @@ def _add_transformed_constraint(self, cpm_expr):
                 else:
                     scip_cons = lin_expr >= rhs
                 if isinstance(cond, NegBoolView):
-                    self.scip_model.addConsIndicator(scip_cons, binvar=self.solver_var(cond._bv), activeone=False)
+                    return self.scip_model.addConsIndicator(scip_cons, binvar=self.solver_var(cond._bv), activeone=False, name=name or "")
                 else:
-                    self.scip_model.addConsIndicator(scip_cons, binvar=self.solver_var(cond), activeone=True)
+                    return self.scip_model.addConsIndicator(scip_cons, binvar=self.solver_var(cond), activeone=True, name=name or "")
             elif sub_expr.name == "==":
-                self._add_transformed_constraint(cond.implies(lhs <= rhs))
-                self._add_transformed_constraint(cond.implies(lhs >= rhs))
+                return [
+                    self._add_transformed_constraint(cond.implies(lhs <= rhs), name=f"{name}_le" if name else None),
+                    self._add_transformed_constraint(cond.implies(lhs >= rhs), name=f"{name}_ge" if name else None),
+                ]
             else:
                 raise Exception(f"Unknown linear expression {sub_expr} name")
 
@@ -361,7 +372,7 @@ def _add_transformed_constraint(self, cpm_expr):
                         scip_args.append(self.solver_var(arg))
 
                 # post constraint (note: `addConsXor` is tested to work for empty lists)
-                self.scip_model.addConsXor(scip_args, rhsvar)
+                return self.scip_model.addConsXor(scip_args, rhsvar, name=name or "")
             else:
                 raise NotImplementedError(
                     f"SCIP does not translate global constraint '{cpm_expr.name}' natively; "
@@ -370,19 +381,125 @@ def _add_transformed_constraint(self, cpm_expr):
                 )
         elif isinstance(cpm_expr, BoolVal):
             if cpm_expr.args[0] is False:
-                self.scip_model.addConsXor([], True)  # easiest way to post False to SCIP (e.g. 0 <= -1 is not allowed, bv <= -1 requires adding a dummy variables, ...)
+                return self.scip_model.addConsXor([], True, name=name or "")  # easiest way to post False to SCIP (e.g. 0 <= -1 is not allowed, bv <= -1 requires adding a dummy variables, ...)
+            return None
 
         elif isinstance(cpm_expr, DirectConstraint):
             cpm_expr.callSolver(self, self.scip_model)
+            return None
 
 
         else:
             raise NotImplementedError(cpm_expr)
 
+    @staticmethod
+    def _flatten_scip_cons(scip_cons):
+        if scip_cons is None:
+            return []
+        if isinstance(scip_cons, list):
+            return [con for sub in scip_cons for con in CPM_scip._flatten_scip_cons(sub)]
+        return [scip_cons]
+
+    @staticmethod
+    def _posts_multiple_scip_constraints(cpm_expr):
+        return (
+            isinstance(cpm_expr, Operator)
+            and cpm_expr.name == "->"
+            and isinstance(cpm_expr.args[1], Comparison)
+            and cpm_expr.args[1].name == "=="
+        )
+
+    @classmethod
+    def mus_native(cls, soft, hard=[]):
+        """
+        Compute a MUS using SCIP's native IIS (Irreducible Infeasible Subsystem) algorithm.
+
+        SCIP minimizes over native constraints, while CPMpy expects a MUS over user-level
+        constraints. If a CPMpy soft constraint transforms to multiple SCIP constraints, we
+        represent it by one activation constraint `a >= 1` and post `a -> transformed_constraint`
+        as hard constraints. The returned IIS can then be mapped back through the single
+        activation constraint.
+        """
+        soft_cons = toplevel_list(soft, merge_and=False)
+        hard_cons = toplevel_list(hard, merge_and=False)
+
+        s = cls()
+        native_soft_names = []
+        captured_hard_cons = []
+
+        def post_transformed(cpm_con):
+            posted = []
+            for tf_con in s.transform(cpm_con):
+                posted.extend(s._flatten_scip_cons(
+                    s._add_transformed_constraint(tf_con)
+                ))
+            return posted
+
+        def post_hard(cpm_con):
+            posted = post_transformed(cpm_con)
+            for scip_con in posted:
+                s.native_model.captureCons(scip_con)
+                captured_hard_cons.append(scip_con)
+            return posted
+
+        # Captured constraints are kept hard by SCIP's IIS extractor.
+        for hard_con in hard_cons:
+            post_hard(hard_con)
+
+        for soft_con in soft_cons:
+            soft_con_tf = s.transform(soft_con)
+
+            if len(soft_con_tf) == 0:
+                native_soft_names.append([])
+                continue
+            elif len(soft_con_tf) == 1 and not s._posts_multiple_scip_constraints(soft_con_tf[0]):
+                # One CPMpy soft constraint maps to one transformed SCIP-level constraint:
+                # post it directly so the IIS can minimize it natively.
+                scip_cons = s._flatten_scip_cons(
+                    s._add_transformed_constraint(soft_con_tf[0])
+                )
+                native_soft_names.append([con.name for con in scip_cons])
+            else:
+                # One CPMpy soft constraint maps to a group of SCIP constraints. Guard the
+                # group by one activation variable, then make that activation the only soft
+                # native constraint.
+                assumption = cp.boolvar()
+                guarded = assumption.implies(cp.all(soft_con_tf))
+                post_hard(guarded)
+
+                scip_cons = s._flatten_scip_cons(
+                    s._add_transformed_constraint(assumption >= 1)
+                )
+                native_soft_names.append([con.name for con in scip_cons])
+
+        try:
+            # `generateIIS()` solves the model if needed and raises if the model is feasible.
+            try:
+                iis = s.native_model.generateIIS()
+            except Exception as e:
+                # Keep the user-facing contract consistent with other native MUS extractors.
+                status = s.native_model.getStatus()
+                if status not in ("infeasible", "inforunbd"):
+                    raise AssertionError("MUS: model must be UNSAT") from e
+                raise
+
+            subscip = iis.getSubscip()
+            iis_names = {con.name for con in subscip.getConss()}
+
+            candidate = [
+                soft_con
+                for soft_con, scip_names in zip(soft_cons, native_soft_names)
+                if any(name in iis_names for name in scip_names)
+            ]
+
+            return candidate
+        finally:
+            for scip_con in captured_hard_cons:
+                s.native_model.releaseCons(scip_con)
+
     def solveAll(self, display=None, time_limit=None, solution_limit=None, call_from_model=False, **kwargs):
         warnings.warn("Solution enumeration is not implemented in PySCIPOpt, defaulting to CPMpy's naive implementation")
         # Issues to track for future reference:
         # - https://github.com/scipopt/PySCIPOpt/issues/549 and
         # - https://github.com/scipopt/PySCIPOpt/issues/248
         return super().solveAll(display, time_limit, solution_limit, call_from_model, **kwargs)
-

cpmpy/tools/explain/mus.py

@@ -7,6 +7,7 @@
     - Native MUS for given solvers:
         - Exact: deletion-based MUS extraction
         - Gurobi: IIS-based MUS extraction
+        - SCIP: IIS-based MUS extraction
 """
 import warnings
 import numpy as np
@@ -69,7 +70,7 @@ def mus_native(soft, hard=[], solver="exact"):
 
     :param soft: soft constraints, list of expressions
     :param hard: hard constraints, optional, list of expressions
-    :param solver: which solver to use (`exact` or `gurobi`)
+    :param solver: which solver to use (`exact`, `gurobi`, `highs`, or `scip`)
     """
 
     # get solver class
@@ -346,4 +347,3 @@ def optimal_mus_naive(soft, hard=[], weights=None, solver="ortools", hs_solver="
         Naive implementation of `optimal_mus` without assumption variables and incremental solving
     """
     return ocus_naive(soft, hard, weights, meta_constraint=True, solver=solver, hs_solver=hs_solver)
-