search

SciML / JumpProcesses.jl

Build and simulate jump equations like Gillespie simulations and jump diffusions with constant and state-dependent rates and mix with differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/JumpProcesses/stable/
Other
140 stars 35 forks source link

Tau-Leaping and beyond #8

Closed ChrisRackauckas closed 4 years ago

ChrisRackauckas commented 7 years ago

We should add more discrete aggeregators! Explicit (Euler) Tau-leaping is the obvious next step, with implicit Euler leaping and Trapezoid as the next steps (from one will come the others).

To implement these, one just needs to implement a new DiscreteAggregator. A new type needs to be declared:

https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/master/src/aggregators/aggregators.jl

And the logic for building the callback needs to be implemented. For example, here's the direct:

https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/master/src/aggregators/direct.jl

The DiscreteCallback for tau-leaping should work similarly. It should store a tstop for the next timepoint where it will be called. For tau-leaping, this is just the $$\Delta t$$. This could then keep the 

@inline function (p::DirectJumpAggregation)(t,u,integrator) # condition
  p.next_jump==t
end

line. But the affect! function will have to change. But that change would be to add some Poissons to the values dependent on the rates. The constructor will be very similar to the one for direct.

Without adapting the time to the next jump, a first implementation should be very short.

elevien commented 7 years ago

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

ChrisRackauckas commented 7 years ago

@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.

gabrielgellner commented 7 years ago

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.

gabrielgellner commented 7 years ago

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?

ChrisRackauckas commented 7 years ago

Yes, use the Petzold methods. Other then them, the other paper to know are the Barrage Poisson Runge Kutta methods, of which these tau leaping methods are a subset. But the general implementation can come later.

On Mar 30, 2017 1:19 PM, "Gabriel Gellner" notifications@github.com wrote:

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?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaDiffEq/DiffEqJump.jl/issues/8#issuecomment-290532702, or mute the thread https://github.com/notifications/unsubscribe-auth/ABuunoQWQRxkY2hb5rRw7ZO4eIYNKEYeks5rrA5YgaJpZM4Mqt8Y .

gabrielgellner commented 7 years ago

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

ChrisRackauckas commented 7 years ago

Might as well throw down my roadmap:

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:

https://github.com/JuliaDiffEq/DiffEqBiological.jl/issues/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 https://github.com/JuliaDiffEq/DiffEqJump.jl/issues/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?

sdwfrost commented 7 years ago

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

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.

ChrisRackauckas commented 7 years ago

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.

ChrisRackauckas commented 6 years ago

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:

https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/250

ChrisRackauckas commented 6 years ago

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.

ChrisRackauckas commented 4 years ago

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.

ChrisRackauckas commented 4 years ago

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!