Demo and cookbook
[1]:
import comod
import numpy as np
Defining a model
Models are defined enumerating the states, parameters and transition rules.
Each of this rules is given by the origin compartment, the destination compartment, and the proportionality coeffient with describes the flow per individual in the origin to the destination.
[2]:
model = comod.Model(
["S", "I", "R"], # States
["beta", "gamma"], # Coefficients
[ # Rules in the form (origin, destination, coefficient)
("S", "I", "beta I / N"), # N is a special state with the total population
("I", "R", "gamma"),
],
# Special state names can be set with the following options:
sum_state="N", # Total population. Can be used in the coefficients, but not as origin/destination.
nihil_state="!", # The nothingness (?) from which one is born and to which one dies.
# Can be used as origin/destination, but not in the cofficients.
)
Plotting the model
Plotting with tikz
Plotting with tikz might be convenient for the generation of high quaility plots of small models.
[3]:
model.plot_tikz(layout={"S": (1, 0), "I": (2, 0), "R": (3, 0)})
[3]:
The generating tex file can be exported for further tweaking.
[4]:
model.plot_tikz("export_demo.tex", layout={"S": (1, 0), "I": (2, 0), "R": (3, 0)})
[5]:
!cat export_demo.tex
\documentclass{standalone}
\usepackage{tikz-network}
\begin{document}
\begin{tikzpicture}
\clip (0,0) rectangle (6,6);
\Vertex[x=0.350,y=3.000,color=blue,label=$S$]{S}
\Vertex[x=3.000,y=3.000,color=blue,label=$I$]{I}
\Vertex[x=5.650,y=3.000,color=blue,label=$R$]{R}
\Edge[,label=$\beta I / N$,Direct](S)(I)
\Edge[,label=$\gamma$,Direct](I)(R)
\end{tikzpicture}
\end{document}
Plotting with igraph
Plotting python igraph might be convenient for larger graphs due to its layout capabilities.
[6]:
model.plot_igraph(layout="circle", bbox=(300, 300 / 1.618))
[6]:
Output the differential equations in LaTeX
[7]:
latex = model.to_latex()
print(latex)
\begin{array}{lcl} \dot{S} &=& -\beta I / N S\\
\dot{I} &=& \beta I / N S - \gamma I\\
\dot{R} &=& \gamma I
\end{array}
[8]:
from IPython.display import display, Math
display(Math(latex))
Solve and plot
Basic solution
The model can be numerically solved given the values of the parameters and the initial conditions.
[9]:
t = np.linspace(0, 150, 150)
S, I, R = model.solve(
(999, 1, 0), # Initial state
[0.3, 0.1], # Coefficient values
t, # Time mesh
method="RK45", # see scipy.integrate.solve_ivp
)
[10]:
import matplotlib.pyplot as plt
plt.plot(t, S, "b", label="S")
plt.plot(t, I, "r", label="I")
plt.plot(t, R, "g", label="R")
plt.xlabel("Time [days]")
plt.ylabel("Number of individuals")
plt.legend()
plt.show()
Dynamic solution using ipywidgets
Interactive visualization using ipywidgets. Only available when running as a notebook.
[11]:
import ipywidgets
[12]:
@ipywidgets.interact
def plot(
beta=ipywidgets.FloatSlider(min=0.0, max=2, value=0.3, description="$\\beta$"),
gamma=ipywidgets.FloatSlider(
min=0.0, max=0.3, value=0.1, step=0.01, description="$\\gamma$"
),
):
t = np.linspace(0, 150, 150)
S, I, R = model.solve((999, 1, 0), [beta, gamma], t)
plt.plot(t, S, "b", label="S")
plt.plot(t, I, "r", label="I")
plt.plot(t, R, "g", label="R")
plt.xlabel("Time [days]")
plt.ylabel("Number of individuals")
plt.legend()
plt.show()
Time-dependent parameter solution
The model can also be solved if the parameters are functions of time
[13]:
def beta_time(t):
return max(0.1, 0.3 - t / 30 * 0.1)
[14]:
plt.plot([beta_time(x) for x in t])
plt.show()
[15]:
solution_time = model.solve_time(
(999, 1, 0), # Initial state
[beta_time, lambda x: 0.1], # Coefficient values
t, # Time mesh
)
[16]:
import matplotlib.pyplot as plt
plt.plot(t, solution_time[0], "b", label="S")
plt.plot(t, solution_time[1], "r", label="I")
plt.plot(t, solution_time[2], "g", label="R")
plt.xlabel("Time [days]")
plt.ylabel("Number of individuals")
plt.legend()
plt.show()
Best fit
Use the previously generated data to best-fit a model.
[17]:
fit_pars = model.best_fit(
np.asarray([S, I, R]), t, [0.2, 0.2] # Existing data # Time mesh # Initial guess
).x
print("Fitted parameters:")
for par, value in zip(model.parameters, fit_pars):
print("%s: %.3f" % (par, value))
Fitted parameters:
beta: 0.298
gamma: 0.101
Best fit in time windows
Fits can be performed using sliding windows.
First, let’s try with the constant model:
[18]:
results = model.best_sliding_fit(
np.asarray([S, I, R]), # Existing data
t, # Time mesh
[0.2, 0.2], # Initial guess
window_size=20, # Window size
step_size=25, # Number of time steps between windows
target="y", # The target to minimize can be "y" (values) or "dy" (their derivatives)
)
results
[18]:
beta | gamma | |
---|---|---|
9.563758 | 0.300000 | 0.100001 |
34.731544 | 0.300037 | 0.100028 |
59.899329 | 0.301793 | 0.100076 |
85.067114 | 0.298463 | 0.099781 |
110.234899 | 0.299341 | 0.099746 |
135.402685 | 0.299556 | 0.099773 |
[19]:
plt.plot(results)
plt.legend(["$\\beta$", "$\\gamma$"])
plt.show()
Now let’s try the time-dependent coefficient model
[20]:
results = model.best_sliding_fit(
solution_time, # Existing data
t, # Time mesh
[0.2, 0.2], # Initial guess
window_size=10, # Window size
step_size=15, # Number of time steps between windows
target="y", # The target to minimize can be "y" (values) or "dy" (their derivatives)
time_criterion="mean", # Determines how a time is assigned to each window in the solution
)
results
[20]:
beta | gamma | |
---|---|---|
4.530201 | 0.288520 | 0.102111 |
19.630872 | 0.238503 | 0.101901 |
34.731544 | 0.188711 | 0.101453 |
49.832215 | 0.138761 | 0.100800 |
64.932886 | 0.098525 | 0.100100 |
80.033557 | 0.099614 | 0.099815 |
95.134228 | 0.099563 | 0.099749 |
110.234899 | 0.100074 | 0.100113 |
125.335570 | 0.100256 | 0.100254 |
140.436242 | 0.099812 | 0.099812 |
[21]:
plt.plot(results)
plt.plot([beta_time(x) for x in t], "C0--")
plt.legend(["$\\beta$", "$\\gamma$", "$\\beta$ (theory)"])
plt.show()
Symbolic calculation of fixed points
A sympy represenation of the model can automatically be obtained if the package is installed:
[22]:
model.get_sympy()
[22]:
An automatic calculation of fixed points can be performed for simple models.
[23]:
for point, jacobian in model.get_fixed_points():
print("Fixed point:")
display(point)
print("Jacobian:")
display(jacobian)
Fixed point:
(S, 0.0, R)
Jacobian:
[24]:
jacobian.eigenvals()
[24]:
{-(R*gamma - S*beta + S*gamma)/(R + S): 1, 0: 2}
Stability analysis is left as an exercise for the reader.