CaoTauLeaping Implementation

[sivasathyaseeelan]

using StochasticDiffEq, JumpProcesses, DiffEqBase, Statistics, Plots
using Test, LinearAlgebra

function regular_rate(out,u,p,t)
    out[1] = (0.1/1000.0)*u[1]*u[2]
    out[2] = 0.01u[2]
end

const dc = zeros(3, 2)
dc[1,1] = -1
dc[2,1] = 1
dc[2,2] = -1
dc[3,2] = 1

function regular_c(du,u,p,t,counts,mark)
    mul!(du,dc,counts)
end

rj = RegularJump(regular_rate,regular_c,2)
jumps = JumpSet(rj)
iip_prob = DiscreteProblem([999.0,1,0],(0.0,250.0))
jump_iipprob = JumpProblem(iip_prob,Direct(),rj)

N = 100
sol = solve(EnsembleProblem(jump_iipprob),CaoTauLeaping();dt=1.0,trajectories = N)
plot(sol)

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[ChrisRackauckas]

what are these doing here?

[ChrisRackauckas]

these don't belong here

[ChrisRackauckas]

why is this needed?

[ChrisRackauckas]

This is fine for now, but we really should be storing that differently in the jumps. @isaacsas this is worth thinking about.

[isaacsas]

I think it would be worthwhile to make leaping work with individually defined jumps too. Split-step methods require that anyways.

But that wouldn’t solve this issue for general non-mass action jumps.

[ChrisRackauckas]

and we'd need two implementations of each leaping? I'm thinking the interface should maybe require more information in the jump building, under the assumption users will be using something like Catalyst rather than building most jumps by hand. To do some of this optimally we have to shift been a few different representations, where for a user this would be redundant information to the solver it's different views that have different optimal uses.

[isaacsas]

For methods that work on the full vector at once, can’t we can manually build the needed function from the list of jumps? So those methods can stick with the current interface, but then we could support methods that really need jump-by-jump rate/affect access too.

But yes, Catalyst should build the needed inputs for different representations.

[ChrisRackauckas]

For methods that work on the full vector at once, can’t we can manually build the needed function from the list of jumps?

You get a lot of compile time trying to inline all of the small functions into the large one. Also there's a lack of fusion. So it's quite slow to do it that way.

I think the point is, some times you need jump-by-jump information, and sometimes you want the aggregate information (all tau leap steppers will use an aggregate function a lot of the time), and those are the same information in two very different computational representations, and building one from the other isn't easy to do from the implementation but it is trivial to do symbolically.

[isaacsas]

Definitely, so it should be preferred that methods that work on the full vector get a single input function, but having a fallback would be useful. Catalyst/MTK, should just generate both inputs (or whichever is most appropriate for the selected integrator -- we can add a function to query the preferred input type).

[ChrisRackauckas]

The formula uses the derivative of the rate function, not the rate itself.

[ChrisRackauckas]

This is supposed to use the sum of rates term, not u[i], also why the 1.0?

[ChrisRackauckas]

Oh I see.

[ChrisRackauckas]

this is missing the g_i factor

[ChrisRackauckas]

There are 3 schemes in the paper https://people.cs.vt.edu/~ycao/publication/newstepsize.pdf, can you implement them as different algorithms (reusing most of the code, just different controller dispatches) so we can compare them?