opencobra / cobrapy

COBRApy is a package for constraint-based modeling of metabolic networks.
http://opencobra.github.io/cobrapy/
GNU General Public License v2.0
465 stars 218 forks source link

Error in loopless FVA #865

Closed Midnighter closed 1 year ago

Midnighter commented 5 years ago

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
cdiener commented 5 years ago

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.

Midnighter commented 5 years ago

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: