.travis.yml

@@ -29,7 +29,7 @@ install:
   - pip install --upgrade pip setuptools wheel tox tox-travis
 
 script:
-  - tox -- --benchmark-skip --cov-report xml --cov-report term
+  - travis_wait tox -- --benchmark-skip --cov-report xml --cov-report term
   - bash <(curl -s https://codecov.io/bash)
 
 stages:
@@ -41,30 +41,31 @@ stages:
 #before_deploy:
 #  - source scripts/prepare_notes.sh
 
+jobs:
+  include:
+    - stage: deploy
+      python: "3.6"
+      install:
+        - pip install --upgrade pip setuptools wheel
+      script:
+        - echo "Deploying..."
       deploy:
         - provider: pypi
           skip_cleanup: true
           distributions: sdist bdist_wheel
-    install: skip
-    script:
-      - echo "Deploying COBRApy to PyPI."
+          user: cobrapy
+          password: $PYPI_PASSWORD
           on:
             tags: true
             repo: $GITHUB_REPO
-    user: cobrapy
-    password: $PYPI_PASSWORD
         - provider: script
           skip_cleanup: true
-    install: skip
           script: scripts/deploy_website.sh
           on:
             tags: true
             repo: $GITHUB_REPO
         - provider: releases
           skip_cleanup: true
-    install: skip
-    script:
-    - echo "Releasing COBRApy to GitHub."
           api_key: $GITHUB_TOKEN
           body: "Please see https://github.com/opencobra/cobrapy/tree/${TRAVIS_TAG}/release-notes/${TRAVIS_TAG}.md for the full release notes."
           on:

appveyor.yml

@@ -4,57 +4,32 @@ branches:
  - devel
  - /^[0-9]+\.[0-9]+\.[0-9]+[.0-9ab]*$/
 
-environment:
-
-  global:
-    PIP_CACHE_DIR: "pip_cache"
-    PIP_DISABLE_PIP_VERSION_CHECK: "yes"
+clone_depth: 2
 
+environment:
   matrix:
-    - PYTHON: "C:\\Miniconda-x64"
-      CONDA: true
-
     - PYTHON: "C:\\Miniconda36-x64"
-      CONDA: true
-
-    - PYTHON: "C:\\Python27-x64"
-      TOXENV: py27
-
+      TOXENV: py36
+    - PYTHON: "C:\\Miniconda37-x64"
+      TOXENV: py37
     - PYTHON: "C:\\Python36-x64"
       TOXENV: py36
 
 matrix:
   fast_finish: true
 
-clone_depth: 2
-
-init:
-  - "ECHO %PYTHON%"
-  - "set PATH=%PYTHON%;%PYTHON%\\Scripts;%PATH%"
-  - python -V
-
 cache:
-    - pip_cache -> appveyor.yml
-
+  - '%LOCALAPPDATA%\pip\Cache'
 
-install:
-  - ps: |
-      if ($env:CONDA -eq "true") {
-        conda config --set always_yes yes --set changeps1 no;
-        conda config --add channels conda-forge
-        conda install -q pip setuptools wheel numpy scipy;
-        pip install -r develop-requirements.txt;
-        pip install .
-      } else {
-        pip install --upgrade setuptools wheel tox
-      }
+platform: x64
 
 build: off
 
+install:
+  - set PATH=%PYTHON%;%PYTHON%\Scripts;%PYTHON%\Library\bin;%PATH%"
+  - "ECHO %PYTHON%"
+  - python --version
+  - python -m pip install --upgrade --disable-pip-version-check pip setuptools wheel tox
+
 test_script:
-  - ps: |
-      if ($env:CONDA -eq "true") {
-        pytest --benchmark-skip cobra/test
-      } else {
-        tox
-      }
+  - tox

cobra/__init__.py

@@ -13,9 +13,10 @@ from cobra.core import (
     Species)
 from cobra import flux_analysis
 from cobra import io
+from cobra import sampling
 from cobra.util import show_versions
 
-__version__ = "0.14.2"
+__version__ = "0.15.4"
 
 # set the warning format to be prettier and fit on one line
 _cobra_path = _dirname(_abspath(__file__))

cobra/core/__init__.py

@@ -9,5 +9,6 @@ from cobra.core.metabolite import Metabolite
 from cobra.core.model import Model
 from cobra.core.object import Object
 from cobra.core.reaction import Reaction
-from cobra.core.solution import LegacySolution, Solution, get_solution
+from cobra.core.group import Group
+from cobra.core.solution import Solution, LegacySolution, get_solution
 from cobra.core.species import Species

cobra/core/configuration.py

@@ -34,6 +34,8 @@ class BaseConfiguration(object):
     solver : {"glpk", "cplex", "gurobi"}
         The default solver for new models. The solver choices are the ones
         provided by `optlang` and solvers installed in your environment.
+    tolerance: float
+        The default tolerance for the solver being used (default 1E-09).
     lower_bound : float
         The standard lower bound for reversible reactions (default -1000).
     upper_bound : float
@@ -50,8 +52,10 @@ class BaseConfiguration(object):
 
     def __init__(self):
         self._solver = None
+        self.tolerance = 1E-07
         self.lower_bound = None
         self.upper_bound = None
+
         # Set the default solver from a preferred order.
         for name in ["gurobi", "cplex", "glpk"]:
             try:
@@ -60,7 +64,9 @@ class BaseConfiguration(object):
                 continue
             else:
                 break
+
         self.bounds = -1000.0, 1000.0
+
         try:
             self.processes = cpu_count()
         except NotImplementedError:
@@ -99,10 +105,49 @@ class BaseConfiguration(object):
         self.upper_bound = bounds[1]
 
     def __repr__(self):
-        return "solver: {}\nlower_bound: {}\nupper_bound: {}\n" \
-            "processes: {}".format(
-                interface_to_str(self.solver), self.lower_bound,
-                self.upper_bound, self.processes)
+        return """
+        solver: {solver}
+        solver tolerance: {tolerance}
+        lower_bound: {lower_bound}
+        upper_bound: {upper_bound}
+        processes: {processes}""".format(
+            solver=interface_to_str(self.solver),
+            tolerance=self.tolerance,
+            lower_bound=self.lower_bound,
+            upper_bound=self.upper_bound,
+            processes=self.processes,
+        )
+
+    def _repr_html_(self):
+        return """
+        <table>
+            <tr>
+                <td><strong>Solver</strong></td>
+                <td>{solver}</td>
+            </tr>
+            <tr>
+                <td><strong>Solver tolerance</strong></td>
+                <td>{tolerance}</td>
+            </tr>
+            <tr>
+                <td><strong>Lower bound</strong></td>
+                <td>{lower_bound}</td>
+            </tr>
+            <tr>
+                <td><strong>Upper bound</strong></td>
+                <td>{upper_bound}</td>
+            </tr>
+            <tr>
+                <td><strong>Processes</strong></td>
+                <td>{processes}</td>
+            </tr>
+        </table>""".format(
+            solver=interface_to_str(self.solver),
+            tolerance=self.tolerance,
+            lower_bound=self.lower_bound,
+            upper_bound=self.upper_bound,
+            processes=self.processes,
+        )
 
 
 class Configuration(with_metaclass(Singleton, BaseConfiguration)):

cobra/core/gene.py

@@ -249,6 +249,11 @@ class Gene(Species):
         the_gene_re = re.compile('(^|(?<=( |\()))%s(?=( |\)|$))' %
                                  re.escape(self.id))
 
+        # remove reference to the gene in all groups
+        associated_groups = self._model.get_associated_groups(self)
+        for group in associated_groups:
+            group.remove_members(self)
+
         self._model.genes.remove(self)
         self._model = None
 

cobra/core/group.py

+# -*- coding: utf-8 -*-
+
+"""Define the group class."""
+
+from __future__ import absolute_import
+
+from warnings import warn
+
+from six import string_types
+
+from cobra.core.object import Object
+
+
+class Group(Object):
+    """
+    Manage groups via this implementation of the SBML group specification.
+
+    `Group` is a class for holding information regarding a pathways,
+    subsystems, or other custom groupings of objects within a cobra.Model
+    object.
+
+    Parameters
+    ----------
+    id : str
+        The identifier to associate with this group
+    name : str, optional
+        A human readable name for the group
+    members : iterable, optional
+        A list object containing references to cobra.Model-associated objects
+        that belong to the group.
+    kind : {"collection", "classification", "partonomy"}, optional
+        The kind of group, as specified for the Groups feature in the SBML
+        level 3 package specification. Can be any of "classification",
+        "partonomy", or "collection". The default is "collection".
+        Please consult the SBML level 3 package specification to ensure you
+        are using the proper value for kind. In short, members of a
+        "classification" group should have an "is-a" relationship to the group
+        (e.g. member is-a polar compound, or member is-a transporter).
+        Members of a "partonomy" group should have a "part-of" relationship
+        (e.g. member is part-of glycolysis). Members of a "collection" group
+        do not have an implied relationship between the members, so use this
+        value for kind when in doubt (e.g. member is a gap-filled reaction,
+        or member is involved in a disease phenotype).
+    """
+    KIND_TYPES = ("collection", "classification", "partonomy")
+
+    def __init__(self, id, name='', members=None, kind=None):
+        Object.__init__(self, id, name)
+
+        self._members = set() if members is None else set(members)
+        self._kind = None
+        self.kind = "collection" if kind is None else kind
+        # self.model is None or refers to the cobra.Model that
+        # contains self
+        self._model = None
+
+    def __len__(self):
+        return len(self._members)
+
+    # read-only
+    @property
+    def members(self):
+        return self._members
+
+    @property
+    def kind(self):
+        return self._kind
+
+    @kind.setter
+    def kind(self, kind):
+        kind = kind.lower()
+        if kind in self.KIND_TYPES:
+            self._kind = kind
+        else:
+            raise ValueError(
+                "Kind can only by one of: {}.".format(", ".join(
+                    self.KIND_TYPES)))
+
+    def add_members(self, new_members):
+        """
+        Add objects to the group.
+
+        Parameters
+        ----------
+        new_members : list
+            A list of cobrapy objects to add to the group.
+
+        """
+
+        if isinstance(new_members, string_types) or \
+                hasattr(new_members, "id"):
+            warn("need to pass in a list")
+            new_members = [new_members]
+
+        self._members.update(new_members)
+
+    def remove_members(self, to_remove):
+        """
+        Remove objects from the group.
+
+        Parameters
+        ----------
+        to_remove : list
+            A list of cobra objects to remove from the group
+        """
+
+        if isinstance(to_remove, string_types) or \
+                hasattr(to_remove, "id"):
+            warn("need to pass in a list")
+            to_remove = [to_remove]
+
+        self._members.difference_update(to_remove)

cobra/core/metabolite.py

@@ -79,7 +79,7 @@ class Metabolite(Species):
         # necessary for some old pickles which use the deprecated
         # Formula class
         tmp_formula = str(self.formula)
-        # commonly occuring characters in incorrectly constructed formulas
+        # commonly occurring characters in incorrectly constructed formulas
         if "*" in tmp_formula:
             warn("invalid character '*' found in formula '%s'" % self.formula)
             tmp_formula = tmp_formula.replace("*", "")

cobra/core/model.py

@@ -15,6 +15,9 @@ from six import iteritems, string_types
 
 from cobra.core.configuration import Configuration
 from cobra.core.dictlist import DictList
+from cobra.core.gene import Gene
+from cobra.core.group import Group
+from cobra.core.metabolite import Metabolite
 from cobra.core.object import Object
 from cobra.core.reaction import Reaction
 from cobra.core.solution import get_solution
@@ -55,12 +58,17 @@ class Model(Object):
     genes : DictList
         A DictList where the key is the gene identifier and the value a
         Gene
+    groups : DictList
+        A DictList where the key is the group identifier and the value a
+        Group
     solution : Solution
         The last obtained solution from optimizing the model.
+
     """
 
     def __setstate__(self, state):
-        """Make sure all cobra.Objects in the model point to the model."""
+        """Make sure all cobra.Objects in the model point to the model.
+        """
         self.__dict__.update(state)
         for y in ['reactions', 'genes', 'metabolites']:
             for x in getattr(self, y):
@@ -93,6 +101,7 @@ class Model(Object):
             self.genes = DictList()
             self.reactions = DictList()  # A list of cobra.Reactions
             self.metabolites = DictList()  # A list of cobra.Metabolites
+            self.groups = DictList()  # A list of cobra.Groups
             # genes based on their ids {Gene.id: Gene}
             self._compartments = {}
             self._contexts = []
@@ -101,11 +110,15 @@ class Model(Object):
 
             # if not hasattr(self, '_solver'):  # backwards compatibility
             # with older cobrapy pickles?
+
             interface = CONFIGURATION.solver
             self._solver = interface.Model()
             self._solver.objective = interface.Objective(Zero)
             self._populate_solver(self.reactions, self.metabolites)
 
+            self._tolerance = None
+            self.tolerance = CONFIGURATION.tolerance
+
     @property
     def solver(self):
         """Get or set the attached solver instance.
@@ -149,6 +162,34 @@ class Model(Object):
             return
         self._solver = interface.Model.clone(self._solver)
 
+    @property
+    def tolerance(self):
+        return self._tolerance
+
+    @tolerance.setter
+    def tolerance(self, value):
+        solver_tolerances = self._solver.configuration.tolerances
+
+        try:
+            solver_tolerances.feasibility = value
+        except AttributeError:
+            LOGGER.info("The current solver doesn't allow setting"
+                        "feasibility tolerance.")
+
+        try:
+            solver_tolerances.optimality = value
+        except AttributeError:
+            LOGGER.info("The current solver doesn't allow setting"
+                        "optimality tolerance.")
+
+        try:
+            solver_tolerances.integrality = value
+        except AttributeError:
+            LOGGER.info("The current solver doesn't allow setting"
+                        "integrality tolerance.")
+
+        self._tolerance = value
+
     @property
     def description(self):
         warn("description deprecated", DeprecationWarning)
@@ -281,7 +322,7 @@ class Model(Object):
         """
         new = self.__class__()
         do_not_copy_by_ref = {"metabolites", "reactions", "genes", "notes",
-                              "annotation"}
+                              "annotation", "groups"}
         for attr in self.__dict__:
             if attr not in do_not_copy_by_ref:
                 new.__dict__[attr] = self.__dict__[attr]
@@ -327,6 +368,39 @@ class Model(Object):
                 new_gene = new.genes.get_by_id(gene.id)
                 new_reaction._genes.add(new_gene)
                 new_gene._reaction.add(new_reaction)
+
+        new.groups = DictList()
+        do_not_copy_by_ref = {"_model", "_members"}
+        # Groups can be members of other groups. We initialize them first and
+        # then update their members.
+        for group in self.groups:
+            new_group = group.__class__(group.id)
+            for attr, value in iteritems(group.__dict__):
+                if attr not in do_not_copy_by_ref:
+                    new_group.__dict__[attr] = copy(value)
+            new_group._model = new
+            new.groups.append(new_group)
+        for group in self.groups:
+            new_group = new.groups.get_by_id(group.id)
+            # update awareness, as in the reaction copies
+            new_objects = []
+            for member in group.members:
+                if isinstance(member, Metabolite):
+                    new_object = new.metabolites.get_by_id(member.id)
+                elif isinstance(member, Reaction):
+                    new_object = new.reactions.get_by_id(member.id)
+                elif isinstance(member, Gene):
+                    new_object = new.genes.get_by_id(member.id)
+                elif isinstance(member, Group):
+                    new_object = new.genes.get_by_id(member.id)
+                else:
+                    raise TypeError(
+                        "The group member {!r} is unexpectedly not a "
+                        "metabolite, reaction, gene, nor another "
+                        "group.".format(member))
+                new_objects.append(new_object)
+            new_group.add_members(new_objects)
+
         try:
             new._solver = deepcopy(self.solver)
             # Cplex has an issue with deep copies
@@ -409,6 +483,11 @@ class Model(Object):
         for x in metabolite_list:
             x._model = None
 
+            # remove reference to the metabolite in all groups
+            associated_groups = self.get_associated_groups(x)
+            for group in associated_groups:
+                group.remove_members(x)
+
             if not destructive:
                 for the_reaction in list(x._reaction):
                     the_coefficient = the_reaction._metabolites[x]
@@ -505,6 +584,7 @@ class Model(Object):
         (0, 1000.0)
         >>> demand.build_reaction_string()
         'atp_c --> '
+
         """
         ub = CONFIGURATION.upper_bound if ub is None else ub
         lb = CONFIGURATION.lower_bound if lb is None else lb
@@ -684,6 +764,100 @@ class Model(Object):
                             if context:
                                 context(partial(self.genes.add, gene))
 
+                # remove reference to the reaction in all groups
+                associated_groups = self.get_associated_groups(reaction)
+                for group in associated_groups:
+                    group.remove_members(reaction)
+
+    def add_groups(self, group_list):
+        """Add groups to the model.
+
+        Groups with identifiers identical to a group already in the model are
+        ignored.
+
+        If any group contains members that are not in the model, these members
+        are added to the model as well. Only metabolites, reactions, and genes
+        can have groups.
+
+        Parameters
+        ----------
+        group_list : list
+            A list of `cobra.Group` objects to add to the model.
+        """
+
+        def existing_filter(group):
+            if group.id in self.groups:
+                LOGGER.warning(
+                    "Ignoring group '%s' since it already exists.", group.id)
+                return False
+            return True
+
+        if isinstance(group_list, string_types) or \
+                hasattr(group_list, "id"):
+            warn("need to pass in a list")
+            group_list = [group_list]
+
+        pruned = DictList(filter(existing_filter, group_list))
+
+        for group in pruned:
+            group._model = self
+            for member in group.members:
+                # If the member is not associated with the model, add it
+                if isinstance(member, Metabolite):
+                    if member not in self.metabolites:
+                        self.add_metabolites([member])
+                if isinstance(member, Reaction):
+                    if member not in self.reactions:
+                        self.add_reactions([member])
+                # TODO(midnighter): `add_genes` method does not exist.
+                # if isinstance(member, Gene):
+                #     if member not in self.genes:
+                #         self.add_genes([member])
+
+            self.groups += [group]
+
+    def remove_groups(self, group_list):
+        """Remove groups from the model.
+
+        Members of each group are not removed
+        from the model (i.e. metabolites, reactions, and genes in the group
+        stay in the model after any groups containing them are removed).
+
+        Parameters
+        ----------
+        group_list : list
+            A list of `cobra.Group` objects to remove from the model.
+        """
+
+        if isinstance(group_list, string_types) or \
+                hasattr(group_list, "id"):
+            warn("need to pass in a list")
+            group_list = [group_list]
+
+        for group in group_list:
+            # make sure the group is in the model
+            if group.id not in self.groups:
+                LOGGER.warning("%r not in %r. Ignored.", group, self)
+            else:
+                self.groups.remove(group)
+                group._model = None
+
+    def get_associated_groups(self, element):
+        """Returns a list of groups that an element (reaction, metabolite, gene)
+        is associated with.
+
+        Parameters
+        ----------
+        element: `cobra.Reaction`, `cobra.Metabolite`, or `cobra.Gene`
+
+        Returns
+        -------
+        list of `cobra.Group`
+            All groups that the provided object is a member of
+        """
+        # check whether the element is associated with the model
+        return [g for g in self.groups if element in g.members]
+
     def add_cons_vars(self, what, **kwargs):
         """Add constraints and variables to the model's mathematical problem.
 
@@ -922,6 +1096,7 @@ class Model(Object):
             self.reactions._generate_index()
             self.metabolites._generate_index()
             self.genes._generate_index()
+            self.groups._generate_index()
         if rebuild_relationships:
             for met in self.metabolites:
                 met._reaction.clear()
@@ -932,8 +1107,9 @@ class Model(Object):
                     met._reaction.add(rxn)
                 for gene in rxn._genes:
                     gene._reaction.add(rxn)
+
         # point _model to self
-        for l in (self.reactions, self.genes, self.metabolites):
+        for l in (self.reactions, self.genes, self.metabolites, self.groups):
             for e in l:
                 e._model = self
 
@@ -1108,6 +1284,9 @@ class Model(Object):
             </tr><tr>
                 <td><strong>Number of reactions</strong></td>
                 <td>{num_reactions}</td>
+            </tr><tr>
+                <td><strong>Number of groups</strong></td>
+                <td>{num_groups}</td>
             </tr><tr>
                 <td><strong>Objective expression</strong></td>
                 <td>{objective}</td>
@@ -1120,6 +1299,7 @@ class Model(Object):
             address='0x0%x' % id(self),
             num_metabolites=len(self.metabolites),
             num_reactions=len(self.reactions),
+            num_groups=len(self.groups),
             objective=format_long_string(str(self.objective.expression), 100),
             compartments=", ".join(
                 v if v else k for k, v in iteritems(self.compartments)

cobra/core/reaction.py

@@ -25,7 +25,7 @@ from cobra.util.solver import (
 from cobra.util.util import format_long_string
 
 
-CONFIGURATION = Configuration()
+config = Configuration()
 
 # precompiled regular expressions
 # Matches and/or in a gene reaction rule
@@ -61,7 +61,7 @@ class Reaction(Object):
     """
 
     def __init__(self, id=None, name='', subsystem='', lower_bound=0.0,
-                 upper_bound=None, objective_coefficient=0.0):
+                 upper_bound=None):
         Object.__init__(self, id, name)
         self._gene_reaction_rule = ''
         self.subsystem = subsystem
