cobrapy icon indicating copy to clipboard operation
cobrapy copied to clipboard

Error in loopless FVA

Open Midnighter opened this issue 5 years ago • 2 comments

Problem description

It seems like the loopless FVA code is not robust against infeasible solutions when it comes to the Gurobi solver. It does not happen with the GLPK solver although many solutions are infeasible or the status is undefined. CPLEX seems also fine.

Code Sample

Download the yeast-GEM 8.3.3.

import cobra
from cobra.io import read_sbml_model
from cobra.flux_analysis import flux_variability_analysis


config = cobra.Configuration()
config.solver = "gurobi"
model = read_sbml_model("yeastGEM.xml")
fva_result = flux_variability_analysis(model, loopless=True)

Actual Output

cobra/util/solver.py:416 UserWarning: solver status is 'infeasible'

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<timed exec> in <module>

~/.virtualenvs/carrot/lib/python3.6/site-packages/cobra/flux_analysis/variability.py in flux_variability_analysis(model, reaction_list, loopless, fraction_of_optimum, pfba_factor, processes)
    195             else:
    196                 _init_worker(model, loopless, what[:3])
--> 197                 for rxn_id, value in map(_fva_step, reaction_ids):
    198                     fva_result.at[rxn_id, what] = value
    199 

~/.virtualenvs/carrot/lib/python3.6/site-packages/cobra/flux_analysis/variability.py in _fva_step(reaction_id)
     47     sutil.check_solver_status(_model.solver.status)
     48     if _loopless:
---> 49         value = loopless_fva_iter(_model, rxn)
     50     else:
     51         value = _model.solver.objective.value

~/.virtualenvs/carrot/lib/python3.6/site-packages/cobra/flux_analysis/loopless.py in loopless_fva_iter(model, reaction, solution, zero_cutoff)
    221         # If the previous optimum is maintained in the loopless solution it was
    222         # loopless and we are done
--> 223         if abs(reaction.flux - current) < zero_cutoff:
    224             if solution:
    225                 return sol

TypeError: unsupported operand type(s) for -: 'float' and 'NoneType'

Expected Output

I expect to receive a valid FVA result.

Dependency Information

System Information

OS Linux OS-release 4.15.0-52-generic Python 3.6.7

Package Versions

cobra 0.15.3 depinfo 1.5.1 future 0.17.1 numpy 1.16.4 optlang 1.4.4 pandas 0.24.2 pbr 5.1.3 pip 19.1.1 python-libsbml-experimental 5.18.0 ruamel.yaml 0.15.97 setuptools 41.0.1 six 1.12.0 swiglpk 4.65.0 tabulate 0.8.3 wheel 0.33.4

Midnighter avatar Jun 27 '19 11:06 Midnighter

Actually, it's more disconcerting to me that the other solvers don't raise an error. Loopless FVA will definitely not work when your model is infeasible and any result you would get from that is probably meaningless. Also getting infeasible solutions if your initial model was feasible is a bug per se but may happen due to numerical instabilities. The original solution is always a feasible solutions on all optimizations in the iteration so there should be no infeasible solutions unless your initial objective optimization is infeasible.

I think the behavioral difference is due to whether solvers invalidate the primal values upon non-optimal solution states and that behavior is managed completely in optlang.

I think the correct handling for that would be to mark those unstable reactions with NaN. If you agree I can start a PR.

cdiener avatar Jun 27 '19 18:06 cdiener

I suppose issuing a warning and maintaining the loop is the safest option?

If you agree I can start a PR.

Please go ahead :slightly_smiling_face:

Midnighter avatar Jun 28 '19 08:06 Midnighter