Implement lazy callbacks & two-phase gap

miplearn/components/component.py

@@ -24,3 +24,6 @@ class Component(ABC):
 
     def after_iteration(self, solver, instance, model):
         return False
+
+    def on_lazy_callback(self, solver, instance, model):
+        return

miplearn/components/lazy_static.py

@@ -21,14 +21,29 @@ class LazyConstraint:
 class StaticLazyConstraintsComponent(Component):
     def __init__(self,
                  classifier=CountingClassifier(),
-                 threshold=0.05):
+                 threshold=0.05,
+                 use_two_phase_gap=True,
+                 large_gap=1e-2,
+                 violation_tolerance=-0.5,
+                ):
         self.threshold = threshold
         self.classifier_prototype = classifier
         self.classifiers = {}
         self.pool = []
+        self.original_gap = None
+        self.large_gap = large_gap
+        self.is_gap_large = False
+        self.use_two_phase_gap = use_two_phase_gap
+        self.violation_tolerance = violation_tolerance
 
     def before_solve(self, solver, instance, model):
         self.pool = []
+        if not solver.use_lazy_cb and self.use_two_phase_gap:
+            logger.info("Increasing gap tolerance to %f", self.large_gap)
+            self.original_gap = solver.gap_tolerance
+            self.is_gap_large = True
+            solver.internal_solver.set_gap_tolerance(self.large_gap)
+
         instance.found_violated_lazy_constraints = []
         if instance.has_static_lazy_constraints():
             self._extract_and_predict_static(solver, instance)
@@ -37,21 +52,39 @@ class StaticLazyConstraintsComponent(Component):
         pass
 
     def after_iteration(self, solver, instance, model):
-        logger.info("Finding violated lazy constraints...")
+        if solver.use_lazy_cb:
+            return False
+        else:
+            should_repeat = self._check_and_add(instance, solver)
+            if should_repeat:
+                return True
+            else:
+                if self.is_gap_large:
+                    logger.info("Restoring gap tolerance to %f", self.original_gap)
+                    solver.internal_solver.set_gap_tolerance(self.original_gap)
+                    self.is_gap_large = False
+                    return True
+                else:
+                    return False
+
+    def on_lazy_callback(self, solver, instance, model):
+        self._check_and_add(instance, solver)
+
+    def _check_and_add(self, instance, solver):
+        logger.debug("Finding violated lazy constraints...")
         constraints_to_add = []
         for c in self.pool:
-            if not solver.internal_solver.is_constraint_satisfied(c.obj):
+            if not solver.internal_solver.is_constraint_satisfied(c.obj,
+                                                                  tol=self.violation_tolerance):
                 constraints_to_add.append(c)
         for c in constraints_to_add:
             self.pool.remove(c)
             solver.internal_solver.add_constraint(c.obj)
             instance.found_violated_lazy_constraints += [c.cid]
         if len(constraints_to_add) > 0:
-            logger.info("Added %d lazy constraints back into the model" % len(constraints_to_add))
-            logger.info("Lazy constraint pool has %d constraints" % len(self.pool))
+            logger.info("%8d lazy constraints added %8d in the pool" % (len(constraints_to_add), len(self.pool)))
             return True
         else:
-            logger.info("Found no violated lazy constraints")
             return False
 
     def fit(self, training_instances):
@@ -92,7 +125,7 @@ class StaticLazyConstraintsComponent(Component):
                                    obj=solver.internal_solver.extract_constraint(cid))
                 constraints[category] += [c]
                 self.pool.append(c)
-        logger.info("Extracted %d lazy constraints" % len(self.pool))
+        logger.info("%8d lazy constraints extracted" % len(self.pool))
         logger.info("Predicting required lazy constraints...")
         n_added = 0
         for (category, x_values) in x.items():
@@ -108,8 +141,7 @@ class StaticLazyConstraintsComponent(Component):
                     self.pool.remove(c)
                     solver.internal_solver.add_constraint(c.obj)
                     instance.found_violated_lazy_constraints += [c.cid]
-        logger.info("Added %d lazy constraints back into the model" % n_added)
-        logger.info("Lazy constraint pool has %d constraints" % len(self.pool))
+        logger.info("%8d lazy constraints added %8d in the pool" % (n_added, len(self.pool)))
 
     def _collect_constraints(self, train_instances):
         constraints = {}

miplearn/components/tests/test_lazy_static.py

@@ -13,6 +13,9 @@ from miplearn.classifiers import Classifier
 
 def test_usage_with_solver():
     solver = Mock(spec=LearningSolver)
+    solver.use_lazy_cb = False
+    solver.gap_tolerance = 1e-4
+
     internal = solver.internal_solver = Mock(spec=InternalSolver)
     internal.get_constraint_ids = Mock(return_value=["c1", "c2", "c3", "c4"])
     internal.extract_constraint = Mock(side_effect=lambda cid: "<%s>" % cid)
@@ -37,7 +40,8 @@ def test_usage_with_solver():
         "c4": "type-b",
     }[cid])
 
-    component = StaticLazyConstraintsComponent(threshold=0.90)
+    component = StaticLazyConstraintsComponent(threshold=0.90,
+                                               use_two_phase_gap=False)
     component.classifiers = {
         "type-a": Mock(spec=Classifier),
         "type-b": Mock(spec=Classifier),

miplearn/log.py

@@ -7,40 +7,43 @@ import logging
 import time
 import sys
 
-if sys.stdout.isatty():
-    log_colors = {
-        "green":  '\033[92m',
-        "yellow": '\033[93m',
-        "red":    '\033[91m',
-        "reset":  '\033[0m',
-    }
-else:
-    log_colors = {
-        "green": "",
-        "yellow": "",
-        "red": "",
-        "reset": "",
-    }
-
 class TimeFormatter():
-    def __init__(self, start_time):
+    def __init__(self, start_time, log_colors):
         self.start_time = start_time
+        self.log_colors = log_colors
 
     def format(self, record):
         if record.levelno >= logging.ERROR:
-            color = log_colors["red"]
+            color = self.log_colors["red"]
         elif record.levelno >= logging.WARNING:
-            color = log_colors["yellow"]
+            color = self.log_colors["yellow"]
         else:
-            color = log_colors["green"]
+            color = self.log_colors["green"]
         return "%s[%12.3f]%s %s" % (color,
                                     record.created - self.start_time,
-                                    log_colors["reset"],
+                                    self.log_colors["reset"],
                                     record.getMessage())
 
-def setup_logger(start_time):
+def setup_logger(start_time=None,
+                 force_color=False):
+    if start_time is None:
+        start_time = time.time()
+    if sys.stdout.isatty() or force_color:
+        log_colors = {
+            "green":  '\033[92m',
+            "yellow": '\033[93m',
+            "red":    '\033[91m',
+            "reset":  '\033[0m',
+        }
+    else:
+        log_colors = {
+            "green": "",
+            "yellow": "",
+            "red": "",
+            "reset": "",
+        }        
     handler = logging.StreamHandler()
-    handler.setFormatter(TimeFormatter(start_time))
+    handler.setFormatter(TimeFormatter(start_time, log_colors))
     logging.getLogger().addHandler(handler)
     logging.getLogger("miplearn").setLevel(logging.INFO)
     lg = logging.getLogger("miplearn")

miplearn/solvers/gurobi.py

@@ -13,13 +13,25 @@ logger = logging.getLogger(__name__)
 
 
 class GurobiSolver(InternalSolver):
-    def __init__(self, params=None):
+    def __init__(self,
+                 params=None,
+                 lazy_cb_frequency=1,
+                 ):
+        """
+        An InternalSolver backed by Gurobi's Python API (without Pyomo).
+
+        Parameters
+        ----------
+        params
+            Parameters to pass to Gurobi. For example, params={"MIPGap": 1e-3}
+            sets the gap tolerance to 1e-3.
+        lazy_cb_frequency
+            If 1, calls lazy constraint callbacks whenever an integer solution
+            is found. If 2, calls it also at every node, after solving the
+            LP relaxation of that node.
+        """
         if params is None:
             params = {}
-            # params = {
-            #     "LazyConstraints": 1,
-            #     "PreCrush": 1,
-            # }
         from gurobipy import GRB
         self.GRB = GRB
         self.instance = None
@@ -27,8 +39,16 @@ class GurobiSolver(InternalSolver):
         self.params = params
         self._all_vars = None
         self._bin_vars = None
+        self.cb_where = None
+        assert lazy_cb_frequency in [1, 2]
+        if lazy_cb_frequency == 1:
+            self.lazy_cb_where = [self.GRB.Callback.MIPSOL]
+        else:
+            self.lazy_cb_where = [self.GRB.Callback.MIPSOL,
+                                  self.GRB.Callback.MIPNODE]
 
     def set_instance(self, instance, model=None):
+        self._raise_if_callback()
         if model is None:
             model = instance.to_model()
         self.instance = instance
@@ -36,6 +56,10 @@ class GurobiSolver(InternalSolver):
         self.model.update()
         self._update_vars()
 
+    def _raise_if_callback(self):
+        if self.cb_where is not None:
+            raise Exception("method cannot be called from a callback")
+
     def _update_vars(self):
         self._all_vars = {}
         self._bin_vars = {}
@@ -63,6 +87,7 @@ class GurobiSolver(InternalSolver):
             self.model.setParam(name, value)
 
     def solve_lp(self, tee=False):
+        self._raise_if_callback()
         self._apply_params()
         streams = [StringIO()]
         if tee:
@@ -83,7 +108,24 @@ class GurobiSolver(InternalSolver):
             "Log": log
         }
 
-    def solve(self, tee=False, iteration_cb=None):
+    def solve(self,
+              tee=False,
+              iteration_cb=None,
+              lazy_cb=None):
+        self._raise_if_callback()
+
+        def cb_wrapper(cb_model, cb_where):
+            try:
+                self.cb_where = cb_where
+                if cb_where in self.lazy_cb_where:
+                    lazy_cb(self, self.model)
+            except:
+                logger.exception("callback error")
+            finally:
+                self.cb_where = None
+
+        if lazy_cb:
+            self.params["LazyConstraints"] = 1
         self._apply_params()
         total_wallclock_time = 0
         total_nodes = 0
@@ -95,7 +137,10 @@ class GurobiSolver(InternalSolver):
         while True:
             logger.debug("Solving MIP...")
             with RedirectOutput(streams):
-                self.model.optimize()
+                if lazy_cb is None:
+                    self.model.optimize()
+                else:
+                    self.model.optimize(cb_wrapper)
             total_wallclock_time += self.model.runtime
             total_nodes += int(self.model.nodeCount)
             should_repeat = iteration_cb()
@@ -114,6 +159,8 @@ class GurobiSolver(InternalSolver):
         }
 
     def get_solution(self):
+        self._raise_if_callback()
+
         solution = {}
         for (varname, vardict) in self._all_vars.items():
             solution[varname] = {}
@@ -121,7 +168,22 @@ class GurobiSolver(InternalSolver):
                 solution[varname][idx] = var.x
         return solution
 
+    def get_value(self, var_name, index):
+        var = self._all_vars[var_name][index]
+        return self._get_value(var)
+
+    def _get_value(self, var):
+        if self.cb_where == self.GRB.Callback.MIPSOL:
+            return self.model.cbGetSolution(var)
+        elif self.cb_where == self.GRB.Callback.MIPNODE:
+            return self.model.cbGetNodeRel(var)
+        elif self.cb_where is None:
+            return var.x
+        else:
+            raise Exception("get_value cannot be called from cb_where=%s" % self.cb_where)
+
     def get_variables(self):
+        self._raise_if_callback()
         variables = {}
         for (varname, vardict) in self._all_vars.items():
             variables[varname] = []
@@ -132,12 +194,18 @@ class GurobiSolver(InternalSolver):
     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)
