SIR test gives funny answers

[sdwfrost]

Hi @ChrisRackauckas

Your version of my SIR model in test/ gives funny answers: the epidemic never seems to take off, and it doesn't give plots similar to that in the introduction, even over multiple simulations.

[sdwfrost]

In addition, the different model spec given in the docs throws an error:

using DifferentialEquations

prob = DiscreteProblem([999,1,0],(0.0,250.0))
rate = (t,u) -> (0.1/1000.0)*u[1]*u[2]
affect! = function (integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
jump = ConstantRateJump(rate,affect!)

rate = (t,u) -> 0.01*u[2]
affect! = function (integrator)
  integrator.u[2] -= 1
  integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,affect!)
jump_prob = JumpProblem(prob,jump,jump2)
sol = solve(jump_prob,Discrete(apply_map=false))

MethodError: no method matching DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}(::DiffEqBase.DiscreteProblem{Array{Int64,1},Float64,true,DiffEqBase.##130#131}, ::DiffEqJump.ConstantRateJump{##17#18,##19#20}, ::DiffEqJump.ConstantRateJump{##21#22,##23#24})
Closest candidates are:
  DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}(::Any, ::DiffEqJump.AbstractAggregatorAlgorithm, ::DiffEqJump.ConstantRateJump{F1,F2}; kwargs...) at /home/simon/juliapro/JuliaPro-0.5.1.1/JuliaPro/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:9
  DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}(::Any, ::DiffEqJump.AbstractAggregatorAlgorithm, ::Any...; kwargs...) at /home/simon/juliapro/JuliaPro-0.5.1.1/JuliaPro/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:11
  DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}(::P, ::A, ::J<:Union{DiffEqJump.AbstractJumpAggregator,Void}, ::C, ::J2) at /home/simon/juliapro/JuliaPro-0.5.1.1/JuliaPro/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:2

[ChrisRackauckas]

Second problem, that needs a jump aggregation method:

jump_prob = JumpProblem(prob,Direct(),jump,jump2)

For the first problem, the first set of random numbers should be good. But, hopefully this version of the loop clarifies what's up:

nums = Int[]
@time for i in 1:1000
  jump_prob = JumpProblem(prob,Direct(),jump,jump2)
  sol = solve(jump_prob,Discrete(apply_map=false))
  push!(nums,sol[end][3])
end
mean(nums)

(Note: I checked the means against your Gillespie.jl, so unless something has changed, it should end up with the same distribution)

Since the first jump is only calculated when the jump_prob is constructed, a funny thing is happening. First of all, if it It's using the last states of the previous run, causing it to want to die out more. But secondly, it's not going to "redo" the 1st-jump calculation unless you call jump_prob = JumpProblem(prob,Direct(),jump,jump2) . So in order to test different runs, you need to do:

jump_prob = JumpProblem(prob,Direct(),jump,jump2)
sol = solve(jump_prob,Discrete(apply_map=false))

You need to run both of those lines! If you do that, it should work.

I don't know when this got introduced, but this is indeed awful behavior. This is because the jump's code isn't ever touched until it reaches the first timepoint of the jump, and so there isn't a way to input "setup the first jump again". But, if there was an initialization phase for callbacks:

https://github.com/JuliaDiffEq/DiffEqBase.jl/issues/19

then this would be fixed.

[ChrisRackauckas]

And thank you for testing this out! I would like to hear some feedback from you. This setup using the DiscreteCallbacks in the ODE solver is more flexible because it can be attached to any DE, but it is still not as efficient as doing the straight loop. I think I got like 4x from Gillespie.jl last time I checked? The only-jumps version might need to be pulled out and not use the ODE solver logic in order to get more efficient. But I am also wondering if the fact that I am using tuples of tuples is making a difference here. For small equations, tuples of tuples to store the reaction rates probably slows things down considerably over using a dense matrix like Gillespie.jl, but should have an advantage in large sparse models. I don't have a good way of testing how much of a difference is directly due to that choice though.

[ChrisRackauckas]

For reference, https://github.com/JuliaDiffEq/DiffEqJump.jl/commit/590cf3f48cc54faa4d2d302e19e6dbac81dec67b pulled out the initialization from the solve command. But that should go directly into the DE solvers. I'll see if I can get this done right now.

[sdwfrost]

Adding in the initialization step worked, but the latest build throws an error.

using DiffEqJump, DiffEqBase, OrdinaryDiffEq
using Base.Test

rate = (t,u) -> (0.1/1000.0)*u[1]*u[2]
affect! = function (integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
jump = ConstantRateJump(rate,affect!;save_positions=(false,true))

rate = (t,u) -> 0.01u[2]
affect! = function (integrator)
  integrator.u[2] -= 1
  integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,affect!;save_positions=(false,true))

prob = DiscreteProblem([999.0,1.0,0.0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
sol = solve(jump_prob,Discrete(apply_map=false))

nums = Int[]
@time for i in 1:1000
  sol = solve(jump_prob,Discrete(apply_map=false))
  push!(nums,sol[end][1])
end
mean(nums)

MethodError: no method matching DiffEqBase.DiscreteCallback{F1,F2,F3}(::DiffEqJump.DirectJumpAggregation{Float64,Tuple{##1#2,##5#6},Tuple{##3#4,##7#8}}, ::DiffEqJump.DirectJumpAggregation{Float64,Tuple{##1#2,##5#6},Tuple{##3#4,##7#8}}, ::Tuple{Bool,Bool})
Closest candidates are:
  DiffEqBase.DiscreteCallback{F1,F2,F3}{F1,F2,F3}(::F1, ::F2, ::F3, ::Tuple{Bool,Bool}) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqBase/src/callbacks.jl:45
  DiffEqBase.DiscreteCallback{F1,F2,F3}(::Any, ::Any; initialize, save_positions) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqBase/src/callbacks.jl:50
  DiffEqBase.DiscreteCallback{F1,F2,F3}(::Any, ::Any, ::DiffEqJump.ConstantRateJump{F1,F2}, ::Any) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqJump/src/callbacks.jl:31
  ...

 in #JumpProblem#10(::Tuple{Bool,Bool}, ::Type{T}, ::DiffEqBase.DiscreteProblem{Array{Float64,1},Float64,true,DiffEqBase.##130#131}, ::DiffEqJump.Direct, ::DiffEqJump.JumpSet{Tuple{},Tuple{DiffEqJump.ConstantRateJump{##1#2,##3#4},DiffEqJump.ConstantRateJump{##5#6,##7#8}}}) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:23
 in #JumpProblem#9(::Array{Any,1}, ::Type{T}, ::DiffEqBase.DiscreteProblem{Array{Float64,1},Float64,true,DiffEqBase.##130#131}, ::DiffEqJump.Direct, ::DiffEqJump.ConstantRateJump{##1#2,##3#4}, ::Vararg{Any,N}) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:11
 in DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}(::DiffEqBase.DiscreteProblem{Array{Float64,1},Float64,true,DiffEqBase.##130#131}, ::DiffEqJump.Direct, ::DiffEqJump.ConstantRateJump{##1#2,##3#4}, ::DiffEqJump.ConstantRateJump{##5#6,##7#8}) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:11

[sdwfrost]

In addition, with the old commit (I get the same error as above with the latest build) I get different results using the Reaction syntax.

using DifferentialEquations

infection = Reaction(0.1/1000,[1,2],[(1,-1),(2,1)])
recovery = Reaction(0.01,[2],[(2,-1),(3,1)])

srand(1234)
nums = Int[]
@time for i in 1:1000
  sir_prob = DiscreteProblem([999,1,0],(0.0,250.0))
  sir_jump_prob = GillespieProblem(sir_prob,Direct(),infection,recovery)
  sir_sol = solve(sir_jump_prob,Discrete())
  push!(nums,sir_sol[end][3])
end
mean(nums)

[ChrisRackauckas]

Right now that master branch is:

Pkg.checkout("OrdinaryDiffEq")
Pkg.checkout("DiffEqJump")

(along with Stochastic and Delay DiffEq, all for the initialize change, along with retcodes). I'll be pushing tags when tests pass.

My test case is:

using DifferentialEquations

infection = Reaction(0.1/1000,[1,2],[(1,-1),(2,1)])
recovery = Reaction(0.01,[2],[(2,-1),(3,1)])
sir_prob = DiscreteProblem([999,1,0],(0.0,250.0))
sir_jump_prob = GillespieProblem(sir_prob,Direct(),infection,recovery)

sir_sol = solve(sir_jump_prob,Discrete())

using Plots; plot(sir_sol)

srand(1234)
nums = Int[]
@time for i in 1:100000
  sir_sol = solve(sir_jump_prob,Discrete())
  push!(nums,sir_sol[end][3])
end
println("Reaction DSL: $(mean(nums))")

using DiffEqJump, DiffEqBase, OrdinaryDiffEq
using Base.Test

rate = (t,u) -> (0.1/1000.0)*u[1]*u[2]
affect! = function (integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
jump = ConstantRateJump(rate,affect!;save_positions=(false,true))

rate = (t,u) -> 0.01u[2]
affect! = function (integrator)
  integrator.u[2] -= 1
  integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,affect!;save_positions=(false,true))

prob = DiscreteProblem([999.0,1.0,0.0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
sol = solve(jump_prob,Discrete(apply_map=false))

using Plots; plot(sol)

nums = Int[]
@time for i in 1:100000
  sol = solve(jump_prob,Discrete(apply_map=false))
  push!(nums,sol[end][3])
end
println("Direct Jumps: $(mean(nums))")

using Gillespie

function F(x,parms)
  (S,I,R) = x
  (beta,gamma) = parms
  infection = beta*S*I
  recovery = gamma*I
  [infection,recovery]
end

x0 = [999,1,0]
nu = [[-1 1 0];[0 -1 1]]
parms = [0.1/1000.0,0.01]
tf = 250.0
srand(1234)

nums = Int[]
@time for i in 1:100000
  result = ssa(x0,F,nu,parms,tf)
  data = ssa_data(result)
  push!(nums,data[:x3][end])
end
println("Gillespie: $(mean(nums))")

and I get for the means:

Reaction DSL: 725.36465
Direct Jumps: 725.03103
Gillepsie: 725.53742

For the record, the times are:

Reaction DSL: 137.620577 seconds (1.97 G allocations: 87.976 GB, 23.05% gc time)
Direct Jumps: 116.954541 seconds (1.81 G allocations: 63.778 GB, 16.63% gc time)
Gillespie: 25.244438 seconds (1.48 G allocations: 57.974 GB, 34.18% gc time)

so there seems to be a performance regression somewhere changing 2x-4x to 5x-6x. Performance will be another issue though. A lot of it might have to do with the splatting penalty, which seems to keep growing:

https://github.com/JuliaLang/julia/issues/13359

The easiest way around all of this would be generated functions, but in theory there's a way without them

[sdwfrost]

Dear @ChrisRackauckas

Interesting. Looking forward to the tagged releases; still, it's good to know that my implementation gets something for being less friendly. May I suggest using a fixed random number seed and testing on the basis of the output (which is what I do in Gillespie.jl to avoid messing things up)?

[ChrisRackauckas]

May I suggest using a fixed random number seed and testing on the basis of the output (which is what I do in Gillespie.jl to avoid messing things up)?

Yeah, that is something a little finicky too. Juno actually uses some random numbers when setting up its progress bar, so many times when I am benchmarking I get subtle differences when I go to a large problem and turn the progress bar on.

Interesting. Looking forward to the tagged releases; still, it's good to know that my implementation gets something for being less friendly.

Your implementation is much easier to check the correctness of, and it's designed as a straight loop and so it's quite easy to see that it's "nearly optimal". I'm really glad you made it, because it's the perfect test case!

https://github.com/JuliaDiffEq/DiffEqJump.jl/issues/10

I think that in the end we can have something more extensive here, but that also means that the barrier to entry for "how this works" is much higher, and much more testing has to be done to make it robust.

[ChrisRackauckas]

I think all of the releases are tagged now. Those test case should run.

[sdwfrost]

Dear @ChrisRackauckas

Looks OK now...thanks!

[ChrisRackauckas]

Cool. Double checking: what performance numbers do you get? Similar to what I posted above?

[sdwfrost]

Yes, pretty similar - 3-4x performance hit on a Mac Pro too.

[ChrisRackauckas]

Good to know. Sounds like a good universal test case to work on. Thanks for your input. This should make the whole jump much easier to use.