Problem with toy SIR example

[sdwfrost]

Hi, I'm having problems using a toy example from my repository of SIR models. Here's the code, where I'm trying to initialize the system with S=990, I=10:

using FiniteStateProjection
using OrdinaryDiffEq
using Plots

rn = @reaction_network SIR begin
    β, S + I --> 2I
    γ, I --> 0
end β γ

p = [0.0005, 0.25]
u0 = [990.0, 10.0]
tspan = (0.0, 40.0)
solver = Tsit5()

prob_ode = ODEProblem(rn, u0, tspan, p)
sol_ode = solve(prob_ode, solver)
plot(sol_ode)

sys_fsp = FSPSystem(rn)
u0f = zeros(2, 1001)
u0f[1,991] = 1.0
u0f[2,11] = 1.0
prob_fsp = convert(ODEProblem, sys_fsp, u0f, tspan, p)
sol_fsp = solve(prob_fsp, Tsit5())

The model runs quite happily, but gives incorrect answers (a mode of zero at the final timepoint):

julia> sol_fsp.u[end]
2×1001 Matrix{Float64}:
 0.975804  0.0429964    0.00096567  1.44262e-5   1.61612e-7   …  1.75808e-303  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.979782  0.000436904  8.76707e-8  1.04251e-11  8.13531e-16     0.0           0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

I thought it may have something to do with the 'nothing', on the RHS of the second reaction. However, modifying this results in an error on solve:

using FiniteStateProjection
using OrdinaryDiffEq
using Plots

rn = @reaction_network SIR begin
    β, S + I --> 2I
    γ, I --> R
end β γ

p = [0.0005, 0.25]
u0 = [990.0, 10.0, 0.0]
tspan = (0.0, 40.0)
solver = Tsit5()

prob_ode = ODEProblem(rn, u0, tspan, p)
sol_ode = solve(prob_ode, solver)
plot(sol_ode)

sys_fsp = FSPSystem(rn)
u0f = zeros(3, 1001)
u0f[1,991] = 1.0
u0f[2,11] = 1.0
prob_fsp = convert(ODEProblem, sys_fsp, u0f, tspan, p)
sol_fsp = solve(prob_fsp, solver)

Output:

ERROR: MethodError: no method matching pairedindices(::DefaultIndexHandler, ::Matrix{Float64}, ::CartesianIndex{3})

[sdwfrost]

My bad! I thought (based on my poor knowledge of the second example) that I could initialize the variables separately rather than jointly (which also explains why my 3-state model fails). Here's a working version, scaled down to run a bit faster.

using FiniteStateProjection
using OrdinaryDiffEq
using Plots

rn = @reaction_network SIR begin
    β, S + I --> 2I
    γ, I --> 0
end β γ

p = [0.005, 0.25]
u0 = [99, 1]
δt = 1.0
tspan = (0.0, 40.0)
solver = Tsit5()

prob_ode = ODEProblem(rn, u0, tspan, p)
sol_ode = solve(prob_ode, solver)
plot(sol_ode)

sys_fsp = FSPSystem(rn)
u0f = zeros(101, 101)
u0f[100,2] = 1.0
prob_fsp = convert(ODEProblem, sys_fsp, u0f, tspan, p)
sol_fsp = solve(prob_fsp, solver, dense=false, saveat=δt)

p1=plot(0:1:100, sum(sol_fsp.u[1],dims=2),label="S",title="t=0", ylims=(0,1))
plot!(p1, 0:1:100, sum(sol_fsp.u[1],dims=1)',label="I")
p2 = plot(0:1:100, sum(sol_fsp.u[21],dims=2),label="S",title="t="*string(sol_fsp.t[21]),ylims=(0,1))
plot!(p2, 0:1:100, sum(sol_fsp.u[21],dims=1)',label="I")
p3 = plot(0:1:100, sum(sol_fsp.u[41],dims=2),label="S",title="t="*string(sol_fsp.t[41]),ylims=(0,1))
plot!(p3, 0:1:100, sum(sol_fsp.u[41],dims=1)',label="I")
l = @layout [a b c]
plot(p1, p2, p3, layout=l)

[kaandocal]

Thanks for opening this issue! You ran into a common pitfall when trying to solve the CME - I added a rough bit of documentation and a troubleshooting section that may help others avoid it in the future, and you should now get a more informative error message in your second example.