GNARWHAL

947 days ago by ethan.kerzner

This code is modified from http://wiki.deductivethinking.com/wiki/Python_Programs_for_Modelling_Infectious_Diseases_book.  The license permits copying, modifying, and distributing changes.  See the bottom of this worksheet for the license.

import scipy.integrate import matplotlib.pyplot as plt import numpy N = 2100000 time_step=.05 time_end=100 betaK=1 # contacts of iK with sU betaU=22600# contacts of iU with sU omegaK=2 # contacts of iU with sE omegaU=1 # contacts of iK with sE phi=.2154/4 # rate of education gamma=.0336 # rate of becoming aware of infection #initial starting populations iU0 = 200000*.69 sU0 = 1900000*.75 sE0 = 1900000*.25 iK0 = 200000*.31 BRiU = .0241*.1 BRsU = .0241*.9 BRsE = .0241*.9 BRiK = .0241*.1 DRiU = .63*.0222 DRsU = .37*.0222 DRsE = .37*.0222 DRiK = .63*.0222 # First, define the function that takes the previous time populations # and returns the new populations def differential_eqs(previous_state,t): sU,sE,iU,iK = previous_state # susceptible uneducated sU_new = - (betaU / N) * (iU) * (sU) - phi * sU + (BRsU * sU + BRsE * sE - DRsU * sU) #susceptible educated sE_new = phi * sU - (DRsE*sE) #infected unknown iU_new = (betaU / N) * (iU) * (sU) - (phi + gamma) * iU + (BRiU - DRiU)* iU #infected known iK_new = (phi + gamma) * iU + (BRiK - DRiK) return sU_new, sE_new, iU_new, iK_new # Now solve the differential equations to get the actual populations start_populations = (sU0, sE0, iU0, iK0) time_range=[0.0..time_end, step=time_step] populations = scipy.integrate.odeint(differential_eqs,start_populations,time_range) # Extract the data for each group of people susceptible_U = populations[:,0] # first column susceptible_E = populations[:,1] # second column infected_U = populations[:,2] # third column infected_K = populations[:,3] #fourth column #Plot the results plt.figure() plt.plot(susceptible_U, '-g', label='Susceptible - Uneducated') plt.plot(susceptible_E, '-b', label='Susceptible - Educated') plt.plot(infected_U, '-r', label='Infected Unknown') plt.plot(infected_K, '-y', label='Infected Known') plt.legend() plt.title('S-I Model with Education') plt.xlabel('Time') plt.ylabel('population') plt.axis([0,70,-2,N]) plt.savefig('test.png') 
       
import scipy.integrate import matplotlib.pyplot as plt import numpy N = 100.0 # the number of people time_step=0.05 time_end=100 betaK=50 # contacts of iK with sU betaU=50 # contacts of iU with sU omegaK=0 # contacts of iU with sE omegaU=0 # contacts of iK with sE phi=1/2 # rate of education gamma=1/9 # rate of becoming aware of infection #initial starting populations iU0 = 1 sU0 = 70 sE0 = 20 iK0 = 9 BRiU = 0 BRsU = 0 BRsE = 0 BRiK = 0 DRiU = .04 DRsU = 0 DRsE = 0 DRiK = .04 # First, define the function that takes the previous time populations # and returns the new populations def differential_eqs(previous_state,t): sU,sE,iU,iK = previous_state # susceptible uneducated sU_new = (-betaK / N) * (iK) * (sU) - (betaU / N) * (iU) * (sU) - phi * sU + (BRsU * sU + BRsE * sE - DRsU * sU) #susceptible educated sE_new = (-omegaK / N) * (iK) * (sE) - (omegaU / N) * (iU) * (sE) + phi * sU - (DRsE*sE) #infected unknown iU_new = (betaK / N) * (iK) * (sU) + (betaU / N) * (iU) * (sU) - (phi + gamma) * iU + (BRiU * iU - DRiU) #infected known iK_new = (omegaK / N) * (iK) * (sE) + (omegaU / N) * (iU) * (sE) + (phi + gamma) * iU + (BRiK * iK - DRiK) return sU_new, sE_new, iU_new, iK_new # Now solve the differential equations to get the actual populations start_populations = (sU0, sE0, iU0, iK0) time_range=[0.0..time_end, step=time_step] populations = scipy.integrate.odeint(differential_eqs,start_populations,time_range) # Extract the data for each group of people susceptible_U = populations[:,0] # first column susceptible_E = populations[:,1] # second column infected_U = populations[:,2] # third column infected_K = populations[:,3] #fourth column total_pop = susceptible_U + susceptible_E + infected_U + infected_K percent_sU = susceptible_U/total_pop*100 percent_sE = susceptible_E/total_pop*100 percent_iU = infected_U/total_pop*100 percent_iK = infected_K/total_pop*100 #print total_pop #print percent_sU, percent_sE, percent_iU, percent_iK #Plot the results plt.figure() plt.plot(percent_sU, '-g', label='Susceptible - Uneducated') plt.plot(percent_sE, '-b', label='Susceptible - Educated') plt.plot(percent_iU, '-r', label='Infected Unknown') plt.plot(percent_iK, '-y', label='Infected Known') plt.legend() plt.title('Population Percentage') plt.xlabel('Time') plt.ylabel('percent population') plt.axis([0,100,0,N]) plt.savefig('test.png') 
       

