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

@@ -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)

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())

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()

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