Error of FVA!

[mahdie hadi]

Hi everybody ,

I have tested FVA for a model , it is a large model and it has about 7000 reactions, 

after testing FVA for  this model , I have got these errors:

???  In an assignment  A(I) = B, the number of elements in B and
I must be the same.

Error in ==> fluxVariability at 163
       maxFlux(i) = LPsolution.full(LPproblem.c~=0);

Since I have tested FVA with the same commands for some other models and it work, I guess the problem is with my model, but I can not detect any special problem in my model , could you please tell me what could be the problem ?

Best,

Mahdieh 

[Ben]

Hi Mahdieh -

I'm getting this error (or a similar error?) while I'm using the iMM904 model, and it's very puzzling because I get it inconsistently when I use the findBlockedReaction function

To troubleshoot, I added the following following line 161 in fluxVariablility.m:

        if length(LPproblem.c) ~= length(LPsolution.full)
            i
        end

When I run findBlockedReaction, which calls fluxVariability thus:
[minMax(:,1),minMax(:,2)] = fluxVariability(model,0);

(and I'll admit, I'm confused by the use of 0 for the optPercentage argument),

it breaks on different iterations, and sometimes doesn't break at all:

Checking for blocked reactions...
i =
   800
Index exceeds matrix dimensions.
Error in fluxVariability (line 166)
        maxFlux(i) = LPsolution.full(LPproblem.c~=0);
Error in testYeastModel_glu_min>findBlockedReaction_ks (line 753)
    [minMax(:,1),minMax(:,2)] = fluxVariability(model,0);

Checking for blocked reactions...
i =
   440
Index exceeds matrix dimensions.
Error in fluxVariability (line 166)
        maxFlux(i) = LPsolution.full(LPproblem.c~=0);
Error in testYeastModel_glu_min>findBlockedReaction_ks (line 753)
    [minMax(:,1),minMax(:,2)] = fluxVariability(model,0);

Checking for blocked reactions...
i =
   330
Index exceeds matrix dimensions.
Error in fluxVariability (line 166)
        maxFlux(i) = LPsolution.full(LPproblem.c~=0);
Error in testYeastModel_glu_min>findBlockedReaction_ks (line 753)
    [minMax(:,1),minMax(:,2)] = fluxVariability(model,0);

etc...

It's a puzzle!

Has anyone else seen this or made any progress on troubleshooting or fixing it?

- Ben

[Ben]

p.s. Since this is using LPsolution.full, it is probably worth mentioning that I'm using Gurobi 5.0.

-b

On Tuesday, November 13, 2012 1:37:05 PM UTC-8, M Hadi wrote:

[Ben]

This error appears to be related to Gurobi, and remains unsettling to me because I do not how to troubleshoot it further, nor the best fix.

With Gurobi, it appears that LPsolution.full (which is returned by solveCobraLP) may sometimes be a vector of length 0 (such as when the problem is infeasible). In such cases, line 163 of fluxVariability,

maxFlux(i) = LPsolution.full(LPproblem.c~=0);

would fail and give the error I've been seeing.

When I use glpk, I don't get the error (at least, not on the half dozen iterations I've tried so far). I suspect (but have not confirmed) that glpk returns a LPsolution.full of 0s in the infeasible case, unlike Gurobi. If so, such behavior would allow the fluxVariability code to continue without error.

It may be worth adding a try/catch statement around line 163 for cases where length(LPproblem.c) ~= length(LPsolution.full), but then I'm not sure how the code should respond if the test fails. Two options I've thought of are:

1) skip that FVA iteration
2) retry that FVA iteration - the intermittent nature of the error seems to suggest that Gurobi returns different LPsolution vectors when it retries the  (same?) problem again. This is of concern, but I'm not sure how to address it.
3) require use of a non-Gurobi solver for fluxVariability (but then how do we test them, and what should the solver be doing, and why is Gurobi behaving differently and seemingly inconsistently?)

What do you all think?

- Ben

[J H]