+            if self.cb_where in [self.GRB.Callback.MIPSOL, self.GRB.Callback.MIPNODE]:
+                self.model.cbLazy(lhs, sense, rhs)
+            else:
+                self.model.addConstr(lhs, sense, rhs, name)
         else:
-            self.model.addConstr(constraint, name=name)
+            if self.cb_where in [self.GRB.Callback.MIPSOL, self.GRB.Callback.MIPNODE]:
+                self.model.cbLazy(constraint)
+            else:
+                self.model.addConstr(constraint, name=name)
 
     def set_warm_start(self, solution):
+        self._raise_if_callback()
         count_fixed, count_total = 0, 0
         for (varname, vardict) in solution.items():
             for (idx, value) in vardict.items():
@@ -149,11 +217,13 @@ class GurobiSolver(InternalSolver):
                     (count_fixed, count_total))
 
     def clear_warm_start(self):
+        self._raise_if_callback()
         for (varname, vardict) in self._all_vars:
             for (idx, var) in vardict.items():
                 var[idx].start = self.GRB.UNDEFINED
 
     def fix(self, solution):
+        self._raise_if_callback()
         for (varname, vardict) in solution.items():
             for (idx, value) in vardict.items():
                 if value is None:
@@ -164,10 +234,12 @@ class GurobiSolver(InternalSolver):
                 var.ub = value
 
     def get_constraint_ids(self):
+        self._raise_if_callback()
         self.model.update()
         return [c.ConstrName for c in self.model.getConstrs()]
 
     def extract_constraint(self, cid):
