benchmark <- function(dry.run=FALSE) { if (dry.run) cat("****************************************\n", "Dry run - no models will actually run...\n", "****************************************\n", sep="") source("dev.sourceAll.R") dev.sourceAll() ## # Define a series of models that will be used for the benchmarking (basically # the demo models but without any plotting and only minimum of concole output) ## Rprof("benchmark.Rprof") #=============================================================================== # Decaying-Dimerization Reaction Set (Gillespie, 2001) #=============================================================================== cat("Running 1/8...\n") flush.console() parms <- c(c1=1.0, c2=0.002, c3=0.5, c4=0.04) # Define parameters x0 <- c(s1=10000, s2=0, s3=0) # Initial state vector nu <- matrix(c(-1, -2, +2, 0, # State-change matrix 0, +1, -1, -1, 0, 0, 0, +1), nrow=3,byrow=TRUE) a <- c("c1*s1", "c2*s1*s1", "c3*s2", "c4*s2") # Propensity vector tf <- 10 # Final time simName <- "Decaying-Dimerizing Reaction Set" set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="D",simName,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="ETL",simName,tau=0.003,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="BTL",simName,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="OTL",simName,verbose=FALSE,consoleInterval=1) #=============================================================================== # SIRS metapopulation model (Pineda-Krch, 2007) #=============================================================================== cat("Running 2/8...\n") flush.console() patchPopSize <- 500 # Patch size U <- 20 # Number of patches # Define parameters parms <- c(beta = 0.001, # Transmission rate gamma = 0.1, # Recovery rate rho = 0.005, # Loss of immunity rate epsilon = 0.01, # Proportion inter-patch transmissions N = patchPopSize) # Patch population size (constant) # Create the named initial state vector for the U-patch system. The structure of # x0 is as follows (assuming a patchsize of 500 individuals), # x0 <- c(499,1,500,0,500,0 ... ) x0 <- c(S1=(patchPopSize-1), I1=1) if (U>1) { for (i in (seq(2,U))) { patch <- c(patchPopSize, 0) names(patch) <- c(paste("S",i,sep=""), paste("I",i,sep="")) x0 <- c(x0, patch) } } # State-change matrix # Define the state change matix for a single patch # Reaction 1 2 3 4 State nu <- matrix(c(-1, -1, 0, +1, # S +1, +1, -1, 0), # I nrow=2,byrow=TRUE) # Create the propensity vector a <- NULL for (patch in (seq(U))) { i <- patch # Intra-patch index if (patch==1) j <- U # Inter-patch index else j <- patch-1 # Construct the propensity functions for the current patch # Reaction: a_patch <- c(paste("(1-epsilon)*beta*S",i,"*I",i,"",sep=""), # 1. Intra-patch infection paste("epsilon*beta*S",i,"*I",j,"",sep=""), # 2. Inter-patch infection paste("gamma*I",i,"",sep=""), # 3. Recovery from infection paste("rho*(N-S",i,"-I",i,")",sep="")) # 4. Loss of immunity a <- c(a, a_patch) } # for() set.seed(1) if (!dry.run) out <- ssa(x0, a, nu, parms, tf=50, method = "D", simName="SIRS metapopulation model", verbose=FALSE, consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0, a, nu, parms, tf=50, method = "ETL", simName="SIRS metapopulation model", verbose=FALSE, consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0, a, nu, parms, tf=50, method = "BTL", simName="SIRS metapopulation model", verbose=FALSE, consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0, a, nu, parms, tf=50, method = "OTL", simName="SIRS metapopulation model", hor=rep(2,length(x0)), verbose=FALSE, consoleInterval=1) #=============================================================================== # Linear Chain System (Cao et al., 2004) #=============================================================================== cat("Running 3/8...\n") flush.console() # Rate parameter parms <- c(c=1) # Number of chain reactions M <- 50 # Initial state vector x0 <- c(1000, rep(0,M)) names(x0) <- paste("x",seq(M+1),sep="") # State-change matrix nu <- matrix(rep(0,(M*(M+1))),ncol=M) diag(nu) <- -1 diag(nu[2:M,]) <- +1 nu[M+1,M] <- +1 # Propensity vector a <- as.vector(paste("c*x",seq(M),"",sep="")) tf <- 5 # Final time simName <- "Linear Chain System" set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="D",simName,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="ETL",simName,tau=0.1,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="BTL",simName,f=50,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="OTL",simName,verbose=FALSE,consoleInterval=1) #=============================================================================== # Logistic growth (Pearl-Verhulst model) (Kot, 2001) #=============================================================================== cat("Running 4/8...\n") flush.console() parms <- c(b=2, d=1, K=1000) # Parameters x0 <- c(N=500) # Initial state vector nu <- matrix(c(+1,-1),ncol=2) # State-change matrix a <- c("b*N", "(d+(b-d)*N/K)*N") # Propensity vector tf <- 10 # Final time simName <- "Logistic growth" set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="D",simName,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="ETL",simName,tau=0.03,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="BTL",simName,f=5,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="OTL",simName,verbose=FALSE,consoleInterval=1) #=============================================================================== # Lotka predator-prey model (Gillespie, 1977; Kot, 2001) #=============================================================================== cat("Running 5/8...\n") flush.console() # Define parameters parms <- c(c1=10, c2=.01, c3=10) # Define system x0 <- c(Y1=1000, Y2=1000) # Initial state vector nu <- matrix(c(+1, -1, 0, 0, 1, -1),nrow=2,byrow=T) # State-change matrix a <- c("c1*Y1", "c2*Y1*Y2","c3*Y2") # Propensity vector tf <- 2 # Final time simName <- "Lotka predator-prey model" set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="D",simName="Lotka predator-prey model",verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="ETL",simName,tau=0.002,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="BTL",simName,f=100,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="OTL",simName,epsilon=0.1,verbose=FALSE,consoleInterval=1) #=============================================================================== # Radioactive decay model [aka Irreversible isomerization reaction set] # (Gillespie, 1977) #=============================================================================== cat("Running 6/8...\n") flush.console() parms <- c(k=0.5) # Define parameters x0 <- c(N=1000) # Initial state vector nu <- matrix(c(-1),nrow=1,byrow=TRUE) # State change matrix a <- c("k*N") # Propensity vector tf <- 20 # Final time simName <- "Radioactive decay model" if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="D",simName,verbose=FALSE) #=============================================================================== # Rosenzweig-MacArthur predator-prey model # (Pineda-Krch et al., 2007, Pineda-Krch, [submitted JSS manuscript]) #=============================================================================== cat("Running 7/8...\n") flush.console() # Define parameters # (B in Figure 1 in Pineda-Krch et al. 2007) parms <- c(b=2, d=1, K=1000, alpha=0.005, w=0.0025, c=2, g=2) # Parameters x0 <- c(N=500, P=500) # Initial state vector nu <- matrix(c(+1, -1, -1, 0, 0, # State-change matrix 0, 0, 0, +1, -1), nrow=2,byrow=TRUE) a <- c("b*N", # Propensity vector "(d+(b-d)*N/K)*N", "alpha/(1+w*N)*N*P", "c*alpha/(1+w*N)*N*P", "g*P") tf <- 10 simName <- "Rosenzweig-MacArthur predator-prey model" set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="D",simName,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="ETL",simName,tau=0.01,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="BTL",simName,verbose=FALSE,consoleInterval=1) set.seed(1) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="OTL",simName,verbose=FALSE,consoleInterval=1) #=============================================================================== # Kermack-McKendric SIR model (Brown & Rothery, 1993) #=============================================================================== cat("Running 8/8...\n") flush.console() # Define parameters parms <- c(beta=.001, gamma=.100) # Define system x0 <- c(S=500, I=1, R=0) # Initial state vector nu <- matrix(c(-1,0,1,-1,0,1),nrow=3,byrow=T) # State-change matrix a <- c("beta*S*I", "gamma*I") # Propensity vector tf <- 100 # Final time simName <- "Kermack-McKendrick SIR" set.seed(2) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="D",simName,verbose=FALSE,consoleInterval=1) set.seed(2) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="ETL",simName,verbose=FALSE,consoleInterval=1) set.seed(2) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="BTL",simName,verbose=FALSE,consoleInterval=1) set.seed(2) if (!dry.run) out <- ssa(x0,a,nu,parms,tf,method="OTL",simName,verbose=FALSE,consoleInterval=1) Rprof(NULL) #ReadProfileData("benchmark.Rprof") }