## ---- include = FALSE--------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", echo = TRUE, hold = TRUE, fig.width = 7, fig.height = 10.5, eval = TRUE ) # set seed for reproduibility set.seed(seed = 42) ## ---- eval=FALSE-------------------------------------------------------------- # #################### # # Load libraries # #################### # library(MGDrivE) # # #################### # # Output Folder # #################### # outFolder <- "mgdrive" # dir.create(path = outFolder) # # #################### # # Simulation Parameters # #################### # # days to run the simulation # tMax <- 365 # # # number of Monte Carlo iterations # nRep <- 5 # # # each Monte Carlo iteration gets its own folder # folderNames <- file.path(outFolder, # formatC(x = 1:nRep, width = 3, format = "d", flag = "0")) # # # entomological parameters # bioParameters <- list(betaK=20,tEgg=5,tLarva=6,tPupa=4,popGrowth=1.175,muAd=0.09) # # # a 5-node network with 5% per day migration rate # sitesNumber <- 5 # adultPopEquilibrium <- 500 # # # auxiliary function # triDiag <- function(upper, lower){ # # # return matrix # retMat <- matrix(data = 0, nrow = length(upper) + 1, ncol = length(upper) + 1) # # # set index values for upper/lower triangles # indx <- 1:length(upper) # # # set forward/backward migration using matrix access # retMat[cbind(indx+1,indx)] <- lower # retMat[cbind(indx,indx+1)] <- upper # # # set stay probs # diag(x = retMat) <- 1-rowSums(x = retMat) # # return(retMat) # } # # # fill movement matrix # # Remember, rows need to sum to 1. # moveMat <- triDiag(upper = rep.int(x = 0.05, times = sitesNumber-1), # lower = rep.int(x = 0.05, times = sitesNumber-1)) # # # batch migration is disabled by setting the probability to 0 # batchMigration <- basicBatchMigration(batchProbs=0, # sexProbs=c(.5,.5), # numPatches=sitesNumber) # # #################### # # Basic Inheritance pattern # #################### # # Mendelian cube with standard (default) paraameters # cube <- cubeMendelian() # # # #################### # # Setup releases and batch migration # #################### # # set up the empty release vector # # MGDrivE pulls things out by name # patchReleases <- replicate(n=sitesNumber, # expr={list(maleReleases=NULL,femaleReleases=NULL, # eggReleases=NULL,matedFemaleReleases=NULL)}, # simplify=FALSE) # # # choose release parameters # # Releases start at time 25, occur once a week, for 10 weeks. # # There are 100 mosquitoes released every time. # releasesParameters <- list(releasesStart=25, # releasesNumber=10, # releasesInterval=7, # releaseProportion=100) # # # generate male release vector # maleReleasesVector <- generateReleaseVector(driveCube=cube, # releasesParameters=releasesParameters) # # # generate female release vector # femaleReleasesVector <- generateReleaseVector(driveCube=cube, # releasesParameters=releasesParameters) # # # put releases into the proper place in the release list # # This performs the releases in the first patch only # patchReleases[[1]]$maleReleases <- maleReleasesVector # patchReleases[[1]]$femaleReleases <- femaleReleasesVector # # # #################### # # Combine parameters and run! # #################### # # setup parameters for the network. This builds a list of parameters required for # # every population in the network. # netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber, # beta=bioParameters$betaK, muAd=bioParameters$muAd, # popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, # tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, # AdPopEQ=adultPopEquilibrium, inheritanceCube = cube) # # # set MGDrivE to run stochastic # setupMGDrivE(stochasticityON = TRUE, verbose = FALSE) # # # build network prior to run # MGDrivESim <- Network$new(params=netPar, # driveCube=cube, # patchReleases=patchReleases, # migrationMale=moveMat, # migrationFemale=moveMat, # migrationBatch=batchMigration, # directory=folderNames, # verbose = FALSE) # # run simulation # MGDrivESim$multRun(verbose = FALSE) # # # #################### # # Post Analysis # #################### # # First, split output by patch # # Second, aggregate females by their mate choice # for(i in 1:nRep){ # splitOutput(readDir = folderNames[i], remFile = TRUE, verbose = FALSE) # aggregateFemales(readDir = folderNames[i], genotypes = cube$genotypesID, # remFile = TRUE, verbose = FALSE) # } # # # plot output of first run to see effect # plotMGDrivESingle(readDir=folderNames[1],totalPop = TRUE,lwd=3.5,alpha=1) # # # plot all 5 repetitions together # plotMGDrivEMult(readDir=outFolder,lwd=0.35,alpha=0.75) ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) ## ----------------------------------------------------------------------------- #################### # Output Folder #################### outFolder <- "mgdrive" dir.create(path = outFolder) ## ----------------------------------------------------------------------------- #################### # Simulation Parameters #################### # days to run the simulation tMax <- 365 # number of Monte Carlo iterations nRep <- 5 # each Monte Carlo iteration gets its own folder folderNames <- file.path(outFolder, formatC(x = 1:nRep, width = 3, format = "d", flag = "0")) folderNames ## ----------------------------------------------------------------------------- # entomological parameters bioParameters <- list(betaK=20,tEgg=5,tLarva=6,tPupa=4,popGrowth=1.175,muAd=0.09) ## ----------------------------------------------------------------------------- # a 5-node network with 5% per day migration rate sitesNumber <- 5 adultPopEquilibrium <- 500 patchPops <- rep(adultPopEquilibrium,sitesNumber) # patchPops is optional. If all populations are the same size, parameterizeMGDrivE # can take a single number. However, if you desire different population sizes, it # must be a vector of length equal to the number of sites # landscape # auxiliary function triDiag <- function(upper, lower){ # return matrix retMat <- matrix(data = 0, nrow = length(upper) + 1, ncol = length(upper) + 1) # set index values for upper/lower triangles indx <- 1:length(upper) # set forward/backward migration using matrix access retMat[cbind(indx+1,indx)] <- lower retMat[cbind(indx,indx+1)] <- upper # set stay probs diag(x = retMat) <- 1-rowSums(x = retMat) return(retMat) } # fill movement matrix # Remember, rows need to sum to 1. moveMat <- triDiag(upper = rep.int(x = 0.05, times = sitesNumber-1), lower = rep.int(x = 0.05, times = sitesNumber-1)) moveMat ## ----------------------------------------------------------------------------- # batch migration is disabled by setting the probability to 0 batchMigration <- basicBatchMigration(batchProbs=0, numPatches=sitesNumber) ## ----------------------------------------------------------------------------- #################### # Basic Inheritance pattern #################### # Mendelian cube with standard (default) paraameters cube <- cubeMendelian() ## ----------------------------------------------------------------------------- #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 25, occur once a week, for 10 weeks. # There are 100 mosquitoes released every time. releasesParameters <- list(releasesStart=25, releasesNumber=10, releasesInterval=7, releaseProportion=100) ## ----------------------------------------------------------------------------- # generate male release vector maleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) maleReleasesVector[[1]] ## ----------------------------------------------------------------------------- # generate female release vector femaleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- maleReleasesVector patchReleases[[1]]$femaleReleases <- femaleReleasesVector ## ----------------------------------------------------------------------------- #################### # Combine parameters and run! #################### # setup parameters for the network. This builds a list of parameters required for # every population in the network. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=patchPops, inheritanceCube = cube) ## ----------------------------------------------------------------------------- # set MGDrivE to run stochastic setupMGDrivE(stochasticityON = TRUE) ## ----------------------------------------------------------------------------- # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=folderNames, verbose = TRUE) # list folders to show that they have been created list.files(path = outFolder) ## ----------------------------------------------------------------------------- # run simulation MGDrivESim$multRun(verbose = FALSE) ## ----------------------------------------------------------------------------- # first and last repetitions list.files(path = outFolder)[1] list.files(path = list.files(path = outFolder, full.names = TRUE)[1], recursive = TRUE) list.files(path = outFolder)[5] list.files(path = list.files(path = outFolder, full.names = TRUE)[5], recursive = TRUE) ## ----------------------------------------------------------------------------- # read in male and female files mMat <- as.matrix(read.csv(file = list.files(path = outFolder, full.names = TRUE, recursive = TRUE)[2], header = TRUE, sep = ",")) fMat <- as.matrix(read.csv(file = list.files(path = outFolder, full.names = TRUE, recursive = TRUE)[1], header = TRUE, sep = ",")) # look at male file header colnames(mMat) ## ----------------------------------------------------------------------------- head(x = mMat, n = 2*sitesNumber) ## ----------------------------------------------------------------------------- # look at female file header colnames(fMat) ## ----------------------------------------------------------------------------- head(x = fMat, n = 2*sitesNumber) ## ----------------------------------------------------------------------------- # First, split output by patch for(i in 1:nRep){ splitOutput(readDir = folderNames[i], remFile = TRUE, verbose = FALSE) } ## ----------------------------------------------------------------------------- # first and last repetitions list.files(path = outFolder)[1] list.files(path = list.files(path = outFolder, full.names = TRUE)[1], recursive = TRUE) list.files(path = outFolder)[5] list.files(path = list.files(path = outFolder, full.names = TRUE)[5], recursive = TRUE) # read in examples of new male/female files twoFiles <- list.files(path = outFolder, full.names = TRUE, recursive = TRUE)[c(1,6)] # read in male and female files mMat <- as.matrix(read.csv(file = twoFiles[2], header = TRUE, sep = ",")) fMat <- as.matrix(read.csv(file = twoFiles[1], header = TRUE, sep = ",")) ## ----------------------------------------------------------------------------- head(x = mMat, n = 5) ## ----------------------------------------------------------------------------- head(x = fMat, n = 5) ## ----------------------------------------------------------------------------- # Second, aggregate females by their mate choice for(i in 1:nRep){ aggregateFemales(readDir = folderNames[i], genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) } ## ----------------------------------------------------------------------------- # changed female files in first repetition list.files(path = outFolder)[1] list.files(path = list.files(path = outFolder, full.names = TRUE)[1], recursive = TRUE) # read in examples of new female file fMat2 <- as.matrix(read.csv(file = list.files(path = outFolder, recursive = TRUE, full.names = TRUE)[1], header = TRUE, sep = ",")) ## ----------------------------------------------------------------------------- head(x = fMat2, n = 5) ## ----------------------------------------------------------------------------- # show non-aggregated female file split by patch # This is for patch number 1 fMat[(releasesParameters$releasesStart-2):(releasesParameters$releasesStart+2), ] cat("\n") # show aggregated female file # This is for patch number 1 fMat2[(releasesParameters$releasesStart-2):(releasesParameters$releasesStart+2), ] ## ----------------------------------------------------------------------------- # plot the first repetition plotMGDrivESingle(readDir=folderNames[1],totalPop = TRUE,lwd=3.5,alpha=1) ## ----------------------------------------------------------------------------- # plot all 5 repetitions together plotMGDrivEMult(readDir=outFolder,lwd=0.35,alpha=0.75) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls())