SIR test gives funny answers #9
Comments
|
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)) |
|
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 = 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: then this would be fixed. |
|
And thank you for testing this out! I would like to hear some feedback from you. This setup using the |
|
For reference, 590cf3f pulled out the initialization from the |
|
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) |
|
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) |
|
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.53742For 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: The easiest way around all of this would be generated functions, but in theory there's a way without them |
|
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)? |
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.
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! 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. |
|
I think all of the releases are tagged now. Those test case should run. |
|
Dear @ChrisRackauckas Looks OK now...thanks! |
|
Cool. Double checking: what performance numbers do you get? Similar to what I posted above? |
|
Yes, pretty similar - 3-4x performance hit on a Mac Pro too. |
|
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. |
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.The text was updated successfully, but these errors were encountered: