Add QP capability for Gurobi (#149)

* Add gurobi QPs

* add some fixes and QP tests for gurobi

optlang/gurobi_interface.py

@@ -54,6 +54,8 @@
 
 _LP_METHODS = {"auto": -1, "primal": 0, "dual": 1, "barrier": 2, "concurrent": 3, "deterministic_concurrent": 4}
 _REVERSE_LP_METHODS = {v: k for k, v in _LP_METHODS.items()}
+_QP_METHODS = {"auto": -1, "primal": 0, "dual": 1, "barrier": 2}
+_REVERSE_QP_METHODS = {v: k for k, v in _QP_METHODS.items()}
 
 _VTYPE_TO_GUROBI_VTYPE = {'continuous': gurobipy.GRB.CONTINUOUS, 'integer': gurobipy.GRB.INTEGER,
                           'binary': gurobipy.GRB.BINARY}
@@ -329,8 +331,8 @@ class Objective(interface.Objective):
     def __init__(self, expression, sloppy=False, *args, **kwargs):
         super(Objective, self).__init__(expression, *args, sloppy=sloppy, **kwargs)
         self._expression_expired = False
-        if not (sloppy or self.is_Linear):  # or self.is_Quadratic: # QP is not yet supported
-            raise ValueError("The given objective is invalid. Must be linear or quadratic (not yet supported")
+        if not (sloppy or self.is_Linear or self.is_Quadratic):
+            raise ValueError("The given objective is invalid. Must be linear or quadratic.")
 
     @property
     def value(self):
@@ -371,22 +373,35 @@ def get_linear_coefficients(self, variables):
     def _get_expression(self):
         if self.problem is not None and self._expression_expired and len(self.problem._variables) > 0:
             grb_obj = self.problem.problem.getObjective()
-            terms = []
             variables = self.problem._variables
-            for i in range(grb_obj.size()):
-                terms.append(grb_obj.getCoeff(i) * variables[grb_obj.getVar(i).getAttr('VarName')])
-            expression = symbolics.add(terms)
-            # TODO implement quadratic objectives
-            self._expression = expression + getattr(self.problem, "_objective_offset", 0)
+            if self.problem.problem.IsQP:
+                quadratic_expression = symbolics.add(
+                    [symbolics.Real(grb_obj.getCoeff(i)) *
+                     variables[grb_obj.getVar1(i).VarName] *
+                     variables[grb_obj.getVar2(i).VarName]
+                     for i in range(grb_obj.size())])
+                linear_objective = grb_obj.getLinExpr()
+            else:
+                quadratic_expression = symbolics.Real(0.0)
+                linear_objective = grb_obj
+            linear_expression = symbolics.add(
+                [symbolics.Real(linear_objective.getCoeff(i)) *
+                 variables[linear_objective.getVar(i).VarName]
+                 for i in range(linear_objective.size())]
+            )
+            self._expression = (linear_expression + quadratic_expression +
+                                getattr(self.problem, "_objective_offset", 0))
             self._expression_expired = False
         return self._expression
 
 
 @six.add_metaclass(inheritdocstring)
 class Configuration(interface.MathematicalProgrammingConfiguration):
-    def __init__(self, lp_method='primal', presolve=False, verbosity=0, timeout=None, *args, **kwargs):
+    def __init__(self, lp_method='primal', qp_method='primal', presolve=False,
+                 verbosity=0, timeout=None, *args, **kwargs):
         super(Configuration, self).__init__(*args, **kwargs)
         self.lp_method = lp_method
+        self.qp_method = qp_method
         self.presolve = presolve
         self.verbosity = verbosity
         self.timeout = timeout
@@ -402,6 +417,17 @@ def lp_method(self, value):
         if getattr(self, "problem", None) is not None:
             self.problem.problem.params.Method = _LP_METHODS[value]
 
+    @property
+    def qp_method(self):
+        return _REVERSE_QP_METHODS[self.problem.problem.params.Method]
+
+    @qp_method.setter
+    def qp_method(self, value):
+        if value not in _QP_METHODS:
+            raise ValueError("Invalid LP method. Please choose among: " + str(list(_QP_METHODS)))
+        if getattr(self, "problem", None) is not None:
+            self.problem.problem.params.Method = _QP_METHODS[value]
+
     @property
     def presolve(self):
         return self._presolve
@@ -525,13 +551,24 @@ def __init__(self, problem=None, *args, **kwargs):
             super(Model, self)._add_constraints(constraints, sloppy=True)
 
             gurobi_objective = self.problem.getObjective()
+            if self.problem.IsQP:
+                quadratic_expression = symbolics.add(
+                    [symbolics.Real(gurobi_objective.getCoeff(i)) *
+                     self.variables[gurobi_objective.getVar1(i).VarName] *
+                     self.variables[gurobi_objective.getVar2(i).VarName]
+                     for i in range(gurobi_objective.size())])
+                linear_objective = gurobi_objective.getLinExpr()
+            else:
+                quadratic_expression = symbolics.Real(0.0)
+                linear_objective = gurobi_objective
             linear_expression = symbolics.add(
-                [symbolics.Real(gurobi_objective.getCoeff(i)) * self.variables[gurobi_objective.getVar(i).VarName]
-                 for i in range(gurobi_objective.size())]
+                [symbolics.Real(linear_objective.getCoeff(i)) *
+                 self.variables[linear_objective.getVar(i).VarName]
+                 for i in range(linear_objective.size())]
             )
 
             self._objective = Objective(
-                linear_expression,
+                quadratic_expression + linear_expression,
                 problem=self,
                 direction={1: 'min', -1: 'max'}[self.problem.getAttr('ModelSense')]
             )
@@ -604,7 +641,7 @@ def objective(self, value):
                 var1, var2 = key
                 var1 = self.problem.getVarByName(var1.name)
                 var2 = self.problem.getVarByName(var2.name)
-                grb_terms.append(coef, var1, var2)
+                grb_terms.append(coef * var1 * var2)
 
         grb_expression = gurobipy.quicksum(grb_terms)
 

optlang/tests/abstract_test_cases.py

@@ -688,13 +688,13 @@ def test_objective_set_linear_coefficients(self):
         self.assertAlmostEqual(y.primal, 0)
 
         obj.set_linear_coefficients({y: 1})
-        self.assertEqual((obj.expression - (x + y)).expand(), 0)
+        self.assertEqual(float((obj.expression - (x + y)).expand()), 0.0)
         self.assertEqual(model.optimize(), optlang.interface.OPTIMAL)
         self.assertAlmostEqual(x.primal, 2)
         self.assertAlmostEqual(y.primal, 2)
 
         obj.set_linear_coefficients({x: 0})
-        self.assertEqual((obj.expression - y).expand(), 0)
+        self.assertEqual(float((obj.expression - y).expand()), 0.0)
         self.assertEqual(model.optimize(), optlang.interface.OPTIMAL)
         self.assertAlmostEqual(x.primal, 0)
         self.assertAlmostEqual(y.primal, 3)

optlang/tests/test_gurobi_interface.py

@@ -25,12 +25,15 @@ def test_fail(self):
     import pickle
 
     from optlang.gurobi_interface import Variable, Constraint, Model, Objective
+    from gurobipy import GurobiError
     from optlang.tests import abstract_test_cases
     from optlang import gurobi_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')
+    CONVEX_QP_PATH = os.path.join(os.path.dirname(__file__), 'data/qplib_3256.lp')
+    NONCONVEX_QP_PATH = os.path.join(os.path.dirname(__file__), 'data/qplib_1832.lp')
 
 
     class VariableTestCase(abstract_test_cases.AbstractVariableTestCase):
@@ -438,5 +441,75 @@ def test_set_linear_coefficients_constraint(self):
                 if col_name == 'R_Biomass_Ecoli_core_w_GAM':
                     self.assertEqual(row.getCoeff(i), 666.)
 
+
+    class QuadraticProgrammingTestCase(abstract_test_cases.AbstractQuadraticProgrammingTestCase):
+            def setUp(self):
+                self.model = Model()
+                self.x1 = Variable("x1", lb=0)
+                self.x2 = Variable("x2", lb=0)
+                self.c1 = Constraint(self.x1 + self.x2, lb=1)
+                self.model.add([self.x1, self.x2, self.c1])
+
+            def test_convex_obj(self):
+                model = self.model
+                obj = Objective(self.x1 ** 2 + self.x2 ** 2, direction="min")
+                model.objective = obj
+                model.optimize()
+                self.assertAlmostEqual(model.objective.value, 0.5)
+                self.assertAlmostEqual(self.x1.primal, 0.5)
+                self.assertAlmostEqual(self.x2.primal, 0.5)
+
+                obj_2 = Objective(self.x1, direction="min")
+                model.objective = obj_2
+                model.optimize()
+                self.assertAlmostEqual(model.objective.value, 0.0)
+                self.assertAlmostEqual(self.x1.primal, 0.0)
+                self.assertGreaterEqual(self.x2.primal, 1.0)
+
+            def test_non_convex_obj(self):
+                model = self.model
+                obj = Objective(self.x1 * self.x2, direction="min")
+                model.objective = obj
+                self.assertRaises(GurobiError, model.optimize)
+
+                obj_2 = Objective(self.x1, direction="min")
+                model.objective = obj_2
+                model.optimize()
+                self.assertAlmostEqual(model.objective.value, 0.0)
+                self.assertAlmostEqual(self.x1.primal, 0.0)
+                self.assertGreaterEqual(self.x2.primal, 1.0)
+
+            def test_qp_convex(self):
+                model = Model(problem=gurobipy.read(CONVEX_QP_PATH))
+                self.assertEqual(len(model.variables), 651)
+                self.assertEqual(len(model.constraints), 501)
+                for constraint in model.constraints:
+                    self.assertTrue(constraint.is_Linear, "%s should be linear" % (str(constraint.expression)))
+                    self.assertFalse(constraint.is_Quadratic, "%s should not be quadratic" % (str(constraint.expression)))
+
+                self.assertTrue(model.objective.is_Quadratic, "objective should be quadratic")
+                self.assertFalse(model.objective.is_Linear, "objective should not be linear")
+
+                model.optimize()
+                self.assertAlmostEqual(model.objective.value, 32.2291282)
+
+            def test_qp_non_convex(self):
+                model = Model(problem=gurobipy.read(NONCONVEX_QP_PATH))
+                self.assertEqual(len(model.variables), 31)
+                self.assertEqual(len(model.constraints), 1)
+                for constraint in model.constraints:
+                    self.assertTrue(constraint.is_Linear, "%s should be linear" % (str(constraint.expression)))
+                    self.assertFalse(constraint.is_Quadratic, "%s should not be quadratic" % (str(constraint.expression)))
+
+                self.assertTrue(model.objective.is_Quadratic, "objective should be quadratic")
+                self.assertFalse(model.objective.is_Linear, "objective should not be linear")
+
+                self.assertRaises(GurobiError, model.optimize)
+
+            def test_quadratic_objective_expression(self):
+                objective = Objective(self.x1 ** 2 + self.x2 ** 2, direction="min")
+                self.model.objective = objective
+                self.assertEqual((self.model.objective.expression - (self.x1 ** 2 + self.x2 ** 2)).simplify(), 0)
+
 if __name__ == '__main__':
     nose.runmodule()

setup.cfg

@@ -15,4 +15,4 @@ exclude = tests/*,versioneer.py,__init__.py,_version.py,plotting_old.py,cplex/*,
 max-complexity = 20
 
 [bdist_wheel]
-universal = 1
+universal = 1