--- title: "2 - Running a numerical experimental design" data: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{2 - Running a numerical experimental design} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` This is an example of how to run a numerical experimental design with *landsepi*. ```{r, eval=FALSE, message=FALSE} library(landsepi) ``` ## Write the experimental design First create a dataframe containing, for each line, the value of target parameters (i.e. those destined to vary from one simulation to another). The dataframe must also contain columns to store simulation outputs. >For the sake of simplicity, the following numerical design contains only 5 simulations for which 7 target parameters have been arbitrarily chosen to vary and 3 output variable are to be computed. Off course, real numerical designs will include more simulations, in particular those representing factorial experimental plans. ```{r} myDesign <- data.frame(nTSpY = c(120, 110, 100, 90, 80) , pI0 = c(0, 1E-4, 5E-4, 1E-3, 1E-2) , infection_rate = c(0.5, 0.4, 0.3, 0.2, 0.1) , id_landscape = 1:5 , aggreg = c(0.07, 0.07, 0.25, 10, 10) , R_efficiency = c(1.00, 0.90, 0.80, 0.70, 0.60) , growth_rate = c(0.1, 0.2, 0.3, 0.1, 0.1) ## create columns to store outputs , durab_MG1 = NA , durab_MG2 = NA , mean_audpc = NA) ``` Then add an identifier for each simulation and specify a random seed (for stochasticity). ```{r} n <- nrow(myDesign) myDesign <- cbind(simul = 1:n, seed = sample(n*100, n), myDesign) myDesign ``` ## Run the simulations Create the object `simul_params` and a repository to store inputs and outputs of the simulations. ```{r, eval=FALSE} simul_params <- createSimulParams(outputDir = getwd()) ``` Then run the experimental design by updating target parameters with the values stored in `myDesign`. Note that some parameters vary jointly (e.g. landscape and dispersal matrix). *See tutorial on how to [run a simple simulation](O1-run_simple_simul.html) for details on the different functions.* ```{r, eval=FALSE} for (i in 1:n){ print(paste("Running simulation", i, "/", n)) ## Set Nyears and nTSpY simul_params <- setTime(simul_params , Nyears = 6 , nTSpY = myDesign$nTSpY[i]) ## update nTSpY ## set seed (to run stochastic replicates) simul_params <- setSeed(simul_params, myDesign$seed[i]) ## update seed ## Pathogen parameters simul_params@ReproSexProb <- logical(0) ## initialize vector of sexual probability (one for each time step) basic_patho_param <- loadPathogen("rust") basic_patho_param$infection_rate <- myDesign$infection_rate[i] ## update inf. rate simul_params <- setPathogen(simul_params, basic_patho_param) ## Landscape parameters landscape <- loadLandscape(myDesign$id_landscape[i]) ## update landscape simul_params <- setLandscape(simul_params, landscape) ## Dispersal parameters disp_patho_clonal <- loadDispersalPathogen(myDesign$id_landscape[i])[[1]] ## update dispersal simul_params <- setDispersalPathogen(simul_params, disp_patho_clonal) ## Genes gene1 <- loadGene(name = "MG 1", type = "majorGene") gene1$mutation_prob<-1e-4 gene2 <- loadGene(name = "MG 2", type = "majorGene") gene2$mutation_prob<-1e-4 genes <- data.frame(rbind(gene1, gene2), stringsAsFactors = FALSE) genes$efficiency <- myDesign$R_efficiency[i] ## update resistance efficiency simul_params <- setGenes(simul_params, genes) ## Cultivars cultivar1 <- loadCultivar(name = "Susceptible", type = "wheat") cultivar2 <- loadCultivar(name = "Resistant1", type = "wheat") cultivar3 <- loadCultivar(name = "Resistant2", type = "wheat") cultivars <- data.frame(rbind(cultivar1, cultivar2, cultivar3) , stringsAsFactors = FALSE) cultivars$growth_rate <- myDesign$growth_rate[i] ## update growth rate simul_params <- setCultivars(simul_params, cultivars) ## Allocate genes to cultivars simul_params <- allocateCultivarGenes(simul_params, "Resistant1", c("MG 1")) simul_params <- allocateCultivarGenes(simul_params, "Resistant2", c("MG 2")) ## Allocate cultivars to croptypes croptypes <- loadCroptypes(simul_params, names = c("Susceptible crop" , "Resistant crop 1" , "Resistant crop 2")) croptypes <- allocateCroptypeCultivars(croptypes, "Susceptible crop", "Susceptible") croptypes <- allocateCroptypeCultivars(croptypes, "Resistant crop 1", "Resistant1") croptypes <- allocateCroptypeCultivars(croptypes, "Resistant crop 2", "Resistant2") simul_params <- setCroptypes(simul_params, croptypes) ## Allocate croptypes to landscape rotation_sequence <- croptypes$croptypeID ## No rotation: 1 rotation_sequence element rotation_period <- 0 ## same croptypes every years prop <- c(1/3, 1/3, 1/3) ## croptypes proportions aggreg <- myDesign$aggreg[i] simul_params <- allocateLandscapeCroptypes(simul_params, rotation_period = rotation_period, rotation_sequence = rotation_sequence, rotation_realloc = FALSE, prop = prop, aggreg = aggreg, graphic = FALSE) ## Initial conditions simul_params <- setInoculum(simul_params, myDesign$pI0[i]) ## update pI0 #Only susceptible hosts are infected, see vignette1 for a more complex inoculum ## configure outputs outputlist <- loadOutputs(epid_outputs = "audpc", evol_outputs = "durability") simul_params <- setOutputs(simul_params, outputlist) ## Check, (save) and run simulation checkSimulParams(simul_params) # simul_params <- saveDeploymentStrategy(simul_params, overwrite = TRUE) res <- runSimul(simul_params, writeTXT=FALSE, graphic = FALSE) ## Extract outputs myDesign$durab_MG1[i] <- res$evol_outputs$durability[,"MG 1"] myDesign$durab_MG2[i] <- res$evol_outputs$durability[,"MG 2"] myDesign$mean_audpc[i] <- mean(res$epid_outputs$audpc$total) ## Create myDesign.txt at first simulation and then append if (i ==1){ write.table(myDesign[i,], paste(simul_params@OutputDir, "/myDesign.txt", sep="") , append=FALSE, row.names = FALSE, col.names = TRUE) } else { write.table(myDesign[i,], paste(simul_params@OutputDir, "/myDesign.txt", sep="") , append=TRUE, row.names = FALSE, col.names = FALSE) } } myDesign ``` ## Compute your own output variables The argument "keepRawResults" from function `runSimul()` allows to keep a set of binary files after the end of the simulation. >This set of binary files is composed of one file per simulated year and per compartment: - H: healthy hosts, - Hjuv: juvenile healthy hosts (for host reproduction), - L: latently infected hosts, - I: infectious hosts, - R: removed hosts, - P: propagules. Each file indicates for every time-step the number of individuals in each field, and when appropriate for each host and pathogen genotypes. ```{r, eval=FALSE} ## Disable computation of outputs: outputlist <- loadOutputs(epid_outputs = "", evol_outputs = "") simul_params <- setOutputs(simul_params, outputlist) ## Run simulation and keep raw binary files: runSimul(simul_params, writeTXT=FALSE, graphic = FALSE, keepRawResults = TRUE) ``` The following code loads the content of the binary files in lists (named "H", "Hjuv", "P", "L", "I", "R"). Each element of these lists contains a matrix or an array of the number of individuals associated with each field, host genotype and pathogen genotype for a given time step (i.e. the number of elements of each list is the total number of time steps). Then, users can compute their own output variables using these lists. ```{r, eval=FALSE} ## retrieve parameters from the object simul_params path <- simul_params@OutputDir Nyears <- simul_params@TimeParam$Nyears nTSpY <- simul_params@TimeParam$nTSpY nTS <- Nyears * nTSpY ## Total number of time-steps Npoly <- nrow(simul_params@Landscape) Nhost <- nrow(simul_params@Cultivars) Npatho <- prod(simul_params@Genes$Nlevels_aggressiveness) ## Initialise lists H <- as.list(1:nTS) Hjuv <- as.list(1:nTS) P <- as.list(1:nTS) L <- as.list(1:nTS) I <- as.list(1:nTS) R <- as.list(1:nTS) index <- 0 ## Read binary files and store values in the lists as matrices or arrays for (year in 1:Nyears) { binfileH <- file(paste(path, sprintf("/H-%02d", year), ".bin", sep = ""), "rb") H.tmp <- readBin(con = binfileH, what = "int", n = Npoly * Nhost * nTSpY, size = 4 , signed = T, endian = "little") close(binfileH) binfileHjuv = file(paste(path, sprintf("/Hjuv-%02d", year), ".bin",sep=""), "rb") Hjuv.tmp <- readBin(con=binfileHjuv, what="int", n=Npoly*Nhost*nTSpY, size = 4 , signed=T,endian="little") close(binfileHjuv) binfileP <- file(paste(path, sprintf("/P-%02d", year), ".bin", sep = ""), "rb") P.tmp <- readBin(con = binfileP, what = "int", n = Npoly * Npatho * nTSpY, size = 4 , signed = T, endian = "little") close(binfileP) binfileL <- file(paste(path, sprintf("/L-%02d", year), ".bin", sep = ""), "rb") L.tmp <- readBin(con = binfileL, what = "int", n = Npoly * Npatho * Nhost * nTSpY , size = 4 , signed = T, endian = "little") close(binfileL) binfileI <- file(paste(path, sprintf("/I-%02d", year), ".bin", sep = ""), "rb") I.tmp <- readBin(con = binfileI, what = "int", n = Npoly * Npatho * Nhost * nTSpY , size = 4 , signed = T, endian = "little") close(binfileI) binfileR <- file(paste(path, sprintf("/R-%02d", year), ".bin", sep = ""), "rb") R.tmp <- readBin(con = binfileR, what = "int", n = Npoly * Npatho * Nhost * nTSpY , size = 4 , signed = T, endian = "little") close(binfileR) ## Convert vectors in matrices or arrays for (t in 1:nTSpY) { H[[t + index]] <- matrix(H.tmp[((Nhost * Npoly) * (t-1)+1):(t * (Nhost * Npoly))] , ncol = Nhost, byrow = T) Hjuv[[t + index]] <- matrix(Hjuv.tmp[((Nhost*Npoly)*(t-1)+1):(t*(Nhost*Npoly))] , ncol=Nhost,byrow=T) P[[t + index]] <- matrix(P.tmp[((Npatho * Npoly) * (t-1)+1):(t * (Npatho * Npoly))] , ncol = Npatho, byrow = T) L[[t + index]] <- array(data = L.tmp[((Npatho * Npoly * Nhost) * (t-1)+1):(t * (Npatho * Npoly * Nhost))] , dim = c(Nhost, Npatho, Npoly)) I[[t + index]] <- array(data = I.tmp[((Npatho * Npoly * Nhost) * (t-1)+1):(t * (Npatho * Npoly * Nhost))] , dim = c(Nhost, Npatho, Npoly)) R[[t + index]] <- array(data = R.tmp[((Npatho * Npoly * Nhost) * (t-1)+1):(t * (Npatho * Npoly * Nhost))] , dim = c(Nhost, Npatho, Npoly)) } index <- index + nTSpY } ```