Merge pull request #583 from opencobra/refactor-bounds

refactor: change behavior of `find_bounds`

HISTORY.rst

@@ -3,6 +3,8 @@ History
 
 Next Release
 ------------
+* Instead of min/max bounds consider the median bounds for testing (un-)bounded 
+  fluxes.
 * Use a model context for every test case.
 * Fix bug which involved find_transport_reactions to ignore compartments.
 * Internal change to use model context rather than copy.

memote/suite/tests/test_consistency.py

@@ -322,7 +322,6 @@ def test_find_reactions_unbounded_flux_default_condition(model):
     missing cofactors, incorrectly defined transport reactions and more.
 
     """
-    # TODO: Arbitrary threshold right now! Update after meta study!
     ann = test_find_reactions_unbounded_flux_default_condition.annotation
     unbounded_rxn_ids, fraction, _ = \
         consistency.find_reactions_with_unbounded_flux_default_condition(model)
@@ -333,7 +332,8 @@ def test_find_reactions_unbounded_flux_default_condition(model):
         reactions) can carry unbounded flux in the default model
         condition. Unbounded reactions may be involved in
         thermodynamically infeasible cycles: {}""".format(
-            len(ann["data"]), ann["metric"], truncate(ann["data"])
+            ann["metric"], len(ann["data"]), truncate(ann["data"])
         )
     )
+    # TODO: Arbitrary threshold right now! Update after meta study!
     assert ann["metric"] <= 0.1, ann["message"]

memote/support/consistency.py

@@ -570,13 +570,14 @@ def find_reactions_with_unbounded_flux_default_condition(model):
     conditionally_blocked = fva_result.loc[
         fva_result.abs().max(axis=1) < TOLERANCE_THRESHOLD
     ].index.tolist()
-    max_bound = max([rxn.upper_bound for rxn in model.reactions])
-    min_bound = min([rxn.lower_bound for rxn in model.reactions])
-    # Find those reactions whose flux is close to the maximal upper or lower
-    # bound, i.e., appears unconstrained.
+    small, large = helpers.find_bounds(model)
+    # Find those reactions whose flux is close to or outside of the median
+    # upper or lower bound, i.e., appears unconstrained.
     unlimited_flux = fva_result.loc[
-       np.isclose(fva_result["maximum"], max_bound, atol=TOLERANCE_THRESHOLD) |
-       np.isclose(fva_result["minimum"], min_bound, atol=TOLERANCE_THRESHOLD)
+        np.isclose(fva_result["maximum"], large, atol=TOLERANCE_THRESHOLD) |
+        (fva_result["maximum"] > large) |
+        np.isclose(fva_result["minimum"], small, atol=TOLERANCE_THRESHOLD) |
+        (fva_result["minimum"] < small)
     ].index.tolist()
     try:
         fraction = len(unlimited_flux) / \

memote/support/helpers.py

@@ -771,28 +771,30 @@ def find_objective_function(model):
 @lrudecorator(size=2)
 def find_bounds(model):
     """
-    Return the maximal and minimal bounds of the reactions of the model.
+    Return the median upper and lower bound of the metabolic model.
 
     Bounds can vary from model to model. Cobrapy defaults to (-1000, 1000) but
     this may not be the case for merged or autogenerated models. In these
     cases, this function is used to iterate over all the bounds of all the
-    reactions and find the most extreme bound values in the model, which are
-    then used as its default max and min bound values.
+    reactions and find the median bound values in the model, which are
+    then used as the 'most common' bounds.
 
     Parameters
     ----------
     model : cobra.Model
         The metabolic model under investigation.
 
     """
-    try:
-        lower_bound = min(rxn.lower_bound for rxn in model.reactions)
-        upper_bound = max(rxn.upper_bound for rxn in model.reactions)
-    except ValueError:
-        lower_bound = -1000
-        upper_bound = 1000
-
-    if lower_bound != 0 or upper_bound != 0:
-        return lower_bound, upper_bound
-    else:
-        return -1000, 1000
+    lower_bounds = np.asarray([rxn.lower_bound for rxn in model.reactions],
+                              dtype=float)
+    upper_bounds = np.asarray([rxn.upper_bound for rxn in model.reactions],
+                              dtype=float)
+    lower_bound = np.nanmedian(lower_bounds[lower_bounds != 0.0])
+    upper_bound = np.nanmedian(upper_bounds[upper_bounds != 0.0])
+    if np.isnan(lower_bound):
+        LOGGER.warning("Could not identify a median lower bound.")
+        lower_bound = -1000.0
+    if np.isnan(upper_bound):
+        LOGGER.warning("Could not identify a median upper bound.")
+        upper_bound = 1000.0
+    return lower_bound, upper_bound