Merge pull request #37 from twallema/master

twallema · web-flow · commit b902ed8eb9f7 · 2020-05-12T10:53:34.000+02:00
Modification to ideal network parameters in hpc scripts and demo
diff --git a/hpc/calibrate_stochastic.py b/hpc/calibrate_stochastic.py
@@ -24,48 +24,50 @@
 
 # Construct the network G
 # ~~~~~~~~~~~~~~~~~~~~~~~
-numNodes = 30000
-baseGraph    = networkx.barabasi_albert_graph(n=numNodes, m=7)
+numNodes = 90000
+baseGraph    = networkx.barabasi_albert_graph(n=numNodes, m=3)
 # Baseline normal interactions:
-G_norm     = models.custom_exponential_graph(baseGraph, scale=200)
+G_norm     = models.custom_exponential_graph(baseGraph, scale=500)
 models.plot_degree_distn(G_norm, max_degree=40)
 
 # Construct the network G under social distancing
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-numNodes = 30000
-baseGraph    = networkx.barabasi_albert_graph(n=numNodes, m=2)
+numNodes = 90000
+baseGraph    = networkx.barabasi_albert_graph(n=numNodes, m=1)
 # Baseline normal interactions:
-G_dist     = models.custom_exponential_graph(baseGraph, scale=20000)
+G_dist     = models.custom_exponential_graph(baseGraph, scale=200000)
 models.plot_degree_distn(G_dist, max_degree=40)
 
 # Define model
 # ~~~~~~~~~~~~
 model = models.SEIRSNetworkModel(
                                  # network connectivty
                                  G = G_norm,
-                                 p       = 0.6,
+                                 p       = 0.51,
                                  # clinical parameters
-                                 beta    = 0.05, 
-                                 sigma   = 5.2,
+                                 beta    = 0.20, 
+                                 sigma   = 4.0,
+                                 omega   = 1.5,
                                  zeta    = 0,
                                  a = 0.43, # probability of an asymptotic (supermild) infection
                                  m = 1-0.43, # probability of a mild infection
                                  h = 0.20, # probability of hospitalisation for a mild infection
                                  c = 2/3, # probability of hospitalisation in cohort
                                  mi = 1/6, # probability of hospitalisation in midcare
-                                 da = 14, # days of infection when asymptomatic (supermild)
-                                 dm = 14, # days of infection when mild
+                                 da = 6.5, # days of infection when asymptomatic (supermild)
+                                 dm = 6.5, # days of infection when mild
                                  dc = 7,
                                  dmi = 14,
                                  dICU = 14,
                                  dICUrec = 6,
                                  dmirec = 6,
-                                 dhospital = 3, # days before reaching the hospital when heavy or critical
+                                 dhospital = 5, # days before reaching the hospital when heavy or critical
                                  m0 = 0.49, # mortality in ICU
                                  maxICU = 2000,
                                  # testing
                                  theta_S = 0,
                                  theta_E = 0,
+                                 theta_I = 0,
                                  theta_A = 0,
                                  theta_M = 0,
                                  theta_R = 0,
@@ -75,11 +77,13 @@
                                  # back-tracking
                                  phi_S   = 0,
                                  phi_E   = 0,
+                                 phi_I   = 0,
                                  phi_A   = 0,
                                  phi_R   = 0,
                                  # initial condition
                                  initN = 11.43e6, #results are extrapolated to entire population
-                                 initE = 100,
+                                 initE = 0,
+                                 initI = 3,
                                  initA = 0, 
                                  initM = 0,
                                  initC = 0,
@@ -89,13 +93,15 @@
                                  initD = 0,
                                  initSQ = 0,
                                  initEQ = 0,
+                                 initIQ = 0,
                                  initAQ = 0,
                                  initMQ = 0,
                                  initRQ = 0,
                                  # monte-carlo sampling
-                                 monteCarlo = True,
+                                 monteCarlo = False,
                                  repeats = 1
                             )
+
 # Load data
 # ~~~~~~~~~~~~
 #[index,data] = model.obtainData()
@@ -111,19 +117,19 @@
 # vector with dates
 index=pd.date_range('2020-03-13', freq='D', periods=ICUvect.size)
 # data series used to calibrate model must be given to function 'plotFit' as a list
-idx = -26
+idx = -23
 index = index[0:idx]
 data=[np.transpose(ICUvect[:,0:idx]),np.transpose(hospital[:,0:idx])]
 # set optimisation settings
-parNames = ['beta','p'] # must be a list!
+parNames = ['beta'] # must be a list!
 positions = [np.array([6]),np.array([4,5,6])] # must be a list!