I know I was previously using an iterated FVA approach with Gurobi to create bounding planes upon solutions within an area of variation in two parameters, and when a solution was infeasible the returned vectors were full length of [N/A]. I'm not sure if I wrote that into the LPsolution in a similar way to how you've suggested, or if it was in the script assembling surfaces in response to a failed solution, but I agree that adding a try/catch or equivalent somewhere is best. 

--
 
---
You received this message because you are subscribed to the Google Groups "COBRA Toolbox" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cobra-toolbo...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[J H]

It'd also be great if you could determine if this "indeterministic" solution behaviour is actually occurring with Gurobi, or if these are simply infeasible solutions. 

[Ben]

Continuing the troubleshooting, I've found that I can pass the same LP problem to Gurobi 10,000 times, and get the result that it is infeasible in 2 cases and feasible the other 9,998 times. So it appears that Gurobi is behaving inconsistently.

I've posted to the gurobi google group to follow up there - https://groups.google.com/forum/?fromgroups=#!topic/gurobi/W9q29-eqKyw

Meanwhile, my hack will be to re-run the gurobi call in FVA if it returns infeasible the first time - that'll let me continue with my research for now. For now, it may be best to let COBRA's fluxVariablility function continue to fail with this gurobi error, rather than hiding it by including my workaround in the COBRA repository.

- Ben

[Ben]

Final followup, from the gurobi group:

"I noticed you are setting the PreDual parameter. When I remove the
PreDual parameter setting your test passes. There was a bug with
PreDual=2 that was fixed in 5.0.2 that could be the culprit. Please
upgrade to the latest version: Gurobi 5.1.0. "
 - Chris Maes

I've upgraded now, and don't get the error anymore with Gurobi 5.1.0.

- Ben

[Ronan Fleming]

Hi Ben,

Chris Maes is certainly a good source of expertise.

NaN would be the appropriate entry for each row of a vector returned
when the problem seems infeasible. I say "seems" because not all
solvers can really always detect when a problem is infeasible or not
without expert tuning of the solvers parameters, before the solve.

Ronan

--
--
Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.
-----------------------------------------------------------------
Assistant Professor,
Dept. Biochemistry & Molecular Biology,
Faculty of Medicine,
University of Iceland,
&
Center for Systems Biology,
University of Iceland,
Sturlugata 8,
101 Reykjavik,
Iceland.
http://www.hi.is/~rfleming
http://systemsbiology.hi.is/
Ph: +354 618 6245
Skype: ronan.fleming
-----------------------------------------------------------------

[Christopher Maes]

> Chris Maes is certainly a good source of expertise.

I didn't realize until his subsequent email that Ben was working on a
flux balance application. (Ronan and I collaborated on an FBA paper
several years ago.) It's nice to hear that FBA, Gurobi, and
optimization are being discussed on this list. I should probably
subscribe. :)

> NaN would be the appropriate entry for each row of a vector returned
> when the problem seems infeasible. I say "seems" because not all
> solvers can really always detect when a problem is infeasible or not
> without expert tuning of the solvers parameters, before the solve.

With the exception of bugs like the one Ben reported (and models
suffering from numerical troubles), if Gurobi returns an infeasible
status the problem should be infeasible. As Ronan noted, there are a
few parameters settings that are useful to know:

1) Setting the parameter InfUnbdInfo=1 will cause Gurobi to return a
Farkas infeasibility proof when an LP is declared infeasible. This
proof should allow you to verify the LP is infeasible.

2) Gurobi may also return the status code INF_OR_UNBD (infeasible or
unbounded) in some cases. This is because our presolve routines
perform dual reductions to reduce the size of the problem. These
reductions speed up the solve, but make it difficult to distinguish
between the infeasible and unbounded cases. To obtain a more
definitive conclusion, set the DualReductions parameter to 0 and
reoptimize. After the reoptimization you should get a infeasible
status code, if your problem is infeasible.

These parameters are describe in further detail in the Parameters
section of the Gurobi reference manual:
http://www.gurobi.com/documentation/5.1/reference-manual/node775#sec:Parameters

If you have further questions (or bug reports!), please reach out to
us on the Gurobi google group:
http://groups.google.com/group/gurobi

Best of luck with your system biology problems,
Chris