Implement UserCutsComponent; modify TravelingSalesmanInstance

src/python/miplearn/__init__.py

@@ -10,6 +10,7 @@ from .extractors import (SolutionExtractor,
 from .components.component import Component
 from .components.objective import ObjectiveValueComponent
 from .components.lazy import LazyConstraintsComponent
+from .components.cuts import UserCutsComponent
 from .components.primal import PrimalSolutionComponent
 from .components.branching import BranchPriorityComponent, BranchPriorityExtractor
 

src/python/miplearn/components/cuts.py

@@ -0,0 +1,86 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+from copy import deepcopy
+
+from miplearn.classifiers.counting import CountingClassifier
+from miplearn.components import classifier_evaluation_dict
+
+from .component import Component
+from ..extractors import *
+
+logger = logging.getLogger(__name__)
+
+
+class UserCutsComponent(Component):
+    """
+    A component that predicts which user cuts to enforce.
+    """
+    
+    def __init__(self,
+                 classifier=CountingClassifier(),
+                 threshold=0.05):
+        self.violations = set()
+        self.count = {}
+        self.n_samples = 0
+        self.threshold = threshold
+        self.classifier_prototype = classifier
+        self.classifiers = {}
+
+    def before_solve(self, solver, instance, model):
+        logger.info("Predicting violated user cuts...")
+        violations = self.predict(instance)
+        logger.info("Enforcing %d cuts..." % len(violations))
+        for v in violations:
+            cut = instance.build_user_cut(model, v)
+            solver.internal_solver.add_constraint(cut)
+
+    def after_solve(self, solver, instance, model, results):
+        pass
+                
+    def fit(self, training_instances):
+        logger.debug("Fitting...")
+        features = InstanceFeaturesExtractor().extract(training_instances)
+
+        self.classifiers = {}
+        violation_to_instance_idx = {}
+        for (idx, instance) in enumerate(training_instances):
+            for v in instance.found_violated_user_cuts:
+                if v not in self.classifiers:
+                    self.classifiers[v] = deepcopy(self.classifier_prototype)
+                    violation_to_instance_idx[v] = []
+                violation_to_instance_idx[v] += [idx]
+
+        for (v, classifier) in tqdm(self.classifiers.items(), desc="Fit (user cuts)"):
+            logger.debug("Training: %s" % (str(v)))
+            label = np.zeros(len(training_instances))
+            label[violation_to_instance_idx[v]] = 1.0
+            classifier.fit(features, label)
+
+    def predict(self, instance):
+        violations = []
+        features = InstanceFeaturesExtractor().extract([instance])
+        for (v, classifier) in self.classifiers.items():
+            proba = classifier.predict_proba(features)
+            if proba[0][1] > self.threshold:
+                violations += [v]
+        return violations
+
+    def evaluate(self, instances):
+        results = {}
+        all_violations = set()
+        for instance in instances:
+            all_violations |= set(instance.found_violated_user_cuts)
+        for idx in tqdm(range(len(instances)), desc="Evaluate (lazy)"):
+            instance = instances[idx]
+            condition_positive = set(instance.found_violated_user_cuts)
+            condition_negative = all_violations - condition_positive
+            pred_positive = set(self.predict(instance)) & all_violations
+            pred_negative = all_violations - pred_positive
+            tp = len(pred_positive & condition_positive)
+            tn = len(pred_negative & condition_negative)
+            fp = len(pred_positive & condition_negative)
+            fn = len(pred_negative & condition_positive)
+            results[idx] = classifier_evaluation_dict(tp, tn, fp, fn)
+        return results

src/python/miplearn/components/lazy.py

@@ -46,7 +46,7 @@ class LazyConstraintsComponent(Component):
         self.classifiers = {}
         violation_to_instance_idx = {}
         for (idx, instance) in enumerate(training_instances):
-            for v in instance.found_violations:
+            for v in instance.found_violated_lazy_constraints:
                 if v not in self.classifiers:
                     self.classifiers[v] = deepcopy(self.classifier_prototype)
                     violation_to_instance_idx[v] = []
@@ -71,10 +71,10 @@ class LazyConstraintsComponent(Component):
         results = {}
         all_violations = set()
         for instance in instances:
-            all_violations |= set(instance.found_violations)
+            all_violations |= set(instance.found_violated_lazy_constraints)
         for idx in tqdm(range(len(instances)), desc="Evaluate (lazy)"):
             instance = instances[idx]
-            condition_positive = set(instance.found_violations)
+            condition_positive = set(instance.found_violated_lazy_constraints)
             condition_negative = all_violations - condition_positive
             pred_positive = set(self.predict(instance)) & all_violations
             pred_negative = all_violations - pred_positive

src/python/miplearn/components/tests/test_cuts.py

@@ -0,0 +1,31 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+import numpy as np
+import pyomo.environ as pe
+
+from miplearn import Instance, GurobiSolver, LearningSolver
+from miplearn.problems.knapsack import ChallengeA
+
+
+class CutInstance(Instance):
+    def to_model(self):
+        model = pe.ConcreteModel()
+        model.x = x = pe.Var([0, 1], domain=pe.Binary)
+        model.OBJ = pe.Objective(expr=x[0] + x[1], sense=pe.maximize)
+        model.eq = pe.Constraint(expr=2 * x[0] + 2 * x[1] <= 3)
+        return model
+
+    def get_instance_features(self):
+        return np.zeros(0)
+
+    def get_variable_features(self, var, index):
+        return np.zeros(0)
+
+
+def test_cut():
+    challenge = ChallengeA()
+    gurobi = GurobiSolver()
+    solver = LearningSolver(solver=gurobi, time_limit=10)
+    solver.solve(challenge.training_instances[0])
+    # assert False

src/python/miplearn/components/tests/test_lazy.py

@@ -15,8 +15,8 @@ E = 0.1
 
 def test_lazy_fit():
     instances, models = get_training_instances_and_models()
-    instances[0].found_violations = ["a", "b"]
-    instances[1].found_violations = ["b", "c"]
+    instances[0].found_violated_lazy_constraints = ["a", "b"]
+    instances[1].found_violated_lazy_constraints = ["b", "c"]
     classifier = Mock(spec=Classifier)
     component = LazyConstraintsComponent(classifier=classifier)
 
@@ -78,6 +78,7 @@ def test_lazy_before():
     # Should ask internal solver to add generated constraint
     solver.internal_solver.add_constraint.assert_called_once_with("c1")
 
+
 def test_lazy_evaluate():
     instances, models = get_training_instances_and_models()
     component = LazyConstraintsComponent()
@@ -88,8 +89,8 @@ def test_lazy_evaluate():
     component.classifiers["b"].predict_proba = Mock(return_value=[[0.0, 1.0]])
     component.classifiers["c"].predict_proba = Mock(return_value=[[0.0, 1.0]])
     
-    instances[0].found_violations = ["a", "b", "c"]
-    instances[1].found_violations = ["b", "d"]
+    instances[0].found_violated_lazy_constraints = ["a", "b", "c"]
+    instances[1].found_violated_lazy_constraints = ["b", "d"]
     assert component.evaluate(instances) == {
         0: {
             "Accuracy": 0.75,

src/python/miplearn/instance.py

@@ -76,7 +76,7 @@ class Instance(ABC):
         """
         return "default"
 
-    def find_violations(self, model):
+    def find_violated_lazy_constraints(self, model):
         """
         Returns lazy constraint violations found for the current solution.
         
@@ -99,16 +99,22 @@ class Instance(ABC):
         """
         Returns a Pyomo constraint which fixes a given violation.
         
-        This method is typically called immediately after find_violations. The violation object
-        provided to this method is exactly the same object returned earlier by find_violations.
+        This method is typically called immediately after find_violated_lazy_constraints. The violation object
+        provided to this method is exactly the same object returned earlier by find_violated_lazy_constraints.
         After some training, LearningSolver may decide to proactively build some lazy constraints
         at the beginning of the optimization process, before a solution is even available. In this
         case, build_lazy_constraints will be called without a corresponding call to
-        find_violations.
+        find_violated_lazy_constraints.
         
         The implementation should not directly add the constraint to the model. The constraint
         will be added by LearningSolver after the method returns.
         
         For a concrete example, see TravelingSalesmanInstance.
         """
-        pass
+        pass
+
+    def find_violated_user_cuts(self, model):
+        return []
+
+    def build_user_cut(self, model, violation):
+        pass

src/python/miplearn/problems/tests/test_tsp.py

@@ -61,6 +61,8 @@ def test_subtour():
     for solver_name in ['gurobi', 'cplex']:
         solver = LearningSolver(solver=solver_name)
         solver.solve(instance)
+        assert hasattr(instance, "found_violated_lazy_constraints")
+        assert hasattr(instance, "found_violated_user_cuts")
         x = instance.solution["x"]
         assert x[0,1] == 1.0
         assert x[0,4] == 1.0
@@ -69,4 +71,4 @@ def test_subtour():
         assert x[3,5] == 1.0
         assert x[4,5] == 1.0
         solver.fit([instance])
-        solver.solve(instance)
+        solver.solve(instance)

src/python/miplearn/problems/tsp.py

@@ -151,7 +151,7 @@ class TravelingSalesmanInstance(Instance):
     def get_variable_category(self, var_name, index):
         return index
     
-    def find_violations(self, model):
+    def find_violated_lazy_constraints(self, model):
         selected_edges = [e for e in model.edges if model.x[e].value > 0.5]
         graph = nx.Graph()
         graph.add_edges_from(selected_edges)
@@ -167,3 +167,9 @@ class TravelingSalesmanInstance(Instance):
                      if (e[0] in component and e[1] not in component) or
                         (e[0] not in component and e[1] in component)]
         return model.eq_subtour.add(sum(model.x[e] for e in cut_edges) >= 2)
+
+    def find_violated_user_cuts(self, model):
+        return self.find_violated_lazy_constraints(model)
+
+    def build_user_cut(self, model, violation):
+        return self.build_lazy_constraint(model, violation)

src/python/miplearn/solvers/gurobi.py

@@ -49,26 +49,45 @@ class GurobiSolver(PyomoSolver):
         from gurobipy import GRB
 
         def cb(cb_model, cb_opt, cb_where):
-            if cb_where == GRB.Callback.MIPSOL:
-                cb_opt.cbGetSolution(self._all_vars)
-                logger.debug("Finding violated constraints...")
-                violations = self.instance.find_violations(cb_model)
-                self.instance.found_violations += violations
-                logger.debug("    %d violations found" % len(violations))
-                for v in violations:
-                    cut = self.instance.build_lazy_constraint(cb_model, v)
-                    cb_opt.cbLazy(cut)
+            try:
+                # User cuts
+                if cb_where == GRB.Callback.MIPNODE:
+                    logger.debug("Finding violated cutting planes...")
+                    cb_opt.cbGetNodeRel(self._all_vars)
+                    violations = self.instance.find_violated_user_cuts(cb_model)
+                    self.instance.found_violated_user_cuts += violations
+                    logger.debug("    %d found" % len(violations))
+                    for v in violations:
+                        cut = self.instance.build_user_cut(cb_model, v)
+                        cb_opt.cbCut(cut)
+
+                # Lazy constraints
+                if cb_where == GRB.Callback.MIPSOL:
+                    cb_opt.cbGetSolution(self._all_vars)
+                    logger.debug("Finding violated lazy constraints...")
+                    violations = self.instance.find_violated_lazy_constraints(cb_model)
+                    self.instance.found_violated_lazy_constraints += violations
+                    logger.debug("    %d found" % len(violations))
+                    for v in violations:
+                        cut = self.instance.build_lazy_constraint(cb_model, v)
+                        cb_opt.cbLazy(cut)
+            except Exception as e:
+                logger.error(e)
+
+        self._pyomo_solver.options["LazyConstraints"] = 1
+        self._pyomo_solver.options["PreCrush"] = 1
+        self._pyomo_solver.set_callback(cb)
+
+        self.instance.found_violated_lazy_constraints = []
+        self.instance.found_violated_user_cuts = []
 
-        if hasattr(self.instance, "find_violations"):
-            self._pyomo_solver.options["LazyConstraints"] = 1
-            self._pyomo_solver.set_callback(cb)
-            self.instance.found_violations = []
         streams = [StringIO()]
         if tee:
             streams += [sys.stdout]
         with RedirectOutput(streams):
             results = self._pyomo_solver.solve(tee=True,
                                                warmstart=self._is_warm_start_available)
+
         self._pyomo_solver.set_callback(None)
         log = streams[0].getvalue()
         return {

src/python/miplearn/solvers/learning.py

@@ -10,10 +10,10 @@ from p_tqdm import p_map
 
 from .cplex import CPLEXSolver
 from .gurobi import GurobiSolver
-from .internal import InternalSolver
 from .. import (ObjectiveValueComponent,
                 PrimalSolutionComponent,
-                LazyConstraintsComponent)
+                LazyConstraintsComponent,
+                UserCutsComponent)
 
 logger = logging.getLogger(__name__)
 
@@ -31,7 +31,8 @@ def _parallel_solve(instance_idx):
         "Results": results,
         "Solution": instance.solution,
         "LP solution": instance.lp_solution,
-        "Violations": instance.found_violations,
+        "Violated lazy constraints": instance.found_violated_lazy_constraints,
+        "Violated user cuts": instance.found_violated_user_cuts,
     }
 
 
@@ -67,6 +68,7 @@ class LearningSolver:
             self.add(ObjectiveValueComponent())
             self.add(PrimalSolutionComponent())
             self.add(LazyConstraintsComponent())
+            self.add(UserCutsComponent())
 
         assert self.mode in ["exact", "heuristic"]
         for component in self.components.values():
@@ -107,7 +109,7 @@ class LearningSolver:
             - instance.lower_bound
             - instance.upper_bound
             - instance.solution
-            - instance.found_violations
+            - instance.found_violated_lazy_constraints
             - instance.solver_log
         Additional solver components may set additional properties. Please
         see their documentation for more details.
@@ -190,7 +192,8 @@ class LearningSolver:
             instances[idx].lp_value = r["Results"]["LP value"]
             instances[idx].lower_bound = r["Results"]["Lower bound"]
             instances[idx].upper_bound = r["Results"]["Upper bound"]
-            instances[idx].found_violations = r["Violations"]
+            instances[idx].found_violated_lazy_constraints = r["Violated lazy constraints"]
+            instances[idx].found_violated_user_cuts = r["Violated user cuts"]
             instances[idx].solver_log = r["Results"]["Log"]
 
         return results

src/python/miplearn/solvers/pyomo.py

@@ -121,20 +121,19 @@ class PyomoSolver(InternalSolver):
         streams = [StringIO()]
         if tee:
             streams += [sys.stdout]
-        self.instance.found_violations = []
+        self.instance.found_violated_lazy_constraints = []
+        self.instance.found_violated_user_cuts = []
         while True:
             logger.debug("Solving MIP...")
             with RedirectOutput(streams):
                 results = self._pyomo_solver.solve(tee=True,
                                                    warmstart=self._is_warm_start_available)
             total_wallclock_time += results["Solver"][0]["Wallclock time"]
-            if not hasattr(self.instance, "find_violations"):
-                break
             logger.debug("Finding violated constraints...")
-            violations = self.instance.find_violations(self.model)
+            violations = self.instance.find_violated_lazy_constraints(self.model)
             if len(violations) == 0:
                 break
-            self.instance.found_violations += violations
+            self.instance.found_violated_lazy_constraints += violations
             logger.debug("    %d violations found" % len(violations))
             for v in violations:
                 cut = self.instance.build_lazy_constraint(self.model, v)

src/python/miplearn/solvers/tests/test_learning_solver.py

@@ -34,7 +34,8 @@ def test_learning_solver():
             assert round(instance.lp_solution["x"][2], 3) == 1.000
             assert round(instance.lp_solution["x"][3], 3) == 0.000
             assert round(instance.lp_value, 3) == 1287.923
-            assert instance.found_violations == []
+            assert instance.found_violated_lazy_constraints == []
+            assert instance.found_violated_user_cuts == []
             assert len(instance.solver_log) > 100
 
             solver.fit([instance])