Implement GurobiSolver (without Pyomo)

src/python/miplearn/__init__.py

@@ -25,5 +25,6 @@ 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.internal import InternalSolver
 from .solvers.learning import LearningSolver

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

@@ -6,11 +6,11 @@ from unittest.mock import Mock
 import numpy as np
 from miplearn import BranchPriorityComponent, BranchPriorityExtractor
 from miplearn.classifiers import Regressor
-from miplearn.tests import get_training_instances_and_models
+from miplearn.tests import get_test_pyomo_instances
 
 
 def test_branch_extract():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     instances[0].branch_priorities = {"x": {0: 100, 1: 200, 2: 300, 3: 400}}
     instances[1].branch_priorities = {"x": {0: 150, 1: 250, 2: 350, 3: 450}}
     priorities = BranchPriorityExtractor().extract(instances)
@@ -18,7 +18,7 @@ def test_branch_extract():
 
 
 def test_branch_calculate():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     comp = BranchPriorityComponent()
 
     # If instances do not have branch_priority property, fit should compute them
@@ -32,7 +32,7 @@ def test_branch_calculate():
 
 
 def test_branch_x_y_predict():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     instances[0].branch_priorities = {"x": {0: 100, 1: 200, 2: 300, 3: 400}}
     instances[1].branch_priorities = {"x": {0: 150, 1: 250, 2: 350, 3: 450}}
     comp = BranchPriorityComponent()

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

@@ -7,14 +7,14 @@ from unittest.mock import Mock
 import numpy as np
 from miplearn import LazyConstraintsComponent, LearningSolver, InternalSolver
 from miplearn.classifiers import Classifier
-from miplearn.tests import get_training_instances_and_models
+from miplearn.tests import get_test_pyomo_instances
 from numpy.linalg import norm
 
 E = 0.1
 
 
 def test_lazy_fit():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     instances[0].found_violated_lazy_constraints = ["a", "b"]
     instances[1].found_violated_lazy_constraints = ["b", "c"]
     classifier = Mock(spec=Classifier)
@@ -51,7 +51,7 @@ def test_lazy_fit():
 
 
 def test_lazy_before():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     instances[0].build_lazy_constraint = Mock(return_value="c1")
     solver = LearningSolver()
     solver.internal_solver = Mock(spec=InternalSolver)
