Various bug fixes (#126)

- Make gurobi models picklable
- Fix gurobi infinite recursion
- Fix gurobi constraint primal calculation
- Fix a bug where cplex took forever to set a large objective (e.g. 10000 vars)
- Fix some issues with symengine compatibility

optlang/cplex_interface.py

@@ -711,7 +711,9 @@ def objective(self, value):
         if self._objective is not None:  # Reset previous objective
             variables = self.objective.variables
             if len(variables) > 0:
-                self.problem.objective.set_linear([(variable.name, 0.) for variable in variables])
+                name_list = [var.name for var in variables]
+                index_dict = {n: i for n, i in zip(name_list, self._get_variable_indices(name_list))}
+                self.problem.objective.set_linear([(index_dict[variable.name], 0.) for variable in variables])
             if self.problem.objective.get_num_quadratic_variables() > 0:
                 self.problem.objective.set_quadratic([0. for _ in range(self.problem.variables.get_num())])
         super(Model, self.__class__).objective.fset(self, value)
@@ -721,7 +723,9 @@ def objective(self, value):
         # self.problem.objective.set_offset(float(offset)) # Not available prior to 12.6.2
         self._objective_offset = offset
         if linear_coefficients:
-            self.problem.objective.set_linear([var.name, float(coef)] for var, coef in linear_coefficients.items())
+            name_list = [var.name for var in linear_coefficients]
+            index_dict = {n: i for n, i in zip(name_list, self._get_variable_indices(name_list))}
+            self.problem.objective.set_linear([index_dict[var.name], float(coef)] for var, coef in linear_coefficients.items())
 
         for key, coef in quadratic_coeffients.items():
             if len(key) == 1:
@@ -901,3 +905,16 @@ def _get_quadratic_expression(self, quadratic=None):
                 else:
                     pass  # Only look at upper triangle
         return add(terms)
+
+    def _get_variable_indices(self, names):
+        # Cplex does not keep an index of variable names
+        # Getting indices thus takes roughly quadratic time
+        # If many indices are required an alternate and faster method is used, where the full name list must only
+        # be traversed once
+        if len(names) < 1000:
+            return self.problem.variables.get_indices(names)
+        else:
+            name_set = set(names)
+            all_names = self.problem.variables.get_names()
+            indices = {n: i for i, n in enumerate(all_names) if n in name_set}
+            return [indices[n] for n in names]

optlang/duality.py

@@ -13,8 +13,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-from sympy import Add, Mul
-
+import optlang
 
 # This function is very complex. Should maybe be refactored
 def convert_linear_problem_to_dual(model, sloppy=False, infinity=None, maintain_standard_form=True, prefix="dual_", dual_model=None):  # NOQA
@@ -84,6 +83,8 @@ def convert_linear_problem_to_dual(model, sloppy=False, infinity=None, maintain_
             if constraint.lb != 0:
                 dual_objective[const_var] = sign * constraint.lb
             for variable, coef in constraint.expression.as_coefficients_dict().items():
+                if variable == 1:
+                    continue
                 coefficients.setdefault(variable.name, {})[const_var] = sign * coef
         else:
             if constraint.lb is not None:
@@ -105,6 +106,8 @@ def convert_linear_problem_to_dual(model, sloppy=False, infinity=None, maintain_
                 coefficients_dict = {constraint.expression.args[1]: constraint.expression.args[0]}
 
             for variable, coef in coefficients_dict.items():
+                if variable == 1:
+                    continue
                 if constraint.lb is not None:
                     coefficients.setdefault(variable.name, {})[lb_var] = -sign * coef
                 if constraint.ub is not None:
@@ -131,7 +134,7 @@ def convert_linear_problem_to_dual(model, sloppy=False, infinity=None, maintain_
     # Add dual constraints from primal objective
     primal_objective_dict = model.objective.expression.as_coefficients_dict()
     for variable in model.variables:
-        expr = Add(*((coef * dual_var) for dual_var, coef in coefficients[variable.name].items()))
+        expr = optlang.symbolics.add([(coef * dual_var) for dual_var, coef in coefficients[variable.name].items()])
         obj_coef = primal_objective_dict[variable]
         if maximization:
             const = model.interface.Constraint(expr, lb=obj_coef, name=prefix + variable.name)
@@ -140,7 +143,7 @@ def convert_linear_problem_to_dual(model, sloppy=False, infinity=None, maintain_
         dual_model.add(const)
 
     # Make dual objective
-    expr = Add(*((coef * dual_var) for dual_var, coef in dual_objective.items() if coef != 0))
+    expr = optlang.symbolics.add([(coef * dual_var) for dual_var, coef in dual_objective.items() if coef != 0])
     if maximization:
         objective = model.interface.Objective(expr, direction="min")
     else:

optlang/glpk_interface.py

@@ -50,7 +50,7 @@
     glp_get_mat_row, glp_get_row_ub, glp_get_row_type, glp_get_row_lb, glp_get_row_name, glp_get_obj_coef, \
     glp_get_obj_dir, glp_scale_prob, GLP_SF_AUTO, glp_get_num_int, glp_get_num_bin, glp_mip_col_val, \
     glp_mip_obj_val, glp_mip_status, GLP_ETMLIM, glp_adv_basis, glp_read_lp, glp_mip_row_val, \
-    get_col_primals, get_col_duals, get_row_primals, get_row_duals
+    get_col_primals, get_col_duals, get_row_primals, get_row_duals, glp_delete_prob
 
 
 
@@ -579,6 +579,10 @@ def __setstate__(self, repr_dict):
         if repr_dict['glpk_status'] == 'optimal':
             self.optimize()  # since the start is an optimal solution, nothing will happen here
 
+    # def __del__(self):  # To make sure that the glpk problem is deleted when this is garbage collected
+    # Gotcha: When objects with a __del__ method are part of a referencing cycle, the entire cycle is never automatically garbage collected
+    #     glp_delete_prob(self.problem)
+
     @property
     def objective(self):
         return self._objective

optlang/gurobi_interface.py

@@ -27,6 +27,7 @@
 
 import os
 import six
+from functools import partial
 from optlang import interface
 from optlang.util import inheritdocstring, TemporaryFilename
 from optlang.expression_parsing import parse_optimization_expression
@@ -236,7 +237,14 @@ def problem(self, value):
     @property
     def primal(self):
         if self.problem is not None:
-            return self._internal_constraint.Slack
+            aux_var = self.problem.problem.getVarByName(self._internal_constraint.getAttr('ConstrName') + '_aux')
+            if aux_var is not None:
+                aux_val = aux_var.X
+            else:
+                aux_val = 0
+            return (self._internal_constraint.RHS -
+                    self._internal_constraint.Slack +
+                    aux_val)
             # return self._round_primal_to_bounds(primal_from_solver)  # Test assertions fail
             # return primal_from_solver
         else:
@@ -325,9 +333,10 @@ def __init__(self, expression, sloppy=False, *args, **kwargs):
     @property
     def value(self):
         if getattr(self, "problem", None) is not None:
+            gurobi_problem = self.problem.problem
             try:
-                return self.problem.problem.getAttr("ObjVal") + getattr(self.problem, "_objective_offset", 0)
-            except gurobipy.GurobiError:  # TODO: let this check the actual error message
+                return gurobi_problem.getAttr("ObjVal") + getattr(self.problem, "_objective_offset", 0)
+            except (gurobipy.GurobiError, AttributeError):  # TODO: let this check the actual error message
                 return None
         else:
             return None
@@ -359,8 +368,9 @@ 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) * self.problem.variables[grb_obj.getVar(i).getAttr('VarName')])
+                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)
@@ -385,20 +395,22 @@ def lp_method(self):
     def lp_method(self, value):
         if value not in _LP_METHODS:
             raise ValueError("Invalid LP method. Please choose among: " + str(list(_LP_METHODS)))
-        self.problem.problem.params.Method = _LP_METHODS[value]
+        if getattr(self, "problem", None) is not None:
+            self.problem.problem.params.Method = _LP_METHODS[value]
 
     @property
     def presolve(self):
         return self._presolve
 
     @presolve.setter
     def presolve(self, value):
-        if value is True:
-            self.problem.problem.params.Presolve = 2
-        elif value is False:
-            self.problem.problem.params.Presolve = 0
-        else:
-            self.problem.problem.params.Presolve = -1
+        if getattr(self, "problem", None) is not None:
+            if value is True:
+                self.problem.problem.params.Presolve = 2
+            elif value is False:
+                self.problem.problem.params.Presolve = 0
+            else:
+                self.problem.problem.params.Presolve = -1
         self._presolve = value
 
     @property
@@ -427,20 +439,37 @@ def timeout(self, value):
                 self.problem.problem.params.TimeLimit = value
         self._timeout = value
 
+    def __getstate__(self):
+        return {'presolve': self.presolve, 'verbosity': self.verbosity, 'timeout': self.timeout}
+
+    def __setstate__(self, state):
+        self.__init__()
+        for key, val in six.iteritems(state):
+            setattr(self, key, val)
+
+    def _get_feasibility(self):
+        return getattr(self.problem.problem.params, "FeasibilityTol")
+
+    def _set_feasibility(self, value):
+        return setattr(self.problem.problem.params, "FeasibilityTol", value)
+
+    def _get_optimality(self):
+        return getattr(self.problem.problem.params, "OptimalityTol")
+
+    def _set_optimality(self, value):
+        return setattr(self.problem.problem.params, "OptimalityTol", value)
+
+    def _get_integrality(self):
+        return getattr(self.problem.problem.params, "IntFeasTol")
+
+    def _set_integrality(self, value):
+        return setattr(self.problem.problem.params, "IntFeasTol", value)
+
     def _tolerance_functions(self):
         return {
-            "feasibility": (
-                lambda: self.problem.problem.params.FeasibilityTol,
-                lambda x: setattr(self.problem.problem.params, "FeasibilityTol", x)
-            ),
-            "optimality": (
-                lambda: self.problem.problem.params.OptimalityTol,
-                lambda x: setattr(self.problem.problem.params, "OptimalityTol", x)
-            ),
-            "integrality": (
-                lambda: self.problem.problem.params.IntFeasTol,
-                lambda x: setattr(self.problem.problem.params, "IntFeasTol", x)
-            )
+            "feasibility": (self._get_feasibility, self._set_feasibility),
+            "optimality": (self._get_optimality, self._set_optimality),
+            "integrality": (self._get_integrality, self._set_integrality)
         }
 
 

optlang/symbolics.py

@@ -72,6 +72,9 @@ def __repr__(self):
         def __str__(self):
             return self._name
 
+        def __getnewargs__(self):
+            return (self._name, {})
+
     def add(*args):
         if len(args) == 1:
             args = args[0]

optlang/tests/abstract_test_cases.py

@@ -731,7 +731,7 @@ def test_objective_expression_includes_constant(self):
         objective = self.model.objective
         self.model.objective = self.interface.Objective(objective.expression + 3)
         self.model.update()
-        self.assertEqual((self.model.objective.expression - (objective.expression + 3)).expand(), 0)
+        self.assertEqual((self.model.objective.expression - (objective.expression + 3.)).expand(), 0.)
 
     def test_is_integer(self):
         model = self.model
@@ -809,6 +809,19 @@ def test_integer_batch_duals(self):
         self.assertEqual(model.reduced_costs[x.name], 0)
         self.assertEqual(model.shadow_prices[c.name], 1)
 
+    def test_large_objective(self):
+        model = self.interface.Model()
+        model.add([self.interface.Variable(str(i), lb=1) for i in range(1100)])
+        model.optimize()
+
+        obj = self.interface.Objective(
+            optlang.symbolics.add([optlang.symbolics.mul((optlang.symbolics.One, v)) for v in model.variables]),
+            direction="min"
+        )
+        model.objective = obj
+        model.optimize()
+        self.assertAlmostEqual(model.objective.value, len(model.variables))
+
 
 @six.add_metaclass(abc.ABCMeta)
 class AbstractConfigurationTestCase(unittest.TestCase):

optlang/tests/test_expression_parsing.py

@@ -7,6 +7,20 @@
 import unittest
 
 
+def _quad_terms_to_expression(terms):
+    args = []
+    for term, coef in terms.items():
+        term = list(term)
+        args.append(coef * term[0] * term[-1])
+    return sum(args)
+
+
+def _compare_term_dicts(test_case, dict1, dict2):
+    for term, coef1 in dict1.items():
+        coef2 = dict2[term]
+        test_case.assertEqual(coef1 - coef2, 0)
+
+
 class ExpressionParsingTestCase(unittest.TestCase):
     def setUp(self):
         self.vars = [Variable(name) for name in ["x", "y", "z", "v", "w"]]
@@ -22,8 +36,8 @@ def test_parse_linear_expression(self):
 
         self.assertEqual(offset_const, 0)
         self.assertEqual(offset_obj, offset)
-        self.assertEqual(linear_terms_const, target)
-        self.assertEqual(linear_terms_obj, target)
+        _compare_term_dicts(self, linear_terms_const, target)
+        _compare_term_dicts(self, linear_terms_obj, target)
         self.assertEqual(quad_terms_const, {})
         self.assertEqual(quad_terms_obj, {})
 
@@ -41,8 +55,10 @@ def test_parse_quadratic_expression(self):
         self.assertEqual(offset_obj, offset)
         self.assertEqual(linear_terms_const, {})
         self.assertEqual(linear_terms_obj, {})
-        self.assertEqual(quad_terms_const, target)
-        self.assertEqual(quad_terms_obj, target)
+        _compare_term_dicts(self, quad_terms_const, target)
+        _compare_term_dicts(self, quad_terms_obj, target)
+        self.assertEqual((_quad_terms_to_expression(quad_terms_obj) - (expr - offset)).expand(), 0)
+        self.assertEqual((_quad_terms_to_expression(quad_terms_const) - (expr - offset)).expand(), 0)
 
     def test_parse_non_expanded_quadratic_expression(self):
         x, y, z = self.vars[:3]
@@ -57,7 +73,7 @@ def test_parse_non_expanded_quadratic_expression(self):
 
         self.assertEqual(offset_const, -4)
         self.assertEqual(offset_obj, -4 + offset)
-        self.assertEqual(linear_terms_const, linear_target)
-        self.assertEqual(linear_terms_obj, linear_target)
-        self.assertEqual(quad_terms_const, target)
-        self.assertEqual(quad_terms_obj, target)
+        _compare_term_dicts(self, linear_terms_const, linear_target)
+        _compare_term_dicts(self, linear_terms_obj, linear_target)
+        _compare_term_dicts(self, quad_terms_const, target)
+        _compare_term_dicts(self, quad_terms_obj, target)

optlang/tests/test_scipy_interface.py

@@ -315,3 +315,6 @@ def test_integer_constraint_dual(self):
 
     def test_integer_batch_duals(self):
         self.skipTest("No duals with scipy")
+
+    def test_large_objective(self):
+        self.skipTest("Quite slow and not necessary")