Make sure to open in Colab to see the plots!
from scipy.integrate import odeint
from scipy import optimize
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline
#!pip install mpld3
#import mpld3
#mpld3.enable_notebook()
def plotseird(t, S, E, I, R, D=None, L=None, R0=None, Alpha=None, CFR=None):
f, ax = plt.subplots(1,1,figsize=(10,4))
ax.plot(t, S, 'b', alpha=0.7, linewidth=2, label='Susceptible')
ax.plot(t, E, 'y', alpha=0.7, linewidth=2, label='Exposed')
ax.plot(t, I, 'r', alpha=0.7, linewidth=2, label='Infected')
ax.plot(t, R, 'g', alpha=0.7, linewidth=2, label='Recovered')
if D is not None:
ax.plot(t, D, 'k', alpha=0.7, linewidth=2, label='Dead')
ax.plot(t, S+E+I+R+D, 'c--', alpha=0.7, linewidth=2, label='Total')
else:
ax.plot(t, S+E+I+R, 'c--', alpha=0.7, linewidth=2, label='Total')
ax.set_xlabel('Time (days)')
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend(borderpad=2.0)
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
ax.spines[spine].set_visible(False)
if L is not None:
plt.title("Lockdown after {} days".format(L))
plt.show();
if R0 is not None or CFR is not None:
f = plt.figure(figsize=(12,4))
if R0 is not None:
# sp1
ax1 = f.add_subplot(121)
ax1.plot(t, R0, 'b--', alpha=0.7, linewidth=2, label='R_0')
ax1.set_xlabel('Time (days)')
ax1.title.set_text('R_0 over time')
# ax.set_ylabel('Number (1000s)')
# ax.set_ylim(0,1.2)
ax1.yaxis.set_tick_params(length=0)
ax1.xaxis.set_tick_params(length=0)
ax1.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax1.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
ax.spines[spine].set_visible(False)
if Alpha is not None:
# sp2
ax2 = f.add_subplot(122)
ax2.plot(t, Alpha, 'r--', alpha=0.7, linewidth=2, label='alpha')
ax2.set_xlabel('Time (days)')
ax2.title.set_text('fatality rate over time')
# ax.set_ylabel('Number (1000s)')
# ax.set_ylim(0,1.2)
ax2.yaxis.set_tick_params(length=0)
ax2.xaxis.set_tick_params(length=0)
ax2.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax2.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
ax.spines[spine].set_visible(False)
plt.show();
def deriv(y, t, N, beta, gamma, delta):
S, E, I, R = y
dSdt = -beta * S * I / N
dEdt = beta * S * I / N - delta * E
dIdt = delta * E - gamma * I
dRdt = gamma * I
return dSdt, dEdt, dIdt, dRdt
def deriv1(y, t, N, R0):
S, E, I, R = y
dSdt = -R0* 0.1 * S * I / N
dEdt = R0*0.1 * S * I / N - (1/5.5) * E
dIdt = (1/5.5) * E - 0.1 * I
dRdt = 0.1 * I
return dSdt, dEdt, dIdt, dRdt
def deriv2(y, t, N, R0, gamma, delta):
S, E, I, R = y
dSdt = -R0*gamma * S * I / N
dEdt = R0*gamma * S * I / N - delta * E
dIdt = delta * E - gamma * I
dRdt = gamma * I
return dSdt, dEdt, dIdt, dRdt
N = 474
D = 10.0 # infections lasts ten days
gamma = 1.0 / D # fix it 0.1
delta = 1.0 / 5.5 # incubation period of five.five days
R_0 = 3.28
beta = R_0 * gamma # R_0 = beta / gamma, so beta = R_0 * gamma
S0, E0, I0, R0 = N-1, 0, 1, 0 # initial conditions: one exposed
t = np.linspace(0, 100, 100) # Grid of time points (in days)
y0 = S0, E0, I0, R0 # Initial conditions vector
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma, delta))
S, E, I, R = ret.T
# Integrate the SIR equations over the time grid, t.
ret1 = odeint(deriv1, y0, t, args = (N, R0))
S1, E1, I1, R1 = ret1.T
# Integrate the SIR equations over the time grid, t.
ret2 = odeint(deriv2, y0, t, args = (N, R0, gamma, delta))
S2, E2, I2, R2 = ret2.T
plotseird(t, S, E, I, R)
plotseird(t, S1, E1, I1, R1)
plotseird(t, S2, E2, I2, R2)#
# get the real data
# fit the curve
data = pd.read_csv('SEIR_Jaisalmer_new.csv')
print(data)
print(data.shape)
def fit_odeint(x, beta, gamma,delta):
return odeint(deriv, y0, x, args=(N, beta, gamma, delta))[:,2]
def fit_odeint1(x, R0):
return odeint(deriv2, y0, x, args=(N, R0))[:,2]
def optimize_sliced_data(x,y):
# Do the sampling 100 times to get the error estimates
popt, pcov = optimize.curve_fit(fit_odeint, x, y)
fitted = fit_odeint(x, *popt)
return popt, pcov, fitted
def show_parameters(popt, pcov):
print(f'beta: {popt[0]}')
print(f'gamma: {popt[1]}')
print(f'delta: {popt[2]}')
perr = np.sqrt(np.diag(pcov))
print('Covariance Matrix', pcov)
print('Error estimates', perr)
#print(type(popt))
R0 = popt[0]/popt[1]
print("R0 value : %2.4f" %(R0) )
return popt[0], popt[1],perr, pcov, R0
def optimize_sliced_data1(x,y):
# Do the sampling 100 times to get the error estimates
popt, pcov = optimize.curve_fit(fit_odeint1, x, y)
return popt, pcov
x = [i for i in range(0,data.shape[0])]
y = np.asarray(data['I'])
x1 = x[0:20]
y1 = y[0:20]
popt11, pcov11 = optimize_sliced_data1(x1,y1)
#beta1,gamma1, errors1,cov1, R01 = show_parameters(popt1, pcov1)
#print(popt1[0]+1.95596*np.sqrt(pcov1[0,0]))
print(popt11)
x1 = x[0:20]
y1 = y[0:20]
popt1, pcov1, fit = optimize_sliced_data(x1,y1)
beta1,gamma1, errors1,cov1, R01 = show_parameters(popt1, pcov1)
print(popt1[0]+1.95596*np.sqrt(pcov1[0,0]))
x2 = x[20:]
y2 = y[20:]
popt2, pcov2 = optimize_sliced_data(x2,y2)
beta2,gamma2, errors2,cov2, R02 = show_parameters(popt2, pcov2)
#print('95% Confidence interval' popt2[0]-1.95596*np.sqrt(pcov2[0,0]))
#print(fitted)
#plt.plot(x, y, '-')
#plt.plot(x, fitted, '-')
#plt.show()
# prepare confidence level curves
#nstd = 1.0 # to draw 5-sigma intervals
#popt_up = popt + nstd*perr
#popt_dw = popt – nstd*perr
#popt_dw= popt-nstd*perr
fit1 = fit_odeint(x1, *popt1)
fit2= fit_odeint(x2, *popt2)
#fit_up = fit_odeint(x, *popt_up)
#fit_dw = fit_odeint(x, *popt_dw)
#print(fit_dw)
#print(fit_up)
#plot
fig, (ax1,ax2) = plt.subplots(1,2, figsize = (12, 5))
ax1.set_xlabel('Days since evacuation', fontsize=18)
ax1.set_ylabel('No of Infected', fontsize=18)
#ax1.set_title('fit with only Infected' , fontsize=18)
ax1.plot(x1, fit1, 'r', lw=2, label='best fit curve')
ax1.plot(x1, y1, 'k-', lw=2, label='True curve')
#plt.plot(x,y,'k-', lw=2)
#ax.fill_between(x,fit_up, fit_dw, alpha = 0.25)
#ax.fill_between(x, fit_up, fit_dw, alpha=.25, label='5-sigma interval')
ax1.text(2, 15, 'R0 value : %2.2f' %(R01) , style='italic',fontsize=20)
#bbox={'facecolor': 'red', 'alpha': 0.5, 'pad': 10})
ax2.set_xlabel('Days since evacuation', fontsize=18)
ax2.set_ylabel('No of Infected', fontsize=18)
#ax2.set_title('fit with only Infected' , fontsize=18)
ax2.plot(x2, fit2, 'r', lw=2, label='best fit curve')
ax2.plot(x2, y2, 'k-', lw=2, label='True curve')
#plt.plot(x,y,'k-', lw=2)
#ax.fill_between(x,fit_up, fit_dw, alpha = 0.25)
#ax.fill_between(x, fit_up, fit_dw, alpha=.25, label='5-sigma interval')
ax2.legend(loc='upper right',fontsize=18)
ax2.text(29, 4.7, 'R0 value : %2.2f' %(R02) , style='italic',fontsize = 20)
#bbox={'facecolor': 'red', 'alpha': 0.5, 'pad': 10})
plt.tight_layout()
plt.savefig('Initial_and_final.png')
plt.show()
# take all the data
popt, pcov = optimize_sliced_data(x,y)
beta, gamma, errors,cov, R0 = show_parameters(popt, pcov)
fit = fit_odeint(x, *popt)
#fit_up = fit_odeint(x, *popt_up)
#fit_dw = fit_odeint(x, *popt_dw)
#print(fit_dw)
#print(fit_up)
#plot
fig, (ax1) = plt.subplots(1,1, figsize = (8, 5))
ax1.set_xlabel('Days since evacuation', fontsize=18)
ax1.set_ylabel('No of Infected', fontsize=18)
#ax1.set_title('fit with only Infected' , fontsize=18)
ax1.plot(x, fit, 'r', lw=2, label='best fit curve')
ax1.plot(x, y, 'k-', lw=2, label='True curve')
#plt.plot(x,y,'k-', lw=2)
#ax.fill_between(x,fit_up, fit_dw, alpha = 0.25)
#ax.fill_between(x, fit_up, fit_dw, alpha=.25, label='5-sigma interval')
ax1.text(25, 15, 'R0 value : %2.2f' %(R0) , style='italic', fontsize = 20)
#bbox={'facecolor': 'red', 'alpha': 0.5, 'pad': 10})
ax1.legend(loc='upper right',fontsize=18)
plt.tight_layout()
plt.savefig('Overall.png')
plt.show()