Let's make this interactive.

 
       

Phase plots

import scipy.integrate import matplotlib.pyplot as plt import numpy N = 100.0 # the number of people time_step=1.0 time_end=70.0 betaK=2 # contacts of iK with sU betaU=6 # contacts of iU with sU omegaK=2 # contacts of iU with sE omegaU=1 # contacts of iK with sE phi=1/4 # rate of education gamma=1/9 # rate of becoming aware of infection #initial starting populations iU0 = 1 sU0 = 70 sE0 = 20 iK0 = 9 BRiU = .06 BRsU = .06 BRsE = .06 BRiK = .02 DRiU = .11 DRsU = .05 DRsE = .04 DRiK = .08 # First, define the function that takes the previous time populations # and returns the new populations def differential_eqs(previous_state,t): "The main set of equations" sU,sE,iU,iK = previous_state # susceptible uneducated sU_new = (-betaK / N) * (iK) * (sU) - (betaU / N) * (iU) * (sU) - phi * sU + (BRsU * sU + BRsE * sE - DRsU * sU) #susceptible educated sE_new = (-omegaK / N) * (iK) * (sE) - (omegaU / N) * (iU) * (sE) + phi * sU - (DRsE*sE) #infected unknown iU_new = (betaK / N) * (iK) * (sU) + (betaU / N) * (iU) * (sU) - (phi + gamma) * iU + (BRiU * iU - DRiU) #infected known iK_new = (omegaK / N) * (iK) * (sE) + (omegaU / N) * (iU) * (sE) + (phi + gamma) * iU + (BRiK * iK - DRiK) return sU_new, sE_new, iU_new, iK_new # Now solve the differential equations to get the actual populations start_populations = (sU0, sE0, iU0, iK0) time_range=[0.0..time_end, step=time_step] populations = scipy.integrate.odeint(differential_eqs,start_populations,time_range) # Extract the data for each group of people susceptible_U = populations[:,0] # first column susceptible_E = populations[:,1] # second column infected_U = populations[:,2] # third column infected_K = populations[:,3] #fourth column total_pop = susceptible_U + susceptible_E + infected_U + infected_K susceptible = susceptible_U + susceptible_E infected = infected_U + infected_K percent_s = susceptible/total_pop*100 percent_i = infected/total_pop*100 #Plot the results plt.figure() plt.plot(percent_s,percent_i) plt.title('Percent SI phase portrait') plt.xlabel('Percent Susceptible (educated and uneducated)') plt.ylabel('Percent Infected (know and unknown)') plt.axis([0,100,0,100]) plt.savefig('test.png') 
       
