cobra pie 0.2.0b1

cobra/__init__.py

@@ -4,11 +4,20 @@
 if __name == 'java':
     #raise Exception("Experimental modules of numpy/scipy for java that are" +\
     #"not yet ready for prime time.")
-    import oven.danielhyduke.jython.numpy as numpy
-    import oven.danielhyduke.jython.scipy as scipy
-    from core import Object, Formula, Metabolite, Gene, Reaction, Model
+    #import oven.danielhyduke.jython.numpy as numpy
+    #import oven.danielhyduke.jython.scipy as scipy
+    from warnings import warn
+    warn("COBRA for Python is not optimized for JAVA.  If it's slow or crashes consider increasing JVM memory")
+    from core import Object, Formula, Metabolite, Gene, Reaction, Model, DictList
 else:
-    from core import *
+    #from core import *
+    from core import Object, Formula, Metabolite, Gene, Reaction, Model, DictList
+    try:
+        from core import ArrayBasedModel
+    except Exception, e:
+        from warnings import warn
+        warn("cobra.ArrayBasedModel class is unavailable: %s"%repr(e))
+
 del __name
 __doc__ = """
 """

cobra/core/ArrayBasedModel.py

@@ -114,6 +114,7 @@ def add_metabolites(self, metabolite_list,
 
         """
         Model.add_metabolites(metabolite_list)
+        self.constraint_sense = [x._constraint_sense for x in self.metabolites]
         if self.S is not None and expand_stoichiometric_matrix:
             s_expansion = len(self.metabolites) - self.S.shape[0]
             if s_expansion > 0:
@@ -210,6 +211,7 @@ def _update_reaction_vectors(self):
         self.lower_bounds = array(lower_bounds)
         self.upper_bounds = array(upper_bounds)
         self.objective_coefficients = array(objective_coefficients)
+
     def _update_metabolite_vectors(self):
         """regenerates _b and _constraint_sense
 
@@ -223,10 +225,9 @@ def _update_metabolite_vectors(self):
         [(_b.append(x._bound),
           _constraint_sense.append(x._constraint_sense))
          for x in self.metabolites]
-        self._b = array(_b)
-        self._constraint_sense = _constraint_sense
+        self.b = array(_b)
+        self.constraint_sense = _constraint_sense
 
-
     def _update_matrices(self, reaction_list=None):
         """
         reaction_list: None or a list of cobra.Reaction objects that are in

cobra/core/Metabolite.py

@@ -2,6 +2,9 @@
 #######################
 #BEGIN Class Metabolite
 #
+#Is it better to restrict a Metabolite to a single model or
+#should we allow a Metabolite to be associated with multiple models?
+#
 from copy import deepcopy
 import re
 from .Formula import Formula
@@ -10,14 +13,15 @@ class Metabolite(Object):
     """Metabolite is a class for holding information regarding
     a metabolite in a cobra.Reaction object.
 
-    TODO: Clean up.  Allow for creation of empty metabolite
+        
     """
     ## __slots__ = ['id', 'formula', 'name', 'compartment',
     ##              'charge']
     ##def __setstate__(self, the_dict):
     ##    from cobra.core.Metabolite import Metabolite
     ##    Object.__setstate__(self, the_dict)
     ##    [self.__setattr__(k, v) for k, v in the_dict]
+
     def __init__(self, id=None, formula=None,
                  name=None, compartment=None):
         """
@@ -60,6 +64,7 @@ def __getstate__(self):
         state = Object.__getstate__(self)
         state['_reaction'] = set()
         return state
+
     def parse_composition(self):
         """Breaks the chemical formula down by element.
         Useful for making sure Reactions are balanced.'
@@ -80,6 +85,7 @@ def copy(self):
         """
         new_metabolite = deepcopy(self)
         return new_metabolite
+
     def guided_copy(self, the_model):
         """Trying to make a faster copy procedure for cases where large
         numbers of metabolites might be copied.  Such as when copying reactions.
@@ -91,26 +97,51 @@ def guided_copy(self, the_model):
         the_copy._model = the_model
         the_copy._reaction = set()
         return(the_copy)
-
-    def remove_from_model(self, model=None):
+
+    def get_model(self):
+        """Returns the Model object that contain this Object
+
+        """
+        return self._model
+
+    def get_reaction(self):
+        """Returns a list of Reactions that contain this Object
+
+        """
+        return list(self._reaction)
+
+    def remove_from_model(self, model=None, method='subtractive'):
         """Removes the association
 
-        model: :class:`~cobra.core.Model` object.
+        model: None or :class:`~cobra.core.Model` object.
             Remove the reaction from this model.
+
+        method: 'subtractive' or 'destructive'.  If 'subtractive' then
+        the metabolite is removed from all associated reactions.  If
+        'destructive' then all associated reactions are removed from the Model.
+
+    
         
         """
-        # why is model being taken in as a parameter?
+        # why is model being taken in as a parameter? This plays
+        #back to the question of allowing a Metabolite to be associated
+        #with multiple Models
         if model != self._model and model is not None:
             raise Exception('%s not in %s ergo it cannot be removed. (%s)'%(self,
                                                                   model,
                                                                   self._model))
 
         self._model.metabolites.remove(self)
         self._model = None
-        for the_reaction in self._reaction:
-            the_coefficient = the_reaction._metabolites[self]
-            the_reaction.subtract_metabolites({self: the_coefficient})
-
+        if method.lower() == 'subtractive':
+            for the_reaction in list(self._reaction):
+                the_coefficient = the_reaction._metabolites[self]
+                the_reaction.subtract_metabolites({self: the_coefficient})
+        elif method.lower() == 'destructive':
+            [x.remove_from_model() for x in self._reaction()]
+        else:
+            raise Exception("method for remove_from_model must be 'subtractive' " +\
+                            "or 'destructive'.  You entered: '%s'"%method)
 
 #
 #END Class Metabolite

cobra/core/Model.py

@@ -180,7 +180,9 @@ def _update_reaction(self, reaction):
             self.reactions[reaction_index] = the_reaction
 
     def update(self):
-        """
+        """Non functional.  Model.update is moved to ArrayBasedModel.  Please use
+        the to_array_based_model property to create an ArrayBasedModel.
+        
         """
         raise Exception("Model.update is moved to ArrayBasedModel.  Please use \n"
                         "the to_array_based_model property to create an ArrayBasedModel.")

cobra/core/Object.py

@@ -14,13 +14,15 @@ def __init__(self, id=None):
         #formatting notes and annotation
         self.notes = {}
         self.annotation = {}
+
     def __getstate__(self):
         """To prevent excessive replication during deepcopy.
         """
         state = self.__dict__.copy()
         if '_model' in state:
             state['_model'] = None
         return state
+
     def guided_copy(self):
         """Trying to make a faster copy procedure for cases where large
         numbers of metabolites might be copied.  Such as when copying reactions.
@@ -32,6 +34,7 @@ def guided_copy(self):
         [setattr(the_copy, k, v)
          for k, v in self.__dict__.iteritems()]
         return(the_copy)
+
     ## def __setstate__(self, state):
     ##     self.__dict__.update(state)
     ## def __getstate__(self):
@@ -45,6 +48,7 @@ def guided_copy(self):
     #Not the best idea.  This will be removed in the next major
     #release
     #
+
     def __lt__(self, other):
         if hasattr(other, 'id'):
             x = self.id < other.id
@@ -74,6 +78,7 @@ def __ge__(self, other):
         return x
 
     def __ne__(self, other):
+        x = True
         if hasattr(other, 'id'):
             x = self.id != other.id
         elif type(other) == type(self.id):

cobra/core/Reaction.py

@@ -2,6 +2,9 @@
 #######################
 #BEGIN Class Reaction
 #
+#Is it better to restrict a Reaction to a single model or
+#should we allow a Reaction to be associated with multiple models?
+#
 from collections import defaultdict
 import re
 from copy import deepcopy
@@ -62,13 +65,22 @@ def _update_awareness(self):
         [x._reaction.add(self) for x in self._metabolites]
         [x._reaction.add(self) for x in self._genes]
 
+    def get_model(self):
+        """Returns the Model object that this Reaction is associated with.
+
+        """
+        return self._model
+
+
     def remove_from_model(self, model=None):
         """Removes the association
 
         model: cobra.Model object.  remove the reaction from this model.
         
         """
-        # why is model being taken in as a parameter?
+        # why is model being taken in as a parameter? This plays
+        #back to the question of allowing a Metabolite to be associated
+        #with multiple Models
         if model != self._model and model is not None:
             raise Exception('%s not in %s ergo it cannot be removed. (%s)'%(self,
                                                                   model,
@@ -133,6 +145,7 @@ def copy(self):
         new_reaction = deepcopy(self)
         ## self._model = the_model
         return new_reaction
+
     def guided_copy(self, the_model, metabolite_dict, gene_dict=None):
         """Trying to make a faster copy procedure for cases where large
         numbers of metabolites might be copied.  Such as when copying reactions.
@@ -245,6 +258,19 @@ def parse_gene_association(self, the_type='gene'):
             #Make the gene aware that it is involved in this reaction
             [x._reaction.add(self) for x in self._genes]
 
+
+    def add_gene_reaction_rule(self, the_rule):
+        """This adds a gene reaction rule.
+
+        the_rule:  A boolean representation of the gene requirements for this reaction
+        to be active as described in Schellenberger et al 2011 Nature Protocols 6(9):1290-307.
+
+        Note that this method currently replaces any pre-existing rules
+        
+        """
+        self.gene_reaction_rule = the_rule
+        self.parse_gene_association()
+
     def get_reactants(self):
         """Return a list of reactants for the reaction.
 
@@ -259,6 +285,12 @@ def get_products(self):
         return [k for k, v in self._metabolites.items()
                 if v > 0]
 
+    def get_gene(self):
+        """Return a list of reactants for the reaction.
+
+        """
+        return self._genes.keys()
+
 
     def get_coefficient(self, the_metabolite):
         """Return the stoichiometric coefficient for a metabolite in
@@ -494,6 +526,11 @@ def print_values(self):
                          self._metabolites]:
             print repr(the_item)
 
