Nothing
###############################################################################
# __ _____________ _ ______
# / |/ / ____/ __ \_____(_) __/ ____/
# / /|_/ / / __/ / / / ___/ / | / / __/
# / / / / /_/ / /_/ / / / /| |/ / /___
# /_/ /_/\____/_____/_/ /_/ |___/_____/
#
# MGDrivE Auxiliary
# Marshall Lab
# jared_bennett@berkeley.edu
# December 2019
#
###############################################################################
#######################################
# Calculate Parameter Values
#######################################
#' Solve for Omega (additional genotype-specific mortality)
#'
#' Solves for root of equation of geometrically-distributed lifespan for value of omega.
#'
#' @param mu Daily mortality probability (discrete-time hazard, called \code{muAd} in code)
#' @param lifespanReduction Target reduced lifespan, between 0 and 1
#' (target average lifespan will be \eqn{\frac{1}{\mu_{Ad}} \times lifespanReduction})
#'
#' @importFrom stats uniroot
#'
#' @examples
#' # reduce lifespan by 10%
#' # Example mu is an average for Aedes
#' newOmega <- calcOmega(mu = 0.11, lifespanReduction = 0.90)
#'
#' @export
calcOmega <- function(mu, lifespanReduction){
# safety check
if(lifespanReduction > 1 | lifespanReduction < 0){
stop("parameter 'lifespanReduction' must be in interval [0,1]")
}
lifespan <- (1/mu) * lifespanReduction
calcOmega_f <- function(omega, mu, lifespan){
(1/(1-omega+(omega*mu)))-lifespan
}
return(uniroot(f=calcOmega_f, interval=c(0,1), lifespan=lifespan,
mu=mu,maxiter=1e4)$root)
}
#######################################
# Delete files in a directory
#######################################
#' Erase all files in a directory
#'
#' Given a directory path, check that it exists, and if so, delete all its contents.
#'
#' @param directory Directory whose contents will be deleted
#' @param verbose Chatty? Default is TRUE
#'
#' @examples
#' \dontrun{
#' # Path to directory, can tilde expand
#' myPath <- "~/path/to/write/output"
#'
#' # Erase directory
#' # No return value
#' eraseDirectory(directory = myPath)
#' }
#'
#' @export
eraseDirectory <- function(directory, verbose = TRUE){
# check directory exists
if(!dir.exists(directory)){
if(verbose){cat("no such directory exists\n")}
}
# list files/subdirectories to delete
dirFiles = list.files(path = directory)
# delete contents
if(length(dirFiles)>0){
for(i in dirFiles){
if(verbose){cat("removing file: ",file.path(directory, i),"\n", sep = "")}
unlink(x = file.path(directory, i), recursive = TRUE)
}
}
}
#######################################
# Population Objects
#######################################
#' Normalise a Numeric Vector
#'
#' Normalise a numeric vector to sum to one
#'
#' @param vector Numeric vector
#'
normalise <- function(vector){
if(all(vector==0)){
return(vector)
} else {
return(vector/sum(vector))
}
}
#######################################
# Kernel-related
#######################################
#' Calculates the zero-inflation part of a hurdle exponential kernel.
#'
#' Given the probability of an adult mosquito to stay in the same patch throughout
#' its whole lifespan, and its mortality, it calculates the height of the pulse-density
#' part of the hurdle kernel.
#'
#' @param stayThroughLifespanProbability Probability of a mosquito to spend its
#' whole lifespan in the same node
#' @param adultMortality Adult mortality rate
#'
#' @examples
#' # setup distance matrix
#' # two-column matrix with latitude/longitude, in degrees
#' latLong = cbind(runif(n = 5, min = 0, max = 90),
#' runif(n = 5, min = 0, max = 180))
#'
#' # Vincenty Ellipsoid distance formula
#' distMat = calcVinEll(latLongs = latLong)
#'
#' # get hurdle height
#' # Lets assume 80% stay probs and adult mortality of 0.1
#' hHeight <- calcZeroInflation(stayThroughLifespanProbability = 0.80,
#' adultMortality = 0.1)
#'
#' # calculate hurdle exponential distribution over distances
#' kernMat = calcHurdleExpKernel(distMat = distMat, rate = 10, p0 = hHeight)
#'
#' @export
calcZeroInflation <- function(stayThroughLifespanProbability,adultMortality){
return(stayThroughLifespanProbability^(adultMortality))
}
#######################################
# Post-processing of Output
#######################################
#' Split Output by Patch
#'
#' Split output into multiple files by patches.
#'
#' @param readDir Directory where output was written to
#' @param writeDir Directory to write output to. Default is readDir
#' @param remFile Remove original output? Default is TRUE
#' @param verbose Chatty? Default is TRUE
#'
#' @importFrom utils setTxtProgressBar txtProgressBar write.table
#'
#' @examples
#' \dontrun{
#' # This example assumes user has already run MGDrivE and generated output.
#' # If that's untree, see vignette for complete example
#' fPath <- "path/to/data/containing/folder"
#' oPath <- "path/to/write/output"
#'
#' # split data by patch, keep original files
#' # no return value
#' splitOutput(readDir = fPath, writeDir = oPath, remFile = FALSE)
#'
#' # Alternatively, remove the original files and write new ones in their place
#' fPath <- "path/to/data/containing/folder"
#'
#' splitOutput(readDir = fPath, writeDir = NULL, remFile = TRUE)
#' }
#'
#' @export
splitOutput <- function(readDir, writeDir=NULL, remFile=TRUE, verbose=TRUE){
# get all files in the read directory
dirFiles = list.files(path = readDir, pattern = "^[MF]_.*\\.csv$", full.names = TRUE)
# check write directory
if(is.null(writeDir)){writeDir <- readDir}
# initialize text, progress bar below
if(verbose){
# what we're doing
cat(" Splitting", length(dirFiles), "files.\n")
# whether or not we remove files
if(remFile){
cat(" Removing original files.\n")
} else {
cat(" Not removing original files.\n\n")
}
}
# loop over files
for(selectFile in dirFiles){
# Read in files
columnNames <- scan(file = selectFile, what = character(), sep = ",",
nlines = 1, quiet = TRUE)
fileIn = matrix(data = scan(file = selectFile, what = numeric(), sep = ",",
skip = 1, quiet = TRUE),
ncol = length(columnNames), byrow = TRUE,
dimnames = list(NULL, columnNames))
# get stuff for loop
patches <- unique(fileIn[ ,2])
patchRows <- seq.int(from = 1, to = dim(fileIn)[1], by = length(patches))
fileBase <- strsplit(x = basename(path = selectFile), split = ".", fixed = TRUE)[[1]][1]
# progress bar stuff
if(verbose){
pb = txtProgressBar(min = 0,max = length(patches),style = 3)
}
pbVal = 0
# loop over patches, split and write out
for(patch in patches){
# output file name
patchName = file.path(writeDir, paste0(fileBase, "_Patch",
formatC(x = patch, width = 3, format = "d", flag = "0"),
".csv") )
# write output
# subset matrix by patch number, then remove patch label column
write.table(x = fileIn[patchRows, ][ ,-2], file = patchName, sep = ",",
row.names = FALSE, col.names = TRUE)
# increment rows
patchRows[] <- patchRows + 1L
# some indication that it's working
pbVal = pbVal+1
if(verbose){setTxtProgressBar(pb = pb, value = pbVal)}
}
# check if removing original output
if(remFile){file.remove(selectFile)}
} # end loop over files
}
#' Aggregate Female Output by Genotype
#'
#' Aggregate over male mate genotype to convert female matrix output into vector output.
#'
#' @param readDir Directory to read input from
#' @param writeDir Directory to write output to. Default is readDir
#' @param genotypes Character vector of possible genotypes; found in \code{driveCube$genotypesID}
#' @param remFile Boolean flag to remove original (unaggregated) file
#' @param verbose Chatty? Default is TRUE
#'
#'
#' @examples
#' \dontrun{
#' # This example assumes user has already run MGDrivE and generated output.
#' # This also assumes that the user has already split output by patch.
#' # See vignette for complete example.
#'
#' # set read/write directory
#' fPath <- "path/to/data/containing/folder"
#'
#' # Need genotypes from the cube run in the simulation
#' # This is dependent on the simulation run
#' # Using Mendelian cube for this example
#' cube <- cubeMendelian()
#'
#' # no return value from function
#' aggregateFemales(readDir= fPath, writeDir = NULL, genotypes = cube$genotypesID,
#' remFile = TRUE)
#' }
#'
#' @export
aggregateFemales <- function(readDir, writeDir=NULL, genotypes, remFile=TRUE,
verbose=TRUE){
# check write directory
if(is.null(writeDir)){writeDir <- readDir}
#get female files
femaleFiles = list.files(path = readDir, pattern = "^F_.*\\.csv$", full.names = TRUE)
# initialize progress bar and text
if(verbose){
# what we're doing
cat(" Aggregating", length(femaleFiles), "files.\n")
# whether or not we remove shit
if(remFile){
cat(" Removing original files.\n")
} else {
cat(" Not removing original files.\n\n")
}
# progress bar
pb = txtProgressBar(min = 0,max = length(femaleFiles),style = 3)
}
pbVal = 0
# setup aggregation columns
# +1 offset for time column
nGenos <- length(genotypes)
idxList <- lapply(X = 0:(nGenos-1) * nGenos, FUN = "+", 2:(nGenos+1))
#loop over files
for(selectFile in femaleFiles){
#read in file
thisOutput = matrix(data = scan(file = selectFile, what = numeric(), sep = ",", skip = 1, quiet = TRUE),
ncol = 1+nGenos^2, byrow = TRUE)
#aggregate output
aggregateOut = vapply(X = idxList, FUN = function(x){
rowSums(x = thisOutput[ ,x])
}, FUN.VALUE = numeric(nrow(thisOutput)))
#name file and write out
fileName <- file.path(writeDir, sub(pattern = "_",replacement = "_Aggregate_",
x = basename(path = selectFile), fixed = TRUE)
)
# combine time column from original file, with aggregated matrix of female counts
write.table(x = cbind(thisOutput[ ,1],aggregateOut), file = fileName,
sep = ",", row.names = FALSE, col.names = c("Time",genotypes))
# check if removing file
if(remFile){file.remove(selectFile)}
# some indication that it's working
pbVal = pbVal+1
if(verbose){setTxtProgressBar(pb = pb, value = pbVal)}
}#end file loop
}#end program
#' Aggregate Output Over Landscape
#'
#' This function aggregates the output of a run over the entire output, i.e., all
#' of the patches. It writes the output one level above the folder pointed to by
#' readDir, if writeDir is NULL. Output consists of 2 csv files, one for males and
#' one for females, "...M_LandscapeAgg_Run...csv".
#'
#' @usage aggregateOutput(readDir, writeDir=NULL)
#'
#' @param readDir Directory where output was written to
#' @param writeDir Directory to write output to. Default is one level above readDir
#'
#'
#' @examples
#' \dontrun{
#' # This assumes user has run MGDrivE and output is in fPath.
#' # See vignette for examples on how to run MGDrivE
#'
#' # read/write dirs
#' fPath <- "folder/containing/output"
#' oPath <- "folder/to/write/stuff"
#'
#' # first, split output by patch and aggregate females by mate genotype
#' # remember, cube is for example and changes with simulation
#' # landscape aggregation will work if females are not aggregated, but it's slower
#' cube <- cubeMendelian()
#'
#' splitOutput(readDir = fPath, writeDir = NULL, remFile = TRUE)
#' aggregateFemales(readDir= fPath, writeDi = NULL, genotypes = cube$genotypesID,
#' remFile = TRUE)
#'
#' # aggregate mosquitoes over entire landscape
#' # no return value
#' aggregateOutput(readDir = fPath, writeDir = NULL)
#' }
#'
#' @export
aggregateOutput <- function(readDir, writeDir=NULL){
# get all files in the read directory
mFiles = list.files(path = readDir, pattern = "^M_.*\\.csv$", full.names = TRUE)
fFiles = list.files(path = readDir, pattern = "^F_.*\\.csv$", full.names = TRUE)
# check write directory
hold <- strsplit(x = readDir, split = "/", fixed = TRUE)[[1]]
if(is.null(writeDir)){
writeDir <- paste0(hold[-length(hold)], collapse = "/")
}
# set object sizes for males and females
columnNames <- scan(file = mFiles[1], what = character(), sep = ",",
quiet = TRUE, nlines = 1)
mMat <- matrix(data = scan(file = mFiles[1], what = integer(), sep = ",", skip = 1, quiet = TRUE),
ncol = length(columnNames), byrow = TRUE, dimnames = list(NULL,columnNames))
columnNames <- scan(file = fFiles[1], what = character(), sep = ",",
quiet = TRUE, nlines = 1)
fMat <- matrix(data = scan(file = fFiles[1], what = integer(), sep = ",", skip = 1, quiet = TRUE),
ncol = length(columnNames), byrow = TRUE, dimnames = list(NULL,columnNames))
# set things to reuse
nRowM <- dim(mMat)[1]
nColM <- dim(mMat)[2]
nRowF <- dim(fMat)[1]
nColF <- dim(fMat)[2]
nReadM <- nRowM * nColM
nReadF <- nRowF * nColF
simTime <- mMat[ ,"Time"]
# aggregate files
for(sFile in 2:length(mFiles)){
# males
mMat <- mMat + matrix(data = scan(file = mFiles[sFile], what = integer(),
sep = ",", skip = 1, n = nReadM, quiet = TRUE),
nrow = nRowM, ncol = nColM, byrow = TRUE)
# females
fMat <- fMat + matrix(data = scan(file = fFiles[sFile], what = integer(),
sep = ",", skip = 1, n = nReadF, quiet = TRUE),
nrow = nRowF, ncol = nColF, byrow = TRUE)
} # end loop over files
# fix time addition
mMat[ ,1] <- fMat[ ,1] <- simTime
# print males
fileName <- file.path(writeDir, file.path("M_LandscapeAgg_Run", hold[length(hold)], ".csv", fsep = ""))
write.table(x = mMat, file = fileName, sep = ",", row.names = FALSE, col.names = TRUE)
# print females
fileName <- file.path(writeDir, file.path("F_LandscapeAgg_Run", hold[length(hold)], ".csv", fsep = ""))
write.table(x = fMat, file = fileName, sep = ",", row.names = FALSE, col.names = TRUE)
} # end aggregate over landscape
#' Retrieve Output
#'
#' Read in output from directory. The resulting object will be a nested list;
#' outermost nesting dimension indexes runID, within runID elements are split by sex
#' and innermost nesting is over patches.
#'
#' @param readDir Directory where output was written to; must not end in path separator
#' @param verbose Chatty? Default is TRUE
#'
#' @importFrom utils read.csv
#' @importFrom stats setNames
#'
#' @examples
#' \dontrun{
#' # Example assumes user has run and analyzed MGDrivE.
#' # See vignette for examples of how to do that.
#'
#' # set read directory
#' fPath <- "path/to/split/aggregated/output"
#'
#' # read in data as nested lists
#' dataList <- retrieveOutput(readDir = fPath)
#' }
#'
#' @return Nested List
#'
#' @export
retrieveOutput <- function(readDir, verbose=TRUE){
dirFiles = list.files(path = readDir)
runID = strsplit(x = dirFiles,split = "_")
runID = unique(x = grep( pattern = "Run", x = unlist(x = runID), value = TRUE))
output = setNames(object = vector(mode = "list",length = length(runID)), runID)
for(run in runID){
if(verbose){cat("processing ",run,"\n",sep="")}
runFiles = dirFiles[grep(pattern = run,x = dirFiles)]
patches = regmatches(x = runFiles, m = regexpr(pattern = "Patch[0-9]+", text = runFiles))
patches = sort(unique(patches))
output[[run]]$F = setNames(object = vector(mode = "list",length = length(patches)),
nm = patches)
output[[run]]$M = setNames(object = vector(mode = "list",length = length(patches)),
nm = patches)
# retrieve males for this run
males = runFiles[grep(pattern = "M_",x = runFiles, fixed = TRUE, useBytes = TRUE)]
for(male in males){
thisPatch = regmatches(x = male, m = regexpr(pattern = "Patch[0-9]+", text = male))
thisOutput = read.csv(file = file.path(readDir, male))
thisOutput = as.matrix(thisOutput)
rownames(thisOutput) = NULL
output[[run]]$M[[thisPatch]] = thisOutput[,-1]
}
# retrieve females for this run
females = runFiles[grep(pattern = "F_Aggregate",x = runFiles)]
for(female in females){
thisPatch = regmatches(x = female, m = regexpr(pattern = "Patch[0-9]+", text = female))
thisOutput = read.csv(file = file.path(readDir, female))
thisOutput = thisOutput[,-1]
output[[run]]$F[[thisPatch]] = as.matrix(thisOutput)
}
}
return(output)
}
#' Summary Statistics for Stochastic MGDrivE
#'
#' This function reads in all repetitions for each patch and calculates either
#' the mean, quantiles, or both. User chooses the quantiles, up to 4 decimal places,
#' and enters them as a vector. Quantiles are calculated empirically. (order does not matter) \cr
#'
#' Given the readDir, this function assumes the follow file structure: \cr
#' * readDir
#' * repetition 1
#' * patch 1
#' * patch 2
#' * patch 3
#' * repetition 2
#' * patch 1
#' * patch 2
#' * patch 3
#' * repetition 3
#' * repetition 4
#' * ... \cr
#'
#' Output files are *.csv contain the mean or quantile in the file name, i.e.
#' {M/F}_Mean_(patchNum).csv and {M/F}_Quantile_(quantNum)_(patchNum).csv.
#'
#' @param readDir Directory to find repetition folders in
#' @param writeDir Directory to write output
#' @param mean Boolean, calculate mean or not. Default is TRUE
#' @param quantiles Vector of quantiles to calculate. Default is NULL
#' @param verbose Chatty? Default is TRUE
#'
#' @examples
#' \dontrun{
#' # This function assumes network$multRun() has been performed, or several
#' # network$oneRun() have been performed and all of the data has been split
#' # and aggregated.
#'
#' # read/write paths
#' fPath <- "path/to/folder/ofFolders/with/data"
#' oPath <- "my/path/output"
#'
#' # here, only calculate mean, no quantiles
#' # no return value
#' calcQuantiles(readDir = fPath, writeDir = oPath, mean = TRUE,
#' quantiles = NULL)
#'
#' # here, calculate 2.5% and 97.5% quantiles
#' calcQuantiles(readDir = fPath, writeDir = oPath, mean = FALSE,
#' quantiles = c(0.025, 0.975))
#' }
#'
#' @return Writes output to files in writeDir
#' @export
calcQuantiles <- function(readDir, writeDir, mean=TRUE, quantiles=NULL,
verbose=TRUE){
# safety check
if(!mean && is.null(quantiles)){
stop("User needs to specify the mean or which quantiles to calculate. ")
}
# get all files
patchFiles <- lapply(X = list.dirs(path = readDir, full.names = TRUE, recursive = FALSE),
FUN = list.files, pattern = "^[FM].*\\.csv$", full.names = TRUE)
#s ubset females/males
malePatches <- lapply(X = patchFiles, FUN = grep, pattern = "M_", fixed = TRUE, value=TRUE)
femalePatches <- lapply(X = patchFiles, FUN = grep, pattern = "F_Aggregate", fixed = TRUE, value=TRUE)
# generate a list of all patches to run over
patchList = unique(regmatches(x = patchFiles[[1]],
m = regexpr(pattern = "Patch[0-9]+",
text = patchFiles[[1]],
perl = TRUE)))
# read in a file initially to get variables and setup return array
testFile <- read.csv(file = patchFiles[[1]][1], header = TRUE)
# bunch of constants that get used several times
numReps <- length(patchFiles)
columnNames <- names(testFile)
numRow <- dim(testFile)[1]
numCol <- dim(testFile)[2]
numRead <- numRow*numCol
# setup input data holder
popDataMale <- array(data = 0, dim = c(numRow, numReps, numCol-1))
popDataFemale <- array(data = 0, dim = c(numRow, numReps, numCol-1))
# setup output data holder
if(is.null(quantiles)){
outputDataMale <- array(data = 0, dim = c(numRow, numCol, 1),
dimnames = list(NULL, columnNames, NULL))
outputDataFemale <- array(data = 0, dim = c(numRow, numCol, 1),
dimnames = list(NULL, columnNames, NULL))
} else {
outputDataMale <- array(data = 0, dim = c(numRow, numCol, length(quantiles)),
dimnames = list(NULL, columnNames, NULL))
outputDataFemale <- array(data = 0, dim = c(numRow, numCol, length(quantiles)),
dimnames = list(NULL, columnNames, NULL))
}
# set time from testFile
outputDataMale[,1,] <- testFile[ ,"Time"]
outputDataFemale[,1,] <- testFile[ ,"Time"]
# write initial outputs and setup progress bar
if(verbose){
# how many patches and repetitions are we running through
cat(" Patches:", length(patchList), "\n Repetitions:", numReps, "\n\n")
# setup progress bar
pb = txtProgressBar(min = 0,max = length(patchList),style = 3)
}
pbVal = 0
# loop over all patches and do stats.
for(patch in patchList){
# get male and female files, all repetitions of this patch
maleFiles <- vapply(X = malePatches,
FUN = grep, pattern = patch, fixed=TRUE, value=TRUE,
FUN.VALUE = character(length = 1L))
femaleFiles <- vapply(X = femalePatches,
FUN = grep, pattern = patch, fixed=TRUE, value=TRUE,
FUN.VALUE = character(length = 1L))
# read in all repetitions for this patch
# remove time column after reading in.
for(repetition in 1:numReps){
popDataMale[ ,repetition, ] <- matrix(data = scan(file = maleFiles[repetition],
what = numeric(), n = numRead,
sep = ",", skip = 1, quiet = TRUE),
nrow = numRow, ncol = numCol, byrow = TRUE)[ ,-1]
popDataFemale[ ,repetition, ] <- matrix(data = scan(file = femaleFiles[repetition],
what = numeric(), n = numRead,
sep = ",", skip = 1, quiet = TRUE),
nrow = numRow, ncol = numCol, byrow = TRUE)[ ,-1]
}
# if the user wants the mean, calculate and put out.
if(mean){
for(whichCol in 1:(numCol-1)){
outputDataMale[ ,whichCol+1,1] <- .rowMeans(x = popDataMale[ , ,whichCol],
m = numRow, n = numReps)
outputDataFemale[ ,whichCol+1,1] <- .rowMeans(x = popDataFemale[ , ,whichCol],
m = numRow, n = numReps)
}
# write output
maleFileName <- file.path(writeDir,
file.path("M_Mean_", patch, ".csv", fsep = "") )
femaleFileName <- file.path(writeDir,
file.path("F_Mean_", patch, ".csv", fsep = "") )
write.table(x = outputDataMale[ , ,1], file = maleFileName, sep = ",",
row.names = FALSE, col.names =TRUE)
write.table(x = outputDataFemale[ , ,1], file = femaleFileName, sep = ",",
row.names = FALSE, col.names =TRUE)
} # end mean
# if the user wants quantiles, do them and write output
if(!is.null(quantiles)){
for(whichCol in 1:(numCol-1)){
outputDataMale[ ,whichCol+1, ] <- quantileC(Trials = popDataMale[ , ,whichCol],
Probs = quantiles)
outputDataFemale[ ,whichCol+1, ] <- quantileC(Trials = popDataFemale[ , ,whichCol],
Probs = quantiles)
}#end loop to calculate quantiles
# write output
for(whichQuant in 1:length(quantiles)){
# file names
maleFileName <- file.path(writeDir,
file.path("M_Quantile_",
formatC(x = quantiles[whichQuant], digits = 4,
format = "f", decimal.mark = "",
big.mark = NULL),
"_", patch, ".csv", fsep = "") )
femaleFileName <- file.path(writeDir,
file.path("F_Quantile_",
formatC(x = quantiles[whichQuant], digits = 4,
format = "f", decimal.mark = "",
big.mark = NULL),
"_", patch, ".csv", fsep = "") )
# write output
write.table(x = outputDataMale[ , ,whichQuant], file = maleFileName, sep = ",",
row.names = FALSE, col.names =TRUE)
write.table(x = outputDataFemale[ , ,whichQuant], file = femaleFileName, sep = ",",
row.names = FALSE, col.names =TRUE)
} # end loop to write files
} # end quantiles
# some indication that it's working
pbVal = pbVal+1
if(verbose){setTxtProgressBar(pb = pb, value = pbVal)}
} # end loop over patches
} # end function
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.