bifurcation analysis support

[isaacsas]

Discourse is suddenly overflowing with mention of packages that could be useful for steady-state / bifurcation analysis. We should look into creating interfaces to use them.

PseudoArcLengthContinuation.jl Bifurcations.jl HomotopyContinuation.jl

I guess @TorkelE has already started on this to some degree with BRNEquilibrate?

[TorkelE]

Yes, this is very useful stuff we should support. Currently BRNEquilibrate works quite well on my machine and I am having use of it, however there are still a few things I want to adjust before I write a proper documentation. (It is possible one could just integrate it all into DiffEqBiological, dpeending on how much stuff we want to put in here)

A large part of it is just to link up the reaction networks to Homotopy Continuation, and provide some good functionality for what typically is interesting when studying these kinds of system. Hopefully it should all be ready quite soon.

[isaacsas]

I'd argue that this is the last major functionality that is needed to offer many biology modelers a complete package for their day to day use. Having

1. Parameterization methods
2. ODE solvers
3. SS solvers with bifurcation analysis
4. SDE solvers
5. SSAs

covers the full gamut of what I think many modelers need on a day to day basis (and goes beyond any other systems biology modeling package that I'm aware of). At that point I'd say the package should be ready to start advertising more widely (i.e. maybe put together a paper introducing it and showing off the functionality for the biological community).

[TorkelE]

Yes, I will get together all the basic functionality of the stability/bifurcation analysis, then I will put it all up and we help finishing it.

[ChrisRackauckas]

You guys have definitely made this amazing and I agree that a publication is in order.

What we are missing is more tau leaping type of stuff. I'm not totally satisfied by the design there quite yet, but don't wait on me since that'll take a bit (we could have a discussion on this).

[isaacsas]

Yes, I totally agree about tau-leaping. Completely forgot about it.

Ultimately it would be great to reach a point where we can resolve different reactions (or species?) among ODEs, SDEs, tau-leaping and jumps with some algorithms for dynamic shuffling of species/reactions between scales. There is certainly prior work on this, some implemented in Bionetgen, but there is still lots of research to do and little that has been incorporated into existing software.

I'd be happy to chat on another issue about what you were thinking for the tau-leaping interface / what methods you'd like to get implemented. There's interesting issues there too with how to ensure non-negativity of solutions (from the Anderson approach, which is accurate and guarantees positiveness but is expensive, to more simple strategies).

[ChrisRackauckas]

Well right now it's currently implemented with the RegularJumps.

http://docs.juliadiffeq.org/latest/tutorials/discrete_stochastic_example.html#RegularJumps-and-Tau-Leaping-1 http://docs.juliadiffeq.org/latest/types/jump_types.html#Defining-a-Regular-Jump-1 http://docs.juliadiffeq.org/latest/solvers/jump_solve.html#RegularJump-Compatible-Methods-1

But they are weirdly siloed from the rest. The reason was because they want to only be about changing values, and changing a lot of them at once because all of them all calculated in the interval instead of individually. That's the core reason why they are so different and optimized so differently. If it's just building a du array each time, why call 1000 functions?

But then we end up in this situation where tau-leaping exists, but it's not easy to swap between the formulations. Maybe that's not an issue if DiffEqBiological is building the jumps anyways. But this issue is why I haven't built them out much. Otherwise, it would be pretty simple to start adding Binomial and implicit Tau leaping strategies, and I am interested in researching some optimizations for Poisson Runge-Kutta methods (willing to collaborate!). But somehow they shouldn't be so buried.

Oh, and the other thing is that our jumps should all have marking distributions associated with them. It's a common thing and has come up on the Discourse: https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453

[TorkelE]

Just a quick note on the stability/bifurcation analysis.

I earlier split it of from the DiffEqBiological package (with the idea that people not interested in it should not have to load all the required stuff). Then we might merge it later if we want.

However there are a few things that would be good to calculate on in the initial macro, and this is still quite useful, so I thought I'd might as well integrate stuff directly. However, here I got into some weirdness I do not understand. The main tool I use for finding steady states is HomotopyContinuation.jl. When I write this into julia:

@time using DiffEqBiological
@time using HomotopyContinuation
@time using DynamicPolynomials

I receive quite short times (~10, ~3s, ~1s). However, if I try to put either using HomotopyContinuation or using DynamicPolynomials inside DiffEqBiological.jl I run into trouble. E.g. if I have:

__precompile__()

module DiffEqBiological

using Reexport
using DiffEqBase
using DiffEqJump
using SymEngine
using MacroTools
using DataStructures
using Parameters

using HomotopyContinuation

@reexport using DiffEqBase, DiffEqJump
...

Suddenly the first time I run

@time using DiffEqBiological

I get ~250s! (first time I run it, 2nd is still fast. But 1st time should be better than that). Intuitively I thought it maybe it should take ~13s, but definatly not this (the same goes if I add only using DynamicPolynomials). I tried running it, stop Julia, running it again, but still ~250s.

Anyone know what might be going on, and how to fix it?

[isaacsas]

That seems really weird. As an alternative, what info do you want to calculate or need from the macro that you can’t currently get at? We could always save more of the generated expressions so you can work with them in a separate package.

[ChrisRackauckas]

Precompilation might be going on a wild goose chase and doing a bunch of unnecessary combinations. I'll get someone who knows this stuff to comment.

[JeffBezanson]

Could you try enabling #define ENABLE_TIMINGS in src/options.h and running the thing that takes 250 sec? That will give us a bit more info.

[TorkelE]

Where do I find the src/options.h file? I have tried to simply add such a file to the DiffEqBiological src folder, adding a single line define ENABLE_TIMINGS (no effect). I have found such a file (with such an entry) here https://github.com/JuliaLang/julia/blob/master/src/options.h but I have not been able to find the corresponding file on my computer (digging through the .julia folder).

[JeffBezanson]

Yes, it's in the julia source tree. You have to rebuild julia for it to take effect.

[ChrisRackauckas]

We can have @YingboMa give it a shot.

[TorkelE]

Is there a good instruction manual for how to rebuild Julia (which does not assume to much competence from the user)?

My (failed) approach:

• Create a new folder (under home) called "my_julia".
• Make it a git folder by moving into it and writing git init in the terminal.
• Fetch the latest julia version using git clone git://github.com/JuliaLang/julia.git
• Go into the src folder, find the options.h file and change the line // #define ENABLE_TIMINGS to define ENABLE_TIMINGS.
• Go into the julia folder and write make in the terminal.
• Wait for quite a while until it finishes.
• Now I tried writing ./julia in various folder, but always get the response bash: ./julia: No such file or directory. If I restart juno it still says that I am on the same julia version as I had before (indicating it does not take the update into account). Running @time using DiffEqBiological yields no additional input.

I probably did something wrong somewhere, but this is also quite far out of my comfort zone, so I am a bit unsure on what is actually happening.

[YingboMa]

diff --git a/REQUIRE b/REQUIRE
index ff7a248..9b2109c 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -7,3 +7,5 @@ Reexport 0.1.0
 SymEngine 0.4.0
 MacroTools 0.4.0
 Parameters 0.10.3
+HomotopyContinuation
+DynamicPolynomials
diff --git a/src/DiffEqBiological.jl b/src/DiffEqBiological.jl
index 3e2b6d2..bb08ceb 100644
--- a/src/DiffEqBiological.jl
+++ b/src/DiffEqBiological.jl
@@ -9,6 +9,8 @@ using SymEngine
 using MacroTools
 using DataStructures
 using Parameters
+using HomotopyContinuation
+using DynamicPolynomials
 @reexport using DiffEqBase, DiffEqJump

 using Compat

ggnore@banana [06:00:16 PM] [~/.julia/dev/DiffEqBiological] [master *]
-> % rm -rf ~/.julia/compiled/v1.2/DiffEqBiological
ggnore@banana [06:00:21 PM] [~/.julia/dev/DiffEqBiological] [master *]
-> % j --startup-file=no -e '@time using DiffEqBiological'
ROOT                      : 2.18 %   294271331
GC                        : 5.35 %   722614381
LOWERING                  : 1.82 %   245775799
PARSING                   : 0.80 %   107689567
INFERENCE                 : 10.97 %   1481990989
CODEGEN                   : 0.67 %   90757761
METHOD_LOOKUP_SLOW        : 0.04 %   5455198
METHOD_LOOKUP_FAST        : 1.67 %   225919429
LLVM_OPT                  : 5.85 %   789692388
LLVM_MODULE_FINISH        : 0.31 %   41640641
LLVM_EMIT                 : 0.01 %   1281748
METHOD_LOOKUP_COMPILE     : 0.02 %   3066380
METHOD_MATCH              : 18.49 %   2497016508
TYPE_CACHE_LOOKUP         : 1.89 %   255215760
TYPE_CACHE_INSERT         : 0.50 %   68063377
STAGED_FUNCTION           : 0.00 %   515479
MACRO_INVOCATION          : 0.01 %   1260518
AST_COMPRESS              : 0.27 %   35939730
AST_UNCOMPRESS            : 0.57 %   76661775
SYSIMG_LOAD               : 1.78 %   240024280
ADD_METHOD                : 38.66 %   5221285920
LOAD_MODULE               : 7.06 %   952887030
SAVE_MODULE               : 1.00 %   135440285
INIT_MODULE               : 0.08 %   11040573
  9.414693 seconds (8.71 M allocations: 483.035 MiB, 2.05% gc time)
ROOT                      : 52.42 %   13696049747
GC                        : 2.07 %   541188975
LOWERING                  : 0.06 %   16657607
PARSING                   : 0.00 %   293871
INFERENCE                 : 3.12 %   813895693
CODEGEN                   : 0.50 %   131298300
METHOD_LOOKUP_SLOW        : 0.01 %   3745914
METHOD_LOOKUP_FAST        : 0.46 %   119120031
LLVM_OPT                  : 11.43 %   2986597813
LLVM_MODULE_FINISH        : 0.17 %   44925647
LLVM_EMIT                 : 0.01 %   1590317
METHOD_LOOKUP_COMPILE     : 0.01 %   2012762
METHOD_MATCH              : 5.87 %   1533651983
TYPE_CACHE_LOOKUP         : 0.51 %   134082358
TYPE_CACHE_INSERT         : 0.22 %   57587979
STAGED_FUNCTION           : 0.00 %   413925
MACRO_INVOCATION          : 0.00 %   33132
AST_COMPRESS              : 0.08 %   21848957
AST_UNCOMPRESS            : 0.17 %   45139219
SYSIMG_LOAD               : 0.94 %   244619858
ADD_METHOD                : 18.21 %   4758477430
LOAD_MODULE               : 3.69 %   963029209
INIT_MODULE               : 0.04 %   10956227

I couldn't recreate the long precompilation issue on Julia master.

[TorkelE]

So basically this was ready, was hoping to get it up before the weekend. But then git got in-between. Haven't been able to figure it out yet and now I've gotten to sleep to continue. Made a thread about it on StackOverflow: https://stackoverflow.com/questions/54837083/git-problem-error-src-refspec-origin-does-not-match-any

[korsbo]

This is mostly done, See #97. What's currently lacking is documentation.

[isaacsas]

From @TorkelE on Slack, to recover the steady-state support we had we can do

using ModelingToolkit
import HomotopyContinuation
const HC = HomotopyContinuation

# Define a MT nonlinear system

@variables x y z
@parameters σ ρ β

eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])


# convert variables and parametrs from nonlinear system to HC vars

hc_vars = HC.Variable.(Symbol.(ns.states))
hc_params = HC.Variable.(Symbol.(ns.ps))

# Now we just would need to evaluate the nonlinear system at (hc_vars, hc_params)
# However this seems to be possibel without actually compiling a function (sigh...)

nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β], expression=Val{false})[1]

hc_eqs = nlsys_func(hc_vars, hc_params)
HC.System(hc_eqs; variables = hc_vars, parameters = hc_params)

[TorkelE]

I'd say this is solver: https://docs.sciml.ai/Catalyst/dev/tutorials/bifurcation_diagram/