Loopless solutions for optlang [fix #385] #391

cdiener · 2017-02-14T22:51:13Z

Up to now only implements the original method from Schellenberger. First time the optlang interface really shines. I could basically implement the original method as is. This requires less constraints than the prior version in cobrapy and brings a speed up (twice as fast).

Missing:

a posteriori version (should be much faster)
faster loopless FVA using the method from CycleFreeFlux (semi a posteriori)

As expected a posteriori methods are faster. For FVA the loopless version is still slower than the non-loopless one for large problems (kind of expected). However it is much faster than using add_loopless or construct_loopless_model followed by calling "normal" FVA. Faster meaning 2 seconds vs 2 minutes for textbook and a much larger gap (~1000x) for ecoli.

API

This one I am still not sure about. So fire away with your comments

A priori:

with model:
    add_loopless(model)
    model.optimize()
    single_gene_deletion(model)

A posteriori:

with model:
    model.optimize()
    fluxes = loopless_solution(model)

Maybe that one should take a flux distribution as argument alongside with the model? That way you could use it to convert flux samples or saved flux distributions to loopless ones as well...

Benchmarks

test_output.txt

pstjohn · 2017-02-15T01:28:58Z

cobra/flux_analysis/loopless.py

+       Erratum in: Biophys J. 2011 Mar 2;100(5):1381.
+    """
+    internal = [i for i, r in enumerate(model.reactions) if not r.boundary]
+    Sint = model.S[:, np.array(internal)]


Do we think this is an operation we'll do many times? Just thinking that we could speed things up if I had the array-generating function accept an optional argument to only give you certain columns...

In this case though, probably not a speed-limiting factor.

For now it's only used here...

pstjohn · 2017-02-15T01:31:40Z

cobra/flux_analysis/loopless.py

+        coefs = {model.solver.variables[
+            "delta_g_" + model.reactions[ridx].id]: row[i]
+            for i, ridx in enumerate(internal) if abs(row[i]) > zero_cutoff}
+        model.solver.constraints[name].set_linear_coefficients(coefs)


is this something we'll want to context-manage eventually?

Already is due to the add_to_solver above. Had to do it like that because you can't set linear coefficients for a Constraint that is not part of an optlang.Model yet.

pstjohn · 2017-02-15T01:42:41Z

cobra/test/test_flux_analysis.py

+        def _():
+            with test_model:
+                add_loopless(test_model)
+                test_model.optimize(solver="optlang-glpk")


so I'm confused: glpk can handle the MILP formulation?

Yes! Absolutely. GLPK includes a MIP solver. It does not support indicator variables like cplex or gurobi so you have to declare them by hand.

codecov-io · 2017-02-15T15:30:43Z

Codecov Report

Merging #391 into devel-2 will increase coverage by 0.16%.
The diff coverage is 79.72%.

@@             Coverage Diff             @@
##           devel-2     #391      +/-   ##
===========================================
+ Coverage    65.35%   65.52%   +0.16%     
===========================================
  Files           64       64              
  Lines         7716     7822     +106     
  Branches      1338     1359      +21     
===========================================
+ Hits          5043     5125      +82     
- Misses        2442     2459      +17     
- Partials       231      238       +7

Impacted Files	Coverage Δ
cobra/core/model.py	`85.97% <ø> (+0.66%)`	✅
cobra/flux_analysis/init.py	`50% <ø> (ø)`	✅
cobra/flux_analysis/sampling.py	`96.1% <ø> (-0.17%)`	❌
cobra/util/array.py	`100% <100%> (ø)`	✅
cobra/flux_analysis/variability.py	`92.85% <100%> (+0.31%)`	✅
cobra/util/solver.py	`80.9% <100%> (-0.67%)`	❌
cobra/flux_analysis/loopless.py	`81.25% <74.39%> (-18.75%)`	❌
cobra/test/test_flux_analysis.py	`86.81% <80%> (-0.64%)`	❌
cobra/solvers/mosek.py	`0% <ø> (ø)`	✅
cobra/core/formula.py	`18.42% <ø> (ø)`	✅
... and 4 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ac928ef...9ed85d7. Read the comment docs.

Midnighter · 2017-02-15T15:34:05Z

cobra/flux_analysis/loopless.py

+    ----------
+    .. [1] Elimination of thermodynamically infeasible loops in steady-state
+       metabolic models.
+       Schellenberger J, Lewis NE, Palsson BO.


Could this not be fixed by using u""" rather than reverting to ASCII? Internationalization is a good thing in general, we're not in the 80s anymore :)

I think in this case it's not that critical since Palsson is cited like that pretty often. If it would in actual code I would agree. But really have no preference. I can add the u""" if you would like.

You're right, he is cited with O fairly often. Let's leave it then.

hredestig

made some cosmetic changes hope it was ok
Although haven't tried it live so don't not very informed opinion but I like the API you suggest, shall we proceed with it as is? Guess just need to remove the additional benchmarking code or you want to keep it?

hredestig · 2017-02-17T12:08:24Z

cobra/flux_analysis/variability.py

@@ -143,7 +158,7 @@ def _fva_optlang(model, reaction_list, fraction):
        fva_old_objective = prob.Variable(


do you need m.optimize() above instead to save the solution? also, I reverted m.optimize() to raise an exception when doing the pfba, was not intended as a final solution there but work-in-progress, with that change I guess you don't need to raise ValueError here , think that's where the error comes from...

Not really. The thing is that everything involving constructing a Solution (even LazySolution) adds overhead. This is why I use the pure optlang interface and m.solver.objective.value below (which is set by calling m.solver.optimize already). In that instance it does not make so much difference except for small models. However, in the loop below that make quite a difference.

hredestig · 2017-02-17T12:11:33Z

cobra/flux_analysis/loopless.py

+    """Add variables and constraints to make a model thermodynamically
+    feasible. This removes flux loops. The used formulation is described
+    in [1]_.
+


perhaps could give an example of when to use add_loopless versus loopless_solution to help understand the reason for both of them?

Yeah that is a good idea. Also note that zero_cutoff becomes unnecessary if there would a way in optlang to get the tolerance for the solver.

pstjohn · 2017-02-17T17:35:33Z

So here's a random question:

what if I did:

add_loopless(model)
model.add_reaction(my_new_reaction)
model.optimize()

or even

add_moma(model)
model.add_reaction(my_new_reaction)
model.optimize()

My point is that if we want these routines to return/modify model objects, maybe we should be superclassing model, and then either providing add_reaction wrappers that handle the new constraints, or return NotImplementedErrors?

but then we might lose speed and the ability to context manage.

Either way, probably not a huge concern... Maybe a simple docstring warning would be sufficient

cdiener · 2017-02-17T19:17:25Z

@pstjohn That is a good point. For now mentioning it in the docstrings and the future docs is the easiest thing. We could add a check in add_reactions but I wouldn't know for what to check right now...

cdiener · 2017-02-17T19:18:27Z

Returning to WIP status since the changes in the Solution class broke the loopless functions.

@hredestig I had to rebase due to some problems. I tried to integrate all your suggestions and some more. Could you check whether it's fine?

cdiener · 2017-02-17T16:15:04Z

cobra/flux_analysis/variability.py

@@ -143,7 +158,7 @@ def _fva_optlang(model, reaction_list, fraction):
        fva_old_objective = prob.Variable(


Not really. The thing is that everything involving constructing a Solution (even LazySolution) adds overhead. This is why I use the pure optlang interface and m.solver.objective.value below (which is set by calling m.solver.optimize already). In that instance it does not make so much difference except for small models. However, in the loop below that make quite a difference.

cdiener · 2017-02-17T16:17:00Z

cobra/flux_analysis/loopless.py

+    """Add variables and constraints to make a model thermodynamically
+    feasible. This removes flux loops. The used formulation is described
+    in [1]_.
+


Yeah that is a good idea. Also note that zero_cutoff becomes unnecessary if there would a way in optlang to get the tolerance for the solver.

cdiener · 2017-02-17T16:20:59Z

cobra/flux_analysis/loopless.py

    """
    internal = [i for i, r in enumerate(model.reactions) if not r.boundary]
-    Sint = model.S[:, numpy.array(internal)]


Actually forgot to throw an error here if there is no numpy :O

cdiener · 2017-02-18T03:09:12Z

So here is the gist of the last commits:

Model.objective still had the @resettable decorator. That together with the new set_objective caused some bugs filling up the contexts. Taking @resettable off Model.objective uncovered some bugs with set_objective where it would first zero the objective and than read the reset value (which is surprisingly zero in this case :) ).

hredestig · 2017-02-20T08:01:36Z

Looks great now, is it still wip?

cdiener · 2017-02-20T16:27:56Z

@hredestig Still figuring out whether loopless solution should have an argument fluxes that takes a dict of fluxes.

Might be useful for flux sampling (convert sample to a loopless solution) or saved flux solutions. The problem there is how to get the objective value from a dict alone. Easy if Model.objective_coefficients gives the correct objective, but what happens if I set an objective in the optlang interface directy (non-linear, uses other variables, etc.)?

* Remove unimplemented/incomplete optlang solvers from get_solver_name * Add test_choose_solver fix from #391 * fix typo

still work in progress

cdiener · 2017-02-22T19:52:52Z

Ok, should be fine now. Added the fluxes arg as well.

Note that this also fixes two bugs I found on the way (#416 and another one in util.solver_valid_atoms). Let me know if you want those in a separate PR.

hredestig · 2017-02-24T17:00:47Z

cobra/flux_analysis/loopless.py

+                model.objective.set_linear_coefficients(
+                    {rxn.forward_variable: -1, rxn.reverse_variable: 1})
+        model.solver.objective.direction = "min"
+        model.optimize(objective_sense=None)


this raises a SolveError if status is not optimal - be in try/catch?

- implement review comment from hredestig - fix some errors in loopless_optlang when using fluxes arg - add a test for fluxes arg

cdiener · 2017-02-24T20:34:04Z

Okay rebased onto devel-2 😅 and added fixes and tests. Also implemented the review from @hredestig. For now this PR also fixes #416 and #418.

hredestig

Very nice!

cdiener added the WIP work in progress label Feb 14, 2017

hredestig added ready and removed WIP work in progress labels Feb 14, 2017

pstjohn reviewed Feb 15, 2017

View reviewed changes

Midnighter reviewed Feb 15, 2017

View reviewed changes

cdiener force-pushed the feature/loopless_optlang branch from 55fd021 to 3b2c241 Compare February 15, 2017 18:32

cdiener changed the title ~~Loopless solutions for optlang [WIP]~~ Loopless solutions for optlang Feb 15, 2017

cdiener mentioned this pull request Feb 15, 2017

refactor loopless model #385

Closed

cdiener force-pushed the feature/loopless_optlang branch from cbbe3d1 to 4d96055 Compare February 16, 2017 17:39

hredestig suggested changes Feb 17, 2017

View reviewed changes

cdiener force-pushed the feature/loopless_optlang branch from dda4101 to 0225482 Compare February 17, 2017 18:31

cdiener changed the title ~~Loopless solutions for optlang~~ Loopless solutions for optlang [WIP] Feb 17, 2017

cdiener commented Feb 18, 2017

View reviewed changes

hredestig previously approved these changes Feb 20, 2017

View reviewed changes

cdiener mentioned this pull request Feb 20, 2017

On devel-2 TestHelpers.test_choose_solver is failing #407

Closed

cdiener added a commit to cdiener/cobrapy that referenced this pull request Feb 20, 2017

Add test_choose_solver fix from opencobra#391

1439dab

hredestig pushed a commit that referenced this pull request Feb 21, 2017

Revise solver lists for optlang (#409)

c0d0205

* Remove unimplemented/incomplete optlang solvers from get_solver_name * Add test_choose_solver fix from #391 * fix typo

cdiener mentioned this pull request Feb 22, 2017

Errors with nested contexts and objectives #416

Closed

cdiener force-pushed the feature/loopless_optlang branch from 8ad4d00 to 4e63570 Compare February 22, 2017 19:50

cdiener changed the title ~~Loopless solutions for optlang [WIP]~~ Loopless solutions for optlang Feb 22, 2017

cdiener mentioned this pull request Feb 23, 2017

Local test failure TestSolverBasedModel.test_transfer_objective #418

Closed

hredestig changed the title ~~Loopless solutions for optlang~~ Loopless solutions for optlang [fix #385] Feb 24, 2017

hredestig reviewed Feb 24, 2017

View reviewed changes

cdiener added 13 commits February 24, 2017 14:03

Add original method from Schellenberger

a654005

Irony is if B. Pallsson's name breaks your cobrapy CI :)

842804b

Add CycleFreeflux methods

6787775

Handle missing numpy

c31c3a0

More numpy catches

37764c8

Fixes for loopless FVA, more tests

87bf2f9

Fix docstrings and variable names.

e44d087

started adapting for new solution class

1883a16

Fix objective reset and reoptimize bugs

2569897

Fix bugs with dict in set_objective

0a2efdd

Add fluxes arg to loopless_solution.

8454ca8

Fix bug in _valid_atoms

52e4946

Several small fixes

9ed85d7

- implement review comment from hredestig - fix some errors in loopless_optlang when using fluxes arg - add a test for fluxes arg

cdiener force-pushed the feature/loopless_optlang branch from 4e63570 to 9ed85d7 Compare February 24, 2017 20:31

hredestig approved these changes Feb 27, 2017

View reviewed changes

hredestig merged commit 6cf03f8 into opencobra:devel-2 Feb 27, 2017

hredestig removed the ready label Feb 27, 2017

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Loopless solutions for optlang [fix #385] #391

Loopless solutions for optlang [fix #385] #391

cdiener commented Feb 14, 2017 •

edited

pstjohn Feb 15, 2017

cdiener Feb 15, 2017

pstjohn Feb 15, 2017

cdiener Feb 15, 2017

pstjohn Feb 15, 2017

cdiener Feb 15, 2017

codecov-io commented Feb 15, 2017 •

edited

Midnighter Feb 15, 2017

cdiener Feb 15, 2017

Midnighter Feb 15, 2017

hredestig left a comment

hredestig Feb 17, 2017 •

edited

cdiener Feb 17, 2017

hredestig Feb 17, 2017

cdiener Feb 17, 2017

pstjohn commented Feb 17, 2017

cdiener commented Feb 17, 2017

cdiener commented Feb 17, 2017 •

edited

cdiener Feb 17, 2017

cdiener Feb 17, 2017

cdiener Feb 17, 2017

cdiener commented Feb 18, 2017

hredestig commented Feb 20, 2017

cdiener commented Feb 20, 2017

cdiener commented Feb 22, 2017

hredestig Feb 24, 2017

cdiener commented Feb 24, 2017

hredestig left a comment

		@@ -143,7 +158,7 @@ def _fva_optlang(model, reaction_list, fraction):
		fva_old_objective = prob.Variable(

Loopless solutions for optlang [fix #385] #391

Loopless solutions for optlang [fix #385] #391

Conversation

cdiener commented Feb 14, 2017 • edited

Missing:

API

Benchmarks

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov-io commented Feb 15, 2017 • edited

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hredestig left a comment

Choose a reason for hiding this comment

hredestig Feb 17, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pstjohn commented Feb 17, 2017

cdiener commented Feb 17, 2017

cdiener commented Feb 17, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cdiener commented Feb 18, 2017

hredestig commented Feb 20, 2017

cdiener commented Feb 20, 2017

cdiener commented Feb 22, 2017

Choose a reason for hiding this comment

cdiener commented Feb 24, 2017

hredestig left a comment

Choose a reason for hiding this comment

cdiener commented Feb 14, 2017 •

edited

codecov-io commented Feb 15, 2017 •

edited

hredestig Feb 17, 2017 •

edited

cdiener commented Feb 17, 2017 •

edited