diff --git a/src/python/miplearn/__init__.py b/src/python/miplearn/__init__.py index b519bca..b5f81af 100644 --- a/src/python/miplearn/__init__.py +++ b/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 diff --git a/src/python/miplearn/components/tests/test_branching.py b/src/python/miplearn/components/tests/test_branching.py index c0a264b..ef03e58 100644 --- a/src/python/miplearn/components/tests/test_branching.py +++ b/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() diff --git a/src/python/miplearn/components/tests/test_lazy.py b/src/python/miplearn/components/tests/test_lazy.py index 6e179f9..3c6f2a3 100644 --- a/src/python/miplearn/components/tests/test_lazy.py +++ b/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), diff --git a/src/python/miplearn/components/tests/test_objective.py b/src/python/miplearn/components/tests/test_objective.py index 4809f84..6a9bc5d 100644 --- a/src/python/miplearn/components/tests/test_objective.py +++ b/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) diff --git a/src/python/miplearn/components/tests/test_primal.py b/src/python/miplearn/components/tests/test_primal.py index ee828e8..5147547 100644 --- a/src/python/miplearn/components/tests/test_primal.py +++ b/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 diff --git a/src/python/miplearn/problems/knapsack.py b/src/python/miplearn/problems/knapsack.py index b840fd1..aa08d34 100644 --- a/src/python/miplearn/problems/knapsack.py +++ b/src/python/miplearn/problems/knapsack.py @@ -250,4 +250,26 @@ class KnapsackInstance(Instance): return np.array([ self.weights[index], self.prices[index], - ]) \ No newline at end of file + ]) + + +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 diff --git a/src/python/miplearn/solvers/guroby.py b/src/python/miplearn/solvers/guroby.py new file mode 100644 index 0000000..616a88c --- /dev/null +++ b/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 diff --git a/src/python/miplearn/solvers/tests/__init__.py b/src/python/miplearn/solvers/tests/__init__.py index b5d5159..87a453a 100644 --- a/src/python/miplearn/solvers/tests/__init__.py +++ b/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., - ) \ No newline at end of file +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] diff --git a/src/python/miplearn/solvers/tests/test_internal_solver.py b/src/python/miplearn/solvers/tests/test_internal_solver.py index 77b8def..374d391 100644 --- a/src/python/miplearn/solvers/tests/test_internal_solver.py +++ b/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 \ No newline at end of file +# 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 diff --git a/src/python/miplearn/solvers/tests/test_learning_solver.py b/src/python/miplearn/solvers/tests/test_learning_solver.py index ab6de6e..ea121fb 100644 --- a/src/python/miplearn/solvers/tests/test_learning_solver.py +++ b/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(): diff --git a/src/python/miplearn/tests/__init__.py b/src/python/miplearn/tests/__init__.py index 751ae47..8896939 100644 --- a/src/python/miplearn/tests/__init__.py +++ b/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.],