-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathgenFigu08.py
150 lines (140 loc) · 8.01 KB
/
genFigu08.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
from importCOVID19 import COVID19
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 10
plt.rcParams['xtick.labelsize'] = 8
plt.rcParams['ytick.labelsize'] = 8
plt.rc('figure', figsize=(10, 5))
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
# For statistical plots
import seaborn as sns
import scipy.stats as st
from scipy.stats import norm
# Reading parameters fit up to 2020/05/15
fitpars = pd.read_csv('20200525fits.csv')
# Excluded states
toexclude = ['TOT','CMX','MEX','GR0','GR1','GR2','COA','BCS','DUR','ZAC','COL','ROO','AGU']
# Excluding further states for high R0 values
#toexclude.extend(['OAX','SON','GRO','HID','SLP','QUE','CAM']) # From data 15/05
toexclude.extend(['OAX','GRO','SLP','CHH','MOR','BCN','CHP']) # From data 25/05
mainpars = fitpars.loc[(fitpars['ID'] == 'TOT') | (fitpars['ID'] == 'CMX') | (fitpars['ID'] == 'MEX')]
print(mainpars)
fitpars = fitpars[~fitpars['ID'].isin(toexclude)]
# Sorting pars in terms of the maximum infected value
fitpars.sort_values(by=['maxI'], inplace=True,ascending=False)
print('Number of states to plot:', len(fitpars.index))
# Get the parVals for a given fitpars row
def getPars(row,stat=False):
if stat:
return (row[1],row[7]*row[5],row[7],row[6],row[2],row[8])
# row = [ID, beta, eta]
# beta, tau, q, delta, eta, epsilon
q = 1; epsilon = 0.2; delta = 1/10
return (row[1],q*2/3,q,delta,row[2],epsilon)
# General definitions
maxX = 600; startV = 1
xtime = np.linspace(startV, maxX + startV -1,maxX)
xdate = pd.to_datetime(pd.Series(pd.date_range('20200315', periods=maxX)))
rawData = '20200525_Start1503.csv'
randData = '20200526_LAST/'
colors = sns.cubehelix_palette(10) #02
poscol = 6
mxCOVID19 = COVID19('SEIRfull', rawData)
mxCOVID19.readStat(randData)
plt.rc('figure', figsize=(12, 12))
fig, axs = plt.subplots(4, 4, gridspec_kw={'wspace': 0.2, 'hspace' : 0.08}, sharex=True, sharey=True)
useconf = True; conf = 0.95
nstate = 0; axins = []; source = mainpars; maxstates = len(source.index)
outres = []; nscheme = 0;
schemes = [6 for i in range(len(fitpars.index)+3)]
startMonthsName = ['Nov','Sep','Oct','Nov','Jan','Nov','Nov','Nov','Jan','Nov','Jan','Nov','Dic','Dic','Oct','Nov','Oct','Nov','Dic']
monthsDict = {'Sep' : 171, 'Oct' : 201, 'Nov' : 232, 'Dic' : 262, 'Jan' : 293 }
startMonths = [ monthsDict[mnth] for mnth in startMonthsName ]
startMonths = [225,181,215,228,292,234,229,239,307,230,297,226,264,267,208,236,215,241,260] # From data 15/05
startMonths = [240,188,233,294,255,250,261,276,222,249,243,340,297,288,244,240,287,261,240] # From data 25/05
for irow in range(4):
for icol in range(4):
print(irow,icol,nstate,maxstates)
for ps in range(maxstates+1):
if (nstate<(len(source.index)+maxstates)):
if ps==maxstates:
source = fitpars
maxstates = 0
print(ps,maxstates)
row = source.iloc[nstate if maxstates==0 else ps].values
state = row[0]
mxCOVID19 = COVID19('SEIRfullV01', rawData,[schemes[nscheme],20,startMonths[nscheme],-1,1])
mxCOVID19.readStat(randData)
nscheme += 1
print('Computing ',state,'(',irow,',',icol,')')
# Getting the scatter data
data = mxCOVID19.getFitData(state)
# Getting the best fitted curve under the given hypothesis
# Setting pars
mxCOVID19.parVals = getPars(row)
yValsModel = mxCOVID19.getModel(xtime)[:,1]
# Analysis of the stat data
mxCOVID19.getStatData(state)
allPars = mxCOVID19.actStatData.copy()
statData = allPars['eta'].values
limitsETA = st.t.interval(conf, len(statData)-1, loc=np.mean(statData), scale=st.sem(statData))
# Initializing the min and max y limits for the time range
ymin = np.zeros(maxX)+1e8; ymax = np.zeros(maxX)
minR0 = 1e5; maxR0 = 0; cnt = 0
# Cycle over the statpars
for i in range(len(allPars.index)):
statpars = allPars.iloc[i].values
if ((not useconf) or (statpars[2]>=limitsETA[0]) and (statpars[2]<=limitsETA[1])):
minR0 = minR0 if (minR0<statpars[3]) else statpars[3]
maxR0 = maxR0 if (maxR0>statpars[3]) else statpars[3]
mxCOVID19.parVals = getPars(statpars,True)
mxCOVID19.fE0 = statpars[9]
yVals = mxCOVID19.getModel(xtime)[:,1]
# Running over x to update max and min values
cnt+=1
for xindx in range(maxX):
ymin[xindx] = ymin[xindx] if (ymin[xindx]<yVals[xindx]) else yVals[xindx]
ymax[xindx] = ymax[xindx] if (ymax[xindx]>yVals[xindx]) else yVals[xindx]
print(conf*100,'% on eta of '+state+':',limitsETA,'. Min R0=',minR0,', Max R0=',maxR0,' useconf=', useconf,' pars used: ',cnt)
outres.append([state,np.max(ymin),xdate.iloc[np.argmax(ymin, axis=0)],minR0,maxR0])
if (nstate<(len(source.index)+maxstates)):
# The main axis
# Plotting the best fit
axs[irow,icol].plot(xdate,yValsModel,linewidth=0.8,color=colors[poscol])
# Plotting the shadow probabilities
axs[irow,icol].fill_between(xdate, ymin, ymax, alpha=0.3,color=colors[poscol])
axs[irow,icol].tick_params(labelrotation=22)
axs[irow,icol].ticklabel_format(axis="y", style="sci", scilimits=(0,0))
upv = 0.05; starttext = 0.10; sz = 6; xpos = 0.67
axs[irow,icol].text(xpos,starttext+upv*2, r'$\beta$='+'{:.2f}'.format(row[1]), transform=axs[irow,icol].transAxes,size = sz)
axs[irow,icol].text(xpos,starttext+upv*1, r'$\eta$='+'{:.2f}'.format(row[2]), transform=axs[irow,icol].transAxes,size = sz)
axs[irow,icol].text(xpos,starttext+upv*0, r'$R_0$='+'{:.1f}'.format(row[3]), transform=axs[irow,icol].transAxes,size = sz)
axs[irow,icol].text(xpos,starttext+upv*4, r'$I_{m}$='+format(int(row[5]), ',d'), transform=axs[irow,icol].transAxes,size = sz)
axs[irow,icol].text(xpos,starttext+upv*3, r'$D_{m}$='+row[6], transform=axs[irow,icol].transAxes,size = sz)
axs[irow,icol].text(0.45,0.95, state, transform=axs[irow,icol].transAxes)#,color='blue', transform=ax.transAxes)
# The inset axis
axins.append(inset_axes(axs[irow,icol], width="30%", height="40%" ,borderpad=1.1))#, loc='lower left', bbox_to_anchor=(0.8, 0.8, 1, 1)))
# Plotting the data
axins[-1].scatter(mxCOVID19.realDays.values,data['Infected'],color='gray',s=2)
# Plotting the shadow probabilities
axins[-1].fill_between(xtime, ymin, ymax, alpha=0.3,color=colors[poscol])
# Plotting the best fit
axins[-1].plot(xtime,yValsModel,linewidth=0.8,color=colors[poscol])
axins[-1].set_xlim(mxCOVID19.realDays.values[0],mxCOVID19.realDays.values[-1]) #self.fitData.iloc[0].values
axins[-1].set_ylim(0,np.max(data['Infected'].values)) #self.fitData.iloc[0].values
axins[-1].tick_params(direction='out', labelsize = 6)
axins[-1].ticklabel_format(axis="y", style="sci", scilimits=(0,0))
print('Maximum for the most optimistic scenario: ', np.max(ymin), " achieved on ", xdate.iloc[np.argmax(ymin, axis=0)])
# Moving to the next state
nstate += 1
outres = pd.DataFrame(outres,columns = ['ID','MaxInf','Date','R0min','R0max'])
outres.to_csv (rawData[:8]+'_bestScenarioStrategy.csv', index = False, header=True)
# add a big axis, hide frame
fig.add_subplot(111, frameon=False)
# hide tick and tick label of the big axis
plt.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
plt.ylabel('Infected')
plt.tight_layout()
plt.savefig('FigureFITstatesSchemes.png',dpi=600)
plt.show()