1.2.2 (#129)

- Fixed a couple of gurobi interface bugs
- Fixed a bug where cplex took forever to set a large objective (e.g. 10000 vars)
- Experimental support for symengine (requires Python >= 3.5, install with `conda install -c symengine python-symengine`)

optlang/cplex_interface.py

@@ -29,22 +29,18 @@
 import os
 from six.moves import StringIO
 
-import sympy
-from sympy.core.add import _unevaluated_Add
-from sympy.core.mul import _unevaluated_Mul
-from sympy.core.singleton import S
 import cplex
 from cplex.exceptions import CplexSolverError
 
 from optlang import interface
+from optlang import symbolics
 from optlang.util import inheritdocstring, TemporaryFilename
 from optlang.expression_parsing import parse_optimization_expression
 from optlang.exceptions import SolverError
 
 log = logging.getLogger(__name__)
 
-Zero = S.Zero
-One = S.One
+from optlang.symbolics import add, mul, One, Zero
 
 _STATUS_MAP = {
     'MIP_abort_feasible': interface.ABORTED,
@@ -242,8 +238,8 @@ def _get_expression(self):
             cplex_problem = self.problem.problem
             cplex_row = cplex_problem.linear_constraints.get_rows(self.name)
             variables = self.problem._variables
-            expression = sympy.Add._from_args(
-                [sympy.Mul._from_args((sympy.RealNumber(cplex_row.val[i]), variables[ind])) for i, ind in
+            expression = add(
+                [mul((symbolics.Real(cplex_row.val[i]), variables[ind])) for i, ind in
                  enumerate(cplex_row.ind)])
             self._expression = expression
         return self._expression
@@ -364,7 +360,7 @@ def _get_expression(self):
         if self.problem is not None and self._expression_expired and len(self.problem._variables) > 0:
             cplex_problem = self.problem.problem
             coeffs = cplex_problem.objective.get_linear()
-            expression = sympy.Add._from_args([coeff * var for coeff, var in zip(coeffs, self.problem._variables) if coeff != 0.])
+            expression = add([coeff * var for coeff, var in zip(coeffs, self.problem._variables) if coeff != 0.])
             if cplex_problem.objective.get_num_quadratic_nonzeros() > 0:
                 expression += self.problem._get_quadratic_expression(cplex_problem.objective.get_quadratic())
             self._expression = expression + getattr(self.problem, "_objective_offset", 0)
@@ -616,9 +612,9 @@ def __init__(self, problem=None, *args, **kwargs):
                 # lhs = _unevaluated_Add(*[val * variables[i - 1] for i, val in zip(row.ind, row.val)])
                 lhs = 0
                 if isinstance(lhs, int):
-                    lhs = sympy.Integer(lhs)
+                    lhs = symbolics.Integer(lhs)
                 elif isinstance(lhs, float):
-                    lhs = sympy.RealNumber(lhs)
+                    lhs = symbolics.Real(lhs)
                 if sense == 'E':
                     constr = Constraint(lhs, lb=rhs, ub=rhs, name=name, problem=self)
                 elif sense == 'G':
@@ -650,10 +646,9 @@ def __init__(self, problem=None, *args, **kwargs):
                 if 'CPLEX Error  1219:' not in str(e):
                     raise e
             else:
-                linear_expression = _unevaluated_Add(
-                    *[_unevaluated_Mul(sympy.RealNumber(coeff), variables[index]) for index, coeff in
-                      enumerate(self.problem.objective.get_linear()) if coeff != 0.])
-
+                linear_expression = add(
+                    [mul(symbolics.Real(coeff), variables[index]) for index, coeff in enumerate(self.problem.objective.get_linear()) if coeff != 0.]
+                )
                 try:
                     quadratic = self.problem.objective.get_quadratic()
                 except IndexError:
@@ -716,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)
@@ -726,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:
@@ -905,4 +904,17 @@ def _get_quadratic_expression(self, quadratic=None):
                     terms.append(0.5 * val * self._variables[i_name] ** 2)
                 else:
                     pass  # Only look at upper triangle
-        return _unevaluated_Add(*terms)
+        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:  # pragma: no cover  # For symengine
+                    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:  # pragma: no cover  # For symengine
+                    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/expression_parsing.py

@@ -13,6 +13,8 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+from optlang.symbolics import One
+
 
 def parse_optimization_expression(obj, linear=True, quadratic=False, expression=None, **kwargs):
     """
@@ -79,9 +81,10 @@ def _parse_linear_expression(expression, expanded=False, **kwargs):
         coefficients = {}
     else:
         raise ValueError("Expression {} seems to be invalid".format(expression))
+
     for var in coefficients:
         if not (var.is_Symbol):
-            if var == 1:
+            if var == One:
                 constant = var
                 offset = float(coefficients[var])
             elif expanded:

optlang/glpk_interface.py

@@ -30,13 +30,11 @@
 
 import os
 import six
-import sympy
-from sympy.core.add import _unevaluated_Add
-from sympy.core.mul import _unevaluated_Mul
 
 from optlang.util import inheritdocstring, TemporaryFilename
 from optlang.expression_parsing import parse_optimization_expression
 from optlang import interface
+from optlang import symbolics
 
 log = logging.getLogger(__name__)
 
@@ -52,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
 
 
 
@@ -165,8 +163,8 @@ def _get_expression(self):
             nnz = glp_get_mat_row(self.problem.problem, self._index, ia, da)
             constraint_variables = [self.problem._variables[glp_get_col_name(self.problem.problem, ia[i])] for i in
                                     range(1, nnz + 1)]
-            expression = sympy.Add._from_args(
-                [sympy.Mul._from_args((sympy.RealNumber(da[i]), constraint_variables[i - 1])) for i in
+            expression = symbolics.add(
+                [symbolics.mul((symbolics.Real(da[i]), constraint_variables[i - 1])) for i in
                  range(1, nnz + 1)])
             self._expression = expression
         return self._expression
@@ -308,9 +306,9 @@ def term_generator():
                 for index in range(1, glp_get_num_cols(self.problem.problem) + 1):
                     coeff = glp_get_obj_coef(self.problem.problem, index)
                     if coeff != 0.:
-                        yield (sympy.RealNumber(coeff), variables[index - 1])
+                        yield (symbolics.Real(coeff), variables[index - 1])
 
-            expression = sympy.Add._from_args([sympy.Mul._from_args(term) for term in term_generator()])
+            expression = symbolics.add([symbolics.mul(term) for term in term_generator()])
             self._expression = expression + getattr(self.problem, "_objective_offset", 0)
             self._expression_expired = False
         return self._expression
@@ -531,9 +529,9 @@ def __init__(self, problem=None, *args, **kwargs):
                     )
                     log.exception()
                 if isinstance(lhs, int):
-                    lhs = sympy.Integer(lhs)
+                    lhs = symbolics.Integer(lhs)
                 elif isinstance(lhs, float):
-                    lhs = sympy.RealNumber(lhs)
+                    lhs = symbolics.Real(lhs)
                 constraint_id = glp_get_row_name(self.problem, j)
                 for variable in constraint_variables:
                     try:
@@ -551,9 +549,10 @@ def __init__(self, problem=None, *args, **kwargs):
                 for index in range(1, glp_get_num_cols(problem) + 1)
             )
             self._objective = Objective(
-                _unevaluated_Add(
-                    *[_unevaluated_Mul(sympy.RealNumber(term[0]), term[1]) for term in term_generator if
-                      term[0] != 0.]),
+                symbolics.add(
+                    [symbolics.mul((symbolics.Real(term[0]), term[1])) for term in term_generator if
+                     term[0] != 0.]
+                ),
                 problem=self,
                 direction={GLP_MIN: 'min', GLP_MAX: 'max'}[glp_get_obj_dir(self.problem)])
         glp_scale_prob(self.problem, GLP_SF_AUTO)
@@ -580,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