Implement MemorizingCutsComponent; STAB: switch to edge formulation

miplearn/collectors/basic.py

@@ -60,8 +60,7 @@ class BasicCollector:
 
                     # Add lazy constraints to model
                     if model.lazy_enforce is not None:
-                        model.lazy_enforce(model, model.lazy_constrs_)
-                        h5.put_scalar("mip_lazy", repr(model.lazy_constrs_))
+                        model.lazy_enforce(model, model.lazy_)
 
                     # Save MPS file
                     model.write(mps_filename)

miplearn/components/cuts/mem.py

@@ -0,0 +1,105 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+import logging
+from typing import List, Dict, Any, Hashable, Union
+
+import numpy as np
+from sklearn.preprocessing import MultiLabelBinarizer
+
+from miplearn.extractors.abstract import FeaturesExtractor
+from miplearn.h5 import H5File
+from miplearn.solvers.abstract import AbstractModel
+
+logger = logging.getLogger(__name__)
+
+
+class _BaseMemorizingConstrComponent:
+    def __init__(self, clf: Any, extractor: FeaturesExtractor, field: str) -> None:
+        self.clf = clf
+        self.extractor = extractor
+        self.constrs_: List[Hashable] = []
+        self.n_features_: int = 0
+        self.n_targets_: int = 0
+        self.field = field
+
+    def fit(
+        self,
+        train_h5: List[str],
+    ) -> None:
+        logger.info("Reading training data...")
+        n_samples = len(train_h5)
+        x, y, constrs, n_features = [], [], [], None
+        constr_to_idx: Dict[Hashable, int] = {}
+        for h5_filename in train_h5:
+            with H5File(h5_filename, "r") as h5:
+                # Store constraints
+                sample_constrs_str = h5.get_scalar(self.field)
+                assert sample_constrs_str is not None
+                assert isinstance(sample_constrs_str, str)
+                sample_constrs = eval(sample_constrs_str)
+                assert isinstance(sample_constrs, list)
+                y_sample = []
+                for c in sample_constrs:
+                    if c not in constr_to_idx:
+                        constr_to_idx[c] = len(constr_to_idx)
+                        constrs.append(c)
+                    y_sample.append(constr_to_idx[c])
+                y.append(y_sample)
+
+                # Extract features
+                x_sample = self.extractor.get_instance_features(h5)
+                assert len(x_sample.shape) == 1
+                if n_features is None:
+                    n_features = len(x_sample)
+                else:
+                    assert len(x_sample) == n_features
+                x.append(x_sample)
+        logger.info("Constructing matrices...")
+        assert n_features is not None
+        self.n_features_ = n_features
+        self.constrs_ = constrs
+        self.n_targets_ = len(constr_to_idx)
+        x_np = np.vstack(x)
+        assert x_np.shape == (n_samples, n_features)
+        y_np = MultiLabelBinarizer().fit_transform(y)
+        assert y_np.shape == (n_samples, self.n_targets_)
+        logger.info(
+            f"Dataset has {n_samples:,d} samples, "
+            f"{n_features:,d} features and {self.n_targets_:,d} targets"
+        )
+        logger.info("Training classifier...")
+        self.clf.fit(x_np, y_np)
+
+    def predict(
+        self,
+        msg: str,
+        test_h5: str,
+    ) -> List[Hashable]:
+        with H5File(test_h5, "r") as h5:
+            x_sample = self.extractor.get_instance_features(h5)
+        assert x_sample.shape == (self.n_features_,)
+        x_sample = x_sample.reshape(1, -1)
+        logger.info(msg)
+        y = self.clf.predict(x_sample)
+        assert y.shape == (1, self.n_targets_)
+        y = y.reshape(-1)
+        return [self.constrs_[i] for (i, yi) in enumerate(y) if yi > 0.5]
+
+
+class MemorizingCutsComponent(_BaseMemorizingConstrComponent):
+    def __init__(self, clf: Any, extractor: FeaturesExtractor) -> None:
+        super().__init__(clf, extractor, "mip_cuts")
+
+    def before_mip(
+        self,
+        test_h5: str,
+        model: AbstractModel,
+        stats: Dict[str, Any],
+    ) -> None:
+        if model.cuts_enforce is None:
+            return
+        assert self.constrs_ is not None
+        model.cuts_aot_ = self.predict("Predicting cutting planes...", test_h5)
+        stats["Cuts: AOT"] = len(model.cuts_aot_)

miplearn/components/lazy/mem.py

@@ -1,74 +1,22 @@
 #  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
 #  Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved.
 #  Released under the modified BSD license. See COPYING.md for more details.
+
 import logging
 from typing import List, Dict, Any, Hashable
 