@@ -80,17 +80,11 @@ class Reaction(Object):
         # contains self
         self._model = None
 
-        if objective_coefficient != 0:
-            raise NotImplementedError('setting objective coefficient when '
-                                      'creating reaction is no longer '
-                                      'supported. Use the model.objective '
-                                      'setter')
-
         # from cameo ...
         self._lower_bound = lower_bound if lower_bound is not None else \
-            CONFIGURATION.lower_bound
+            config.lower_bound
         self._upper_bound = upper_bound if upper_bound is not None else \
-            CONFIGURATION.upper_bound
+            config.upper_bound
 
     def _set_id_with_model(self, value):
         if value in self.model.reactions:
@@ -1050,14 +1044,12 @@ class Reaction(Object):
         # reversible case
         arrow_match = reversible_arrow_finder.search(reaction_str)
         if arrow_match is not None:
-            self.lower_bound = -1000
-            self.upper_bound = 1000
+            self.bounds = config.lower_bound, config.upper_bound
         else:  # irreversible
             # try forward
             arrow_match = forward_arrow_finder.search(reaction_str)
             if arrow_match is not None:
-                self.upper_bound = 1000
-                self.lower_bound = 0
+                self.bounds = 0, config.upper_bound
             else:
                 # must be reverse
                 arrow_match = reverse_arrow_finder.search(reaction_str)