+        self._raise_if_callback()
         constr = self.model.getConstrByName(cid)
         cobj = (self.model.getRow(constr),
                 constr.sense,
@@ -178,29 +250,41 @@ class GurobiSolver(InternalSolver):
 
     def is_constraint_satisfied(self, cobj, tol=1e-5):
         lhs, sense, rhs, name = cobj
-        lhs_value = lhs.getValue()
+        if self.cb_where is not None:
+            lhs_value = lhs.getConstant()
+            for i in range(lhs.size()):
+                var = lhs.getVar(i)
+                coeff = lhs.getCoeff(i)
+                lhs_value += self._get_value(var) * coeff
+        else:
+            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
+            return abs(rhs - lhs_value) < abs(tol)
         else:
             raise Exception("Unknown sense: %s" % sense)
 
     def set_branching_priorities(self, priorities):
+        self._raise_if_callback()
         logger.warning("set_branching_priorities not implemented")
 
     def set_threads(self, threads):
+        self._raise_if_callback()
         self.params["Threads"] = threads
 
     def set_time_limit(self, time_limit):
+        self._raise_if_callback()
         self.params["TimeLimit"] = time_limit
 
     def set_node_limit(self, node_limit):
+        self._raise_if_callback()
         self.params["NodeLimit"] = node_limit
 
     def set_gap_tolerance(self, gap_tolerance):
+        self._raise_if_callback()
         self.params["MIPGap"] = gap_tolerance
 
     def _extract_warm_start_value(self, log):

miplearn/solvers/internal.py

@@ -119,7 +119,7 @@ class InternalSolver(ABC):
         pass
 
     @abstractmethod
-    def solve(self, tee=False, iteration_cb=None):
+    def solve(self, tee=False, iteration_cb=None, lazy_cb=None):
         """
         Solves the currently loaded instance. After this method finishes,
         the best solution found can be retrieved by calling `get_solution`.
@@ -132,7 +132,15 @@ class InternalSolver(ABC):
             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
+        lazy_cb: (internal_solver, model) -> None
+            This function is called whenever the solver finds a new candidate
+            solution and can be used to add lazy constraints to the model. Only
+            two operations within the callback are allowed:
+                - Querying the value of a variable, through `get_value(var, idx)`
+                - Querying if a constraint is satisfied, through `is_constraint_satisfied(cobj)`
+                - Adding a new constraint to the problem, through `add_constraint`
+            Additional operations may be allowed by specific subclasses.
+        tee: Bool
             If true, prints the solver log to the screen.
 
         Returns
@@ -144,6 +152,13 @@ class InternalSolver(ABC):
         """
         pass
 
+    @abstractmethod
+    def get_value(self, var_name, index):
+        """
+        Returns the current value of a decision variable.
+        """
+        pass
+
     @abstractmethod
     def get_constraint_ids(self):
         """

miplearn/solvers/learning.py

@@ -31,34 +31,61 @@ def _parallel_solve(instance_idx):
         "Solution": instance.solution,
         "LP solution": instance.lp_solution,
         "Violated lazy constraints": instance.found_violated_lazy_constraints,
-        "Violated user cuts": instance.found_violated_user_cuts,
+        #"Violated user cuts": instance.found_violated_user_cuts,
     }
 
 
 class LearningSolver:
-    """
-    Mixed-Integer Linear Programming (MIP) solver that extracts information
-    from previous runs, using Machine Learning methods, to accelerate the
-    solution of new (yet unseen) instances.
-
-    Parameters
-    ----------
-    solve_lp_first: bool
-        If true, solve LP relaxation first, then solve original MILP. This
-        option should be activated if the LP relaxation is not very
-        expensive to solve and if it provides good hints for the integer
-        solution.
-    """
-
     def __init__(self,
                  components=None,
-                 gap_tolerance=None,
+                 gap_tolerance=1e-4,
                  mode="exact",
                  solver="gurobi",
                  threads=None,
                  time_limit=None,
                  node_limit=None,
-                 solve_lp_first=True):
+                 solve_lp_first=True,
+                 use_lazy_cb=False):
+        """
+        Mixed-Integer Linear Programming (MIP) solver that extracts information
+        from previous runs and uses Machine Learning methods to accelerate the
+        solution of new (yet unseen) instances.
+
+        Parameters
+        ----------
+        components
+            Set of components in the solver. By default, includes:
+                - ObjectiveValueComponent
+                - PrimalSolutionComponent
+                - DynamicLazyConstraintsComponent
+                - UserCutsComponent
+        gap_tolerance
+            Relative MIP gap tolerance. By default, 1e-4.
+        mode
+            If "exact", solves problem to optimality, keeping all optimality
+            guarantees provided by the MIP solver. If "heuristic", uses machine
+            learning more agressively, and may return suboptimal solutions.
+        solver
+            The internal MIP solver to use. Can be either "cplex", "gurobi", a
+            solver class such as GurobiSolver, or a solver instance such as
+            GurobiSolver().
+        threads
+            Maximum number of threads to use. If None, uses solver default.
+        time_limit
+            Maximum running time in seconds. If None, uses solver default.
+        node_limit
+            Maximum number of branch-and-bound nodes to explore. If None, uses
+            solver default.
+        use_lazy_cb
+            If True, uses lazy callbacks to enforce lazy constraints, instead of
+            a simple solver loop. This functionality may not supported by
+            all internal MIP solvers.
+        solve_lp_first: bool
+            If true, solve LP relaxation first, then solve original MILP. This
+            option should be activated if the LP relaxation is not very
+            expensive to solve and if it provides good hints for the integer
+            solution.
+        """
         self.components = {}
         self.mode = mode
         self.internal_solver = None
@@ -69,6 +96,7 @@ class LearningSolver:
         self.tee = False
         self.node_limit = node_limit
         self.solve_lp_first = solve_lp_first
+        self.use_lazy_cb = use_lazy_cb
 
         if components is not None:
             for comp in components:
@@ -122,14 +150,12 @@ class LearningSolver:
             - instance.lower_bound
             - instance.upper_bound
             - instance.solution
-            - instance.found_violated_lazy_constraints
             - instance.solver_log
-
         Additional solver components may set additional properties. Please
         see their documentation for more details.
 
-        If `solve_lp_first` is False, the properties lp_solution and lp_value
-        will be set to dummy values.
+        If `solver.solve_lp_first` is False, the properties lp_solution and
+        lp_value will be set to dummy values.
 
         Parameters
         ----------
@@ -175,13 +201,23 @@ class LearningSolver:
 
         def iteration_cb():
             should_repeat = False
-            for component in self.components.values():
-                if component.after_iteration(self, instance, model):
+            for comp in self.components.values():
+                if comp.after_iteration(self, instance, model):
                     should_repeat = True
             return should_repeat
 
+        def lazy_cb_wrapper(cb_solver, cb_model):
+            for comp in self.components.values():
+                comp.on_lazy_callback(self, instance, model)
+
+        lazy_cb = None
+        if self.use_lazy_cb:
+            lazy_cb = lazy_cb_wrapper
+
         logger.info("Solving MILP...")
-        results = self.internal_solver.solve(tee=tee, iteration_cb=iteration_cb)
+        results = self.internal_solver.solve(tee=tee,
+                                             iteration_cb=iteration_cb,
+                                             lazy_cb=lazy_cb)
         results["LP value"] = instance.lp_value
 
         # Read MIP solution and bounds
@@ -217,7 +253,7 @@ class LearningSolver:
             instances[idx].lower_bound = r["Results"]["Lower bound"]
             instances[idx].upper_bound = r["Results"]["Upper bound"]
             instances[idx].found_violated_lazy_constraints = r["Violated lazy constraints"]
-            instances[idx].found_violated_user_cuts = r["Violated user cuts"]
+            #instances[idx].found_violated_user_cuts = r["Violated user cuts"]
             instances[idx].solver_log = r["Results"]["Log"]
 
         return results

miplearn/solvers/pyomo/base.py

@@ -57,6 +57,10 @@ class BasePyomoSolver(InternalSolver):
                 solution[str(var)][index] = var[index].value
         return solution
 
+    def get_value(self, var_name, index):
+        var = self._varname_to_var[var_name]
+        return var[index].value
+
     def get_variables(self):
         variables = {}
         for var in self.model.component_objects(Var):
@@ -137,7 +141,12 @@ class BasePyomoSolver(InternalSolver):
         self._pyomo_solver.add_constraint(constraint)
         self._update_constrs()
 
-    def solve(self, tee=False, iteration_cb=None):
+    def solve(self,
+              tee=False,
+              iteration_cb=None,
+              lazy_cb=None):
+        if lazy_cb is not None:
+            raise Exception("lazy callback not supported")
         total_wallclock_time = 0
         streams = [StringIO()]
         if tee:

miplearn/solvers/tests/__init__.py

@@ -7,13 +7,13 @@ from miplearn.problems.knapsack import KnapsackInstance, GurobiKnapsackInstance
 
 
 def _get_instance(solver):
-    if issubclass(solver, BasePyomoSolver):
+    if issubclass(solver, BasePyomoSolver) or isinstance(solver, BasePyomoSolver):
         return KnapsackInstance(
             weights=[23., 26., 20., 18.],
             prices=[505., 352., 458., 220.],
             capacity=67.,
         )
-    if issubclass(solver, GurobiSolver):
+    if issubclass(solver, GurobiSolver) or isinstance(solver, GurobiSolver):
         return GurobiKnapsackInstance(
             weights=[23., 26., 20., 18.],
             prices=[505., 352., 458., 220.],

miplearn/solvers/tests/test_lazy_cb.py

@@ -0,0 +1,27 @@
+#  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 logging
+
+from . import _get_instance
+from ... import GurobiSolver
+
+logger = logging.getLogger(__name__)
+
+
+def test_lazy_cb():
+    solver = GurobiSolver()
+    instance = _get_instance(solver)
+    model = instance.to_model()
+
+    def lazy_cb(cb_solver, cb_model):
+        logger.info("x[0] = %.f" % cb_solver.get_value("x", 0))
+        cobj = (cb_model.getVarByName("x[0]") * 1.0, "<", 0.0, "cut")
+        if not cb_solver.is_constraint_satisfied(cobj):
+            cb_solver.add_constraint(cobj)
+
+    solver.set_instance(instance, model)
+    solver.solve(lazy_cb=lazy_cb)
+    solution = solver.get_solution()
+    assert solution["x"][0] == 0.0