+    def get_compartments(self):
+        """
+        """
+        return(list(set([x.compartment for x in self._metabolites])))
+
 #DEPRECATED SECTION
     def process_prefixed_reaction(self):
         """Deal with reaction names that have a prefix.

cobra/core/__init__.py

@@ -8,6 +8,9 @@
 from .Model import Model 
 from os import name as __name
 if __name != 'java':
-    from .ArrayBasedModel import ArrayBasedModel 
-
+    try:
+        from .ArrayBasedModel import ArrayBasedModel 
+    except Exception, e:
+        from warnings import warn
+        warn("ArrayBasedModel is not accessible: %s"%e)
 

cobra/examples/01_create_model.py

@@ -41,10 +41,13 @@
 #Print the reaction to make sure it worked:
 print reaction.reaction
 
-reaction.gene_reaction_rule = '( STM2378  or  STM1197 )'
+#A boolean representation of the gene requirements for this reaction
+#to be active as described in Schellenberger et al 2011 Nature Protocols 6(9):1290-307.
+gene_reaction_rule = '( STM2378  or  STM1197 )'
+
 #The next step will create cobra.Gene objects from the gene reaction rule, which will
 #be used later by the cobra.Model to modulate reaction bounds after deleting genes.
-reaction.parse_gene_association()
+reaction.add_gene_reaction_rule(gene_reaction_rule)
 
 #The model's initially empty:
 print '%i reactions in model'%len(cobra_model.reactions)
@@ -66,5 +69,5 @@
     print '%s has formula %s'%(x,x.formula)
 for x in cobra_model.genes:
     print 'gene %s is associated with reactions:'%x
-    for y in x._reaction:
+    for y in x.get_reaction():
         print '\t%s'%y

cobra/examples/07_ice_cream_milp.py

@@ -2,21 +2,8 @@
 #
 # Advanced user example showing how to set up and solve an MILP
 #
-
 from cobra import Model, Metabolite, Reaction
-try:
-    import glpk
-    solver = 'glpk'
-except:
-    try:
-        from cplex import Cplex
-        solver = 'cplex'
-    except:
-        try:
-            from gurobipy import Model as gurobi_model
-            solver = 'gurobi'
-        except:
-            raise Exception("Couldn't import glpk, cplex, or gurobi")
+solver = 'cplex'  #With libglpk-java there is an untraced memory bug.
 
 cone_selling_price = 7.
 cone_production_cost = 3.
@@ -110,22 +97,18 @@
 Cone_production.objective_coefficient = -1*cone_production_cost
 Popsicle_production.objective_coefficient = -1*popsicle_production_cost
 
-reaction_dict = dict([(x.id, x) for x in cobra_model.reactions])
 
 production_capacity_constraint = Metabolite(id='production_capacity_constraint')
 production_capacity_constraint._constraint_sense = 'L'
 production_capacity_constraint._bound = starting_budget;
 
-the_reaction = reaction_dict['Cone_production']
-the_reaction.add_metabolites({production_capacity_constraint: cone_production_cost })
-
-the_reaction = reaction_dict['Popsicle_production']
-the_reaction.add_metabolites({production_capacity_constraint: popsicle_production_cost })
+Cone_production.add_metabolites({production_capacity_constraint: cone_production_cost })
+Popsicle_production.add_metabolites({production_capacity_constraint: popsicle_production_cost })
 
 print
 print('Here is what happens in the continuous (LP) case...')
 
-the_program = cobra_model.optimize(objective_sense='maximize',solver=solver)
+the_program = cobra_model.optimize(objective_sense='maximize', solver=solver)
 print
 print('Status is: %s'%cobra_model.solution.status)
 print('Objective value is: %1.2f'%cobra_model.solution.f)
@@ -142,9 +125,7 @@
 print
 print('Who wants 1/3 of a cone, WTF???  Cones and popsicles are units aka integers, reformulate as MILP')
 Cone_production.variable_kind = 'integer'
-Cone_production.upper_bound = 50000;
 Popsicle_production.variable_kind = 'integer'
-Popsicle_production.upper_bound = 50000;
 
 the_program = cobra_model.optimize(objective_sense='maximize', solver=solver)
 print