import scipy.integrate import matplotlib.pyplot as plt import numpy N = 100.0 # the number of people time_step=1.0 time_end=70.0 betaK=2 # contacts of iK with sU betaU=6 # contacts of iU with sU omegaK=2 # contacts of iU with sE omegaU=1 # contacts of iK with sE phi=1/4 # rate of education gamma=1/9 # rate of becoming aware of infection #initial starting populations iU0 = 1 sU0 = 70 sE0 = 20 iK0 = 9 BRiU = .06 BRsU = .06 BRsE = .06 BRiK = .02 DRiU = .11 DRsU = .05 DRsE = .04 DRiK = .08 # First, define the function that takes the previous time populations # and returns the new populations def differential_eqs(previous_state,t): "The main set of equations" sU,sE,iU,iK = previous_state # susceptible uneducated sU_new = (-betaK / N) * (iK) * (sU) - (betaU / N) * (iU) * (sU) - phi * sU + (BRsU * sU + BRsE * sE - DRsU * sU) #susceptible educated sE_new = (-omegaK / N) * (iK) * (sE) - (omegaU / N) * (iU) * (sE) + phi * sU - (DRsE*sE) #infected unknown iU_new = (betaK / N) * (iK) * (sU) + (betaU / N) * (iU) * (sU) - (phi + gamma) * iU + (BRiU * iU - DRiU) #infected known iK_new = (omegaK / N) * (iK) * (sE) + (omegaU / N) * (iU) * (sE) + (phi + gamma) * iU + (BRiK * iK - DRiK) return sU_new, sE_new, iU_new, iK_new # Now solve the differential equations to get the actual populations start_populations = (sU0, sE0, iU0, iK0) time_range=[0.0..time_end, step=time_step] populations = scipy.integrate.odeint(differential_eqs,start_populations,time_range) # Extract the data for each group of people susceptible_U = populations[:,0] # first column susceptible_E = populations[:,1] # second column infected_U = populations[:,2] # third column infected_K = populations[:,3] #fourth column susceptible = susceptible_U + susceptible_E infected = infected_U + infected_K #Plot the results plt.figure() plt.plot(susceptible,infected) plt.title('SI phase portrait') plt.xlabel('Susceptible (educated and uneducated)') plt.ylabel('Infected (know and unknown') plt.axis([0,N,0,N]) plt.savefig('test.png') 
       
import scipy.integrate import matplotlib.pyplot as plt import numpy N = 100.0 # the number of people time_step=1.0 time_end=70.0 betaK=2 # contacts of iK with sU betaU=6 # contacts of iU with sU omegaK=2 # contacts of iU with sE omegaU=1 # contacts of iK with sE phi=1/4 # rate of education gamma=1/9 # rate of becoming aware of infection #initial starting populations iU0 = 1 sU0 = 70 sE0 = 20 iK0 = 9 BRiU = .06 BRsU = .06 BRsE = .06 BRiK = .02 DRiU = .11 DRsU = .05 DRsE = .04 DRiK = .08 # First, define the function that takes the previous time populations # and returns the new populations def differential_eqs(previous_state,t): "The main set of equations" sU,sE,iU,iK = previous_state # susceptible uneducated sU_new = (-betaK / N) * (iK) * (sU) - (betaU / N) * (iU) * (sU) - phi * sU + (BRsU * sU + BRsE * sE - DRsU * sU) #susceptible educated sE_new = (-omegaK / N) * (iK) * (sE) - (omegaU / N) * (iU) * (sE) + phi * sU - (DRsE*sE) #infected unknown iU_new = (betaK / N) * (iK) * (sU) + (betaU / N) * (iU) * (sU) - (phi + gamma) * iU + (BRiU * iU - DRiU) #infected known iK_new = (omegaK / N) * (iK) * (sE) + (omegaU / N) * (iU) * (sE) + (phi + gamma) * iU + (BRiK * iK - DRiK) return sU_new, sE_new, iU_new, iK_new # Now solve the differential equations to get the actual populations start_populations = (sU0, sE0, iU0, iK0) time_range=[0.0..time_end, step=time_step] populations = scipy.integrate.odeint(differential_eqs,start_populations,time_range) # Extract the data for each group of people susceptible_U = populations[:,0] # first column susceptible_E = populations[:,1] # second column infected_U = populations[:,2] # third column infected_K = populations[:,3] #fourth column educated = infected_K + susceptible_E uneducated = susceptible_U + infected_K total_pop = susceptible_U + susceptible_E + infected_U + infected_K percent_educated = educated/total_pop*100 percent_uneducated = uneducated/total_pop*100 #Plot the results plt.figure() plt.plot(percent_educated,percent_uneducated) plt.title('Education phase portrait') plt.xlabel('Educated') plt.ylabel('Uneducated') plt.axis([0,N,0,N]) plt.savefig('test.png') 
       
 
       
 
       

SIR with births and deaths (constant population)

 
       
 
       
 
       

 

##########################################################################
### Copyright (C) <2008> Ilias Soumpasis                                 #
### ilias.soumpasis@deductivethinking.com                                #
### ilias.soumpasis@gmail.com                                             #
###                                                                      #
### This program is free software: you can redistribute it and/or modify #
### it under the terms of the GNU General Public License as published by #
### the Free Software Foundation, version 3.                             #
###                                                                      #
### This program is distributed in the hope that it will be useful,      #
### but WITHOUT ANY WARRANTY; without even the implied warranty of       #
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        #
### GNU General Public License for more details.                         #
###                                                                      #
### You should find a copy of the GNU General Public License at          #
###the Copyrights section or, see http://www.gnu.org/licenses.           #
##########################################################################