-import numpy as np
-from sklearn.preprocessing import MultiLabelBinarizer
-
+from miplearn.components.cuts.mem import (
+    _BaseMemorizingConstrComponent,
+)
 from miplearn.extractors.abstract import FeaturesExtractor
-from miplearn.h5 import H5File
 from miplearn.solvers.abstract import AbstractModel
 
 logger = logging.getLogger(__name__)
 
 
-class MemorizingLazyConstrComponent:
+class MemorizingLazyComponent(_BaseMemorizingConstrComponent):
     def __init__(self, clf: Any, extractor: FeaturesExtractor) -> None:
-        self.clf = clf
-        self.extractor = extractor
-        self.constrs_: List[Hashable] = []
-        self.n_features_: int = 0
-        self.n_targets_: int = 0
-
-    def fit(self, train_h5: List[str]) -> None:
-        logger.info("Reading training data...")
-        n_samples = len(train_h5)
-        x, y, constrs, n_features = [], [], [], None
-        constr_to_idx: Dict[Hashable, int] = {}
-        for h5_filename in train_h5:
-            with H5File(h5_filename, "r") as h5:
-
-                # Store lazy constraints
-                sample_constrs_str = h5.get_scalar("mip_lazy")
-                assert sample_constrs_str is not None
-                assert isinstance(sample_constrs_str, str)
-                sample_constrs = eval(sample_constrs_str)
-                assert isinstance(sample_constrs, list)
-                y_sample = []
-                for c in sample_constrs:
-                    if c not in constr_to_idx:
-                        constr_to_idx[c] = len(constr_to_idx)
-                        constrs.append(c)
-                    y_sample.append(constr_to_idx[c])
-                y.append(y_sample)
-
-                # Extract features
-                x_sample = self.extractor.get_instance_features(h5)
-                assert len(x_sample.shape) == 1
-                if n_features is None:
-                    n_features = len(x_sample)
-                else:
-                    assert len(x_sample) == n_features
-                x.append(x_sample)
-
-        logger.info("Constructing matrices...")
-        assert n_features is not None
-        self.n_features_ = n_features
-        self.constrs_ = constrs
-        self.n_targets_ = len(constr_to_idx)
-        x_np = np.vstack(x)
-        assert x_np.shape == (n_samples, n_features)
-        y_np = MultiLabelBinarizer().fit_transform(y)
-        assert y_np.shape == (n_samples, self.n_targets_)
-        logger.info(
-            f"Dataset has {n_samples:,d} samples, "
-            f"{n_features:,d} features and {self.n_targets_:,d} targets"
-        )
-
-        logger.info("Training classifier...")
-        self.clf.fit(x_np, y_np)
+        super().__init__(clf, extractor, "mip_lazy")
 
     def before_mip(
         self,
@@ -78,23 +26,8 @@ class MemorizingLazyConstrComponent:
     ) -> None:
         if model.lazy_enforce is None:
             return
-
         assert self.constrs_ is not None
-
-        # Read features
-        with H5File(test_h5, "r") as h5:
-            x_sample = self.extractor.get_instance_features(h5)
-        assert x_sample.shape == (self.n_features_,)
-        x_sample = x_sample.reshape(1, -1)
-
-        # Predict violated constraints
-        logger.info("Predicting violated lazy constraints...")
-        y = self.clf.predict(x_sample)
-        assert y.shape == (1, self.n_targets_)
-        y = y.reshape(-1)
-
-        # Enforce constraints
-        violations = [self.constrs_[i] for (i, yi) in enumerate(y) if yi > 0.5]
+        violations = self.predict("Predicting violated lazy constraints...", test_h5)
         logger.info(f"Enforcing {len(violations)} constraints ahead-of-time...")
         model.lazy_enforce(model, violations)
         stats["Lazy Constraints: AOT"] = len(violations)

miplearn/problems/stab.py

@@ -1,14 +1,13 @@
 #  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
 #  Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved.
 #  Released under the modified BSD license. See COPYING.md for more details.
-
+import logging
 from dataclasses import dataclass
-from typing import List, Union
+from typing import List, Union, Any, Hashable
 
 import gurobipy as gp
 import networkx as nx
 import numpy as np
-import pyomo.environ as pe
 from gurobipy import GRB, quicksum
 from networkx import Graph
 from scipy.stats import uniform, randint
@@ -16,7 +15,8 @@ from scipy.stats.distributions import rv_frozen
 
 from miplearn.io import read_pkl_gz
 from miplearn.solvers.gurobi import GurobiModel
-from miplearn.solvers.pyomo import PyomoModel
+
+logger = logging.getLogger(__name__)
 
 
 @dataclass
