diff --git a/Analysis/SEIRDS.pdf b/Analysis/SEIRDS.pdf new file mode 100644 index 0000000..07d8447 Binary files /dev/null and b/Analysis/SEIRDS.pdf differ diff --git a/Analysis/SERIDS.pdf b/Analysis/SERIDS.pdf deleted file mode 100644 index d8cd1e2..0000000 Binary files a/Analysis/SERIDS.pdf and /dev/null differ diff --git a/Analysis/ode.R b/Analysis/ode.R index 1092930..75ac6b1 100644 --- a/Analysis/ode.R +++ b/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)) \ No newline at end of file diff --git a/Analysis/prep_model.py b/Analysis/prep_model.py index 5f30a7c..05194f9 100644 --- a/Analysis/prep_model.py +++ b/Analysis/prep_model.py @@ -11,104 +11,136 @@ print(os.getcwd()) #os.chdir(os.path.dirname(sys.argv[0])) #print(os.getcwd()) -# AGraph preserves attributes, networkx Graph does not. -# Many of the desired functions are in networkx. -# So import AGraph to keep attributes, then convert to Networkx. -A = nx.drawing.nx_agraph.to_agraph(nx.drawing.nx_pydot.read_dot("./1_mo_color_DOTFILE.dot")) -A.layout('dot') -#A.draw('tree.png') -A.remove_node('\\n') # Remove "newline" node from newline end of dot file -G=nx.DiGraph(A) - -color_map = [] -color_d = {} -node_pos = {} # used for drawing/graphing - -# Compartments -S = 0 -I_R = 0 -I_D = 0 -E = 0 -R = 0 -D = 0 - -ep_tmp = 0 # counter for epsilon - -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(',') - x = coords[0] # layout for draw function - y = coords[1] - node_pos[node] = float(x), float(y) - if color is None or color == '': - color_map.append("white") - color_d[node] = color - in_edges = list(G.in_edges(node)) - tmp_S = 1 - for source in in_edges: - tmp_S = 1 - # If previous node was infected, then we are recovered +def prep_seirds(weighted): + print("Prepping the SEIRDS model, using trivial weighting=", weighted) + # AGraph preserves attributes, networkx Graph does not. + # Many of the desired functions are in networkx. + # So import AGraph to keep attributes, then convert to Networkx. + A = nx.drawing.nx_agraph.to_agraph(nx.drawing.nx_pydot.read_dot("./1_mo_color_DOTFILE.dot")) + A.layout('dot') + #A.draw('tree.png') + A.remove_node('\\n') # Remove "newline" node from newline end of dot file + G=nx.DiGraph(A) + + color_map = [] + color_d = {} + node_pos = {} # used for drawing/graphing + + # Compartments + S = 0 + I_R = 0 + I_D = 0 + E = 0 + R = 0 + D = 0 + + ep_tmp = 0 # counter for epsilon + inf_ct = 0 # infection rate counter + + 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(',') + x = coords[0] # layout for draw function + y = coords[1] + node_pos[node] = float(x), float(y) + if color is None or color == '': + 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 + if (color_d[source[0]] == 'red'): + 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'): - R = R + 1 - tmp_S = 0 - break # No need to check the other nodes - S = S + tmp_S - #G[source[0]][node]['weight'] = 3 - elif color == 'yellow': - color_map.append(color) - color_d[node] = color - in_edges = list(G.in_edges(node)) - tmp_E = 1 - 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_E = 0 - break # No need to check the other nodes - E = E + tmp_E - else: - color_map.append(color) - color_d[node] = color - # Check if node dies - out_edges = list(G.out_edges(node)) - if not out_edges: - D = D + 1 - I_D = I_D + 1 - else: - I_R = I_R + 1 - # Check if imported - in_edges = list(G.in_edges(node)) - if not in_edges: - ep_tmp = ep_tmp + 1 - -# Params -beta = (I_R+I_D)/len(A) # rate of infec (I/total?) -#delta = E/len(A) # symptom appearance rate (E/total?) -delta = 1 # incubation period -gamma_r = R/len(A) # recov rate (R/total?) -gamma_d = D/len(A) # death rate (D/total?) -mu = D/(I_R+I_D) # fatality ratio (D/I) -epsilon = ep_tmp/len(A) # infected import rate -omega = 1 # waning immunity rate - - -print("Model Compartments:") -print("S:", str(S)) -print("E:", str(E)) -print("I_R:", str(I_R)) -print("I_D:", str(I_D)) -print("R:", str(R)) -print("D:", str(D)) -print("\n") - -print("Model Parameters:") -print("beta:", str(beta)) -print("delta:", str(delta)) -print("gamma_r:", str(gamma_r)) -print("gamma_d:", str(gamma_d)) -print("mu:", str(mu)) -print("epsilon:", str(epsilon)) -print("omega:", str(omega)) + 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': + color_map.append(color) + 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: + # If previous node was infected, then we are recovered + if (color_d[source[0]] == 'red'): + tmp_R = 1 + tmp_E = 0 + 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 + # Check if node dies + out_edges = list(G.out_edges(node)) + if not out_edges: + D = D + 1 + I_D = I_D + 1 + else: + I_R = I_R + 1 + # Check if imported + 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 + beta = (inf_ct)/len(A) # rate of infec + delta = 1 # incubation period + gamma_r = R/len(A) # recov rate + gamma_d = D/len(A) # death rate + mu = D/(I_R+I_D) # fatality ratio + epsilon = ep_tmp/len(A) # infected import rate + omega = 1 # waning immunity rate + + + print("Model Compartments:") + print("S:", str(S)) + print("E:", str(E)) + print("I_R:", str(I_R)) + print("I_D:", str(I_D)) + print("R:", str(R)) + print("D:", str(D)) + print("\n") + + print("Model Parameters:") + print("infect counter:", str(inf_ct)) + print("beta:", str(beta)) + print("delta:", str(delta)) + print("gamma_r:", str(gamma_r)) + print("gamma_d:", str(gamma_d)) + print("mu:", str(mu)) + 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)