Merge pull request #14 from biosustain/FVA_infeasibility

Fva infeasibility
biosustain · Oct 24, 2014 · c65ea31 · c65ea31
2 parents 4955212 + 53c0364
commit c65ea31
Show file tree

Hide file tree

Showing 2 changed files with 26 additions and 5 deletions.
diff --git a/cameo/flux_analysis/analysis.py b/cameo/flux_analysis/analysis.py
@@ -25,6 +25,10 @@
 from cameo.flux_analysis.simulation import _cycle_free_flux
 import pandas
 
+import logging
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
 
 def flux_variability_analysis(model, reactions=None, fraction_of_optimum=0., remove_cycles=True, view=None):
     """Flux variability analysis.
@@ -57,7 +61,7 @@ def flux_variability_analysis(model, reactions=None, fraction_of_optimum=0., rem
             try:
                 obj_val = model.solve().f
             except SolveError as e:
-                print "flux_variability_analyis was not able to determine an optimal solution for objective %s" % model.objective
+                logger.debug("flux_variability_analyis was not able to determine an optimal solution for objective %s" % model.objective)
                 raise e
             if model.objective.direction == 'max':
                 fix_obj_constraint = model.solver.interface.Constraint(model.objective.expression,
@@ -153,6 +157,7 @@ def _flux_variability_analysis(model, reactions=None):
     else:
         reactions = model._ids_to_reactions(reactions)
     fva_sol = OrderedDict()
+    [lb_flag, ub_flag] = [False, False]
     for reaction in reactions:
         fva_sol[reaction.id] = dict()
         model.objective = reaction
@@ -163,21 +168,36 @@ def _flux_variability_analysis(model, reactions=None):
         except Unbounded:
             fva_sol[reaction.id]['lower_bound'] = -numpy.inf
         except Infeasible:
-            fva_sol[reaction.id]['lower_bound'] = 0
-    for reaction in reactions:
-        model.objective = reaction
+            lb_flag = True
+
         model.objective.direction = 'max'
         try:
             solution = model.solve()
             fva_sol[reaction.id]['upper_bound'] = solution.f
         except Unbounded:
             fva_sol[reaction.id]['upper_bound'] = numpy.inf
         except Infeasible:
+            ub_flag = True
+
+        if lb_flag is True and ub_flag is True:
+            fva_sol[reaction.id]['lower_bound'] = 0
             fva_sol[reaction.id]['upper_bound'] = 0
+            [lb_flag, ub_flag] = [False, False]
+        elif lb_flag is True and ub_flag is False:
+            fva_sol[reaction.id]['lower_bound'] = fva_sol[reaction.id]['upper_bound']
+            lb_flag = False
+        elif lb_flag is False and ub_flag is True:
+            fva_sol[reaction.id]['upper_bound'] = fva_sol[reaction.id]['lower_bound']
+            ub_flag = False
+
+
     model.objective = original_objective
     df = pandas.DataFrame.from_dict(fva_sol, orient='index')
     lb_higher_ub = df[df.lower_bound > df.upper_bound]
-    assert ((lb_higher_ub.lower_bound - lb_higher_ub.upper_bound) < 1e-6).all()  # Assert that these cases really only numerical artifacts
+    try: # this is an alternative solution to what I did above with flags
+        assert ((lb_higher_ub.lower_bound - lb_higher_ub.upper_bound) < 1e-6).all()  # Assert that these cases really only numerical artifacts
+    except AssertionError as e:
+        logger.debug(zip(model.reactions, (lb_higher_ub.lower_bound - lb_higher_ub.upper_bound) < 1e-6))
     df.lower_bound[lb_higher_ub.index] = df.upper_bound[lb_higher_ub.index]
     return df
 

diff --git a/cameo/solver_based_model.py b/cameo/solver_based_model.py
@@ -29,6 +29,7 @@
 
 from cobra.core import Solution
 from cobra.core.Reaction import Reaction as OriginalReaction
+from cobra.core.Metabolite import Metabolite
 from cobra.core.Model import Model
 from cobra.core.DictList import DictList
 from cobra.manipulation.delete import find_gene_knockout_reactions