running python through R with reticulate

Analysis/ode.R

@@ -1,11 +1,12 @@
+## Set Working Directory to file directory - RStudio approach
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+
 seirds.f <- function(t, y, k) {
   S <- y[1]
   E <- y[2]
-  I_R <- y[3]
+  I <- y[3]
-  I_D <- y[4]
+  R <- y[4]
-  R <- y[5]
+  D <- y[5]
-  D <- y[6]
-  N <- y[7]
 
   beta = k[1]
   delta = k[2]
@@ -15,64 +16,64 @@ seirds.f <- function(t, y, k) {
   epsilon = k[6]
   waning = k[7]
 
-  dS <- S - beta*(I_R+I_D)/N + waning*R
+  # Saying infec rate of S in contact with E same as contact with I
-  dE <- E + beta*S*(I_R+I_D)/N - delta*E + epsilon
+  dS <- epsilon - (beta*E + beta*I)*S + waning*R - mu*S
-  dI_R <- I_R + delta*(1-mu)*E - gamma_r*I_R + epsilon
+  dE <- beta*S*E + beta*S*I - (delta+mu)*E
-  dI_D <- I_D + delta*mu*E - gamma_d*I_D + epsilon
+  dI <- delta*E - mu*I - gamma_r*I
-  dR <- R + gamma_r*I_R - waning*R
+  dR <- gamma_r*I - (waning+mu)*R
-  dD <- D + gamma_d*I_D
+  dD <- gamma_d*I
-  return(as.matrix(c(dS, dE, dI_R, dI_D, dR, dD, N)))
+  
+  return(as.matrix(c(dS, dE, dI, dR, dD)))
 }
 
+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')
+
 # a) a=0.5, b=1, S(0)=0.9, I(0)=0.1, R(0)=0
 seirds.params <- c(0.7316455696202532,    # beta
                    0.020253164556962026,  # delta
                    0.09113924050632911,   # gamma_r
                    0.002531645569620253,  # gamma_d
                    0.0034602076124567475, # mu
-                   0.0)                   #epsilon
+                   0.0,                   # epsilon
-
+                   0.0)                   # waning
-seirds.params <- c(1,    # beta
-                   1,  # delta
-                   1,   # gamma_r
-                   1,  # gamma_d
-                   1, # mu
-                   1)                   #epsilon
-
 
 S <- 62
 E <- 8
-I_R <- 288
+I <- 288
-I_D <- 1
 R <- 36
 D <- 1
 
-N <- S+E+I_R+I_D+R+D
+N <- S+E+I+R+D
 S <- S/N
 E <- E/N
-I_R <- I_R/N
+I <- I/N
-I_D <- I_D/N
 R <- R/N
 D <- D/N
 
 tmin <- 0
-tmax <- 20
+tmax <- 7
 
-install.packages("pracma")
 library(pracma)
 seirds.ode.sol <- ode45(f = function(t,y){seirds.f(t,y,k=seirds.params)},
-                     y = c(S, E, I_R, I_D, R, D, N),
+                     y = c(S, E, I, R, D),
                      t0 = tmin, tfinal = tmax)
+
 # plot
+library(reshape2)
+library(ggplot2)
 plot.seirds <- function(sol, method){
   # - pred/prey columns melted to 1 column called "value"
   sol.melted <- melt(data.frame(time=sol$t,
                                 S=sol$y[,1],
                                 E=sol$y[,2],
-                                I_R=sol$y[,3],
+                                I=sol$y[,3],
-                                I_D=sol$y[,4],
+                                R=sol$y[,4],
-                                R=sol$y[,5],
+                                D=sol$y[,5]),
-                                D=sol$y[,6],),
                      id = "time") # melt based on time
   # New column created with variable names, called “variable”
   colnames(sol.melted)[2] <- "Group" # used in legend
@@ -84,4 +85,6 @@ plot.seirds <- function(sol, method){
   show(g)
 }
 
-plot.seirds(seirds.ode.sol, "ODE45")
+plot.seirds(seirds.ode.sol, "ODE45")
+
+print(seirds.ode.sol$y[tmax,1]+seirds.ode.sol$y[tmax,2]+seirds.ode.sol$y[tmax,3]+seirds.ode.sol$y[tmax,4]+seirds.ode.sol$y[tmax,5])

Analysis/prep_model.py

@@ -7,7 +7,9 @@ from operator import getitem
 import itertools, os, sys
 
 # Change dir to location of this python file
-os.chdir(os.path.dirname(sys.argv[0]))
+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.
@@ -106,4 +108,4 @@ 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("epsilon:", str(epsilon))