@@ -1065,8 +1057,7 @@ class Reaction(Object):
                     raise ValueError("no suitable arrow found in '%s'" %
                                      reaction_str)
                 else:
-                    self.upper_bound = 0
-                    self.lower_bound = -1000
+                    self.bounds = config.lower_bound, 0
         reactant_str = reaction_str[:arrow_match.start()].strip()
         product_str = reaction_str[arrow_match.end():].strip()
 

cobra/flux_analysis/__init__.py

@@ -3,6 +3,7 @@
 from cobra.flux_analysis.deletion import (
     double_gene_deletion, double_reaction_deletion, single_gene_deletion,
     single_reaction_deletion)
+from cobra.flux_analysis.fastcc import fastcc
 from cobra.flux_analysis.gapfilling import gapfill
 from cobra.flux_analysis.geometric import geometric_fba
 from cobra.flux_analysis.loopless import (loopless_solution, add_loopless)
@@ -13,4 +14,3 @@ from cobra.flux_analysis.variability import (
     flux_variability_analysis)
 from cobra.flux_analysis.phenotype_phase_plane import production_envelope
 from cobra.flux_analysis.room import add_room, room
-from cobra.flux_analysis.sampling import sample

cobra/flux_analysis/fastcc.py

+# -*- coding: utf-8 -*-
+
+"""Provide an implementation of FASTCC."""
+
+from __future__ import absolute_import
+
+import logging
+
+from optlang.symbolics import Zero
+
+from cobra.flux_analysis.helpers import normalize_cutoff
+
+
+LOGGER = logging.getLogger(__name__)
+
+
+def fastcc(model, flux_threshold=1.0, zero_cutoff=None):
+    r"""
+    Check consistency of a metabolic network using FASTCC [1]_.
+
+    FASTCC (Fast Consistency Check) is an algorithm for rapid and
+    efficient consistency check in metabolic networks. FASTCC is
+    a pure LP implementation and is low on computation resource
+    demand. FASTCC also circumvents the problem associated with
+    reversible reactions for the purpose. Given a global model,
+    it will generate a consistent global model i.e., remove
+    blocked reactions. For more details on FASTCC, please
+    check [1]_.
+
+    Parameters
+    ----------
+    model: cobra.Model
+        The constraint-based model to operate on.
+    flux_threshold: float, optional (default 1.0)
+        The flux threshold to consider.
+    zero_cutoff: float, optional
+        The cutoff to consider for zero flux (default model.tolerance).
+
+    Returns
+    -------
+    cobra.Model
+        The consistent constraint-based model.
+
+    Notes
+    -----
+    The LP used for FASTCC is like so:
+    maximize: \sum_{i \in J} z_i
+    s.t.    : z_i \in [0, \varepsilon] \forall i \in J, z_i \in \mathbb{R}_+
+              v_i \ge z_i \forall i \in J
+              Sv = 0 v \in B
+
+    References
+    ----------
+    .. [1] Vlassis N, Pacheco MP, Sauter T (2014)
+           Fast Reconstruction of Compact Context-Specific Metabolic Network
+           Models.
+           PLoS Comput Biol 10(1): e1003424. doi:10.1371/journal.pcbi.1003424
+
+    """
+    zero_cutoff = normalize_cutoff(model, zero_cutoff)
+
+    with model:
+        obj_vars = []
+        vars_and_cons = []
+        prob = model.problem
+
+        for rxn in model.reactions:
+            var = prob.Variable("auxiliary_{}".format(rxn.id),
+                                lb=0.0, ub=flux_threshold)
+            const = prob.Constraint(rxn.forward_variable +
+                                    rxn.reverse_variable -
+                                    var, name="constraint_{}".format(rxn.id),
+                                    lb=0.0)
+            vars_and_cons.extend([var, const])
+            obj_vars.append(var)
+        model.add_cons_vars(vars_and_cons)
+
+        model.objective = prob.Objective(Zero, sloppy=True, direction="max")
+        model.objective.set_linear_coefficients({v: 1.0 for v in obj_vars})
+
+        sol = model.optimize()
+
+    rxns_to_remove = sol.fluxes[sol.fluxes.abs() < zero_cutoff].index
+    consistent_model = model.copy()
+    consistent_model.remove_reactions(rxns_to_remove, remove_orphans=True)
+    return consistent_model

cobra/flux_analysis/geometric.py

@@ -15,7 +15,7 @@ from cobra.flux_analysis.variability import flux_variability_analysis
 LOGGER = logging.getLogger(__name__)
 
 
-def geometric_fba(model, epsilon=1E-06, max_tries=200):
+def geometric_fba(model, epsilon=1E-06, max_tries=200, processes=None):
     """
     Perform geometric FBA to obtain a unique, centered flux distribution.
 
@@ -32,6 +32,9 @@ def geometric_fba(model, epsilon=1E-06, max_tries=200):
         The convergence tolerance of the model (default 1E-06).
     max_tries: int, optional
         Maximum number of iterations (default 200).
+    processes : int, optional
+        The number of parallel processes to run. If not explicitly passed,
+        will be set from the global configuration singleton.
 
     Returns
     -------
@@ -58,7 +61,7 @@ def geometric_fba(model, epsilon=1E-06, max_tries=200):
         prob = model.problem
         add_pfba(model)  # Minimize the solution space to a convex hull.
         model.optimize()
-        fva_sol = flux_variability_analysis(model)
+        fva_sol = flux_variability_analysis(model, processes=processes)
         mean_flux = (fva_sol["maximum"] + fva_sol["minimum"]).abs() / 2
 
         # Set the gFBA constraints.
@@ -84,7 +87,7 @@ def geometric_fba(model, epsilon=1E-06, max_tries=200):
         model.objective.set_linear_coefficients({v: 1.0 for v in obj_vars})
         # Update loop variables.
         sol = model.optimize()
-        fva_sol = flux_variability_analysis(model)
+        fva_sol = flux_variability_analysis(model, processes=processes)
         mean_flux = (fva_sol["maximum"] + fva_sol["minimum"]).abs() / 2
         delta = (fva_sol["maximum"] - fva_sol["minimum"]).max()
         count = 1
@@ -92,14 +95,14 @@ def geometric_fba(model, epsilon=1E-06, max_tries=200):
                      count, delta, sol.status)
 
         # Following iterations that minimize the distance below threshold.
-        while delta > epsilon and count <= max_tries:
+        while delta > epsilon and count < max_tries:
             for rxn_id, var, u_c, l_c in updating_vars_cons:
                 var.ub = mean_flux[rxn_id]
                 u_c.ub = mean_flux[rxn_id]
                 l_c.lb = fva_sol.at[rxn_id, "minimum"]
             # Update loop variables.
             sol = model.optimize()
-            fva_sol = flux_variability_analysis(model)
+            fva_sol = flux_variability_analysis(model, processes=processes)
             mean_flux = (fva_sol["maximum"] + fva_sol["minimum"]).abs() / 2
             delta = (fva_sol["maximum"] - fva_sol["minimum"]).max()
             count += 1

cobra/flux_analysis/helpers.py

+# -*- coding: utf-8 -*-
+
+"""Helper functions for all flux analysis methods."""
+
+from __future__ import absolute_import
+
+import logging
+
+
+LOGGER = logging.getLogger(__name__)
+
+
+def normalize_cutoff(model, zero_cutoff=None):
+    """Return a valid zero cutoff value."""
+    if zero_cutoff is None:
+        return model.tolerance
+    else:
+        if zero_cutoff < model.tolerance:
+            raise ValueError(
+                "The chosen zero cutoff cannot be less than the model's "
+                "tolerance value."
+            )
+        else:
+            return zero_cutoff

cobra/flux_analysis/loopless.py

@@ -4,14 +4,20 @@
 
 from __future__ import absolute_import
 
+import logging
+
 import numpy
 from optlang.symbolics import Zero
 
 from cobra.core import get_solution
+from cobra.flux_analysis.helpers import normalize_cutoff
 from cobra.util import create_stoichiometric_matrix, nullspace
 
 
-def add_loopless(model, zero_cutoff=1e-12):
+LOGGER = logging.getLogger(__name__)
+
+
+def add_loopless(model, zero_cutoff=None):
     """Modify a model so all feasible flux distributions are loopless.
 
     In most cases you probably want to use the much faster `loopless_solution`.
@@ -28,7 +34,7 @@ def add_loopless(model, zero_cutoff=1e-12):
         The model to which to add the constraints.
     zero_cutoff : positive float, optional
         Cutoff used for null space. Coefficients with an absolute value smaller
-        than `zero_cutoff` are considered to be zero.
+        than `zero_cutoff` are considered to be zero (default model.tolerance).
 
     Returns
     -------
@@ -41,6 +47,8 @@ def add_loopless(model, zero_cutoff=1e-12):
        2011 Feb 2;100(3):544-53. doi: 10.1016/j.bpj.2010.12.3707. Erratum
        in: Biophys J. 2011 Mar 2;100(5):1381.
     """
+    zero_cutoff = normalize_cutoff(model, zero_cutoff)
+
     internal = [i for i, r in enumerate(model.reactions) if not r.boundary]
     s_int = create_stoichiometric_matrix(model)[:, numpy.array(internal)]
     n_int = nullspace(s_int).T
@@ -163,7 +171,7 @@ def loopless_solution(model, fluxes=None):
     return solution
 
 
-def loopless_fva_iter(model, reaction, solution=False, zero_cutoff=1e-6):
+def loopless_fva_iter(model, reaction, solution=False, zero_cutoff=None):
     """Plugin to get a loopless FVA solution from single FVA iteration.
 
     Assumes the following about `model` and `reaction`:
@@ -184,7 +192,7 @@ def loopless_fva_iter(model, reaction, solution=False, zero_cutoff=1e-6):
         `reaction`.
     zero_cutoff : positive float, optional
         Cutoff used for loop removal. Fluxes with an absolute value smaller
-        than `zero_cutoff` are considered to be zero.
+        than `zero_cutoff` are considered to be zero (default model.tolerance).
 
     Returns
     -------
@@ -193,6 +201,8 @@ def loopless_fva_iter(model, reaction, solution=False, zero_cutoff=1e-6):
         all_fluxes == False (default). Otherwise returns a loopless flux
         solution containing the minimum/maximum flux for `reaction`.
     """
+    zero_cutoff = normalize_cutoff(model, zero_cutoff)
+
     current = model.objective.value
     sol = get_solution(model)
     objective_dir = model.objective.direction

cobra/flux_analysis/phenotype_phase_plane.py

@@ -2,6 +2,7 @@
 
 from __future__ import absolute_import, division
 
+import logging
 from itertools import product
 
 import pandas as pd
@@ -12,10 +13,14 @@ from six import iteritems
 import cobra.util.solver as sutil
 from cobra.exceptions import OptimizationError
 from cobra.flux_analysis import flux_variability_analysis as fva
+from cobra.flux_analysis.helpers import normalize_cutoff
+
+
+LOGGER = logging.getLogger(__name__)
 
 
 def production_envelope(model, reactions, objective=None, carbon_sources=None,
-                        points=20, threshold=1e-7):
+                        points=20, threshold=None):
     """Calculate the objective value conditioned on all combinations of
     fluxes for a set of chosen reactions
 
