|
24 | 24 |
|
25 | 25 | # Construct the network G |
26 | 26 | # ~~~~~~~~~~~~~~~~~~~~~~~ |
27 | | -numNodes = 30000 |
28 | | -baseGraph = networkx.barabasi_albert_graph(n=numNodes, m=7) |
| 27 | +numNodes = 90000 |
| 28 | +baseGraph = networkx.barabasi_albert_graph(n=numNodes, m=3) |
29 | 29 | # Baseline normal interactions: |
30 | | -G_norm = models.custom_exponential_graph(baseGraph, scale=200) |
| 30 | +G_norm = models.custom_exponential_graph(baseGraph, scale=500) |
31 | 31 | models.plot_degree_distn(G_norm, max_degree=40) |
32 | 32 |
|
33 | 33 | # Construct the network G under social distancing |
34 | 34 | # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
35 | | -numNodes = 30000 |
36 | | -baseGraph = networkx.barabasi_albert_graph(n=numNodes, m=2) |
| 35 | +numNodes = 90000 |
| 36 | +baseGraph = networkx.barabasi_albert_graph(n=numNodes, m=1) |
37 | 37 | # Baseline normal interactions: |
38 | | -G_dist = models.custom_exponential_graph(baseGraph, scale=20000) |
| 38 | +G_dist = models.custom_exponential_graph(baseGraph, scale=200000) |
39 | 39 | models.plot_degree_distn(G_dist, max_degree=40) |
40 | 40 |
|
41 | 41 | # Define model |
42 | 42 | # ~~~~~~~~~~~~ |
43 | 43 | model = models.SEIRSNetworkModel( |
44 | 44 | # network connectivty |
45 | 45 | G = G_norm, |
46 | | - p = 0.6, |
| 46 | + p = 0.51, |
47 | 47 | # clinical parameters |
48 | | - beta = 0.05, |
49 | | - sigma = 5.2, |
| 48 | + beta = 0.20, |
| 49 | + sigma = 4.0, |
| 50 | + omega = 1.5, |
50 | 51 | zeta = 0, |
51 | 52 | a = 0.43, # probability of an asymptotic (supermild) infection |
52 | 53 | m = 1-0.43, # probability of a mild infection |
53 | 54 | h = 0.20, # probability of hospitalisation for a mild infection |
54 | 55 | c = 2/3, # probability of hospitalisation in cohort |
55 | 56 | mi = 1/6, # probability of hospitalisation in midcare |
56 | | - da = 14, # days of infection when asymptomatic (supermild) |
57 | | - dm = 14, # days of infection when mild |
| 57 | + da = 6.5, # days of infection when asymptomatic (supermild) |
| 58 | + dm = 6.5, # days of infection when mild |
58 | 59 | dc = 7, |
59 | 60 | dmi = 14, |
60 | 61 | dICU = 14, |
61 | 62 | dICUrec = 6, |
62 | 63 | dmirec = 6, |
63 | | - dhospital = 3, # days before reaching the hospital when heavy or critical |
| 64 | + dhospital = 5, # days before reaching the hospital when heavy or critical |
64 | 65 | m0 = 0.49, # mortality in ICU |
65 | 66 | maxICU = 2000, |
66 | 67 | # testing |
67 | 68 | theta_S = 0, |
68 | 69 | theta_E = 0, |
| 70 | + theta_I = 0, |
69 | 71 | theta_A = 0, |
70 | 72 | theta_M = 0, |
71 | 73 | theta_R = 0, |
|
75 | 77 | # back-tracking |
76 | 78 | phi_S = 0, |
77 | 79 | phi_E = 0, |
| 80 | + phi_I = 0, |
78 | 81 | phi_A = 0, |
79 | 82 | phi_R = 0, |
80 | 83 | # initial condition |
81 | 84 | initN = 11.43e6, #results are extrapolated to entire population |
82 | | - initE = 100, |
| 85 | + initE = 0, |
| 86 | + initI = 3, |
83 | 87 | initA = 0, |
84 | 88 | initM = 0, |
85 | 89 | initC = 0, |
|
89 | 93 | initD = 0, |
90 | 94 | initSQ = 0, |
91 | 95 | initEQ = 0, |
| 96 | + initIQ = 0, |
92 | 97 | initAQ = 0, |
93 | 98 | initMQ = 0, |
94 | 99 | initRQ = 0, |
95 | 100 | # monte-carlo sampling |
96 | | - monteCarlo = True, |
| 101 | + monteCarlo = False, |
97 | 102 | repeats = 1 |
98 | 103 | ) |
| 104 | + |
99 | 105 | # Load data |
100 | 106 | # ~~~~~~~~~~~~ |
101 | 107 | #[index,data] = model.obtainData() |
|
111 | 117 | # vector with dates |
112 | 118 | index=pd.date_range('2020-03-13', freq='D', periods=ICUvect.size) |
113 | 119 | # data series used to calibrate model must be given to function 'plotFit' as a list |
114 | | -idx = -26 |
| 120 | +idx = -23 |
115 | 121 | index = index[0:idx] |
116 | 122 | data=[np.transpose(ICUvect[:,0:idx]),np.transpose(hospital[:,0:idx])] |
117 | 123 | # set optimisation settings |
118 | | -parNames = ['beta','p'] # must be a list! |
| 124 | +parNames = ['beta'] # must be a list! |
119 | 125 | positions = [np.array([6]),np.array([4,5,6])] # must be a list! |
120 | | -bounds=[(10,100),(0.1,0.5),(0.4,0.8)] # must be a list! |
| 126 | +bounds=[(10,100),(0.25,0.60)] # must be a list! |
121 | 127 | weights = np.array([0,1]) |
122 | 128 | # run optimisation |
123 | | -theta = model.fit(data,parNames,positions,bounds,weights,setvar=True,maxiter=10,popsize=5) |
| 129 | +theta = model.fit(data,parNames,positions,bounds,weights,setvar=True,maxiter=15,popsize=multiprocessing.cpu_count()) |
124 | 130 | stop = timeit.default_timer() |
125 | 131 | print('Required time: ', stop - start,' seconds') |
126 | 132 |
|
127 | 133 | # Make a graphical representation of results |
128 | 134 | # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
129 | | -model.plotFit(index,data,positions,modelClr=['red','orange'],legendText=('ICU (model)','ICU (data)','Hospital (model)','Hospital (data)'),titleText='Belgium',filename-'calibration.svg') |
| 135 | +model.plotFit(index,data,positions,modelClr=['red','orange'],legendText=('ICU (model)','ICU (data)','Hospital (model)','Hospital (data)'),titleText='Belgium',filename-'calibration90K.svg') |
0 commit comments