@@ -82,35 +82,43 @@ class MaxWeightStableSetGenerator:
         return nx.generators.random_graphs.binomial_graph(self.n.rvs(), self.p.rvs())
 
 
-def build_stab_model_gurobipy(data: MaxWeightStableSetData) -> GurobiModel:
-    data = _read_stab_data(data)
-    model = gp.Model()
-    nodes = list(data.graph.nodes)
-    x = model.addVars(nodes, vtype=GRB.BINARY, name="x")
-    model.setObjective(quicksum(-data.weights[i] * x[i] for i in nodes))
-    for clique in nx.find_cliques(data.graph):
-        model.addConstr(quicksum(x[i] for i in clique) <= 1)
-    model.update()
-    return GurobiModel(model)
-
-
-def build_stab_model_pyomo(
-    data: MaxWeightStableSetData,
-    solver: str = "gurobi_persistent",
-) -> PyomoModel:
-    data = _read_stab_data(data)
-    model = pe.ConcreteModel()
-    nodes = pe.Set(initialize=list(data.graph.nodes))
-    model.x = pe.Var(nodes, domain=pe.Boolean, name="x")
-    model.obj = pe.Objective(expr=sum([-data.weights[i] * model.x[i] for i in nodes]))
-    model.clique_eqs = pe.ConstraintList()
-    for clique in nx.find_cliques(data.graph):
-        model.clique_eqs.add(expr=sum(model.x[i] for i in clique) <= 1)
-    return PyomoModel(model, solver)
-
-
-def _read_stab_data(data: Union[str, MaxWeightStableSetData]) -> MaxWeightStableSetData:
+def build_stab_model(data: MaxWeightStableSetData) -> GurobiModel:
     if isinstance(data, str):
         data = read_pkl_gz(data)
     assert isinstance(data, MaxWeightStableSetData)
-    return data
+
+    model = gp.Model()
+    nodes = list(data.graph.nodes)
+
+    # Variables and objective function
+    x = model.addVars(nodes, vtype=GRB.BINARY, name="x")
+    model.setObjective(quicksum(-data.weights[i] * x[i] for i in nodes))
+
+    # Edge inequalities
+    for (i1, i2) in data.graph.edges:
+        model.addConstr(x[i1] + x[i2] <= 1)
+
+    def cuts_separate(m: GurobiModel) -> List[Hashable]:
+        # Retrieve optimal fractional solution
+        x_val = m.inner.cbGetNodeRel(x)
+
+        # Check that we selected at most one vertex for each
+        # clique in the graph (sum <= 1)
+        violations: List[Hashable] = []
+        for clique in nx.find_cliques(data.graph):
+            if sum(x_val[i] for i in clique) > 1.0001:
+                violations.append(tuple(sorted(clique)))
+        return violations
+
+    def cuts_enforce(m: GurobiModel, violations: List[Any]) -> None:
+        logger.info(f"Adding {len(violations)} clique cuts...")
+        for clique in violations:
+            m.add_constr(quicksum(x[i] for i in clique) <= 1)
+
+    model.update()
+
+    return GurobiModel(
+        model,
+        cuts_separate=cuts_separate,
+        cuts_enforce=cuts_enforce,
+    )

miplearn/solvers/abstract.py

@@ -23,7 +23,11 @@ class AbstractModel(ABC):
     def __init__(self) -> None:
         self.lazy_enforce: Optional[Callable] = None
         self.lazy_separate: Optional[Callable] = None
-        self.lazy_constrs_: Optional[List[Any]] = None
+        self.lazy_: Optional[List[Any]] = None
+        self.cuts_enforce: Optional[Callable] = None
+        self.cuts_separate: Optional[Callable] = None
+        self.cuts_: Optional[List[Any]] = None
+        self.cuts_aot_: Optional[List[Any]] = None
         self.where = self.WHERE_DEFAULT
 
     @abstractmethod

miplearn/solvers/gurobi.py

@@ -1,6 +1,7 @@
 #  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
 #  Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved.
 #  Released under the modified BSD license. See COPYING.md for more details.
+import logging
 from typing import Dict, Optional, Callable, Any, List
 
 import gurobipy as gp
@@ -11,16 +12,40 @@ from scipy.sparse import lil_matrix
 from miplearn.h5 import H5File
 from miplearn.solvers.abstract import AbstractModel
 
+logger = logging.getLogger(__name__)
 
-def _gurobi_callback(model: AbstractModel, where: int) -> None:
-    assert model.lazy_separate is not None
-    assert model.lazy_enforce is not None
-    assert model.lazy_constrs_ is not None
-    if where == GRB.Callback.MIPSOL:
-        model.where = model.WHERE_LAZY
-        violations = model.lazy_separate(model)
-        model.lazy_constrs_.extend(violations)
-        model.lazy_enforce(model, violations)
+
+def _gurobi_callback(model: AbstractModel, gp_model: gp.Model, where: int) -> None:
+    # Lazy constraints
+    if model.lazy_separate is not None:
+        assert model.lazy_enforce is not None
+        assert model.lazy_ is not None
+        if where == GRB.Callback.MIPSOL:
+            model.where = model.WHERE_LAZY
+            violations = model.lazy_separate(model)
+            if len(violations) > 0:
+                model.lazy_.extend(violations)
+                model.lazy_enforce(model, violations)
+
+    # User cuts
+    if model.cuts_separate is not None:
+        assert model.cuts_enforce is not None
+        assert model.cuts_ is not None
+        if where == GRB.Callback.MIPNODE:
+            status = gp_model.cbGet(GRB.Callback.MIPNODE_STATUS)
+            if status == GRB.OPTIMAL:
+                model.where = model.WHERE_CUTS
+                if model.cuts_aot_ is not None:
+                    violations = model.cuts_aot_
+                    model.cuts_aot_ = None
+                    logger.info(f"Enforcing {len(violations)} cuts ahead-of-time...")
+                else:
+                    violations = model.cuts_separate(model)
+                if len(violations) > 0:
+                    model.cuts_.extend(violations)
+                    model.cuts_enforce(model, violations)
+
+    # Cleanup
     model.where = model.WHERE_DEFAULT
 
 
@@ -44,10 +69,14 @@ class GurobiModel(AbstractModel):
         inner: gp.Model,
         lazy_separate: Optional[Callable] = None,
         lazy_enforce: Optional[Callable] = None,
+        cuts_separate: Optional[Callable] = None,
+        cuts_enforce: Optional[Callable] = None,
     ) -> None:
         super().__init__()
         self.lazy_separate = lazy_separate
         self.lazy_enforce = lazy_enforce
+        self.cuts_separate = cuts_separate
+        self.cuts_enforce = cuts_enforce
         self.inner = inner
 
     def add_constrs(
@@ -125,6 +154,10 @@ class GurobiModel(AbstractModel):
         except AttributeError:
             pass
         self._extract_after_mip_solution_pool(h5)
+        if self.lazy_ is not None:
+            h5.put_scalar("mip_lazy", repr(self.lazy_))
+        if self.cuts_ is not None:
+            h5.put_scalar("mip_cuts", repr(self.cuts_))
 
     def fix_variables(
         self,
@@ -149,14 +182,22 @@ class GurobiModel(AbstractModel):
             stats["Fixed variables"] = n_fixed
 
     def optimize(self) -> None:
-        self.lazy_constrs_ = []
+        self.lazy_ = []
+        self.cuts_ = []
 
         def callback(_: gp.Model, where: int) -> None:
-            _gurobi_callback(self, where)
+            _gurobi_callback(self, self.inner, where)
 
+        # Required parameters for lazy constraints
         if self.lazy_enforce is not None:
             self.inner.setParam("PreCrush", 1)
             self.inner.setParam("LazyConstraints", 1)
+
+        # Required parameters for user cuts
+        if self.cuts_enforce is not None:
+            self.inner.setParam("PreCrush", 1)
+
+        if self.lazy_enforce is not None or self.cuts_enforce is not None:
             self.inner.optimize(callback)
         else:
             self.inner.optimize()

miplearn/solvers/pyomo.py

@@ -36,7 +36,6 @@ class PyomoModel(AbstractModel):
         self._is_warm_start_available = False
         self.lazy_separate = lazy_separate
         self.lazy_enforce = lazy_enforce
-        self.lazy_constrs_: Optional[List[Any]] = None
         if not hasattr(self.inner, "dual"):
             self.inner.dual = Suffix(direction=Suffix.IMPORT)
             self.inner.rc = Suffix(direction=Suffix.IMPORT)
@@ -131,15 +130,14 @@ class PyomoModel(AbstractModel):
                 self.solver.update_var(var)
 
     def optimize(self) -> None:
-        self.lazy_constrs_ = []
-
+        self.lazy_ = []
         if self.lazy_separate is not None:
             assert (
                 self.solver_name == "gurobi_persistent"
             ), "Callbacks are currently only supported on gurobi_persistent"
 
             def callback(_: Any, __: Any, where: int) -> None:
-                _gurobi_callback(self, where)
+                _gurobi_callback(self, self.solver, where)
 
             self.solver.set_gurobi_param("PreCrush", 1)
             self.solver.set_gurobi_param("LazyConstraints", 1)