-bounds=[(10,100),(0.1,0.5),(0.4,0.8)] # must be a list!
+bounds=[(10,100),(0.25,0.60)] # must be a list!
 weights = np.array([0,1])
 # run optimisation
-theta = model.fit(data,parNames,positions,bounds,weights,setvar=True,maxiter=10,popsize=5)
+theta = model.fit(data,parNames,positions,bounds,weights,setvar=True,maxiter=15,popsize=multiprocessing.cpu_count())
 stop = timeit.default_timer()
 print('Required time: ', stop - start,' seconds')
 
 # Make a graphical representation of results
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-model.plotFit(index,data,positions,modelClr=['red','orange'],legendText=('ICU (model)','ICU (data)','Hospital (model)','Hospital (data)'),titleText='Belgium',filename-'calibration.svg')
+model.plotFit(index,data,positions,modelClr=['red','orange'],legendText=('ICU (model)','ICU (data)','Hospital (model)','Hospital (data)'),titleText='Belgium',filename-'calibration90K.svg')
diff --git a/hpc/sim_stochastic.py b/hpc/sim_stochastic.py
@@ -24,43 +24,44 @@
 
 # Construct the network G
 # ~~~~~~~~~~~~~~~~~~~~~~~
-numNodes = 60000
-baseGraph    = networkx.barabasi_albert_graph(n=numNodes, m=7)
+numNodes = 90000
+baseGraph    = networkx.barabasi_albert_graph(n=numNodes, m=3)
 # Baseline normal interactions:
-G_norm     = models.custom_exponential_graph(baseGraph, scale=200)
+G_norm     = models.custom_exponential_graph(baseGraph, scale=500)
 models.plot_degree_distn(G_norm, max_degree=40)
 
 # Construct the network G under social distancing
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-numNodes = 60000
-baseGraph    = networkx.barabasi_albert_graph(n=numNodes, m=2)
+numNodes = 90000
+baseGraph    = networkx.barabasi_albert_graph(n=numNodes, m=1)
 # Baseline normal interactions:
-G_dist     = models.custom_exponential_graph(baseGraph, scale=20000)
+G_dist     = models.custom_exponential_graph(baseGraph, scale=200000)
 models.plot_degree_distn(G_dist, max_degree=40)
 
 # Define model
 # ~~~~~~~~~~~~
 model = models.SEIRSNetworkModel(
                                  # network connectivty
                                  G = G_norm,
-                                 p       = 0.6,
+                                 p       = 0.51,
                                  # clinical parameters
-                                 beta    = 0.05, 
-                                 sigma   = 5.2,
+                                 beta    = 0.03, 
+                                 sigma   = 4.0,
+                                 omega   = 1.5,
                                  zeta    = 0,
                                  a = 0.43, # probability of an asymptotic (supermild) infection
                                  m = 1-0.43, # probability of a mild infection
                                  h = 0.20, # probability of hospitalisation for a mild infection
                                  c = 2/3, # probability of hospitalisation in cohort
                                  mi = 1/6, # probability of hospitalisation in midcare
-                                 da = 14, # days of infection when asymptomatic (supermild)
-                                 dm = 14, # days of infection when mild
+                                 da = 6.5, # days of infection when asymptomatic (supermild)
+                                 dm = 6.5, # days of infection when mild
                                  dc = 7,
                                  dmi = 14,
                                  dICU = 14,
                                  dICUrec = 6,
                                  dmirec = 6,
-                                 dhospital = 3, # days before reaching the hospital when heavy or critical
+                                 dhospital = 5, # days before reaching the hospital when heavy or critical
                                  m0 = 0.49, # mortality in ICU
                                  maxICU = 2000,
                                  # testing
@@ -75,11 +76,13 @@
                                  # back-tracking
                                  phi_S   = 0,
                                  phi_E   = 0,
+                                 phi_I   = 0,
                                  phi_A   = 0,
                                  phi_R   = 0,
                                  # initial condition
                                  initN = 11.43e6, #results are extrapolated to entire population
-                                 initE = 100,
+                                 initE = 0,
+                                 initI = 1,
                                  initA = 0, 
                                  initM = 0,
                                  initC = 0,
@@ -89,12 +92,13 @@
                                  initD = 0,
                                  initSQ = 0,
                                  initEQ = 0,
+                                 initIQ = 0,
                                  initAQ = 0,
                                  initMQ = 0,
                                  initRQ = 0,
                                  # monte-carlo sampling
-                                 monteCarlo = True,
-                                 repeats = 10
+                                 monteCarlo = False,
+                                 repeats = 3
                             )
 
 # Create checkpoints dictionary
diff --git a/src/SEIRSNetworkModel_demo.ipynb b/src/SEIRSNetworkModel_demo.ipynb
