cobrapy
cobrapy copied to clipboard
Error in loopless FVA
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
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.
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: