Temporarily remove native solver callbacks; add iteration_cb

miplearn/__init__.py

@@ -24,7 +24,7 @@ from .instance import Instance
 from .solvers.pyomo.base import BasePyomoSolver
 from .solvers.pyomo.cplex import CplexPyomoSolver
 from .solvers.pyomo.gurobi import GurobiPyomoSolver
-from .solvers.guroby import GurobiSolver
+from .solvers.gurobi import GurobiSolver
 from .solvers.internal import InternalSolver
 from .solvers.learning import LearningSolver
 

miplearn/components/cuts.py

@@ -30,6 +30,7 @@ class UserCutsComponent(Component):
         self.classifiers = {}
 
     def before_solve(self, solver, instance, model):
+        instance.found_violated_user_cuts = []
         logger.info("Predicting violated user cuts...")
         violations = self.predict(instance)
         logger.info("Enforcing %d user cuts..." % len(violations))

miplearn/components/lazy.py

@@ -30,6 +30,7 @@ class LazyConstraintsComponent(Component):
         self.classifiers = {}
 
     def before_solve(self, solver, instance, model):
+        instance.found_violated_lazy_constraints = []
         logger.info("Predicting violated lazy constraints...")
         violations = self.predict(instance)
         logger.info("Enforcing %d lazy constraints..." % len(violations))

miplearn/problems/knapsack.py

@@ -269,7 +269,8 @@ class GurobiKnapsackInstance(KnapsackInstance):
         n = len(self.weights)
         x = model.addVars(n, vtype=GRB.BINARY, name="x")
         model.addConstr(gp.quicksum(x[i] * self.weights[i]
-                                    for i in range(n)) <= self.capacity)
+                                    for i in range(n)) <= self.capacity,
+                        "eq_capacity")
         model.setObjective(gp.quicksum(x[i] * self.prices[i]
                                        for i in range(n)), GRB.MAXIMIZE)
         return model

miplearn/problems/tests/test_tsp.py

@@ -23,52 +23,52 @@ def test_generator():
     assert np.std(d) > 0
 
     
-def test_instance():
-    n_cities = 4
-    distances = np.array([
-        [0., 1., 2., 1.],
-        [1., 0., 1., 2.],
-        [2., 1., 0., 1.],
-        [1., 2., 1., 0.],
-    ])
-    instance = TravelingSalesmanInstance(n_cities, distances)
-    for solver_name in ['gurobi', 'cplex']:
-        solver = LearningSolver(solver=solver_name)
-        solver.solve(instance)
-        x = instance.solution["x"]
-        assert x[0,1] == 1.0
-        assert x[0,2] == 0.0
-        assert x[0,3] == 1.0
-        assert x[1,2] == 1.0
-        assert x[1,3] == 0.0
-        assert x[2,3] == 1.0
-        assert instance.lower_bound == 4.0
-        assert instance.upper_bound == 4.0
-
-    
-def test_subtour():
-    n_cities = 6
-    cities = np.array([
-        [0., 0.],
-        [1., 0.],
-        [2., 0.],
-        [3., 0.],
-        [0., 1.],
-        [3., 1.],
-    ])
-    distances = squareform(pdist(cities))
-    instance = TravelingSalesmanInstance(n_cities, distances)
-    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
-        assert x[1,2] == 1.0
-        assert x[2,3] == 1.0
-        assert x[3,5] == 1.0
-        assert x[4,5] == 1.0
-        solver.fit([instance])
-        solver.solve(instance)
+# def test_instance():
+#     n_cities = 4
+#     distances = np.array([
+#         [0., 1., 2., 1.],
+#         [1., 0., 1., 2.],
+#         [2., 1., 0., 1.],
+#         [1., 2., 1., 0.],
+#     ])
+#     instance = TravelingSalesmanInstance(n_cities, distances)
+#     for solver_name in ['gurobi', 'cplex']:
+#         solver = LearningSolver(solver=solver_name)
+#         solver.solve(instance)
+#         x = instance.solution["x"]
+#         assert x[0,1] == 1.0
+#         assert x[0,2] == 0.0
+#         assert x[0,3] == 1.0
+#         assert x[1,2] == 1.0
+#         assert x[1,3] == 0.0
+#         assert x[2,3] == 1.0
+#         assert instance.lower_bound == 4.0
+#         assert instance.upper_bound == 4.0
+#
+#
+# def test_subtour():
+#     n_cities = 6
+#     cities = np.array([
+#         [0., 0.],
+#         [1., 0.],
+#         [2., 0.],
+#         [3., 0.],
+#         [0., 1.],
+#         [3., 1.],
+#     ])
+#     distances = squareform(pdist(cities))
+#     instance = TravelingSalesmanInstance(n_cities, distances)
+#     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
+#         assert x[1,2] == 1.0
+#         assert x[2,3] == 1.0
+#         assert x[3,5] == 1.0
+#         assert x[4,5] == 1.0
+#         solver.fit([instance])
+#         solver.solve(instance)

miplearn/solvers/gurobi.py

@@ -26,7 +26,6 @@ class GurobiSolver(InternalSolver):
         self.params = params
         self._all_vars = None
         self._bin_vars = None
-        self._varname_to_var = None
 
     def set_instance(self, instance, model=None):
         if model is None:
@@ -83,45 +82,30 @@ class GurobiSolver(InternalSolver):
             "Log": log
         }
 
-    def solve(self, tee=False):
-        self.instance.found_violated_lazy_constraints = []
-        self.instance.found_violated_user_cuts = []
+    def solve(self, tee=False, iteration_cb=None):
+        total_wallclock_time = 0
+        total_nodes = 0
         streams = [StringIO()]
         if tee:
             streams += [sys.stdout]
+        if iteration_cb is None:
+            iteration_cb = lambda : False
+        while True:
+            logger.debug("Solving MIP...")
+            with RedirectOutput(streams):
+                self.model.optimize()
+            total_wallclock_time += self.model.runtime
+            total_nodes += int(self.model.nodeCount)
+            should_repeat = iteration_cb()
+            if not should_repeat:
+                break
 
-        def cb(cb_model, cb_where):
-            try:
-                # User cuts
-                if cb_where == self.GRB.Callback.MIPNODE:
-                    logger.debug("Finding violated cutting planes...")
-                    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_model.cbCut(cut)
-
-                # Lazy constraints
-                if cb_where == self.GRB.Callback.MIPSOL:
-                    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_model.cbLazy(cut)
-            except Exception as e:
-                logger.error(e)
-
-        with RedirectOutput(streams):
-            self.model.optimize(cb)
         log = streams[0].getvalue()
         return {
             "Lower bound": self.model.objVal,
             "Upper bound": self.model.objBound,
-            "Wallclock time": self.model.runtime,
-            "Nodes": int(self.model.nodeCount),
+            "Wallclock time": total_wallclock_time,
+            "Nodes": total_nodes,
             "Sense": ("min" if self.model.modelSense == 1 else "max"),
             "Log": log,
             "Warm start value": self._extract_warm_start_value(log),
@@ -143,8 +127,13 @@ class GurobiSolver(InternalSolver):
                 variables[varname] += [idx]
         return variables
 
-    def add_constraint(self, constraint):
-        self.model.addConstr(constraint)
+    def add_constraint(self, constraint, name=""):
+        if type(constraint) is tuple:
+            lhs, sense, rhs, name = constraint
+            logger.debug(lhs, sense, rhs)
+            self.model.addConstr(lhs, sense, rhs, name)
+        else:
+            self.model.addConstr(constraint, name=name)
 
     def set_warm_start(self, solution):
         count_fixed, count_total = 0, 0
@@ -172,6 +161,31 @@ class GurobiSolver(InternalSolver):
                 var.lb = value
                 var.ub = value
 
+    def get_constraints_ids(self):
+        self.model.update()
+        return [c.ConstrName for c in self.model.getConstrs()]
+
+    def extract_constraint(self, cid):
+        constr = self.model.getConstrByName(cid)
+        cobj = (self.model.getRow(constr),
+                constr.sense,
+                constr.RHS,
+                constr.ConstrName)
+        self.model.remove(constr)
+        return cobj
+
+    def is_constraint_satisfied(self, cobj, tol=1e-5):
+        lhs, sense, rhs, name = cobj
+        lhs_value = lhs.getValue()
+        if sense == "<":
+            return lhs_value <= rhs + tol
+        elif sense == ">":
+            return lhs_value >= rhs - tol
+        elif sense == "=":
+            return abs(rhs - lhs_value) < tol
+        else:
+            raise Exception("Unknown sense: %s" % sense)
+
     def set_branching_priorities(self, priorities):
         logger.warning("set_branching_priorities not implemented")
 

miplearn/solvers/internal.py

@@ -119,13 +119,19 @@ class InternalSolver(ABC):
         pass
 
     @abstractmethod
-    def solve(self, tee=False):
+    def solve(self, tee=False, iteration_cb=None):
         """
         Solves the currently loaded instance. After this method finishes,
         the best solution found can be retrieved by calling `get_solution`.
 
         Parameters
         ----------
+        iteration_cb: function
+            By default, InternalSolver makes a single call to the native `solve`
+            method and returns the result. If an iteration callback is provided
+            instead, InternalSolver enters a loop, where `solve` and `iteration_cb`
+            are called alternatively. To stop the loop, `iteration_cb` should
+            return False. Any other result causes the solver to loop again.
         tee: bool
             If true, prints the solver log to the screen.
 
@@ -138,25 +144,25 @@ class InternalSolver(ABC):
         """
         pass
 
-    # @abstractmethod
-    def get_constraint_names(self):
+    @abstractmethod
+    def get_constraints_ids(self):
         """
-        Returns a list of strings, containing the name of each constraint in the
-        model.
+        Returns a list of ids, which uniquely identify each constraint in the model.
         """
         pass
 
-    # @abstractmethod
-    def extract_constraint(self, cname):
+    @abstractmethod
+    def extract_constraint(self, cid):
         """
-        Removes a given constraint from the model and returns an object `c` which
+        Removes a given constraint from the model and returns an object `cobj` which
         can be used to verify if the removed constraint is still satisfied by
-        the current solution, using `is_constraint_satisfied(c)`, and can potentially
-        be re-added to the model using `add_constraint(c)`.
+        the current solution, using `is_constraint_satisfied(cobj)`, and can potentially
+        be re-added to the model using `add_constraint(cobj)`.
         """
         pass
 
-    def is_constraint_satisfied(self, c):
+    @abstractmethod
+    def is_constraint_satisfied(self, cobj):
         pass
 
     @abstractmethod

miplearn/solvers/pyomo/base.py

@@ -9,7 +9,7 @@ import pyomo
 from abc import abstractmethod
 from io import StringIO
 from pyomo import environ as pe
-from pyomo.core import Var
+from pyomo.core import Var, Constraint
 
 from .. import RedirectOutput
 from ..internal import InternalSolver
@@ -32,6 +32,7 @@ class BasePyomoSolver(InternalSolver):
         self._pyomo_solver = None
         self._obj_sense = None
         self._varname_to_var = {}
+        self._cname_to_constr = {}
 
     def solve_lp(self, tee=False):
         for var in self._bin_vars:
@@ -93,23 +94,31 @@ class BasePyomoSolver(InternalSolver):
         self.instance = instance
         self.model = model
         self._pyomo_solver.set_instance(model)
+        self._update_obj()
+        self._update_vars()
+        self._update_constrs()
 
-        # Update objective sense
+    def _update_obj(self):
         self._obj_sense = "max"
         if self._pyomo_solver._objective.sense == pyomo.core.kernel.objective.minimize:
             self._obj_sense = "min"
 
-        # Update variables
+    def _update_vars(self):
         self._all_vars = []
         self._bin_vars = []
         self._varname_to_var = {}
-        for var in model.component_objects(Var):
+        for var in self.model.component_objects(Var):
             self._varname_to_var[var.name] = var
             for idx in var:
                 self._all_vars += [var[idx]]
                 if var[idx].domain == pyomo.core.base.set_types.Binary:
                     self._bin_vars += [var[idx]]
 
+    def _update_constrs(self):
+        self._cname_to_constr = {}
+        for constr in self.model.component_objects(Constraint):
+            self._cname_to_constr[constr.name] = constr
+
     def fix(self, solution):
         count_total, count_fixed = 0, 0
         for varname in solution:
@@ -126,12 +135,15 @@ class BasePyomoSolver(InternalSolver):
 
     def add_constraint(self, constraint):
         self._pyomo_solver.add_constraint(constraint)
+        self._update_constrs()
 
-    def solve(self, tee=False):
+    def solve(self, tee=False, iteration_cb=None):
         total_wallclock_time = 0
         streams = [StringIO()]
         if tee:
             streams += [sys.stdout]
+        if iteration_cb is None:
+            iteration_cb = lambda: False
         self.instance.found_violated_lazy_constraints = []
         self.instance.found_violated_user_cuts = []
         while True:
@@ -140,16 +152,9 @@ class BasePyomoSolver(InternalSolver):
                 results = self._pyomo_solver.solve(tee=True,
                                                    warmstart=self._is_warm_start_available)
             total_wallclock_time += results["Solver"][0]["Wallclock time"]
-            logger.debug("Finding violated constraints...")
-            violations = self.instance.find_violated_lazy_constraints(self.model)
-            if len(violations) == 0:
+            should_repeat = iteration_cb()
+            if not should_repeat:
                 break
-            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)
-                self.add_constraint(cut)
-
         log = streams[0].getvalue()
         return {
             "Lower bound": results["Problem"][0]["Lower bound"],
@@ -198,6 +203,15 @@ class BasePyomoSolver(InternalSolver):
         key = self._get_gap_tolerance_option_name()
         self._pyomo_solver.options[key] = gap_tolerance
 
+    def get_constraints_ids(self):
+        return list(self._cname_to_constr.keys())
+
+    def extract_constraint(self, cid):
+        raise Exception("Not implemented")
+
+    def is_constraint_satisfied(self, cobj):
+        raise Exception("Not implemented")
+
     @abstractmethod
     def _get_warm_start_regexp(self):
         pass

miplearn/solvers/pyomo/gurobi.py

@@ -16,89 +16,23 @@ logger = logging.getLogger(__name__)
 
 class GurobiPyomoSolver(BasePyomoSolver):
     def __init__(self,
-                 use_lazy_callbacks=True,
                  options=None):
         """
         Creates a new Gurobi solver, accessed through Pyomo.
 
         Parameters
         ----------
-        use_lazy_callbacks: bool
-            If true, lazy constraints will be enforced via lazy callbacks.
-            Otherwise, they will be enforced via a simple solve-check loop.
         options: dict
             Dictionary of options to pass to the Pyomo solver. For example,
             {"Threads": 4} to set the number of threads.
         """
         super().__init__()
-        self._use_lazy_callbacks = use_lazy_callbacks
         self._pyomo_solver = pe.SolverFactory('gurobi_persistent')
         self._pyomo_solver.options["Seed"] = randint(low=0, high=1000).rvs()
         if options is not None:
             for (key, value) in options.items():
                 self._pyomo_solver.options[key] = value
 
-    def solve(self, tee=False):
-        if self._use_lazy_callbacks:
-            return self._solve_with_callbacks(tee)
-        else:
-            return super().solve(tee)
-
-    def _solve_with_callbacks(self, tee):
-        from gurobipy import GRB
-
-        def cb(cb_model, cb_opt, cb_where):
-            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 = []
-
-        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 {
-            "Lower bound": results["Problem"][0]["Lower bound"],
-            "Upper bound": results["Problem"][0]["Upper bound"],
-            "Wallclock time": results["Solver"][0]["Wallclock time"],
-            "Nodes": self._extract_node_count(log),
-            "Sense": self._obj_sense,
-            "Log": log,
-            "Warm start value": self._extract_warm_start_value(log),
-        }
-
     def _extract_node_count(self, log):
         return max(1, int(self._pyomo_solver._solver_model.getAttr("NodeCount")))
 

miplearn/solvers/tests/test_internal_solver.py

@@ -6,10 +6,9 @@ import logging
 from io import StringIO
 
 import pyomo.environ as pe
-from miplearn import BasePyomoSolver
-from miplearn.problems.knapsack import ChallengeA
-from miplearn.solvers import RedirectOutput
 
+from miplearn import BasePyomoSolver, GurobiSolver
+from miplearn.solvers import RedirectOutput
 from . import _get_instance, _get_internal_solvers
 
 logger = logging.getLogger(__name__)
@@ -99,17 +98,59 @@ def test_internal_solver():
         assert solution["x"][2] == 1.0
         assert solution["x"][3] == 1.0
 
+        # Add a brand new constraint
         if isinstance(solver, BasePyomoSolver):
-            model.cut = pe.Constraint(expr=model.x[0] <= 0.5)
+            model.cut = pe.Constraint(expr=model.x[0] <= 0.0, name="cut")
             solver.add_constraint(model.cut)
-            solver.solve_lp()
-            assert model.x[0].value == 0.5
+        elif isinstance(solver, GurobiSolver):
+            x = model.getVarByName("x[0]")
+            solver.add_constraint(x <= 0.0, name="cut")
+        else:
+            raise Exception("Illegal state")
+
+        # New constraint should affect solution and should be listed in
+        # constraint ids
+        assert solver.get_constraints_ids() == ["eq_capacity", "cut"]
+        stats = solver.solve()
+        assert stats["Lower bound"] == 1030.0
+
+        if isinstance(solver, GurobiSolver):
+            # Extract new constraint
+            cobj = solver.extract_constraint("cut")
+
+            # New constraint should no longer affect solution and should no longer
+            # be listed in constraint ids
+            assert solver.get_constraints_ids() == ["eq_capacity"]
+            stats = solver.solve()
+            assert stats["Lower bound"] == 1183.0
+
+            # New constraint should not be satisfied by current solution
+            assert not solver.is_constraint_satisfied(cobj)
+
+            # Re-add constraint
+            solver.add_constraint(cobj)
+
+            # Constraint should affect solution again
+            assert solver.get_constraints_ids() == ["eq_capacity", "cut"]
+            stats = solver.solve()
+            assert stats["Lower bound"] == 1030.0
+
+            # New constraint should now be satisfied
+            assert solver.is_constraint_satisfied(cobj)
+
+
+def test_iteration_cb():
+    for solver_class in _get_internal_solvers():
+        logger.info("Solver: %s" % solver_class)
+        instance = _get_instance(solver_class)
+        solver = solver_class()
+        solver.set_instance(instance)
+        count = 0
 
+        def custom_iteration_cb():
+            nonlocal count
+            count += 1
+            return count < 5
 
-# def test_node_count():
-#     for solver in _get_internal_solvers():
-#         challenge = ChallengeA()
-#         solver.set_time_limit(1)
-#         solver.set_instance(challenge.test_instances[0])
-#         stats = solver.solve(tee=True)
-#         assert stats["Nodes"] > 1
+        solver.solve(iteration_cb=custom_iteration_cb)
+        assert count == 5