Browse Source

Added support for SteadyStateProblems and JumpProblems (#29)

* Added support for SteadStateProblems and JumpProblems

* Fixed compat entries so packages are upgradable
pull/31/head v1.1.0
Micah Halter 10 months ago
committed by GitHub
parent
commit
6feac0e654
No known key found for this signature in database GPG Key ID: 4AEE18F83AFDEB23
  1. 20
      Project.toml
  2. 1
      docs/Project.toml
  3. 1
      examples/Project.toml
  4. 40
      examples/epidemiology.jl
  5. 56
      src/solvers.jl
  6. 106
      test/solvers.jl

20
Project.toml

@ -1,25 +1,31 @@
name = "Petri"
uuid = "4259d249-1051-49fa-8328-3f8ab9391c33"
authors = ["Micah Halter <micah.halter@gtri.gatech.edu>", "James Fairbanks <james.fairbanks@gtri.gatech.edu>"]
version = "1.0.0"
version = "1.1.0"
[deps]
AutoHashEquals = "15f4f7f2-30c1-5605-9d31-71845cf9641f"
Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
DiffEqJump = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
[compat]
AutoHashEquals = "<1.0"
Catlab = "<1.0"
OrdinaryDiffEq = "<6.0"
StochasticDiffEq = "<7.0"
AutoHashEquals = "^0.2.0"
Catlab = "^0.7.0"
DiffEqBase = "^6.40.0"
DiffEqJump = "^6.9.0"
OrdinaryDiffEq = "^5.41.0"
SteadyStateDiffEq = "^1.5.0"
StochasticDiffEq = "^6.24.0"
julia = "1.0"
[extras]
LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[targets]
test = ["LabelledArrays", "OrdinaryDiffEq", "Test"]
test = ["LabelledArrays", "Random", "Test"]

1
docs/Project.toml

@ -1,4 +1,5 @@
[deps]
DiffEqJump = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"

1
examples/Project.toml

@ -1,4 +1,5 @@
[deps]
DiffEqJump = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Petri = "4259d249-1051-49fa-8328-3f8ab9391c33"

40
examples/epidemiology.jl

@ -4,9 +4,10 @@
using Petri
using LabelledArrays
using Plots
using DiffEqJump
using StochasticDiffEq
using OrdinaryDiffEq
using Plots
# ### SIR Model
#
@ -29,22 +30,29 @@ Graph(sir)
# Once a model is defined, we can define out initial parameters `u0`, a time
# span `tspan`, and the transition rates of the interactions `β`
u0 = LVector(S=10.0, I=1.0, R=0.0)
tspan = (0.0,7.5)
β = LVector(inf=0.4, rec=0.4);
u0 = LVector(S=990.0, I=10.0, R=0.0)
tspan = (0.0,40.0)
β = LVector(inf=0.5/sum(u0), rec=0.25);
# Petri.jl provides interfaces to StochasticDiffEq.jl and OrdinaryDiffEq.jl
# Here, we call the `SDEProblem` function that returns an StochasticDiffEq
# problem object and an appropriate Callback set that can be passed to the
# StochasticDiffEq solver which can then be plotted and visualized
# Petri.jl provides interfaces to StochasticDiffEq.jl, DiffEqJump.jl, and
# OrdinaryDiffEq.jl Here, we call the `JumpProblem` function that returns an
# DiffEqJump problem object that can be passed to the DiffEqJump solver which
# can then be plotted and visualized
prob, cb = SDEProblem(sir, u0, tspan, β)
prob = JumpProblem(sir, u0, tspan, β)
sol = DiffEqJump.solve(prob,SSAStepper())
plot(sol)
# Similarly, we can generated `SDEProblem` statements that can be used with
# StochasticDiffEq solvers
prob, cb = SDEProblem(sir, u0, tspan, β)
sol = StochasticDiffEq.solve(prob,SRA1(),callback=cb)
plot(sol)
# Similarly, we can generated `ODEProblem` statements that can be used with
# Lastly, we can generated `ODEProblem` statements that can be used with
# OrdinOrdinaryDiffEq solvers
prob = ODEProblem(sir, u0, tspan, β)
@ -64,9 +72,9 @@ seir = Petri.Model(S, Δ)
Graph(seir)
#-
u0 = LVector(S=10.0, E=1.0, I=0.0, R=0.0)
tspan = (0.0,15.0)
β = LVector(exp=0.9, inf=0.2, rec=0.5)
u0 = LVector(S=990.0, E=10.0, I=0.0, R=0.0)
tspan = (0.0,40.0)
β = LVector(exp=0.7/sum(u0), inf=0.5, rec=0.25)
prob, cb = SDEProblem(seir, u0, tspan, β)
sol = StochasticDiffEq.solve(prob,SRA1(),callback=cb)
@ -86,9 +94,9 @@ seird = Petri.Model(S, Δ)
Graph(seird)
#-
u0 = LVector(S=10.0, E=1.0, I=0.0, R=0.0, D=0.0)
tspan = (0.0,15.0)
β = LVector(exp=0.9, inf=0.2, rec=0.5, die=0.1)
u0 = LVector(S=990.0, E=10.0, I=0.0, R=0.0, D=0.0)
tspan = (0.0,40.0)
β = LVector(exp=0.9/sum(u0), inf=0.9, rec=0.25, die=0.03)
prob, cb = SDEProblem(seird, u0, tspan, β)
sol = StochasticDiffEq.solve(prob,SRA1(),callback=cb)

56
src/solvers.jl

@ -1,12 +1,20 @@
using DiffEqBase
using OrdinaryDiffEq
using SteadyStateDiffEq
using StochasticDiffEq
using DiffEqJump
import OrdinaryDiffEq: ODEProblem
import SteadyStateDiffEq: SteadyStateProblem
import StochasticDiffEq: SDEProblem
import DiffEqJump: JumpProblem
funcindex!(list, key, f, vals...) = list[key] = f(list[key],vals...)
valueat(x::Number, t) = x
valueat(f::Function, t) = f(t)
transitionrate(S, T, k, rate, t) = exp(reduce((x,y)->x+log(S[y] <= 0 ? 0 : S[y]),
keys(first(T[k]));
init=log(valueat(rate[k],t))))
""" vectorfield(m::Model)
@ -19,20 +27,18 @@ function vectorfield(m::Model)
ϕ = Dict()
f(du, u, p, t) = begin
for k in keys(T)
ins = first(T[k])
ϕ[k] = reduce((x,y)->x*u[y]/ins[y], keys(ins); init=valueat(p[k],t))
ϕ[k] = transitionrate(u, T, k, p, t)
end
for s in S
du[s] = 0
end
for k in keys(T)
ins = first(T[k])
outs = last(T[k])
for s in keys(ins)
funcindex!(du, s, -, ϕ[k] * ins[s])
l,r = T[k]
for s in keys(l)
funcindex!(du, s, -, ϕ[k] * l[s])
end
for s in keys(outs)
funcindex!(du, s, +, ϕ[k] * outs[s])
for s in keys(r)
funcindex!(du, s, +, ϕ[k] * r[s])
end
end
return du
@ -46,6 +52,12 @@ Generate an OrdinaryDiffEq ODEProblem
"""
ODEProblem(m::Model, u0, tspan, β) = ODEProblem(vectorfield(m), u0, tspan, β)
""" SteadyStateProblem(m::Model, u0, tspan, β)
Generate an SteadyStateDiffEq SteadyStateProblem
"""
SteadyStateProblem(m::Model, u0, tspan, β) = SteadyStateProblem(ODEProblem(m, u0, tspan, β))
function statecb(s)
cond = (u,t,integrator) -> u[s]
aff = (integrator) -> integrator.u[s] = 0.0
@ -76,7 +88,7 @@ function SDEProblem(m::Model, u0, tspan, β)
sum_u = sum(u)
for k in keys(T)
ins = first(T[k])
ϕ[k] = reduce((x,y)->x*u[y]/(sum_u*ins[y]), keys(ins); init=valueat(p[k],t))
ϕ[k] = transitionrate(u, T, k, p, t)
end
for k in keys(T)
@ -94,3 +106,29 @@ function SDEProblem(m::Model, u0, tspan, β)
return SDEProblem(vectorfield(m),noise,u0,tspan,β,noise_rate_prototype=nu),
CallbackSet([statecb(s) for s in S]...)
end
jumpTransitionRate(T, k) = (u,p,t) -> transitionrate(u, T, k, p, t)
function jumpTransitionFunction(t)
return (integrator) -> begin
l,r = t
for i in keys(l)
integrator.u[i] -= l[i]
end
for i in keys(r)
integrator.u[i] += r[i]
end
end
end
""" JumpProblem(m::Model, u0, tspan, β)
Generate an DiffEqJump JumpProblem
"""
function JumpProblem(m::Model, u0, tspan, β)
T = m.Δ
prob = DiscreteProblem(u0, tspan, β)
return JumpProblem(prob, Direct(), [ConstantRateJump(jumpTransitionRate(T, k),
jumpTransitionFunction(T[k])
) for k in keys(T)]...)
end

106
test/solvers.jl

@ -3,55 +3,34 @@ using Petri
using LabelledArrays
using OrdinaryDiffEq
using StochasticDiffEq
using DiffEqJump
using Random
Random.seed!(1234);
@testset "Generation of ODE formulas" begin
sir = Petri.Model([:S,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
(LVector(I=1), LVector(R=1))])
sirf = vectorfield(sir)
@test sirf(LVector(S=0.0, I=0.0, R=0.0), LVector(S=100.0,I=1.0,R=0.0), [0.35,0.05],0.0) == LVector(S=-35.0,I=34.95,R=0.05)
rates = sirf(LVector(S=0.0, I=0.0, R=0.0), LVector(S=100.0,I=1.0,R=0.0), [0.5/101,0.25],0.0)
@test rates.S -.495 atol=1e-3
@test rates.I .24505 atol=1e-5
@test rates.R .25 atol=1e-2
end
@testset "Generation of ODE solutions" begin
@testset "SIR Solving" begin
sir = Petri.Model([:S,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
(LVector(I=1), LVector(R=1))])
u0 = LVector(S=100.0,I=1.0,R=0.0)
p = [0.35, 0.05]
prob = ODEProblem(sir,u0,(0.0,365.0),p)
sol = OrdinaryDiffEq.solve(prob,Tsit5())
@test sol[end].S < 1
@test sol[end].I < 1
@test sol[end].R > 99
end
@testset "SEIR Solving" begin
seir = Petri.Model([:S,:E,:I,:R],LVector(exp=(LVector(S=1,I=1), LVector(I=1,E=1)),
inf=(LVector(E=1), LVector(I=1)),
rec=(LVector(I=1), LVector(R=1))))
u0 = LVector(S=100.0, E=1.0, I=0.0, R=0.0)
p = LVector(exp=0.35, inf=0.05, rec=0.05)
prob = ODEProblem(seir,u0,(0.0,365.0),p)
u0 = LVector(S=990.0,I=10.0,R=0.0)
p = [0.5/sum(u0), 0.25]
prob = ODEProblem(sir,u0,(0.0,40.0),p)
sol = OrdinaryDiffEq.solve(prob,Tsit5())
@test sol[end].S < 1
@test sol[end].E < 1
@test sol[end].I < 1
@test sol[end].R > 99
end
@testset "SEIRS Solving" begin
seirs = Petri.Model([:S,:E,:I,:R],LVector(exp=(LVector(S=1,I=1), LVector(I=1,E=1)),
inf=(LVector(E=1), LVector(I=1)),
rec=(LVector(I=1), LVector(R=1)),
deg=(LVector(R=1), LVector(S=1))))
u0 = LVector(S=100.0, E=1.0, I=0.0, R=0.0)
p = LVector(exp=0.3, inf=0.4, rec=0.01, deg=0.01)
prob = ODEProblem(seirs,u0,(0.0,60.0),p)
sol = OrdinaryDiffEq.solve(prob,Tsit5())
@test sol[end].S < 1
@test sol[end].E < 1
@test sol[end].I > 60
@test sol[end].R > 30
@test sum(sol[end]) 1000 atol=1e-3
@test sol[end].S 209.843 atol=1e-3
@test sol[end].I 14.474 atol=1e-3
@test sol[end].R 775.684 atol=1e-3
end
end
@ -59,41 +38,28 @@ end
@testset "SIR Solving" begin
sir = Petri.Model([:S,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
(LVector(I=1), LVector(R=1))])
u0 = LVector(S=100.0,I=1.0,R=0.0)
p = [0.35, 0.05]
prob,cb = SDEProblem(sir,u0,(0.0,365.0),p)
sol = StochasticDiffEq.solve(prob,SRA1(),callback=cb)
@test sol[end].S < 1.1
@test sol[end].I < 1.1
@test sol[end].R > 98.9
end
@testset "SEIR Solving" begin
seir = Petri.Model([:S,:E,:I,:R],LVector(exp=(LVector(S=1,I=1), LVector(I=1,E=1)),
inf=(LVector(E=1), LVector(I=1)),
rec=(LVector(I=1), LVector(R=1))))
u0 = LVector(S=100.0, E=1.0, I=0.0, R=0.0)
p = LVector(exp=0.35, inf=0.05, rec=0.05)
prob,cb = SDEProblem(seir,u0,(0.0,365.0),p)
u0 = LVector(S=990.0,I=10.0,R=0.0)
p = [0.5/sum(u0), 0.25]
prob,cb = SDEProblem(sir,u0,(0.0,40.0),p)
sol = StochasticDiffEq.solve(prob,SRA1(),callback=cb)
@test sol[end].S < 1.1
@test sol[end].E < 1.1
@test sol[end].I < 1.1
@test sol[end].R > 98.9
@test sum(sol[end]) 1000 atol=1e-3
@test sol[end].S 186.054 atol=1e-3
@test sol[end].I 13.447 atol=1e-3
@test sol[end].R 800.499 atol=1e-3
end
end
@testset "SEIRS Solving" begin
seirs = Petri.Model([:S,:E,:I,:R],LVector(exp=(LVector(S=1,I=1), LVector(I=1,E=1)),
inf=(LVector(E=1), LVector(I=1)),
rec=(LVector(I=1), LVector(R=1)),
deg=(LVector(R=1), LVector(S=1))))
u0 = LVector(S=100.0, E=1.0, I=0.0, R=0.0)
p = LVector(exp=0.3, inf=0.4, rec=0.01, deg=0.01)
prob,cb = SDEProblem(seirs,u0,(0.0,60.0),p)
sol = StochasticDiffEq.solve(prob,SRA1(),callback=cb)
@test sol[end].S < 1.1
@test sol[end].E < 1.1
@test sol[end].I > 59.9
@test sol[end].R > 29.9
@testset "Generation of Jump solutions" begin
@testset "SIR Solving" begin
sir = Petri.Model([:S,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
(LVector(I=1), LVector(R=1))])
u0 = LVector(S=990.0,I=10.0,R=0.0)
p = [0.5/sum(u0), 0.25]
prob = JumpProblem(sir,u0,(0.0,40.0),p)
sol = DiffEqJump.solve(prob,SSAStepper())
@test sum(sol[end]) == 1000
@test sol[end].S == 263
@test sol[end].I == 10
@test sol[end].R == 727
end
end
end
Loading…
Cancel
Save