Browse Source

[WIP] Add a theory for modeling SEIRD diseases w/ primitive transmision model (covid?) (#241)

* Add a theory for modeling SEIRD diseases primitive transmision model

* update covid example with operadic composition and SEIRD model

* Update catlab dependency

* Fix typo

* Updated dependency versions

* Updated catlab syntax and petri net syntax throughout files

* Updated packages

* attempt to update malaria example. Fixes problems with OpenPetri.model.model

* More fixes to malaria.ipynb

* fix the mk_function calls in malaria.ipynb

* update wd render to orientation from direction. update malaria example

* Add some TODOs to the CT api

* attempt to use the new CT API for petri composition.

* add implementation of Covid model without CT interface.

* figured out how to use the new API

* introduce PetriCospan instance, tests pass

* WIP saving state

* fix test failure on Death Sq

* fix compose_pushout attempt 1

* remove compose_pushout debug stmts

* collecting tests at the bottom of the file

* move covid example to bottom of file

* remove unnecessarily tight type constraint

* add tests for building SEI and SEIR models

* comment out drawing to make testing easier

* WIP building intercity flow model

* add more PetriCospan tests

* add covid model and parameter exploration

* refactor covid odemodel

* move covid images to a `img` folder

* ignore images in covid example directory

* add SEIRD case for COVID

* fix plotting

* add peakgap analysis of SEIRD multicity

* delete unused test file for catlab

* add Pseird city model to covid example

* update project.toml and tests

* add comments to odemodel.jl

* integrate PetriCospans into SemanticModels add tests

* ammend previous commit

* add estimate delay to odemodel.jl

* add SEIR vs SEIRD comparison based on gap

* add validation data generation script

Co-authored-by: Micah Halter <micah.halter@gtri.gatech.edu>
pull/247/head
James 2 years ago committed by GitHub
parent
commit
5c0e9c9ba8
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
  1. 1
      .gitignore
  2. 24
      Project.toml
  3. 1
      REQUIRE
  4. 4
      examples/decorations/decoration-demo.jl
  5. 6
      examples/decorations/relolog.jl
  6. 4
      examples/epicookbook/src/SEIRmodel.jl
  7. 6
      examples/petri/covid/Project.toml
  8. 242
      examples/petri/covid/covid.jl
  9. 31
      examples/petri/covid/fallback.jl
  10. 466
      examples/petri/covid/odemodel.jl
  11. 169
      examples/petri/covid/validation.jl
  12. 6984
      examples/petri/malaria.ipynb
  13. 1
      examples/petri/malaria.jl
  14. 6
      examples/petri/petri_demo.jl
  15. 92
      examples/petri/rewrite.jl
  16. 73
      examples/petri/wiring_ode.jl
  17. 20
      examples/petri/wiring_petri.jl
  18. 5
      src/modeltools/CategoryTheory.jl
  19. 1
      src/modeltools/ModelTools.jl
  20. 14
      src/modeltools/OpenModels.jl
  21. 52
      src/modeltools/OpenPetris.jl
  22. 260
      src/modeltools/PetriCospans.jl
  23. 42
      src/modeltools/PetriModels.jl
  24. 8
      src/modeltools/WiringDiagrams.jl
  25. 226
      src/petri.jl
  26. 329
      test/catlab.jl
  27. 215
      test/modeltools/epidemics.jl
  28. 210
      test/modeltools/petricospans.jl
  29. 2
      test/odegraft.jl
  30. 11
      test/runtests.jl
  31. 4
      test/transform/ode.jl

1
.gitignore vendored

@ -277,3 +277,4 @@ examples/*.dot @@ -277,3 +277,4 @@ examples/*.dot
examples/*.dot.svg
doc/src/img/petri/*.png
/doc/workflow.slides.html
/examples/petri/covid/img/

24
Project.toml

@ -8,15 +8,17 @@ BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" @@ -8,15 +8,17 @@ BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
Cassette = "7057c7e9-c182-5462-911a-8362d720325c"
Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
Compose = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
GraphDataFrameBridge = "3c71623a-a715-5176-9801-629b201a4880"
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@ -24,26 +26,25 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" @@ -24,26 +26,25 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Petri = "4259d249-1051-49fa-8328-3f8ab9391c33"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
TikzPictures = "37f6aa50-8035-52d0-81c2-5a1d08754b2d"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
[compat]
julia = "1.0"
BoundaryValueDiffEq = "2.3.0"
CSV = "<1.0"
Cassette = "<1.0"
Catlab = "0.3"
Colors = "<1.0"
Catlab = "<1.0"
DataFrames = "<1.0"
DataFramesMeta = "<1.0"
DelimitedFiles = "<2.0"
Distributions = "<1.0"
Distributions = "0.23.0"
Documenter = "<1.0"
GLM = "<2.0"
GeneralizedGenerated = "0.1.4"
GraphDataFrameBridge = "<1.0"
JSON = "<1.0"
Latexify = "<1.0"
@ -52,22 +53,21 @@ LinearAlgebra = "<2.0" @@ -52,22 +53,21 @@ LinearAlgebra = "<2.0"
Logging = "<2.0"
MacroTools = "<1.0"
MetaGraphs = "<1.0"
ModelingToolkit = "0.6"
Petri = "<1.0"
Petri = "0.1.3"
Plots = "<1.0"
Random = "<2.0"
StaticArrays = "0.11"
StochasticDiffEq = "6.10.0"
Unitful = "<1.0"
julia = "1.0"
[extras]
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[targets]
test = ["DifferentialEquations", "LsqFit", "Polynomials", "Printf", "Statistics", "Test"]
test = ["OrdinaryDiffEq", "LsqFit", "Polynomials", "Printf", "Statistics", "Test"]

1
REQUIRE

@ -1,7 +1,6 @@ @@ -1,7 +1,6 @@
julia 1.0.0
CSV
Cassette
Colors
Distributions
Documenter
LightGraphs

4
examples/decorations/decoration-demo.jl

@ -45,7 +45,7 @@ Petri.Graph(devar(sir_petri)) @@ -45,7 +45,7 @@ Petri.Graph(devar(sir_petri))
S::Ob
I::Ob
R::Ob
infects::Hom(I, S)
recovers_to::Hom(I, R)
end
@ -68,7 +68,7 @@ Petri.Graph(devar(sii_petri)) @@ -68,7 +68,7 @@ Petri.Graph(devar(sii_petri))
S::Ob
I::Ob
I_p::Ob
infects::Hom(I, S)
infects_p::Hom(I_p, S)
end

6
examples/decorations/relolog.jl

@ -3,7 +3,7 @@ using Catlab @@ -3,7 +3,7 @@ using Catlab
using Catlab.Doctrines
import TikzPictures
import Catlab.Graphics: to_tikz
import SemanticModels.ModelTools.RelOlogModels: RelOlogModel, model, ,
import SemanticModels.ModelTools.RelOlogModels: RelOlogModel, model,
using SemanticModels.ModelTools.CategoryTheory
import SemanticModels.ModelTools.CategoryTheory: , FinSetMorph
@ -15,12 +15,12 @@ import SemanticModels.ModelTools.CategoryTheory: ⊔, FinSetMorph @@ -15,12 +15,12 @@ import SemanticModels.ModelTools.CategoryTheory: ⊔, FinSetMorph
E::Ob
I::Ob
R::Ob
exposes::Hom(I, S)
becomes_exposed::Hom(S, E)
falls_ill::Hom(E, I)
recovers_to::Hom(I, R)
illness := compose(becomes_exposed, falls_ill)
exposure := compose(exposes, becomes_exposed)
end

4
examples/epicookbook/src/SEIRmodel.jl

@ -14,7 +14,7 @@ @@ -14,7 +14,7 @@
# ---
module SEIRmodel
using DifferentialEquations
using OrdinaryDiffEq
#Susceptible-exposed-infected-recovered model function
function seir_ode(dY,Y,p,t)
@ -50,7 +50,7 @@ function main() @@ -50,7 +50,7 @@ function main()
seir_prob = ODEProblem(seir_ode,init,tspan,pram)
sol=solve(seir_prob);
sol=solve(seir_prob, alg=Tsit5());
return sol
end

6
examples/petri/covid/Project.toml

@ -0,0 +1,6 @@ @@ -0,0 +1,6 @@
[deps]
Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
Compose = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
TikzPictures = "37f6aa50-8035-52d0-81c2-5a1d08754b2d"

242
examples/petri/covid/covid.jl

@ -0,0 +1,242 @@ @@ -0,0 +1,242 @@
using Catlab
using Catlab.Doctrines
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.Programs
import Base.Multimedia: display
import Catlab.Graphics: to_graphviz, LeftToRight
println("Done Importing Catlab")
import Base: (==), length, show
using Test
using Petri
using SemanticModels.ModelTools.CategoryTheory
import SemanticModels.ModelTools.CategoryTheory: undecorate,
using SemanticModels.ModelTools.PetriModels
using SemanticModels.ModelTools.PetriCospans
import SemanticModels.ModelTools.PetriCospans: otimes_ipm, compose_pushout
println("Done Importing SemanticModels")
import Catlab.Doctrines:
Ob, Hom, dom, codom, compose, , , id, oplus, otimes, , , munit, mzero, braid,
dagger, dunit, dcounit, mcopy, Δ, delete, , mmerge, , create, ,
plus, zero, coplus, cozero, meet, top, join, bottom
wd = to_wiring_diagram
draw(d::WiringDiagram) = to_graphviz(add_junctions(d), orientation=LeftToRight, labels=true)
draw(d::HomExpr) = draw(wd(d))
println("Making COVID Model")
@present Disease(FreeEpidemiology) begin
S::Ob
E::Ob
I::Ob
R::Ob
end
S,E,I,R = generators(Disease)
# draw(spontaneous(E,I)⋅spontaneous(I,R)⊗(spontaneous(E,I)))
sei = compose(exposure(S,I,E), otimes(spontaneous(E,I), id(I)), mmerge(I))
# draw(sei)
seir = sei⋅Δ(I)(id(I)spontaneous(I, R))
# draw(seir)
seir2 = compose(mcopy(S)id(I), id(S)seir)
# draw(seir2)
# draw(death(S))
d = @program Disease (s::S, e::E, i::I) begin
e1, i1 = exposure{S,I,E}(s,i)
i2 = spontaneous{E,I}(e1)
e = [e, e1]
e_out = spontaneous{E,E}(e)
i1 = [i1, i2]
r = spontaneous{I,R}(i1)
s_out = spontaneous{S,S}(s)
return s_out, e_out, spontaneous{I,I}(i1)
end
# draw(d)
# draw(d⋅d)
seirdef = to_hom_expr(FreeEpidemiology, d)
try
add_definition!(Disease, :seir, seirdef)
catch
println("INFO: definition already added.")
end
# if the disease is fatal, we need to add a death component
seird = @program Disease (s::S, e::E, i::I) begin
e1, i1 = exposure{S,I,E}(s,i)
i2 = spontaneous{E,I}(e1)
e = [e, e1]
e_out = spontaneous{E,E}(e)
i1 = [i1, i2]
r = spontaneous{I,R}(i1)
s_out = spontaneous{S,S}(s)
death{I}(i1)
return s_out, e_out, spontaneous{I,I}(i1)
end
#TODO: This does not get translated correctly, bug?
seirddef = to_hom_expr(FreeEpidemiology, seird)
try
add_definition!(Disease, :seird, seirddef)
catch
println("INFO: definition already added.")
end
seirgen = generator(Disease, :seir)
seirdgen = generator(Disease, :seird)
ncities(city,n::Int) = compose([city for i in 1:n]...)
city³ = ncities(seirgen, 3)
# draw(city³)
dcity³ = wd(city³)
dc3 = substitute(dcity³, box_ids(dcity³), [d,d,d])
@show dc3 == ncities(d, 3)
# draw(dc3)
import Base: repeat
repeat(d::WiringDiagram, n::Int) = compose([d for i in 1:n]...)
repeat(d::FreeEpidemiology.Hom, n::Int) = compose([d for i in 1:n]...)
# draw(ncities(seirdgen, 3))
# draw(repeat(seird, 3))
# draw(seirddef)
#
# using TikzPictures
# using Catlab.Graphics.TikZWiringDiagrams
# using Convex
# using SCS
#
# to_tikz(seirddef, labels=true)
X = FinSet(1)
F(ex) = functor((FinSet, FinSetMorph),
ex,
generators=Dict(S=>X, I=>X, E=>X, R=>X))
Fseir = compose(exposure(X,X,X),otimes(spontaneous(X,X),id(X)),mmerge(X),mcopy(X),otimes(id(X),spontaneous(X,X)))
# states are [S, I, E, R]
seir_petri = left(Fseir.f).d[1]
f = FinSetMorph(1:4, [1, 2, 3])
g = FinSetMorph(1:4, [1, 2, 3])
Fseir′ = PetriCospan(Cospan(Decorated(f,seir_petri), Decorated(g, seir_petri)))
s = spontaneous(X,X)
Fflow = otimes(s, s, s)
Fcity = Fseir′ Fflow
Fcity³ = Fcity⋅Fcity⋅Fcity
@test left(Fcity³.f).d[1].model.S == 1:15
@test left(Fcity³.f).d[1].model.Δ == [([1, 2], [3, 2]),
([3], [2]),
([2], [4]),
#outflow 1
([1], [5]),
([2], [6]),
([3], [7]),
# SEIR 2
([5, 6], [7, 6]),
([7], [6]),
([6], [8]),
#outflow 2
([5], [9]),
([6], [10]),
([7], [11]),
# SEIR 3
([9, 10], [11, 10]),
([11], [10]),
([10], [12]),
#outflow 3
([9], [13]),
([10], [14]),
([11], [15])
]
Fcity₀ = Fseir′ Fflow
Fcity₁ = Fflow Fseir′ Fflow
Fcityₑ = Fflow Fseir′
Fcity³ = Fcity₀⋅Fcity₁⋅Fcityₑ
@test left(Fcity³.f).d[1].model.S == 1:18
@test left(Fcity³.f).d[1].model.Δ == [
# SEIR 1
([1, 2], [3, 2]),
([3], [2]),
([2], [4]),
# outflow 1→2
([1], [5]),
([2], [6]),
([3], [7]),
# inflow 1→2
([5], [8]),
([6], [9]),
([7], [10]),
# SEIR 2
([8, 9], [10, 9]),
([10], [9]),
([9], [11]),
# outflow 2→3
([8], [12]),
([9], [13]),
([10], [14]),
# inflow 2→3
([12], [15]),
([13], [16]),
([14], [17]),
# SEIR 3
([15, 16], [17, 16]),
([17], [16]),
([16], [18])
]
Fcity₀ = Fseir′ Fflow
Fcity₁ = Fseir′ Fflow
Fcityₑ = Fseir′
Fcity³ = Fcity₀⋅Fcity₁⋅Fcityₑ
@test left(Fcity³.f).d[1].model.S == 1:12
@test left(Fcity³.f).d[1].model.Δ == [
# City 1 SEIR
([1, 2], [3, 2]),
([3], [2]), # E→I
([2], [4]), # I→R
# outflow 1→2
([1], [5]), # S→S′
([2], [6]), # I→I′
([3], [7]), # E→E′
# City 2 SEIR
([5, 6], [7, 6]),
([7], [6]), # E→I
([6], [8]), # I→R
# outflow 2→3
([5], [9]), # S→S′
([6], [10]),# I→I′
([7], [11]),# E→E′
# City 3 SEIR
([9, 10], [11, 10]),
([11], [10]), # E→I
([10], [12]) # I→R
]
Pseird = PetriModel(
Petri.Model(1:5,[
([1,2],[3,2]), # exposure
([3],[2]), # onset
([2],[4]), # recovery
([2],[5]), # death
], missing, missing))
inputs = FinSetMorph(1:5, [1,2,3])
outputs = FinSetMorph(1:5, [1,2,3])
Fcityd = PetriCospan(Cospan(Decorated(inputs, Pseird),
Decorated(outputs, Pseird)))
Fcity₀ = Fcityd Fflow
Fcity₁ = Fcityd Fflow
Fcityₑ = Fcityd
Fcity³ = Fcity₀⋅Fcity₁⋅Fcityₑ
left(Fcity³.f).d[1].model.S == 1:15
length(left(Fcity³.f).d[1].model.Δ) == 18

31
examples/petri/covid/fallback.jl

@ -0,0 +1,31 @@ @@ -0,0 +1,31 @@
using ModelingToolkit
using SemanticModels.ModelTools
using SemanticModels.ModelTools.PetriModels
using SemanticModels.ModelTools.OpenPetris
using SemanticModels.ModelTools.OpenModels
using Petri
MAX_STATES = 100
X = @variables(X[1:MAX_STATES])[1]
STATELOOKUP = Dict(s.op.name=>i for (i,s) in enumerate(X))
#const OpenPetri{V} = OpenModel{V, PetriModel}
Petri.N(s) = 1
SEIRD = Petri.Model(1:5, [
(X[1]+X[2], X[3]+X[2]),
(X[3], X[2]),
(X[2], X[4]),
(X[2], X[5])
], missing, missing)
flows = Petri.Model(1:3, [
(X[1], X[4]),
(X[2], X[5]),
(X[3], X[6])], missing, missing
)
seirdpetri = OpenModel([1,2,3], PetriModel(SEIRD), [1,2,3])
openflows = OpenModel([1,2,3], PetriModel(flows), [4,5,6])
city = OpenPetris.compose(seirdpetri, openflows)
city³ = compose(city, city, city)

466
examples/petri/covid/odemodel.jl

@ -0,0 +1,466 @@ @@ -0,0 +1,466 @@
using OrdinaryDiffEq
import OrdinaryDiffEq: ODEProblem
using Petri
using Test
#this function was removed in Julia 0.7
""" linreg(x,y)
perform a simple linear regression for scalar case. Coefficients are returned
as (intercept, scalar)
"""
linreg(x, y) = hcat(fill!(similar(x), 1), x) \ y
""" peakgap(s1, s2)
identify the time delay between the peak in two time series. Time series should
be stored as (x(t), t) for consistency with the tuples(sol) from OrdinaryDiffEq.
"""
function peakgap(s1, s2)
@show p1, i1 = findmax(first.(s1))
@show p2, i2 = findmax(first.(s2))
return s2[i2][end] - s1[i1][end]
end
""" peakgap(sol, i::Int, j::Int)
identify the time delay between the peak in two dimensions of an ODE solution structure.
"""
function peakgap(sol, i::Int, j::Int)
s1 = [(u[i],t) for (u,t) in tuples(sol)]
s2 = [(u[j],t) for (u,t) in tuples(sol)]
return peakgap(s1,s2)
end
""" paramsweep(f::Function, m::Petri.Model, initials, tspan, params)
solve an ODEProblem for each parameter setting in params and apply the function f to the solution.
Usage:
paramsweep(m,initials,tspan,[p1,p2,p3]) do sol
return f(sol)
end
"""
function paramsweep(f::Function, m::Petri.Model, initials, tspan, params)
map(params) do p
prob = ODEProblem(fluxes(m), initials, tspan, p)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
return f(sol)
end
end
""" fluxes(m::Petri.Model)
a PetriNet interpreter that computes the mass action kinetics of a petri net.
Usage: pass this to ODEProblem to set up an ODE for a given model.
"""
function fluxes(m::Petri.Model)
S = m.S
T = m.Δ
nS = length(S)
nT = length(T)
ϕ = zeros(Float64, nT)
f(du, u, p, t) = begin
for (i, t) in enumerate(T)
ins = t[1]
# TODO: accomodate multiplicites here
ϕ[i] = p[i]*prod(u[ins])
end
for i in S
du[i] = 0
end
for (i, t) in enumerate(T)
ins = t[1]
out = t[2]
for s in ins
# TODO: accomodate multiplicites here
du[s] -= ϕ[i]
end
for s in out
# TODO: accomodate multiplicites here
du[s] += ϕ[i]
end
end
return du
end
return f
end
u₀(m::Petri.Model) = begin
zeros(Float64, length(m.S))
end
u₀(m::Petri.Model, initialS, initialI=1) = begin
u0=zeros(Float64, length(m.S))
u0[1] = initialS
u0[2] = initialI
return u0
end
function makeplots_seir(sol, prefix)
mkpath(dirname(prefix))
p1 = plot(sol,vars=[1,2,3,4], xlabel="", ylabel="people", linewidth=3,title="Cities", legend=false)
p2 = plot(sol,vars=[5,6,7,8], xlabel="", ylabel="people", linewidth=3, legend=false)
p3 = plot(sol,vars=[9,10,11,12], xlabel="time", ylabel="people", linewidth=3, legend=false)
p4 = plot(sol,vars=[2,6,10], xlabel="", linewidth=3, labels=["i1" "i2" "i3"], legend=true)
p5 = plot(sol,vars=[3,7,11], xlabel="", linewidth=3,title="Populations", labels=["e1" "e2" "e3"], legend=true)
p6 = plot(sol,vars=[4,8,12], xlabel="time", linewidth=3, labels=["r1" "r2" "r3"], legend=true)
p = plot(p1, p5, p2, p4, p3, p6, layout=(3,2), linewidth=3, link=:both)
savefig(p, "$(prefix)combined.pdf")
p
end
function makeplots_seird(sol, prefix)
mkpath(dirname(prefix))
p1 = plot(sol,vars=[1,2,3,4,5], xlabel="", ylabel="people", linewidth=3,title="Cities", legend=false)
p2 = plot(sol,vars=[6,7,8,9,10], xlabel="", ylabel="people", linewidth=3, legend=false)
p3 = plot(sol,vars=[11,12,13,14,15], xlabel="time", ylabel="people", linewidth=3, legend=false)
p4 = plot(sol,vars=[2,7,12], xlabel="", linewidth=3, labels=["i1" "i2" "i3"], legend=true)
p5 = plot(sol,vars=[3,8,13], xlabel="", linewidth=3,title="Populations", labels=["e1" "e2" "e3"], legend=true)
p6 = plot(sol,vars=[5,10,15], xlabel="time", linewidth=3, labels=["d1" "d2" "d3"], legend=true)
p = plot(p1, p5, p2, p4, p3, p6, layout=(3,2), linewidth=3, link=:both)
savefig(p, "$(prefix)combined.pdf")
p
end
T = [
# City 1 SEIR
([1, 2], [3, 2]),
([3], [2]), # E→I
([2], [4]), # I→R
# outflow 1→2
([1], [5]), # S→S′
([2], [6]), # I→I′
([3], [7]), # E→E′
# City 2 SEIR
([5, 6], [7, 6]),
([7], [6]), # E→I
([6], [8]), # I→R
# outflow 2→3
([5], [9]), # S→S′
([6], [10]),# I→I′
([7], [11]),# E→E′
# City 3 SEIR
([9, 10], [11, 10]),
([11], [10]), # E→I
([10], [12]) # I→R
]
m = Petri.Model(1:12, T, missing, missing)
nS = length(m.S)
β = ones(Float64, length(T))
@test fluxes(m)(zeros(Float64, nS), ones(Float64, nS), β, 1) |> length == length(m.S)
tspan = (0.0,60.0)
prob = ODEProblem(fluxes(m), u₀(m, 1,0), tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@test sol.u[end][end-3] > .90
@test sol.u[end][end-2] == 0
@test sol.u[end][end-1] == 0
@test sol.u[end][end] == 0
u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[5] = 10000
u0[9] = 10000
u0[2] = 1
u0
βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/2, 1/2, 1/2]/1000
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
@show β
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]
using Plots
makeplots_seir(sol, "img/seir/baseline")
βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/2, 1/2, 1/2]/100
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]
makeplots_seir(sol, "img/seir/travel")
βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/2, 1/200, 1/2]/100
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]
makeplots_seir(sol, "img/seir/screening")
βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/200, 1/200, 1/200]/100000
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]
makeplots_seir(sol, "img/seir/shutdown")
println("Exploring the SEIRD^3 example")
S = 1:15
T = [([1, 2], [3, 2]),
([3], [2]),
([2], [4]),
([2], [5]),
([1], [6]),
([2], [7]),
([3], [8]),
([6, 7], [8, 7]),
([8], [7]),
([7], [9]),
([7], [10]),
([6], [11]),
([7], [12]),
([8], [13]),
([11, 12], [13, 12]),
([13], [12]),
([12], [14]),
([12], [15])]
m = Petri.Model(S, T, missing, missing)
u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[6] = 10000
u0[11] = 10000
u0[2] = 1
u0
βseir = [10/sum(u0), 1/2, 1/5, 1/10]
βtravel = [1/2, 1/2, 1/2]/1000
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
@show β
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]
makeplots_seird(sol, "img/seird/baseline")
function paramsweep(f::Function, m::Petri.Model, initials, tspan, params)
map(params) do p
prob = ODEProblem(fluxes(m), initials, tspan, p)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
return f(sol)
end
end
packparams(βlocal, βflow) = vcat(βlocal, βflow, βlocal, βflow, βlocal)
# The βseir parameters are the mass action kinetics rates of
# 1. S+I --> E+I
# 2. E --> I
# 3. I --> R
# The βtravel parameters are the mass action kinetics rates of outflow to the next
# city in the transit network.
# 1. S --> S⁺
# 2. I --> I⁺
# 3. E --> E⁺
# Baseline scenario
βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/2, 1/2, 1/2]/100
β1 = packparams(βseir, βtravel)
# Travel scenario
βtravel = [1/2, 1/2, 1/2]/10
β2 = packparams(βseir, βtravel)
# Screening scenario
βtravel = [1/2, 1/20, 1/2]/100
β3 = packparams(βseir, βtravel)
# Shutdown scenario
βtravel = [1/2, 1/2, 1/2]/1000
β4 = packparams(βseir, βtravel)
# Total and Complete Shutdown scenario
βtravel = [1/2, 1/2, 1/2]/10000
β5 = packparams(βseir, βtravel)
βs = [β1, β2, β3, β4, β5]
T = [
# City 1 SEIR
([1, 2], [3, 2]),
([3], [2]), # E→I
([2], [4]), # I→R
# outflow 1→2
([1], [5]), # S→S′
([2], [6]), # I→I′
([3], [7]), # E→E′
# City 2 SEIR
([5, 6], [7, 6]),
([7], [6]), # E→I
([6], [8]), # I→R
# outflow 2→3
([5], [9]), # S→S′
([6], [10]),# I→I′
([7], [11]),# E→E′
# City 3 SEIR
([9, 10], [11, 10]),
([11], [10]), # E→I
([10], [12]) # I→R
]
# u0 is the initial configuration of the system
# u0[[1,5,9]] = number of initially Susceptible people in each population
# u0[2] = number of Infected individuals
u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[5] = 10000
u0[9] = 10000
u0[2] = 1
u0
m = Petri.Model(1:12, T, missing, missing)
gaps = paramsweep(sol->peakgap(sol, 2, 10), m, u0, tspan, βs)
# prob = ODEProblem(fluxes(m), u0, tspan, βs[1])
# sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
# makeplots_seir(sol, "img/debug_1")
travelparams(k::Number) = begin
βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/2, 1/2, 1/2]/(100k)
β = packparams(βseir, βtravel)
return β
end
travelsweep(r::AbstractVector) = begin
return travelparams.(2 .^ r)
end
βrange = 1:8
gaps = paramsweep(sol->peakgap(sol, 2, 10), m, u0, tspan, travelsweep(βrange))
coeffs = linreg(βrange, gaps)
@test sum(abs.((coeffs[1] .+ (coeffs[2]*(βrange))) .- gaps))/length(βrange)<= 1e-1
# struct GapResult{P,G}
# param::P
# gap::G
# end
#
# GapResult(t::Tuple) = GapResult(t...)
#
# struct GapStudy
# result::Vector
# coeffs::Vector
# end
#
# slope(gs::GapStudy) = gs.coeffs[end]
function estimatedelay(r::AbstractVector)
gaps = paramsweep(sol->peakgap(sol, 2, 10), m, u0, tspan, travelsweep(r))
coeffs = linreg(βrange, gaps)
return (results=zip(r, gaps), coeffs=coeffs)
end
function estimatedelay(m::Petri.Model, u0, tspan, parameters)
gaps = paramsweep(sol->peakgap(sol, 2, 10), m, u0, tspan, parameters)
coeffs = linreg(βrange, gaps)
return (results=zip(parameters, gaps), coeffs=coeffs)
end
slope(gs) = gs.coeffs[end]
gs = estimatedelay(βrange)
@test slope(gs) > 1
collect(gs.results)
gs.coeffs
using Catlab
using Catlab.Doctrines
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.Programs
using SemanticModels.ModelTools.CategoryTheory
import SemanticModels.ModelTools.CategoryTheory: undecorate,
using SemanticModels.ModelTools.PetriModels
using SemanticModels.ModelTools.PetriCospans
import SemanticModels.ModelTools: model
model(c::PetriCospan) = left(c.f).d[1].model
X = FinSet(1)
Fseir = compose(exposure(X,X,X),otimes(spontaneous(X,X),id(X)),mmerge(X),mcopy(X),otimes(id(X),spontaneous(X,X)))
# states are [S, I, E, R]
seir_petri = left(Fseir.f).d[1]
f = FinSetMorph(1:4, [1, 2, 3])
g = FinSetMorph(1:4, [1, 2, 3])
Fseir′ = PetriCospan(Cospan(Decorated(f,seir_petri), Decorated(g, seir_petri)))
s = spontaneous(X,X)
Fflow = otimes(s, s, s)
@test length(model(Fseir′).S) == 4
@test length(model(Fseir′).Δ) == 3
m = model(Fseir′ Fflow Fseir′ Fflow Fseir′)
u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[5] = 10000
u0[9] = 10000
u0[2] = 1
res = estimatedelay(m::Petri.Model, u0, tspan, travelsweep(1:8))
seirgaps = collect(map(last, res.results))
@test 1.35 <= slope(res) <= 1.45
Pseird = PetriModel(
Petri.Model(1:5,[
([1,2],[3,2]), # exposure
([3],[2]), # onset
([2],[4]), # recovery
([2],[5]), # death
], missing, missing))
inputs = FinSetMorph(1:5, [1,2,3])
outputs = FinSetMorph(1:5, [1,2,3])
Fcityd = PetriCospan(Cospan(Decorated(inputs, Pseird),
Decorated(outputs, Pseird)))
s = spontaneous(X,X)
Fflow = otimes(s, s, s)
Fcity₀ = Fcityd Fflow
Fcity₁ = Fcityd Fflow
Fcityₑ = Fcityd
Fcity³ = Fcity₀⋅Fcity₁⋅Fcityₑ
m = model(Fcity³)
u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[6] = 10000
u0[11] = 10000
u0[2] = 1
seirdparams(k) = begin
βseird = [10/sum(u0), 1/2, 1/5, 1/16]
βtravel = [1/2, 1/2, 1/2]/(100k)
β = vcat(βseird, βtravel)
return vcat(β, β, β)
end
gaps = paramsweep(sol->peakgap(sol, 2, 12), m,
u0, tspan, seirdparams.(2 .^ (1:8)))
prob = ODEProblem(fluxes(model(Fcity³)), u0, tspan, seirdparams(2))
sol = OrdinaryDiffEq.solve(prob, alg=Tsit5())
seirdgaps = gaps
@test model(Fcity³).S == 1:15
@test length(model(Fcity³).Δ) == 18
makeplots_seird(sol, "tmp")
gp = plot([seirgaps seirdgaps], marker=:square, labels=[:seir :seird], xlabel="2^x fold reduction in travel", ylabel="Gap between first and last epidemic", legend=:bottomright, title="Effect of model assumptions on travel restrictions")
savefig(gp, "gapplot.pdf")

169
examples/petri/covid/validation.jl

@ -0,0 +1,169 @@ @@ -0,0 +1,169 @@
using OrdinaryDiffEq
import OrdinaryDiffEq: ODEProblem
using Petri
using Test
using Catlab
using Catlab.Doctrines
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.Programs
using SemanticModels.ModelTools.CategoryTheory
import SemanticModels.ModelTools.CategoryTheory: undecorate,
using SemanticModels.ModelTools.PetriModels
using SemanticModels.ModelTools.PetriCospans
import SemanticModels.ModelTools: model
""" fluxes(m::Petri.Model)
a PetriNet interpreter that computes the mass action kinetics of a petri net.
Usage: pass this to ODEProblem to set up an ODE for a given model.
"""
function fluxes(m::Petri.Model)
S = m.S
T = m.Δ
nS = length(S)
nT = length(T)
ϕ = zeros(Float64, nT)
f(du, u, p, t) = begin
for (i, t) in enumerate(T)
ins = t[1]
# TODO: accomodate multiplicites here
ϕ[i] = p[i]*prod(u[ins])
end
for i in S
du[i] = 0
end
for (i, t) in enumerate(T)
ins = t[1]
out = t[2]
for s in ins
# TODO: accomodate multiplicites here
du[s] -= ϕ[i]
end
for s in out
# TODO: accomodate multiplicites here
du[s] += ϕ[i]
end
end
return du
end
return f
end
u₀(m::Petri.Model) = begin
zeros(Float64, length(m.S))
end
u₀(m::Petri.Model, initialS, initialI=1) = begin
u0=zeros(Float64, length(m.S))
u0[1] = initialS
u0[2] = initialI
return u0
end
function savedata(sol::ODESolution, colnames, fname::String)
open(fname, "w") do fp
println(fp, "time, $colnames")
map(tuples(sol)) do (u,t)
print(fp, "$t")
map(u) do x
print(fp, ",$x")
end
println(fp, "")
end
end
end
model(c::PetriCospan) = left(c.f).d[1].model
X = FinSet(1)
# Fseir = compose(exposure(X,X,X),otimes(spontaneous(X,X),id(X)),mmerge(X),mcopy(X),otimes(id(X),spontaneous(X,X)))
# states are [S, I, E, R]
# m_seir = left(Fseir.f).d[1]
Psir = PetriModel(
Petri.Model(1:3,[
([1,2],[2,2]), # exposure
([2],[3]), # recovery
], missing, missing))
m = Psir.model
u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[2] = 1
β = [10/sum(u0), 1/5]
tspan = (0,100.0)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob, alg=Tsit5())
savedata(sol, "S,I,R", "sirdata.csv")
Psird = PetriModel(
Petri.Model(1:3,[
([1,2],[2,2]), # exposure
([2],[3]), # recovery
([2],[4]), # death
], missing, missing))
m = Psir.model
u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[2] = 1
β = [10/sum(u0), 1/5, 1/10]
tspan = (0,100.0)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob, alg=Tsit5())
savedata(sol, "S,I,R,D", "sirddata.csv")
Pseir = PetriModel(
Petri.Model(1:5,[
([1,2],[3,2]), # exposure
([3],[2]), # onset
([2],[4]), # recovery
], missing, missing))
m = Pseir.model
u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[2] = 1
β = [10/sum(u0), 1/2, 1/5]
tspan = (0,100.0)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob, alg=Tsit5())
savedata(sol, "S,I,E,R", "seirdata.csv")
Pseird = PetriModel(
Petri.Model(1:5,[
([1,2],[3,2]), # exposure
([3],[2]), # onset
([2],[4]), # recovery
([2],[5]), # death
], missing, missing))
m = Pseird.model
u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[2] = 1
seirdparams() = begin
βseird = [10/sum(u0), 1/2, 1/5, 1/16]
end
tspan = (0,100.0)
prob = ODEProblem(fluxes(m), u0, tspan, seirdparams())
sol = OrdinaryDiffEq.solve(prob, alg=Tsit5())
savedata(sol, "S,I,E,R,D", "seirddata.csv")

6984
examples/petri/malaria.ipynb

File diff suppressed because one or more lines are too long

1
examples/petri/malaria.jl

@ -1,3 +1,4 @@ @@ -1,3 +1,4 @@
# -*- coding: utf-8 -*-
module Malaria
using Petri
using MacroTools

6
examples/petri/petri_demo.jl

@ -134,7 +134,7 @@ p, p2, p3 = main() @@ -134,7 +134,7 @@ p, p2, p3 = main()
Petri.solve(p)
@time Petri.solve(p)
mf = Petri.eval(Petri.funckit(p))
mf = Petri.evaluate(Petri.funckit(p))
pf = Petri.Problem(mf, SIRState(100, 1, 0, 0.5, 0.15, 0.05), 150)
Petri.solve(pf)
@time Petri.solve(pf)
@ -145,7 +145,7 @@ Petri.solve(pf) @@ -145,7 +145,7 @@ Petri.solve(pf)
Petri.solve(p2)
@time Petri.solve(p2)
mf2 = Petri.eval(Petri.funckit(p2))
mf2 = Petri.evaluate(Petri.funckit(p2))
pf2 = Petri.Problem(mf2, SEIRState(100, 1, 0, 0.5, 0.15, 0.05, 0, 0.12), 150)
Petri.solve(pf2)
@time Petri.solve(pf2)
@ -156,7 +156,7 @@ Petri.solve(pf2) @@ -156,7 +156,7 @@ Petri.solve(pf2)
Petri.solve(p3)
@time Petri.solve(p3)
mf3 = Petri.eval(Petri.funckit(p3))
mf3 = Petri.evaluate(Petri.funckit(p3))
pf3 = Petri.Problem(mf3, SEIRDState(100, 1, 0, 0.5, 0.15, 0.05, 0, 0.12, 0, 0.1), 150)
Petri.solve(pf3)
@time Petri.solve(pf3)

92
examples/petri/rewrite.jl

@ -22,11 +22,12 @@ import Base: ∘ @@ -22,11 +22,12 @@ import Base: ∘
(a::WiringDiagram, b::WiringDiagram) = compose(b, a)
(a,b) = b a
S, E, I, R, D= Ob(FreeSymmetricMonoidalCategory, :S, :E, :I, :R, :D)
wd(args...) = to_wiring_diagram(args...)
inf = WiringDiagram(Hom(:infection, S I, I⊗I))
expo = WiringDiagram(Hom(:exposure, S⊗I, E⊗I))
rec = WiringDiagram(Hom(:recovery, I, R))
wan = WiringDiagram(Hom(:waning, R, S))
inf = wd(Hom(:infection, S I, I⊗I))
expo = wd(Hom(:exposure, S⊗I, E⊗I))
rec = wd(Hom(:recovery, I, R))
wan = wd(Hom(:waning, R, S))
sir_wire = inf (rec rec)
@ -41,37 +42,41 @@ dump(S) @@ -41,37 +42,41 @@ dump(S)
# -
sir = Petri.Model([S, I, R], [(I, R), (S+I, 2I)])
sir_model = model(PetriModel, sir)
ir = Petri.Model([S, I, R], [(I, R)])
ir_model = model(PetriModel, ir)
seir = Petri.Model([S, I, R], [(I, R), (S+I, I+E), (E, I)])
seir_model = model(PetriModel, seir)
# +
irs = Petri.Model([S, I, R], [(I, R), (R, S)])
irs_model = model(PetriModel, irs)
dump(sir)
dump(irs)
# -
sirs = Petri.Model([S, I, R], [(I, R), (S+I, 2*I), (R, S)])
rule = PetriModels.Span(sir, ir, seir)
sirs′ = PetriModels.pushout(irs, sir)
@test sirs′.Δ == sirs.Δ
seirs = PetriModels.solve(PetriModels.DPOProblem(rule, irs))
@test all(Set(seirs.Δ) .== Set([(S+I, I+E),
(E, I),
(I, R),
(R, S)]))
l = sir
c = ir
r = seir
c′ = irs
sirs_model = model(PetriModel, sirs)
rule = PetriModels.PetriSpan(sir_model, ir_model, seir_model)
sirs_model′ = PetriModels.pushout(irs_model, sir_model)
@test sirs_model′.model.Δ == sirs_model.model.Δ
seirs = PetriModels.solve(PetriModels.DPOProblem(rule, irs_model))
@test all(Set(seirs.model.Δ) .== Set([(S+I, I+E),
(E, I),
(I, R),
(R, S)]))
l = sir_model
c = ir_model
r = seir_model
c′ = irs_model
l′ = PetriModels.pushout(l, c′)
@test l′.Δ == sirs.Δ
@test PetriModels.dropdown(l,c,l′).Δ == c′.Δ
@test PetriModels.pushout(r, c′).Δ == seirs.Δ
@test l′.model.Δ == sirs_model.model.Δ
@test PetriModels.dropdown(l,c,l′).model.Δ == c′.model.Δ
@test PetriModels.pushout(r, c′).model.Δ == seirs.model.Δ
function Δ(m::Petri.Model, ctx=:state)
function updateblock(exp, sym)
@ -173,7 +178,7 @@ function Δ(m::Petri.Model, ctx=:state) @@ -173,7 +178,7 @@ function Δ(m::Petri.Model, ctx=:state)
end
@show "SIR"
Δ(l, :state)
Δ(l.model, :state)
# @show "IR"
# funckit(c, :state)
@show "SEIR"
@ -229,35 +234,35 @@ function test_1() @@ -229,35 +234,35 @@ function test_1()
no_transitions = Tuple{Operation, Operation}[]
@variables A, B, C, D
states = [A, B, C, D]
l = Petri.Model(states, [(A, B)])
c = Petri.Model(states, no_transitions)
r = Petri.Model(states, [(A, B + C)])
rule = PetriModels.Span(l, c, r)
c′ = Petri.Model(states, [(B, A)])
l = model(PetriModel, Petri.Model(states, [(A, B)]))
c = model(PetriModel, Petri.Model(states, no_transitions))
r = model(PetriModel, Petri.Model(states, [(A, B + C)]))
rule = PetriModels.PetriSpan(l, c, r)
c′ = model(PetriModel, Petri.Model(states, [(B, A)]))
r′ = PetriModels.solve(PetriModels.DPOProblem(rule, c′))
@test r′.Δ == [(A, B+C), (B, A)]
@test r′.model.Δ == [(A, B+C), (B, A)]
l′ = PetriModels.pushout(l, c′)
@test l′.Δ == [(A, B), (B, A)]
@test PetriModels.dropdown(l,c,l′).Δ == [(B, A)]
@test l′.model.Δ == [(A, B), (B, A)]
@test PetriModels.dropdown(l,c,l′).model.Δ == [(B, A)]
end
function test_2()
no_transitions = Tuple{Operation, Operation}[]
@variables A, B, C, D
states = [A, B, C, D]
l = Petri.Model(states, [(A, B), (B,C)])
c = Petri.Model(states, no_transitions)
r = Petri.Model(states, [(A, C)])
rule = PetriModels.Span(l, c, r)
c′ = Petri.Model(states, [(C, D)])
l = model(PetriModel, Petri.Model(states, [(A, B), (B,C)]))
c = model(PetriModel, Petri.Model(states, no_transitions))
r = model(PetriModel, Petri.Model(states, [(A, C)]))
rule = PetriModels.PetriSpan(l, c, r)
c′ = model(PetriModel, Petri.Model(states, [(C, D)]))
r′ = PetriModels.solve(PetriModels.DPOProblem(rule, c′))
@test r′.Δ == [(A, C), (C, D)]
@test r′.model.Δ == [(A, C), (C, D)]
l′ = PetriModels.pushout(l, c′)
@test l′.Δ == [(A, B), (B, C), (C, D)]
@test PetriModels.dropdown(l,c,l′).Δ == [(C, D)]
@test PetriModels.pushout(rule.r, c′).Δ == [(A, C), (C, D)]
@test l′.model.Δ == [(A, B), (B, C), (C, D)]
@test PetriModels.dropdown(l,c,l′).model.Δ == [(C, D)]
@test PetriModels.pushout(rule.r, c′).model.Δ == [(A, C), (C, D)]
end
test_1()
test_2()
@ -322,7 +327,7 @@ end @@ -322,7 +327,7 @@ end
sir′ = funckit(sir)
m = sir′
m′ = eval(m)
m′ = evaluate(m)
@show m′
@show typeof(m′)
p = Petri.Problem(m′, ParamSIR(100, 1, 0, [0.15, 0.55/101]), 250)
@ -330,7 +335,7 @@ soln = Petri.solve(p) @@ -330,7 +335,7 @@ soln = Petri.solve(p)
@test soln.S <= 10
@test soln.S + soln.I + soln.R == 101
m′ = eval(funckit(sirs))
m′ = evaluate(funckit(sirs))
@show m′
@show typeof(m′)
p = Petri.Problem(m′, ParamSIR(100, 1, 0, [0.15, 0.55/101, 0.15]), 250)
@ -339,6 +344,3 @@ sirs_soln = Petri.solve(p) @@ -339,6 +344,3 @@ sirs_soln = Petri.solve(p)
@test sirs_soln.S + sirs_soln.I + sirs_soln.R == 101
@show soln
@show sirs_soln
m′ = Petri.eval(funckit(seirs))
# @code_native m′.Δ[1](ParamSIR(100, 1, 0, [ 0.15, 0.55/101, 0.15, 0.1 ]))

73
examples/petri/wiring_ode.jl

@ -10,21 +10,22 @@ import Catlab.Doctrines.⊗ @@ -10,21 +10,22 @@ import Catlab.Doctrines.⊗
import Base:
(a::WiringDiagram, b::WiringDiagram) = compose(b, a)
(a,b) = b a
?WiringDiagram
wd(args...) = to_wiring_diagram(args...)
# +
# Generators
A, B, C, D = Ob(FreeSymmetricMonoidalCategory, :A, :B, :C, :D)
f = WiringDiagram(Hom(:f,A,B))
g = WiringDiagram(Hom(:g,B,A))
f = wd(Hom(:f,A,B))
g = wd(Hom(:g,B,A))
boxes(f)
@test nboxes(f) == 1
@test boxes(f) == [ Box(Hom(:f,A,B)) ]
@test boxes(f) == AbstractBox[ Box(:f,[:A],[:B]) ]
@test nwires(f) == 2
@test Ports([:A]) == Ports([:A])
@test WiringDiagram(Hom(:f,A,B)) == WiringDiagram(Hom(:f,A,B))
@test wd(Hom(:f,A,B)) == wd(Hom(:f,A,B))