Implement MemorizingLazyConstrComponent

miplearn/collectors/basic.py

@@ -62,9 +62,7 @@ class BasicCollector:
                         and model.fix_violations is not None
                     ):
                         model.fix_violations(model, model.violations_, "aot")
-                        h5.put_scalar(
+                        h5.put_scalar("mip_constr_violations", repr(model.violations_))
-                            "mip_constr_violations", json.dumps(model.violations_)
-                        )
 
                     # Save MPS file
                     model.write(mps_filename)

miplearn/collectors/lazy.py

@@ -1,117 +0,0 @@
-#  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.
-
-from io import StringIO
-from typing import Callable
-
-import gurobipy as gp
-import numpy as np
-from gurobipy import GRB, LinExpr
-
-from ..h5 import H5File
-from ..io import _RedirectOutput
-
-
-class LazyCollector:
-    def __init__(
-        self,
-        min_constrs: int = 100_000,
-        time_limit: float = 900,
-    ) -> None:
-        self.min_constrs = min_constrs
-        self.time_limit = time_limit
-
-    def collect(
-        self, data_filename: str, build_model: Callable, tol: float = 1e-6
-    ) -> None:
-        h5_filename = f"{data_filename}.h5"
-        with H5File(h5_filename, "r+") as h5:
-            streams = [StringIO()]
-            lazy = None
-            with _RedirectOutput(streams):
-                slacks = h5.get_array("mip_constr_slacks")
-                assert slacks is not None
-
-                # Check minimum problem size
-                if len(slacks) < self.min_constrs:
-                    print("Problem is too small. Skipping.")
-                    h5.put_array("mip_constr_lazy", np.zeros(len(slacks)))
-                    return
-
-                # Load model
-                print("Loading model...")
-                model = build_model(data_filename)
-                model.params.LazyConstraints = True
-                model.params.timeLimit = self.time_limit
-                gp_constrs = np.array(model.getConstrs())
-                gp_vars = np.array(model.getVars())
-
-                # Load constraints
-                lhs = h5.get_sparse("static_constr_lhs")
-                rhs = h5.get_array("static_constr_rhs")
-                sense = h5.get_array("static_constr_sense")
-                assert lhs is not None
-                assert rhs is not None
-                assert sense is not None
-                lhs_csr = lhs.tocsr()
-                lhs_csc = lhs.tocsc()
-                constr_idx = np.array(range(len(rhs)))
-                lazy = np.zeros(len(rhs))
-
-                # Drop loose constraints
-                selected = (slacks > 0) & ((sense == b"<") | (sense == b">"))
-                loose_constrs = gp_constrs[selected]
-                print(
-                    f"Removing {len(loose_constrs):,d} constraints (out of {len(rhs):,d})..."
-                )
-                model.remove(list(loose_constrs))
-
-                # Filter to constraints that were dropped
-                lhs_csr = lhs_csr[selected, :]
-                lhs_csc = lhs_csc[selected, :]
-                rhs = rhs[selected]
-                sense = sense[selected]
-                constr_idx = constr_idx[selected]
-                lazy[selected] = 1
-
-                # Load warm start
-                var_names = h5.get_array("static_var_names")
-                var_values = h5.get_array("mip_var_values")
-                assert var_values is not None
-                assert var_names is not None
-                for (var_idx, var_name) in enumerate(var_names):
-                    var = model.getVarByName(var_name.decode())
-                    var.start = var_values[var_idx]
-
-                print("Solving MIP with lazy constraints callback...")
-
-                def callback(model: gp.Model, where: int) -> None:
-                    assert rhs is not None
-                    assert lazy is not None
-                    assert sense is not None
-
-                    if where == GRB.Callback.MIPSOL:
-                        x_val = np.array(model.cbGetSolution(model.getVars()))
-                        slack = lhs_csc * x_val - rhs
-                        slack[sense == b">"] *= -1
-                        is_violated = slack > tol
-
-                        for (j, rhs_j) in enumerate(rhs):
-                            if is_violated[j]:
-                                lazy[constr_idx[j]] = 0
-                                expr = LinExpr(
-                                    lhs_csr[j, :].data, gp_vars[lhs_csr[j, :].indices]
-                                )
-                                if sense[j] == b"<":
-                                    model.cbLazy(expr <= rhs_j)
-                                elif sense[j] == b">":
-                                    model.cbLazy(expr >= rhs_j)
-                                else:
-                                    raise RuntimeError(f"Unknown sense: {sense[j]}")
-
-                model.optimize(callback)
-                print(f"Marking {lazy.sum():,.0f} constraints as lazy...")
-
-            h5.put_array("mip_constr_lazy", lazy)
-            h5.put_scalar("mip_constr_lazy_log", streams[0].getvalue())

miplearn/components/lazy.py

@@ -1,43 +0,0 @@
-#  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 json
-from typing import Any, Dict, List
-
-import gurobipy as gp
-
-from ..h5 import H5File
-
-
-class ExpertLazyComponent:
-    def __init__(self) -> None:
-        pass
-
-    def fit(self, train_h5: List[str]) -> None:
-        pass
-
-    def before_mip(self, test_h5: str, model: gp.Model, stats: Dict[str, Any]) -> None:
-        with H5File(test_h5, "r") as h5:
-            constr_names = h5.get_array("static_constr_names")
-            constr_lazy = h5.get_array("mip_constr_lazy")
-            constr_violations = h5.get_scalar("mip_constr_violations")
-
-            assert constr_names is not None
-            assert constr_violations is not None
-
-            # Static lazy constraints
-            n_static_lazy = 0
-            if constr_lazy is not None:
-                for (constr_idx, constr_name) in enumerate(constr_names):
-                    if constr_lazy[constr_idx]:
-                        constr = model.getConstrByName(constr_name.decode())
-                        constr.lazy = 3
-                        n_static_lazy += 1
-            stats.update({"Static lazy constraints": n_static_lazy})
-
-            # Dynamic lazy constraints
-            if hasattr(model, "_fix_violations"):
-                violations = json.loads(constr_violations)
-                model._fix_violations(model, violations, "aot")
-                stats.update({"Dynamic lazy constraints": len(violations)})

miplearn/components/lazy/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
+
+import numpy as np
+from sklearn.preprocessing import MultiLabelBinarizer
+
+from miplearn.extractors.abstract import FeaturesExtractor
+from miplearn.h5 import H5File
+from miplearn.solvers.gurobi import GurobiModel
+
+logger = logging.getLogger(__name__)
+
+# TODO: Replace GurobiModel by AbstractModel
+# TODO: fix_violations: remove model.inner
+# TODO: fix_violations: remove `where`
+# TODO: Write documentation
+# TODO: Implement ExpertLazyConstrComponent
+
+
+class MemorizingLazyConstrComponent:
+    def __init__(self, clf: Any, extractor: FeaturesExtractor) -> None:
+        self.clf = clf
+        self.extractor = extractor
+        self.violations_: 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, violations, n_features = [], [], [], None
+        violation_to_idx: Dict[Hashable, int] = {}
+        for h5_filename in train_h5:
+            with H5File(h5_filename, "r") as h5:
+
+                # Store lazy constraints
+                sample_violations_str = h5.get_scalar("mip_constr_violations")
+                assert sample_violations_str is not None
+                assert isinstance(sample_violations_str, str)
+                sample_violations = eval(sample_violations_str)
+                assert isinstance(sample_violations, list)
+                y_sample = []
+                for v in sample_violations:
+                    if v not in violation_to_idx:
+                        violation_to_idx[v] = len(violation_to_idx)
+                        violations.append(v)
+                    y_sample.append(violation_to_idx[v])
+                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.violations_ = violations
+        self.n_targets_ = len(violation_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 before_mip(
+        self,
+        test_h5: str,
+        model: GurobiModel,
+        stats: Dict[str, Any],
+    ) -> None:
+        assert self.violations_ is not None
+        if model.fix_violations is None:
+            return
+
+        # 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.violations_[i] for (i, yi) in enumerate(y) if yi > 0.5]
+        logger.info(f"Enforcing {len(violations)} constraints ahead-of-time...")
+        model.fix_violations(model, violations, "aot")
+        stats["Lazy Constraints: AOT"] = len(violations)

miplearn/problems/tsp.py

@@ -150,12 +150,12 @@ def build_tsp_model(data: Union[str, TravelingSalesmanData]) -> GurobiModel:
         graph.add_edges_from(selected_edges)
         for component in list(nx.connected_components(graph)):
             if len(component) < model.inner._n_cities:
-                cut_edges = [
+                cut_edges = tuple(
-                    e
+                    (e[0], e[1])
                     for e in model.inner._edges
                     if (e[0] in component and e[1] not in component)
                     or (e[0] not in component and e[1] in component)
-                ]
+                )
                 violations.append(cut_edges)
         return violations
 

tests/components/lazy/test_mem.py

@@ -0,0 +1,62 @@
+#  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.
+
+from typing import List, Dict, Any
+from unittest.mock import Mock
+
+from sklearn.dummy import DummyClassifier
+from sklearn.neighbors import KNeighborsClassifier
+
+from miplearn.components.lazy.mem import MemorizingLazyConstrComponent
+from miplearn.extractors.abstract import FeaturesExtractor
+from miplearn.problems.tsp import build_tsp_model
+from miplearn.solvers.learning import LearningSolver
+
+
+def test_mem_component(
+    tsp_h5: List[str],
+    default_extractor: FeaturesExtractor,
+) -> None:
+    clf = Mock(wraps=DummyClassifier())
+    comp = MemorizingLazyConstrComponent(clf=clf, extractor=default_extractor)
+    comp.fit(tsp_h5)
+
+    # Should call fit method with correct arguments
+    clf.fit.assert_called()
+    x, y = clf.fit.call_args.args
+    assert x.shape == (3, 190)
+    assert y.tolist() == [
+        [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
+        [1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0],
+        [1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1],
+    ]
+
+    # Should store violations
+    assert comp.violations_ is not None
+    assert comp.n_features_ == 190
+    assert comp.n_targets_ == 22
+    assert len(comp.violations_) == 22
+
+    # Call before-mip
+    stats: Dict[str, Any] = {}
+    model = Mock()
+    comp.before_mip(tsp_h5[0], model, stats)
+
+    # Should call predict with correct args
+    clf.predict.assert_called()
+    (x_test,) = clf.predict.call_args.args
+    assert x_test.shape == (1, 190)
+
+
+def test_usage_tsp(
+    tsp_h5: List[str],
+    default_extractor: FeaturesExtractor,
+) -> None:
+    # Should not crash
+    data_filenames = [f.replace(".h5", ".pkl.gz") for f in tsp_h5]
+    clf = KNeighborsClassifier(n_neighbors=1)
+    comp = MemorizingLazyConstrComponent(clf=clf, extractor=default_extractor)
+    solver = LearningSolver(components=[comp])
+    solver.fit(data_filenames)
+    solver.optimize(data_filenames[0], build_tsp_model)

tests/components/primal/test_mem.py

@@ -20,7 +20,8 @@ logger = logging.getLogger(__name__)
 
 
 def test_mem_component(
-    multiknapsack_h5: List[str], default_extractor: FeaturesExtractor
+    multiknapsack_h5: List[str],
+    default_extractor: FeaturesExtractor,
 ) -> None:
     # Create mock classifier
     clf = Mock(wraps=DummyClassifier())

tests/conftest.py

@@ -8,13 +8,18 @@ from typing import List
 
 import pytest
 
-from miplearn.extractors.fields import H5FieldsExtractor
 from miplearn.extractors.abstract import FeaturesExtractor
+from miplearn.extractors.fields import H5FieldsExtractor
 
 
 @pytest.fixture()
 def multiknapsack_h5() -> List[str]:
-    return sorted(glob(f"{dirname(__file__)}/fixtures/multiknapsack*.h5"))
+    return sorted(glob(f"{dirname(__file__)}/fixtures/multiknapsack-n100*.h5"))
+
+
+@pytest.fixture()
+def tsp_h5() -> List[str]:
+    return sorted(glob(f"{dirname(__file__)}/fixtures/tsp-n20*.h5"))
 
 
 @pytest.fixture()

tests/fixtures/gen_tsp.py

@@ -0,0 +1,22 @@
+from os.path import dirname
+
+import numpy as np
+from scipy.stats import uniform, randint
+
+from miplearn.collectors.basic import BasicCollector
+from miplearn.io import write_pkl_gz
+from miplearn.problems.tsp import TravelingSalesmanGenerator, build_tsp_model
+
+np.random.seed(42)
+gen = TravelingSalesmanGenerator(
+    x=uniform(loc=0.0, scale=1000.0),
+    y=uniform(loc=0.0, scale=1000.0),
+    n=randint(low=20, high=21),
+    gamma=uniform(loc=1.0, scale=0.25),
+    fix_cities=True,
+    round=True,
+)
+data = gen.generate(3)
+data_filenames = write_pkl_gz(data, dirname(__file__), prefix="tsp-n20-")
+collector = BasicCollector()
+collector.collect(data_filenames, build_tsp_model)