From f1c8bb37c5a2cab4a0b9d8e72de58c2509390ee5 Mon Sep 17 00:00:00 2001 From: Nikolaus Sonnenschein Date: Tue, 7 Nov 2017 15:22:44 +0100 Subject: [PATCH] Add an interface to GLPK's exact solver (#145) * feat: add GLPK exact solver interface * fix: don't monkey patch the glpk_exact_interface test suite * style: flake8 --- optlang/glpk_exact_interface.py | 150 ++++++++++++++++++ optlang/glpk_interface.py | 2 +- optlang/tests/test_glpk_exact_interface.py | 38 +++++ optlang/tests/test_glpk_interface.py | 63 ++++---- .../tests/test_netlib_glpk_exact_interface.py | 130 +++++++++++++++ 5 files changed, 350 insertions(+), 33 deletions(-) create mode 100644 optlang/glpk_exact_interface.py create mode 100644 optlang/tests/test_glpk_exact_interface.py create mode 100644 optlang/tests/test_netlib_glpk_exact_interface.py diff --git a/optlang/glpk_exact_interface.py b/optlang/glpk_exact_interface.py new file mode 100644 index 00000000..1e3bd991 --- /dev/null +++ b/optlang/glpk_exact_interface.py @@ -0,0 +1,150 @@ +# Copyright 2017 Novo Nordisk Foundation Center for Biosustainability, +# Technical University of Denmark. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +""" +Interface for the GNU Linear Programming Kit (GLPK) + +GLPK is an open source LP solver, with MILP capabilities. This interface exposes its GLPK's exact solver. +To use GLPK you need to install the 'swiglpk' python package (with pip or from http://github.com/biosustain/swiglpk) +and make sure that 'import swiglpk' runs without error. +""" + +import logging + +import six + +from optlang.util import inheritdocstring +from optlang import interface +from optlang import glpk_interface +from optlang.glpk_interface import _GLPK_STATUS_TO_STATUS + +log = logging.getLogger(__name__) + +from swiglpk import glp_exact, glp_create_prob, glp_get_status, \ + GLP_SF_AUTO, GLP_ETMLIM, glp_adv_basis, glp_read_lp, glp_scale_prob + + +@six.add_metaclass(inheritdocstring) +class Variable(glpk_interface.Variable): + pass + + +@six.add_metaclass(inheritdocstring) +class Constraint(glpk_interface.Constraint): + pass + + +@six.add_metaclass(inheritdocstring) +class Objective(glpk_interface.Objective): + pass + + +@six.add_metaclass(inheritdocstring) +class Configuration(glpk_interface.Configuration): + pass + + +@six.add_metaclass(inheritdocstring) +class Model(glpk_interface.Model): + def _run_glp_exact(self): + return_value = glp_exact(self.problem, self.configuration._smcp) + glpk_status = glp_get_status(self.problem) + if return_value == 0: + status = _GLPK_STATUS_TO_STATUS[glpk_status] + elif return_value == GLP_ETMLIM: + status = interface.TIME_LIMIT + else: + status = _GLPK_STATUS_TO_STATUS[glpk_status] + if status == interface.UNDEFINED: + log.debug("Status undefined. GLPK status code returned by glp_simplex was %d" % return_value) + return status + + def _optimize(self): + # Solving inexact first per GLPK manual + # Computations in exact arithmetic are very time consuming, so solving LP + # problem with the routine glp_exact from the very beginning is not a good + # idea. It is much better at first to find an optimal basis with the routine + # glp_simplex and only then to call glp_exact, in which case only a few + # simplex iterations need to be performed in exact arithmetic. + status = super(Model, self)._optimize() + if status != interface.OPTIMAL: + return status + else: + status = self._run_glp_exact() + + # Sometimes GLPK gets itself stuck with an invalid basis. This will help it get rid of it. + if status == interface.UNDEFINED and self.configuration.presolve is not True: + glp_adv_basis(self.problem, 0) + status = self._run_glp_exact() + + if status == interface.UNDEFINED and self.configuration.presolve is True: + # If presolve is on, status will be undefined if not optimal + self.configuration.presolve = False + status = self._run_glp_exact() + self.configuration.presolve = True + if self._glpk_is_mip(): + status = self._run_glp_mip() + if status == 'undefined' or status == 'infeasible': + # Let's see if the presolver and some scaling can fix this issue + glp_scale_prob(self.problem, GLP_SF_AUTO) + original_presolve_setting = self.configuration.presolve + self.configuration.presolve = True + status = self._run_glp_mip() + self.configuration.presolve = original_presolve_setting + return status + + +if __name__ == '__main__': + import pickle + + x1 = Variable('x1', lb=0) + x2 = Variable('x2', lb=0) + x3 = Variable('x3', lb=0, ub=1, type='binary') + c1 = Constraint(x1 + x2 + x3, lb=-100, ub=100, name='c1') + c2 = Constraint(10 * x1 + 4 * x2 + 5 * x3, ub=600, name='c2') + c3 = Constraint(2 * x1 + 2 * x2 + 6 * x3, ub=300, name='c3') + obj = Objective(10 * x1 + 6 * x2 + 4 * x3, direction='max') + model = Model(name='Simple model') + model.objective = obj + model.add([c1, c2, c3]) + model.configuration.verbosity = 3 + status = model.optimize() + print("status:", model.status) + print("objective value:", model.objective.value) + + for var_name, var in model.variables.items(): + print(var_name, "=", var.primal) + + print(model) + + problem = glp_create_prob() + glp_read_lp(problem, None, "tests/data/model.lp") + + solver = Model(problem=problem) + print(solver.optimize()) + print(solver.objective) + + import time + + t1 = time.time() + print("pickling") + pickle_string = pickle.dumps(solver) + resurrected_solver = pickle.loads(pickle_string) + t2 = time.time() + print("Execution time: %s" % (t2 - t1)) + + resurrected_solver.optimize() + print("Halelujah!", resurrected_solver.objective.value) diff --git a/optlang/glpk_interface.py b/optlang/glpk_interface.py index 6eb2562b..633433f7 100644 --- a/optlang/glpk_interface.py +++ b/optlang/glpk_interface.py @@ -839,7 +839,7 @@ def _remove_constraints(self, constraints): print(model) problem = glp_create_prob() - glp_read_lp(problem, None, "../tests/data/model.lp") + glp_read_lp(problem, None, "tests/data/model.lp") solver = Model(problem=problem) print(solver.optimize()) diff --git a/optlang/tests/test_glpk_exact_interface.py b/optlang/tests/test_glpk_exact_interface.py new file mode 100644 index 00000000..24218dd1 --- /dev/null +++ b/optlang/tests/test_glpk_exact_interface.py @@ -0,0 +1,38 @@ +# Copyright (c) 2013 Novo Nordisk Foundation Center for Biosustainability, DTU. +# See LICENSE for details. + +import os +import random + +import nose + +from optlang import glpk_exact_interface +from optlang.tests import test_glpk_interface + +random.seed(666) +TESTMODELPATH = os.path.join(os.path.dirname(__file__), 'data/model.lp') +TESTMILPMODELPATH = os.path.join(os.path.dirname(__file__), 'data/simple_milp.lp') + + +class VariableTestCase(test_glpk_interface.VariableTestCase): + interface = glpk_exact_interface + + +class ConstraintTestCase(test_glpk_interface.ConstraintTestCase): + interface = glpk_exact_interface + + +class ObjectiveTestCase(test_glpk_interface.ObjectiveTestCase): + interface = glpk_exact_interface + + +class ConfigurationTestCase(test_glpk_interface.ConfigurationTestCase): + interface = glpk_exact_interface + + +class ModelTestCase(test_glpk_interface.ModelTestCase): + interface = glpk_exact_interface + + +if __name__ == '__main__': + nose.runmodule() diff --git a/optlang/tests/test_glpk_interface.py b/optlang/tests/test_glpk_interface.py index 0b00e18d..85aedf09 100644 --- a/optlang/tests/test_glpk_interface.py +++ b/optlang/tests/test_glpk_interface.py @@ -9,13 +9,12 @@ import nose import optlang from optlang import glpk_interface -from optlang.glpk_interface import Variable, Constraint, Model, Objective from optlang.tests import abstract_test_cases from optlang.util import glpk_read_cplex from swiglpk import glp_get_num_rows, glp_get_col_name, glp_get_num_cols, glp_get_prob_name, glp_get_row_name, \ glp_get_col_kind, glp_find_col, intArray, doubleArray, glp_get_mat_row, glp_get_row_type, glp_get_row_lb, \ glp_get_row_ub, glp_get_obj_coef, GLP_UP, GLP_DB, GLP_LO, GLP_CV, GLP_IV, GLP_FX, GLP_FR, glp_get_col_lb, \ - glp_get_col_ub + glp_get_col_ub, glp_get_obj_dir, GLP_MIN, GLP_MAX random.seed(666) TESTMODELPATH = os.path.join(os.path.dirname(__file__), 'data/model.lp') @@ -30,7 +29,7 @@ def test_variable_without_problem_returns_None_index(self): def test_get_primal(self): self.assertEqual(self.var.primal, None) - model = Model(problem=glpk_read_cplex(TESTMODELPATH)) + model = glpk_interface.Model(problem=glpk_read_cplex(TESTMODELPATH)) model.optimize() for i, j in zip([var.primal for var in model.variables], [0.8739215069684306, -16.023526143167608, 16.023526143167604, -14.71613956874283, @@ -52,7 +51,7 @@ def test_get_primal(self): self.assertAlmostEqual(i, j) def test_changing_variable_names_is_reflected_in_the_solver(self): - model = Model(problem=glpk_read_cplex(TESTMODELPATH)) + model = self.interface.Model(problem=glpk_read_cplex(TESTMODELPATH)) for i, variable in enumerate(model.variables): variable.name = "var" + str(i) self.assertEqual(variable.name, "var" + str(i)) @@ -65,10 +64,10 @@ def test_glpk_setting_bounds(self): model = self.model var.lb = 1 self.assertEqual(var.lb, 1) - self.assertEqual(glpk_interface.glp_get_col_lb(model.problem, var._index), 1) + self.assertEqual(glp_get_col_lb(model.problem, var._index), 1) var.ub = 2 self.assertEqual(var.ub, 2) - self.assertEqual(glpk_interface.glp_get_col_ub(model.problem, var._index), 2) + self.assertEqual(glp_get_col_ub(model.problem, var._index), 2) class ConstraintTestCase(abstract_test_cases.AbstractConstraintTestCase): @@ -98,17 +97,17 @@ class ObjectiveTestCase(abstract_test_cases.AbstractObjectiveTestCase): interface = glpk_interface def setUp(self): - self.model = Model(problem=glpk_read_cplex(TESTMODELPATH)) + self.model = self.interface.Model(problem=glpk_read_cplex(TESTMODELPATH)) self.obj = self.model.objective def test_change_direction(self): self.obj.direction = "min" self.assertEqual(self.obj.direction, "min") - self.assertEqual(glpk_interface.glp_get_obj_dir(self.model.problem), glpk_interface.GLP_MIN) + self.assertEqual(glp_get_obj_dir(self.model.problem), GLP_MIN) self.obj.direction = "max" self.assertEqual(self.obj.direction, "max") - self.assertEqual(glpk_interface.glp_get_obj_dir(self.model.problem), glpk_interface.GLP_MAX) + self.assertEqual(glp_get_obj_dir(self.model.problem), GLP_MAX) class ConfigurationTestCase(abstract_test_cases.AbstractConfigurationTestCase): @@ -119,7 +118,7 @@ class ModelTestCase(abstract_test_cases.AbstractModelTestCase): interface = glpk_interface def test_glpk_create_empty_model(self): - model = Model(name="empty_problem") + model = self.interface.Model(name="empty_problem") self.assertEqual(glp_get_prob_name(model.problem), "empty_problem") def test_pickle_ability(self): @@ -150,7 +149,7 @@ def test_init_from_existing_problem(self): [glp_get_row_name(inner_prob, j) for j in range(1, glp_get_num_rows(inner_prob) + 1)]) def test_add_non_cplex_conform_variable(self): - var = Variable('12x!!@#5_3', lb=-666, ub=666) + var = self.interface.Variable('12x!!@#5_3', lb=-666, ub=666) self.assertEqual(var._index, None) self.model.add(var) self.assertTrue(var in self.model.variables.values()) @@ -194,14 +193,14 @@ def test_glpk_remove_variable(self): self.assertEqual(var.problem, None) def test_add_constraints(self): - x = Variable('x', lb=0, ub=1, type='binary') - y = Variable('y', lb=-181133.3, ub=12000., type='continuous') - z = Variable('z', lb=0., ub=10., type='integer') - constr1 = Constraint(0.3 * x + 0.4 * y + 66. * z, lb=-100, ub=0., name='test') - constr2 = Constraint(2.333 * x + y + 3.333, ub=100.33, name='test2') - constr3 = Constraint(2.333 * x + y + z, lb=-300) - constr4 = Constraint(x, lb=-300, ub=-300) - constr5 = Constraint(3 * x) + x = self.interface.Variable('x', lb=0, ub=1, type='binary') + y = self.interface.Variable('y', lb=-181133.3, ub=12000., type='continuous') + z = self.interface.Variable('z', lb=0., ub=10., type='integer') + constr1 = self.interface.Constraint(0.3 * x + 0.4 * y + 66. * z, lb=-100, ub=0., name='test') + constr2 = self.interface.Constraint(2.333 * x + y + 3.333, ub=100.33, name='test2') + constr3 = self.interface.Constraint(2.333 * x + y + z, lb=-300) + constr4 = self.interface.Constraint(x, lb=-300, ub=-300) + constr5 = self.interface.Constraint(3 * x) self.model.add(constr1) self.model.add(constr2) self.model.add(constr3) @@ -269,9 +268,9 @@ def test_add_constraints(self): self.assertGreater(glp_get_row_ub(self.model.problem, constr5._index), 1e30) def test_change_of_constraint_is_reflected_in_low_level_solver(self): - x = Variable('x', lb=-83.3, ub=1324422.) - y = Variable('y', lb=-181133.3, ub=12000.) - constraint = Constraint(0.3 * x + 0.4 * y, lb=-100, name='test') + x = self.interface.Variable('x', lb=-83.3, ub=1324422.) + y = self.interface.Variable('y', lb=-181133.3, ub=12000.) + constraint = self.interface.Constraint(0.3 * x + 0.4 * y, lb=-100, name='test') self.assertEqual(constraint._index, None) self.model.add(constraint) self.assertEqual( @@ -280,7 +279,7 @@ def test_change_of_constraint_is_reflected_in_low_level_solver(self): self.assertEqual(self.model.constraints["test"].lb, -100) self.assertEqual(constraint._index, 73) - z = Variable('z', lb=3, ub=10, type='integer') + z = self.interface.Variable('z', lb=3, ub=10, type='integer') self.assertEqual(z._index, None) constraint += 77. * z self.assertEqual(z._index, 98) @@ -293,11 +292,11 @@ def test_change_of_constraint_is_reflected_in_low_level_solver(self): self.assertEqual(constraint._index, 73) def test_constraint_set_problem_to_None_caches_the_latest_expression_from_solver_instance(self): - x = Variable('x', lb=-83.3, ub=1324422.) - y = Variable('y', lb=-181133.3, ub=12000.) - constraint = Constraint(0.3 * x + 0.4 * y, lb=-100, name='test') + x = self.interface.Variable('x', lb=-83.3, ub=1324422.) + y = self.interface.Variable('y', lb=-181133.3, ub=12000.) + constraint = self.interface.Constraint(0.3 * x + 0.4 * y, lb=-100, name='test') self.model.add(constraint) - z = Variable('z', lb=2, ub=5, type='integer') + z = self.interface.Variable('z', lb=2, ub=5, type='integer') constraint += 77. * z self.model.remove(constraint) self.assertEqual( @@ -307,9 +306,9 @@ def test_constraint_set_problem_to_None_caches_the_latest_expression_from_solver self.assertEqual(constraint.ub, None) def test_change_of_objective_is_reflected_in_low_level_solver(self): - x = Variable('x', lb=-83.3, ub=1324422.) - y = Variable('y', lb=-181133.3, ub=12000.) - objective = Objective(0.3 * x + 0.4 * y, name='test', direction='max') + x = self.interface.Variable('x', lb=-83.3, ub=1324422.) + y = self.interface.Variable('y', lb=-181133.3, ub=12000.) + objective = self.interface.Objective(0.3 * x + 0.4 * y, name='test', direction='max') self.model.objective = objective self.assertEqual( @@ -322,7 +321,7 @@ def test_change_of_objective_is_reflected_in_low_level_solver(self): for i in range(1, glp_get_num_cols(self.model.problem) + 1): if i != x._index and i != y._index: self.assertEqual(glp_get_obj_coef(self.model.problem, i), 0) - z = Variable('z', lb=4, ub=4, type='integer') + z = self.interface.Variable('z', lb=4, ub=4, type='integer') self.model.objective += 77. * z self.assertEqual( @@ -416,7 +415,7 @@ def test_set_linear_coefficients_objective(self): self.assertEqual(glp_get_obj_coef(self.model.problem, self.model.variables.R_TPI._index), 666.) def test_instantiating_model_with_different_solver_problem_raises(self): - self.assertRaises(TypeError, Model, problem='Chicken soup') + self.assertRaises(TypeError, self.interface.Model, problem='Chicken soup') def test_set_linear_coefficients_constraint(self): constraint = self.model.constraints.M_atp_c diff --git a/optlang/tests/test_netlib_glpk_exact_interface.py b/optlang/tests/test_netlib_glpk_exact_interface.py new file mode 100644 index 00000000..56257f7c --- /dev/null +++ b/optlang/tests/test_netlib_glpk_exact_interface.py @@ -0,0 +1,130 @@ +# Copyright (c) 2013 Novo Nordisk Foundation Center for Biosustainability, DTU. +# See LICENSE for details. + +import glob +import os +import pickle +import tarfile +import tempfile +from functools import partial + +import nose +import six +from optlang.glpk_exact_interface import Model +from swiglpk import glp_create_prob, GLP_MPS_DECK, glp_read_mps, glp_get_num_cols, glp_smcp, glp_exact, \ + glp_get_status, \ + GLP_OPT, glp_get_obj_val + + +def test_netlib(netlib_tar_path=os.path.join(os.path.dirname(__file__), 'data/netlib_lp_problems.tar.gz')): + """ + Test netlib with glpk interface + """ + if six.PY3: + nose.SkipTest('Skipping because py3') + elif os.getenv('TRAVIS'): + nose.SkipTest('Skipping because travis (solving problems exactly takes very long.)') + else: + with open(os.path.join(os.path.dirname(__file__), 'data/the_final_netlib_results.pcl'), 'rb') as fhandle: + THE_FINAL_NETLIB_RESULTS = pickle.load(fhandle) + + # noinspection PyShadowingNames + def read_netlib_sif_glpk(fhandle): + tmp_file = tempfile.mktemp(suffix='.mps') + with open(tmp_file, 'w') as tmp_handle: + content = ''.join([str(s) for s in fhandle if str(s.strip())]) + tmp_handle.write(content) + fhandle.close() + problem = glp_create_prob() + glp_read_mps(problem, GLP_MPS_DECK, None, tmp_file) + # glp_read_mps(problem, GLP_MPS_FILE, None, tmp_file) + return problem + + def check_dimensions(glpk_problem, model): + """ + Tests that the glpk problem and the interface model have the same + number of rows (constraints) and columns (variables). + """ + assert glp_get_num_cols(glpk_problem) == len(model.variables) + + def check_objval_against_the_final_netlib_results(netlib_id, model_objval): + print(float(THE_FINAL_NETLIB_RESULTS[netlib_id]['Objvalue']), model_objval) + relative_error = abs(1 - (model_objval / float(THE_FINAL_NETLIB_RESULTS[netlib_id]['Objvalue']))) + print(relative_error) + nose.tools.assert_true(relative_error < 0.01) + # nose.tools.assert_almost_equal(model_objval, float(THE_FINAL_NETLIB_RESULTS[netlib_id]['Objvalue']), places=4) + + tar = tarfile.open(netlib_tar_path) + model_paths_in_tar = glob.fnmatch.filter(tar.getnames(), '*.SIF') + + for model_path_in_tar in model_paths_in_tar: + netlib_id = os.path.basename(model_path_in_tar).replace('.SIF', '') + # TODO: get the following problems to work + # E226 seems to be a MPS related problem, see http://lists.gnu.org/archive/html/bug-glpk/2003-01/msg00003.html + if netlib_id in ('AGG', 'E226', 'SCSD6', 'DFL001'): + # def test_skip(netlib_id): + # raise SkipTest('Skipping netlib problem %s ...' % netlib_id) + # test_skip(netlib_id) + # class TestWeirdNetlibProblems(unittest.TestCase): + + # @unittest.skip('Skipping netlib problem') + # def test_fail(): + # pass + continue + # TODO: For now, test only models that are covered by the final netlib results + else: + if netlib_id not in THE_FINAL_NETLIB_RESULTS.keys(): + continue + fhandle = tar.extractfile(model_path_in_tar) + glpk_problem = read_netlib_sif_glpk(fhandle) + model = Model(problem=glpk_problem) + model.configuration.presolve = True + model.configuration.verbosity = 3 + func = partial(check_dimensions, glpk_problem, model) + func.description = "test_netlib_check_dimensions_%s (%s)" % (netlib_id, os.path.basename(str(__file__))) + yield func + + model.optimize() + if model.status == 'optimal': + model_objval = model.objective.value + else: + raise Exception('No optimal solution found for netlib model %s' % netlib_id) + + func = partial(check_objval_against_the_final_netlib_results, netlib_id, model_objval) + func.description = "test_netlib_check_objective_value__against_the_final_netlib_results_%s (%s)" % ( + netlib_id, os.path.basename(str(__file__))) + yield func + + if not os.getenv('TRAVIS', False): + # check that a cloned model also gives the correct result + model = Model.clone(model, use_json=False, use_lp=False) + model.optimize() + if model.status == 'optimal': + model_objval = model.objective.value + else: + raise Exception('No optimal solution found for netlib model %s' % netlib_id) + + if not os.getenv('TRAVIS', False): + func = partial(check_objval_against_the_final_netlib_results, netlib_id, model_objval) + func.description = "test_netlib_check_objective_value__against_the_final_netlib_results_after_cloning_%s (%s)" % ( + netlib_id, os.path.basename(str(__file__))) + yield func + + +if __name__ == '__main__': + # tar = tarfile.open('data/netlib_lp_problems.tar.gz') + # model_paths_in_tar = glob.fnmatch.filter(tar.getnames(), '*.SIF') + # fhandle = tar.extractfile('./netlib/ADLITTLE.SIF') + # glpk_problem = read_netlib_sif_glpk(fhandle) + # glp_simplex(glpk_problem, None) + # print glp_get_obj_val(glpk_problem) + # print glpk_problem + # fhandle = tar.extractfile('./netlib/ADLITTLE.SIF') + # glpk_problem = read_netlib_sif_glpk(fhandle) + # model = Model(problem=glpk_problem) + # glp_simplex(glpk_problem, None) + # model.optimize() + # print model.objective.value + # print model + # test_netlib().next() + nose.runmodule()