Infection count edge weighting

Analysis/ode.R

@@ -1,6 +1,12 @@
 ## Set Working Directory to file directory - RStudio approach
 setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
 
+weighted <- 'False'
+#conda_install(envname = "r-reticulate", packages="networkx")
+#conda_install(envname = "r-reticulate", packages="matplotlib")
+#conda_install(envname = "r-reticulate", packages="pydot")
+#conda_install(envname = "r-reticulate", packages="pygraphviz")
+
 seirds.f <- function(t, y, k) {
   S <- y[1]
   E <- y[2]
@@ -18,7 +24,7 @@ seirds.f <- function(t, y, k) {
 
   # Saying infec rate of S in contact with E same as contact with I
   dS <- epsilon - (beta*E + beta*I)*S + waning*R - gamma_d*S
-  dE <- beta*S*E + beta*S*I - (delta+gamma_d)*E
+  dE <- (beta*E + beta*I)*S - (delta+gamma_d)*E
   dI <- delta*E - (1+gamma_d)*I
   dR <- (1-mu)*I - (waning+gamma_d)*R
   dD <- mu*I
@@ -27,11 +33,22 @@ seirds.f <- function(t, y, k) {
 }
 
 library(reticulate)
-#conda_install(envname = "r-reticulate", packages="networkx")
-#conda_install(envname = "r-reticulate", packages="matplotlib")
-#conda_install(envname = "r-reticulate", packages="pydot")
-#conda_install(envname = "r-reticulate", packages="pygraphviz")
 source_python('prep_model.py')
+model_data <- prep_seirds(weighted)
+
+S <- unlist(model_data)[1]
+E <- unlist(model_data)[2]
+I <- unlist(model_data)[3]
+R <- unlist(model_data)[4]
+D <- unlist(model_data)[5]
+
+beta <- unlist(model_data)[6]
+delta <- unlist(model_data)[7]
+gamma_r <- unlist(model_data)[8]
+gamma_d <- unlist(model_data)[9]
+mu <- unlist(model_data)[10]
+epsilon <- unlist(model_data)[11]
+omega <- unlist(model_data)[12]
 
 # Obtained from prep_model.py
 seirds.params <- c(beta,     # beta
@@ -79,6 +96,6 @@ plot.seirds <- function(sol, method){
 }
 
 plot.seirds(seirds.ode.sol, "ODE45")
-ggsave("SERIDS.pdf")
+ggsave("SEIRDS.pdf")
 # Sanity check: Make sure sums to ~1.0
 sum(tail(seirds.ode.sol$y,1))

Analysis/prep_model.py

@@ -11,30 +11,33 @@ print(os.getcwd())
 #os.chdir(os.path.dirname(sys.argv[0]))
 #print(os.getcwd())
 
-# AGraph preserves attributes, networkx Graph does not.
+def prep_seirds(weighted):
-# Many of the desired functions are in networkx.
+  print("Prepping the SEIRDS model, using trivial weighting=", weighted)
-# So import AGraph to keep attributes, then convert to Networkx.
+  # AGraph preserves attributes, networkx Graph does not.
-A = nx.drawing.nx_agraph.to_agraph(nx.drawing.nx_pydot.read_dot("./1_mo_color_DOTFILE.dot"))
+  # Many of the desired functions are in networkx.
-A.layout('dot')
+  # So import AGraph to keep attributes, then convert to Networkx.
-#A.draw('tree.png')
+  A = nx.drawing.nx_agraph.to_agraph(nx.drawing.nx_pydot.read_dot("./1_mo_color_DOTFILE.dot"))
-A.remove_node('\\n') # Remove "newline" node from newline end of dot file
+  A.layout('dot')
-G=nx.DiGraph(A)
+  #A.draw('tree.png')
+  A.remove_node('\\n') # Remove "newline" node from newline end of dot file
+  G=nx.DiGraph(A)
   
-color_map = []
+  color_map = []
-color_d = {}
+  color_d = {}
-node_pos = {} # used for drawing/graphing
+  node_pos = {} # used for drawing/graphing
   
-# Compartments
+  # Compartments
-S = 0
+  S = 0
-I_R = 0
+  I_R = 0
-I_D = 0
+  I_D = 0
-E = 0
+  E = 0
-R = 0
+  R = 0
-D = 0
+  D = 0
   
-ep_tmp = 0 # counter for epsilon
+  ep_tmp = 0 # counter for epsilon
+  inf_ct = 0 # infection rate counter
   
-for node in A:
+  for node in A:
       color = A.get_node(node).attr.to_dict()['fillcolor']
       str_pos = A.get_node(node).attr.to_dict()['pos']
       coords = str_pos.split(',')
@@ -45,7 +48,11 @@ for node in A:
           color_map.append("white")
           color_d[node] = color
           in_edges = list(G.in_edges(node))
+          out_edges = list(G.out_edges(node))
+          
+          
           tmp_S = 1
+          tmp_inf = 0
           for source in in_edges:
               tmp_S = 1
               # If previous node was infected, then we are recovered
@@ -53,6 +60,13 @@ for node in A:
                   R = R + 1
                   tmp_S = 0
                   break # No need to check the other nodes
+          for source in out_edges:
+            if (color_d[source[0]] == 'red'):
+             if(weighted == 'False' or not out_edges):
+                tmp_inf = 1
+            else:
+                inf_ct = inf_ct + 1/len(out_edges) # trivial weighting
+          inf_ct = inf_ct + tmp_inf
           S = S + tmp_S
               #G[source[0]][node]['weight'] = 3 
       elif color == 'yellow':
@@ -60,14 +74,21 @@ for node in A:
           color_d[node] = color
           in_edges = list(G.in_edges(node))
           tmp_E = 1
+          tmp_R = 0
+          tmp_inf = 0
           for source in in_edges:
-            tmp_E = 1
               # If previous node was infected, then we are recovered
               if (color_d[source[0]] == 'red'):
-                R = R + 1
+                  tmp_R = 1
                   tmp_E = 0
-                break # No need to check the other nodes
+              if (color_d[source[0]] == '' or color_d[source[0]] == 'white'):
+                  if(weighted == 'False' or not in_edges):
+                    tmp_inf = 1 # add 1 for the inf counter
+                  else:
+                    inf_ct = inf_ct + 1/len(in_edges) # trivial weighting
           E = E + tmp_E
+          R = R + tmp_R
+          inf_ct = inf_ct + tmp_inf
       else:
           color_map.append(color)
           color_d[node] = color
@@ -82,33 +103,44 @@ for node in A:
           in_edges = list(G.in_edges(node))
           if not in_edges:
               ep_tmp = ep_tmp + 1
+          else:
+            tmp_inf = 0
+            for source in in_edges:
+              if (color_d[source[0]] == '' or color_d[source[0]] == 'white'):
+                if(weighted == 'False' or not in_edges):
+                  tmp_inf = 1
+                else:
+                  inf_ct = inf_ct + 1/len(in_edges) # trivial weightin
+            inf_ct = inf_ct + tmp_inf
   
-# Params
+  # Params
-beta = (I_R+I_D)/len(A) # rate of infec (I/total?)
+  beta = (inf_ct)/len(A) # rate of infec
-#delta = E/len(A) # symptom appearance rate (E/total?)
+  delta = 1 # incubation period
-delta = 1 # incubation period
+  gamma_r = R/len(A) # recov rate
-gamma_r = R/len(A) # recov rate (R/total?)
+  gamma_d = D/len(A) # death rate
-gamma_d = D/len(A) # death rate (D/total?)
+  mu = D/(I_R+I_D) # fatality ratio
-mu = D/(I_R+I_D) # fatality ratio (D/I)
+  epsilon = ep_tmp/len(A) # infected import rate
-epsilon = ep_tmp/len(A) # infected import rate
+  omega = 1 # waning immunity rate
-omega = 1 # waning immunity rate
   
   
-print("Model Compartments:")
+  print("Model Compartments:")
-print("S:", str(S))
+  print("S:", str(S))
-print("E:", str(E))
+  print("E:", str(E))
-print("I_R:", str(I_R))
+  print("I_R:", str(I_R))
-print("I_D:", str(I_D))
+  print("I_D:", str(I_D))
-print("R:", str(R))
+  print("R:", str(R))
-print("D:", str(D))
+  print("D:", str(D))
-print("\n")
+  print("\n")
   
-print("Model Parameters:")
+  print("Model Parameters:")
-print("beta:", str(beta))
+  print("infect counter:", str(inf_ct))
-print("delta:", str(delta))
+  print("beta:", str(beta))
-print("gamma_r:", str(gamma_r))
+  print("delta:", str(delta))
-print("gamma_d:", str(gamma_d))
+  print("gamma_r:", str(gamma_r))
-print("mu:", str(mu))
+  print("gamma_d:", str(gamma_d))
-print("epsilon:", str(epsilon))
+  print("mu:", str(mu))
-print("omega:", str(omega))
+  print("epsilon:", str(epsilon))
+  print("omega:", str(omega))
+  
+  return (S, E, I_R+I_D, R, D, beta, delta, gamma_r, gamma_d, mu, epsilon, omega)