.github/ISSUE_TEMPLATE.md

@@ -8,8 +8,7 @@ Please explain:
 
 #### Code Sample
 
-Create a [minimal, complete, verifiable example
-](https://stackoverflow.com/help/mcve).
+Create a [minimal, complete, verifiable example](https://stackoverflow.com/help/mcve).
 
 ```python
 # Paste your code here.
@@ -20,9 +19,13 @@ Create a [minimal, complete, verifiable example
 
 #### Expected Output
 
-#### Output of `cobra.show_versions()`
+#### Dependency Information
+
+Please paste the output of `python -c "from cobra import show_versions;
+show_versions()"` in the details block below.
 
 <details>
-# Paste the output of `import cobra;cobra.show_versions()` here.
+
+[Paste output here.]
 
 </details>

.travis.yml

@@ -2,8 +2,7 @@ language: python
 python: 3.5
 sudo: false
 cache:
-  directories:
-    - $HOME/.cache/pip
+  pip: true
 git:
   depth: 2
 
@@ -15,8 +14,6 @@ branches:
 
 env:
   global:
-    - secure: "hkKBaGLvoDVgktSKR3BmX+mYlGzHw9EO11MRHtiH8D9BbdygOR9p9aSV/OxkaRWhnkSP5/0SXqVgBrvU1g5OsR6cc85UQSpJ5H5jVnLoWelIbTxMCikjxDSkZlseD7ZEWrKZjRo/ZN2qym0HRWpsir3qLpl8W25xHRv/sK7Z6g8="
-    - secure: "DflyBz+QiyhlhBxn4wN00xu248EJUMjKTxUZQN6wq22qV55xO3ToGJTy9i4D6OBfZGAlSXxjjKCJ2+0sAjsghBSDEK56ud3EEg/08TIo7/T8ex/C58FsGoGFz3yDBATmquClEWN8vAMrLdxwniHmQVCBZCP/phdt5dct0AUuDc8="
     - GITHUB_REPO=opencobra/cobrapy
 
 matrix:
@@ -80,7 +77,7 @@ deploy:
       repo: $GITHUB_REPO
       python: '3.6'
       condition: $TRAVIS_OS_NAME == "linux"
-    user: $PYPI_USERNAME
+    user: cobrapy
     password: $PYPI_PASSWORD
   - provider: script
     skip_cleanup: true

appveyor.yml

@@ -8,23 +8,18 @@ environment:
 
   global:
     PIP_CACHE_DIR: "pip_cache"
+    PIP_DISABLE_PIP_VERSION_CHECK: "yes"
 
   matrix:
     - PYTHON: "C:\\Miniconda-x64"
       CONDA: true
 
-    - PYTHON: "C:\\Miniconda35-x64"
-      CONDA: true
-
     - PYTHON: "C:\\Miniconda36-x64"
       CONDA: true
 
     - PYTHON: "C:\\Python27-x64"
       TOXENV: py27
 
-    - PYTHON: "C:\\Python35-x64"
-      TOXENV: py35
-
     - PYTHON: "C:\\Python36-x64"
       TOXENV: py36
 
@@ -46,6 +41,7 @@ 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 .

benchmarks/geometric_fba.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import logging
+```
+
+%% Cell type:code id: tags:
+
+``` python
+logging.basicConfig(level="DEBUG")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+from cobra.test import create_test_model
+```
+
+%% Output
+
+    DEBUG:optlang.util:Gurobi python bindings found at /home/moritz/.virtualenvs/cobra/lib/python3.5/site-packages/gurobipy
+    DEBUG:optlang.util:GLPK python bindings found at /home/moritz/.virtualenvs/cobra/lib/python3.5/site-packages/swiglpk
+    DEBUG:optlang.util:Mosek python bindings not available.
+    DEBUG:optlang.util:CPLEX python bindings found at /home/moritz/.virtualenvs/cobra/lib/python3.5/site-packages/cplex
+    DEBUG:optlang.util:Scipy solver not available
+    DEBUG:pip._internal.utils.misc:lzma module is not available
+    DEBUG:pip._internal.vcs:Registered VCS backend: git
+    DEBUG:pip._internal.vcs:Registered VCS backend: hg
+    DEBUG:pip._internal.vcs:Registered VCS backend: svn
+    DEBUG:pip._internal.vcs:Registered VCS backend: bzr
+
+%% Cell type:code id: tags:
+
+``` python
+from cobra.flux_analysis import geometric_fba
+```
+
+%% Cell type:code id: tags:
+
+``` python
+model = create_test_model("textbook")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "cplex"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+logging.getLogger().setLevel(logging.DEBUG)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%time geometric_fba(model)
+```
+
+%% Output
+
+    DEBUG:cobra.flux_analysis.geometric:Iteration: 1; delta: 4e-13; status: optimal.
+
+    CPU times: user 481 ms, sys: 0 ns, total: 481 ms
+    Wall time: 480 ms
+
+    <Solution 0.000 at 0x7f0c62d3f6a0>
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "gurobi"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%time geometric_fba(model)
+```
+
+%% Output
+
+    DEBUG:cobra.flux_analysis.geometric:Iteration: 1; delta: 6.47e-12; status: optimal.
+
+    CPU times: user 504 ms, sys: 2.91 ms, total: 507 ms
+    Wall time: 505 ms
+
+    <Solution 0.000 at 0x7f0c6288be10>
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "glpk"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%time geometric_fba(model)
+```
+
+%% Output
+
+    DEBUG:optlang.glpk_interface:Status undefined. GLPK status code returned by glp_simplex was 1
+    DEBUG:cobra.flux_analysis.geometric:Iteration: 1; delta: 2.96e-13; status: optimal.
+
+    CPU times: user 263 ms, sys: 3.46 ms, total: 267 ms
+    Wall time: 263 ms
+
+    <Solution -0.000 at 0x7f0c6283d9b0>

benchmarks/objective_benchmark.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from cobra.io import read_sbml_model
+```
+
+%% Cell type:code id: tags:
+
+``` python
+from optlang.symbolics import Zero, One, add, mul, Add, Mul, Real
+```
+
+%% Cell type:code id: tags:
+
+``` python
+one = Real(1.0)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+model = read_sbml_model("/home/moritz/Projects/memote/Models/iJO1366.xml.gz")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+import cobra.util.solver as sutil
+```
+
+%% Cell type:code id: tags:
+
+``` python
+from itertools import chain
+```
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "glpk"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+reaction_variables = ((rxn.forward_variable, rxn.reverse_variable)
+                      for rxn in model.reactions)
+variables = list(chain(*reaction_variables))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%%timeit
+with model:
+    model.objective = model.problem.Objective(
+        Zero, direction='min', sloppy=True,
+        name="_pfba_objective")
+    model.objective.set_linear_coefficients({v: 1.0 for v in variables})
+```
+
+%% Output
+
+    146 ms ± 3.94 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
+
+%% Cell type:code id: tags:
+
+``` python
+%%timeit
+with model:
+    model.objective = model.problem.Objective(
+        add([mul(one, v) for v in variables]),
+        direction='min', sloppy=True, name="_pfba_objective")
+```
+
+%% Output
+
+    149 ms ± 2.13 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "gurobi"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+reaction_variables = ((rxn.forward_variable, rxn.reverse_variable)
+                      for rxn in model.reactions)
+variables = list(chain(*reaction_variables))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%%timeit
+with model:
+    model.objective = model.problem.Objective(
+        Zero, direction='min', sloppy=True,
+        name="_pfba_objective")
+    model.objective.set_linear_coefficients({v: 1.0 for v in variables})
+```
+
+%% Output
+
+    13.7 ms ± 269 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
+
+%% Cell type:code id: tags:
+
+``` python
+%%timeit
+with model:
+    model.objective = model.problem.Objective(
+        add([mul(one, v) for v in variables]),
+        direction='min', sloppy=True, name="_pfba_objective")
+```
+
+%% Output
+
+    218 ms ± 2.55 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "cplex"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+reaction_variables = ((rxn.forward_variable, rxn.reverse_variable)
+                      for rxn in model.reactions)
+variables = list(chain(*reaction_variables))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%%timeit
+with model:
+    model.objective = model.problem.Objective(
+        Zero, direction='min', sloppy=True,
+        name="_pfba_objective")
+    model.objective.set_linear_coefficients({v: 1.0 for v in variables})
+```
+
+%% Output
+
+    2.04 s ± 32.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
+
+%% Cell type:code id: tags:
+
+``` python
+%%timeit
+with model:
+    model.objective = model.problem.Objective(
+        add([mul(one, v) for v in variables]),
+        direction='min', sloppy=True, name="_pfba_objective")
+```
+
+%% Output
+
+    183 ms ± 5.28 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

benchmarks/single_gene_deletion_linear_room.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import logging
+```
+
+%% Cell type:code id: tags:
+
+``` python
+logging.basicConfig(level="DEBUG")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+from cobra.test import create_test_model
+```
+
+%% Output
+
+    DEBUG:optlang.util:Gurobi python bindings found at /home/moritz/.virtualenvs/cobra/lib/python3.5/site-packages/gurobipy
+    DEBUG:optlang.util:GLPK python bindings found at /home/moritz/.virtualenvs/cobra/lib/python3.5/site-packages/swiglpk
+    DEBUG:optlang.util:Mosek python bindings not available.
+    DEBUG:optlang.util:CPLEX python bindings found at /home/moritz/.virtualenvs/cobra/lib/python3.5/site-packages/cplex
+    DEBUG:optlang.util:Scipy solver not available
+    DEBUG:pip._internal.utils.misc:lzma module is not available
+    DEBUG:pip._internal.vcs:Registered VCS backend: git
+    DEBUG:pip._internal.vcs:Registered VCS backend: hg
+    DEBUG:pip._internal.vcs:Registered VCS backend: svn
+    DEBUG:pip._internal.vcs:Registered VCS backend: bzr
+
+%% Cell type:code id: tags:
+
+``` python
+from cobra.flux_analysis import single_gene_deletion
+```
+
+%% Cell type:code id: tags:
+
+``` python
+model = create_test_model("textbook")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+genes = ['b0008', 'b0114', 'b2276', 'b1779']
+```
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "cplex"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+logging.getLogger().setLevel(logging.DEBUG)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%time single_gene_deletion(model, method="linear room", gene_list=genes, processes=1)
+```
+
+%% Output
+
+    CPU times: user 308 ms, sys: 0 ns, total: 308 ms
+    Wall time: 307 ms
+
+                   growth   status
+    ids
+    (b1779)  5.034815e+00  optimal
+    (b2276)  5.475871e+00  optimal
+    (b0114)  1.224159e+00  optimal
+    (b0008)  3.028591e-15  optimal
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "gurobi"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%time single_gene_deletion(model, method="linear room", gene_list=genes, processes=1)
+```
+
+%% Output
+
+    CPU times: user 271 ms, sys: 0 ns, total: 271 ms
+    Wall time: 265 ms
+
+               growth   status
+    ids
+    (b1779)  5.034815  optimal
+    (b2276)  5.475871  optimal
+    (b0114)  1.224159  optimal
+    (b0008)  0.000000  optimal
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "glpk"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver.configuration.timeout = 20
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%time single_gene_deletion(model, method="linear room", gene_list=genes, processes=1)
+```
+
+%% Output
+
+    DEBUG:optlang.glpk_interface:Status undefined. GLPK status code returned by glp_simplex was 1
+
+    CPU times: user 205 ms, sys: 4.99 ms, total: 210 ms
+    Wall time: 203 ms
+
+                   growth   status
+    ids
+    (b1779)  5.034815e+00  optimal
+    (b2276)  5.475871e+00  optimal
+    (b0114)  1.224159e+00  optimal
+    (b0008) -5.857507e-17  optimal
+
+%% Cell type:code id: tags:
+
+``` python
+```

benchmarks/single_gene_deletion_room.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import logging
+```
+
+%% Cell type:code id: tags:
+
+``` python
+logging.basicConfig(level="DEBUG")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+from cobra.test import create_test_model
+```
+
+%% Output
+
+    DEBUG:optlang.util:Gurobi python bindings found at /home/moritz/.virtualenvs/cobra/lib/python3.5/site-packages/gurobipy
+    DEBUG:optlang.util:GLPK python bindings found at /home/moritz/.virtualenvs/cobra/lib/python3.5/site-packages/swiglpk
+    DEBUG:optlang.util:Mosek python bindings not available.
+    DEBUG:optlang.util:CPLEX python bindings found at /home/moritz/.virtualenvs/cobra/lib/python3.5/site-packages/cplex
+    DEBUG:optlang.util:Scipy solver not available
+    DEBUG:pip._internal.utils.misc:lzma module is not available
+    DEBUG:pip._internal.vcs:Registered VCS backend: git
+    DEBUG:pip._internal.vcs:Registered VCS backend: hg
+    DEBUG:pip._internal.vcs:Registered VCS backend: svn
+    DEBUG:pip._internal.vcs:Registered VCS backend: bzr
+
+%% Cell type:code id: tags:
+
+``` python
+from cobra.flux_analysis import single_gene_deletion
+```
+
+%% Cell type:code id: tags:
+
+``` python
+model = create_test_model("textbook")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+genes = ['b0008', 'b0114', 'b2276', 'b1779']
+```
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "cplex"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+logging.getLogger().setLevel(logging.DEBUG)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%time single_gene_deletion(model, method="room", gene_list=genes, processes=1)
+```
+
+%% Output
+
+    WARNING:root:Warning:  No solution found from 6 MIP starts.
+    
+
+    CPU times: user 21.4 s, sys: 84.2 ms, total: 21.5 s
+    Wall time: 3.23 s
+
+             growth   status
+    ids
+    (b1779)    45.0  optimal
+    (b0008)     0.0  optimal
+    (b0114)    27.0  optimal
+    (b2276)    46.0  optimal
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "gurobi"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%time single_gene_deletion(model, method="room", gene_list=genes, processes=1)
+```
+
+%% Output
+
+    CPU times: user 1.97 s, sys: 4.23 ms, total: 1.98 s
+    Wall time: 1.97 s
+
+             growth   status
+    ids
+    (b1779)    45.0  optimal
+    (b0008)     0.0  optimal
+    (b0114)    27.0  optimal
+    (b2276)    46.0  optimal
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver = "glpk"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+model.solver.configuration.timeout = 20
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%time single_gene_deletion(model, method="room", gene_list=genes, processes=1)
+```
+
+%% Output
+
+    DEBUG:optlang.glpk_interface:Status undefined. GLPK status code returned by glp_simplex was 1
+
+    CPU times: user 1min 21s, sys: 53.1 ms, total: 1min 21s
+    Wall time: 1min 21s
+
+             growth      status
+    ids
+    (b1779)    45.0     optimal
+    (b0008)     0.0     optimal
+    (b0114)     NaN  time_limit
+    (b2276)     NaN  time_limit

cobra/__init__.py

@@ -11,9 +11,9 @@ from os.path import dirname as _dirname
 from cobra import flux_analysis, io
 from cobra.core import (
     DictList, Gene, Metabolite, Model, Object, Reaction, Species)
-from cobra.util.version_info import show_versions
+from cobra.util import show_versions
 
-__version__ = "0.11.3"
+__version__ = "0.13.3"
 
 # set the warning format to be prettier and fit on one line
 _cobra_path = _dirname(_abspath(__file__))

cobra/core/metabolite.py

@@ -126,7 +126,7 @@ class Metabolite(Species):
             return sum([count * elements_and_molecular_weights[element]
                         for element, count in self.elements.items()])
         except KeyError as e:
-            warn("The element %s does not appear in the peridic table" % e)
+            warn("The element %s does not appear in the periodic table" % e)
 
     @property
     def y(self):
@@ -207,38 +207,40 @@ class Metabolite(Species):
         """
         self._model.remove_metabolites(self, destructive)
 
-    def summary(self, solution=None, threshold=0.01, fva=False,
+    def summary(self, solution=None, threshold=0.01, fva=None, names=False,
                 floatfmt='.3g'):
-        """Print a summary of the reactions which produce and consume this
-        metabolite.
+        """
+        Print a summary of the production and consumption fluxes.
 
         This method requires the model for which this metabolite is a part
         to be solved.
 
         Parameters
         ----------
-        solution : cobra.core.Solution
+        solution : cobra.Solution, optional
             A previously solved model solution to use for generating the
             summary. If none provided (default), the summary method will
             resolve the model. Note that the solution object must match the
             model, i.e., changes to the model such as changed bounds,
             added or removed reactions are not taken into account by this
             method.
-
-        threshold : float
-            a value below which to ignore reaction fluxes
-
-        fva : float (0->1), or None
+        threshold : float, optional
+            Threshold below which fluxes are not reported.
+        fva : pandas.DataFrame, float or None, optional
             Whether or not to include flux variability analysis in the output.
-            If given, fva should be a float between 0 and 1, representing the
+            If given, fva should either be a previous FVA solution matching
+            the model or a float between 0 and 1 representing the
             fraction of the optimum objective to be searched.
+        names : bool, optional
+            Emit reaction and metabolite names rather than identifiers (default
+            False).
+        floatfmt : string, optional
+            Format string for floats (default '.3g').
 
-        floatfmt : string
-            format method for floats, passed to tabulate. Default is '.3g'.
         """
         from cobra.flux_analysis.summary import metabolite_summary
         return metabolite_summary(self, solution=solution, threshold=threshold,
-                                  fva=fva, floatfmt=floatfmt)
+                                  fva=fva, names=names, floatfmt=floatfmt)
 
     def _repr_html_(self):
         return """

cobra/core/model.py

@@ -23,6 +23,7 @@ from cobra.util.solver import (
     get_solver_name, interface_to_str, set_objective, solvers,
     add_cons_vars_to_problem, remove_cons_vars_from_problem, assert_optimal)
 from cobra.util.util import AutoVivification, format_long_string
+from cobra.medium import find_boundary_types
 
 LOGGER = logging.getLogger(__name__)
 
@@ -233,16 +234,20 @@ class Model(Object):
 
         # Set the given media bounds
         media_rxns = list()
+        exchange_rxns = frozenset(self.exchanges)
         for rxn_id, bound in iteritems(medium):
             rxn = self.reactions.get_by_id(rxn_id)
+            if rxn not in exchange_rxns:
+                LOGGER.warn("%s does not seem to be an"
+                            " an exchange reaction. Applying bounds anyway.",
+                            rxn.id)
             media_rxns.append(rxn)
             set_active_bound(rxn, bound)
 
-        boundary_rxns = set(self.exchanges)
-        media_rxns = set(media_rxns)
+        media_rxns = frozenset(media_rxns)
 
         # Turn off reactions not present in media
-        for rxn in (boundary_rxns - media_rxns):
+        for rxn in (exchange_rxns - media_rxns):
             set_active_bound(rxn, 0)
 
     def __add__(self, other_model):
@@ -725,13 +730,36 @@ class Model(Object):
         """
         return self.solver.constraints
 
+    @property
+    def boundary(self):
+        """Boundary reactions in the model.
+        Reactions that either have no substrate or product.
+        """
+        return [rxn for rxn in self.reactions if rxn.boundary]
+
     @property
     def exchanges(self):
         """Exchange reactions in model.
+        Reactions that exchange mass with the exterior. Uses annotations
+        and heuristics to exclude non-exchanges such as sink reactions.
+        """
+        return find_boundary_types(self, "exchange", None)
 
-        Reactions that either don't have products or substrates.
+    @property
+    def demands(self):
+        """Demand reactions in model.
+        Irreversible reactions that accumulate or consume a metabolite in
+        the inside of the model.
         """
-        return [rxn for rxn in self.reactions if rxn.boundary]
+        return find_boundary_types(self, "demand", None)
+
+    @property
+    def sinks(self):
+        """Sink reactions in model.
+        Reversible reactions that accumulate or consume a metabolite in
+        the inside of the model.
+        """
+        return find_boundary_types(self, "sink", None)
 
     def _populate_solver(self, reaction_list, metabolite_list=None):
         """Populate attached solver with constraints and variables that
@@ -934,34 +962,37 @@ class Model(Object):
         else:
             raise ValueError("Unknown objective direction '{}'.".format(value))
 
-    def summary(self, solution=None, threshold=1E-8, fva=None, floatfmt='.3g'):
-        """Print a summary of the input and output fluxes of the model. This
-        method requires the model to have been previously solved.
+    def summary(self, solution=None, threshold=1E-06, fva=None, names=False,
+                floatfmt='.3g'):
+        """
+        Print a summary of the input and output fluxes of the model.
 
         Parameters
         ----------
-        solution: cobra.core.Solution
+        solution: cobra.Solution, optional
             A previously solved model solution to use for generating the
             summary. If none provided (default), the summary method will
             resolve the model. Note that the solution object must match the
             model, i.e., changes to the model such as changed bounds,
             added or removed reactions are not taken into account by this
             method.
-
-        threshold : float
-            tolerance for determining if a flux is zero (not printed)
-
-        fva : int or None
-            Whether or not to calculate and report flux variability in the
-            output summary
-
-        floatfmt : string
-            format method for floats, passed to tabulate. Default is '.3g'.
+        threshold : float, optional
+            Threshold below which fluxes are not reported.
+        fva : pandas.DataFrame, float or None, optional
+            Whether or not to include flux variability analysis in the output.
+            If given, fva should either be a previous FVA solution matching
+            the model or a float between 0 and 1 representing the
+            fraction of the optimum objective to be searched.
+        names : bool, optional
+            Emit reaction and metabolite names rather than identifiers (default
+            False).
+        floatfmt : string, optional
+            Format string for floats (default '.3g').
 
         """
         from cobra.flux_analysis.summary import model_summary
         return model_summary(self, solution=solution, threshold=threshold,
-                             fva=fva, floatfmt=floatfmt)
+                             fva=fva, names=names, floatfmt=floatfmt)
 
     def __enter__(self):
         """Record all future changes to the model, undoing them when a call to

cobra/core/reaction.py

@@ -739,6 +739,14 @@ class Reaction(Object):
         _id_to_metabolites = dict([(x.id, x) for x in self._metabolites])
 
         for metabolite, coefficient in iteritems(metabolites_to_add):
+
+            # Make sure metabolites being added belong to the same model, or
+            # else copy them.
+            if isinstance(metabolite, Metabolite):
+                if ((metabolite.model is not None) and
+                        (metabolite.model is not self._model)):
+                    metabolite = metabolite.copy()
+
             met_id = str(metabolite)
             # If a metabolite already exists in the reaction then
             # just add them.

cobra/flux_analysis/__init__.py

 # -*- coding: utf-8 -*-
 
 from cobra.flux_analysis.gapfilling import gapfill
+from cobra.flux_analysis.geometric import geometric_fba
 from cobra.flux_analysis.loopless import (
     construct_loopless_model, loopless_solution, add_loopless)
 from cobra.flux_analysis.parsimonious import pfba
+from cobra.flux_analysis.moma import moma, add_moma
+from cobra.flux_analysis.room import room, add_room
 from cobra.flux_analysis.deletion import (
     single_gene_deletion, single_reaction_deletion)
 from cobra.flux_analysis.variability import (

cobra/flux_analysis/deletion.py

@@ -13,6 +13,7 @@ import pandas as pd
 from cobra.manipulation.delete import find_gene_knockout_reactions
 import cobra.util.solver as sutil
 from cobra.flux_analysis.moma import add_moma
+from cobra.flux_analysis.room import add_room
 
 LOGGER = logging.getLogger(__name__)
 
@@ -72,7 +73,7 @@ def _init_worker(model):
 
 
 def _multi_deletion(model, entity, element_lists, method="fba",
-                    processes=None):
+                    solution=None, processes=None, **kwargs):
     """
     Provide a common interface for single or multiple knockouts.
 
@@ -80,21 +81,21 @@ def _multi_deletion(model, entity, element_lists, method="fba",
     ----------
     model : cobra.Model
         The metabolic model to perform deletions in.
-
     entity : 'gene' or 'reaction'
         The entity to knockout (``cobra.Gene`` or ``cobra.Reaction``).
-
     element_lists : list
         List of iterables ``cobra.Reaction``s or ``cobra.Gene``s (or their IDs)
         to be deleted.
-
-    method: {"fba", "moma", "linear moma"}, optional
+    method: {"fba", "moma", "linear moma", "room", "linear room"}, optional
         Method used to predict the growth rate.
-
+    solution : cobra.Solution, optional
+        A previous solution to use as a reference for (linear) MOMA or ROOM.
     processes : int, optional
         The number of parallel processes to run. Can speed up the computations
         if the number of knockouts to perform is large. If not passed,
         will be set to the number of CPUs found.
+    kwargs :
+        Passed on to underlying simulation functions.
 
     Returns
     -------
@@ -110,7 +111,7 @@ def _multi_deletion(model, entity, element_lists, method="fba",
             The solution's status.
     """
     solver = sutil.interface_to_str(model.problem.__name__)
-    if "moma" in method and solver not in sutil.qp_solvers:
+    if method == "moma" and solver not in sutil.qp_solvers:
         raise RuntimeError(
             "Cannot use MOMA since '{}' is not QP-capable."
             "Please choose a different solver or use FBA only.".format(solver))
@@ -124,7 +125,10 @@ def _multi_deletion(model, entity, element_lists, method="fba",
 
     with model:
         if "moma" in method:
-            add_moma(model, linear="linear" in method)
+            add_moma(model, solution=solution, linear="linear" in method)
+        elif "room" in method:
+            add_room(model, solution=solution, linear="linear" in method,
+                     **kwargs)
 
         args = set([frozenset(comb) for comb in product(*element_lists)])
         processes = min(processes, len(args))
@@ -180,7 +184,7 @@ def _element_lists(entities, *ids):
 
 
 def single_reaction_deletion(model, reaction_list=None, method="fba",
-                             processes=None):
+                             solution=None, processes=None, **kwargs):
     """
     Knock out each reaction from a given list.
 
@@ -188,18 +192,20 @@ def single_reaction_deletion(model, reaction_list=None, method="fba",
     ----------
     model : cobra.Model
         The metabolic model to perform deletions in.
-
-    reaction_list : iterable
+    reaction_list : iterable, optional
         ``cobra.Reaction``s to be deleted. If not passed,
         all the reactions from the model are used.
-
-    method: {"fba", "moma", "linear moma"}, optional
+    method: {"fba", "moma", "linear moma", "room", "linear room"}, optional
         Method used to predict the growth rate.
-
+    solution : cobra.Solution, optional
+        A previous solution to use as a reference for (linear) MOMA or ROOM.
     processes : int, optional
         The number of parallel processes to run. Can speed up the computations
         if the number of knockouts to perform is large. If not passed,
         will be set to the number of CPUs found.
+    kwargs :
+        Keyword arguments are passed on to underlying simulation functions
+        such as ``add_room``.
 
     Returns
     -------
@@ -218,10 +224,11 @@ def single_reaction_deletion(model, reaction_list=None, method="fba",
     return _multi_deletion(
         model, 'reaction',
         element_lists=_element_lists(model.reactions, reaction_list),
-        method=method, processes=processes)
+        method=method, solution=solution, processes=processes, **kwargs)
 
 
-def single_gene_deletion(model, gene_list=None, method="fba", processes=None):
+def single_gene_deletion(model, gene_list=None, method="fba", solution=None,
+                         processes=None, **kwargs):
     """
     Knock out each gene from a given list.
 
@@ -229,18 +236,20 @@ def single_gene_deletion(model, gene_list=None, method="fba", processes=None):
     ----------
     model : cobra.Model
         The metabolic model to perform deletions in.
-
     gene_list : iterable
         ``cobra.Gene``s to be deleted. If not passed,
         all the genes from the model are used.
-
-    method: {"fba", "moma", "linear moma"}, optional
+    method: {"fba", "moma", "linear moma", "room", "linear room"}, optional
         Method used to predict the growth rate.
-
+    solution : cobra.Solution, optional
+        A previous solution to use as a reference for (linear) MOMA or ROOM.
     processes : int, optional
         The number of parallel processes to run. Can speed up the computations
         if the number of knockouts to perform is large. If not passed,
         will be set to the number of CPUs found.
+    kwargs :
+        Keyword arguments are passed on to underlying simulation functions
+        such as ``add_room``.
 
     Returns
     -------
@@ -258,11 +267,12 @@ def single_gene_deletion(model, gene_list=None, method="fba", processes=None):
     """
     return _multi_deletion(
         model, 'gene', element_lists=_element_lists(model.genes, gene_list),
-        method=method, processes=processes)
+        method=method, solution=solution, processes=processes, **kwargs)
 
 
 def double_reaction_deletion(model, reaction_list1=None, reaction_list2=None,
-                             method="fba", processes=None):
+                             method="fba", solution=None, processes=None,
+                             **kwargs):
     """
     Knock out each reaction pair from the combinations of two given lists.
 
@@ -272,22 +282,23 @@ def double_reaction_deletion(model, reaction_list1=None, reaction_list2=None,
     ----------
     model : cobra.Model
         The metabolic model to perform deletions in.
-
     reaction_list1 : iterable, optional
         First iterable of ``cobra.Reaction``s to be deleted. If not passed,
         all the reactions from the model are used.
-
     reaction_list2 : iterable, optional
         Second iterable of ``cobra.Reaction``s to be deleted. If not passed,
         all the reactions from the model are used.
-
-    method: {"fba", "moma", "linear moma"}, optional
+    method: {"fba", "moma", "linear moma", "room", "linear room"}, optional
         Method used to predict the growth rate.
-
+    solution : cobra.Solution, optional
+        A previous solution to use as a reference for (linear) MOMA or ROOM.
     processes : int, optional
         The number of parallel processes to run. Can speed up the computations
         if the number of knockouts to perform is large. If not passed,
         will be set to the number of CPUs found.
+    kwargs :
+        Keyword arguments are passed on to underlying simulation functions
+        such as ``add_room``.
 
     Returns
     -------
@@ -309,11 +320,11 @@ def double_reaction_deletion(model, reaction_list1=None, reaction_list2=None,
                                                     reaction_list2)
     return _multi_deletion(
         model, 'reaction', element_lists=[reaction_list1, reaction_list2],
-        method=method, processes=processes)
+        method=method, solution=solution, processes=processes, **kwargs)
 
 
-def double_gene_deletion(model, gene_list1=None, gene_list2=None,
-                         method="fba", processes=None):
+def double_gene_deletion(model, gene_list1=None, gene_list2=None, method="fba",
+                         solution=None, processes=None, **kwargs):
     """
     Knock out each gene pair from the combination of two given lists.
 
@@ -323,22 +334,23 @@ def double_gene_deletion(model, gene_list1=None, gene_list2=None,
     ----------
     model : cobra.Model
         The metabolic model to perform deletions in.
-
     gene_list1 : iterable, optional
         First iterable of ``cobra.Gene``s to be deleted. If not passed,
         all the genes from the model are used.
-
     gene_list2 : iterable, optional
         Second iterable of ``cobra.Gene``s to be deleted. If not passed,
         all the genes from the model are used.
-
-    method: {"fba", "moma", "linear moma"}, optional
+    method: {"fba", "moma", "linear moma", "room", "linear room"}, optional
         Method used to predict the growth rate.
-
+    solution : cobra.Solution, optional
+        A previous solution to use as a reference for (linear) MOMA or ROOM.
     processes : int, optional
         The number of parallel processes to run. Can speed up the computations
         if the number of knockouts to perform is large. If not passed,
         will be set to the number of CPUs found.
+    kwargs :
+        Keyword arguments are passed on to underlying simulation functions
+        such as ``add_room``.
 
     Returns
     -------
@@ -359,4 +371,4 @@ def double_gene_deletion(model, gene_list1=None, gene_list2=None,
                                             gene_list2)
     return _multi_deletion(
         model, 'gene', element_lists=[gene_list1, gene_list2],
-        method=method, processes=processes)
+        method=method, solution=solution, processes=processes, **kwargs)

cobra/flux_analysis/gapfilling.py

@@ -2,7 +2,7 @@
 
 from __future__ import absolute_import
 
-from optlang.symbolics import Add
+from optlang.symbolics import add, Zero
 
 from optlang.interface import OPTIMAL
 from cobra.core import Model
@@ -183,7 +183,9 @@ class GapFiller(object):
         self.model.add_cons_vars(self.indicators)
         self.model.add_cons_vars(constraints, sloppy=True)
         self.model.objective = prob.Objective(
-            Add(*self.indicators), direction='min')
+            Zero, direction='min', sloppy=True)
+        self.model.objective.set_linear_coefficients({
+            i: 1 for i in self.indicators})
         self.update_costs()
 
     def fill(self, iterations=1):

cobra/flux_analysis/geometric.py

+# -*- coding: utf8 -*-
+
+"""Provide an implementation of geometric FBA."""
+
+from __future__ import absolute_import, division
+
+import logging
+
+from optlang.symbolics import Zero
+
+from cobra.flux_analysis.parsimonious import add_pfba
+from cobra.flux_analysis.variability import flux_variability_analysis
+
+LOGGER = logging.getLogger(__name__)
+
+
+def geometric_fba(model, epsilon=1E-06, max_tries=200):
+    """
+    Perform geometric FBA to obtain a unique, centered flux distribution.
+
+    Geometric FBA [1]_ formulates the problem as a polyhedron and
+    then solves it by bounding the convex hull of the polyhedron.
+    The bounding forms a box around the convex hull which reduces
+    with every iteration and extracts a unique solution in this way.
+
+    Parameters
+    ----------
+    model: cobra.Model
+        The model to perform geometric FBA on.
+    epsilon: float, optional
+        The convergence tolerance of the model (default 1E-06).
+    max_tries: int, optional
+        Maximum number of iterations (default 200).
+
+    Returns
+    -------
+    cobra.Solution
+        The solution object containing all the constraints required
+        for geometric FBA.
+
+    References
+    ----------
+    .. [1] Smallbone, Kieran & Simeonidis, Vangelis. (2009).
+           Flux balance analysis: A geometric perspective.
+           Journal of theoretical biology.258. 311-5.
+           10.1016/j.jtbi.2009.01.027.
+
+    """
+
+    with model:
+        # Variables' and constraints' storage variables.
+        consts = []
+        obj_vars = []
+        updating_vars_cons = []
+
+        # The first iteration.
+        prob = model.problem
+        add_pfba(model)  # Minimize the solution space to a convex hull.
+        model.optimize()
+        fva_sol = flux_variability_analysis(model)
+        mean_flux = (fva_sol["maximum"] + fva_sol["minimum"]).abs() / 2
+
+        # Set the gFBA constraints.
+        for rxn in model.reactions:
+            var = prob.Variable("geometric_fba_" + rxn.id,
+                                lb=0,
+                                ub=mean_flux[rxn.id])
+            upper_const = prob.Constraint(rxn.flux_expression - var,
+                                          ub=mean_flux[rxn.id],
+                                          name="geometric_fba_upper_const_" +
+                                          rxn.id)
+            lower_const = prob.Constraint(rxn.flux_expression + var,
+                                          lb=fva_sol.at[rxn.id, "minimum"],
+                                          name="geometric_fba_lower_const_" +
+                                          rxn.id)
+            updating_vars_cons.append((rxn.id, var, upper_const, lower_const))
+            consts.extend([var, upper_const, lower_const])
+            obj_vars.append(var)
+        model.add_cons_vars(consts)
+
+        # Minimize the distance between the flux distribution and center.
+        model.objective = prob.Objective(Zero, sloppy=True, direction="min")
+        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)
+        mean_flux = (fva_sol["maximum"] + fva_sol["minimum"]).abs() / 2
+        delta = (fva_sol["maximum"] - fva_sol["minimum"]).max()
+        count = 1
+        LOGGER.debug("Iteration: %d; delta: %.3g; status: %s.",
+                     count, delta, sol.status)
+
+        # Following iterations that minimize the distance below threshold.
+        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)
+            mean_flux = (fva_sol["maximum"] + fva_sol["minimum"]).abs() / 2
+            delta = (fva_sol["maximum"] - fva_sol["minimum"]).max()
+            count += 1
+            LOGGER.debug("Iteration: %d; delta: %.3g; status: %s.",
+                         count, delta, sol.status)
+
+        if count == max_tries:
+            raise RuntimeError(
+                "The iterations have exceeded the maximum value of {}. "
+                "This is probably due to the increased complexity of the "
+                "model and can lead to inaccurate results. Please set a "
+                "different convergence tolerance and/or increase the "
+                "maximum iterations".format(max_tries)
+            )
+
+    return sol

cobra/flux_analysis/loopless.py

@@ -82,7 +82,8 @@ def add_loopless(model, zero_cutoff=1e-12):
 
 def _add_cycle_free(model, fluxes):
     """Add constraints for CycleFreeFlux."""
-    model.objective = Zero
+    model.objective = model.solver.interface.Objective(
+        Zero, direction="min", sloppy=True)
     objective_vars = []
     for rxn in model.reactions:
         flux = fluxes[rxn.id]
@@ -90,14 +91,13 @@ def _add_cycle_free(model, fluxes):
             rxn.bounds = (flux, flux)
             continue
         if flux >= 0:
-            rxn.lower_bound = max(0, rxn.lower_bound)
+            rxn.bounds = max(0, rxn.lower_bound), max(flux, rxn.upper_bound)
             objective_vars.append(rxn.forward_variable)
         else:
-            rxn.upper_bound = min(0, rxn.upper_bound)
+            rxn.bounds = min(flux, rxn.lower_bound), min(0, rxn.upper_bound)
             objective_vars.append(rxn.reverse_variable)
 
-    model.objective.set_linear_coefficients(dict.fromkeys(objective_vars, 1.0))
-    model.objective.direction = "min"
+    model.objective.set_linear_coefficients({v: 1.0 for v in objective_vars})
 
 
 def loopless_solution(model, fluxes=None):
@@ -208,11 +208,11 @@ def loopless_fva_iter(model, reaction, solution=False, zero_cutoff=1e-6):
 
     with model:
         _add_cycle_free(model, sol.fluxes)
-        flux = model.slim_optimize()
+        model.slim_optimize()
 
         # If the previous optimum is maintained in the loopless solution it was
         # loopless and we are done
-        if abs(flux - current) < zero_cutoff:
+        if abs(reaction.flux - current) < zero_cutoff:
             if solution:
                 return sol
             return current
@@ -221,11 +221,11 @@ def loopless_fva_iter(model, reaction, solution=False, zero_cutoff=1e-6):
         # almost loopless solution containing only loops including the current
         # reaction. Than remove all of those loops.
         ll_sol = get_solution(model).fluxes
-        bounds = reaction.bounds
         reaction.bounds = (current, current)
         model.slim_optimize()
         almost_ll_sol = get_solution(model).fluxes
-        reaction.bounds = bounds
+
+    with model:
         # find the reactions with loops using the current reaction and remove
         # the loops
         for rxn in model.reactions:

cobra/flux_analysis/moma.py

 # -*- coding: utf-8 -*-
 
-"""Contains functions to run minimization of metabolic adjustment (MOMA)."""
+"""Provide minimization of metabolic adjustment (MOMA)."""
 
 from __future__ import absolute_import
 
-from optlang.symbolics import Zero
+from optlang.symbolics import Zero, add
 
 import cobra.util.solver as sutil
 
+from cobra.flux_analysis.parsimonious import pfba
 
-def add_moma(model, solution=None, linear=False):
+
+def moma(model, solution=None, linear=True):
+    """
+    Compute a single solution based on (linear) MOMA.
+
+    Compute a new flux distribution that is at a minimal distance to a
+    previous reference solution. Minimization of metabolic adjustment (MOMA) is
+    generally used to assess the impact
+    of knock-outs. Thus the typical usage is to provide a wildtype flux
+    distribution as reference and a model in knock-out state.
+
+    Parameters
+    ----------
+    model : cobra.Model
+        The model state to compute a MOMA-based solution for.
+    solution : cobra.Solution, optional
+        A (wildtype) reference solution.
+    linear : bool, optional
+        Whether to use the linear MOMA formulation or not (default True).
+
+    Returns
+    -------
+    cobra.Solution
+        A flux distribution that is at a minimal distance compared to the
+        reference solution.
+
+    See Also
+    --------
+    add_moma : add MOMA constraints and objective
+
+    """
+    with model:
+        add_moma(model=model, solution=solution, linear=linear)
+        solution = model.optimize()
+    return solution
+
+
+def add_moma(model, solution=None, linear=True):
     r"""Add constraints and objective representing for MOMA.
 
     This adds variables and constraints for the minimization of metabolic
@@ -19,18 +57,15 @@ def add_moma(model, solution=None, linear=False):
     ----------
     model : cobra.Model
         The model to add MOMA constraints and objective to.
-    solution : cobra.Solution
-        A previous solution to use as a reference.
-    linear : bool
-        Whether to use the linear MOMA formulation or not.
-
-    Returns
-    -------
-    Nothing.
+    solution : cobra.Solution, optional
+        A previous solution to use as a reference. If no solution is given,
+        one will be computed using pFBA.
+    linear : bool, optional
+        Whether to use the linear MOMA formulation or not (default True).
 
     Notes
     -----
-    In the original MOMA specification one looks for the flux distribution
+    In the original MOMA [1]_ specification one looks for the flux distribution
     of the deletion (v^d) closest to the fluxes without the deletion (v).
     In math this means:
 
@@ -46,20 +81,35 @@ def add_moma(model, solution=None, linear=False):
          v^t = v^d_i - v_i
          lb_i <= v^d_i <= ub_i
 
-    So basically we just re-center the flux space at the old solution and than
+    So basically we just re-center the flux space at the old solution and then
     find the flux distribution closest to the new zero (center). This is the
     same strategy as used in cameo.
 
-    In the case of linear MOMA, we instead minimize \sum_i abs(v^t_i). The
+    In the case of linear MOMA [2]_, we instead minimize \sum_i abs(v^t_i). The
     linear MOMA is typically significantly faster. Also quadratic MOMA tends
     to give flux distributions in which all fluxes deviate from the reference
     fluxes a little bit whereas linear MOMA tends to give flux distributions
-    where the majority of fluxes are the same reference which few fluxes
+    where the majority of fluxes are the same reference with few fluxes
     deviating a lot (typical effect of L2 norm vs L1 norm).
 
     The former objective function is saved in the optlang solver interface as
-    "moma_old_objective" and this can be used to immediately extract the value
-    of the former objective after MOMA optimization.
+    ``"moma_old_objective"`` and this can be used to immediately extract the
+    value of the former objective after MOMA optimization.
+
+    See Also
+    --------
+    pfba : parsimonious FBA
+
+    References
+    ----------
+    .. [1] Segrè, Daniel, Dennis Vitkup, and George M. Church. “Analysis of
+           Optimality in Natural and Perturbed Metabolic Networks.”
+           Proceedings of the National Academy of Sciences 99, no. 23
+           (November 12, 2002): 15112. https://doi.org/10.1073/pnas.232349399.
+    .. [2] Becker, Scott A, Adam M Feist, Monica L Mo, Gregory Hannum,
+           Bernhard Ø Palsson, and Markus J Herrgard. “Quantitative
+           Prediction of Cellular Metabolism with Constraint-Based Models:
+           The COBRA Toolbox.” Nature Protocols 2 (March 29, 2007): 727.
     """
     if 'moma_old_objective' in model.solver.variables:
         raise ValueError('model is already adjusted for MOMA')
@@ -69,13 +119,14 @@ def add_moma(model, solution=None, linear=False):
         model.solver = sutil.choose_solver(model, qp=True)
 
     if solution is None:
-        solution = model.optimize()
+        solution = pfba(model)
     prob = model.problem
     v = prob.Variable("moma_old_objective")
     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 = Zero
+    model.objective = prob.Objective(Zero, direction="min", sloppy=True)
+    obj_vars = []
     for r in model.reactions:
         flux = solution.fluxes[r.id]
         if linear:
@@ -83,12 +134,16 @@ def add_moma(model, solution=None, linear=False):
                 model, r.flux_expression, name="moma_dist_" + r.id,
                 difference=flux, add=False)
             to_add.extend(components)
-            new_obj += components.variable
+            obj_vars.append(components.variable)
         else:
             dist = prob.Variable("moma_dist_" + r.id)
             const = prob.Constraint(r.flux_expression - dist, lb=flux, ub=flux,
                                     name="moma_constraint_" + r.id)
             to_add.extend([dist, const])
-            new_obj += dist**2
+            obj_vars.append(dist ** 2)
     model.add_cons_vars(to_add)
-    model.objective = prob.Objective(new_obj, direction='min')
+    if linear:
+        model.objective.set_linear_coefficients({v: 1.0 for v in obj_vars})
+    else:
+        model.objective = prob.Objective(
+            add(obj_vars), direction="min", sloppy=True)

cobra/flux_analysis/parsimonious.py

@@ -100,6 +100,5 @@ def add_pfba(model, objective=None, fraction_of_optimum=1.0):
                           for rxn in model.reactions)
     variables = chain(*reaction_variables)
     model.objective = model.problem.Objective(
-        Zero, direction='min', sloppy=True,
-        name="_pfba_objective")
-    model.objective.set_linear_coefficients(dict.fromkeys(variables, 1.0))
+        Zero, direction='min', sloppy=True, name="_pfba_objective")
+    model.objective.set_linear_coefficients({v: 1.0 for v in variables})

cobra/flux_analysis/room.py

+# -*- coding: utf-8 -*-
+
+"""Provide regulatory on/off minimization (ROOM)."""
+
+from __future__ import absolute_import, division
+
+from optlang.symbolics import Zero
+
+from cobra.flux_analysis.parsimonious import pfba
+
+
+def room(model, solution=None, linear=False, delta=0.03, epsilon=1E-03):
+    """
+    Compute a single solution based on regulatory on/off minimization (ROOM).
+
+    Compute a new flux distribution that minimizes the number of active
+    reactions needed to accommodate a previous reference solution.
+    Regulatory on/off minimization (ROOM) is generally used to assess the
+    impact of knock-outs. Thus the typical usage is to provide a wildtype flux
+    distribution as reference and a model in knock-out state.
+
+    Parameters
+    ----------
+    model : cobra.Model
+        The model state to compute a ROOM-based solution for.
+    solution : cobra.Solution, optional
+        A (wildtype) reference solution.
+    linear : bool, optional
+        Whether to use the linear ROOM formulation or not (default False).
+    delta: float, optional
+        The relative tolerance range (additive) (default 0.03).
+    epsilon: float, optional
+        The absolute tolerance range (multiplicative) (default 0.001).
+
+    Returns
+    -------
+    cobra.Solution
+        A flux distribution with minimal active reaction changes compared to
+        the reference.
+
+    See Also
+    --------
+    add_room : add ROOM constraints and objective
+
+    """
+    with model:
+        add_room(model=model, solution=solution, linear=linear, delta=delta,
+                 epsilon=epsilon)
+        solution = model.optimize()
+    return solution
+
+
+def add_room(model, solution=None, linear=False, delta=0.03, epsilon=1E-03):
+    r"""
+    Add constraints and objective for ROOM.
+
+    This function adds variables and constraints for applying regulatory
+    on/off minimization (ROOM) to the model.
+
+    Parameters
+    ----------
+    model : cobra.Model
+        The model to add ROOM constraints and objective to.
+    solution : cobra.Solution, optional
+        A previous solution to use as a reference. If no solution is given,
+        one will be computed using pFBA.
+    linear : bool, optional
+        Whether to use the linear ROOM formulation or not (default False).
+    delta: float, optional
+        The relative tolerance range which is additive in nature
+        (default 0.03).
+    epsilon: float, optional
+        The absolute range of tolerance which is multiplicative
+        (default 0.001).
+
+    Notes
+    -----
+    The formulation used here is the same as stated in the original paper [1]_.
+    The mathematical expression is given below:
+
+    minimize \sum_{i=1}^m y^i
+    s.t. Sv = 0
+         v_min <= v <= v_max
+         v_j = 0
+         j ∈ A
+         for 1 <= i <= m
+         v_i - y_i(v_{max,i} - w_i^u) <= w_i^u      (1)
+         v_i - y_i(v_{min,i} - w_i^l) <= w_i^l      (2)
+         y_i ∈ {0,1}                                (3)
+         w_i^u = w_i + \delta|w_i| + \epsilon
+         w_i^l = w_i - \delta|w_i| - \epsilon
+
+    So, for the linear version of the ROOM , constraint (3) is relaxed to
+    0 <= y_i <= 1.
+
+    See Also
+    --------
+    pfba : parsimonious FBA
+
+    References
+    ----------
+    .. [1] Tomer Shlomi, Omer Berkman and Eytan Ruppin, "Regulatory on/off
+     minimization of metabolic flux changes after genetic perturbations",
+     PNAS 2005 102 (21) 7695-7700; doi:10.1073/pnas.0406346102
+
+    """
+
+    if 'room_old_objective' in model.solver.variables:
+        raise ValueError('model is already adjusted for ROOM')
+
+    # optimizes if no reference solution is provided
+    if solution is None:
+        solution = pfba(model)
+
+    prob = model.problem
+    variable = prob.Variable("room_old_objective", ub=solution.objective_value)
+    constraint = prob.Constraint(
+        model.solver.objective.expression - variable,
+        ub=0.0,
+        lb=0.0,
+        name="room_old_objective_constraint"
+    )
+    model.objective = prob.Objective(Zero, direction="min", sloppy=True)
+    vars_and_cons = [variable, constraint]
+    obj_vars = []
+
+    for rxn in model.reactions:
+        flux = solution.fluxes[rxn.id]
+
+        if linear:
+            y = prob.Variable("y_" + rxn.id, lb=0, ub=1)
+            delta = epsilon = 0.0
+        else:
+            y = prob.Variable("y_" + rxn.id, type="binary")
+
+        # upper constraint
+        w_u = flux + (delta * abs(flux)) + epsilon
+        upper_const = prob.Constraint(
+            rxn.flux_expression - y * (rxn.upper_bound - w_u),
+            ub=w_u, name="room_constraint_upper_" + rxn.id)
+        # lower constraint
+        w_l = flux - (delta * abs(flux)) - epsilon
+        lower_const = prob.Constraint(
+            rxn.flux_expression - y * (rxn.lower_bound - w_l),
+            lb=w_l, name="room_constraint_lower_" + rxn.id)
+        vars_and_cons.extend([y, upper_const, lower_const])
+        obj_vars.append(y)
+
+    model.add_cons_vars(vars_and_cons)
+    model.objective.set_linear_coefficients({v: 1.0 for v in obj_vars})

cobra/flux_analysis/summary.py

 # -*- coding: utf-8 -*-
 
-from __future__ import absolute_import
+from __future__ import absolute_import, division
+
+import logging
+from operator import attrgetter
 
 import pandas as pd
 from numpy import zeros
@@ -13,68 +16,93 @@ from cobra.util.solver import linear_reaction_coefficients
 from cobra.util.util import format_long_string
 from cobra.core import get_solution
 
+LOGGER = logging.getLogger(__name__)
 
-def metabolite_summary(met, solution=None, threshold=0.01, fva=False,
-                       floatfmt='.3g'):
-    """Print a summary of the reactions which produce and consume this
-    metabolite
 
-    solution : cobra.core.Solution
-        A previously solved model solution to use for generating the
-        summary. If none provided (default), the summary method will resolve
-        the model. Note that the solution object must match the model, i.e.,
-        changes to the model such as changed bounds, added or removed
-        reactions are not taken into account by this method.
+def metabolite_summary(met, solution=None, threshold=0.01, fva=False,
+                       names=False, floatfmt='.3g'):
+    """
+    Print a summary of the production and consumption fluxes.
 
-    threshold : float
-        a value below which to ignore reaction fluxes
+    This method requires the model for which this metabolite is a part
+    to be solved.
 
-    fva : float (0->1), or None
+    Parameters
+    ----------
+    solution : cobra.Solution, optional
+        A previously solved model solution to use for generating the
+        summary. If none provided (default), the summary method will
+        resolve the model. Note that the solution object must match the
+        model, i.e., changes to the model such as changed bounds,
+        added or removed reactions are not taken into account by this
+        method.
+    threshold : float, optional
+        Threshold below which fluxes are not reported.
+    fva : pandas.DataFrame, float or None, optional
         Whether or not to include flux variability analysis in the output.
-        If given, fva should be a float between 0 and 1, representing the
+        If given, fva should either be a previous FVA solution matching
+        the model or a float between 0 and 1 representing the
         fraction of the optimum objective to be searched.
-
-    floatfmt : string
-        format method for floats, passed to tabulate. Default is '.3g'.
+    names : bool, optional
+        Emit reaction and metabolite names rather than identifiers (default
+        False).
+    floatfmt : string, optional
+        Format string for floats (default '.3g').
 
     """
+    if names:
+        emit = attrgetter('name')
+    else:
+        emit = attrgetter('id')
     if solution is None:
         met.model.slim_optimize(error_value=None)
         solution = get_solution(met.model, reactions=met.reactions)
 
+    rxns = sorted(met.reactions, key=attrgetter("id"))
     rxn_id = list()
+    rxn_name = list()
     flux = list()
     reaction = list()
-    for rxn in met.reactions:
-        rxn_id.append(format_long_string(rxn.id, 10))
-        flux.append(solution.fluxes[rxn.id] * rxn.metabolites[met])
-        reaction.append(format_long_string(rxn.reaction, 40 if fva else 50))
-
-    flux_summary = pd.DataFrame(data={
-        "id": rxn_id, "flux": flux, "reaction": reaction})
-
-    if fva:
+    for rxn in rxns:
+        rxn_id.append(rxn.id)
+        rxn_name.append(format_long_string(emit(rxn), 10))
+        flux.append(solution[rxn.id] * rxn.metabolites[met])
+        txt = rxn.build_reaction_string(use_metabolite_names=names)
+        reaction.append(format_long_string(txt, 40 if fva is not None else 50))
+
+    flux_summary = pd.DataFrame({
+        "id": rxn_name,
+        "flux": flux,
+        "reaction": reaction
+    }, index=rxn_id)
+
+    if fva is not None:
+        if hasattr(fva, 'columns'):
+            fva_results = fva
+        else:
             fva_results = flux_variability_analysis(
-            met.model, met.reactions, fraction_of_optimum=fva)
-
-        flux_summary.index = flux_summary["id"]
-        flux_summary["maximum"] = zeros(len(rxn_id))
-        flux_summary["minimum"] = zeros(len(rxn_id))
-        for rid, rxn in zip(rxn_id, met.reactions):
-            imax = rxn.metabolites[met] * fva_results.loc[rxn.id, "maximum"]
-            imin = rxn.metabolites[met] * fva_results.loc[rxn.id, "minimum"]
-            flux_summary.loc[rid, "fmax"] = (imax if abs(imin) <= abs(imax)
-                                             else imin)
-            flux_summary.loc[rid, "fmin"] = (imin if abs(imin) <= abs(imax)
-                                             else imax)
+                met.model, list(met.reactions), fraction_of_optimum=fva)
+
+        flux_summary["maximum"] = zeros(len(rxn_id), dtype=float)
+        flux_summary["minimum"] = zeros(len(rxn_id), dtype=float)
+        for rxn in rxns:
+            fmax = rxn.metabolites[met] * fva_results.at[rxn.id, "maximum"]
+            fmin = rxn.metabolites[met] * fva_results.at[rxn.id, "minimum"]
+            if abs(fmin) <= abs(fmax):
+                flux_summary.at[rxn.id, "fmax"] = fmax
+                flux_summary.at[rxn.id, "fmin"] = fmin
+            else:
+                # Reverse fluxes.
+                flux_summary.at[rxn.id, "fmax"] = fmin
+                flux_summary.at[rxn.id, "fmin"] = fmax
 
-    assert flux_summary.flux.sum() < 1E-6, "Error in flux balance"
+    assert flux_summary["flux"].sum() < 1E-6, "Error in flux balance"
 
     flux_summary = _process_flux_dataframe(flux_summary, fva, threshold,
                                            floatfmt)
 
     flux_summary['percent'] = 0
-    total_flux = flux_summary[flux_summary.is_input].flux.sum()
+    total_flux = flux_summary.loc[flux_summary.is_input, "flux"].sum()
 
     flux_summary.loc[flux_summary.is_input, 'percent'] = \
         flux_summary.loc[flux_summary.is_input, 'flux'] / total_flux
@@ -84,7 +112,7 @@ def metabolite_summary(met, solution=None, threshold=0.01, fva=False,
     flux_summary['percent'] = flux_summary.percent.apply(
         lambda x: '{:.0%}'.format(x))
 
-    if fva:
+    if fva is not None:
         flux_table = tabulate(
             flux_summary.loc[:, ['percent', 'flux', 'fva_fmt', 'id',
                                  'reaction']].values, floatfmt=floatfmt,
@@ -115,28 +143,38 @@ def metabolite_summary(met, solution=None, threshold=0.01, fva=False,
         pd.np.array(flux_table[2:])[~flux_summary.is_input.values]))
 
 
-def model_summary(model, solution=None, threshold=1E-8, fva=None,
+def model_summary(model, solution=None, threshold=0.01, fva=None, names=False,
                   floatfmt='.3g'):
-    """Print a summary of the input and output fluxes of the model.
+    """
+    Print a summary of the input and output fluxes of the model.
 
-    solution : cobra.core.Solution
+    Parameters
+    ----------
+    solution: cobra.Solution, optional
         A previously solved model solution to use for generating the
-        summary. If none provided (default), the summary method will resolve
-        the model. Note that the solution object must match the model, i.e.,
-        changes to the model such as changed bounds, added or removed
-        reactions are not taken into account by this method.
-
-    threshold : float
-        tolerance for determining if a flux is zero (not printed)
-
-    fva : int or None
-        Whether or not to calculate and report flux variability in the
-        output summary
-
-    floatfmt : string
-        format method for floats, passed to tabulate. Default is '.3g'.
+        summary. If none provided (default), the summary method will
+        resolve the model. Note that the solution object must match the
+        model, i.e., changes to the model such as changed bounds,
+        added or removed reactions are not taken into account by this
+        method.
+    threshold : float, optional
+        Threshold below which fluxes are not reported.
+    fva : pandas.DataFrame, float or None, optional
+        Whether or not to include flux variability analysis in the output.
+        If given, fva should either be a previous FVA solution matching
+        the model or a float between 0 and 1 representing the
+        fraction of the optimum objective to be searched.
+    names : bool, optional
+        Emit reaction and metabolite names rather than identifiers (default
+        False).
+    floatfmt : string, optional
+        Format string for floats (default '.3g').
 
     """
+    if names:
+        emit = attrgetter('name')
+    else:
+        emit = attrgetter('id')
     objective_reactions = linear_reaction_coefficients(model)
     boundary_reactions = model.exchanges
     summary_rxns = set(objective_reactions.keys()).union(boundary_reactions)
@@ -146,58 +184,70 @@ def model_summary(model, solution=None, threshold=1E-8, fva=None,
         solution = get_solution(model, reactions=summary_rxns)
 
     # Create a dataframe of objective fluxes
-    obj_fluxes = pd.DataFrame({key: solution.fluxes[key.id] * value for key,
+    obj_fluxes = pd.DataFrame({key: solution[key.id] * value for key,
                                value in iteritems(objective_reactions)},
                               index=['flux']).T
     obj_fluxes['id'] = obj_fluxes.apply(
         lambda x: format_long_string(x.name.id, 15), 1)
 
     # Build a dictionary of metabolite production from the boundary reactions
-    metabolite_fluxes = {}
+    metabolites = {m for r in boundary_reactions for m in r.metabolites}
+    index = sorted(metabolites, key=attrgetter('id'))
+    metabolite_fluxes = pd.DataFrame({
+        'id': [format_long_string(emit(m), 15) for m in index],
+        'flux': zeros(len(index), dtype=float)
+    }, index=[m.id for m in index])
     for rxn in boundary_reactions:
         for met, stoich in iteritems(rxn.metabolites):
-            metabolite_fluxes[met] = {
-                'id': format_long_string(met.id, 15),
-                'flux': stoich * solution.fluxes[rxn.id]}
+            metabolite_fluxes.at[met.id, 'flux'] += stoich * solution[rxn.id]
 
     # Calculate FVA results if requested
-    if fva:
+    if fva is not None:
+        if len(index) != len(boundary_reactions):
+            LOGGER.warning(
+                "There exists more than one boundary reaction per metabolite. "
+                "Please be careful when evaluating flux ranges.")
+        metabolite_fluxes['fmin'] = zeros(len(index), dtype=float)
+        metabolite_fluxes['fmax'] = zeros(len(index), dtype=float)
+        if hasattr(fva, 'columns'):
+            fva_results = fva
+        else:
             fva_results = flux_variability_analysis(
-            model, reaction_list=boundary_reactions, fraction_of_optimum=fva)
+                model, reaction_list=boundary_reactions,
+                fraction_of_optimum=fva)
 
         for rxn in boundary_reactions:
             for met, stoich in iteritems(rxn.metabolites):
-                imin = stoich * fva_results.loc[rxn.id]['minimum']
-                imax = stoich * fva_results.loc[rxn.id]['maximum']
+                fmin = stoich * fva_results.at[rxn.id, 'minimum']
+                fmax = stoich * fva_results.at[rxn.id, 'maximum']
                 # Correct 'max' and 'min' for negative values
-                metabolite_fluxes[met].update({
-                    'fmin': imin if abs(imin) <= abs(imax) else imax,
-                    'fmax': imax if abs(imin) <= abs(imax) else imin,
-                })
+                if abs(fmin) <= abs(fmax):
+                    metabolite_fluxes.at[met.id, 'fmin'] += fmin
+                    metabolite_fluxes.at[met.id, 'fmax'] += fmax
+                else:
+                    metabolite_fluxes.at[met.id, 'fmin'] += fmax
+                    metabolite_fluxes.at[met.id, 'fmax'] += fmin
 
     # Generate a dataframe of boundary fluxes
-    metabolite_fluxes = pd.DataFrame(metabolite_fluxes).T
     metabolite_fluxes = _process_flux_dataframe(
         metabolite_fluxes, fva, threshold, floatfmt)
 
     # Begin building string output table
     def get_str_table(species_df, fva=False):
         """Formats a string table for each column"""
-
-        if not fva:
-            return tabulate(species_df.loc[:, ['id', 'flux']].values,
-                            floatfmt=floatfmt, tablefmt='plain').split('\n')
-
-        else:
+        if fva:
             return tabulate(
                 species_df.loc[:, ['id', 'flux', 'fva_fmt']].values,
                 floatfmt=floatfmt, tablefmt='simple',
                 headers=['id', 'Flux', 'Range']).split('\n')
+        else:
+            return tabulate(species_df.loc[:, ['id', 'flux']].values,
+                            floatfmt=floatfmt, tablefmt='plain').split('\n')
 
     in_table = get_str_table(
-        metabolite_fluxes[metabolite_fluxes.is_input], fva=fva)
+        metabolite_fluxes[metabolite_fluxes['is_input']], fva=fva is not None)
     out_table = get_str_table(
-        metabolite_fluxes[~metabolite_fluxes.is_input], fva=fva)
+        metabolite_fluxes[~metabolite_fluxes['is_input']], fva=fva is not None)
     obj_table = get_str_table(obj_fluxes, fva=False)
 
     # Print nested output table
@@ -210,23 +260,27 @@ def _process_flux_dataframe(flux_dataframe, fva, threshold, floatfmt):
     """Some common methods for processing a database of flux information into
     print-ready formats. Used in both model_summary and metabolite_summary. """
 
+    abs_flux = flux_dataframe['flux'].abs()
+    flux_threshold = threshold * abs_flux.max()
+
     # Drop unused boundary fluxes
-    if not fva:
-        flux_dataframe = flux_dataframe[
-            flux_dataframe.flux.abs() > threshold].copy()
+    if fva is None:
+        flux_dataframe = flux_dataframe.loc[
+            abs_flux >= flux_threshold, :].copy()
     else:
-        flux_dataframe = flux_dataframe[
-            (flux_dataframe.flux.abs() > threshold) |
-            (flux_dataframe.fmin.abs() > threshold) |
-            (flux_dataframe.fmax.abs() > threshold)].copy()
+        flux_dataframe = flux_dataframe.loc[
+            (abs_flux >= flux_threshold) |
+            (flux_dataframe['fmin'].abs() >= flux_threshold) |
+            (flux_dataframe['fmax'].abs() >= flux_threshold), :].copy()
 
-        flux_dataframe.loc[flux_dataframe.flux.abs() < threshold, 'flux'] = 0
+        # Why set to zero? If included show true value?
+        # flux_dataframe.loc[
+        #     flux_dataframe['flux'].abs() < flux_threshold, 'flux'] = 0
 
     # Make all fluxes positive
-    if not fva:
-        flux_dataframe['is_input'] = flux_dataframe.flux >= 0
-        flux_dataframe.flux = \
-            flux_dataframe.flux.abs().astype('float')
+    if fva is None:
+        flux_dataframe['is_input'] = (flux_dataframe['flux'] >= 0)
+        flux_dataframe['flux'] = flux_dataframe['flux'].abs()
     else:
 
         def get_direction(flux, fmin, fmax):
@@ -258,7 +312,7 @@ def _process_flux_dataframe(flux_dataframe, fva, threshold, floatfmt):
             flux_dataframe.loc[:, ['flux', 'fmin', 'fmax']].applymap(
                 lambda x: x if abs(x) > 1E-6 else 0)
 
-    if fva:
+    if fva is not None:
         flux_dataframe['fva_fmt'] = flux_dataframe.apply(
             lambda x: ("[{0.fmin:" + floatfmt + "}, {0.fmax:" +
                        floatfmt + "}]").format(x), 1)