@@ -46,7 +51,8 @@ def production_envelope(model, reactions, objective=None, carbon_sources=None,
     points : int, optional
        The number of points to calculate production for.
     threshold : float, optional
-        A cut-off under which flux values will be considered to be zero.
+        A cut-off under which flux values will be considered to be zero
+        (default model.tolerance).
 
     Returns
     -------
@@ -71,14 +77,18 @@ def production_envelope(model, reactions, objective=None, carbon_sources=None,
     >>> from cobra.flux_analysis import production_envelope
     >>> model = cobra.test.create_test_model("textbook")
     >>> production_envelope(model, ["EX_glc__D_e", "EX_o2_e"])
+
     """
+
     reactions = model.reactions.get_by_any(reactions)
     objective = model.solver.objective if objective is None else objective
     data = dict()
+
     if carbon_sources is None:
         c_input = find_carbon_sources(model)
     else:
         c_input = model.reactions.get_by_any(carbon_sources)
+
     if c_input is None:
         data['carbon_source'] = None
     elif hasattr(c_input, 'id'):
@@ -86,17 +96,23 @@ def production_envelope(model, reactions, objective=None, carbon_sources=None,
     else:
         data['carbon_source'] = ', '.join(rxn.id for rxn in c_input)
 
+    threshold = normalize_cutoff(model, threshold)
+
     size = points ** len(reactions)
+
     for direction in ('minimum', 'maximum'):
         data['flux_{}'.format(direction)] = full(size, nan, dtype=float)
         data['carbon_yield_{}'.format(direction)] = full(
             size, nan, dtype=float)
         data['mass_yield_{}'.format(direction)] = full(
             size, nan, dtype=float)
+
     grid = pd.DataFrame(data)
+
     with model:
         model.objective = objective
         objective_reactions = list(sutil.linear_reaction_coefficients(model))
+
         if len(objective_reactions) != 1:
             raise ValueError('cannot calculate yields for objectives with '
                              'multiple reactions')
@@ -110,6 +126,7 @@ def production_envelope(model, reactions, objective=None, carbon_sources=None,
         tmp = pd.DataFrame(points, columns=[rxn.id for rxn in reactions])
         grid = pd.concat([grid, tmp], axis=1, copy=False)
         add_envelope(model, reactions, grid, c_input, c_output, threshold)
+
     return grid
 
 
@@ -117,6 +134,7 @@ def add_envelope(model, reactions, grid, c_input, c_output, threshold):
     if c_input is not None:
         input_components = [reaction_elements(rxn) for rxn in c_input]
         output_components = reaction_elements(c_output)
+
         try:
             input_weights = [reaction_weight(rxn) for rxn in c_input]
             output_weight = reaction_weight(c_output)
@@ -132,16 +150,20 @@ def add_envelope(model, reactions, grid, c_input, c_output, threshold):
     for direction in ('minimum', 'maximum'):
         with model:
             model.objective_direction = direction
+
             for i in range(len(grid)):
                 with model:
                     for rxn in reactions:
                         point = grid.at[i, rxn.id]
                         rxn.bounds = point, point
                     obj_val = model.slim_optimize()
+
                     if model.solver.status != OPTIMAL:
                         continue
+
                     grid.at[i, 'flux_{}'.format(direction)] = \
                         0.0 if abs(obj_val) < threshold else obj_val
+
                     if c_input is not None:
                         grid.at[i, 'carbon_yield_{}'.format(direction)] = \
                             total_yield([rxn.flux for rxn in c_input],
@@ -213,10 +235,13 @@ def reaction_elements(reaction):
 
 def reaction_weight(reaction):
     """Return the metabolite weight times its stoichiometric coefficient."""
+
     if len(reaction.metabolites) != 1:
         raise ValueError('Reaction weight is only defined for single '
                          'metabolite products or educts.')
+
     met, coeff = next(iteritems(reaction.metabolites))
+
     return [coeff * met.formula_weight]
 
 
@@ -234,8 +259,10 @@ def total_components_flux(flux, components, consumption=True):
         Whether to sum up consumption or production fluxes.
 
     """
+
     direction = 1 if consumption else -1
     c_flux = [elem * flux * direction for elem in components]
+
     return sum([flux for flux in c_flux if flux > 0])
 
 
@@ -254,6 +281,7 @@ def find_carbon_sources(model):
        The medium reactions with carbon input flux.
 
     """
+
     try:
         model.slim_optimize(error_value=None)
     except OptimizationError:

cobra/flux_analysis/variability.py

@@ -14,6 +14,7 @@ from pandas import DataFrame
 from cobra.core import Configuration, get_solution
 from cobra.flux_analysis.deletion import (
     single_gene_deletion, single_reaction_deletion)
+from cobra.flux_analysis.helpers import normalize_cutoff
 from cobra.flux_analysis.loopless import loopless_fva_iter
 from cobra.flux_analysis.parsimonious import add_pfba
 from cobra.util import solver as sutil
@@ -201,7 +202,7 @@ def flux_variability_analysis(model, reaction_list=None, loopless=False,
 
 def find_blocked_reactions(model,
                            reaction_list=None,
-                           zero_cutoff=1e-9,
+                           zero_cutoff=None,
                            open_exchanges=False,
                            processes=None):
     """
@@ -223,7 +224,8 @@ def find_blocked_reactions(model,
         List of reactions to consider, the default includes all model
         reactions.
     zero_cutoff : float, optional
-        Flux value which is considered to effectively be zero.
+        Flux value which is considered to effectively be zero
+        (default model.tolerance).
     open_exchanges : bool, optional
         Whether or not to open all exchange reactions to very high flux ranges.
     processes : int, optional
@@ -235,7 +237,10 @@ def find_blocked_reactions(model,
     -------
     list
         List with the identifiers of blocked reactions.
+
     """
+    zero_cutoff = normalize_cutoff(model, zero_cutoff)
+
     with model:
         if open_exchanges:
             for reaction in model.exchanges:

cobra/io/__init__.py

@@ -5,7 +5,6 @@ from __future__ import absolute_import
 from cobra.io.dict import model_from_dict, model_to_dict
 from cobra.io.json import from_json, load_json_model, save_json_model, to_json
 from cobra.io.mat import load_matlab_model, save_matlab_model
-from cobra.io.sbml import read_legacy_sbml
-from cobra.io.sbml import write_cobra_model_to_sbml_file as write_legacy_sbml
-from cobra.io.sbml3 import read_sbml_model, write_sbml_model
+from cobra.io.sbml import read_sbml_model, write_sbml_model, \
+    validate_sbml_model
 from cobra.io.yaml import from_yaml, load_yaml_model, save_yaml_model, to_yaml