Open Circular Port Graph Examples using AlgebraicDynamics.DWDDynam
using AlgebraicDynamics.CPortGraphDynam
using AlgebraicDynamics.CPortGraphDynam: draw, barbell, gridpath, grid, meshpath

using Catlab
using Catlab.WiringDiagrams
using Catlab.WiringDiagrams.CPortGraphs
using Catlab.Theories
using Catlab.CategoricalAlgebra

using OrdinaryDiffEq
using Plots, Plots.PlotMeasures

using PrettyTables
WARNING: using Plots.grid in module Main conflicts with an existing identifier.

SIR Epidemiology Model

An SIR epidemiology model has three types of people: susceptible, infected, and recovered. When a susceptible person interacts with infected person, the susceptible person also becomes infected. Over time infected people recover. Transition rates determine the how frequently susceptible people come into contact with infected people and how fast infected people recover. The system evolves according to the law of mass action.

In a multi-city SIR model, each city has susceptible, infected, and recovered populations. To see the spread of the disease we will consider both susceptible and infected people moving between cities. To define a multi-city SIR model, we can compose multiple single-city SIR models using the composition syntax of open CPGs. The composition pattern will consist of three boxes each of which will be filled by a single-city SIR model. Ports expose the susceptible and infected populations of each city. One set of wires connect the susceptible and infected populations of cities 1 and 2. A second set of wires connect the susceptible and infected popuation of cities 2 and 3.

# Define the composition pattern
d₀ = OpenCPortGraph()
d₁ = barbell(2)
F = ACSetTransformation((Box=,), d₀, d₁)
G = ACSetTransformation((Box=,), d₀, d₁)
d₂ = apex(pushout(F,G))

# Define the primitive systems
β, μ, α₁, α₂ = 0.4, 0.4, 0.01, 0.01

sirfuncb = (u,x,p,t)->[-β*u*u - α₁*(u-x), # Ṡ
β*u*u - μ*u - α₂*(u-x), #İ
μ*u # Ṙ
]
sirfuncm = (u,x,p,t)->[-β*u*u - α₁*(u-(x+x)/2),
β*u*u - μ*u - α₂*(u-(x+x)/2),
μ*u
]

boundary  = ContinuousMachine{Float64}(2,3,sirfuncb, (u,p,t)->u[1:2])
middle    = ContinuousMachine{Float64}(4,3, sirfuncm, (u,p,t)->u[[1,2,1,2]])

# Compose
threecity = oapply(d₂, [boundary,middle,boundary])

First, we will approximate the solution to the three city SIR model using Euler's method. The initial condition has 100 susceptible people in each city a single infected person in the first city. We show the infected populations in each city over time.

u0 = [100,1,0,100,0,0,100,0,0.0]

h = 0.01
tspan = (0.0, 1.0)
threecity_approx = euler_approx(threecity, h)
prob = DiscreteProblem(threecity_approx, u0, [], tspan, nothing)
sol = solve(prob, FunctionMap(); dt = h)

