Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Failing to reproduce a model #16

Open
sdwfrost opened this issue Jun 19, 2024 · 0 comments
Open

Failing to reproduce a model #16

sdwfrost opened this issue Jun 19, 2024 · 0 comments

Comments

@sdwfrost
Copy link

sdwfrost commented Jun 19, 2024

I'm trying to reproduce a 'standard' SIR model in ReactiveDynamics, but I get very different results (my hand-coded discrete problem gives the expected result). Can you point me in the direction of where I'm going wrong in defining my model with ReactiveDynamics?

using ReactiveDynamics
using OrdinaryDiffEq
using Random
using Distributions
using Plots

function sir_markov(u,p,t)
    (S, I, R) = u
    (β, γ, N, δt) = p
    infrate = abs*I/N*S*δt)
    recrate = abs*I*δt)
    infection = rand(Poisson(infrate))
    recovery = rand(Poisson(recrate))
    S = S-infection
    I = I+infection-recovery
    R = R+recovery
    [S, I, R]
end

tspan = (0.0,40.0)
u0 = [990, 10, 0] # S, I, R
β = 0.5
γ = 0.25
N = 1000.0
δt = 0.1
p = [β, γ, N, δt] # β, γ, N
seed = 1234
Random.seed!(seed)

prob = DiscreteProblem(sir_markov, u0, tspan, p, dt=δt)
sol = solve(prob, FunctionMap())
plot(sol)

# model dynamics
sir_acs = @ReactionNetwork begin
        β*S*I, S+I --> 2I, name=>infection
        γ*I, I --> R, name=>recovery 
end

# simulation parameters
## initial values
@prob_init sir_acs S=990 I=10 R=0
## uncertainty in initial values (Gaussian)
@prob_uncertainty sir_acs S=0.0 I=0.0
## parameters
@prob_params sir_acs β=0.5 γ=0.25
## other arguments passed to the solver
@prob_meta sir_acs tspan=40 dt=0.01
# turn model into a discrete? problem
sir_prob = @problematize sir_acs
# solve the problem over multiple trajectories
sir_sol = @solve sir_prob
# plot the solution
@plot sir_sol
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant