Tau-Leaping and beyond about jumpprocesses.jl HOT 13 CLOSED

Might as well throw down my roadmap:

• Explicit tau-leap
• Implicit Euler tau-leap
• Trapezoid tau-leap
• Burrage's Poisson Runge-Kutta methods (generalized tau-leap)
• Order calculation tools for Poisson Runge-Kutta methods (for methods research, I'm really interested in building a high order tau-leap and this will be the tool)
• Binomial tau-leap (https://www.ncbi.nlm.nih.gov/pubmed/15638577)
• Other implementations of direct (our current one is a little more optimized than the Optimized Direct method of Cao Petzold, mostly because there is a bisection where they do a linear search. Should add a citation to the citation page for this https://engineering.ucsb.edu/~cse/Files/optSSA042.pdf). There's like, First Reaction method and something else.

I would like to add a bunch of tests to DiffEqBenchmarks.jl. Then hybrid methods. Method which tag some variables to tau-leap, and others to do direct on. Those can be implemented directly here without anything extra. Of course, benchmark here.

Then full hybrid methods. The design is similar, except there's 3 modes. There's tau mode, direct mode, and external mode. It's like this, if rare, do direct on this reaction. If semi-rare, do some sort of tau-leap. If not rare, then flip some boolean. Then, have this array of booleans accessible in the ODE/SDE solver function, allowing users to do something like:

if reaction_is_ode[10]
du[10] = ... # this is the ODE function
end

inside of the ODE. This can be implemented pretty simply with a smart use of pointers. Then I want to round it out with:

SciML/Catalyst.jl#2

that is, allow the reaction DSL to build both the jumps and the ODE/SDE functions, and build a DSL macro. This should give us everything that's needed to research new methods and benchmark... oh, and I guess we could actually do some biology.

Am I missing any methods? I am not an expert in this domain (yet! 👍) so I don't know if I've read all of the methods yet. My main goal here is to have a whole suite of methods to start testing all of this hybrid stuff, and have something to start building hybrid RDME+SPDE models. Like the ODE/SDE solvers, I am aiming for having a bunch of methods and offering recommendations, since this suite I would like to offer both performance and a methods research suite.

(Then the other thing is testing all of the addons for compatibility. For example, parameter estimation routines won't work here yet, but that should get fixed up this summer. We need #11 as well. @elevien already setup split coupling, which is nice! I need to find out how to have that generalize to more cases, or if that can.)

@sdwfrost you have anything to add?

from jumpprocesses.jl.

Is anyone else planning on doing this? If not, I'll probably put in a PR this weekend.

from jumpprocesses.jl.

from jumpprocesses.jl.

Adaptive tau-leaping with post leap estimates is now TauLeaping in StochasticDiffEq.jl, and it supports mixing with events and adapted jumps. Next will be a version that allows SDEs and ODEs, but that's not too far away given that it's already mixed into the SDE solver.

from jumpprocesses.jl.

@gabrielgellner is. At least he claimed it in the chat channel awhile ago. I'll let him post to see if he's still up to it.

from jumpprocesses.jl.

I have been working on a PR. But if you feel you might get it done in a weekend your are welcome to bump my efforts. If you lose interest I will continue my work anyway as it is helping understand off the DiffEq stuff works, I can file if you don't.

from jumpprocesses.jl.

Does anyone know if the Cao, Gillespie and Petzold (2006) is the most recent word on how to pick \tau for the basic leaping algorithm?

from jumpprocesses.jl.

Perfect. I want to use the simplest, but not old fashioned version first!

from jumpprocesses.jl.

Hey, @ChrisRackauckas , if you're going crazy on a roadmap, how about:

• Uniformization: slows things down, but very handy for plugging into MCMC for inference. I just added this to Gillespie.jl
• Retrodictive SSA, described here. It has funny behaviour for some models thanks to the choice of a uniform prior for previous states, but it works OK for some models.
• GPU acceleration
• Simulations with events forced at specific times (see TJ McKinley's paper)

I recall a paper that argued that the first reaction method doesn't really speed things up. I'll try and track it down, unless you're already aware.

from jumpprocesses.jl.

Thanks for the suggestions! Yup, we're going crazy here.

GPU acceleration

Does this actually help? These algorithms are essentially serial. If we get a version of the ODE solvers with events working through CUDANative.jl (I hope to do that), then we can do the Monte Carlo over the GPU. But I don't think it would necessarily help most SSA problems. Is there some form of collocation SSA that I don't know about? I don't think you can collocate a Poisson process like you can a Wiener process (Hermite), which is essential to parallelism-in-time.

Simulations with events forced at specific times (see TJ McKinley's paper)

We can already do this if the user only has events on variables which only show up in the rate equations for jumps which are VariableRateJumps. I'll have to see if the paper adds anything more.

I recall a paper that argued that the first reaction method doesn't really speed things up. I'll try and track it down, unless you're already aware.

That's the Petzold paper 2003 I linked for the optimized direct.

from jumpprocesses.jl.

I realized that we need to constrain the jumps we're using a bit to do this. Allowing a general affect! can be difficult to deal with here (if not impossible in some cases). This issue defines a new jump form which is constrained to better handle this case:

SciML/DifferentialEquations.jl#250

from jumpprocesses.jl.

A basic non-adaptive tau-leaping is implemented which shows off the new RegularJump.

https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/master/src/simple_regular_solve.jl#L4

Then there's a RegularJump SSA right below it. That show show how to do all of the pure-jump regular and hybrid methods. Adding the others is just details now. Some adaptive solvers which require pullback will be build on top of the StochasticDiffEq.jl so that way they also support jump diffusions with regular stepping.

from jumpprocesses.jl.

A lot of this has been implemented in StochasticDiffEq, with details remaining, but getting their own issues and their own progress, so I'm closing this. This issue definitely jump started the jump diffusion and leaping process so thanks to everyone involved!

from jumpprocesses.jl.