map(sol) do u
return (i1=u, i2=u, i3=u)
end |> pretty_table
┌─────────┬─────────────┬─────────────┐
│      i1 │          i2 │          i3 │
│ Float64 │     Float64 │     Float64 │
├─────────┼─────────────┼─────────────┤
│     1.0 │         0.0 │         0.0 │
│  1.3959 │      5.0e-5 │         0.0 │
│  1.9463 │  0.00013959 │      5.0e-9 │
│  2.7094 │ 0.000292169 │  2.09385e-8 │
│ 3.76334 │ 0.000543308 │  5.84449e-8 │
│ 5.21121 │  0.00094657 │  1.35914e-7 │
│ 7.18564 │  0.00158188 │   2.8438e-7 │
│ 9.85078 │  0.00256741 │  5.55153e-7 │
│ 13.3983 │  0.00407637 │  1.03168e-6 │
│ 18.0311 │  0.00636006 │  1.84776e-6 │
│ 23.9277 │   0.0097794 │  3.21529e-6 │
│ 31.1815 │   0.0148471 │  5.46617e-6 │
│ 39.7176 │   0.0222833 │  9.11493e-6 │
│ 49.2146 │   0.0330891 │  1.49519e-5 │
│ 59.0815 │   0.0486455 │  2.41802e-5 │
│ 68.5483 │   0.0708489 │  3.86177e-5 │
│ 76.8715 │    0.102305 │  6.09913e-5 │
│ 83.5617 │     0.14661 │  9.53683e-5 │
│ 88.4952 │    0.208744 │ 0.000147786 │
│ 91.8553 │    0.295635 │ 0.000227168 │
│ 93.9784 │    0.416918 │ 0.000346667 │
│ 95.2144 │    0.585974 │ 0.000525604 │
│ 95.8529 │    0.821336 │ 0.000792287 │
│ 96.1049 │     1.14857 │  0.00118809 │
│ 96.1134 │     1.60277 │   0.0017733 │
│ 95.9711 │     2.23173 │  0.00263561 │
│ 95.7361 │     3.09995 │   0.0039022 │
│  95.445 │     4.29318 │  0.00575701 │
│ 95.1199 │     5.92314 │   0.0084654 │
│ 94.7747 │     8.13108 │   0.0124089 │
│ 94.4178 │     11.0877 │    0.018134 │
│ 94.0544 │      14.985 │   0.0264207 │
│ 93.6876 │     20.0143 │   0.0383764 │
│ 93.3194 │     26.3226 │   0.0555651 │
│ 92.9511 │     33.9449 │   0.0801831 │
│ 92.5834 │     42.7238 │    0.115296 │
│ 92.2167 │     52.2483 │     0.16516 │
│ 91.8514 │     61.8688 │    0.235662 │
│ 91.4875 │     70.8272 │    0.334923 │
│ 91.1249 │     78.4739 │    0.474149 │
│ 90.7636 │     84.4564 │    0.668803 │
│ 90.4034 │     88.7674 │     0.94022 │
│ 90.0444 │     91.6474 │     1.31776 │
│ 89.6865 │     93.4347 │      1.8416 │
│ 89.3296 │     94.4517 │     2.56633 │
│ 88.9739 │     94.9542 │     3.56516 │
│ 88.6194 │      95.125 │     4.93473 │
│ 88.2662 │     95.0865 │     6.79946 │
│ 87.9142 │      94.918 │     9.31398 │
│ 87.5635 │     94.6694 │     12.6602 │
│ 87.2141 │     94.3719 │     17.0337 │
│ 86.8661 │     94.0449 │     22.6133 │
│ 86.5195 │     93.7002 │     29.5066 │
│ 86.1742 │     93.3454 │      37.674 │
│ 85.8302 │     92.9848 │     46.8507 │
│ 85.4877 │     92.6213 │     56.5124 │
│ 85.1464 │     92.2566 │      65.938 │
│ 84.8066 │     91.8916 │     74.3883 │
│ 84.4681 │      91.527 │     81.3273 │
│ 84.1309 │     91.1632 │     86.5583 │
│ 83.7951 │     90.8002 │     90.2011 │
│ 83.4606 │     90.4383 │     92.5573 │
│ 83.1275 │     90.0774 │     93.9686 │
│ 82.7957 │     89.7178 │     94.7313 │
│ 82.4652 │     89.3594 │     95.0685 │
│  82.136 │     89.0024 │     95.1342 │
│ 81.8081 │     88.6466 │     95.0299 │
│ 81.4816 │     88.2922 │     94.8204 │
│ 81.1564 │     87.9391 │     94.5464 │
│ 80.8324 │     87.5874 │     94.2331 │
│ 80.5098 │     87.2372 │     93.8962 │
│ 80.1884 │     86.8882 │     93.5451 │
│ 79.8683 │     86.5407 │      93.186 │
│ 79.5495 │     86.1946 │     92.8223 │
│  79.232 │     85.8498 │     92.4565 │
│ 78.9157 │     85.5064 │     92.0898 │
│ 78.6007 │     85.1644 │     91.7232 │
│ 78.2869 │     84.8237 │     91.3572 │
│ 77.9745 │     84.4844 │      90.992 │
│ 77.6632 │     84.1465 │      90.628 │
│ 77.3532 │     83.8099 │     90.2652 │
│ 77.0444 │     83.4747 │     89.9038 │
│ 76.7369 │     83.1408 │     89.5437 │
│ 76.4306 │     82.8082 │      89.185 │
│ 76.1255 │      82.477 │     88.8276 │
│ 75.8216 │     82.1471 │     88.4717 │
│  75.519 │     81.8185 │     88.1172 │
│ 75.2175 │     81.4912 │     87.7642 │
│ 74.9173 │     81.1653 │     87.4125 │
│ 74.6183 │     80.8406 │     87.0622 │
│ 74.3204 │     80.5172 │     86.7134 │
│ 74.0237 │     80.1952 │     86.3659 │
│ 73.7283 │     79.8744 │     86.0198 │
│  73.434 │     79.5549 │     85.6751 │
│ 73.1408 │     79.2367 │     85.3318 │
│ 72.8489 │     78.9197 │     84.9899 │
│ 72.5581 │      78.604 │     84.6493 │
│ 72.2685 │     78.2896 │     84.3101 │
│   71.98 │     77.9765 │     83.9723 │
│ 71.6927 │     77.6646 │     83.6358 │
│ 71.4065 │     77.3539 │     83.3006 │
└─────────┴─────────────┴─────────────┘

Next, we will solve the continuous system and plot the results. Over time the infected populations increase and the susceptible populations decrease. The delays in the plots illustrate how the disease spreads from city 1 to city 2 and then from city 2 to city 3.

# Solve and plot
prob = ODEProblem(threecity, u0, tspan)
sol = solve(prob, Tsit5(); dtmax = 0.01)

plot(sol, lw=2, title = "SIR Epidemiology Model", bottom_margin=10mm, left_margin=10mm,
label=["S" "I" "R"])

Cellular automata

Circular port graphs are particularly useful for modeling systems where the composition pattern is given by a grid and where the dynamics are repetative. In the case of cellular automata the composition pattern is a row of $n$ cells each of which is connected to its two neighbors. The primitive systems are identical machines whose discrete dynamics are a specified rule. See here for a complete set of rules and the patterns they generate.

function Rule(k::Int)
(left_neighbor, x, right_neighbor) ->
Bool(digits(k, base=2, pad=8)[1 + right_neighbor + 2*x + 4*left_neighbor])
end

# Define the composition pattern
n = 100
row = apex(gridpath(n, 1))

# Define the primitive system which will be repeated for each cell
rule = DiscreteMachine{Bool}(2, 1, 2,
(u, x, p, t)->Rule(p)(x, u, x),
(u,p,t)->[u, u])

# Compose
automaton = oapply(row, rule)

# Solve and plot
u0 = zeros(Int, n); u0[Int(n/2)] = 1

rule_number = 126
tspan = (0.0, 100.0)
prob = DiscreteProblem(automaton, u0, [0,0], tspan, rule_number)
sol = solve(prob, FunctionMap())
spy(transpose(reduce(hcat, sol.u)))