@@ -80,7 +80,7 @@ def test_lazy_before():
 
 
 def test_lazy_evaluate():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     component = LazyConstraintsComponent()
     component.classifiers = {"a": Mock(spec=Classifier),
                              "b": Mock(spec=Classifier),

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

@@ -7,11 +7,11 @@ from unittest.mock import Mock
 import numpy as np
 from miplearn import ObjectiveValueComponent
 from miplearn.classifiers import Regressor
-from miplearn.tests import get_training_instances_and_models
+from miplearn.tests import get_test_pyomo_instances
 
 
 def test_usage():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     comp = ObjectiveValueComponent()
     comp.fit(instances)
     assert instances[0].lower_bound == 1183.0
@@ -21,7 +21,7 @@ def test_usage():
 
 
 def test_obj_evaluate():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     reg = Mock(spec=Regressor)
     reg.predict = Mock(return_value=np.array([1000.0, 1000.0]))
     comp = ObjectiveValueComponent(regressor=reg)

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

@@ -7,11 +7,11 @@ from unittest.mock import Mock
 import numpy as np
 from miplearn import PrimalSolutionComponent
 from miplearn.classifiers import Classifier
-from miplearn.tests import get_training_instances_and_models
+from miplearn.tests import get_test_pyomo_instances
 
 
 def test_predict():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     comp = PrimalSolutionComponent()
     comp.fit(instances)
     solution = comp.predict(instances[0])
@@ -23,7 +23,7 @@ def test_predict():
 
 
 def test_evaluate():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     clf_zero = Mock(spec=Classifier)
     clf_zero.predict_proba = Mock(return_value=np.array([
         [0., 1.],  # x[0]
@@ -93,7 +93,7 @@ def test_evaluate():
 
 
 def test_primal_parallel_fit():
-    instances, models = get_training_instances_and_models()
+    instances, models = get_test_pyomo_instances()
     comp = PrimalSolutionComponent()
     comp.fit(instances, n_jobs=2)
     assert len(comp.classifiers) == 2

src/python/miplearn/problems/knapsack.py

@@ -250,4 +250,26 @@ class KnapsackInstance(Instance):
         return np.array([
             self.weights[index],
             self.prices[index],
-        ])
+        ])
+
+
+class GurobiKnapsackInstance(KnapsackInstance):
+    """
+    Simpler (one-dimensional) knapsack instance, implemented directly in Gurobi
+    instead of Pyomo, used for testing.
+    """
+    def __init__(self, weights, prices, capacity):
+        super().__init__(weights, prices, capacity)
+
+    def to_model(self):
+        import gurobipy as gp
+        from gurobipy import GRB
+
+        model = gp.Model("Knapsack")
+        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)
+        model.setObjective(gp.quicksum(x[i] * self.prices[i]
+                                       for i in range(n)), GRB.MAXIMIZE)
+        return model

src/python/miplearn/solvers/guroby.py

@@ -0,0 +1,202 @@
+#  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 re
+import sys
+import logging
+from io import StringIO
+
+from . import RedirectOutput
+from .internal import InternalSolver
+
+logger = logging.getLogger(__name__)
+
+
+class GurobiSolver(InternalSolver):
+    def __init__(self, params=None):
+        if params is None:
+            params = {
+                "LazyConstraints": 1,
+                "PreCrush": 1,
+            }
+        from gurobipy import GRB
+        self.GRB = GRB
+        self.instance = None
+        self.model = None
+        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:
+            model = instance.to_model()
+        self.instance = instance
+        self.model = model
+        self.model.update()
+        self._update_vars()
+
+    def _update_vars(self):
+        self._all_vars = {}
+        self._bin_vars = {}
+        for var in self.model.getVars():
+            m = re.search(r"([^[]*)\[(.*)\]", var.varName)
+            if m is None:
+                name = var.varName
+                idx = [0]
+            else:
+                name = m.group(1)
+                idx = tuple(int(k) if k.isdecimal() else k
+                            for k in m.group(2).split(","))
+            if len(idx) == 1:
+                idx = idx[0]
+            if name not in self._all_vars:
+                self._all_vars[name] = {}
+            self._all_vars[name][idx] = var
+            if var.vtype != 'C':
+                if name not in self._bin_vars:
+                    self._bin_vars[name] = {}
+                self._bin_vars[name][idx] = var
+
+    def _apply_params(self):
+        for (name, value) in self.params.items():
+            self.model.setParam(name, value)
+
+    def solve_lp(self, tee=False):
+        self._apply_params()
+        streams = [StringIO()]
+        if tee:
+            streams += [sys.stdout]
+        for (varname, vardict) in self._bin_vars.items():
+            for (idx, var) in vardict.items():
+                var.vtype = self.GRB.CONTINUOUS
+                var.lb = 0.0
+                var.ub = 1.0
+        with RedirectOutput(streams):
+            self.model.optimize()
+        for (varname, vardict) in self._bin_vars.items():
+            for (idx, var) in vardict.items():
+                var.vtype = self.GRB.BINARY
+        log = streams[0].getvalue()
+        return {
+            "Optimal value": self.model.objVal,
+            "Log": log
+        }
+
+    def solve(self, tee=False):
+        all_vars = self.model.getVars()
+        self.instance.found_violated_lazy_constraints = []
+        self.instance.found_violated_user_cuts = []
+        streams = [StringIO()]
+        if tee:
+            streams += [sys.stdout]
+
+        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),
+            "Sense": ("min" if self.model.modelSense == 1 else "max"),
+            "Log": log,
+            "Warm start value": self._extract_warm_start_value(log),
+        }
+
+    def get_solution(self):
+        solution = {}
+        for (varname, vardict) in self._all_vars.items():
+            solution[varname] = {}
+            for (idx, var) in vardict.items():
+                solution[varname][idx] = var.x
+        return solution
+
+    def add_constraint(self, constraint):
+        self.model.addConstr(constraint)
+
+    def set_warm_start(self, solution):
+        count_fixed, count_total = 0, 0
+        for (varname, vardict) in solution.items():
+            for (idx, value) in vardict.items():
+                count_total += 1
+                if value is not None:
+                    count_fixed += 1
+                    self._all_vars[varname][idx].start = value
+        logger.info("Setting start values for %d variables (out of %d)" %
+                    (count_fixed, count_total))
+
+    def clear_warm_start(self):
+        for (varname, vardict) in self._all_vars:
+            for (idx, var) in vardict.items():
+                var[idx].start = self.GRB.UNDEFINED
+
+    def fix(self, solution):
+        for (varname, vardict) in solution.items():
+            for (idx, value) in vardict.items():
+                if value is None:
+                    continue
+                var = self._all_vars[varname][idx]
+                var.vtype = self.GRB.CONTINUOUS
+                var.lb = value
+                var.ub = value
+
+    def set_branching_priorities(self, priorities):
+        logger.warning("set_branching_priorities not implemented")
+
+    def set_threads(self, threads):
+        self.params["Threads"] = threads
+
+    def set_time_limit(self, time_limit):
+        self.params["TimeLimit"] = time_limit
+
+    def set_node_limit(self, node_limit):
+        self.params["NodeLimit"] = node_limit
+
+    def set_gap_tolerance(self, gap_tolerance):
+        self.params["MIPGap"] = gap_tolerance
+
+    def _extract_warm_start_value(self, log):
+        ws = self.__extract(log, "MIP start with objective ([0-9.e+-]*)")
+        if ws is not None:
+            ws = float(ws)
+        return ws
+
+    def __extract(self, log, regexp, default=None):
+        value = default
+        for line in log.splitlines():
+            matches = re.findall(regexp, line)
+            if len(matches) == 0:
+                continue
+            value = matches[0]
+        return value
+
+    def __getstate__(self):
+        return self.params
+
+    def __setstate__(self, state):
+        self.params = state

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

@@ -2,12 +2,25 @@
 #  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 #  Released under the modified BSD license. See COPYING.md for more details.
 
-from miplearn.problems.knapsack import KnapsackInstance
+from miplearn import BasePyomoSolver, GurobiSolver, GurobiPyomoSolver, CplexPyomoSolver
+from miplearn.problems.knapsack import KnapsackInstance, GurobiKnapsackInstance
 
 
-def _get_instance():
-    return KnapsackInstance(
-        weights=[23., 26., 20., 18.],
-        prices=[505., 352., 458., 220.],
-        capacity=67.,
-    )
+def _get_instance(solver):
+    if issubclass(solver, BasePyomoSolver):
+        return KnapsackInstance(
+            weights=[23., 26., 20., 18.],
+            prices=[505., 352., 458., 220.],
+            capacity=67.,
+        )
+    if issubclass(solver, GurobiSolver):
+        return GurobiKnapsackInstance(
+            weights=[23., 26., 20., 18.],
+            prices=[505., 352., 458., 220.],
+            capacity=67.,
+        )
+    assert False
+
+
+def _get_internal_solvers():
+    return [GurobiPyomoSolver, CplexPyomoSolver, GurobiSolver]

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

@@ -1,14 +1,18 @@
 #  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 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 CplexPyomoSolver, GurobiPyomoSolver
 
-from . import _get_instance
-from ...problems.knapsack import ChallengeA
+from . import _get_instance, _get_internal_solvers
+
+logger = logging.getLogger(__name__)
 
 
 def test_redirect_output():
@@ -22,10 +26,11 @@ def test_redirect_output():
 
 
 def test_internal_solver_warm_starts():
-    for solver in [GurobiPyomoSolver(), CplexPyomoSolver()]:
-        instance = _get_instance()
+    for solver_class in _get_internal_solvers():
+        logger.info("Solver: %s" % solver_class)
+        instance = _get_instance(solver_class)
         model = instance.to_model()
