CoCalc Shared Filesprojects / Project 2.ipynb
Authors: Owen Dray, Eamon Ito-Fisher
Views : 46
In [1]:
# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import functions from the modsim.py module
from modsim import *

In [2]:
def make_system():
"""Make a system object for the HIV model.

returns: System object
"""
gamma=1.36
tau=0.2
mu=1.36e-3
beta=0.00027
ro=0.1
alpha=3.6e-2
delta=0.33
pi=100
sigma=2
t_0 = 0

t_end = 120
dT=0.06

return System(gamm=gamma, tau=tau, mu=mu,
beta=beta, ro=ro, alpha=alpha,delta=delta, pi=pi, sigma=sigma,t_0=t_0,t_end=t_end,dT=dT)

In [3]:
def update(state, time, system):
R, L, E, V = state
gamma, tau, mu, beta, ro, alpha, delta, pi, sigma, t_0, t_end, dT = system

dRdT = gamma * tau - mu * R - beta * R * V
dLdT = ro * beta * R * V - mu * L - alpha * L
dEdT = (1 - ro) * beta * R * V + alpha * L - delta * E
dVdT = pi * E - sigma * V
R += dRdT * dT
L += dLdT * dT
E += dEdT * dT
V += dVdT * dT

return (State(R=R, L=L, E=E, V=V))

In [4]:
def run_simulation(system, update_func):
"""Runs a simulation of the system.

Add a TimeFrame to the System: results

system: System object
"""
init = State(R=200, L=0, E=0, V=4e-7)
gamma, tau, mu, beta, ro, alpha, delta, pi, sigma, t_0, t_end, dT = system

frame = TimeFrame(columns=init.index)
frame.row[t_0] = init
ts = linrange(t_0, t_end, dT)

for t in ts:
frame.row[t+dT] = update_func(frame.row[t], t, system)

return frame


In [5]:
sys = make_system()
frame = run_simulation(sys, update)

R L E V
0.00 200.000000 0.000000e+00 0.000000e+00 4.000000e-07
0.06 200.000000 1.296000e-10 1.166400e-09 3.520000e-07
0.12 200.000000 2.433575e-10 2.170017e-09 3.167584e-07
0.18 200.000000 3.454417e-10 3.051244e-09 2.917675e-07
0.24 200.000000 4.392000e-10 3.842370e-09 2.750629e-07
... ... ... ... ...
119.76 17.993435 6.626104e-01 2.389403e-01 1.201070e+01
119.82 18.004785 6.614752e-01 2.387915e-01 1.200306e+01
119.88 18.016135 6.603425e-01 2.386431e-01 1.199544e+01
119.94 18.027484 6.592124e-01 2.384952e-01 1.198784e+01
120.00 18.038832 6.580848e-01 2.383478e-01 1.198027e+01

2001 rows × 4 columns

In [15]:
# plot(frame["V"])
# plot(frame["R"])
# plot(frame["L"])
# plot(frame["E"])
fig = plt.figure(figsize=(8, 10))

line = ax.plot(frame["V"],label='V')
ax.set_yscale('log')
ax.set_ylim([0.1,10e3])
legend()
line = ax.plot(frame["R"],label='R')
legend()
line = ax.plot(frame["L"],label='L')
legend()
line = ax.plot(frame["E"],label='E')
title
legend()

In [ ]:


In [ ]:


In [ ]:


In [ ]: