Make cobrapy compatible with symengine (#613)

* Add symengine support

cobra/core/model.py

@@ -9,10 +9,9 @@
 from warnings import warn
 
 import optlang
+from optlang.symbolics import Basic, Zero
 import six
-import sympy
 from six import iteritems, string_types
-from sympy import S
 
 from cobra.core.dictlist import DictList
 from cobra.core.object import Object
@@ -101,7 +100,7 @@ def __init__(self, id_or_model=None, name=None):
             # with older cobrapy pickles?
             interface = solvers[get_solver_name()]
             self._solver = interface.Model()
-            self._solver.objective = interface.Objective(S.Zero)
+            self._solver.objective = interface.Objective(Zero)
             self._populate_solver(self.reactions, self.metabolites)
 
     @property
@@ -368,7 +367,7 @@ def add_metabolites(self, metabolite_list):
         for met in metabolite_list:
             if met.id not in self.constraints:
                 constraint = self.problem.Constraint(
-                    S.Zero, name=met.id, lb=0, ub=0)
+                    Zero, name=met.id, lb=0, ub=0)
                 to_add += [constraint]
 
         self.add_cons_vars(to_add)
@@ -743,7 +742,7 @@ def _populate_solver(self, reaction_list, metabolite_list=None):
         if metabolite_list is not None:
             for met in metabolite_list:
                 to_add += [self.problem.Constraint(
-                    S.Zero, name=met.id, lb=0, ub=0)]
+                    Zero, name=met.id, lb=0, ub=0)]
         self.add_cons_vars(to_add)
 
         for reaction in reaction_list:
@@ -764,7 +763,7 @@ def _populate_solver(self, reaction_list, metabolite_list=None):
                     constraint = self.constraints[metabolite.id]
                 else:
                     constraint = self.problem.Constraint(
-                        S.Zero,
+                        Zero,
                         name=metabolite.id,
                         lb=0, ub=0)
                     self.add_cons_vars(constraint, sloppy=True)
@@ -943,7 +942,7 @@ def objective(self):
 
     @objective.setter
     def objective(self, value):
-        if isinstance(value, sympy.Basic):
+        if isinstance(value, Basic):
             value = self.problem.Objective(value, sloppy=False)
         if not isinstance(value, (dict, optlang.interface.Objective)):
             try:

cobra/core/reaction.py

@@ -118,7 +118,7 @@ def flux_expression(self):
         Returns
         -------
         sympy expression
-            The expression represeenting the the forward flux (if associated
+            The expression representing the the forward flux (if associated
             with model), otherwise None. Representing the net flux if
             model.reversible_encoding == 'unsplit' or None if reaction is
             not associated with a model """

cobra/flux_analysis/gapfilling.py

@@ -2,7 +2,7 @@
 
 from __future__ import absolute_import
 
-from sympy import Add
+from optlang.symbolics import Add
 from warnings import warn
 
 from optlang.interface import OPTIMAL

cobra/flux_analysis/loopless.py

@@ -6,7 +6,7 @@
 
 import numpy
 from six import iteritems
-from sympy.core.singleton import S
+from optlang.symbolics import Zero
 
 from cobra.core import Metabolite, Reaction, get_solution
 from cobra.util import (linear_reaction_coefficients,
@@ -72,7 +72,7 @@ def add_loopless(model, zero_cutoff=1e-12):
     # Add nullspace constraints for G_i
     for i, row in enumerate(n_int):
         name = "nullspace_constraint_" + str(i)
-        nullspace_constraint = prob.Constraint(S.Zero, lb=0, ub=0, name=name)
+        nullspace_constraint = prob.Constraint(Zero, lb=0, ub=0, name=name)
         model.add_cons_vars([nullspace_constraint])
         coefs = {model.variables[
                  "delta_g_" + model.reactions[ridx].id]: row[i]
@@ -83,7 +83,7 @@ def add_loopless(model, zero_cutoff=1e-12):
 
 def _add_cycle_free(model, fluxes):
     """Add constraints for CycleFreeFlux."""
-    model.objective = S.Zero
+    model.objective = Zero
     for rxn in model.reactions:
         flux = fluxes[rxn.id]
         if rxn.boundary:

cobra/flux_analysis/moma.py

@@ -5,7 +5,7 @@
 from __future__ import absolute_import
 
 from scipy.sparse import dok_matrix
-from sympy.core.singleton import S
+from optlang.symbolics import Zero
 
 import cobra.util.solver as sutil
 from cobra.solvers import get_solver_name, solver_dict
@@ -77,7 +77,7 @@ def add_moma(model, solution=None, linear=False):
     c = prob.Constraint(model.solver.objective.expression - v,
                         lb=0.0, ub=0.0, name="moma_old_objective_constraint")
     to_add = [v, c]
-    new_obj = S.Zero
+    new_obj = Zero
     for r in model.reactions:
         flux = solution.fluxes[r.id]
         if linear:

cobra/flux_analysis/parsimonious.py

@@ -6,16 +6,14 @@
 from warnings import warn
 from itertools import chain
 
-import sympy
+from optlang.symbolics import Zero
 
 from cobra.util import solver as sutil
 from cobra.manipulation.modify import (
     convert_to_irreversible, revert_to_reversible)
 from cobra.util import linear_reaction_coefficients, set_objective
 from cobra.core.solution import get_solution
 
-add = sympy.Add._from_args
-mul = sympy.Mul._from_args
 LOGGER = logging.getLogger(__name__)
 
 
@@ -128,11 +126,10 @@ def add_pfba(model, objective=None, fraction_of_optimum=1.0):
     reaction_variables = ((rxn.forward_variable, rxn.reverse_variable)
                           for rxn in model.reactions)
     variables = chain(*reaction_variables)
-    pfba_objective = model.problem.Objective(add(
-        [mul((sympy.singleton.S.One, variable))
-         for variable in variables]), direction='min', sloppy=True,
+    model.objective = model.problem.Objective(
+        Zero, direction='min', sloppy=True,
         name="_pfba_objective")
-    set_objective(model, pfba_objective)
+    model.objective.set_linear_coefficients(dict.fromkeys(variables, 1.0))
 
 
 def _pfba_optlang(model, objective=None, reactions=None,

cobra/flux_analysis/sampling.py

@@ -17,7 +17,7 @@
 import numpy as np
 import pandas
 from optlang.interface import OPTIMAL
-from sympy.core.singleton import S
+from optlang.symbolics import Zero
 from cobra.util import (create_stoichiometric_matrix, constraint_matrices,
                         nullspace)
 
@@ -254,7 +254,7 @@ def generate_fva_warmup(self):
         self.n_warmup = 0
         idx = np.hstack([self.fwd_idx, self.rev_idx])
         self.warmup = np.zeros((len(idx), len(self.model.variables)))
-        self.model.objective = S.Zero
+        self.model.objective = Zero
         self.model.objective.direction = "max"
         variables = self.model.variables
         for i in idx:

cobra/flux_analysis/variability.py

@@ -3,7 +3,7 @@
 from __future__ import absolute_import
 
 import pandas
-from sympy.core.singleton import S
+from optlang.symbolics import Zero
 from warnings import warn
 from itertools import chain
 
@@ -211,7 +211,7 @@ def _fva_optlang(model, reaction_list, fraction, loopless, pfba_factor):
                     name="flux_sum_constraint")
             m.add_cons_vars([flux_sum, flux_sum_constraint])
 
-        m.objective = S.Zero  # This will trigger the reset as well
+        m.objective = Zero  # This will trigger the reset as well
         for what in ("minimum", "maximum"):
             sense = "min" if what == "minimum" else "max"
             for rxn in reaction_list:

cobra/io/sbml3.py

@@ -48,7 +48,7 @@
     from cobra.io.sbml import write_cobra_model_to_sbml_file as write_sbml2
 
 try:
-    from sympy import Basic
+    from optlang.symbolics import Basic
 except ImportError:
     class Basic:
         pass

cobra/solvers/cglpk.pyx

@@ -14,7 +14,7 @@ from warnings import warn as _warn
 
 from six import StringIO, iteritems
 try:
-    from sympy import Basic, Number
+    from optlang.symbolics import Basic, Number
 except:
     class Basic:
         pass

cobra/solvers/cplex_solver.py

@@ -15,7 +15,7 @@
 from cobra.core.solution import LegacySolution
 
 try:
-    from sympy import Basic, Number
+    from optlang.symbolics import Basic, Number
 except:
     class Basic:
         pass

cobra/solvers/gurobi_solver.py

@@ -40,7 +40,7 @@ def test_import():
 
 
 try:
-    from sympy import Basic, Number
+    from optlang.symbolics import Basic, Number
 except:
     class Basic:
         pass

cobra/test/test_model.py

@@ -8,7 +8,7 @@
 from math import isnan
 import pytest
 import pandas as pd
-from sympy import S
+from optlang.symbolics import Zero
 
 import cobra.util.solver as su
 from cobra.core import Metabolite, Model, Reaction
@@ -754,8 +754,9 @@ def test_merge_models(self, model, tiny_toy_model):
                 assert 'constraint' in model.constraints
                 assert 'foo' in model.variables
                 assert 'tiny_EX_glc__D_e' in model.reactions
-                assert model.objective.expression == model.reactions.get_by_id(
-                    'Biomass_Ecoli_core').flux_expression
+                assert (model.objective.expression.simplify() ==
+                        model.reactions.get_by_id(
+                            'Biomass_Ecoli_core').flux_expression.simplify())
             assert 'ex1' not in model.reactions
             assert 'constraint' not in model.constraints
             assert 'foo' not in model.variables
@@ -853,7 +854,7 @@ def test_change_objective(self, model):
 
     def test_problem_properties(self, model):
         new_variable = model.problem.Variable("test_variable")
-        new_constraint = model.problem.Constraint(S.Zero,
+        new_constraint = model.problem.Constraint(Zero,
                                                   name="test_constraint")
         model.add_cons_vars([new_variable, new_constraint])
         assert "test_variable" in model.variables.keys()