-
+        solver = solver_class()
         solver.set_instance(instance, model)
         solver.set_warm_start({
             "x": {
@@ -63,11 +68,23 @@ def test_internal_solver_warm_starts():
 
 
 def test_internal_solver():
-    for solver in [GurobiPyomoSolver(), CplexPyomoSolver()]:
-        instance = _get_instance()
-        model = instance.to_model()
+    for solver_class in _get_internal_solvers():
+        logger.info("Solver: %s" % solver_class)
 
+        instance = _get_instance(solver_class)
+        model = instance.to_model()
+        solver = solver_class()
         solver.set_instance(instance, model)
+
+        stats = solver.solve_lp()
+        assert round(stats["Optimal value"], 3) == 1287.923
+
+        solution = solver.get_solution()
+        assert round(solution["x"][0], 3) == 1.000
+        assert round(solution["x"][1], 3) == 0.923
+        assert round(solution["x"][2], 3) == 1.000
+        assert round(solution["x"][3], 3) == 0.000
+
         stats = solver.solve(tee=True)
         assert len(stats["Log"]) > 100
         assert stats["Lower bound"] == 1183.0
@@ -82,27 +99,17 @@ def test_internal_solver():
         assert solution["x"][2] == 1.0
         assert solution["x"][3] == 1.0
 
-        stats = solver.solve_lp()
-        assert round(stats["Optimal value"], 3) == 1287.923
-
-        solution = solver.get_solution()
-        assert round(solution["x"][0], 3) == 1.000
-        assert round(solution["x"][1], 3) == 0.923
-        assert round(solution["x"][2], 3) == 1.000
-        assert round(solution["x"][3], 3) == 0.000
-
-        model.cut = pe.Constraint(expr=model.x[0] <= 0.5)
-        solver.add_constraint(model.cut)
-        solver.solve_lp()
-        assert model.x[0].value == 0.5
+        if isinstance(solver, BasePyomoSolver):
+            model.cut = pe.Constraint(expr=model.x[0] <= 0.5)
+            solver.add_constraint(model.cut)
+            solver.solve_lp()
+            assert model.x[0].value == 0.5
 
 
-def test_node_count():
-    for solver in [GurobiPyomoSolver(),
-                   GurobiPyomoSolver(use_lazy_callbacks=False),
-                   CplexPyomoSolver()]:
-        challenge = ChallengeA()
-        solver.set_time_limit(1)
-        solver.set_instance(challenge.test_instances[0])
-        stats = solver.solve(tee=True)
-        assert stats["Nodes"] > 1
+# 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

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

@@ -2,19 +2,23 @@
 #  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 #  Released under the modified BSD license. See COPYING.md for more details.
 
+import logging
 import pickle
 import tempfile
 
-from miplearn import BranchPriorityComponent, GurobiPyomoSolver
+from miplearn import BranchPriorityComponent
 from miplearn import LearningSolver
 
-from . import _get_instance
+from . import _get_instance, _get_internal_solvers
+
+logger = logging.getLogger(__name__)
 
 
 def test_learning_solver():
-    instance = _get_instance()
     for mode in ["exact", "heuristic"]:
-        for internal_solver in ["cplex", "gurobi", GurobiPyomoSolver]:
+        for internal_solver in _get_internal_solvers():
+            logger.info("Solver: %s" % internal_solver)
+            instance = _get_instance(internal_solver)
             solver = LearningSolver(time_limit=300,
                                     gap_tolerance=1e-3,
                                     threads=1,
@@ -46,12 +50,13 @@ def test_learning_solver():
 
 
 def test_parallel_solve():
-    instances = [_get_instance() for _ in range(10)]
-    solver = LearningSolver()
-    results = solver.parallel_solve(instances, n_jobs=3)
-    assert len(results) == 10
-    for instance in instances:
-        assert len(instance.solution["x"].keys()) == 4
+    for internal_solver in _get_internal_solvers():
+        instances = [_get_instance(internal_solver) for _ in range(10)]
+        solver = LearningSolver(solver=internal_solver)
+        results = solver.parallel_solve(instances, n_jobs=3)
+        assert len(results) == 10
+        for instance in instances:
+            assert len(instance.solution["x"].keys()) == 4
 
 
 def test_add_components():

src/python/miplearn/tests/__init__.py

@@ -5,7 +5,7 @@ from miplearn import LearningSolver
 from miplearn.problems.knapsack import KnapsackInstance
 
 
-def get_training_instances_and_models():
+def get_test_pyomo_instances():
     instances = [
         KnapsackInstance(
             weights=[23., 26., 20., 18.],