diff --git a/.gitignore b/.gitignore index c8dd2d6a..0428c571 100644 --- a/.gitignore +++ b/.gitignore @@ -13,6 +13,8 @@ # Output files from R CMD build /*.tar.gz +src/*.o +src/*.so # Output files from R CMD check /*.Rcheck/ diff --git a/DESCRIPTION b/DESCRIPTION index 5ac2bc7d..b2842785 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: SIMplyBee Type: Package Title: 'AlphaSimR' Extension for Simulating Honeybee Populations and Breeding Programmes -Version: 0.4.1 +Version: 0.5.0 Authors@R: c( person("Jana", "Obšteter", email = "obsteter.jana@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-1511-3916")), @@ -26,8 +26,8 @@ URL: https://github.com/HighlanderLab/SIMplyBee License: MIT + file LICENSE Encoding: UTF-8 LazyData: true -Imports: methods, R6, stats, utils, extraDistr (>= 1.9.1), RANN, Rcpp (>= 0.12.7) -Depends: R (>= 3.3.0), AlphaSimR (>= 1.5.3) +Imports: methods, R6, stats, utils, extraDistr (>= 1.9.1), RANN, Rcpp (>= 0.12.7), foreach +Depends: R (>= 3.3.0), AlphaSimR (>= 2.0.0) LinkingTo: Rcpp, RcppArmadillo (>= 0.7.500.0.0), BH RoxygenNote: 7.3.2 Suggests: diff --git a/NAMESPACE b/NAMESPACE index 231169a0..a61e2857 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,6 +2,7 @@ export(SimParamBee) export(addCastePop) +export(addCastePop_internal) export(addDrones) export(addVirginQueens) export(addWorkers) @@ -69,13 +70,11 @@ export(getGv) export(getIbdHaplo) export(getId) export(getLocation) -export(getMisc) export(getPheno) export(getPooledGeno) export(getQtlGeno) export(getQtlHaplo) export(getQueen) -export(getQueenAge) export(getQueenCsdAlleles) export(getQueenCsdGeno) export(getQueenGv) @@ -87,7 +86,6 @@ export(getQueenSegSiteGeno) export(getQueenSegSiteHaplo) export(getQueenSnpGeno) export(getQueenSnpHaplo) -export(getQueenYearOfBirth) export(getSegSiteGeno) export(getSegSiteHaplo) export(getSnpGeno) @@ -192,8 +190,6 @@ export(replaceWorkers) export(resetEvents) export(selectColonies) export(setLocation) -export(setMisc) -export(setQueensYearOfBirth) export(simulateHoneyBeeGenomes) export(split) export(splitPColonyStrength) @@ -207,6 +203,9 @@ import(AlphaSimR) import(Rcpp) importFrom(R6,R6Class) importFrom(extraDistr,rtpois) +importFrom(foreach,"%do%") +importFrom(foreach,"%dopar%") +importFrom(foreach,foreach) importFrom(methods,"slot<-") importFrom(methods,classLabel) importFrom(methods,is) diff --git a/NEWS.md b/NEWS.md index d72a41ac..838297ad 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,33 @@ editor_options: wrap: 72 --- +# SIMplyBee version 0.5.0 + +- 2025-11-27 + +## Major changes + +- swarm/split/supersede do no longer store the year of the queen + +- colonies with high inbreeding that do not produce a viable virgin + queens in max(10, SP\$nVirginQueens) attempts are removed in + swarm/supersede + +- split no longer creates virgin queens in the split colonies but + returns colonies with workers and meta data, but no virgin queens + +- createMultiColony() no longer creates an empty apiary, but it adds + empty colonies with IDs + +## New features + +- parallelised all the major functions (so they run on + simParamBee\$nThreads cores) with PSOCK system. Since the parallelisation setup within functions + takes additional time, we recommend using a single threads for a small number of colonies + +## Bug fixes + + # SIMplyBee version 0.4.1 - 2024-09-19 @@ -57,8 +84,6 @@ which caused an error. We now read in the locations from a csv file. - Bug fix - get\*Haplo() functions were returning diploid drones when input was a Pop-class -- - # SIMplyBee version 0.3.0 - 2022-12-05 First public/CRAN version of the package diff --git a/R/Class-Colony.R b/R/Class-Colony.R index 7bdd2595..d5fcff2b 100644 --- a/R/Class-Colony.R +++ b/R/Class-Colony.R @@ -94,7 +94,8 @@ setClassUnion("ColonyOrNULL", c("Colony", "NULL")) setValidity(Class = "Colony", method = function(object) { errors <- character() - if ((ifelse(test = !is.null(slot(object, name = "queen")), yes = nInd(slot(object, name = "queen")), no = 0)) > 1) { #Don't use nQueen because of the SP problem + test <- !is.null(slot(object, name = "queen")) + if ((ifelse(test, yes = nInd(slot(object, name = "queen")), no = 0)) > 1) { #Don't use nQueen because of the SP problem errors <- c(errors, "There can be only one queen per colony!") } if (length(errors) == 0) { diff --git a/R/Class-SimParamBee.R b/R/Class-SimParamBee.R index 5ca066c3..762541b5 100644 --- a/R/Class-SimParamBee.R +++ b/R/Class-SimParamBee.R @@ -425,13 +425,112 @@ SimParamBee <- R6Class( invisible(self) }, + #' @description For internal use only. + #' + #' @param nNewInd Number of newly created individuals + #' @param id the name of each individual + #' @param mother vector of mother iids + #' @param father vector of father iids + #' @param isDH indicator for DH lines + addToBeePed = function(nNewInd,id,mother,father,isDH) { + stopifnot(nNewInd>0) + if(length(isDH)==1) isDH = rep(isDH,nNewInd) + mother = as.integer(mother) + father = as.integer(father) + isDH = as.integer(isDH) + stopifnot(length(mother)==nNewInd, + length(father)==nNewInd, + length(isDH)==nNewInd) + tmp = cbind(mother,father,isDH) + rownames(tmp) = id + private$.pedigree = rbind(private$.pedigree,tmp) + private$.lastId = private$.lastId + as.integer(nNewInd) + invisible(self) + }, + + + #' @description For internal use only. + #' + #' @param nNewInd Number of newly created individuals + #' @param id the name of each individual + #' @param mother vector of mother iids + #' @param father vector of father iids + #' @param isDH indicator for DH lines + #' @param hist new recombination history + #' @param ploidy ploidy level + addToBeeRec = function(nNewInd,id,mother,father,isDH, + hist,ploidy){ + stopifnot(nNewInd>0) + if(length(isDH)==1) isDH = rep(isDH,nNewInd) + mother = as.integer(mother) + father = as.integer(father) + isDH = as.integer(isDH) + stopifnot(length(mother)==nNewInd, + length(father)==nNewInd, + length(isDH)==nNewInd) + tmp = cbind(mother,father,isDH) + rownames(tmp) = id + if(is.null(hist)){ + newRecHist = vector("list",nNewInd) + tmpLastHaplo = private$.lastHaplo + if(all(isDH==1L)){ + for(i in seq_len(nNewInd)){ + tmpLastHaplo = tmpLastHaplo + 1L + newRecHist[[i]] = rep(tmpLastHaplo, ploidy) + } + }else{ + for(i in seq_len(nNewInd)){ + newRecHist[[i]] = (tmpLastHaplo+1L):(tmpLastHaplo+ploidy) + tmpLastHaplo = tmpLastHaplo + ploidy + } + } + private$.hasHap = c(private$.hasHap, rep(FALSE, nNewInd)) + private$.isFounder = c(private$.isFounder, rep(TRUE, nNewInd)) + #names(newRecHist) = id + private$.recHist = c(private$.recHist, newRecHist) + private$.lastHaplo = tmpLastHaplo + }else{ + # Add hist to recombination history + private$.hasHap = c(private$.hasHap, rep(FALSE, nNewInd)) + private$.isFounder = c(private$.isFounder, rep(FALSE, nNewInd)) + #names(hist) = id + private$.recHist = c(private$.recHist, hist) + } + private$.pedigree = rbind(private$.pedigree, tmp) + private$.lastId = private$.lastId + as.integer(nNewInd) + + invisible(self) + }, + + #' @description A function to update the caste + #' For internal use only. + #' + #' @param caste vector, named vector of castes to be added + updateCaste = function(caste) { + private$.caste = c(private$.caste, caste) + invisible(self) + }, + + #' @description A function to update the last + #' ID everytime we create an individual + #' For internal use in SIMplyBee only. + #' + #' @param lastId integer, last colony ID assigned + #' @param n integer, how many individuals to add + updateLastBeeId = function(n = 1L) { + private$.lastId = private$.lastId + as.integer(n) + invisible(self) + }, + #' @description A function to update the colony last #' ID everytime we create a Colony-class with createColony. #' For internal use only. #' #' @param lastColonyId integer, last colony ID assigned - updateLastColonyId = function() { - private$.lastColonyId = private$.lastColonyId + 1L + #' @param n integer, how many colonies to add + updateLastColonyId = function(n = 1) { + n = as.integer(n) + private$.lastColonyId = private$.lastColonyId + n invisible(self) } ), @@ -459,7 +558,7 @@ SimParamBee <- R6Class( #' created caste = function(value) { if (missing(value)) { - x = private$.caste + x = private$.caste } else { stop("`$caste` is read only", call. = FALSE) } @@ -469,7 +568,7 @@ SimParamBee <- R6Class( #' created with \code{\link[SIMplyBee]{createColony}} lastColonyId = function(value) { if (missing(value)) { - private$.lastColonyId + private$.lastColonyId } else { stop("`$lastColonyId` is read only", call. = FALSE) } @@ -571,16 +670,15 @@ isSimParamBee <- function(x) { # nFunctions ---- -#' @rdname nWorkersFun -#' @title Sample a number of workers +#' @rdname nCasteFun +#' @title Sample a number of caste members (workers, drones, virgin queens) #' -#' @description Sample a number of workers - used when \code{nInd = NULL} -#' (see \code{\link[SIMplyBee]{SimParamBee}$nWorkers}). +#' @description Sample a number of caste member - used when \code{nInd = NULL} #' #' This is just an example. You can provide your own functions that satisfy #' your needs! #' -#' @param colony \code{\link[SIMplyBee]{Colony-class}} +#' @param x \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} #' @param n integer, number of samples #' @param average numeric, average number of workers #' @param lowerLimit numeric, returned numbers will be above this value @@ -596,35 +694,35 @@ isSimParamBee <- function(x) { #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters #' @param ... other arguments of \code{\link[SIMplyBee]{mapCasteToColonyPheno}} #' -#' @details \code{nWorkersPoisson} samples from a Poisson distribution with a -#' given average, which can return a value 0. \code{nDronesTruncPoisson} +#' @details \code{nCastePoisson} samples from a Poisson distribution with a +#' given average, which can return a value 0. \code{nCasteTruncPoisson} #' samples from a zero truncated Poisson distribution. #' -#' \code{nWorkersColonyPhenotype} returns a number (above \code{lowerLimit}) +#' \code{nCasteColonyPhenotype} returns a number (above \code{lowerLimit}) #' as a function of colony phenotype, say queen's fecundity. Colony phenotype #' is provided by \code{\link[SIMplyBee]{mapCasteToColonyPheno}}. You need to set up #' traits influencing the colony phenotype and their parameters (mean and #' variances) via \code{\link[SIMplyBee]{SimParamBee}} (see examples). #' -#' @seealso \code{\link[SIMplyBee]{SimParamBee}} field \code{nWorkers} and +#' @seealso \code{\link[SIMplyBee]{SimParamBee}} field \code{nWorkers}, \code{nDrones} and #' \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")} #' -#' @return numeric, number of workers +#' @return numeric, number of caste members #' #' @examples -#' nWorkersPoisson() -#' nWorkersPoisson() -#' n <- nWorkersPoisson(n = 1000) +#' nCastePoisson() +#' nCastePoisson() +#' n <- nCastePoisson(n = 1000) #' hist(n, breaks = seq(from = min(n), to = max(n)), xlim = c(0, 200)) #' table(n) #' -#' nWorkersTruncPoisson() -#' nWorkersTruncPoisson() -#' n <- nWorkersTruncPoisson(n = 1000) +#' nCasteTruncPoisson() +#' nCasteTruncPoisson() +#' n <- nCasteTruncPoisson(n = 1000) #' hist(n, breaks = seq(from = min(n), to = max(n)), xlim = c(0, 200)) #' table(n) #' -#' # Example for nWorkersColonyPhenotype() +#' # Example for nCasteColonyPhenotype() #' founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100) #' SP <- SimParamBee$new(founderGenomes) #' \dontshow{SP$nThreads = 1L} @@ -641,327 +739,129 @@ isSimParamBee <- function(x) { #' colony2 <- cross(colony2, drones = droneGroups[[2]]) #' colony1@queen@pheno #' colony2@queen@pheno -#' createWorkers(colony1, nInd = nWorkersColonyPhenotype) -#' createWorkers(colony2, nInd = nWorkersColonyPhenotype) +#' createWorkers(colony1, nInd = nCasteColonyPhenotype) +#' createWorkers(colony2, nInd = nCasteColonyPhenotype) #' @export -nWorkersPoisson <- function(colony, n = 1, average = 100) { +nCastePoisson <- function(x, n = 1, average = 100) { + # We keep the x because for nCasteColonyPhenotype we need colony/multicolony access + # These are used inside other functions when these n functions are called + if (isColony(x)) { + n <- 1 + } else if (isMultiColony(x)) { + n <- nColonies(x) + } return(rpois(n = n, lambda = average)) } -#' @describeIn nWorkersFun Sample a non-zero number of workers +#' @describeIn nCastePoisson #' @export -nWorkersTruncPoisson <- function(colony, n = 1, average = 100, lowerLimit = 0) { - return(extraDistr::rtpois(n = n, lambda = average, a = lowerLimit)) +nVirginQueensPoisson <- function(x, n = 1, average = 10) { + nCastePoisson(x = x, n = n, average = average) +} +#' @describeIn nCastePoisson +#' @export +nFathersPoisson <- function(x, n = 1, average = 15) { + nCastePoisson(x = x, n = n, average = average) +} +#' @describeIn nCastePoisson +#' @export +nWorkersPoisson <- function(x, n = 1, average = 100) { + nCastePoisson(x = x, n = n, average = average) +} +#' @describeIn nCastePoisson +#' @export +nDronesPoisson <- function(x, n = 1, average = 100) { + nCastePoisson(x = x, n = n, average = average) } -#' @describeIn nWorkersFun Sample a non-zero number of workers based on -#' colony phenotype, say queen's fecundity +#' @describeIn nCasteFun Sample a non-zero number of caste individuals #' @export -nWorkersColonyPhenotype <- function(colony, queenTrait = 1, workersTrait = NULL, - checkProduction = FALSE, lowerLimit = 0, - simParamBee = NULL, - ...) { - if (is.null(simParamBee)) { - simParamBee <- get(x = "SP", envir = .GlobalEnv) +nCasteTruncPoisson <- function(x, n = 1, average = 100, lowerLimit = 0) { + if (isColony(x)) { + n <- 1 + } else if (isMultiColony(x)) { + n <- nColonies(x) } - ret <- round(mapCasteToColonyPheno( - colony = colony, - queenTrait = queenTrait, - workersTrait = workersTrait, - checkProduction = checkProduction, - simParamBee = simParamBee, - ... - )) - if (ret < (lowerLimit + 1)) { - ret <- lowerLimit + 1 - } - return(ret) + return(extraDistr::rtpois(n = n, lambda = average, a = lowerLimit)) } -#' @rdname nDronesFun -#' @title Sample a number of drones -#' -#' @description Sample a number of drones - used when \code{nDrones = NULL} -#' (see \code{\link[SIMplyBee]{SimParamBee}$nDrones}). -#' -#' This is just an example. You can provide your own functions that satisfy -#' your needs! -#' -#' @param x \code{\link[AlphaSimR]{Pop-class}} or \code{\link[SIMplyBee]{Colony-class}} -#' @param n integer, number of samples -#' @param average numeric, average number of drones -#' @param lowerLimit numeric, returned numbers will be above this value -#' @param queenTrait numeric (column position) or character (column name), trait -#' that represents queen's effect on the colony phenotype (defined in -#' \code{\link[SIMplyBee]{SimParamBee}} - see examples); if \code{0} then this effect is 0 -#' @param workersTrait numeric (column position) or character (column name), trait -#' that represents workers's effect on the colony phenotype (defined in -#' \code{\link[SIMplyBee]{SimParamBee}} - see examples); if \code{0} then this effect is 0 -#' @param checkProduction logical, does the phenotype depend on the production -#' status of colony; if yes and production is not \code{TRUE}, the result is -#' above \code{lowerLimit} -#' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters -#' @param ... other arguments of \code{\link[SIMplyBee]{mapCasteToColonyPheno}} -#' -#' @details \code{nDronesPoisson} samples from a Poisson distribution with a -#' given average, which can return a value 0. -#' -#' \code{nDronesTruncPoisson} samples from a zero truncated Poisson -#' distribution. -#' -#' \code{nDronesColonyPhenotype} returns a number (above \code{lowerLimit}) as -#' a function of colony phenotype, say queen's fecundity. Colony phenotype is -#' provided by \code{\link[SIMplyBee]{mapCasteToColonyPheno}}. You need to set up -#' traits influencing the colony phenotype and their parameters (mean and -#' variances) via \code{\link[SIMplyBee]{SimParamBee}} (see examples). -#' -#' When \code{x} is \code{\link[AlphaSimR]{Pop-class}}, only \code{workersTrait} is not -#' used, that is, only \code{queenTrait} is used. -#' -#' @seealso \code{\link[SIMplyBee]{SimParamBee}} field \code{nDrones} and -#' \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")} -#' -#' @return numeric, number of drones -#' -#' @examples -#' nDronesPoisson() -#' nDronesPoisson() -#' n <- nDronesPoisson(n = 1000) -#' hist(n, breaks = seq(from = min(n), to = max(n)), xlim = c(0, 200)) -#' table(n) -#' -#' nDronesTruncPoisson() -#' nDronesTruncPoisson() -#' n <- nDronesTruncPoisson(n = 1000) -#' hist(n, breaks = seq(from = min(n), to = max(n)), xlim = c(0, 200)) -#' table(n) -#' -#' # Example for nDronesColonyPhenotype() -#' founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100) -#' SP <- SimParamBee$new(founderGenomes) -#' \dontshow{SP$nThreads = 1L} -#' average <- 100 -#' h2 <- 0.1 -#' SP$addTraitA(nQtlPerChr = 100, mean = average, var = average * h2) -#' SP$setVarE(varE = average * (1 - h2)) -#' basePop <- createVirginQueens(founderGenomes) -#' drones <- createDrones(x = basePop[1], nInd = 50) -#' droneGroups <- pullDroneGroupsFromDCA(drones, n = 2, nDrones = 15) -#' colony1 <- createColony(x = basePop[2]) -#' colony2 <- createColony(x = basePop[3]) -#' colony1 <- cross(colony1, drones = droneGroups[[1]]) -#' colony2 <- cross(colony2, drones = droneGroups[[2]]) -#' colony1@queen@pheno -#' colony2@queen@pheno -#' createDrones(colony1, nInd = nDronesColonyPhenotype) -#' createDrones(colony2, nInd = nDronesColonyPhenotype) +#' @describeIn nCasteTruncPoisson #' @export -nDronesPoisson <- function(x, n = 1, average = 100) { - return(rpois(n = n, lambda = average)) +nVirginQueensTruncPoisson <- function(x, n = 1, average = 10, lowerLimit = 0) { + nCasteTruncPoisson(x = x, n = n, average = average, lowerLimit = lowerLimit) } - -#' @describeIn nDronesFun Sample a non-zero number of drones +#' @describeIn nCasteTruncPoisson +#' @export +nFathersTruncPoisson <- function(x, n = 1, average = 15, lowerLimit = 0) { + nCasteTruncPoisson(x = x, n = n, average = average, lowerLimit = lowerLimit) +} +#' @describeIn nCasteTruncPoisson +#' @export +nWorkersTruncPoisson <- function(x, n = 1, average = 100, lowerLimit = 0) { + nCasteTruncPoisson(x = x, n = n, average = average, lowerLimit = lowerLimit) +} +#' @describeIn nCasteTruncPoisson #' @export nDronesTruncPoisson <- function(x, n = 1, average = 100, lowerLimit = 0) { - return(extraDistr::rtpois(n = n, lambda = average, a = lowerLimit)) + nCasteTruncPoisson(x = x, n = n, average = average, lowerLimit = lowerLimit) } -#' @describeIn nDronesFun Sample a non-zero number of drones based on +#' @describeIn nCasteFun Sample a non-zero number of caste individuals based on #' colony phenotype, say queen's fecundity #' @export -nDronesColonyPhenotype <- function(x, queenTrait = 1, workersTrait = NULL, - checkProduction = FALSE, lowerLimit = 0, - simParamBee = NULL, - ...) { +nCasteColonyPhenotype <- function(x, queenTrait = 1, workersTrait = NULL, + checkProduction = FALSE, lowerLimit = 0, + simParamBee = NULL, + ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } - # This one is special because we cater drone production from base population - # virgin queens and colonies if (isPop(x)) { ret <- round(x@pheno[, queenTrait]) } else { - ret <- round(mapCasteToColonyPheno( - colony = x, + ret <- mapCasteToColonyPheno( + x = x, queenTrait = queenTrait, workersTrait = workersTrait, checkProduction = checkProduction, simParamBee = simParamBee, ... - )) - } - if (ret < (lowerLimit + 1)) { - ret <- lowerLimit + 1 + ) + ret <- sapply(ret, FUN = function(x) round(x)) + test <- ret < (lowerLimit + 1) + if (any(test)) { + ret[test] <- lowerLimit + 1 + } } return(ret) } -#' @rdname nVirginQueensFun -#' @title Sample a number of virgin queens -#' -#' @description Sample a number of virgin queens - used when -#' \code{nFathers = NULL} (see \code{\link[SIMplyBee]{SimParamBee}$nVirginQueens}). -#' -#' This is just an example. You can provide your own functions that satisfy -#' your needs! -#' -#' @param colony \code{\link[SIMplyBee]{Colony-class}} -#' @param n integer, number of samples -#' @param average numeric, average number of virgin queens -#' @param lowerLimit numeric, returned numbers will be above this value -#' @param queenTrait numeric (column position) or character (column name), trait -#' that represents queen's effect on the colony phenotype (defined in -#' \code{\link[SIMplyBee]{SimParamBee}} - see examples); if \code{NULL} then this effect is 0 -#' @param workersTrait numeric (column position) or character (column name), trait -#' that represents workers's effect on the colony phenotype (defined in -#' \code{\link[SIMplyBee]{SimParamBee}} - see examples); if \code{NULL} then this effect is 0 -#' @param checkProduction logical, does the phenotype depend on the production -#' status of colony; if yes and production is not \code{TRUE}, the result is -#' above \code{lowerLimit} -#' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters -#' @param ... other arguments of \code{\link[SIMplyBee]{mapCasteToColonyPheno}} -#' -#' @details \code{nVirginQueensPoisson} samples from a Poisson distribution, -#' which can return a value 0 (that would mean a colony will fail to raise a -#' single virgin queen after the queen swarms or dies). -#' -#' \code{nVirginQueensTruncPoisson} samples from a truncated Poisson -#' distribution (truncated at zero) to avoid failure. -#' -#' \code{nVirginQueensColonyPhenotype} returns a number (above -#' \code{lowerLimit}) as a function of colony phenotype, say swarming -#' tendency. Colony phenotype is provided by -#' \code{\link[SIMplyBee]{mapCasteToColonyPheno}}. You need to set up traits -#' influencing the colony phenotype and their parameters (mean and variances) -#' via \code{\link[SIMplyBee]{SimParamBee}} (see examples). -#' -#' @seealso \code{\link[SIMplyBee]{SimParamBee}} field \code{nVirginQueens} and -#' \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")} -#' -#' @return numeric, number of virgin queens -#' -#' @examples -#' nVirginQueensPoisson() -#' nVirginQueensPoisson() -#' n <- nVirginQueensPoisson(n = 1000) -#' hist(n, breaks = seq(from = min(n), to = max(n)), xlim = c(0, 30)) -#' table(n) -#' -#' nVirginQueensTruncPoisson() -#' nVirginQueensTruncPoisson() -#' n <- nVirginQueensTruncPoisson(n = 1000) -#' hist(n, breaks = seq(from = min(n), to = max(n)), xlim = c(0, 30)) -#' table(n) -#' -#' # Example for nVirginQueensColonyPhenotype() -#' founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100) -#' SP <- SimParamBee$new(founderGenomes) -#' \dontshow{SP$nThreads = 1L} -#' # Setting trait scale such that mean is 10 split into queen and workers effects -#' meanP <- c(5, 5 / SP$nWorkers) -#' # setup variances such that the total phenotype variance will match the mean -#' varA <- c(3 / 2, 3 / 2 / SP$nWorkers) -#' corA <- matrix(data = c( -#' 1.0, -0.5, -#' -0.5, 1.0 -#' ), nrow = 2, byrow = TRUE) -#' varE <- c(7 / 2, 7 / 2 / SP$nWorkers) -#' varA / (varA + varE) -#' varP <- varA + varE -#' varP[1] + varP[2] * SP$nWorkers -#' SP$addTraitA(nQtlPerChr = 100, mean = meanP, var = varA, corA = corA) -#' SP$setVarE(varE = varE) -#' basePop <- createVirginQueens(founderGenomes) -#' drones <- createDrones(x = basePop[1], nInd = 50) -#' droneGroups <- pullDroneGroupsFromDCA(drones, n = 2, nDrones = 15) -#' colony1 <- createColony(x = basePop[2]) -#' colony2 <- createColony(x = basePop[3]) -#' colony1 <- cross(colony1, drones = droneGroups[[1]]) -#' colony2 <- cross(colony2, drones = droneGroups[[2]]) -#' colony1 <- buildUp(colony1) -#' colony2 <- buildUp(colony2) -#' nVirginQueensColonyPhenotype(colony1) -#' nVirginQueensColonyPhenotype(colony2) -#' @export -nVirginQueensPoisson <- function(colony, n = 1, average = 10) { - return(rpois(n = n, lambda = average)) -} - -#' @describeIn nVirginQueensFun Sample a non-zero number of virgin queens +#' @describeIn nCasteColonyPhenotype #' @export -nVirginQueensTruncPoisson <- function(colony, n = 1, average = 10, lowerLimit = 0) { - return(extraDistr::rtpois(n = n, lambda = average, a = lowerLimit)) -} - -#' @describeIn nVirginQueensFun Sample a non-zero number of virgin queens -#' based on colony's phenotype, say, swarming tendency -#' @export -nVirginQueensColonyPhenotype <- function(colony, queenTrait = 1, - workersTrait = 2, - checkProduction = FALSE, - lowerLimit = 0, - simParamBee = NULL, - ...) { - if (is.null(simParamBee)) { - simParamBee <- get(x = "SP", envir = .GlobalEnv) - } - ret <- round(mapCasteToColonyPheno( - colony = colony, - queenTrait = queenTrait, - workersTrait = workersTrait, - simParamBee = simParamBee, - ... - )) - if (ret < (lowerLimit + 1)) { - ret <- lowerLimit + 1 - } - return(ret) +nVirginQueensColonyPhenotype <- function(x, queenTrait = 1, workersTrait = NULL, + checkProduction = FALSE, lowerLimit = 0, + simParamBee = NULL, + ...) { + nCasteColonyPhenotype(x = x, queenTrait = queenTrait, workersTrait = workersTrait, + checkProduction = checkProduction, lowerLimit = lowerLimit, + simParamBee = simParamBee, + ...) } - -#' @rdname nFathersFun -#' @title Sample a number of fathers -#' -#' @description Sample a number of fathers - use when \code{nFathers = NULL} -#' (see \code{\link[SIMplyBee]{SimParamBee}$nFathers}). -#' -#' This is just an example. You can provide your own functions that satisfy -#' your needs! -#' -#' @param n integer, number of samples -#' @param average numeric, average number of fathers -#' @param lowerLimit numeric, returned numbers will be above this value -#' -#' @details \code{nFathersPoisson} samples from a Poisson distribution, which -#' can return a value 0 (that would mean a failed queen mating). -#' -#' \code{nFathersTruncPoisson} samples from a truncated Poisson distribution -#' (truncated at zero) to avoid failed matings. -#' -#' @seealso \code{\link[SIMplyBee]{SimParamBee}} field \code{nFathers} -#' -#' @return numeric, number of fathers -#' -#' @examples -#' nFathersPoisson() -#' nFathersPoisson() -#' n <- nFathersPoisson(n = 1000) -#' hist(n, breaks = seq(from = min(n), to = max(n)), xlim = c(0, 40)) -#' table(n) -#' -#' nFathersTruncPoisson() -#' nFathersTruncPoisson() -#' n <- nFathersTruncPoisson(n = 1000) -#' hist(n, breaks = seq(from = min(n), to = max(n)), xlim = c(0, 40)) -#' table(n) +#' @describeIn nCasteColonyPhenotype #' @export -nFathersPoisson <- function(n = 1, average = 15) { - return(rpois(n = n, lambda = average)) +nWorkersColonyPhenotype <- function(x, n = 1, average = 100, lowerLimit = 0) { + nCasteColonyPhenotype(x = x, queenTrait = queenTrait, workersTrait = workersTrait, + checkProduction = checkProduction, lowerLimit = lowerLimit, + simParamBee = simParamBee, + ...) } - -#' @describeIn nFathersFun Sample a non-zero number of fathers +#' @describeIn nCasteColonyPhenotype #' @export -nFathersTruncPoisson <- function(n = 1, average = 15, lowerLimit = 0) { - return(extraDistr::rtpois(n = n, lambda = average, a = lowerLimit)) +nDronesColonyPhenotype <- function(x, n = 1, average = 100, lowerLimit = 0) { + nCasteColonyPhenotype(x = x, queenTrait = queenTrait, workersTrait = workersTrait, + checkProduction = checkProduction, lowerLimit = lowerLimit, + simParamBee = simParamBee, + ...) } # pFunctions ---- @@ -975,7 +875,7 @@ nFathersTruncPoisson <- function(n = 1, average = 15, lowerLimit = 0) { #' This is just an example. You can provide your own functions that satisfy #' your needs! #' -#' @param colony \code{\link[SIMplyBee]{Colony-class}} +#' @param x \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} #' @param n integer, number of samples #' @param min numeric, lower limit for \code{swarmPUnif} #' @param max numeric, upper limit for \code{swarmPUnif} @@ -998,7 +898,12 @@ nFathersTruncPoisson <- function(n = 1, average = 15, lowerLimit = 0) { #' p <- swarmPUnif(n = 1000) #' hist(p, breaks = seq(from = 0, to = 1, by = 0.01), xlim = c(0, 1)) #' @export -swarmPUnif <- function(colony, n = 1, min = 0.4, max = 0.6) { +swarmPUnif <- function(x, n = 1, min = 0.4, max = 0.6) { + if (isColony(x)) { + n <- 1 + } else if (isMultiColony(x)) { + n <- nColonies(x) + } return(runif(n = n, min = min, max = max)) } @@ -1013,7 +918,7 @@ swarmPUnif <- function(colony, n = 1, min = 0.4, max = 0.6) { #' This is just an example. You can provide your own functions that satisfy #' your needs! #' -#' @param colony \code{\link[SIMplyBee]{Colony-class}} +#' @param x \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} #' @param n integer, number of samples #' @param min numeric, lower limit for \code{splitPUnif} #' @param max numeric, upper limit for \code{splitPUnif} @@ -1079,14 +984,24 @@ swarmPUnif <- function(colony, n = 1, min = 0.4, max = 0.6) { #' plot(pKeep ~ nWorkers, ylim = c(0, 1)) #' abline(v = nWorkersFull) #' @export -splitPUnif <- function(colony, n = 1, min = 0.2, max = 0.4) { +splitPUnif <- function(x, n = 1, min = 0.2, max = 0.4) { + if (isColony(x)) { + n <- 1 + } else if (isMultiColony(x)) { + n <- nColonies(x) + } return(runif(n = n, min = min, max = max)) } #' @describeIn splitPFun Sample the split proportion - the proportion of #' removed workers in a managed split based on the colony strength #' @export -splitPColonyStrength <- function(colony, n = 1, nWorkersFull = 100, scale = 1) { +splitPColonyStrength <- function(x, n = 1, nWorkersFull = 100, scale = 1) { + if (isColony(x)) { + n <- 1 + } else if (isMultiColony(x)) { + n <- nColonies(x) + } nW <- nWorkers(colony) pKeep <- rbeta( n = n, @@ -1107,7 +1022,7 @@ splitPColonyStrength <- function(colony, n = 1, nWorkersFull = 100, scale = 1) { #' This is just an example. You can provide your own functions that satisfy #' your needs! #' -#' @param colony \code{\link[SIMplyBee]{Colony-class}} +#' @param x \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} #' @param n integer, number of samples #' @param min numeric, lower limit for \code{downsizePUnif} #' @param max numeric, upper limit for \code{downsizePUnif} @@ -1122,7 +1037,12 @@ splitPColonyStrength <- function(colony, n = 1, nWorkersFull = 100, scale = 1) { #' p <- downsizePUnif(n = 1000) #' hist(p, breaks = seq(from = 0, to = 1, by = 0.01), xlim = c(0, 1)) #' @export -downsizePUnif <- function(colony, n = 1, min = 0.8, max = 0.9) { +downsizePUnif <- function(x, n = 1, min = 0.8, max = 0.9) { + if (isColony(x)) { + n <- 1 + } else if (isMultiColony(x)) { + n <- nColonies(x) + } return(runif(n = n, min = min, max = max)) } @@ -1146,21 +1066,24 @@ downsizePUnif <- function(colony, n = 1, min = 0.8, max = 0.9) { #' Note though that you can achieve this impact also via multiple correlated #' traits, such as a queen and a workers trait. #' -#' @param colony \code{\link[SIMplyBee]{Colony-class}} +#' @param x \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} #' @param value character, one of \code{pheno} or \code{gv} #' @param queenTrait numeric (column position) or character (column name), #' trait(s) that represents queen's contribution to colony value(s); if -#' \code{NULL} then this contribution is 0; you can pass more than one trait +#' \code{NULL} or there is no queen present, then this contribution is 0; +#' you can pass more than one trait #' here, but make sure that \code{combineFUN} works with these trait dimensions -#' @param queenFUN function, function that will be applied to queen's value +#' @param queenFUN function, function that will be applied to queen's value. #' @param workersTrait numeric (column position) or character (column name), #' trait(s) that represents workers' contribution to colony value(s); if -#' \code{NULL} then this contribution is 0; you can pass more than one trait +#' \code{NULL} or there are no workers present, then this contribution is 0; +#' you can pass more than one trait #' here, but make sure that \code{combineFUN} works with these trait dimensions #' @param workersFUN function, function that will be applied to workers values #' @param dronesTrait numeric (column position) or character (column name), #' trait(s) that represents drones' contribution to colony value(s); if -#' \code{NULL} then this contribution is 0; you can pass more than one trait +#' \code{NULL} or there are no drones present then this contribution is 0; +#' you can pass more than one trait #' here, but make sure that \code{combineFUN} works with these trait dimensions #' @param dronesFUN function, function that will be applied to drone values #' @param traitName, the name of the colony trait(s), say, honeyYield; you can pass @@ -1191,7 +1114,8 @@ downsizePUnif <- function(colony, n = 1, min = 0.8, max = 0.9) { #' \code{\link[SIMplyBee]{calcColonyValue}}. It only works on a single colony - use #' \code{\link[SIMplyBee]{calcColonyValue}} to get Colony or MultiColony values. #' -#' @return numeric matrix with one value or a row of values +#' @return numeric matrix with one value or a row of values when input is \code{\link[SIMplyBee]{Colony-class}} +#' or list of numeric matrices when input is \code{\link[SIMplyBee]{MultiColony-class}} #' #' @examples #' founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 100) @@ -1245,7 +1169,7 @@ downsizePUnif <- function(colony, n = 1, min = 0.8, max = 0.9) { # https://github.com/HighlanderLab/SIMplyBee/issues/353 # TODO: Develop theory for colony genetic values under non-linearity/non-additivity #403 # https://github.com/HighlanderLab/SIMplyBee/issues/403 -mapCasteToColonyValue <- function(colony, +mapCasteToColonyValue <- function(x, value = "pheno", queenTrait = 1, queenFUN = function(x) x, workersTrait = 2, workersFUN = colSums, @@ -1267,95 +1191,129 @@ mapCasteToColonyValue <- function(colony, if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } - if (is.null(queenTrait)) { - queenEff <- 0 - } else { - if (value %in% c("pheno", "gv")) { - tmp <- valueFUN(colony@queen)[, queenTrait, drop = FALSE] - } else { # bv, dd, and aa: leaving this in for future use! - tmp <- valueFUN(colony@queen, simParam = simParamBee)[, queenTrait, drop = FALSE] - } - queenEff <- queenFUN(tmp) - } - if (is.null(workersTrait)) { - workersEff <- 0 - } else { - if (value %in% c("pheno", "gv")) { - tmp <- valueFUN(colony@workers)[, workersTrait, drop = FALSE] - } else { # bv, dd, and aa - tmp <- valueFUN(colony@workers, simParam = simParamBee)[, workersTrait, drop = FALSE] + + if (isColony(x)) { + if (is.null(queenTrait)) { + queenEff <- 0 + } else { + if (isQueenPresent(x)) { + if (value %in% c("pheno", "gv")) { + tmp <- valueFUN(x@queen)[, queenTrait, drop = FALSE] + } else { # bv, dd, and aa: leaving this in for future use! + tmp <- valueFUN(x@queen, simParam = simParamBee)[, queenTrait, drop = FALSE] + } + queenEff <- queenFUN(tmp) + } else { + queenEff <- 0 + } } - workersEff <- workersFUN(tmp) - } - if (is.null(dronesTrait)) { - dronesEff <- 0 - } else { - if (value %in% c("pheno", "gv")) { - tmp <- valueFUN(colony@drones)[, dronesTrait, drop = FALSE] - } else { # bv, dd, and aa - tmp <- valueFUN(colony@drones, simParam = simParamBee)[, dronesTrait, drop = FALSE] + if (is.null(workersTrait)) { + workersEff <- 0 + } else { + if (nWorkers(x) != 0) { + if (value %in% c("pheno", "gv")) { + tmp <- valueFUN(x@workers)[, workersTrait, drop = FALSE] + } else { # bv, dd, and aa + tmp <- valueFUN(x@workers, simParam = simParamBee)[, workersTrait, drop = FALSE] + } + workersEff <- workersFUN(tmp) + } else { + workersEff <- 0 + } } - dronesEff <- dronesFUN(tmp) - } - colonyValue <- combineFUN(q = queenEff, w = workersEff, d = dronesEff) - nColTrt <- length(colonyValue) - colnames(colonyValue) <- traitName - if (any(checkProduction) && !isProductive(colony)) { - if (length(checkProduction) == 1 && nColTrt != 1) { - checkProduction <- rep(checkProduction, times = nColTrt) + if (is.null(dronesTrait)) { + dronesEff <- 0 + } else { + if (nDrones(x) != 0) { + if (value %in% c("pheno", "gv")) { + tmp <- valueFUN(x@drones)[, dronesTrait, drop = FALSE] + } else { # bv, dd, and aa + tmp <- valueFUN(x@drones, simParam = simParamBee)[, dronesTrait, drop = FALSE] + } + dronesEff <- dronesFUN(tmp) + } else { + dronesEff <- 0 + } } - if (length(notProductiveValue) == 1 && nColTrt != 1) { - notProductiveValue <- rep(notProductiveValue, times = nColTrt) + colonyValue <- combineFUN(q = queenEff, w = workersEff, d = dronesEff) + nColTrt <- length(colonyValue) + colnames(colonyValue) <- traitName + if (any(checkProduction) && !isProductive(x)) { + if (length(checkProduction) == 1 && nColTrt != 1) { + checkProduction <- rep(checkProduction, times = nColTrt) + } + if (length(notProductiveValue) == 1 && nColTrt != 1) { + notProductiveValue <- rep(notProductiveValue, times = nColTrt) + } + if (length(checkProduction) != nColTrt) { + stop("Dimension of checkProduction does not match the number of traits from combineFUN()!") + } + if (length(checkProduction) != length(notProductiveValue)) { + stop("Dimensions of checkProduction and notProductiveValue must match!") + } + colonyValue[checkProduction] <- notProductiveValue[checkProduction] } - if (length(checkProduction) != nColTrt) { - stop("Dimension of checkProduction does not match the number of traits from combineFUN()!") + } else if (isMultiColony(x)) { + nCol <- nColonies(x) + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") } - if (length(checkProduction) != length(notProductiveValue)) { - stop("Dimensions of checkProduction and notProductiveValue must match!") + colonyValue <- vector(mode = "list", length = nCol) + names(colonyValue) <- getId(x) + for (colony in 1:nCol) { + colonyValue[[colony]] <- mapCasteToColonyValue(x[[colony]], + value = value, + queenTrait = queenTrait, queenFUN = queenFUN, + workersTrait = workersTrait, workersFUN = workersFUN, + dronesTrait = dronesTrait, dronesFUN = dronesFUN, + traitName = traitName, + combineFUN = combineFUN, + checkProduction = checkProduction, + notProductiveValue = notProductiveValue, + simParamBee = simParamBee) } - colonyValue[checkProduction] <- notProductiveValue[checkProduction] } return(colonyValue) } #' @describeIn mapCasteToColonyValue Map caste member (individual) phenotype values to a colony phenotype value #' @export -mapCasteToColonyPheno <- function(colony, simParamBee = NULL, ...) { +mapCasteToColonyPheno <- function(x, simParamBee = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } - mapCasteToColonyValue(colony, value = "pheno", simParamBee = simParamBee, ...) + mapCasteToColonyValue(x, value = "pheno", simParamBee = simParamBee, ...) } #' @describeIn mapCasteToColonyValue Map caste member (individual) genetic values to a colony genetic value #' @export -mapCasteToColonyGv <- function(colony, simParamBee = NULL, ...) { +mapCasteToColonyGv <- function(x, simParamBee = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } - mapCasteToColonyValue(colony, value = "gv", checkProduction = FALSE, simParamBee = simParamBee, ...) + mapCasteToColonyValue(x, value = "gv", checkProduction = FALSE, simParamBee = simParamBee, ...) } #' @describeIn mapCasteToColonyValue Map caste member (individual) breeding values to a colony breeding value -mapCasteToColonyBv <- function(colony, simParamBee = NULL, ...) { +mapCasteToColonyBv <- function(x, simParamBee = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } - mapCasteToColonyValue(colony, value = "bv", checkProduction = FALSE, simParamBee = simParamBee, ...) + mapCasteToColonyValue(x, value = "bv", checkProduction = FALSE, simParamBee = simParamBee, ...) } #' @describeIn mapCasteToColonyValue Map caste member (individual) dominance values to a colony dominance value -mapCasteToColonyDd <- function(colony, simParamBee = NULL, ...) { +mapCasteToColonyDd <- function(x, simParamBee = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } - mapCasteToColonyValue(colony, value = "dd", checkProduction = FALSE, simParamBee = simParamBee, ...) + mapCasteToColonyValue(x, value = "dd", checkProduction = FALSE, simParamBee = simParamBee, ...) } #' @describeIn mapCasteToColonyValue Map caste member (individual) epistasis values to a colony epistasis value -mapCasteToColonyAa <- function(colony, simParamBee = NULL, ...) { +mapCasteToColonyAa <- function(x, simParamBee = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } - mapCasteToColonyValue(colony, value = "aa", checkProduction = FALSE, simParamBee = simParamBee, ...) + mapCasteToColonyValue(x, value = "aa", checkProduction = FALSE, simParamBee = simParamBee, ...) } diff --git a/R/Functions_L0_auxilary.R b/R/Functions_L0_auxilary.R index 28e56689..fbbb20b5 100644 --- a/R/Functions_L0_auxilary.R +++ b/R/Functions_L0_auxilary.R @@ -359,11 +359,11 @@ pHomBrood <- function(x, simParamBee = NULL) { } } } else if (isColony(x)) { - if (is.null(x@queen@misc$pHomBrood[[1]])) { - ret <- NA - } else { - ret <- x@queen@misc$pHomBrood[[1]] - } + if (is.null(x@queen@misc$pHomBrood[[1]])) { + ret <- NA + } else { + ret <- x@queen@misc$pHomBrood[[1]] + } } else if (isMultiColony(x)) { ret <- sapply(X = x@colonies, FUN = pHomBrood) names(ret) <- getId(x) @@ -952,150 +952,6 @@ isNULLColonies <- function(multicolony) { # get (general) ---- -#' @rdname getQueenYearOfBirth -#' @title Access the queen's year of birth -#' -#' @description Level 0 function that returns the queen's year of birth. -#' -#' @param x \code{\link[AlphaSimR]{Pop-class}} (one or more than one queen), -#' \code{\link[SIMplyBee]{Colony-class}} (one colony), or -#' \code{\link[SIMplyBee]{MultiColony-class}} (more colonies) -#' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters -#' -#' @return numeric, the year of birth of the queen(s); named when theres is more -#' than one queen; \code{NA} if queen not present -#' -#' @examples -#' founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) -#' SP <- SimParamBee$new(founderGenomes) -#' \dontshow{SP$nThreads = 1L} -#' basePop <- createVirginQueens(founderGenomes) -#' -#' drones <- createDrones(x = basePop[1], nInd = 1000) -#' droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson) -#' -#' # Create a Colony and a MultiColony class -#' colony <- createColony(x = basePop[2]) -#' colony <- cross(colony, drones = droneGroups[[1]]) -#' -#' apiary <- createMultiColony(basePop[3:4], n = 2) -#' apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) -#' -#' queen <- getQueen(colony) -#' queen <- setQueensYearOfBirth(queen, year = 2022) -#' getQueenYearOfBirth(queen) -#' -#' getQueenYearOfBirth(getQueen(colony)) -#' colony <- setQueensYearOfBirth(colony, year = 2030) -#' getQueenYearOfBirth(colony) -#' -#' apiary <- setQueensYearOfBirth(apiary, year = 2022) -#' getQueenYearOfBirth(apiary) -#' @export -getQueenYearOfBirth <- function(x, simParamBee = NULL) { - if (is.null(simParamBee)) { - simParamBee <- get(x = "SP", envir = .GlobalEnv) - } - if (isPop(x)) { - if (any(!(isVirginQueen(x, simParamBee = simParamBee) | isQueen(x, simParamBee = simParamBee)))) { - stop("Individuals in x must be virgin queens or queens!") - } - nInd <- nInd(x) - ret <- rep(x = NA, times = nInd) - for (ind in seq_len(nInd)) { - if (!is.null(x@misc$yearOfBirth[[ind]])) { - ret[ind] <- x@misc$yearOfBirth[[ind]] - } - } - if (nInd > 1) { - names(ret) <- getId(x) - } - } else if (isColony(x)) { - ret <- ifelse(is.null(x@queen@misc$yearOfBirth[[1]]), NA, x@queen@misc$yearOfBirth[[1]]) - } else if (isMultiColony(x)) { - ret <- sapply(X = x@colonies, FUN = getQueenYearOfBirth, simParamBee = simParamBee) - names(ret) <- getId(x) - } else { - stop("Argument x must be a Pop, Colony, or MultiColony class object!") - } - return(ret) -} - -#' @rdname getQueenAge -#' @title Get (calculate) the queen's age -#' -#' @description Level 0 function that returns the queen's age. -#' -#' @param x \code{\link[AlphaSimR]{Pop-class}}, \code{\link[SIMplyBee]{Colony-class}}, or -#' \code{\link[SIMplyBee]{MultiColony-class}} -#' @param currentYear integer, current year -#' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters -#' -#' @return numeric, the age of the queen(s); named when theres is more -#' than one queen; \code{NA} if queen not present -#' -#' @examples -#' founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) -#' SP <- SimParamBee$new(founderGenomes) -#' \dontshow{SP$nThreads = 1L} -#' basePop <- createVirginQueens(founderGenomes) -#' -#' drones <- createDrones(x = basePop[1], nInd = 1000) -#' droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson) -#' -#' # Create a Colony and a MultiColony class -#' colony <- createColony(x = basePop[2]) -#' colony <- cross(colony, drones = droneGroups[[1]]) -#' apiary <- createMultiColony(basePop[3:4], n = 2) -#' apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) -#' -#' queen <- getQueen(colony) -#' queen <- setQueensYearOfBirth(queen, year = 2020) -#' getQueenAge(queen, currentYear = 2022) -#' -#' colony <- setQueensYearOfBirth(colony, year = 2021) -#' getQueenAge(colony, currentYear = 2022) -#' -#' apiary <- setQueensYearOfBirth(apiary, year = 2018) -#' getQueenAge(apiary, currentYear = 2022) -#' @export -getQueenAge <- function(x, currentYear, simParamBee = NULL) { - if (is.null(simParamBee)) { - simParamBee <- get(x = "SP", envir = .GlobalEnv) - } - if (isPop(x)) { - if (any(!(isVirginQueen(x, simParamBee = simParamBee) | isQueen(x, simParamBee = simParamBee)))) { - stop("Individuals in x must be virgin queens or queens!") - } - nInd <- nInd(x) - ret <- rep(x = NA, times = nInd) - for (ind in seq_len(nInd)) { - if (!is.null(x@misc$yearOfBirth[[ind]])) { - ret[ind] <- currentYear - x@misc$yearOfBirth[[ind]] - } - } - if (nInd > 1) { - names(ret) <- getId(x) - } - } else if (isColony(x)) { - if (isQueenPresent(x, simParamBee = simParamBee)) { - if(packageVersion("AlphaSimR") > package_version("1.5.3")){ - ret <- currentYear - x@queen@misc$yearOfBirth[[1]] - }else{ - ret <- currentYear - x@queen@misc[[1]]$yearOfBirth - } - } else { - ret <- NA - } - } else if (isMultiColony(x)) { - ret <- sapply(X = x@colonies, FUN = getQueenAge, currentYear = currentYear) - names(ret) <- getId(x) - } else { - stop("Argument x must be a Pop, Colony, or MultiColony class object!") - } - return(ret) -} - #' @rdname getId #' @title Get the colony ID #' @@ -1769,7 +1625,7 @@ getEvents <- function(x) { #' colony <- buildUp(x = colony, nWorkers = 6, nDrones = 3) #' colony <- addVirginQueens(colony, nInd = 5) #' -#' apiary <- createMultiColony(basePop[3:4], n = 2) +#' apiary <- createMultiColony(basePop[3:4]) #' apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) #' apiary <- buildUp(x = apiary, nWorkers = 6, nDrones = 3) #' @@ -2342,6 +2198,9 @@ getQueenCsdAlleles <- function(x, allele = "all", unique = FALSE, collapse = FAL if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (!(isColony(x) | isMultiColony(x))) { + stop("Argument x must be a Colony or MultiColony class object!") + } ret <- getCsdAlleles(x, caste = "queen", allele = allele, @@ -2359,6 +2218,9 @@ getFathersCsdAlleles <- function(x, nInd = NULL, allele = "all", dronesHaploid = if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (!(isColony(x) | isMultiColony(x))) { + stop("Argument x must be a Colony or MultiColony class object!") + } ret <- getCsdAlleles(x, caste = "fathers", nInd = nInd, @@ -2378,6 +2240,9 @@ getVirginQueensCsdAlleles <- function(x, nInd = NULL, allele = "all", if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (!(isColony(x) | isMultiColony(x))) { + stop("Argument x must be a Colony or MultiColony class object!") + } ret <- getCsdAlleles(x, caste = "virginQueens", nInd = nInd, @@ -2396,6 +2261,9 @@ getWorkersCsdAlleles <- function(x, nInd = NULL, allele = "all", if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (!(isColony(x) | isMultiColony(x))) { + stop("Argument x must be a Colony or MultiColony class object!") + } ret <- getCsdAlleles(x, caste = "workers", nInd = nInd, @@ -2414,6 +2282,9 @@ getDronesCsdAlleles <- function(x, nInd = NULL, allele = "all", dronesHaploid = if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (!(isColony(x) | isMultiColony(x))) { + stop("Argument x must be a Colony or MultiColony class object!") + } ret <- getCsdAlleles(x, caste = "drones", nInd = nInd, @@ -2545,7 +2416,7 @@ getCsdGeno <- function(x, caste = NULL, nInd = NULL, dronesHaploid = TRUE, } else { ret <- getCsdGeno( x = getCastePop(x, caste, simParamBee = simParamBee), nInd = nInd, - dronesHaploid = dronesHaploid, simParamBee = simParamBee + dronesHaploid = dronesHaploid, simParamBee = simParamBee ) } } else if (isMultiColony(x)) { @@ -2581,6 +2452,9 @@ getQueenCsdGeno <- function(x, collapse = FALSE, simParamBee = NULL) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (!(isColony(x) | isMultiColony(x))) { + stop("Argument x must be a Colony or MultiColony class object!") + } ret <- getCsdAlleles(x, caste = "queen", collapse = collapse, @@ -2596,6 +2470,9 @@ getFathersCsdGeno <- function(x, nInd = NULL, dronesHaploid = TRUE, collapse = F if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (!(isColony(x) | isMultiColony(x))) { + stop("Argument x must be a Colony or MultiColony class object!") + } ret <- getCsdAlleles(x, nInd = nInd, caste = "fathers", @@ -2613,6 +2490,9 @@ getVirginQueensCsdGeno <- function(x, nInd = NULL, collapse = FALSE, if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (!(isColony(x) | isMultiColony(x))) { + stop("Argument x must be a Colony or MultiColony class object!") + } ret <- getCsdAlleles(x, nInd = nInd, caste = "virginQueens", @@ -2629,6 +2509,9 @@ getWorkersCsdGeno <- function(x, nInd = NULL, collapse = FALSE, if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (!(isColony(x) | isMultiColony(x))) { + stop("Argument x must be a Colony or MultiColony class object!") + } ret <- getCsdAlleles(x, nInd = nInd, caste = "workers", @@ -2646,6 +2529,9 @@ getDronesCsdGeno <- function(x, nInd = NULL, dronesHaploid = TRUE, if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (!(isColony(x) | isMultiColony(x))) { + stop("Argument x must be a Colony or MultiColony class object!") + } ret <- getCsdAlleles(x, nInd = nInd, caste = "drones", @@ -2679,7 +2565,8 @@ isGenoHeterozygous <- function(x) { #' #' @description Level 0 function that returns if individuals of a population are #' heterozygous at the csd locus. See \code{\link[SIMplyBee]{SimParamBee}} for more -#' information about the csd locus. +#' information about the csd locus. The function also return \code{TRUE} for drones to +#' mark their viability, although they are haploid. #' #' @param pop \code{\link[AlphaSimR]{Pop-class}} #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters @@ -2881,10 +2768,10 @@ nCsdAlleles <- function(x, collapse = FALSE, simParamBee = NULL) { #' \code{\link[SIMplyBee]{MultiColony-class}} #' #' @examples -#' founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 50) -#' SP <- SimParamBee$new(founderGenomes) +#' founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 5) +#' SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 4) #' \dontshow{SP$nThreads = 1L} -#' SP$setTrackRec(TRUE) +#' SP$setTrackRec(isTrackRec = TRUE) #' SP$setTrackPed(isTrackPed = TRUE) #' basePop <- createVirginQueens(founderGenomes) #' @@ -2894,13 +2781,13 @@ nCsdAlleles <- function(x, collapse = FALSE, simParamBee = NULL) { #' # Create a Colony and a MultiColony class #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) -#' colony <- buildUp(x = colony, nWorkers = 6, nDrones = 3) -#' colony <- addVirginQueens(x = colony, nInd = 5) +#' colony <- buildUp(x = colony, nWorkers = 3, nDrones = 2) +#' colony <- addVirginQueens(x = colony, nInd = 2) #' -#' apiary <- createMultiColony(basePop[3:4], n = 2) +#' apiary <- createMultiColony(basePop[3:4]) #' apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) -#' apiary <- buildUp(x = apiary, nWorkers = 6, nDrones = 3) -#' apiary <- addVirginQueens(x = apiary, nInd = 5) +#' apiary <- buildUp(x = apiary, nWorkers = 3, nDrones = 2) +#' apiary <- addVirginQueens(x = apiary, nInd = 2) #' #' # Input is a population #' getIbdHaplo(x = getQueen(colony)) @@ -2912,6 +2799,8 @@ nCsdAlleles <- function(x, collapse = FALSE, simParamBee = NULL) { #' getQueenIbdHaplo(colony) #' #' getIbdHaplo(colony, caste = "workers", nInd = 3) +#' getIbdHaplo(colony, caste = "virginQueens") +#' getIbdHaplo(colony, caste = "drones") #' getWorkersIbdHaplo(colony) #' # Same aliases exist for all castes! #' @@ -2926,6 +2815,9 @@ nCsdAlleles <- function(x, collapse = FALSE, simParamBee = NULL) { #' # Or collapse all the haplotypes into a single matrix #' getQueenIbdHaplo(apiary, collapse = TRUE) #' +#' +#' getIbdHaplo(x = apiary, caste = "workers") +#' getIbdHaplo(x = apiary, caste = "drones") #' # Get the haplotypes of all individuals either by colony or in a single matrix #' getIbdHaplo(apiary, caste = "all") #' getIbdHaplo(apiary, caste = "all", collapse = TRUE) @@ -5226,14 +5118,14 @@ calcColonyValue <- function(x, FUN = NULL, simParamBee = NULL, ...) { stop("You must provide FUN or set it in the SimParamBee object!") } if (isColony(x)) { - ret <- FUN(colony = x, ...) + ret <- FUN(x = x, ...) } else if (isMultiColony(x)) { nCol <- nColonies(x) # We could create a matrix output container here, BUT we don't know the output # dimension of FUN() so we create list and row bind the list nodes later ret <- vector(mode = "list", length = nCol) for (colony in seq_len(nCol)) { - ret[[colony]] <- FUN(colony = x[[colony]], ...) + ret[[colony]] <- FUN(x = x[[colony]], ...) } ret <- do.call("rbind", ret) rownames(ret) <- getId(x) @@ -6338,7 +6230,7 @@ editCsdLocus <- function(pop, alleles = NULL, simParamBee = NULL) { alleles <- expand.grid(as.data.frame(matrix(rep(0:1, length(csdSites)), nrow = 2, byrow = FALSE))) # Sample two different alleles (without replacement) for each individual nAlleles <- simParamBee$nCsdAlleles - alleles <- sapply(seq_len(pop@nInd), FUN = function(x) list(alleles[sample(nAlleles, size = 2, replace = F), ])) + alleles <- sapply(seq_len(pop@nInd), FUN = function(x) list(alleles[sample(nAlleles, size = 2, replace = FALSE), ])) } if (pop@nInd != length(alleles)) { @@ -6406,10 +6298,6 @@ editCsdLocus <- function(pop, alleles = NULL, simParamBee = NULL) { #' virginColonies2 <- setLocation(virginColonies2, #' location = Map(c, runif(30, 0, 2*pi), #' runif(30, 0, 2*pi))) -#' virginColonies3 <- createMultiColony(basePop[61:90]) -#' virginColonies3 <- setLocation(virginColonies3, -#' location = Map(c, runif(30, 0, 2*pi), -#' runif(30, 0, 2*pi))) #' #' # Create drone colonies #' droneColonies <- createMultiColony(basePop[121:200]) @@ -6429,58 +6317,23 @@ editCsdLocus <- function(pop, alleles = NULL, simParamBee = NULL) { #' nDrones = nFathersPoisson, #' crossPlan = randomCrossPlan) #' -#' # Plot the colonies in space -#' virginLocations <- as.data.frame(getLocation(c(virginColonies1, virginColonies2, virginColonies3), -#' collapse= TRUE)) -#' virginLocations$Type <- "Virgin" -#' droneLocations <- as.data.frame(getLocation(droneColonies, collapse= TRUE)) -#' droneLocations$Type <- "Drone" -#' locations <- rbind(virginLocations, droneLocations) -#' -#' plot(x = locations$V1, y = locations$V2, -#' col = c("red", "blue")[as.numeric(as.factor(locations$Type))]) -#' -#' # Cross according to a spatial cross plan according to the colonies' locations -#' crossPlanSpatial <- createCrossPlan(x = virginColonies1, -#' droneColonies = droneColonies, -#' nDrones = nFathersPoisson, -#' spatial = TRUE, -#' radius = 1.5) -#' -#' # Plot the crossing for the first colony in the crossPlan -#' virginLocations1 <- as.data.frame(getLocation(virginColonies1, collapse= TRUE)) -#' virginLocations1$Type <- "Virgin" -#' droneLocations <- as.data.frame(getLocation(droneColonies, collapse= TRUE)) -#' droneLocations$Type <- "Drone" -#' locations1 <- rbind(virginLocations1, droneLocations) -#' -#' # Blue marks the target virgin colony and blue marks the drone colonies in the chosen radius -#' plot(x = locations1$V1, y = locations1$V2, pch = c(1, 2)[as.numeric(as.factor(locations1$Type))], -#' col = ifelse(rownames(locations1) %in% crossPlanSpatial[[1]], -#' "red", -#' ifelse(rownames(locations1) == names(crossPlanSpatial)[[1]], -#' "blue", "black"))) -#' -#' colonies1 <- cross(x = virginColonies1, -#' crossPlan = crossPlanSpatial, -#' droneColonies = droneColonies, -#' nDrones = nFathersPoisson) -#' nFathers(colonies1) -#' #' # Cross according to a cross plan that is created internally within the cross function #' # The cross plan is created at random, regardless the location of the colonies -#' colonies2 <- cross(x = virginColonies2, +#' colonies1 <- cross(x = virginColonies1, #' droneColonies = droneColonies, #' nDrones = nFathersPoisson, #' crossPlan = "create") #' -#' # Mate spatially with cross plan created internally by the cross function -#' colonies3 <- cross(x = virginColonies3, -#' droneColonies = droneColonies, +#' +#' # Cross according to a spatial cross plan created internally according to the colonies' locations +#' colonies2 <- cross(x = virginColonies2, #' crossPlan = "create", -#' checkCross = "warning", +#' droneColonies = droneColonies, +#' nDrones = nFathersPoisson, #' spatial = TRUE, -#' radius = 1) +#' radius = 1.5) +#' nFathers(colonies2) +#' #' #' @export createCrossPlan <- function(x, @@ -6548,93 +6401,6 @@ createCrossPlan <- function(x, return(crossPlan) } -# Misc helpers -# These functions replace the defunct functions of the same name in AlphaSimR - -#' @rdname setMisc -#' @title Set miscellaneous information in a population -#' -#' @description Set miscellaneous information in a population -#' -#' @param x \code{\link[AlphaSimR]{Pop-class}} -#' @param node character, name of the node to set within the \code{x@misc} slot -#' @param value, value to be saved into \code{x@misc[[*]][[node]]}; length of -#' \code{value} should be equal to \code{nInd(x)}; if its length is 1, then -#' it is repeated using \code{rep} (see examples) -#' -#' @details A \code{NULL} in \code{value} is ignored -#' -#' @return \code{\link[AlphaSimR]{Pop-class}} -#' -#' @export -setMisc <- function(x, node = NULL, value = NULL) { - if (isPop(x)) { - if (is.null(node)) { - stop("Argument node must be provided!") - } - if (is.null(value)) { - stop("Argument value must be provided!") - } - n <- nInd(x) - if (length(value) == 1 && n > 1) { - value <- rep(x = value, times = n) - } - if (length(value) != n) { - stop("Argument value must be of length 1 or nInd(x)!") - } - - # Check current AlphaSimR version for new or legacy misc slot - if(packageVersion("AlphaSimR") > package_version("1.5.3")){ - # New misc slot - x@misc[[node]] = value - }else{ - # Legacy misc slot - names(value) = rep(x = node, times = n) - inode = match(names(x@misc[[1]]),node) - inode = inode[!is.na(inode)] - if(length(inode) == 0){ - x@misc = sapply(seq_len(n),function(ind){ - c(x@misc[[ind]],value[ind]) - },simplify = FALSE) - }else{ - x@misc = sapply(seq_len(n),function(ind){ - c(x@misc[[ind]],value[ind])[-inode] - },simplify = FALSE) - } - } - - } - - return(x) -} - -#' @rdname getMisc -#' @title Get miscellaneous information in a population -#' -#' @description Get miscellaneous information in a population -#' -#' @param x \code{\link[AlphaSimR]{Pop-class}} -#' @param node character, name of the node to get from the \code{x@misc} slot; -#' if \code{NULL} the whole \code{x@misc} slot is returned -#' -#' @return The \code{x@misc} slot or its nodes \code{x@misc[[*]][[node]]} -#' -#' @export -getMisc <- function(x, node = NULL) { - if (isPop(x)) { - if (is.null(node)) { - ret <- x@misc - } else { - # Check current AlphaSimR version for new or legacy misc slot - ret = x@misc[[node]] - } - } else { - stop("Argument x must be a Pop class object!") - } - return(ret) -} - - #' @rdname mapLoci #' @title Finds loci on a genetic map and return a list of positions #' diff --git a/R/Functions_L1_Pop.R b/R/Functions_L1_Pop.R index 005fe4fd..2503270f 100644 --- a/R/Functions_L1_Pop.R +++ b/R/Functions_L1_Pop.R @@ -1,4 +1,8 @@ # ---- Level 1 Pop Functions ---- +utils::globalVariables("colony") +utils::globalVariables("i") +utils::globalVariables("cl") + #' @rdname getCastePop #' @title Access individuals of a caste @@ -218,18 +222,18 @@ getFathers <- function(x, nInd = NULL, use = "rand", collapse = FALSE, simParamB if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } - if (isPop(x)) { # DO WE WANT TO PUT THIS IN getCastePop??? - ret = lapply(X = x@misc$fathers, - FUN = function(z){ - if(is.null(z)){ - ret = NULL - }else{ - if (is.null(nInd)) { - n <- nInd(z) - } - ret <- selectInd(pop = z, nInd = n, use = use, simParam = simParamBee) - } - } + if (isPop(x)) { # TODO: DO WE WANT TO PUT THIS IN getCastePop??? + ret <- lapply(X = x@misc$fathers, + FUN = function(z){ + if(is.null(z)){ + ret = NULL + }else{ + if (is.null(nInd)) { + n <- nInd(z) + } + ret <- selectInd(pop = z, nInd = n, use = use, simParam = simParamBee) + } + } ) if (nInd(x) == 1) { ret <- ret[[1]] @@ -295,11 +299,6 @@ getVirginQueens <- function(x, nInd = NULL, use = "rand", collapse = FALSE, simP #' only used when \code{x} is \code{\link[SIMplyBee]{Colony-class}} or #' \code{\link[SIMplyBee]{MultiColony-class}}, when \code{x} is \code{link[AlphaSimR]{MapPop-class}} #' all individuals in \code{x} are converted into virgin queens -#' @param exact logical, only relevant when creating workers, -#' if the csd locus is active and exact is \code{TRUE}, -#' create the exactly specified number of viable workers (heterozygous on the -#' csd locus) -#' @param year numeric, year of birth for virgin queens #' @param editCsd logical (only active when \code{x} is \code{link[AlphaSimR]{MapPop-class}}), #' whether the csd locus should be edited to ensure heterozygosity at the csd #' locus (to get viable virgin queens); see \code{csdAlleles} @@ -313,6 +312,10 @@ getVirginQueens <- function(x, nInd = NULL, use = "rand", collapse = FALSE, simP #' in \code{\link[SIMplyBee]{SimParamBee}}. The two csd alleles must be different to #' ensure heterozygosity at the csd locus. #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters +#' @param returnSP logical, whether to return the pedigree, caste, and recHist information +#' for each created population (used internally for parallel computing) +#' @param ids character, IDs of the individuals that are going to be created (used internally +#' for parallel computing) #' @param ... additional arguments passed to \code{nInd} when this argument is a function #' #' @return when \code{x} is \code{\link[AlphaSimR]{MapPop-class}} returns @@ -327,7 +330,7 @@ getVirginQueens <- function(x, nInd = NULL, use = "rand", collapse = FALSE, simP #' founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 50) #' SP <- SimParamBee$new(founderGenomes) #' \dontshow{SP$nThreads = 1L} -#' SP$setTrackRec(TRUE) +#' SP$setTrackRec(isTrackRec = TRUE) #' SP$setTrackPed(isTrackPed = TRUE) #' #' # Create virgin queens on a MapPop @@ -345,7 +348,7 @@ getVirginQueens <- function(x, nInd = NULL, use = "rand", collapse = FALSE, simP #' # Create a Colony and a MultiColony class #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) -#' apiary <- createMultiColony(basePop[3:4], n = 2) +#' apiary <- createMultiColony(basePop[3:4]) #' apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) #' #' # Using default nInd in SP @@ -398,13 +401,15 @@ getVirginQueens <- function(x, nInd = NULL, use = "rand", collapse = FALSE, simP # patrilines # https://github.com/HighlanderLab/SIMplyBee/issues/78 createCastePop <- function(x, caste = NULL, nInd = NULL, - exact = TRUE, year = NULL, editCsd = TRUE, csdAlleles = NULL, simParamBee = NULL, + returnSP = FALSE, + ids = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (is.null(nInd)) { if (caste == "virginQueens") { nInd <- simParamBee$nVirginQueens @@ -417,6 +422,9 @@ createCastePop <- function(x, caste = NULL, nInd = NULL, if (is.function(nInd)) { nInd <- nInd(x, ...) } + if (any(nInd == 0)) { + stop("nInd set to 0, should be > 0!") + } # doing "if (is.function(nInd))" below if (isMapPop(x)) { if (caste != "virginQueens") { # Creating virgin queens if input is a MapPop @@ -430,9 +438,7 @@ createCastePop <- function(x, caste = NULL, nInd = NULL, } ret@sex[] <- "F" simParamBee$changeCaste(id = ret@id, caste = "virginQueens") - if (!is.null(year)) { - ret <- setQueensYearOfBirth(x = ret, year = year, simParamBee = simParamBee) - } + } else if (isPop(x)) { if (caste != "drones") { # Creating drones if input is a Pop stop("Pop-class can only be used to create drones!") @@ -453,73 +459,151 @@ createCastePop <- function(x, caste = NULL, nInd = NULL, } ret <- list() for (virginQueen in 1:nInd(x)) { - ret[[virginQueen]] <- makeDH(pop = x[virginQueen], nDH = nInd[virginQueen], keepParents = FALSE, simParam = simParamBee) + ret[[virginQueen]] <- makeDH(pop = x[virginQueen], nDH = nInd[virginQueen], + keepParents = FALSE, simParam = simParamBee) } ret <- mergePops(ret) } ret@sex[] <- "M" simParamBee$addToCaste(id = ret@id, caste = "drones") } else if (isColony(x)) { - if (!isQueenPresent(x, simParamBee = simParamBee)) { - stop("Missing queen!") - } + originalThreads <- simParamBee$nThreads + simParamBee$nThreads <- 1 + if (length(nInd) > 1) { warning("More than one value in the nInd argument, taking only the first value!") nInd <- nInd[1] } - if (caste == "workers") { - ret <- vector(mode = "list", length = 2) - names(ret) <- c("workers", "nHomBrood") - workers <- combineBeeGametes( - queen = getQueen(x, simParamBee = simParamBee), drones = getFathers(x, simParamBee = simParamBee), - nProgeny = nInd, simParamBee = simParamBee - ) - if (isCsdActive(simParamBee = simParamBee)) { - sel <- isCsdHeterozygous(pop = workers, simParamBee = simParamBee) - ret$workers <- workers[sel] - ret$nHomBrood <- nInd - sum(sel) - if (exact) { - if (nInd(ret$workers) < nInd) { - nMiss <- nInd - nInd(ret$workers) - while (0 < nMiss) { - workers <- combineBeeGametes( - queen = getQueen(x, simParamBee = simParamBee), - drones = getFathers(x, simParamBee = simParamBee), - nProgeny = nMiss, - simParamBee = simParamBee - ) - sel <- isCsdHeterozygous(pop = workers, simParamBee = simParamBee) - ret$workers <- c(ret$workers, workers[sel]) - ret$nHomBrood <- ret$nHomBrood + sum(!sel) - nMiss <- nInd - nInd(ret$workers) + if (!isQueenPresent(x, simParamBee = simParamBee)) { + stop("Missing queen!") + } + + if (nInd > 0) { + if (caste == "workers") { + if (!returnSP) { + ret <- vector(mode = "list", length = 2) + names(ret) <- c("workers", "nHomBrood") + } else { + ret <- vector(mode = "list", length = 5) + names(ret) <- c("workers", "nHomBrood", "pedigree", "caste", "recHist") + } + simParamBee$nThreads <- 1 + ret$workers <- combineBeeGametes( + queen = getQueen(x, simParamBee = simParamBee), + drones = getFathers(x, simParamBee = simParamBee), + nProgeny = nInd, + simParamBee = simParamBee + ) + + simParamBee$addToCaste(id = ret$workers@id, caste = "workers") + ret$workers@sex[] <- "F" + + if (returnSP) { + ret$caste <- simParamBee$caste[ret$workers@id, drop = FALSE] + if (simParamBee$isTrackPed) { + ret$pedigree <- simParamBee$pedigree[ret$workers@id, , drop = FALSE] + } + if (simParamBee$isTrackRec) { + ret$recHist <- simParamBee$recHist[ret$workers@iid] + } + } + + if (!is.null(ids)) { + if (nInd(ret$workers) > length(ids)) { + stop("Not enough IDs provided") + } + if (nInd(ret$workers) < length(ids)) { + stop("Too many IDs provided!") + } + ret$workers@id <- as.character(ids) + ret$workers@iid <- as.integer(ids) + if (returnSP) { + names(ret$caste) <- as.character(ids) + if (simParamBee$isTrackPed) { + rownames(ret$pedigree) <- ids + } + if (simParamBee$isTrackRec) { + names(ret$recHist) <- ids } } } - } else { - ret$workers <- workers - ret$nHomBrood <- NA - } - ret$workers@sex[] <- "F" - simParamBee$addToCaste(id = ret$workers@id, caste = "workers") - } else if (caste == "virginQueens") { # Creating virgin queens if input is a Colony - ret <- createCastePop(x = x, caste = "workers", - nInd = nInd, exact = TRUE, simParamBee = simParamBee)$workers - ret@sex[] <- "F" - simParamBee$changeCaste(id = ret@id, caste = "virginQueens") - if (!is.null(year)) { - ret <- setQueensYearOfBirth(x = ret, year = year, simParamBee = simParamBee) + + if (isCsdActive(simParamBee = simParamBee)) { + sel <- isCsdHeterozygous(pop = ret$workers, simParamBee = simParamBee) + ret$nHomBrood <- nInd(ret$workers) - sum(sel) + ret$workers <- ret$workers[sel] + } else { + ret$nHomBrood <- NA + } + + } else if (caste == "virginQueens") { + ret <- createCastePop(x = x, caste = "workers", + nInd = nInd, simParamBee = simParamBee, + returnSP = returnSP, ids = ids, ...) + ret$caste = rep("virginQueens", length(ret$caste)) + names(ret$caste) = ids + simParamBee$changeCaste(id = ret$workers@id, caste = "virginQueens") + if (!returnSP) { + ret <- ret$workers + } + + } else if (caste == "drones") { + drones <- makeDH( + pop = getQueen(x, simParamBee = simParamBee), nDH = nInd, keepParents = FALSE, + simParam = simParamBee + ) + + drones@sex[] <- "M" + simParamBee$addToCaste(id = drones@id, caste = "drones") + + if (returnSP) { + ret <- vector(mode = "list", length = 4) + names(ret) <- c("drones", "pedigree", "caste", "recHist") + ret$caste <- simParamBee$caste[drones@id, drop = FALSE] + if (simParamBee$isTrackPed) { + ret$pedigree <- simParamBee$pedigree[drones@id, , drop = FALSE] + } + if (simParamBee$isTrackRec) { + ret$recHist <- simParamBee$recHist[drones@iid] + } + } + + if (!is.null(ids)) { + if (nInd(drones) != length(ids)) { + stop("Not enough IDs provided") + } + drones@id = ids + drones@iid = as.integer(ids) + if (returnSP) { + names(ret$caste) <- ids + if (simParamBee$isTrackPed) { + rownames(ret$pedigree) <- ids + } + if (simParamBee$isTrackRec) { + names(ret$recHist) <- ids + } + } + } + + if (returnSP) { + ret$drones <-drones + } else { + ret <- drones + } } - } else if (caste == "drones") { # Creating drones if input is a Colony - ret <- makeDH( - pop = getQueen(x, simParamBee = simParamBee), nDH = nInd, keepParents = FALSE, - simParam = simParamBee - ) - ret@sex[] <- "M" - simParamBee$addToCaste(id = ret@id, caste = "drones") + } else { + ret <- NULL } + simParamBee$nThreads <- originalThreads } else if (isMultiColony(x)) { + if (is.null(nInd)) { + string = paste0("n", toupper(substr(caste, 1, 1)), substr(caste, 2, nchar(caste))) + nInd <- simParamBee[[string]] + } + nCol <- nColonies(x) nNInd <- length(nInd) + if (nNInd > 1 && nNInd < nCol) { stop("Too few values in the nInd argument!") } @@ -527,56 +611,147 @@ createCastePop <- function(x, caste = NULL, nInd = NULL, warning(paste0("Too many values in the nInd argument, taking only the first ", nCol, "values!")) nInd <- nInd[1:nCol] } - ret <- vector(mode = "list", length = nCol) - for (colony in seq_len(nCol)) { - if (is.null(nInd)) { - nIndColony <- NULL + + nNInd <- length(nInd) + totalNInd = ifelse(nNInd == 1, nInd * nCol, sum(nInd)) + if (totalNInd == 0) { + stop("Nothing to create.") + } + + lastId = simParamBee$lastId + ids = (lastId+1):(lastId+totalNInd) + + combine_list <- function(a, b) { + if (!is.null(names(a))) { + "Combine first" + c(list(a), list(b)) } else { - nIndColony <- ifelse(nNInd == 1, nInd, nInd[colony]) + if ((is.null(a) | is.null(b)) & !(is.null(a) & is.null(b))) { + c(a, list(b)) + } else if (is.null(a) & is.null(b)) { + c(list(a), list(b)) + } else { + c(a, list(b)) + } } - ret[[colony]] <- createCastePop( - x = x[[colony]], caste = caste, - nInd = nIndColony, - exact = exact, - year = year, - editCsd = TRUE, csdAlleles = NULL, - simParamBee = simParamBee, ... - ) + } + + ret <- foreach(colony = seq_len(nCol), .combine=combine_list, .packages = c("SIMplyBee")) %dopar% { + nIndColony <- ifelse(nNInd == 1, nInd, nInd[colony]) + if (nIndColony > 0) { + if (nNInd == 1) { + colonyIds = ids[((colony-1)*nIndColony+1):(colony*nIndColony)] + } else { + colonyIds = base::split(ids, rep(seq_along(nInd), nInd))[[as.character(colony)]] + } + createCastePop( + x = x[[colony]], caste = caste, + nInd = nIndColony, + editCsd = TRUE, csdAlleles = NULL, + simParamBee = simParamBee, + returnSP = TRUE, + ids = as.character(colonyIds) + ) + } else { + NULL + } + } + + if (nCol == 1) { + ret <- list(ret) } names(ret) <- getId(x) - } else { + # Add to simParamBee: pedigree, caste, recHist + notNull = sapply(ret, FUN = function(x) !is.null(x)) + + if (!simParamBee$isTrackPed) { + simParamBee$updateLastBeeId(n = totalNInd) + } else if (simParamBee$isTrackPed) { + Pedigree <- do.call("rbind", lapply(ret[notNull], '[[', "pedigree")) + if (!simParamBee$isTrackRec) { + simParamBee$addToBeePed(nNewInd = totalNInd, id = rownames(Pedigree), + mother = Pedigree[, 'mother'], father = Pedigree[, 'father'], + isDH = Pedigree[, 'isDH']) + #simParamBee$updatePedigree(pedigree = Pedigree) + } else { + RecHist = do.call("c", lapply(ret[notNull], '[[', "recHist")) + if (caste == "drones") { + ploidy = rep(1, totalNInd) + } else { + ploidy = rep(2, totalNInd) + } + simParamBee$addToBeeRec(nNewInd = totalNInd, id = rownames(Pedigree), + mother = Pedigree[, 'mother'], father = Pedigree[, 'father'], + isDH = Pedigree[, 'isDH'], + hist = RecHist, ploidy = ploidy) + } + } + # Extend caste + Caste <- do.call("c", lapply(ret[notNull], '[[', "caste")) + Names <- do.call("c", lapply(ret[notNull], function(x) names(x$caste))) + names(Caste) <- Names + simParamBee$updateCaste(caste = Caste) + + + if (!returnSP) { + if (caste %in% c("drones", "virginQueens")) { + ret <- lapply(ret, FUN = function(x) { + if (is.null(x)) return(NULL) # Return NULL if the element is NULL + x[!names(x) %in% c("pedigree", "caste", "recHist")][[1]] + }) + } else { + ret <- lapply(ret, FUN = function(x) { + if (is.null(x)) return(NULL) + x[!names(x) %in% c("pedigree", "caste", "recHist")] + }) + } + } + } + else { stop("Argument x must be a Map-Pop (only for virgin queens), Pop (only for drones), Colony, or MultiColony class object!") } + return(ret) } #' @describeIn createCastePop Create workers from a colony #' @export -createWorkers <- function(x, nInd = NULL, exact = FALSE, simParamBee = NULL, ...) { +createWorkers <- function(x, nInd = NULL, simParamBee = NULL, + returnSP = FALSE, + ids = NULL, ...) { ret <- createCastePop(x, caste = "workers", nInd = nInd, - exact = exact, simParamBee = simParamBee, ...) + simParamBee = simParamBee, + returnSP = returnSP, + ids = ids, ...) return(ret) } #' @describeIn createCastePop Create drones from a colony #' @export -createDrones <- function(x, nInd = NULL, simParamBee = NULL, ...) { +createDrones <- function(x, nInd = NULL, simParamBee = NULL, + returnSP = FALSE, + ids = NULL, ...) { ret <- createCastePop(x, caste = "drones", nInd = nInd, - simParamBee = simParamBee, ...) + simParamBee = simParamBee, + returnSP = returnSP, + ids = ids, ...) return(ret) } #' @describeIn createCastePop Create virgin queens from a colony #' @export createVirginQueens <- function(x, nInd = NULL, - year = NULL, editCsd = TRUE, csdAlleles = NULL, simParamBee = NULL, + returnSP = FALSE, + ids = NULL, ...) { ret <- createCastePop(x, caste = "virginQueens", nInd = nInd, - year = year, editCsd = editCsd, - csdAlleles = csdAlleles, simParamBee = simParamBee, ...) + editCsd = editCsd, + csdAlleles = csdAlleles, simParamBee = simParamBee, + returnSP = returnSP, + ids = ids, ...) return(ret) } @@ -817,7 +992,7 @@ createDCA <- function(x, nInd = NULL, removeFathers = TRUE, simParamBee = NULL) #' SP <- SimParamBee$new(founderGenomes) #' \dontshow{SP$nThreads = 1L} #' basePop <- createVirginQueens(founderGenomes) -#' drones <- createDrones(basePop[1], n = 1000) +#' drones <- createDrones(basePop[1], nInd = 1000) #' droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10) #' #' # Create a colony and cross it @@ -1011,12 +1186,12 @@ pullDroneGroupsFromDCA <- function(DCA, n, nDrones = NULL, #' # Create a Colony and a MultiColony class #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) -#' colony <- buildUp(x = colony, nWorkers = 100, nDrones = 10, exact = TRUE) +#' colony <- buildUp(x = colony, nWorkers = 100, nDrones = 10) #' colony <- addVirginQueens(x = colony, nInd = 3) #' #' apiary <- createMultiColony(basePop[3:4], n = 2) #' apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) -#' apiary <- buildUp(x = apiary, nWorkers = 100, nDrones = 10, exact = TRUE) +#' apiary <- buildUp(x = apiary, nWorkers = 100, nDrones = 10) #' apiary <- addVirginQueens(x = apiary, nInd = 3) #' #' # pullCastePop on Colony class @@ -1092,24 +1267,25 @@ pullCastePop <- function(x, caste, nInd = NULL, use = "rand", ret$pulled <- vector(mode = "list", length = nCol) names(ret$pulled) <- getId(x) ret$remnant <- x - for (colony in seq_len(nCol)) { + + tmp = foreach(colony = seq_len(nCol), .packages = c("SIMplyBee")) %dopar% { if (is.null(nInd)) { nIndColony <- NULL } else { nIndColony <- ifelse(nNInd == 1, nInd, nInd[colony]) } - tmp <- pullCastePop(x = x[[colony]], - caste = caste, - nInd = nIndColony, - use = use, - removeFathers = removeFathers, - collapse = collapse, - simParamBee = simParamBee) - if (!is.null(tmp$pulled)) { - ret$pulled[[colony]] <- tmp$pulled - } - ret$remnant[[colony]] <- tmp$remnant + pullCastePop(x = x[[colony]], + caste = caste, + nInd = nIndColony, + use = use, + removeFathers = removeFathers, + collapse = collapse, + simParamBee = simParamBee) } + + ret$pulled <- lapply(tmp, '[[', "pulled") + ret$remnant@colonies <- lapply(tmp, '[[', "remnant") + if (collapse) { ret$pulled <- mergePops(ret$pulled) } @@ -1195,7 +1371,8 @@ pullVirginQueens <- function(x, nInd = NULL, use = "rand", collapse = FALSE, sim #' to their distance from the virgin colony (that is, in a radius) #' @param radius numeric, the radius around the virgin colony in which to sample mating partners, #' only needed when \code{spatial = TRUE} -#' @param checkCross character, throw a warning (when \code{checkCross = "warning"}), +#' @param checkCross character, throw a warning (when \code{checkCross = "warning"}). +#' This will also remove the unmated queens and return only the mated ones. #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters #' @param ... other arguments for \code{nDrones}, when \code{nDrones} is a function #' @@ -1224,7 +1401,7 @@ pullVirginQueens <- function(x, nInd = NULL, use = "rand", collapse = FALSE, sim #' @examples #' founderGenomes <- quickHaplo(nInd = 30, nChr = 1, segSites = 100) #' SP <- SimParamBee$new(founderGenomes) -#' \dontshow{SP$nThreads = 1L} +#' SP$nThreads = 1L #' basePop <- createVirginQueens(founderGenomes) #' #' drones <- createDrones(x = basePop[1], nInd = 1000) @@ -1347,14 +1524,20 @@ cross <- function(x, if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + + if (isPop(x)) { + type = "Pop" + } else if (isColony(x)) { + type = "Colony" + } else if (isMultiColony(x)) { + type = "MultiColony" + } else { + stop("Input x must be a Pop-class, Colony-class, or MultiColony-class!") + } + if (is.null(nDrones)) { nDrones <- simParamBee$nFathers } - if (is.function(nDrones)) { - nD <- nDrones(...) - } else { - nD <- nDrones - } IDs <- as.character(getId(x)) oneColony <- (isPop(drones)) && (length(IDs) == 1) && (is.null(crossPlan)) @@ -1366,6 +1549,10 @@ cross <- function(x, # Do all the tests here to simplify the function + if (is.null(crossPlan) & (length(IDs) > 1) & isPop(drones)) { + stop("When supplying drones as a single population for mating multiple virgin queens, + crossPlan argument must be set to 'create' to internally create a mating plan!") + } if (crossPlan_droneID && !isPop(drones)) { stop("When using a cross plan, drones must be supplied as a single Pop-class!") } @@ -1382,6 +1569,9 @@ cross <- function(x, stop("The argument drones must be a Pop-class or a list of drone Pop-class objects!") } + if (isPop(drones) && nInd(drones) == 0) { + stop("Argument drones is a Pop-class with 0 individuals!") + } if (crossPlan_given && !is.null(drones) && !all(unlist(crossPlan) %in% drones@id)) { stop("Some drones from the crossPlan are missing in the drones population!") } @@ -1408,195 +1598,205 @@ cross <- function(x, } } + if (crossPlan_create | crossPlan_given) { + if (crossPlan_create) { + crossPlan <- createCrossPlan(x = x, + drones = drones, + droneColonies = droneColonies, + nDrones = nDrones, + spatial = spatial, + radius = radius, + simParamBee = simParamBee) + } - if (crossPlan_create) { - crossPlan <- createCrossPlan(x = x, - drones = drones, - droneColonies = droneColonies, - nDrones = nDrones, - spatial = spatial, - radius = radius, - simParamBee = simParamBee) noMatches <- sapply(crossPlan, FUN = length) + if (all(noMatches == 0)) { + stop("All crossings failed!") + } if (0 %in% noMatches) { - message("There are no potential crosses for some colonies! The cross() will fail - unless argument checkCross is set to 'warning'.") + if (checkCross == "warning") { + message("Crossing failed, unmated virgin queens will be removed!") + ret <- x + } else if (checkCross == "error") { + stop("Crossing failed!") + } } } - if (isPop(x) | isColony(x)) { - ret <- list() - for (virgin in seq_len(length(IDs))) { - virginID <- IDs[virgin] - if (oneColony) { - virginQueenDrones <- drones - } else if (dronePackages) { - virginQueenDrones <- drones[[virgin]] - } else if (crossPlan_given | crossPlan_create) { - if (crossPlan_droneID) { - virginQueenDrones <- drones[crossPlan[[virginID]]] - } else if (crossPlan_colonyID) { - virginMatches <- crossPlan[[virginID]] - if (length(virginMatches) > 0) { - nD <- ifelse(is.function(nDrones), nDrones(...), nDrones) - selectedDPQ <- table(sample(virginMatches, size = nD, replace = TRUE)) - virginQueenDrones <- mergePops(createDrones(droneColonies[names(selectedDPQ)], - nInd = selectedDPQ, simParamBee = simParamBee)) - } else { - virginQueenDrones <- new("Pop") - } - } - } - if (any((virginQueenDrones@nInd == 0), (length(virginQueenDrones@nInd) == 0))) { - msg <- "Crossing failed!" - if (checkCross == "warning") { - message(msg) - ret <- x - } else if (checkCross == "error") { - stop(msg) - } - } else if (virginQueenDrones@nInd > 0) { - if (!all(isDrone(virginQueenDrones, simParamBee = simParamBee))) { - stop("Individuals in drones must be drones!") - } - if (isPop(x)) { - virginQueen <- x[virgin] - } else if (isColony(x)) { - virginQueen <- selectInd(x@virginQueens, nInd = 1, use = "rand", simParam = simParamBee) - } + # Convert everything to a Pop + if (isColony(x) | isMultiColony(x)) { + inputId <- getId(x) + if (isColony(x)) { + colony <- x + x <- pullCastePop(x, caste = "virginQueens", nInd = 1, simParamBee = simParamBee)$pulled + ID_by_input <- data.frame(inputId = inputId, + virginId = getId(x)) + } else if (isMultiColony(x)) { + multicolony <- x + x <- pullCastePop(x, caste = "virginQueens", nInd = 1, simParamBee = simParamBee)$pulled + ID_by_input <- data.frame(inputId = inputId, + virginId = unlist(sapply(x, FUN = function(y) getId(y)))) + x <- mergePops(x) + } + # Rename crossPlan + if (crossPlan_create | crossPlan_given) { + names(crossPlan) <- ID_by_input$virginId[match(ID_by_input$inputId, names(crossPlan))] + } + } - virginQueen@misc$fathers[[1]] <- virginQueenDrones + IDs <- as.character(getId(x)) + #Now x is always a Pop + ret <- list() + nVirgin = nInd(x) - simParamBee$changeCaste(id = virginQueen@id, caste = "queen") - simParamBee$changeCaste(id = virginQueenDrones@id, caste = "fathers") - virginQueen <- setMisc(x = virginQueen, node = "nWorkers", value = 0) - virginQueen <- setMisc(x = virginQueen, node = "nDrones", value = 0) - virginQueen <- setMisc(x = virginQueen, node = "nHomBrood", value = 0) - if (isCsdActive(simParamBee = simParamBee)) { - val <- calcQueensPHomBrood(x = virginQueen, simParamBee = simParamBee) - } else { - val <- NA - } - virginQueen <- setMisc(x = virginQueen, node = "pHomBrood", value = val) - - if (isPop(x)) { - ret[[virgin]] <- virginQueen - } else if (isColony(x)) { - x <- reQueen(x = x, queen = virginQueen, simParamBee = simParamBee) - x <- removeVirginQueens(x, simParamBee = simParamBee) - ret <- x - } - } - } - if (isPop(x)) { - ret <- mergePops(ret) - } - } else if (isMultiColony(x)) { - nCol <- nColonies(x) - if (nCol == 0) { - ret <- createMultiColony(simParamBee = simParamBee) + + if (is.function(nDrones)) { + nD = nDrones(n = nVirgin, ...) + } else { + nD = nDrones + } + + if (length(IDs) > 0 & length(nD) == 1) { + nD = rep(nD, length(IDs)) + } + + + combine_list <- function(a, b) { + if (isPop(a)) { + c(list(a), list(b)) } else { - ret <- createMultiColony(n = nCol, simParamBee = simParamBee) - for (colony in seq_len(nCol)) { - if (oneColony) { - colonyDrones <- drones - } else if (dronePackages) { - colonyDrones <- drones[[colony]] - } else { - if (crossPlan_colonyID) { - colonyDrones <- NULL - } else if(crossPlan_droneID) { - colonyDrones <- drones - } - } - ret[[colony]] <- cross( - x = x[[colony]], - drones = colonyDrones, - crossPlan = crossPlan, - droneColonies = droneColonies, - nDrones = nDrones, - spatial = spatial, - radius = radius, - checkCross = checkCross, - simParamBee = simParamBee - ) + c(a, list(b)) + } + } + + if (crossPlan_given | crossPlan_create) { + if (crossPlan_colonyID) { # WHAT IF ONE ELEMENT IS EMPTY + # This is the crossPlan - for spatial, these are all DPCs found in a radius + crossPlanDF <- data.frame(virginID = rep(names(crossPlan), unlist(sapply(crossPlan, length))), + DPC = unlist(crossPlan)) + # If some of the crossing would fail, we only return the queens that mated successfully + IDs = IDs[IDs %in% crossPlanDF$virginID] + x = x[IDs] + if (type == "MultiColony") { + multicolony <- multicolony[getId(getVirginQueens(multicolony, collapse=TRUE)) %in% IDs] + } + # Here we sample from the DPC in the cross plan to get the needed number of drones (nD) + crossPlanDF_sample <- do.call("rbind", lapply(IDs, + FUN = function(x) { + data.frame(virginID = x, DPC = sample(crossPlan[[x]], size = nD[which(x == IDs)], replace = TRUE)) + } )) + crossPlanDF_sample <- crossPlanDF_sample[order(as.integer(crossPlanDF_sample$DPC)),] + # Here I gather how many drones each DPC needs to produce + crossPlanDF_DPCtable <- as.data.frame(table(crossPlanDF_sample$DPC)) + crossPlanDF_DPCtable <- crossPlanDF_DPCtable[order(as.integer(as.character(crossPlanDF_DPCtable$Var1))),] + colnames(crossPlanDF_DPCtable) <- c("DPC", "noDrones") + # Here I select only the DPCs that have been sampled to produce drones + selectedDPC = selectColonies(droneColonies, ID = as.character(crossPlanDF_DPCtable$DPC)) + # And here I create the drones + dronesByDPC <- createCastePop(selectedDPC, caste = "drones", + nInd = as.integer(crossPlanDF_DPCtable$noDrones), + simParamBee = simParamBee) + # This is where I link the drone ID to the DPC ID + dronesByDPC_DF <- data.frame(DPC = rep(names(dronesByDPC), as.vector(crossPlanDF_DPCtable$noDrones)), + droneID = unlist(sapply(dronesByDPC, FUN = function(x) getId(x)))) + dronesByDPC_DF <- dronesByDPC_DF[order(as.integer(dronesByDPC_DF$DPC)),] + dronePop = mergePops(dronesByDPC) + + if (any(!crossPlanDF_sample$DPC == dronesByDPC_DF$DPC)) { + stop("Something went wrong with cross plan - drone matching!") } + + dronesByVirgin_DF <- cbind(dronesByDPC_DF, crossPlanDF_sample[, c("virginID"), drop = FALSE]) + dronesByVirgin_DF <- dronesByVirgin_DF[order(as.integer(dronesByVirgin_DF$virginID)),] + dronesByVirgin_list <- lapply(IDs, + FUN = function(x) dronesByVirgin_DF$droneID[dronesByVirgin_DF$virginID == x]) + names(dronesByVirgin_list) <- IDs + + dronesByVirgin <- foreach(i = IDs, .combine = combine_list, + .packages = c("SIMplyBee")) %do% { + dronePop[as.character(dronesByVirgin_list[[i]])] + } + } else if (crossPlan_droneID) { + dronesByVirgin <- foreach(i = IDs, .combine = combine_list, + .packages = c("SIMplyBee")) %do% { + drones[as.character(crossPlan[[i]])] + } } } - validObject(ret) - return(ret) -} + # At this point, x is a Pop and dronesByVirgin are a list (so are they if they come if via drone packages) + if (oneColony) { + dronesByVirgin <- list(drones) + } + if (dronePackages) { + dronesByVirgin <- drones + } -#' @rdname setQueensYearOfBirth -#' @title Set the queen's year of birth -#' -#' @description Level 1 function that sets the queen's year of birth. -#' -#' @param x \code{\link[AlphaSimR]{Pop-class}} (one or more than one queen), -#' \code{\link[SIMplyBee]{Colony-class}} (one colony), or -#' \code{\link[SIMplyBee]{MultiColony-class}} (more colonies) -#' @param year integer, the year of the birth of the queen -#' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters -#' -#' @return \code{\link[AlphaSimR]{Pop-class}}, \code{\link[SIMplyBee]{Colony-class}}, or -#' \code{\link[SIMplyBee]{MultiColony-class}} with queens having the year of birth set -#' -#' @examples -#' founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) -#' SP <- SimParamBee$new(founderGenomes) -#' \dontshow{SP$nThreads = 1L} -#' basePop <- createVirginQueens(founderGenomes) -#' -#' drones <- createDrones(x = basePop[1], nInd = 1000) -#' droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson) -#' -#' # Create a Colony and a MultiColony class -#' colony <- createColony(x = basePop[2]) -#' colony <- cross(x = colony, drones = droneGroups[[1]]) -#' apiary <- createMultiColony(basePop[3:4], n = 2) -#' apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) -#' -#' # Example on Colony class -#' getQueenYearOfBirth(colony) -#' getQueenYearOfBirth(apiary) -#' -#' queen1 <- getQueen(colony) -#' queen1 <- setQueensYearOfBirth(queen1, year = 2022) -#' getQueenYearOfBirth(queen1) -#' -#' colony <- setQueensYearOfBirth(colony, year = 2022) -#' getQueenYearOfBirth(colony) -#' -#' apiary <- setQueensYearOfBirth(apiary, year = 2022) -#' getQueenYearOfBirth(apiary) -#' @export -setQueensYearOfBirth <- function(x, year, simParamBee = NULL) { - if (is.null(simParamBee)) { - simParamBee <- get(x = "SP", envir = .GlobalEnv) + names(dronesByVirgin) <- IDs + nDronesByVirgin <- sapply(dronesByVirgin, FUN = function(x) nInd(x)) + + + #if (all(nDronesByVirgin > 0)) { #There was a mistake here - if the message is warning, this still needs to happen + if (!all(sapply(dronesByVirgin, + FUN = function(x) all(isDrone(x, simParamBee = simParamBee))))) { + stop("Individuals in drones must be drones!") } - if (isPop(x)) { - if (any(!(isVirginQueen(x, simParamBee = simParamBee) | isQueen(x, simParamBee = simParamBee)))) { - stop("Individuals in x must be virgin queens or queens!") + + if (nInd(x) != length(dronesByVirgin)) { + stop("Number of virgin queens does not match the length of the assigned drones!") + } + + for (id in IDs) { + simParamBee$changeCaste(id = id, caste = "queen") + } + + for (id in as.vector(Reduce("c", sapply(dronesByVirgin, FUN = function(x) getId(x))))) { + simParamBee$changeCaste(id = id, caste = "fathers") + } + + # All of the input has been transformed to a Pop + crossVirginQueen <- function(virginQueen, virginQueenDrones, simParamBee = NULL) { + virginQueen@misc$fathers[[1]] <- virginQueenDrones + virginQueen@misc[["nWorkers"]] <- 0 + virginQueen@misc[["nDrones"]] <- 0 + virginQueen@misc[["nHomBrood"]] <- 0 + + if (isCsdActive(simParamBee = simParamBee)) { #This does still not work it the CSD is turned on + val <- calcQueensPHomBrood(x = virginQueen, simParamBee = simParamBee) + } else { + val <- NA } - nInd <- nInd(x) - x <- setMisc(x = x, node = "yearOfBirth", value = year) - } else if (isColony(x)) { - if (isQueenPresent(x, simParamBee = simParamBee)) { - x@queen <- setMisc(x = x@queen, node = "yearOfBirth", value = year) + + virginQueen@misc[["pHomBrood"]] <- val + return(virginQueen) + } + + # Add drones in the queens father slot + + x <- foreach(i = 1:length(IDs), .combine = combine_list, .packages = "SIMplyBee") %dopar% { + crossVirginQueen(virginQueen = x[i], + virginQueenDrones = dronesByVirgin[[i]], + simParamBee = simParamBee) + } + + if (type == "Pop") { + if (length(x) == 1) { + ret <- x } else { - stop("Missing queen!") + ret <- mergePops(x) } - } else if (isMultiColony(x)) { - nCol <- nColonies(x) - for (colony in seq_len(nCol)) { - x[[colony]]@queen <- setMisc( - x = x[[colony]]@queen, node = "yearOfBirth", - value = year - ) + } else if (type == "Colony") { + ret <- reQueen(x = colony, queen = x[1], simParamBee = simParamBee) + ret <- removeVirginQueens(ret, simParamBee = simParamBee) + } else if (type == "MultiColony") { + if (length(IDs) > 1) { + x <- mergePops(x) } - } else { - stop("Argument x must be a Pop, Colony or MultiColony class object!") + ret <- reQueen(x = multicolony, queen = x, simParamBee = simParamBee) + ret <- removeCastePop(ret, caste = "virginQueens", simParamBee = simParamBee) } - return(x) + + validObject(ret) + return(ret) } + diff --git a/R/Functions_L2_Colony.R b/R/Functions_L2_Colony.R index 390e76f5..7320c3f6 100644 --- a/R/Functions_L2_Colony.R +++ b/R/Functions_L2_Colony.R @@ -8,6 +8,7 @@ #' #' @param x \code{\link[AlphaSimR]{Pop-class}}, one queen or virgin queen(s) #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters +#' @param id character, ID of the colony that is going to be created (used internally for parallel computing) #' #' @return new \code{\link[SIMplyBee]{Colony-class}} #' @@ -31,15 +32,19 @@ #' colony1 <- cross(colony1, drones = drones) #' colony1 #' @export -createColony <- function(x = NULL, simParamBee = NULL) { +createColony <- function(x = NULL, simParamBee = NULL, id = NULL) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } simParamBee$updateLastColonyId() + if (is.null(id)) { + id <- simParamBee$lastColonyId + } + if (is.null(x)) { colony <- new( Class = "Colony", - id = simParamBee$lastColonyId + id = id ) } else { if (!isPop(x)) { @@ -60,13 +65,13 @@ createColony <- function(x = NULL, simParamBee = NULL) { colony <- new( Class = "Colony", - id = simParamBee$lastColonyId, + id = id, queen = queen, location = c(0, 0), virginQueens = virginQueens ) } - colony <- resetEvents(colony) + colony <- resetEvents(colony, simParamBee = simParamBee) validObject(colony) return(colony) } @@ -107,7 +112,7 @@ createColony <- function(x = NULL, simParamBee = NULL) { #' # Create and cross Colony and MultiColony class #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) -#' apiary <- createMultiColony(basePop[3:4], n = 2) +#' apiary <- createMultiColony(basePop[3:4]) #' apiary <- cross(apiary, drones = droneGroups[2:3]) #' #' # Check queen and virgin queens IDs @@ -158,16 +163,21 @@ reQueen <- function(x, queen, removeVirginQueens = TRUE, simParamBee = NULL) { x <- removeVirginQueens(x, simParamBee = simParamBee) } } else { - x <- removeQueen(x, addVirginQueens = FALSE, simParamBee = simParamBee) + x <- removeQueen(x, simParamBee = simParamBee) x@virginQueens <- queen } } else if (isMultiColony(x)) { nCol <- nColonies(x) + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") + } if (nInd(queen) < nCol) { stop("Not enough queens provided!") } - for (colony in seq_len(nCol)) { - x[[colony]] <- reQueen( + + + x@colonies = foreach(colony = seq_len(nCol), .packages = c("SIMplyBee")) %dopar% { + reQueen( x = x[[colony]], queen = queen[colony], simParamBee = simParamBee @@ -180,6 +190,33 @@ reQueen <- function(x, queen, removeVirginQueens = TRUE, simParamBee = NULL) { return(x) } +#' @rdname addCastePop_internal +#' @title An internal function to add a population in a caste slot of the colony +#' +#' @description Helper function that returns a colony to allow parallelisation, +#' only for internal use. +#' +#' @param colony \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} +#' @param pop \code{\link[AlphaSimR]{Pop-class}} with one or many individual +#' @param caste character +#' @param new logical +#' +#' @return \code{\link[SIMplyBee]{Colony-class}} +#' @export +addCastePop_internal <- function(pop, colony, caste, new = FALSE) { + if (!is.null(pop)) { + if (caste == "queen" & nInd(pop) > 1) { + stop("Cannot add more than one queen!") + } + } + if (is.null(slot(colony, caste)) | new) { + slot(colony, caste) <- pop + } else { + slot(colony, caste) <- c(slot(colony, caste), pop) + } + return(colony) +} + #' @rdname addCastePop #' @title Add caste individuals to the colony #' @@ -197,10 +234,6 @@ reQueen <- function(x, queen, removeVirginQueens = TRUE, simParamBee = NULL) { #' a single value is provided, the same value will be used for all the colonies. #' @param new logical, should the number of individuals be added to the caste population #' anew or should we only top-up the existing number of individuals to \code{nInd} -#' @param exact logical, only relevant when adding workers - if the csd locus is turned -#' on and exact is \code{TRUE}, we add the exact specified number of viable workers -#' (heterozygous at the csd locus) -#' @param year numeric, only relevant when adding virgin queens - year of birth for virgin queens #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters #' @param ... additional arguments passed to \code{nInd} when this argument is a function #' @@ -222,7 +255,7 @@ reQueen <- function(x, queen, removeVirginQueens = TRUE, simParamBee = NULL) { #' # Create and cross Colony and MultiColony class #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) -#' apiary <- createMultiColony(basePop[4:5], n = 2) +#' apiary <- createMultiColony(basePop[4:5]) #' apiary <- cross(apiary, drones = droneGroups[3:4]) #' #' #Here we show an example for workers, but same holds for drones and virgin queens! @@ -249,7 +282,7 @@ reQueen <- function(x, queen, removeVirginQueens = TRUE, simParamBee = NULL) { #' # nVirginQueens/nWorkers/nDrones will vary between function calls when a function is used #' #' # Queen's counters -#' getMisc(getQueen(addWorkers(colony))) +#' getQueen(addWorkers(colony))@misc #' #' # Add individuals to a MultiColony object #' apiary <- addWorkers(apiary) @@ -258,7 +291,7 @@ reQueen <- function(x, queen, removeVirginQueens = TRUE, simParamBee = NULL) { #' #' @export addCastePop <- function(x, caste = NULL, nInd = NULL, new = FALSE, - exact = FALSE, year = NULL, simParamBee = NULL, ...) { + simParamBee = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } @@ -293,16 +326,22 @@ addCastePop <- function(x, caste = NULL, nInd = NULL, new = FALSE, stop("nInd must be non-negative or NULL!") } } + if (length(nInd) > 1) { + warning("More than one value in the nInd argument, taking only the first value!") + nInd <- nInd[1] + } if (0 < nInd) { newInds <- createCastePop(x, nInd, - caste = caste, exact = exact, - year = year, simParamBee = simParamBee + caste = caste, + simParamBee = simParamBee ) if (caste == "workers") { homInds <- newInds$nHomBrood newInds <- newInds$workers x@queen@misc$nWorkers[[1]] <- x@queen@misc$nWorkers[[1]] + nInd(newInds) - x@queen@misc$nHomBrood[[1]] <- x@queen@misc$nHomBrood[[1]] + homInds + if (isCsdActive(simParamBee = simParamBee)) { + x@queen@misc$nHomBrood[[1]] <- x@queen@misc$nHomBrood[[1]] + homInds + } } if (caste == "drones") { x@queen@misc$nDrones[[1]] <- x@queen@misc$nDrones[[1]] + nInd(newInds) @@ -317,27 +356,48 @@ addCastePop <- function(x, caste = NULL, nInd = NULL, new = FALSE, } } else if (isMultiColony(x)) { nCol <- nColonies(x) - nNInd <- length(nInd) - if (nNInd > 1 && nNInd < nCol) { - stop("Too few values in the nInd argument!") + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") } - if (nNInd > 1 && nNInd > nCol) { - warning(paste0("Too many values in the nInd argument, taking only the first ", nCol, "values!")) - nInd <- nInd[1:nCol] + if (any(hasCollapsed(x))) { + stop(paste0("The colony ", getId(x), " collapsed, hence you can not add individuals (from the queen) to it!")) } - for (colony in seq_len(nCol)) { - if (is.null(nInd)) { - nIndColony <- NULL + + newInds <- createCastePop(x, nInd, + caste = caste, + simParamBee = simParamBee, + returnSP = FALSE, ...) + + + if (caste == "workers") { + homInds <- lapply(newInds, function(x) { + if (is.null(x)) return(NULL) + x[['nHomBrood']] + }) + newInds <- lapply(newInds, function(x) { + if (is.null(x)) return(NULL) + x[["workers"]] + }) + } + nInds <- lapply(newInds, function(x) { + if (is.null(x)) return(NULL) + nInd(x) + }) + + x@colonies <- foreach(colony = seq_len(nCol), .packages = c("SIMplyBee")) %dopar% { + if (!is.null(nInds[[colony]])) { + if (caste == "workers") { + x[[colony]]@queen@misc$nWorkers[[1]] <- x[[colony]]@queen@misc$nWorkers[[1]] + nInds[[colony]] + x[[colony]]@queen@misc$nHomBrood[[1]] <- x[[colony]]@queen@misc$nHomBrood[[1]] + ifelse(is.null(homInds[[colony]]), 0, homInds[[colony]]) + } else if (caste == "drones") { + x[[colony]]@queen@misc$nDrones[[1]] <- x[[colony]]@queen@misc$nDrones[[1]] + nInds[[colony]] + } + addCastePop_internal(colony = x[[colony]], pop = newInds[[colony]], caste = caste, new = new) } else { - nIndColony <- ifelse(nNInd == 1, nInd, nInd[colony]) + x[[colony]] } - x[[colony]] <- addCastePop( - x = x[[colony]], caste = caste, - nInd = nIndColony, - new = new, - exact = exact, simParamBee = simParamBee, ... - ) } + } else { stop("Argument x must be a Colony or MultiColony class object!") } @@ -348,17 +408,18 @@ addCastePop <- function(x, caste = NULL, nInd = NULL, new = FALSE, #' @describeIn addCastePop Add workers to a colony #' @export addWorkers <- function(x, nInd = NULL, new = FALSE, - exact = FALSE, simParamBee = NULL, ...) { + simParamBee = NULL, ...) { ret <- addCastePop( x = x, caste = "workers", nInd = nInd, new = new, - exact = exact, simParamBee = simParamBee, ... + simParamBee = simParamBee, ... ) return(ret) } #' @describeIn addCastePop Add drones to a colony #' @export -addDrones <- function(x, nInd = NULL, new = FALSE, simParamBee = NULL, ...) { +addDrones <- function(x, nInd = NULL, new = FALSE, + simParamBee = NULL, ...) { ret <- addCastePop( x = x, caste = "drones", nInd = nInd, new = new, simParamBee = simParamBee, ... @@ -369,10 +430,10 @@ addDrones <- function(x, nInd = NULL, new = FALSE, simParamBee = NULL, ...) { #' @describeIn addCastePop Add virgin queens to a colony #' @export addVirginQueens <- function(x, nInd = NULL, new = FALSE, - year = NULL, simParamBee = NULL, ...) { + simParamBee = NULL, ...) { ret <- addCastePop( x = x, caste = "virginQueens", nInd = nInd, new = new, - year = year, simParamBee = simParamBee, ... + simParamBee = simParamBee, ... ) return(ret) } @@ -398,9 +459,6 @@ addVirginQueens <- function(x, nInd = NULL, new = FALSE, #' @param new logical, should the number of workers and drones be added anew or #' should we only top-up the existing number of workers and drones to #' \code{nWorkers} and \code{nDrones} (see details) -#' @param exact logical, if the csd locus is turned on and exact is \code{TRUE}, -#' create the exact specified number of only viable workers (heterozygous on -#' the csd locus) #' @param resetEvents logical, call \code{\link[SIMplyBee]{resetEvents}} as part of the #' build up #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters @@ -440,7 +498,7 @@ addVirginQueens <- function(x, nInd = NULL, new = FALSE, #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) #' isProductive(colony) -#' apiary <- createMultiColony(basePop[3:4], n = 2) +#' apiary <- createMultiColony(basePop[3:4]) #' apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) #' isProductive(apiary) #' @@ -472,11 +530,12 @@ addVirginQueens <- function(x, nInd = NULL, new = FALSE, #' nDrones(apiary) #' #' # Queen's counters -#' getMisc(getQueen(buildUp(colony))) +#' getQueen(buildUp(colony))@misc +#' #' @export buildUp <- function(x, nWorkers = NULL, nDrones = NULL, - new = TRUE, exact = FALSE, resetEvents = FALSE, - simParamBee = NULL, ...) { + new = TRUE, resetEvents = FALSE, + simParamBee = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } @@ -494,7 +553,7 @@ buildUp <- function(x, nWorkers = NULL, nDrones = NULL, if (isColony(x)) { if (is.function(nWorkers)) { - nWorkers <- nWorkers(colony = x,...) + nWorkers <- nWorkers(x = x,...) } if (hasCollapsed(x)) { stop(paste0("The colony ", getId(x), " collapsed, hence you can not build it up!")) @@ -512,7 +571,7 @@ buildUp <- function(x, nWorkers = NULL, nDrones = NULL, if (0 < n) { x <- addWorkers( x = x, nInd = n, new = new, - exact = exact, simParamBee = simParamBee) + simParamBee = simParamBee) } else if (n < 0) { x@workers <- getWorkers(x, nInd = nWorkers, simParamBee = simParamBee) } @@ -543,7 +602,14 @@ buildUp <- function(x, nWorkers = NULL, nDrones = NULL, } x@production <- TRUE } else if (isMultiColony(x)) { + + if (any(hasCollapsed(x))) { + stop(paste0("Some colonies are collapsed, hence you can not build it up!")) + } nCol <- nColonies(x) + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") + } nNWorkers <- length(nWorkers) nNDrones <- length(nDrones) if (nNWorkers > 1 && nNWorkers < nCol) { @@ -560,27 +626,32 @@ buildUp <- function(x, nWorkers = NULL, nDrones = NULL, warning(paste0("Too many values in the nDrones argument, taking only the first ", nCol, "values!")) nNDrones <- nNDrones[1:nCol] } - for (colony in seq_len(nCol)) { - if (is.null(nWorkers)) { - nWorkersColony <- NULL - } else { - nWorkersColony <- ifelse(nNWorkers == 1, nWorkers, nWorkers[colony]) - } - if (is.null(nDrones)) { - nDronesColony <- NULL - } else { - nDronesColony <- ifelse(nNDrones == 1, nDrones, nDrones[colony]) - } - x[[colony]] <- buildUp( - x = x[[colony]], - nWorkers = nWorkersColony, - nDrones = nDronesColony, - new = new, - exact = exact, - resetEvents = resetEvents, - simParamBee = simParamBee, ... - ) + + if (is.function(nWorkers)) { + nWorkers <- nWorkers(x = x,...) + } + + if (new) { + n <- nWorkers + } else { + n <- nWorkers - nWorkers(x, simParamBee = simParamBee) + } + + if (sum(nWorkers) > 0) { + x <- addWorkers( + x = x, nInd = n, new = new, + simParamBee = simParamBee) + } + if (sum(nDrones) > 0) { + x <- addDrones( + x = x, nInd = n, new = new, + simParamBee = simParamBee) } + x <- setEvents(x, slot = "production", value = TRUE, simParamBee = simParamBee) + if (resetEvents) { + x <- resetEvents(x, simParamBee = simParamBee) + } + } else { stop("Argument x must be a Colony or MultiColony class object!") } @@ -589,7 +660,6 @@ buildUp <- function(x, nWorkers = NULL, nDrones = NULL, return(x) } - #' @rdname downsize #' @title Reduce number of workers and remove all drones and virgin queens from #' a Colony or MultiColony object @@ -627,7 +697,7 @@ buildUp <- function(x, nWorkers = NULL, nDrones = NULL, #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) #' colony <- buildUp(colony) -#' apiary <- createMultiColony(basePop[3:4], n = 2) +#' apiary <- createMultiColony(basePop[3:4]) #' apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) #' apiary <- buildUp(apiary) #' @@ -642,6 +712,7 @@ buildUp <- function(x, nWorkers = NULL, nDrones = NULL, #' apiary <- downsize(x = apiary, p = c(0.5, 0.1), new = TRUE, use = "rand") #' nWorkers(apiary); nDrones(apiary) #' @export +#' downsize <- function(x, p = NULL, use = "rand", new = FALSE, simParamBee = NULL, ...) { if (is.null(simParamBee)) { @@ -681,6 +752,19 @@ downsize <- function(x, p = NULL, use = "rand", new = FALSE, } else if (isMultiColony(x)) { nCol <- nColonies(x) nP <- length(p) + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") + } + + if (any(hasCollapsed(x))) { + stop("Some of the colonies have collapsed, hence you can not downsize them!") + } + if (is.null(p)) { + p <- simParamBee$downsizeP + } + if (is.function(p)) { + p <- p(x, ...) + } if (nP > 1 && nP < nCol) { stop("Too few values in the p argument!") } @@ -688,20 +772,20 @@ downsize <- function(x, p = NULL, use = "rand", new = FALSE, warning(paste0("Too many values in the p argument, taking only the first ", nCol, "values!")) p <- p[1:nCol] } - for (colony in seq_len(nCol)) { - if (is.null(p)) { - pColony <- NULL - } else { - pColony <- ifelse(nP == 1, p, p[colony]) - } - x[[colony]] <- downsize( - x = x[[colony]], - p = pColony, - use = use, - new = new, - simParamBee = simParamBee, ... - ) + if (new == TRUE) { + n <- round(nWorkers(x, simParamBee = simParamBee) * (1 - p)) + x <- addWorkers(x = x, nInd = n, new = TRUE, + simParamBee = simParamBee) + } else { + x <- removeWorkers(x = x, p = p, use = use, + simParamBee = simParamBee) + } + x <- removeDrones(x = x, p = 1, simParamBee = simParamBee) + x <- removeVirginQueens(x = x, p = 1, simParamBee = simParamBee) + for (colony in 1:nCol) { + x[[colony]]@production <- FALSE } + } else { stop("Argument x must be a Colony or MultiColony class object!") } @@ -710,6 +794,8 @@ downsize <- function(x, p = NULL, use = "rand", new = FALSE, return(x) } + + #' @rdname replaceCastePop #' @title Replace a proportion of caste individuals with new ones #' @@ -726,12 +812,6 @@ downsize <- function(x, p = NULL, use = "rand", new = FALSE, #' a single value is provided, the same value will be applied to all the colonies #' @param use character, all the options provided by \code{\link[AlphaSimR]{selectInd}} - #' guides selection of caste individuals that stay when \code{p < 1} -#' @param exact logical, only relevant when adding workers - if the csd locus is turned -#' on and exact is \code{TRUE}, we replace the exact specified number of viable workers -#' (heterozygous at the csd locus). You probably want this set to TRUE since you want to -#' replace with the same number of workers. -#' @param year numeric, only relevant when replacing virgin queens, -#' year of birth for virgin queens #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters #' #' @return \code{\link[SIMplyBee]{Colony-class}} or or \code{\link[SIMplyBee]{MultiColony-class}} with @@ -749,7 +829,7 @@ downsize <- function(x, p = NULL, use = "rand", new = FALSE, #' # Create and cross Colony and MultiColony class #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) -#' apiary <- createMultiColony(basePop[4:5], n = 2) +#' apiary <- createMultiColony(basePop[4:5]) #' apiary <- cross(apiary, drones = droneGroups[3:4]) #' #' # Add individuals @@ -768,8 +848,8 @@ downsize <- function(x, p = NULL, use = "rand", new = FALSE, #' apiary <- replaceWorkers(apiary, p = 0.5) #' getCasteId(apiary, caste="workers") #' @export -replaceCastePop <- function(x, caste = NULL, p = 1, use = "rand", exact = TRUE, - year = NULL, simParamBee = NULL) { +replaceCastePop <- function(x, caste = NULL, p = 1, use = "rand", + simParamBee = NULL) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } @@ -781,72 +861,50 @@ replaceCastePop <- function(x, caste = NULL, p = 1, use = "rand", exact = TRUE, } else if (any(p < 0)) { stop("p must not be less than 0!") } - if (isColony(x)) { - if (hasCollapsed(x)) { - stop(paste0("The colony ", getId(x), " collapsed, hence you can not replace individuals in it!")) + if (isColony(x) | isMultiColony(x)) { + nP <- length(p) + if (isColony(x)) { + nCol <- 1 + } else if (isMultiColony(x)) { + nCol <- nColonies(x) } - if (!isQueenPresent(x, simParamBee = simParamBee)) { - stop("Missing queen!") + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") } - if (length(p) > 1) { - warning("More than one value in the p argument, taking only the first value!") - p <- p[1] + if (any(hasCollapsed(x))) { + stop(paste0("The colony or some of the colonies have collapsed, hence you can not replace individuals in it!")) + } + if (any(!isQueenPresent(x, simParamBee = simParamBee))) { + stop("Missing queen in at least one colony!") + } + if (nP > 1 && nP < nCol) { + stop("Too few values in the p argument!") + } + if (length(p) > nCol) { + warning(paste0("More than one value in the p argument, taking only the first ", nCol, " values!")) + p <- p[nCol] } nInd <- nCaste(x, caste, simParamBee = simParamBee) - if (nInd > 0) { + if (any(nInd > 0)) { nIndReplaced <- round(nInd * p) - if (nIndReplaced < nInd) { - nIndStay <- nInd - nIndReplaced - if (nIndReplaced > 0) { - tmp <- createCastePop(x, - caste = caste, - nInd = nIndReplaced, exact = exact, - year = year, simParamBee = simParamBee - ) - if (caste == "workers") { - x@queen@misc$nWorkers[[1]] <- x@queen@misc$nWorkers[[1]] + nIndReplaced - x@queen@misc$nHomBrood[[1]] <- x@queen@misc$nHomBrood[[1]] + tmp$nHomBrood - tmp <- tmp$workers - } - if (caste == "drones") { - x@queen@misc$nDrones[[1]] <- x@queen@misc$nDrones[[1]] + nIndReplaced - } - - slot(x, caste) <- c( - selectInd(slot(x, caste), nInd = nIndStay, use = use, simParam = simParamBee), - tmp - ) - } + if (any(nIndReplaced < nInd)) { + + x <- removeCastePop(x, + caste = caste, + p = p, simParamBee = simParamBee) + nIndAdd <- nInd - nCaste(x, caste, simParamBee = simParamBee) + x <- addCastePop(x, + caste = caste, + nInd = nIndAdd, + simParamBee = simParamBee + ) } else { x <- addCastePop( x = x, caste = caste, nInd = nIndReplaced, new = TRUE, - year = year, simParamBee = simParamBee + simParamBee = simParamBee ) } } - } else if (isMultiColony(x)) { - nCol <- nColonies(x) - nP <- length(p) - if (nP > 1 && nP < nCol) { - stop("Too few values in the p argument!") - } - if (nP > 1 && nP > nCol) { - warning(paste0("Too many values in the p argument, taking only the first ", nCol, "values!")) - p <- p[1:nCol] - } - for (colony in seq_len(nCol)) { - if (is.null(p)) { - pColony <- NULL - } else { - pColony <- ifelse(nP == 1, p, p[colony]) - } - x[[colony]] <- replaceCastePop( - x = x[[colony]], caste = caste, - p = pColony, - use = use, year = year, - simParamBee = simParamBee - ) - } } else { stop("Argument x must be a Colony or MultiColony class object!") } @@ -856,10 +914,10 @@ replaceCastePop <- function(x, caste = NULL, p = 1, use = "rand", exact = TRUE, #' @describeIn replaceCastePop Replaces some workers in a colony #' @export -replaceWorkers <- function(x, p = 1, use = "rand", exact = TRUE, simParamBee = NULL) { +replaceWorkers <- function(x, p = 1, use = "rand", simParamBee = NULL) { ret <- replaceCastePop( x = x, caste = "workers", p = p, - use = use, exact = exact, + use = use, simParamBee = simParamBee ) return(ret) @@ -885,6 +943,7 @@ replaceVirginQueens <- function(x, p = 1, use = "rand", simParamBee = NULL) { return(ret) } + #' @rdname removeCastePop #' @title Remove a proportion of caste individuals from a colony #' @@ -898,13 +957,7 @@ replaceVirginQueens <- function(x, p = 1, use = "rand", simParamBee = NULL) { #' a single value is provided, the same value will be applied to all the colonies #' @param use character, all the options provided by \code{\link[AlphaSimR]{selectInd}} - #' guides selection of virgins queens that will stay when \code{p < 1} -#' @param addVirginQueens logical, whether virgin queens should be added; only -#' used when removing the queen from the colony -#' @param nVirginQueens integer, the number of virgin queens to be created in the -#' colony; only used when removing the queen from the colony. If \code{0}, no virgin -#' queens are added; If \code{NULL}, the value from \code{simParamBee$nVirginQueens} -#' is used -#' @param year numeric, only relevant when adding virgin queens - year of birth for virgin queens + #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters #' #' @return \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} without virgin queens @@ -922,7 +975,7 @@ replaceVirginQueens <- function(x, p = 1, use = "rand", simParamBee = NULL) { #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) #' colony <- buildUp(colony) -#' apiary <- createMultiColony(basePop[4:5], n = 2) +#' apiary <- createMultiColony(basePop[4:5]) #' apiary <- cross(apiary, drones = droneGroups[3:4]) #' apiary <- buildUp(apiary) #' @@ -943,8 +996,7 @@ replaceVirginQueens <- function(x, p = 1, use = "rand", simParamBee = NULL) { #' nWorkers(removeWorkers(apiary, p = c(0.1, 0.5))) #' @export removeCastePop <- function(x, caste = NULL, p = 1, use = "rand", - addVirginQueens = FALSE, nVirginQueens = NULL, - year = NULL, simParamBee = NULL) { + simParamBee = NULL) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } @@ -964,14 +1016,6 @@ removeCastePop <- function(x, caste = NULL, p = 1, use = "rand", warning("More than one value in the p argument, taking only the first value!") p <- p[1] } - if (caste == "queen") { - if (addVirginQueens) { - if (is.null(nVirginQueens)) { - nVirginQueens <- simParamBee$nVirginQueens - } - x <- addVirginQueens(x, nInd = nVirginQueens, year = year, simParamBee = simParamBee) - } - } if (p == 1) { slot(x, caste) <- NULL } else { @@ -990,6 +1034,9 @@ removeCastePop <- function(x, caste = NULL, p = 1, use = "rand", } else if (isMultiColony(x)) { nCol <- nColonies(x) nP <- length(p) + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") + } if (nP > 1 && nP < nCol) { stop("Too few values in the p argument!") } @@ -997,19 +1044,21 @@ removeCastePop <- function(x, caste = NULL, p = 1, use = "rand", warning(paste0("Too many values in the p argument, taking only the first ", nCol, "values!")) p <- p[1:nCol] } - for (colony in seq_len(nCol)) { + + x@colonies <- foreach(colony = seq_len(nCol), .packages = c("SIMplyBee")) %dopar% { if (is.null(p)) { pColony <- NULL } else { pColony <- ifelse(nP == 1, p, p[colony]) } - x[[colony]] <- removeCastePop( + removeCastePop( x = x[[colony]], caste = caste, p = pColony, use = use, simParamBee = simParamBee ) } + } else { stop("Argument x must be a Colony or MultiColony class object!") } @@ -1019,12 +1068,13 @@ removeCastePop <- function(x, caste = NULL, p = 1, use = "rand", #' @describeIn removeCastePop Remove queen from a colony #' @export -removeQueen <- function(x, addVirginQueens = FALSE, nVirginQueens = NULL, year = NULL, simParamBee = NULL) { - ret <- removeCastePop(x = x, caste = "queen", p = 1, addVirginQueens = addVirginQueens, - nVirginQueens = nVirginQueens, year = year, simParamBee = simParamBee) +#' +removeQueen <- function(x, simParamBee = NULL) { + ret <- removeCastePop(x = x, caste = "queen", p = 1, simParamBee = simParamBee) return(ret) } + #' @describeIn removeCastePop Remove workers from a colony #' @export removeWorkers <- function(x, p = 1, use = "rand", simParamBee = NULL) { @@ -1032,6 +1082,8 @@ removeWorkers <- function(x, p = 1, use = "rand", simParamBee = NULL) { return(ret) } + + #' @describeIn removeCastePop Remove workers from a colony #' @export removeDrones <- function(x, p = 1, use = "rand", simParamBee = NULL) { @@ -1039,6 +1091,8 @@ removeDrones <- function(x, p = 1, use = "rand", simParamBee = NULL) { return(ret) } + + #' @describeIn removeCastePop Remove virgin queens from a colony #' @export removeVirginQueens <- function(x, p = 1, use = "rand", simParamBee = NULL) { @@ -1059,6 +1113,7 @@ removeVirginQueens <- function(x, p = 1, use = "rand", simParamBee = NULL) { #' up a new colony, which the default of \code{NULL} caters for; otherwise, a #' collapsed colony should be left collapsed forever, unless you force #' resetting this event with \code{collapse = TRUE}) +#' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters #' #' @return \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} with #' events reset @@ -1075,7 +1130,7 @@ removeVirginQueens <- function(x, p = 1, use = "rand", simParamBee = NULL) { #' # Create and cross Colony and MultiColony class #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) -#' apiary <- createMultiColony(basePop[4:5], n = 2) +#' apiary <- createMultiColony(basePop[4:5]) #' apiary <- cross(apiary, drones = droneGroups[3:4]) #' #' # Build-up - this sets Productive to TRUE @@ -1125,7 +1180,10 @@ removeVirginQueens <- function(x, p = 1, use = "rand", simParamBee = NULL) { #' hasSplit(remnants[[1]]) #' resetEvents(remnants)[[1]] #' @export -resetEvents <- function(x, collapse = NULL) { +resetEvents <- function(x, collapse = NULL, simParamBee = NULL) { + if (is.null(simParamBee)) { + simParamBee <- get(x = "SP", envir = .GlobalEnv) + } if (isColony(x)) { x@swarm <- FALSE x@split <- FALSE @@ -1142,10 +1200,15 @@ resetEvents <- function(x, collapse = NULL) { validObject(x) } else if (isMultiColony(x)) { nCol <- nColonies(x) - for (colony in seq_len(nCol)) { - x[[colony]] <- resetEvents( + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") + } + + x@colonies <- foreach(colony = seq_len(nCol), .packages = c("SIMplyBee")) %dopar% { + resetEvents( x = x[[colony]], - collapse = collapse + collapse = collapse, + simParamBee = simParamBee ) } validObject(x) @@ -1166,6 +1229,8 @@ resetEvents <- function(x, collapse = NULL) { #' #' @return \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} with the collapse #' event set to \code{TRUE} +#' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters +#' #' #' @details You should use this function in an edge-case when you #' want to indicate that the colony has collapsed, but you still want to @@ -1178,13 +1243,13 @@ resetEvents <- function(x, collapse = NULL) { #' SP <- SimParamBee$new(founderGenomes) #' \dontshow{SP$nThreads = 1L} #' basePop <- createVirginQueens(founderGenomes) -#' drones <- createDrones(basePop[1], n = 1000) +#' drones <- createDrones(basePop[1], nInd = 1000) #' droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10) #' #' # Create Colony and MultiColony class #' colony <- createColony(x = basePop[1]) #' colony <- cross(colony, drones = droneGroups[[1]]) -#' apiary <- createMultiColony(x = basePop[2:10], n = 9) +#' apiary <- createMultiColony(x = basePop[2:10]) #' apiary <- cross(apiary, drones = droneGroups[2:10]) #' #' # Collapse @@ -1200,14 +1265,22 @@ resetEvents <- function(x, collapse = NULL) { #' apiaryLeft <- tmp$remnant #' hasCollapsed(apiaryLeft) #' @export -collapse <- function(x) { +collapse <- function(x, simParamBee = NULL) { + if (is.null(simParamBee)) { + simParamBee <- get(x = "SP", envir = .GlobalEnv) + } if (isColony(x)) { x@collapse <- TRUE x@production <- FALSE } else if (isMultiColony(x)) { nCol <- nColonies(x) - for (colony in seq_len(nCol)) { - x[[colony]] <- collapse(x = x[[colony]]) + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") + } + + x@colonies <- foreach(colony = seq_len(nCol), .packages = c("SIMplyBee")) %dopar% { + collapse(x = x[[colony]], + simParamBee = simParamBee) } } else { stop("Argument x must be a Colony or MultiColony class object!") @@ -1223,7 +1296,11 @@ collapse <- function(x) { #' an event where the queen #' leaves with a proportion of workers to create a new colony (the swarm). The #' remnant colony retains the other proportion of workers and all drones, and -#' the workers raise virgin queens, of which only one prevails. Location of +#' the workers raise virgin queens, of which only one prevails. The function +#' will create either 10 or \code{SimParamBee$nVirginQueens} virgin queens, +#' whichever is higher, and select one at random. In case of high inbreeding, +#' it could be that none of the virgin queens are viable. In that case, you might +#' want to increase \code{SimParamBee$nVirginQueens} or discard the colony. Location of #' the swarm is the same as for the remnant or sampled as deviation from the #' remnant. #' @@ -1233,11 +1310,6 @@ collapse <- function(x) { #' If input is \code{\link[SIMplyBee]{MultiColony-class}}, #' the input could also be a vector of the same length as the number of colonies. If #' a single value is provided, the same value will be applied to all the colonies -#' @param year numeric, year of birth for virgin queens -#' @param nVirginQueens integer, the number of virgin queens to be created in the -#' colony; of these one is randomly selected as the new virgin queen of the -#' remnant colony. If \code{NULL}, the value from \code{simParamBee$nVirginQueens} -#' is used #' @param sampleLocation logical, sample location of the swarm by taking #' the current colony location and adding deviates to each coordinate using #' \code{\link[SIMplyBee]{rcircle}} @@ -1255,6 +1327,8 @@ collapse <- function(x) { #' @examples #' founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 50) #' SP <- SimParamBee$new(founderGenomes) +#' SP$setTrackPed(TRUE) +#' SP$setTrackRec(TRUE) #' \dontshow{SP$nThreads = 1L} #' basePop <- createVirginQueens(founderGenomes) #' drones <- createDrones(basePop[1], n = 1000) @@ -1264,7 +1338,7 @@ collapse <- function(x) { #' colony <- createColony(x = basePop[2]) #' colony <- cross(colony, drones = droneGroups[[1]]) #' (colony <- buildUp(colony, nWorkers = 100)) -#' apiary <- createMultiColony(basePop[3:8], n = 6) +#' apiary <- createMultiColony(basePop[3:8]) #' apiary <- cross(apiary, drones = droneGroups[2:7]) #' apiary <- buildUp(apiary, nWorkers = 100) #' @@ -1288,150 +1362,172 @@ collapse <- function(x) { #' # Swarm only the pulled colonies #' (swarm(tmp$pulled, p = 0.6)) #' @export -swarm <- function(x, p = NULL, year = NULL, nVirginQueens = NULL, +swarm <- function(x, p = NULL, sampleLocation = TRUE, radius = NULL, simParamBee = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (isMultiColony(x)) { + parallel <- TRUE + } if (is.null(p)) { p <- simParamBee$swarmP } if (is.null(radius)) { radius <- simParamBee$swarmRadius } - if (is.null(nVirginQueens)) { - nVirginQueens <- simParamBee$nVirginQueens - } - if (isColony(x)) { - if (hasCollapsed(x)) { - stop(paste0("The colony ", getId(x), " collapsed, hence it can not swarm!")) + if (isColony(x) | isMultiColony(x)) { + if (isColony(x)) { + nCol <- 1 + } else if (isMultiColony(x)) { + nCol <- nColonies(x) } - if (!isQueenPresent(x, simParamBee = simParamBee)) { - stop("No queen present in the colony!") + nP <- length(p) + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") } - if (!isWorkersPresent(x, simParamBee = simParamBee)) { - stop("No workers present in the colony!") + + if (any(hasCollapsed(x))) { + stop(paste0("One of the collonies is collapsed, hence you can not split it!")) + } + if (any(!isQueenPresent(x, simParamBee = simParamBee))) { + stop("No queen present in one of the colonies!") + } + if (any(!isWorkersPresent(x, simParamBee = simParamBee))) { + stop("No workers present in one of the colonies!") } if (is.function(p)) { p <- p(x, ...) } else { - if (p < 0 | 1 < p) { + if (any(p < 0) | any(1 < p)) { stop("p must be between 0 and 1 (inclusive)!") } - if (length(p) > 1) { + if (length(p) > nCol) { warning("More than one value in the p argument, taking only the first value!") - p <- p[1] + p <- p[nCol] + } + if (nP > 1 && nP < nCol) { + stop("Too few values in the p argument!") } - } - if (is.function(nVirginQueens)) { - nVirginQueens <- nVirginQueens(x, ...) } nWorkers <- nWorkers(x, simParamBee = simParamBee) nWorkersSwarm <- round(nWorkers * p) # TODO: Add use="something" to select pWorkers that swarm # https://github.com/HighlanderLab/SIMplyBee/issues/160 - tmp <- pullWorkers(x = x, nInd = nWorkersSwarm, simParamBee = simParamBee) + + tmpVirginQueens <- createCastePop( + x = x, nInd = max(10, simParamBee$nVirginQueens), + caste = "virginQueens", + simParamBee = simParamBee + ) + + if (isColony(x)) { + homCol = nInd(tmpVirginQueens) == 0 + } else if (isMultiColony(x)) { + homCol = lapply(tmpVirginQueens, nInd) == 0 + } + + if (sum(homCol) > 0) { + if (isColony(x)) { + stop("Colony too inbred to produce any virgin queens!") + } else if (isMultiColony(x)) { + warning(paste0(sum(homCol), " colonies produced 0 virgin queens due to high colony homozygosity, removing these colonies!")) + tmpVirginQueens <- tmpVirginQueens[!homCol] + x = x[!homCol] + location = location[!homCol] + nWorkersSwarm = nWorkersSwarm[!homCol] + nCol = nColonies(x) + } + } + + tmp <- pullCastePop(x = x, caste = "workers", + nInd = nWorkersSwarm, simParamBee = simParamBee) + remnantColony <- tmp$remnant + remnantColony <- removeQueen(remnantColony, simParamBee = simParamBee) + if (isColony(x)) { + remnantColony <- reQueen(remnantColony, + queen = selectInd(tmpVirginQueens, nInd = 1, use = "rand", simParam = simParamBee), + simParamBee = simParamBee) + } else { + tmpVirginQueens <- lapply(tmpVirginQueens, FUN = function(x) selectInd(x, nInd = 1, use = "rand", simParam = simParamBee)) + remnantColony <- reQueen(remnantColony, + queen = mergePops(tmpVirginQueens), + simParamBee = simParamBee) + } currentLocation <- getLocation(x) + if (sampleLocation) { - newLocation <- c(currentLocation + rcircle(radius = radius)) + newLocation <- lapply(1:nCol, function(x) currentLocation[[x]] + rcircle(n = nCol, radius = radius)[x,]) # c() to convert row-matrix to a numeric vector } else { newLocation <- currentLocation } - swarmColony <- createColony(x = x@queen, simParamBee = simParamBee) - # It's not re-queening, but the function also sets the colony id - swarmColony@workers <- tmp$pulled - swarmColony <- setLocation(x = swarmColony, location = newLocation) + if (isColony(x)) { + swarmColony <- createColony(x = x@queen, simParamBee = simParamBee) + # It's not re-queening, but the function also sets the colony id - tmpVirginQueen <- createVirginQueens( - x = x, nInd = nVirginQueens, - year = year, - simParamBee = simParamBee - ) - tmpVirginQueen <- selectInd(tmpVirginQueen, nInd = 1, use = "rand", simParam = simParamBee) + swarmColony@workers <- tmp$pulled + swarmColony <- setLocation(x = swarmColony, location = newLocation[[1]], simParamBee = simParamBee) - remnantColony <- createColony(x = tmpVirginQueen, simParamBee = simParamBee) - remnantColony@workers <- getWorkers(tmp$remnant, simParamBee = simParamBee) - remnantColony@drones <- getDrones(x, simParamBee = simParamBee) - # Workers raise virgin queens from eggs laid by the queen and one random - # virgin queen prevails, so we create just one - # Could consider that a non-random one prevails (say the more aggressive one), - # by creating many virgin queens and then picking the one with highest - # gv/pheno for competition or some other criteria (patri-lineage) + remnantColony <- setLocation(x = remnantColony, location = currentLocation, simParamBee = simParamBee) - remnantColony <- setLocation(x = remnantColony, location = currentLocation) + remnantColony@swarm <- TRUE + swarmColony@swarm <- TRUE - remnantColony@swarm <- TRUE - swarmColony@swarm <- TRUE - remnantColony@production <- FALSE - swarmColony@production <- FALSE + remnantColony@production <- FALSE + swarmColony@production <- FALSE + + ret <- list(swarm = swarmColony, remnant = remnantColony) + } else if (isMultiColony(x)) { + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") + } - ret <- list(swarm = swarmColony, remnant = remnantColony) - } else if (isMultiColony(x)) { - nCol <- nColonies(x) - nP <- length(p) - if (nP > 1 && nP < nCol) { - stop("Too few values in the p argument!") - } - if (nP > 1 && nP > nCol) { - warning(paste0("Too many values in the p argument, taking only the first ", nCol, "values!")) - p <- p[1:nCol] - } - if (nCol == 0) { - ret <- list( - swarm = createMultiColony(simParamBee = simParamBee), - remnant = createMultiColony(simParamBee = simParamBee) - ) - } else { ret <- list( - swarm = createMultiColony(n = nCol, simParamBee = simParamBee), - remnant = createMultiColony(n = nCol, simParamBee = simParamBee) + swarm = createMultiColony(x = getQueen(x, collapse = TRUE, simParamBee = simParamBee), + simParamBee = simParamBee), + remnant = remnantColony ) - for (colony in seq_len(nCol)) { - if (is.null(p)) { - pColony <- NULL - } else { - pColony <- ifelse(nP == 1, p, p[colony]) - } - tmp <- swarm(x[[colony]], - p = pColony, - year = year, - nVirginQueens = nVirginQueens, - sampleLocation = sampleLocation, - radius = radius, - simParamBee = simParamBee, ... - ) - ret$swarm[[colony]] <- tmp$swarm - ret$remnant[[colony]] <- tmp$remnant + + ret$swarm@colonies <- foreach(colony = seq_len(nCol), .packages = c("SIMplyBee")) %dopar% { + addCastePop_internal(colony = ret$swarm@colonies[[colony]], + pop = tmp$pulled[[colony]], caste = "workers") } + + + ret$remnant <- setEvents(ret$remnant, slot = "swarm", value = TRUE, simParamBee = simParamBee) + ret$swarm <- setEvents(ret$swarm, slot = "swarm", value = TRUE, simParamBee = simParamBee) + ret$swarm <- setEvents(ret$swarm, slot = "production", value = FALSE, simParamBee = simParamBee) + ret$remnant <- setEvents(ret$remnant, slot = "production", value = FALSE, simParamBee = simParamBee) } - } else { + } + else { stop("Argument x must be a Colony or MultiColony class object!") } - validObject(ret$swarmColony) validObject(ret$remnantColony) return(ret) } + + #' @rdname supersede #' @title Supersede #' #' @description Level 2 function that supersedes a Colony or MultiColony object - #' an event where the #' queen dies. The workers and drones stay unchanged, but workers raise virgin -#' queens, of which only one prevails. +#' queens, of which only one prevails. The function +#' will create either 10 or \code{SimParamBee$nVirginQueens} virgin queens, +#' whichever is higher, and select one at random In case of high inbreeding, +#' it could be that none of the virgin queens are viable.In that case, you might +#' want to increase \code{SimParamBee$nVirginQueens} or discard the colony. #' #' @param x \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} -#' @param year numeric, year of birth for virgin queens -#' @param nVirginQueens integer, the number of virgin queens to be created in the -#' colony; of these one is randomly selected as the new virgin queen of the -#' remnant colony. If \code{NULL}, the value from \code{simParamBee$nVirginQueens} -#' is used #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters #' @param ... additional arguments passed to \code{nVirginQueens} when this #' argument is a function @@ -1474,26 +1570,61 @@ swarm <- function(x, p = NULL, year = NULL, nVirginQueens = NULL, #' # Swarm only the pulled colonies #' (supersede(tmp$pulled)) #' @export -supersede <- function(x, year = NULL, nVirginQueens = NULL, simParamBee = NULL, ...) { +supersede <- function(x, simParamBee = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (isColony(x)) { + parallel <- FALSE + } else if (isMultiColony(x)) { + parallel <- TRUE + } if (is.null(nVirginQueens)) { nVirginQueens <- simParamBee$nVirginQueens } + + if (any(hasCollapsed(x))) { + stop(paste0("One of the collonies is collapsed, hence you can not split it!")) + } + if (any(!isQueenPresent(x, simParamBee = simParamBee))) { + stop("No queen present in one of the colonies!") + } + if (is.function(nVirginQueens)) { + nVirginQueens <- nVirginQueens(x, simParamBee = simParamBee, ...) + } + + # Do this because some colonies might not produce a viable virgin queen + tmpVirginQueens <- createCastePop( + x = x, nInd = max(10, simParamBee$nVirginQueens), + caste = "virginQueens", + simParamBee = simParamBee + ) + if (isColony(x)) { - if (hasCollapsed(x)) { - stop(paste0("The colony ", getId(x), " collapsed, hence it can not supresede!")) - } - if (!isQueenPresent(x, simParamBee = simParamBee)) { - stop("No queen present in the colony!") + homCol = nInd(tmpVirginQueens) == 0 + } else if (isMultiColony(x)) { + homCol = sapply(tmpVirginQueens, nInd) == 0 + } + + if (sum(homCol) > 0) { + if (isColony(x)) { + print("X is colony") + print(class(x)) + stop("Colony to inbred to produce any virgin queens!") + } else if (isMultiColony(x)) { + warning(paste0(sum(homCol), " colonies produced 0 virgin queens due to high colony homozygosity, removing these colonies!")) + tmpVirginQueens <- tmpVirginQueens[!homCol] + x = x[!homCol] + nCol = nColonies(x) } - if (is.function(nVirginQueens)) { - nVirginQueens <- nVirginQueens(x, ...) + } + + if (isColony(x)) { + if (!parallel) { + x <- addCastePop_internal(selectInd(tmpVirginQueens, nInd = 1, use = "rand", simParam = simParamBee), + colony = x, caste = "virginQueens") } - x <- removeQueen(x, addVirginQueens = TRUE, nVirginQueens = nVirginQueens, - year = year, simParamBee = simParamBee) - x@virginQueens <- selectInd(x@virginQueens, nInd = 1, use = "rand", simParam = simParamBee) + x <- removeQueen(x, simParamBee = simParamBee) # TODO: We could consider that a non-random virgin queen prevails (say the most # aggressive one), by creating many virgin queens and then picking the # one with highest pheno for competition or some other criteria @@ -1502,17 +1633,25 @@ supersede <- function(x, year = NULL, nVirginQueens = NULL, simParamBee = NULL, } else if (isMultiColony(x)) { nCol <- nColonies(x) if (nCol == 0) { - x <- createMultiColony(simParamBee = simParamBee) - } else { - for (colony in seq_len(nCol)) { - x[[colony]] <- supersede(x[[colony]], - year = year, - nVirginQueens = nVirginQueens, - simParamBee = simParamBee, ... - ) + stop("The Multicolony contains 0 colonies!") + } + tmpVirginQueens <- lapply(tmpVirginQueens, FUN = function(x) selectInd(x, nInd = 1, use = "rand", simParam = simParamBee)) + + combine_list <- function(a, b) { + if (length(a) == 1) { + c(list(a), list(b)) + } else { + c(a, list(b)) } } - } else { + + x@colonies <- foreach(colony = seq_len(nColonies(x)), .packages = c("SIMplyBee")) %dopar% { + addCastePop_internal(colony = removeQueen(x[[colony]], simParamBee = simParamBee), + pop = tmpVirginQueens[[colony]], caste = "virginQueens") + } + x = setEvents(x, slot = "supersedure", value = TRUE, simParamBee = simParamBee) + } + else { stop("Argument x must be a Colony or MultiColony class object!") } validObject(x) @@ -1526,8 +1665,9 @@ supersede <- function(x, year = NULL, nVirginQueens = NULL, simParamBee = NULL, #' into two new colonies to #' prevent swarming (in managed situation). The remnant colony retains the #' queen and a proportion of the workers and all drones. The split colony gets -#' the other part of the workers, which raise virgin queens, of which only one -#' prevails. Location of the split is the same as for the remnant. +#' the other part of the workers, but note that it is queenless, since the beekeepers +#' would normally requeen with a different queen. +#' Location of the split is the same as for the remnant. #' #' @param x \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} #' @param p numeric, proportion of workers that will go to the split colony; if @@ -1535,7 +1675,6 @@ supersede <- function(x, year = NULL, nVirginQueens = NULL, simParamBee = NULL, #' If input is \code{\link[SIMplyBee]{MultiColony-class}}, #' the input could also be a vector of the same length as the number of colonies. If #' a single value is provided, the same value will be applied to all the colonies -#' @param year numeric, year of birth for virgin queens #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters #' @param ... additional arguments passed to \code{p} when this argument is a #' function @@ -1580,109 +1719,160 @@ supersede <- function(x, year = NULL, nVirginQueens = NULL, simParamBee = NULL, #' # Split only the pulled colonies #' (split(tmp$pulled, p = 0.5)) #' @export -split <- function(x, p = NULL, year = NULL, simParamBee = NULL, ...) { +split <- function(x, p = NULL, simParamBee = NULL, ...) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } if (is.null(p)) { p <- simParamBee$splitP } - if (isColony(x)) { - if (hasCollapsed(x)) { - stop(paste0("The colony ", getId(x), " collapsed, hence you can not split it!")) + if (isMultiColony(x)) { + parallel <- TRUE + } + + if (isColony(x) | isMultiColony(x)) { + if (isColony(x)) { + nCol <- 1 + } else if (isMultiColony(x)) { + nCol <- nColonies(x) } - if (!isQueenPresent(x, simParamBee = simParamBee)) { - stop("No queen present in the colony!") + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") + } + nP <- length(p) + + location <- getLocation(x) + if (any(hasCollapsed(x))) { + stop(paste0("One of the collonies is collapsed, hence you can not split it!")) + } + if (any(!isQueenPresent(x, simParamBee = simParamBee))) { + stop("No queen present in one of the colonies!") } - if (!isWorkersPresent(x, simParamBee = simParamBee)) { - stop("No workers present in the colony!") + if (any(!isWorkersPresent(x, simParamBee = simParamBee))) { + stop("No workers present in one of the colonies!") } if (is.function(p)) { p <- p(x, ...) } else { - if (p < 0 | 1 < p) { + if (any(p < 0) | any(1 < p)) { stop("p must be between 0 and 1 (inclusive)!") } - if (length(p) > 1) { + if (length(p) > nCol) { warning("More than one value in the p argument, taking only the first value!") - p <- p[1] + p <- p[nCol] + } + if (nP > 1 && nP < nCol) { + stop("Too few values in the p argument!") } } + nWorkers <- nWorkers(x, simParamBee = simParamBee) nWorkersSplit <- round(nWorkers * p) # TODO: Split colony at random by default, but we could make it as a # function of some parameters # https://github.com/HighlanderLab/SIMplyBee/issues/179 - tmp <- pullWorkers(x = x, nInd = nWorkersSplit, simParamBee = simParamBee) + tmp <- pullCastePop(x = x, caste = "workers", nInd = nWorkersSplit, simParamBee = simParamBee) remnantColony <- tmp$remnant - tmpVirginQueens <- createVirginQueens( - x = x, nInd = 1, - year = year, - simParamBee = simParamBee - ) - # Workers raise virgin queens from eggs laid by the queen (assuming) that - # a frame of brood is also provided to the split and then one random virgin - # queen prevails, so we create just one - # TODO: Could consider that a non-random one prevails (say the most aggressive - # one), by creating many virgin queens and then picking the one with - # highest pheno for competition or some other criteria - # https://github.com/HighlanderLab/SIMplyBee/issues/239 - splitColony <- createColony(x = tmpVirginQueens, simParamBee = simParamBee) - splitColony@workers <- tmp$pulled - splitColony <- setLocation(x = splitColony, location = getLocation(splitColony)) + if (isColony(x)) { - remnantColony@split <- TRUE - splitColony@split <- TRUE + # Workers raise virgin queens from eggs laid by the queen (assuming) that + # a frame of brood is also provided to the split and then one random virgin + # queen prevails, so we create just one + # TODO: Could consider that a non-random one prevails (say the most aggressive + # one), by creating many virgin queens and then picking the one with + # highest pheno for competition or some other criteria + # https://github.com/HighlanderLab/SIMplyBee/issues/239 - remnantColony@production <- TRUE - splitColony@production <- FALSE + splitColony <- createColony(simParamBee = simParamBee) + splitColony <- setLocation(x = splitColony, location = location, simParamBee = simParamBee) + + splitColony@workers <- tmp$pulled + + remnantColony@split <- TRUE + splitColony@split <- TRUE + + remnantColony@production <- TRUE + splitColony@production <- FALSE + + ret <- list(split = splitColony, remnant = remnantColony) + } else if (isMultiColony(x)) { + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") + } - ret <- list(split = splitColony, remnant = remnantColony) - } else if (isMultiColony(x)) { - nCol <- nColonies(x) - nP <- length(p) - if (nP > 1 && nP < nCol) { - stop("Too few values in the p argument!") - } - if (nP > 1 && nP > nCol) { - warning(paste0("Too many values in the nInd argument, taking only the first ", nCol, "values!")) - p <- p[1:nCol] - } - if (nCol == 0) { - ret <- list( - split = createMultiColony(simParamBee = simParamBee), - remnant = createMultiColony(simParamBee = simParamBee) - ) - } else { ret <- list( - split = createMultiColony(n = nCol, simParamBee = simParamBee), - remnant = createMultiColony(n = nCol, simParamBee = simParamBee) + split = createMultiColony(n = nCol, + simParamBee = simParamBee, + populateColonies = TRUE), + remnant = remnantColony + ) - for (colony in seq_len(nCol)) { - if (is.null(p)) { - pColony <- NULL - } else { - pColony <- ifelse(nP == 1, p, p[colony]) - } - tmp <- split(x[[colony]], - p = pColony, - year = year, - simParamBee = simParamBee, ... - ) - ret$split[[colony]] <- tmp$split - ret$remnant[[colony]] <- tmp$remnant + ret$split <- setLocation(x = ret$split, location = location, simParamBee = simParamBee) + + ret$split@colonies <- foreach(colony = seq_len(nCol), .packages = c("SIMplyBee")) %dopar% { + addCastePop_internal(colony = ret$split@colonies[[colony]], + pop = tmp$pulled[[colony]], caste = "workers") } + ret$split <- setEvents(ret$split, slot = "split", value = TRUE, simParamBee = simParamBee) + ret$remnant <- setEvents(ret$remnant, slot = "split", value = TRUE, simParamBee = simParamBee) + ret$split <- setEvents(ret$split, slot = "production", value = FALSE, simParamBee = simParamBee) + ret$remnant <- setEvents(ret$remnant, slot = "production", value = TRUE, simParamBee = simParamBee) } - } else { + } + else { stop("Argument x must be a Colony or MultiColony class object!") } - validObject(ret$splitColony) - validObject(ret$remnantColony) + validObject(ret$split) + validObject(ret$remnant) return(ret) } +#' @rdname setEvents +#' @title Set colony events +#' +#' @description Helper Level 2 function that populates the events slot. Not interded +#' for external use, intended for internal use in parallel computing +#' +#' @param x \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} +#' @param slot character, which event to set +#' @param value logical, the value for the event +#' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters +#' +#' @return \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} with +#' events reset +#' +#' @examples +#' founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 50) +#' SP <- SimParamBee$new(founderGenomes) +#' \dontshow{SP$nThreads = 1L} +#' basePop <- createVirginQueens(founderGenomes) +#' +#' drones <- createDrones(x = basePop[1], nInd = 100) +#' droneGroups <- pullDroneGroupsFromDCA(drones, n = 5, nDrones = nFathersPoisson) +#' +#' # Create and cross Colony and MultiColony class +#' colony <- createColony(x = basePop[2]) +#' colony <- cross(colony, drones = droneGroups[[1]]) +#' apiary <- createMultiColony(basePop[4:5]) +#' SIMplyBee:::setEvents(apiary, slot = "swarm", value = c(TRUE, TRUE)) +setEvents <- function(x, slot, value, simParamBee = NULL) { + if (is.null(simParamBee)) { + simParamBee <- get(x = "SP", envir = .GlobalEnv) + } + if (isColony(x)) { + slot(x, slot) <- value + } + if (isMultiColony(x)) { + x@colonies <- foreach(colony = seq_len(nColonies(x)), .packages = c("SIMplyBee")) %dopar% { + setEvents(x[[colony]], slot, value, simParamBee = simParamBee) + } + } + return(x) +} + + #' @rdname combine #' @title Combine two colony objects #' @@ -1695,6 +1885,7 @@ split <- function(x, p = NULL, year = NULL, simParamBee = NULL, ...) { #' #' @param strong \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} #' @param weak \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} +#' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters #' #' @return a combined \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} #' @@ -1702,8 +1893,9 @@ split <- function(x, p = NULL, year = NULL, simParamBee = NULL, ...) { #' founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 50) #' SP <- SimParamBee$new(founderGenomes) #' \dontshow{SP$nThreads = 1L} +#' print(SP$nThreads) #' basePop <- createVirginQueens(founderGenomes) -#' drones <- createDrones(basePop[1], n = 1000) +#' drones <- createDrones(basePop[1], nInd = 1000) #' droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10) #' #' # Create weak and strong Colony and MultiColony class @@ -1732,26 +1924,33 @@ split <- function(x, p = NULL, year = NULL, simParamBee = NULL, ...) { #' #' nWorkers(apiary1); nWorkers(apiary2) #' nDrones(apiary1); nDrones(apiary2) -#' apiary1 <- combine(strong = apiary1, weak = apiary2) +#' apiary1 <- combine(strong = apiary1, weak = apiary2, simParamBee = SP) #' nWorkers(apiary1); nWorkers(apiary2) #' nDrones(apiary1); nDrones(apiary2) #' rm(apiary2) #' @export -combine <- function(strong, weak) { +combine <- function(strong, weak, simParamBee = NULL) { + if (is.null(simParamBee)) { + simParamBee <- get(x = "SP", envir = .GlobalEnv) + } + if (any(hasCollapsed(strong))) { + stop(paste0("Some of the strong colonies have collapsed, hence you can not combine it!")) + } + if (any(hasCollapsed(weak))) { + stop(paste0("Some of the weak colonies have collapsed, hence you can not combine it!")) + } if (isColony(strong) & isColony(weak)) { - if (hasCollapsed(strong)) { - stop(paste0("The colony ", getId(strong), " (strong) has collapsed, hence you can not combine it!")) - } - if (hasCollapsed(weak)) { - stop(paste0("The colony ", getId(weak), " (weak) has collapsed, hence you can not combine it!")) - } strong@workers <- c(strong@workers, weak@workers) strong@drones <- c(strong@drones, weak@drones) } else if (isMultiColony(strong) & isMultiColony(weak)) { + if (nColonies(weak) == nColonies(strong)) { nCol <- nColonies(weak) - for (colony in seq_len(nCol)) { - strong[[colony]] <- combine(strong = strong[[colony]], weak = weak[[colony]]) + + strong@colonies <- foreach(colony = seq_len(nCol), .packages = c("SIMplyBee")) %dopar% { + combine(strong = strong[[colony]], + weak = weak[[colony]], + simParamBee = simParamBee) } } else { stop("Weak and strong MultiColony objects must be of the same length!") @@ -1762,6 +1961,7 @@ combine <- function(strong, weak) { return(strong) } + #' @rdname setLocation #' @title Set colony location #' @@ -1774,6 +1974,7 @@ combine <- function(strong, weak) { #' \code{c(x1, y1)} (the same location set to all colonies), #' \code{list(c(x1, y1), c(x2, y2))}, or #' \code{data.frame(x = c(x1, x2), y = c(y1, y2))} +#' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters #' #' @return \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} with set #' location @@ -1812,7 +2013,10 @@ combine <- function(strong, weak) { #' apiary <- setLocation(apiary, location = locDF) #' getLocation(apiary) #' @export -setLocation <- function(x, location = c(0, 0)) { +setLocation <- function(x, location = c(0, 0), simParamBee = NULL) { + if (is.null(simParamBee)) { + simParamBee <- get(x = "SP", envir = .GlobalEnv) + } if (isColony(x)) { if (is.list(location)) { # is.list() captures also is.data.frame() stop("Argument location must be numeric, when x is a Colony class object!") @@ -1822,21 +2026,24 @@ setLocation <- function(x, location = c(0, 0)) { } x@location <- location } else if (isMultiColony(x)) { - n <- nColonies(x) + nCol <- nColonies(x) + if (nCol == 0) { + stop("The Multicolony contains 0 colonies!") + } if (!is.null(location)) { if (is.numeric(location)) { if (length(location) != 2) { stop("When argument location is a numeric, it must be of length 2!") } } else if (is.data.frame(location)) { - if (nrow(location) != n) { + if (nrow(location) != nCol) { stop("When argument location is a data.frame, it must have as many rows as the number of colonies!") } if (ncol(location) != 2) { stop("When argument location is a data.frame, it must have 2 columns!") } } else if (is.list(location)) { - if (length(location) != n) { + if (length(location) != nCol) { stop("When argument location is a list, it must be of length equal to the number of colonies!") } tmp <- sapply(X = location, FUN = length) @@ -1851,7 +2058,15 @@ setLocation <- function(x, location = c(0, 0)) { stop("Argument location must be numeric, list, or data.frame!") } } - for (colony in seq_len(n)) { + combine_list <- function(a, b) { + if (length(a) == 1) { + c(list(a), list(b)) + } else { + c(a, list(b)) + } + } + + tmp <- foreach(colony = seq_len(nCol), .combine = combine_list, .packages = c("SIMplyBee")) %dopar% { if (is.data.frame(location)) { loc <- location[colony, ] loc <- c(loc$x, loc$y) @@ -1860,9 +2075,18 @@ setLocation <- function(x, location = c(0, 0)) { } else { loc <- location } + if (!is.null(x[[colony]])) { x[[colony]]@location <- loc } + + x[[colony]] + } + + if (nCol == 1) { + x@colonies = list(tmp) + } else { + x@colonies = tmp } } else { stop("Argument x must be a Colony or MultiColony class object!") diff --git a/R/Functions_L3_Colonies.R b/R/Functions_L3_Colonies.R index b8dbc191..1c8d7ff3 100644 --- a/R/Functions_L3_Colonies.R +++ b/R/Functions_L3_Colonies.R @@ -13,6 +13,7 @@ #' given then \code{\link[SIMplyBee]{MultiColony-class}} is created with \code{n} #' \code{NULL}) individual colony - this is mostly useful for programming) #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters +#' @param populateColonies boolean, whether to create n empty Colony objects within with assigned ID #' #' @details When both \code{x} and \code{n} are \code{NULL}, then a #' \code{\link[SIMplyBee]{MultiColony-class}} with 0 colonies is created. @@ -40,7 +41,7 @@ #' #' # Create mated colonies by crossing #' apiary <- createMultiColony(x = basePop[1:2], n = 2) -#' drones <- createDrones(x = basePop[3], n = 30) +#' drones <- createDrones(x = basePop[3], nInd = 30) #' droneGroups <- pullDroneGroupsFromDCA(drones, n = 2, nDrones = 15) #' apiary <- cross(apiary, drones = droneGroups) #' apiary @@ -48,15 +49,26 @@ #' apiary[[2]] #' #' @export -createMultiColony <- function(x = NULL, n = NULL, simParamBee = NULL) { +createMultiColony <- function(x = NULL, n = NULL, simParamBee = NULL, populateColonies = FALSE) { if (is.null(simParamBee)) { simParamBee <- get(x = "SP", envir = .GlobalEnv) } + if (is.null(x)) { if (is.null(n)) { ret <- new(Class = "MultiColony") } else { ret <- new(Class = "MultiColony", colonies = vector(mode = "list", length = n)) + if (populateColonies) { + ids <- (simParamBee$lastColonyId+1):(simParamBee$lastColonyId + n) + + ret@colonies <- foreach(colony = seq_len(n), .packages = c("SIMplyBee")) %dopar% { + createColony(simParamBee = simParamBee, id = ids[colony]) + } + simParamBee$updateLastColonyId(n = n) + } else { + + } } } else { if (!isPop(x)) { @@ -72,10 +84,14 @@ createMultiColony <- function(x = NULL, n = NULL, simParamBee = NULL) { stop("Not enough individuals in the x to create n colonies!") } ret <- new(Class = "MultiColony", colonies = vector(mode = "list", length = n)) - for (colony in seq_len(n)) { - ret[[colony]] <- createColony(x = x[colony], simParamBee = simParamBee) + ids <- (simParamBee$lastColonyId+1):(simParamBee$lastColonyId + n) + + ret@colonies <- foreach(colony = seq_len(n), .packages = c("SIMplyBee")) %dopar% { + createColony(x = x[colony], simParamBee = simParamBee, id = ids[colony]) } + simParamBee$updateLastColonyId(n = n) } + validObject(ret) return(ret) } diff --git a/R/SIMplyBee.R b/R/SIMplyBee.R index ce9260f8..04ad88bb 100644 --- a/R/SIMplyBee.R +++ b/R/SIMplyBee.R @@ -7,6 +7,7 @@ #' @importFrom stats rnorm rbeta runif rpois na.omit #' @importFrom extraDistr rtpois #' @importFrom utils packageVersion +#' @importFrom foreach foreach "%dopar%" "%do%" # see https://r-pkgs.org/namespace.html on description what to import/depend/... #' @description diff --git a/man/SimParamBee.Rd b/man/SimParamBee.Rd index 664475f4..00f6e5f5 100644 --- a/man/SimParamBee.Rd +++ b/man/SimParamBee.Rd @@ -317,6 +317,10 @@ generate this object} \item \href{#method-SimParamBee-new}{\code{SimParamBee$new()}} \item \href{#method-SimParamBee-addToCaste}{\code{SimParamBee$addToCaste()}} \item \href{#method-SimParamBee-changeCaste}{\code{SimParamBee$changeCaste()}} +\item \href{#method-SimParamBee-addToBeePed}{\code{SimParamBee$addToBeePed()}} +\item \href{#method-SimParamBee-addToBeeRec}{\code{SimParamBee$addToBeeRec()}} +\item \href{#method-SimParamBee-updateCaste}{\code{SimParamBee$updateCaste()}} +\item \href{#method-SimParamBee-updateLastBeeId}{\code{SimParamBee$updateLastBeeId()}} \item \href{#method-SimParamBee-updateLastColonyId}{\code{SimParamBee$updateLastColonyId()}} \item \href{#method-SimParamBee-clone}{\code{SimParamBee$clone()}} } @@ -532,6 +536,99 @@ SP$caste } +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParamBee-addToBeePed}{}}} +\subsection{Method \code{addToBeePed()}}{ +For internal use only. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{SimParamBee$addToBeePed(nNewInd, id, mother, father, isDH)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{nNewInd}}{Number of newly created individuals} + +\item{\code{id}}{the name of each individual} + +\item{\code{mother}}{vector of mother iids} + +\item{\code{father}}{vector of father iids} + +\item{\code{isDH}}{indicator for DH lines} +} +\if{html}{\out{
}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParamBee-addToBeeRec}{}}} +\subsection{Method \code{addToBeeRec()}}{ +For internal use only. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{SimParamBee$addToBeeRec(nNewInd, id, mother, father, isDH, hist, ploidy)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{nNewInd}}{Number of newly created individuals} + +\item{\code{id}}{the name of each individual} + +\item{\code{mother}}{vector of mother iids} + +\item{\code{father}}{vector of father iids} + +\item{\code{isDH}}{indicator for DH lines} + +\item{\code{hist}}{new recombination history} + +\item{\code{ploidy}}{ploidy level} +} +\if{html}{\out{
}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParamBee-updateCaste}{}}} +\subsection{Method \code{updateCaste()}}{ +A function to update the caste + For internal use only. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{SimParamBee$updateCaste(caste)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{caste}}{vector, named vector of castes to be added} +} +\if{html}{\out{
}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParamBee-updateLastBeeId}{}}} +\subsection{Method \code{updateLastBeeId()}}{ +A function to update the last + ID everytime we create an individual + For internal use in SIMplyBee only. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{SimParamBee$updateLastBeeId(n = 1L)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{n}}{integer, how many individuals to add} + +\item{\code{lastId}}{integer, last colony ID assigned} +} +\if{html}{\out{
}} +} } \if{html}{\out{
}} \if{html}{\out{}} @@ -541,12 +638,14 @@ A function to update the colony last ID everytime we create a Colony-class with createColony. For internal use only. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{SimParamBee$updateLastColonyId()}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{SimParamBee$updateLastColonyId(n = 1)}\if{html}{\out{
}} } \subsection{Arguments}{ \if{html}{\out{
}} \describe{ +\item{\code{n}}{integer, how many colonies to add} + \item{\code{lastColonyId}}{integer, last colony ID assigned} } \if{html}{\out{
}} diff --git a/man/addCastePop.Rd b/man/addCastePop.Rd index 89e896e3..87e1f727 100644 --- a/man/addCastePop.Rd +++ b/man/addCastePop.Rd @@ -7,29 +7,13 @@ \alias{addVirginQueens} \title{Add caste individuals to the colony} \usage{ -addCastePop( - x, - caste = NULL, - nInd = NULL, - new = FALSE, - exact = FALSE, - year = NULL, - simParamBee = NULL, - ... -) - -addWorkers(x, nInd = NULL, new = FALSE, exact = FALSE, simParamBee = NULL, ...) +addCastePop(x, caste = NULL, nInd = NULL, new = FALSE, simParamBee = NULL, ...) + +addWorkers(x, nInd = NULL, new = FALSE, simParamBee = NULL, ...) addDrones(x, nInd = NULL, new = FALSE, simParamBee = NULL, ...) -addVirginQueens( - x, - nInd = NULL, - new = FALSE, - year = NULL, - simParamBee = NULL, - ... -) +addVirginQueens(x, nInd = NULL, new = FALSE, simParamBee = NULL, ...) } \arguments{ \item{x}{\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}} @@ -45,12 +29,6 @@ a single value is provided, the same value will be used for all the colonies.} \item{new}{logical, should the number of individuals be added to the caste population anew or should we only top-up the existing number of individuals to \code{nInd}} -\item{exact}{logical, only relevant when adding workers - if the csd locus is turned -on and exact is \code{TRUE}, we add the exact specified number of viable workers -(heterozygous at the csd locus)} - -\item{year}{numeric, only relevant when adding virgin queens - year of birth for virgin queens} - \item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} \item{...}{additional arguments passed to \code{nInd} when this argument is a function} @@ -90,7 +68,7 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 5, nDrones = nFathersPoisson) # Create and cross Colony and MultiColony class colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) -apiary <- createMultiColony(basePop[4:5], n = 2) +apiary <- createMultiColony(basePop[4:5]) apiary <- cross(apiary, drones = droneGroups[3:4]) #Here we show an example for workers, but same holds for drones and virgin queens! @@ -117,7 +95,7 @@ SP$nWorkers <- nWorkersPoisson # nVirginQueens/nWorkers/nDrones will vary between function calls when a function is used # Queen's counters -getMisc(getQueen(addWorkers(colony))) +getQueen(addWorkers(colony))@misc # Add individuals to a MultiColony object apiary <- addWorkers(apiary) diff --git a/man/addCastePop_internal.Rd b/man/addCastePop_internal.Rd new file mode 100644 index 00000000..bba89e34 --- /dev/null +++ b/man/addCastePop_internal.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Functions_L2_Colony.R +\name{addCastePop_internal} +\alias{addCastePop_internal} +\title{An internal function to add a population in a caste slot of the colony} +\usage{ +addCastePop_internal(pop, colony, caste, new = FALSE) +} +\arguments{ +\item{pop}{\code{\link[AlphaSimR]{Pop-class}} with one or many individual} + +\item{colony}{\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}} + +\item{caste}{character} + +\item{new}{logical} +} +\value{ +\code{\link[SIMplyBee]{Colony-class}} +} +\description{ +Helper function that returns a colony to allow parallelisation, +only for internal use. +} diff --git a/man/buildUp.Rd b/man/buildUp.Rd index 5e280a04..650d8169 100644 --- a/man/buildUp.Rd +++ b/man/buildUp.Rd @@ -9,7 +9,6 @@ buildUp( nWorkers = NULL, nDrones = NULL, new = TRUE, - exact = FALSE, resetEvents = FALSE, simParamBee = NULL, ... @@ -34,10 +33,6 @@ a single value is provided, the same value will be applied to all the colonies.} should we only top-up the existing number of workers and drones to \code{nWorkers} and \code{nDrones} (see details)} -\item{exact}{logical, if the csd locus is turned on and exact is \code{TRUE}, -create the exact specified number of only viable workers (heterozygous on -the csd locus)} - \item{resetEvents}{logical, call \code{\link[SIMplyBee]{resetEvents}} as part of the build up} @@ -86,7 +81,7 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson) colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) isProductive(colony) -apiary <- createMultiColony(basePop[3:4], n = 2) +apiary <- createMultiColony(basePop[3:4]) apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) isProductive(apiary) @@ -118,5 +113,6 @@ nWorkers(apiary) nDrones(apiary) # Queen's counters -getMisc(getQueen(buildUp(colony))) +getQueen(buildUp(colony))@misc + } diff --git a/man/collapse.Rd b/man/collapse.Rd index e00a37b2..34326080 100644 --- a/man/collapse.Rd +++ b/man/collapse.Rd @@ -4,10 +4,12 @@ \alias{collapse} \title{Collapse} \usage{ -collapse(x) +collapse(x, simParamBee = NULL) } \arguments{ \item{x}{\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}} + +\item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} } \value{ \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} with the collapse @@ -30,13 +32,13 @@ founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 50) SP <- SimParamBee$new(founderGenomes) \dontshow{SP$nThreads = 1L} basePop <- createVirginQueens(founderGenomes) -drones <- createDrones(basePop[1], n = 1000) +drones <- createDrones(basePop[1], nInd = 1000) droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10) # Create Colony and MultiColony class colony <- createColony(x = basePop[1]) colony <- cross(colony, drones = droneGroups[[1]]) -apiary <- createMultiColony(x = basePop[2:10], n = 9) +apiary <- createMultiColony(x = basePop[2:10]) apiary <- cross(apiary, drones = droneGroups[2:10]) # Collapse diff --git a/man/combine.Rd b/man/combine.Rd index c14a3a67..9e98b85b 100644 --- a/man/combine.Rd +++ b/man/combine.Rd @@ -4,12 +4,14 @@ \alias{combine} \title{Combine two colony objects} \usage{ -combine(strong, weak) +combine(strong, weak, simParamBee = NULL) } \arguments{ \item{strong}{\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}} \item{weak}{\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}} + +\item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} } \value{ a combined \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} @@ -26,8 +28,9 @@ Level 2 function that combines two Colony or MultiColony objects founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 50) SP <- SimParamBee$new(founderGenomes) \dontshow{SP$nThreads = 1L} +print(SP$nThreads) basePop <- createVirginQueens(founderGenomes) -drones <- createDrones(basePop[1], n = 1000) +drones <- createDrones(basePop[1], nInd = 1000) droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10) # Create weak and strong Colony and MultiColony class @@ -56,7 +59,7 @@ rm(colony2) nWorkers(apiary1); nWorkers(apiary2) nDrones(apiary1); nDrones(apiary2) -apiary1 <- combine(strong = apiary1, weak = apiary2) +apiary1 <- combine(strong = apiary1, weak = apiary2, simParamBee = SP) nWorkers(apiary1); nWorkers(apiary2) nDrones(apiary1); nDrones(apiary2) rm(apiary2) diff --git a/man/createCastePop.Rd b/man/createCastePop.Rd index 4766cbd1..5913ca00 100644 --- a/man/createCastePop.Rd +++ b/man/createCastePop.Rd @@ -11,25 +11,40 @@ createCastePop( x, caste = NULL, nInd = NULL, - exact = TRUE, - year = NULL, editCsd = TRUE, csdAlleles = NULL, simParamBee = NULL, + returnSP = FALSE, + ids = NULL, ... ) -createWorkers(x, nInd = NULL, exact = FALSE, simParamBee = NULL, ...) +createWorkers( + x, + nInd = NULL, + simParamBee = NULL, + returnSP = FALSE, + ids = NULL, + ... +) -createDrones(x, nInd = NULL, simParamBee = NULL, ...) +createDrones( + x, + nInd = NULL, + simParamBee = NULL, + returnSP = FALSE, + ids = NULL, + ... +) createVirginQueens( x, nInd = NULL, - year = NULL, editCsd = TRUE, csdAlleles = NULL, simParamBee = NULL, + returnSP = FALSE, + ids = NULL, ... ) } @@ -47,13 +62,6 @@ only used when \code{x} is \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}, when \code{x} is \code{link[AlphaSimR]{MapPop-class}} all individuals in \code{x} are converted into virgin queens} -\item{exact}{logical, only relevant when creating workers, -if the csd locus is active and exact is \code{TRUE}, -create the exactly specified number of viable workers (heterozygous on the -csd locus)} - -\item{year}{numeric, year of birth for virgin queens} - \item{editCsd}{logical (only active when \code{x} is \code{link[AlphaSimR]{MapPop-class}}), whether the csd locus should be edited to ensure heterozygosity at the csd locus (to get viable virgin queens); see \code{csdAlleles}} @@ -70,6 +78,12 @@ ensure heterozygosity at the csd locus.} \item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} +\item{returnSP}{logical, whether to return the pedigree, caste, and recHist information +for each created population (used internally for parallel computing)} + +\item{ids}{character, IDs of the individuals that are going to be created (used internally +for parallel computing)} + \item{...}{additional arguments passed to \code{nInd} when this argument is a function} } \value{ @@ -100,7 +114,7 @@ Level 1 function that creates the specified number of caste founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 50) SP <- SimParamBee$new(founderGenomes) \dontshow{SP$nThreads = 1L} -SP$setTrackRec(TRUE) +SP$setTrackRec(isTrackRec = TRUE) SP$setTrackPed(isTrackPed = TRUE) # Create virgin queens on a MapPop @@ -118,7 +132,7 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 3, nDrones = nFathersPoisson) # Create a Colony and a MultiColony class colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) -apiary <- createMultiColony(basePop[3:4], n = 2) +apiary <- createMultiColony(basePop[3:4]) apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) # Using default nInd in SP diff --git a/man/createColony.Rd b/man/createColony.Rd index a8a96649..d1707f07 100644 --- a/man/createColony.Rd +++ b/man/createColony.Rd @@ -4,12 +4,14 @@ \alias{createColony} \title{Create a new Colony} \usage{ -createColony(x = NULL, simParamBee = NULL) +createColony(x = NULL, simParamBee = NULL, id = NULL) } \arguments{ \item{x}{\code{\link[AlphaSimR]{Pop-class}}, one queen or virgin queen(s)} \item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} + +\item{id}{character, ID of the colony that is going to be created (used internally for parallel computing)} } \value{ new \code{\link[SIMplyBee]{Colony-class}} diff --git a/man/createCrossPlan.Rd b/man/createCrossPlan.Rd index ae90395f..8db00494 100644 --- a/man/createCrossPlan.Rd +++ b/man/createCrossPlan.Rd @@ -69,10 +69,6 @@ virginColonies2 <- createMultiColony(basePop[31:60]) virginColonies2 <- setLocation(virginColonies2, location = Map(c, runif(30, 0, 2*pi), runif(30, 0, 2*pi))) -virginColonies3 <- createMultiColony(basePop[61:90]) -virginColonies3 <- setLocation(virginColonies3, - location = Map(c, runif(30, 0, 2*pi), - runif(30, 0, 2*pi))) # Create drone colonies droneColonies <- createMultiColony(basePop[121:200]) @@ -92,57 +88,22 @@ droneColonies <- cross(droneColonies, nDrones = nFathersPoisson, crossPlan = randomCrossPlan) -# Plot the colonies in space -virginLocations <- as.data.frame(getLocation(c(virginColonies1, virginColonies2, virginColonies3), - collapse= TRUE)) -virginLocations$Type <- "Virgin" -droneLocations <- as.data.frame(getLocation(droneColonies, collapse= TRUE)) -droneLocations$Type <- "Drone" -locations <- rbind(virginLocations, droneLocations) - -plot(x = locations$V1, y = locations$V2, - col = c("red", "blue")[as.numeric(as.factor(locations$Type))]) - -# Cross according to a spatial cross plan according to the colonies' locations -crossPlanSpatial <- createCrossPlan(x = virginColonies1, - droneColonies = droneColonies, - nDrones = nFathersPoisson, - spatial = TRUE, - radius = 1.5) - -# Plot the crossing for the first colony in the crossPlan -virginLocations1 <- as.data.frame(getLocation(virginColonies1, collapse= TRUE)) -virginLocations1$Type <- "Virgin" -droneLocations <- as.data.frame(getLocation(droneColonies, collapse= TRUE)) -droneLocations$Type <- "Drone" -locations1 <- rbind(virginLocations1, droneLocations) - -# Blue marks the target virgin colony and blue marks the drone colonies in the chosen radius -plot(x = locations1$V1, y = locations1$V2, pch = c(1, 2)[as.numeric(as.factor(locations1$Type))], - col = ifelse(rownames(locations1) \%in\% crossPlanSpatial[[1]], - "red", - ifelse(rownames(locations1) == names(crossPlanSpatial)[[1]], - "blue", "black"))) - -colonies1 <- cross(x = virginColonies1, - crossPlan = crossPlanSpatial, - droneColonies = droneColonies, - nDrones = nFathersPoisson) -nFathers(colonies1) - # Cross according to a cross plan that is created internally within the cross function # The cross plan is created at random, regardless the location of the colonies -colonies2 <- cross(x = virginColonies2, +colonies1 <- cross(x = virginColonies1, droneColonies = droneColonies, nDrones = nFathersPoisson, crossPlan = "create") -# Mate spatially with cross plan created internally by the cross function -colonies3 <- cross(x = virginColonies3, - droneColonies = droneColonies, + +# Cross according to a spatial cross plan created internally according to the colonies' locations +colonies2 <- cross(x = virginColonies2, crossPlan = "create", - checkCross = "warning", + droneColonies = droneColonies, + nDrones = nFathersPoisson, spatial = TRUE, - radius = 1) + radius = 1.5) +nFathers(colonies2) + } diff --git a/man/createMatingStationDCA.Rd b/man/createMatingStationDCA.Rd index 8ccf320d..58c138a7 100644 --- a/man/createMatingStationDCA.Rd +++ b/man/createMatingStationDCA.Rd @@ -35,7 +35,7 @@ founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes) \dontshow{SP$nThreads = 1L} basePop <- createVirginQueens(founderGenomes) -drones <- createDrones(basePop[1], n = 1000) +drones <- createDrones(basePop[1], nInd = 1000) droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10) # Create a colony and cross it diff --git a/man/createMultiColony.Rd b/man/createMultiColony.Rd index 21f2bc14..7951d430 100644 --- a/man/createMultiColony.Rd +++ b/man/createMultiColony.Rd @@ -4,7 +4,12 @@ \alias{createMultiColony} \title{Create MultiColony object} \usage{ -createMultiColony(x = NULL, n = NULL, simParamBee = NULL) +createMultiColony( + x = NULL, + n = NULL, + simParamBee = NULL, + populateColonies = FALSE +) } \arguments{ \item{x}{\code{\link[AlphaSimR]{Pop-class}}, virgin queens or queens for the colonies @@ -16,6 +21,8 @@ given then \code{\link[SIMplyBee]{MultiColony-class}} is created with \code{n} \code{NULL}) individual colony - this is mostly useful for programming)} \item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} + +\item{populateColonies}{boolean, whether to create n empty Colony objects within with assigned ID} } \value{ \code{\link[SIMplyBee]{MultiColony-class}} @@ -49,7 +56,7 @@ apiary[[2]] # Create mated colonies by crossing apiary <- createMultiColony(x = basePop[1:2], n = 2) -drones <- createDrones(x = basePop[3], n = 30) +drones <- createDrones(x = basePop[3], nInd = 30) droneGroups <- pullDroneGroupsFromDCA(drones, n = 2, nDrones = 15) apiary <- cross(apiary, drones = droneGroups) apiary diff --git a/man/cross.Rd b/man/cross.Rd index 8b4ee5fe..8d34e3a6 100644 --- a/man/cross.Rd +++ b/man/cross.Rd @@ -54,7 +54,8 @@ to their distance from the virgin colony (that is, in a radius)} \item{radius}{numeric, the radius around the virgin colony in which to sample mating partners, only needed when \code{spatial = TRUE}} -\item{checkCross}{character, throw a warning (when \code{checkCross = "warning"}),} +\item{checkCross}{character, throw a warning (when \code{checkCross = "warning"}). +This will also remove the unmated queens and return only the mated ones.} \item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} @@ -95,7 +96,7 @@ This function changes caste for the mated drones to fathers, and \examples{ founderGenomes <- quickHaplo(nInd = 30, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes) -\dontshow{SP$nThreads = 1L} +SP$nThreads = 1L basePop <- createVirginQueens(founderGenomes) drones <- createDrones(x = basePop[1], nInd = 1000) diff --git a/man/downsize.Rd b/man/downsize.Rd index e418ad0b..d2acf269 100644 --- a/man/downsize.Rd +++ b/man/downsize.Rd @@ -48,7 +48,7 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 3, nDrones = 12) colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) colony <- buildUp(colony) -apiary <- createMultiColony(basePop[3:4], n = 2) +apiary <- createMultiColony(basePop[3:4]) apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) apiary <- buildUp(apiary) diff --git a/man/getIbdHaplo.Rd b/man/getIbdHaplo.Rd index fd6dfe79..08895a7a 100644 --- a/man/getIbdHaplo.Rd +++ b/man/getIbdHaplo.Rd @@ -115,10 +115,10 @@ Level 0 function that returns IBD (identity by descent) }} \examples{ -founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 50) -SP <- SimParamBee$new(founderGenomes) +founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 5) +SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 4) \dontshow{SP$nThreads = 1L} -SP$setTrackRec(TRUE) +SP$setTrackRec(isTrackRec = TRUE) SP$setTrackPed(isTrackPed = TRUE) basePop <- createVirginQueens(founderGenomes) @@ -128,13 +128,13 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson) # Create a Colony and a MultiColony class colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) -colony <- buildUp(x = colony, nWorkers = 6, nDrones = 3) -colony <- addVirginQueens(x = colony, nInd = 5) +colony <- buildUp(x = colony, nWorkers = 3, nDrones = 2) +colony <- addVirginQueens(x = colony, nInd = 2) -apiary <- createMultiColony(basePop[3:4], n = 2) +apiary <- createMultiColony(basePop[3:4]) apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) -apiary <- buildUp(x = apiary, nWorkers = 6, nDrones = 3) -apiary <- addVirginQueens(x = apiary, nInd = 5) +apiary <- buildUp(x = apiary, nWorkers = 3, nDrones = 2) +apiary <- addVirginQueens(x = apiary, nInd = 2) # Input is a population getIbdHaplo(x = getQueen(colony)) @@ -146,6 +146,8 @@ getIbdHaplo(x = colony, caste = "queen") getQueenIbdHaplo(colony) getIbdHaplo(colony, caste = "workers", nInd = 3) +getIbdHaplo(colony, caste = "virginQueens") +getIbdHaplo(colony, caste = "drones") getWorkersIbdHaplo(colony) # Same aliases exist for all castes! @@ -160,6 +162,9 @@ getQueenIbdHaplo(apiary) # Or collapse all the haplotypes into a single matrix getQueenIbdHaplo(apiary, collapse = TRUE) + +getIbdHaplo(x = apiary, caste = "workers") +getIbdHaplo(x = apiary, caste = "drones") # Get the haplotypes of all individuals either by colony or in a single matrix getIbdHaplo(apiary, caste = "all") getIbdHaplo(apiary, caste = "all", collapse = TRUE) diff --git a/man/getMisc.Rd b/man/getMisc.Rd deleted file mode 100644 index 60b342d8..00000000 --- a/man/getMisc.Rd +++ /dev/null @@ -1,20 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Functions_L0_auxilary.R -\name{getMisc} -\alias{getMisc} -\title{Get miscellaneous information in a population} -\usage{ -getMisc(x, node = NULL) -} -\arguments{ -\item{x}{\code{\link[AlphaSimR]{Pop-class}}} - -\item{node}{character, name of the node to get from the \code{x@misc} slot; -if \code{NULL} the whole \code{x@misc} slot is returned} -} -\value{ -The \code{x@misc} slot or its nodes \code{x@misc[[*]][[node]]} -} -\description{ -Get miscellaneous information in a population -} diff --git a/man/getQueenAge.Rd b/man/getQueenAge.Rd deleted file mode 100644 index dba30d22..00000000 --- a/man/getQueenAge.Rd +++ /dev/null @@ -1,48 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Functions_L0_auxilary.R -\name{getQueenAge} -\alias{getQueenAge} -\title{Get (calculate) the queen's age} -\usage{ -getQueenAge(x, currentYear, simParamBee = NULL) -} -\arguments{ -\item{x}{\code{\link[AlphaSimR]{Pop-class}}, \code{\link[SIMplyBee]{Colony-class}}, or -\code{\link[SIMplyBee]{MultiColony-class}}} - -\item{currentYear}{integer, current year} - -\item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} -} -\value{ -numeric, the age of the queen(s); named when theres is more - than one queen; \code{NA} if queen not present -} -\description{ -Level 0 function that returns the queen's age. -} -\examples{ -founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) -SP <- SimParamBee$new(founderGenomes) -\dontshow{SP$nThreads = 1L} -basePop <- createVirginQueens(founderGenomes) - -drones <- createDrones(x = basePop[1], nInd = 1000) -droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson) - -# Create a Colony and a MultiColony class -colony <- createColony(x = basePop[2]) -colony <- cross(colony, drones = droneGroups[[1]]) -apiary <- createMultiColony(basePop[3:4], n = 2) -apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) - -queen <- getQueen(colony) -queen <- setQueensYearOfBirth(queen, year = 2020) -getQueenAge(queen, currentYear = 2022) - -colony <- setQueensYearOfBirth(colony, year = 2021) -getQueenAge(colony, currentYear = 2022) - -apiary <- setQueensYearOfBirth(apiary, year = 2018) -getQueenAge(apiary, currentYear = 2022) -} diff --git a/man/getQueenYearOfBirth.Rd b/man/getQueenYearOfBirth.Rd deleted file mode 100644 index 0db401d3..00000000 --- a/man/getQueenYearOfBirth.Rd +++ /dev/null @@ -1,49 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Functions_L0_auxilary.R -\name{getQueenYearOfBirth} -\alias{getQueenYearOfBirth} -\title{Access the queen's year of birth} -\usage{ -getQueenYearOfBirth(x, simParamBee = NULL) -} -\arguments{ -\item{x}{\code{\link[AlphaSimR]{Pop-class}} (one or more than one queen), -\code{\link[SIMplyBee]{Colony-class}} (one colony), or -\code{\link[SIMplyBee]{MultiColony-class}} (more colonies)} - -\item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} -} -\value{ -numeric, the year of birth of the queen(s); named when theres is more - than one queen; \code{NA} if queen not present -} -\description{ -Level 0 function that returns the queen's year of birth. -} -\examples{ -founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) -SP <- SimParamBee$new(founderGenomes) -\dontshow{SP$nThreads = 1L} -basePop <- createVirginQueens(founderGenomes) - -drones <- createDrones(x = basePop[1], nInd = 1000) -droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson) - -# Create a Colony and a MultiColony class -colony <- createColony(x = basePop[2]) -colony <- cross(colony, drones = droneGroups[[1]]) - -apiary <- createMultiColony(basePop[3:4], n = 2) -apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) - -queen <- getQueen(colony) -queen <- setQueensYearOfBirth(queen, year = 2022) -getQueenYearOfBirth(queen) - -getQueenYearOfBirth(getQueen(colony)) -colony <- setQueensYearOfBirth(colony, year = 2030) -getQueenYearOfBirth(colony) - -apiary <- setQueensYearOfBirth(apiary, year = 2022) -getQueenYearOfBirth(apiary) -} diff --git a/man/hasSwarmed.Rd b/man/hasSwarmed.Rd index 870d2d14..427d41d3 100644 --- a/man/hasSwarmed.Rd +++ b/man/hasSwarmed.Rd @@ -31,7 +31,7 @@ colony <- cross(colony, drones = droneGroups[[1]]) colony <- buildUp(x = colony, nWorkers = 6, nDrones = 3) colony <- addVirginQueens(colony, nInd = 5) -apiary <- createMultiColony(basePop[3:4], n = 2) +apiary <- createMultiColony(basePop[3:4]) apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) apiary <- buildUp(x = apiary, nWorkers = 6, nDrones = 3) diff --git a/man/isCsdHeterozygous.Rd b/man/isCsdHeterozygous.Rd index 7a5d9d26..1bd9f5ed 100644 --- a/man/isCsdHeterozygous.Rd +++ b/man/isCsdHeterozygous.Rd @@ -17,7 +17,8 @@ logical \description{ Level 0 function that returns if individuals of a population are heterozygous at the csd locus. See \code{\link[SIMplyBee]{SimParamBee}} for more - information about the csd locus. + information about the csd locus. The function also return \code{TRUE} for drones to + mark their viability, although they are haploid. } \details{ We could expand \code{isCsdHeterozygous} to work also with diff --git a/man/mapCasteToColonyValue.Rd b/man/mapCasteToColonyValue.Rd index 290d5e5d..3742a6cb 100644 --- a/man/mapCasteToColonyValue.Rd +++ b/man/mapCasteToColonyValue.Rd @@ -42,21 +42,24 @@ mapCasteToColonyAa(colony, simParamBee = NULL, ...) \item{queenTrait}{numeric (column position) or character (column name), trait(s) that represents queen's contribution to colony value(s); if -\code{NULL} then this contribution is 0; you can pass more than one trait +\code{NULL} or there is no queen present, then this contribution is 0; +you can pass more than one trait here, but make sure that \code{combineFUN} works with these trait dimensions} -\item{queenFUN}{function, function that will be applied to queen's value} +\item{queenFUN}{function, function that will be applied to queen's value.} \item{workersTrait}{numeric (column position) or character (column name), trait(s) that represents workers' contribution to colony value(s); if -\code{NULL} then this contribution is 0; you can pass more than one trait +\code{NULL} or there are no workers present, then this contribution is 0; +you can pass more than one trait here, but make sure that \code{combineFUN} works with these trait dimensions} \item{workersFUN}{function, function that will be applied to workers values} \item{dronesTrait}{numeric (column position) or character (column name), trait(s) that represents drones' contribution to colony value(s); if -\code{NULL} then this contribution is 0; you can pass more than one trait +\code{NULL} or there are no drones present then this contribution is 0; +you can pass more than one trait here, but make sure that \code{combineFUN} works with these trait dimensions} \item{dronesFUN}{function, function that will be applied to drone values} diff --git a/man/pullCastePop.Rd b/man/pullCastePop.Rd index dcf63748..42a890b2 100644 --- a/man/pullCastePop.Rd +++ b/man/pullCastePop.Rd @@ -95,12 +95,12 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson) # Create a Colony and a MultiColony class colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) -colony <- buildUp(x = colony, nWorkers = 100, nDrones = 10, exact = TRUE) +colony <- buildUp(x = colony, nWorkers = 100, nDrones = 10) colony <- addVirginQueens(x = colony, nInd = 3) apiary <- createMultiColony(basePop[3:4], n = 2) apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) -apiary <- buildUp(x = apiary, nWorkers = 100, nDrones = 10, exact = TRUE) +apiary <- buildUp(x = apiary, nWorkers = 100, nDrones = 10) apiary <- addVirginQueens(x = apiary, nInd = 3) # pullCastePop on Colony class diff --git a/man/reQueen.Rd b/man/reQueen.Rd index e90abb52..61e22665 100644 --- a/man/reQueen.Rd +++ b/man/reQueen.Rd @@ -46,7 +46,7 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 7, nDrones = nFathersPoisson) # Create and cross Colony and MultiColony class colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) -apiary <- createMultiColony(basePop[3:4], n = 2) +apiary <- createMultiColony(basePop[3:4]) apiary <- cross(apiary, drones = droneGroups[2:3]) # Check queen and virgin queens IDs diff --git a/man/removeCastePop.Rd b/man/removeCastePop.Rd index 9acecac0..cbe9c2ec 100644 --- a/man/removeCastePop.Rd +++ b/man/removeCastePop.Rd @@ -8,24 +8,9 @@ \alias{removeVirginQueens} \title{Remove a proportion of caste individuals from a colony} \usage{ -removeCastePop( - x, - caste = NULL, - p = 1, - use = "rand", - addVirginQueens = FALSE, - nVirginQueens = NULL, - year = NULL, - simParamBee = NULL -) - -removeQueen( - x, - addVirginQueens = FALSE, - nVirginQueens = NULL, - year = NULL, - simParamBee = NULL -) +removeCastePop(x, caste = NULL, p = 1, use = "rand", simParamBee = NULL) + +removeQueen(x, simParamBee = NULL) removeWorkers(x, p = 1, use = "rand", simParamBee = NULL) @@ -45,16 +30,6 @@ a single value is provided, the same value will be applied to all the colonies} \item{use}{character, all the options provided by \code{\link[AlphaSimR]{selectInd}} - guides selection of virgins queens that will stay when \code{p < 1}} -\item{addVirginQueens}{logical, whether virgin queens should be added; only -used when removing the queen from the colony} - -\item{nVirginQueens}{integer, the number of virgin queens to be created in the -colony; only used when removing the queen from the colony. If \code{0}, no virgin -queens are added; If \code{NULL}, the value from \code{simParamBee$nVirginQueens} -is used} - -\item{year}{numeric, only relevant when adding virgin queens - year of birth for virgin queens} - \item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} } \value{ @@ -88,7 +63,7 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 5, nDrones = nFathersPoisson) colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) colony <- buildUp(colony) -apiary <- createMultiColony(basePop[4:5], n = 2) +apiary <- createMultiColony(basePop[4:5]) apiary <- cross(apiary, drones = droneGroups[3:4]) apiary <- buildUp(apiary) diff --git a/man/replaceCastePop.Rd b/man/replaceCastePop.Rd index bd2c3756..44fd8728 100644 --- a/man/replaceCastePop.Rd +++ b/man/replaceCastePop.Rd @@ -7,17 +7,9 @@ \alias{replaceVirginQueens} \title{Replace a proportion of caste individuals with new ones} \usage{ -replaceCastePop( - x, - caste = NULL, - p = 1, - use = "rand", - exact = TRUE, - year = NULL, - simParamBee = NULL -) +replaceCastePop(x, caste = NULL, p = 1, use = "rand", simParamBee = NULL) -replaceWorkers(x, p = 1, use = "rand", exact = TRUE, simParamBee = NULL) +replaceWorkers(x, p = 1, use = "rand", simParamBee = NULL) replaceDrones(x, p = 1, use = "rand", simParamBee = NULL) @@ -36,14 +28,6 @@ a single value is provided, the same value will be applied to all the colonies} \item{use}{character, all the options provided by \code{\link[AlphaSimR]{selectInd}} - guides selection of caste individuals that stay when \code{p < 1}} -\item{exact}{logical, only relevant when adding workers - if the csd locus is turned -on and exact is \code{TRUE}, we replace the exact specified number of viable workers -(heterozygous at the csd locus). You probably want this set to TRUE since you want to -replace with the same number of workers.} - -\item{year}{numeric, only relevant when replacing virgin queens, -year of birth for virgin queens} - \item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} } \value{ @@ -77,7 +61,7 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 5, nDrones = nFathersPoisson) # Create and cross Colony and MultiColony class colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) -apiary <- createMultiColony(basePop[4:5], n = 2) +apiary <- createMultiColony(basePop[4:5]) apiary <- cross(apiary, drones = droneGroups[3:4]) # Add individuals diff --git a/man/resetEvents.Rd b/man/resetEvents.Rd index 2e8b6642..27836670 100644 --- a/man/resetEvents.Rd +++ b/man/resetEvents.Rd @@ -4,7 +4,7 @@ \alias{resetEvents} \title{Reset colony events} \usage{ -resetEvents(x, collapse = NULL) +resetEvents(x, collapse = NULL, simParamBee = NULL) } \arguments{ \item{x}{\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}} @@ -13,6 +13,8 @@ resetEvents(x, collapse = NULL) up a new colony, which the default of \code{NULL} caters for; otherwise, a collapsed colony should be left collapsed forever, unless you force resetting this event with \code{collapse = TRUE})} + +\item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} } \value{ \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} with @@ -36,7 +38,7 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 5, nDrones = nFathersPoisson) # Create and cross Colony and MultiColony class colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) -apiary <- createMultiColony(basePop[4:5], n = 2) +apiary <- createMultiColony(basePop[4:5]) apiary <- cross(apiary, drones = droneGroups[3:4]) # Build-up - this sets Productive to TRUE diff --git a/man/setEvents.Rd b/man/setEvents.Rd new file mode 100644 index 00000000..5035cd70 --- /dev/null +++ b/man/setEvents.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Functions_L2_Colony.R +\name{setEvents} +\alias{setEvents} +\title{Set colony events} +\usage{ +setEvents(x, slot, value, simParamBee = NULL) +} +\arguments{ +\item{x}{\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}} + +\item{slot}{character, which event to set} + +\item{value}{logical, the value for the event} + +\item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} +} +\value{ +\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} with + events reset +} +\description{ +Helper Level 2 function that populates the events slot. Not interded + for external use, intended for internal use in parallel computing +} +\examples{ +founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 50) +SP <- SimParamBee$new(founderGenomes) +\dontshow{SP$nThreads = 1L} +basePop <- createVirginQueens(founderGenomes) + +drones <- createDrones(x = basePop[1], nInd = 100) +droneGroups <- pullDroneGroupsFromDCA(drones, n = 5, nDrones = nFathersPoisson) + +# Create and cross Colony and MultiColony class +colony <- createColony(x = basePop[2]) +colony <- cross(colony, drones = droneGroups[[1]]) +apiary <- createMultiColony(basePop[4:5]) +SIMplyBee:::setEvents(apiary, slot = "swarm", value = c(TRUE, TRUE)) +} diff --git a/man/setLocation.Rd b/man/setLocation.Rd index a9ee600d..e9596075 100644 --- a/man/setLocation.Rd +++ b/man/setLocation.Rd @@ -4,7 +4,7 @@ \alias{setLocation} \title{Set colony location} \usage{ -setLocation(x, location = c(0, 0)) +setLocation(x, location = c(0, 0), simParamBee = NULL) } \arguments{ \item{x}{\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}} @@ -14,6 +14,8 @@ locations as \code{c(x1, y1)} (the same location set to all colonies), \code{list(c(x1, y1), c(x2, y2))}, or \code{data.frame(x = c(x1, x2), y = c(y1, y2))}} + +\item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} } \value{ \code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}} with set diff --git a/man/setMisc.Rd b/man/setMisc.Rd deleted file mode 100644 index 0b2c0291..00000000 --- a/man/setMisc.Rd +++ /dev/null @@ -1,26 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Functions_L0_auxilary.R -\name{setMisc} -\alias{setMisc} -\title{Set miscellaneous information in a population} -\usage{ -setMisc(x, node = NULL, value = NULL) -} -\arguments{ -\item{x}{\code{\link[AlphaSimR]{Pop-class}}} - -\item{node}{character, name of the node to set within the \code{x@misc} slot} - -\item{value, }{value to be saved into \code{x@misc[[*]][[node]]}; length of -\code{value} should be equal to \code{nInd(x)}; if its length is 1, then -it is repeated using \code{rep} (see examples)} -} -\value{ -\code{\link[AlphaSimR]{Pop-class}} -} -\description{ -Set miscellaneous information in a population -} -\details{ -A \code{NULL} in \code{value} is ignored -} diff --git a/man/setQueensYearOfBirth.Rd b/man/setQueensYearOfBirth.Rd deleted file mode 100644 index 094eea2f..00000000 --- a/man/setQueensYearOfBirth.Rd +++ /dev/null @@ -1,53 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Functions_L1_Pop.R -\name{setQueensYearOfBirth} -\alias{setQueensYearOfBirth} -\title{Set the queen's year of birth} -\usage{ -setQueensYearOfBirth(x, year, simParamBee = NULL) -} -\arguments{ -\item{x}{\code{\link[AlphaSimR]{Pop-class}} (one or more than one queen), -\code{\link[SIMplyBee]{Colony-class}} (one colony), or -\code{\link[SIMplyBee]{MultiColony-class}} (more colonies)} - -\item{year}{integer, the year of the birth of the queen} - -\item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} -} -\value{ -\code{\link[AlphaSimR]{Pop-class}}, \code{\link[SIMplyBee]{Colony-class}}, or - \code{\link[SIMplyBee]{MultiColony-class}} with queens having the year of birth set -} -\description{ -Level 1 function that sets the queen's year of birth. -} -\examples{ -founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) -SP <- SimParamBee$new(founderGenomes) -\dontshow{SP$nThreads = 1L} -basePop <- createVirginQueens(founderGenomes) - -drones <- createDrones(x = basePop[1], nInd = 1000) -droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson) - -# Create a Colony and a MultiColony class -colony <- createColony(x = basePop[2]) -colony <- cross(x = colony, drones = droneGroups[[1]]) -apiary <- createMultiColony(basePop[3:4], n = 2) -apiary <- cross(apiary, drones = droneGroups[c(2, 3)]) - -# Example on Colony class -getQueenYearOfBirth(colony) -getQueenYearOfBirth(apiary) - -queen1 <- getQueen(colony) -queen1 <- setQueensYearOfBirth(queen1, year = 2022) -getQueenYearOfBirth(queen1) - -colony <- setQueensYearOfBirth(colony, year = 2022) -getQueenYearOfBirth(colony) - -apiary <- setQueensYearOfBirth(apiary, year = 2022) -getQueenYearOfBirth(apiary) -} diff --git a/man/split.Rd b/man/split.Rd index 7def12fb..0fab0d75 100644 --- a/man/split.Rd +++ b/man/split.Rd @@ -4,7 +4,7 @@ \alias{split} \title{Split colony in two MultiColony} \usage{ -split(x, p = NULL, year = NULL, simParamBee = NULL, ...) +split(x, p = NULL, simParamBee = NULL, ...) } \arguments{ \item{x}{\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}} @@ -15,8 +15,6 @@ If input is \code{\link[SIMplyBee]{MultiColony-class}}, the input could also be a vector of the same length as the number of colonies. If a single value is provided, the same value will be applied to all the colonies} -\item{year}{numeric, year of birth for virgin queens} - \item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} \item{...}{additional arguments passed to \code{p} when this argument is a @@ -32,8 +30,9 @@ Level 2 function that splits a Colony or MultiColony object into two new colonies to prevent swarming (in managed situation). The remnant colony retains the queen and a proportion of the workers and all drones. The split colony gets - the other part of the workers, which raise virgin queens, of which only one - prevails. Location of the split is the same as for the remnant. + the other part of the workers, but note that it is queenless, since the beekeepers + would normally requeen with a different queen. + Location of the split is the same as for the remnant. } \examples{ founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 50) diff --git a/man/supersede.Rd b/man/supersede.Rd index 04291135..bc1de781 100644 --- a/man/supersede.Rd +++ b/man/supersede.Rd @@ -4,18 +4,11 @@ \alias{supersede} \title{Supersede} \usage{ -supersede(x, year = NULL, nVirginQueens = NULL, simParamBee = NULL, ...) +supersede(x, simParamBee = NULL, ...) } \arguments{ \item{x}{\code{\link[SIMplyBee]{Colony-class}} or \code{\link[SIMplyBee]{MultiColony-class}}} -\item{year}{numeric, year of birth for virgin queens} - -\item{nVirginQueens}{integer, the number of virgin queens to be created in the -colony; of these one is randomly selected as the new virgin queen of the -remnant colony. If \code{NULL}, the value from \code{simParamBee$nVirginQueens} -is used} - \item{simParamBee}{\code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters} \item{...}{additional arguments passed to \code{nVirginQueens} when this @@ -29,7 +22,11 @@ supersede event set to \code{TRUE} Level 2 function that supersedes a Colony or MultiColony object - an event where the queen dies. The workers and drones stay unchanged, but workers raise virgin - queens, of which only one prevails. + queens, of which only one prevails. The function + will create either 10 or \code{SimParamBee$nVirginQueens} virgin queens, + whichever is higher, and select one at random In case of high inbreeding, + it could be that none of the virgin queens are viable.In that case, you might + want to increase \code{SimParamBee$nVirginQueens} or discard the colony. } \examples{ founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 50) diff --git a/man/swarm.Rd b/man/swarm.Rd index e178fe26..6c6293f1 100644 --- a/man/swarm.Rd +++ b/man/swarm.Rd @@ -7,8 +7,6 @@ swarm( x, p = NULL, - year = NULL, - nVirginQueens = NULL, sampleLocation = TRUE, radius = NULL, simParamBee = NULL, @@ -24,13 +22,6 @@ If input is \code{\link[SIMplyBee]{MultiColony-class}}, the input could also be a vector of the same length as the number of colonies. If a single value is provided, the same value will be applied to all the colonies} -\item{year}{numeric, year of birth for virgin queens} - -\item{nVirginQueens}{integer, the number of virgin queens to be created in the -colony; of these one is randomly selected as the new virgin queen of the -remnant colony. If \code{NULL}, the value from \code{simParamBee$nVirginQueens} -is used} - \item{sampleLocation}{logical, sample location of the swarm by taking the current colony location and adding deviates to each coordinate using \code{\link[SIMplyBee]{rcircle}}} @@ -54,13 +45,19 @@ Level 2 function that swarms a Colony or MultiColony object - an event where the queen leaves with a proportion of workers to create a new colony (the swarm). The remnant colony retains the other proportion of workers and all drones, and - the workers raise virgin queens, of which only one prevails. Location of + the workers raise virgin queens, of which only one prevails. The function + will create either 10 or \code{SimParamBee$nVirginQueens} virgin queens, + whichever is higher, and select one at random. In case of high inbreeding, + it could be that none of the virgin queens are viable. In that case, you might + want to increase \code{SimParamBee$nVirginQueens} or discard the colony. Location of the swarm is the same as for the remnant or sampled as deviation from the remnant. } \examples{ founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 50) SP <- SimParamBee$new(founderGenomes) +SP$setTrackPed(TRUE) +SP$setTrackRec(TRUE) \dontshow{SP$nThreads = 1L} basePop <- createVirginQueens(founderGenomes) drones <- createDrones(basePop[1], n = 1000) @@ -70,7 +67,7 @@ droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10) colony <- createColony(x = basePop[2]) colony <- cross(colony, drones = droneGroups[[1]]) (colony <- buildUp(colony, nWorkers = 100)) -apiary <- createMultiColony(basePop[3:8], n = 6) +apiary <- createMultiColony(basePop[3:8]) apiary <- cross(apiary, drones = droneGroups[2:7]) apiary <- buildUp(apiary, nWorkers = 100) diff --git a/tests/testthat/test-L0_auxiliary_functions.R b/tests/testthat/test-L0_auxiliary_functions.R index b9883b0a..5b32280a 100644 --- a/tests/testthat/test-L0_auxiliary_functions.R +++ b/tests/testthat/test-L0_auxiliary_functions.R @@ -1,8 +1,9 @@ # ---- nColonies ---- + test_that("nColonies", { founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) - SP$nThreads = 1L + SP$nThreads <- 1L basePop <- createVirginQueens(founderGenomes, simParamBee = SP) expect_equal(nColonies(createMultiColony(n = 2, simParamBee = SP)), 2) expect_equal(nColonies(createMultiColony(simParamBee = SP)), 0) @@ -14,6 +15,7 @@ test_that("nColonies", { test_that("nCaste", { founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) + SP$nThreads <- 1L basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 45, simParamBee = SP) @@ -26,12 +28,12 @@ test_that("nCaste", { expect_equal(nCaste(colony, caste = "virginQueens", simParamBee = SP), 0) expect_equal(nCaste(colony, caste = "fathers", simParamBee = SP), 10) - apiary <- createMultiColony(basePop[3:4], n = 2, simParamBee = SP) - apiary <- cross(apiary, drones = droneGroups[c(2, 3)], simParamBee = SP) - apiary <- buildUp(apiary, nWorkers = 20, nDrones = 10, simParamBee = SP) - expect_equal(sum(nCaste(apiary, caste = "queen", simParamBee = SP)), 2) - expect_equal(sum(nCaste(apiary, caste = "virginQueens", simParamBee = SP)), 0) - expect_equal(sum(nCaste(apiary, caste = "fathers", simParamBee = SP)), 20) + #apiary <- createMultiColony(basePop[3:4], n = 2, simParamBee = SP) + #apiary <- cross(apiary, drones = droneGroups[c(2, 3)], simParamBee = SP) + #apiary <- buildUp(apiary, nWorkers = 20, nDrones = 10, simParamBee = SP) + #expect_equal(sum(nCaste(apiary, caste = "queen", simParamBee = SP)), 2) + #expect_equal(sum(nCaste(apiary, caste = "virginQueens", simParamBee = SP)), 0) + #expect_equal(sum(nCaste(apiary, caste = "fathers", simParamBee = SP)), 20) }) # ---- nQueens ---- @@ -40,6 +42,7 @@ test_that("nQueens", { founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 1000, simParamBee = SP) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) @@ -60,6 +63,7 @@ test_that("nDrones", { founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 1000, simParamBee = SP) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) @@ -103,6 +107,7 @@ test_that("isCaste", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -136,13 +141,14 @@ test_that("calcQueensPHomBrood", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) # Create a Colony class object - colony <- createColony(x = basePop[2], simParamBee = SP) + colony <- createColony(x = basePop[1], simParamBee = SP) colony <- cross(colony, drones = fatherGroups[[1]], simParamBee = SP) colony <- buildUp(x = colony, nWorkers = 120, nDrones = 20, simParamBee = SP) colony <- addVirginQueens(x = colony, nInd = 1, simParamBee = SP) @@ -150,6 +156,7 @@ test_that("calcQueensPHomBrood", { expect_error(calcQueensPHomBrood(colony@drones, simParamBee = SP)) expect_error(calcQueensPHomBrood(colony@workers, simParamBee = SP)) expect_true(is.numeric(calcQueensPHomBrood(colony@queen, simParamBee = SP))) + expect_true(calcQueensPHomBrood(colony@queen, simParamBee = SP) > 0) colony@queen <- NULL expect_error(calcQueensPHomBrood(colony@queen, simParamBee = SP)) @@ -159,6 +166,8 @@ test_that("calcQueensPHomBrood", { colony@virginQueens <- NULL expect_error(calcQueensPHomBrood(colony, simParamBee = SP)) expect_equal((length(calcQueensPHomBrood(apiary, simParamBee = SP))), 0) + + }) # ---- pHomBrood ---- @@ -167,6 +176,7 @@ test_that("pHomBrood", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -199,6 +209,7 @@ test_that("nHomBrood", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -215,13 +226,12 @@ test_that("nHomBrood", { expect_error(nHomBrood(colony@drones, simParamBee = SP)) expect_true(is.numeric(nHomBrood(colony@queen, simParamBee = SP))) - colony@queen <- NULL - expect_error(nHomBrood(colony@queen, simParamBee = SP)) + apiary <- createMultiColony(simParamBee = SP) colony@workers <- NULL colony@drones <- NULL colony@virginQueens <- NULL - expect_error(nHomBrood(colony, simParamBee = SP)) + expect_error(nHomBrood(removeQueen(colony), simParamBee = SP)) expect_equal(length(nHomBrood(apiary, simParamBee = SP)), 0) }) @@ -231,6 +241,7 @@ test_that("isQueenPresent", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -259,6 +270,7 @@ test_that("isVirginQueensPresent", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -289,6 +301,7 @@ test_that("isProductive", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -323,6 +336,7 @@ test_that("reduceDroneHaplo", { founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 2, simParamBee = SP) virginQueens <- c(basePop[2:3]) @@ -347,6 +361,7 @@ test_that("reduceDroneGeno", { founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 2, simParamBee = SP) virginQueens <- c(basePop[2:3]) @@ -370,9 +385,14 @@ test_that("getCsdAlleles", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) + baseAlleles <- getCsdAlleles(basePop, simParamBee = SP) + expect_equal(nrow(baseAlleles), 8 * 2) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) + dronesAlleles <- getCsdAlleles(drones, simParamBee = SP) + expect_equal(nrow(dronesAlleles), 1000) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) # Create a Colony class @@ -387,6 +407,7 @@ test_that("getCsdAlleles", { rm(SP) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -402,6 +423,8 @@ test_that("getCsdAlleles", { # test unique and colapse rm(SP) SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 5) + SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -428,6 +451,7 @@ test_that("getCsdGeno", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 5) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -447,6 +471,7 @@ test_that("getCsdGeno", { SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -464,12 +489,15 @@ test_that("getCsdGeno", { # ---- isCsdHeterozygous ---- test_that("isCsdHeterozygous", { - founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) + founderGenomes <- quickHaplo(nInd = 50, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 5) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) + expect_true(all(isCsdHeterozygous(basePop, simParamBee = SP))) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) + expect_true(all(isCsdHeterozygous(drones, simParamBee = SP))) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) # Create a Colony class @@ -486,17 +514,9 @@ test_that("isCsdHeterozygous", { # set CSD to NULL SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L - basePop <- createVirginQueens(founderGenomes, simParamBee = SP) - - drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) - fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) - - # Create a Colony class - colony <- createColony(x = basePop[2], simParamBee = SP) - colony <- cross(colony, drones = fatherGroups[[1]], simParamBee = SP) - colony <- buildUp(x = colony, simParamBee = SP) - expect_error(isCsdHeterozygous(colony@queen, simParamBee = SP)) + basePop <- createVirginQueens(founderGenomes[10:15], simParamBee = SP) + expect_error(isCsdHeterozygous(basePop, simParamBee = SP)) }) # ---- nCsdAlleles ---- @@ -505,6 +525,7 @@ test_that("nCsdAlleles", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -523,6 +544,7 @@ test_that("nCsdAlleles", { # set CSD to NULL SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -538,6 +560,8 @@ test_that("nCsdAlleles", { #collapse argument nCsdAlleles <- 5 SP <- SimParamBee$new(founderGenomes, nCsdAlleles = nCsdAlleles) + SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -562,6 +586,7 @@ test_that("calcBeeGRMIbs", { SP$addTraitA(10) SP$addSnpChip(5) + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) @@ -615,6 +640,7 @@ test_that("editCsdLocus", { founderGenomes <- quickHaplo(nInd = 100, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = 1, nCsdAlleles = 8) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, editCsd = FALSE, simParamBee = SP) nrow(getCsdAlleles(basePop, unique = TRUE, simParamBee = SP)) expect_false(all(isCsdHeterozygous(basePop, simParamBee = SP))) @@ -629,6 +655,8 @@ test_that("editCsdLocus", { test_that("emptyNULL", { founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = 1, nCsdAlleles = 8) + SP$nThreads <- 1L + basePop <- createVirginQueens(founderGenomes, editCsd = FALSE, simParamBee = SP) expect_true(isEmpty(new(Class = "Pop"))) @@ -665,6 +693,7 @@ test_that("isDronesPresent", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -690,6 +719,7 @@ test_that("isFathersPresent", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -718,6 +748,7 @@ test_that("isWorkersPresent", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -743,9 +774,18 @@ test_that("isGenoHeterozygous", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) + #Test on a Pop + (tmp <- getCsdGeno(basePop, simParamBee = SP)) + expect_true(all(SIMplyBee:::isGenoHeterozygous(tmp))) + drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) + # Test on drones + (tmp <- getCsdGeno(drones, simParamBee = SP)) + expect_true(all(SIMplyBee:::isGenoHeterozygous(tmp))) + droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson, simParamBee = SP) # Create a Colony and a MultiColony class @@ -776,6 +816,7 @@ test_that("getBV", { SP$nThreads = 1L SP$addTraitA(nQtlPerChr = 10, var = 1) SP$addSnpChip(5) + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -804,6 +845,7 @@ test_that("getDd", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + SP$addTraitAD(nQtlPerChr = 10, meanDD = 0.2, varDD = 0.1) basePop <- createVirginQueens(founderGenomes, simParamBee = SP) @@ -833,6 +875,7 @@ test_that("getAa", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + SP$addTraitADE(nQtlPerChr = 10, meanDD = 0.2, varDD = 0.1, relAA = 0.5) basePop <- createVirginQueens(founderGenomes, simParamBee = SP) @@ -862,6 +905,7 @@ test_that("editCsdLocus", { founderGenomes <- quickHaplo(nInd = 100, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = 1, nCsdAlleles = 8) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, editCsd = FALSE, simParamBee = SP) nrow(getCsdAlleles(basePop, unique = TRUE, simParamBee = SP)) all(isCsdHeterozygous(basePop, simParamBee = SP)) @@ -891,9 +935,9 @@ test_that("getLocation", { expect_equal(getLocation(apiary, collapse = TRUE), tmp) loc <- c(123, 456) - expect_equal(getLocation(setLocation(colony, location = loc)), loc) + expect_equal(getLocation(setLocation(colony, location = loc, simParamBee = SP)), loc) - expect_equal(getLocation(setLocation(apiary, location = loc)), + expect_equal(getLocation(setLocation(apiary, location = loc, simParamBee = SP)), list("2" = loc, "3" = loc)) }) @@ -902,27 +946,32 @@ test_that("createCrossPlan", { founderGenomes <- quickHaplo(nInd = 1000, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) # Create three virgin MultiColony objects with locations virginColonies1 <- createMultiColony(basePop[1:2], simParamBee = SP) virginColonies1 <- setLocation(virginColonies1, location = Map(c, runif(2, 0, 2*pi), - runif(2, 0, 2*pi))) + runif(2, 0, 2*pi)), + simParamBee = SP) virginColonies2 <- createMultiColony(basePop[3:4], simParamBee = SP) virginColonies2 <- setLocation(virginColonies2, location = Map(c, runif(2, 0, 2*pi), - runif(2, 0, 2*pi))) + runif(2, 0, 2*pi)), + simParamBee = SP) virginColonies3 <- createMultiColony(basePop[5:6], simParamBee = SP) virginColonies3 <- setLocation(virginColonies3, location = Map(c, runif(2, 0, 2*pi), - runif(2, 0, 2*pi))) + runif(2, 0, 2*pi)), + simParamBee = SP) # Create drone colonies droneColonies <- createMultiColony(basePop[7:9], simParamBee = SP) droneColonies <- setLocation(droneColonies, location = Map(c, runif(3, 0, 2*pi), - runif(3, 0, 2*pi))) + runif(3, 0, 2*pi)), + simParamBee = SP) # Create some drones to mate initial drone colonies with DCA <- createDrones(basePop[10:12], nInd = 20, simParamBee = SP) @@ -959,6 +1008,7 @@ test_that("getCaste", { founderGenomes <- quickHaplo(nInd = 1000, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) expect_vector(getCaste(basePop, simParamBee = SP), "virginQueens") @@ -1023,5 +1073,59 @@ test_that("getCaste", { +test_that("getIbdHaplo", { + founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 5) + SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 4) + SP$nThreads = 1L + SP$setTrackRec(isTrackRec = TRUE) + SP$setTrackPed(isTrackPed = TRUE) + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) + baseHaplo = getIbdHaplo(basePop, simParamBee = SP) + expect_equal(nrow(baseHaplo), 2*4) + expect_equal(ncol(baseHaplo), 5) + + drones <- createDrones(x = basePop[1], nInd = 200, simParamBee = SP) + expect_equal(nrow(getIbdHaplo(drones, simParamBee = SP)), 200*1) + droneGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson, simParamBee = SP) + + # Create a Colony and a MultiColony class + colony <- createColony(x = basePop[2], simParamBee = SP) + colony <- cross(colony, drones = droneGroups[[1]], simParamBee = SP) + colony <- buildUp(x = colony, nWorkers = 3, nDrones = 2, simParamBee = SP) + colony <- addVirginQueens(x = colony, nInd = 2, simParamBee = SP) + expect_length(getIbdHaplo(colony, simParamBee = SP), 5) + + apiary <- createMultiColony(basePop[3:4], simParamBee = SP) + apiary <- cross(apiary, drones = droneGroups[c(2, 3)], simParamBee = SP) + apiary <- buildUp(x = apiary, nWorkers = 3, nDrones = 2, simParamBee = SP) + apiary <- addVirginQueens(x = apiary, nInd = 2, simParamBee = SP) + expect_length(getIbdHaplo(apiary, simParamBee = SP), 2) + }) + + +test_that("trackingHomozygotes", { + founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) + SP <- SimParamBee$new(founderGenomes) + SP$nThreads = 1L + SP$setTrackPed(T) + SP$setTrackRec(T) + expect_equal(nrow(SP$pedigree), 0) + expect_equal(length(SP$caste), 0) + expect_equal(length(SP$recHist), 0) + + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) + expect_equal(nrow(SP$pedigree), length(SP$caste)) + expect_equal(nrow(SP$pedigree), length(SP$recHist)) + + drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) + expect_equal(nrow(SP$pedigree), length(SP$caste)) + expect_equal(nrow(SP$pedigree), length(SP$recHist)) + fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) + colony <- createColony(x = basePop[1], simParamBee = SP) + colony <- cross(colony, drones = fatherGroups[[1]], simParamBee = SP) + colony <- buildUp(x = colony, nWorkers = 200, nDrones = 50, simParamBee = SP) + expect_equal(nrow(SP$pedigree), length(SP$caste)) + expect_equal(nrow(SP$pedigree), length(SP$recHist)) +}) diff --git a/tests/testthat/test-L1_pop_functions.R b/tests/testthat/test-L1_pop_functions.R index 7cabd2b8..d73e5420 100644 --- a/tests/testthat/test-L1_pop_functions.R +++ b/tests/testthat/test-L1_pop_functions.R @@ -4,6 +4,7 @@ test_that("getCastePop", { founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 15, simParamBee = SP) @@ -39,6 +40,7 @@ test_that("createVirginQueens", { founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + #check that output is virginqueens ? expect_true(all(isVirginQueen(createVirginQueens(founderGenomes, simParamBee = SP), simParamBee = SP))) @@ -90,6 +92,7 @@ test_that("createDrones", { founderGenomes <- quickHaplo(nInd = 6, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + # Error: x can't be a MapPop expect_error(createDrones(founderGenomes, simParamBee = SP)) @@ -140,6 +143,7 @@ test_that("combineBeeGametes", { founderGenomes <- quickHaplo(nInd = 6, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 15, simParamBee = SP) @@ -171,6 +175,7 @@ test_that("pullCastePop", { founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 15, simParamBee = SP) @@ -227,6 +232,7 @@ test_that("cross", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, nInd = 100, simParamBee = SP) drones <- createDrones(basePop[1], n = 1000, simParamBee = SP) dronesGroups <- pullDroneGroupsFromDCA(drones, n = 7, nDrones = 15, simParamBee = SP) @@ -282,47 +288,17 @@ test_that("cross", { expect_error(cross(colony2, drones = dronesGroups[7], simParamBee = SP)) # Message if fathers == 0 "Mating failed" - expect_error(cross(virginQueen2, drones= selectInd(colony@drones,nInd = 0, use = "rand", simParam = SP), simParamBee = SP)) + expect_error(cross(virginQueen2, drones= selectInd(colony@drones, nInd = 0, use = "rand", simParam = SP), simParamBee = SP)) #expect_message(cross(virginQueen2, drones= selectInd(colony@drones,nInd = 0, use = "rand", simParam = SP), checkCross = "warning", simParamBee = SP)) }) -# ---- setQueensYearOfBirth ---- -test_that("setQueensYearOfBirth", { - founderGenomes <- quickHaplo(nInd = 7, nChr = 1, segSites = 100) - SP <- SimParamBee$new(founderGenomes, csdChr = NULL) - SP$nThreads = 1L - basePop <- createVirginQueens(founderGenomes, nInd = 100, simParamBee = SP) - drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) - dronesGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = nFathersPoisson, simParamBee = SP) - - colony <- createColony(x = basePop[2], simParamBee = SP) - colony <- cross(x = colony, drones = dronesGroups[[1]], simParamBee = SP) - colony <- buildUp(colony, simParamBee = SP) - # Error if x = pop, and not a vq or q - expect_error(setQueensYearOfBirth(colony@workers, simParamBee = SP)) - expect_error(setQueensYearOfBirth(colony@drones, simParamBee = SP)) - - colony <- removeQueen(colony, simParamBee = SP) - # Error if x = colony and no queen is present - expect_error(setQueensYearOfBirth(colony, simParamBee = SP)) - - apiary <- createMultiColony(basePop[3:4], n = 2, simParamBee = SP) - apiary <- cross(apiary, drones = dronesGroups[c(2, 3)], simParamBee = SP) - - colony1 <- createColony(x = basePop[5], simParamBee = SP) - colony1 <- cross(colony1, drones = dronesGroups[[4]], simParamBee = SP) - queen1 <- getQueen(colony1, simParamBee = SP) - - expect_s4_class(setQueensYearOfBirth(queen1, year = 2022, simParamBee = SP), "Pop") - expect_s4_class(setQueensYearOfBirth(colony1, year = 2022, simParamBee = SP), "Colony") - expect_s4_class(setQueensYearOfBirth(apiary, year = 2022, simParamBee = SP), "MultiColony") -}) # ---- createDCA ---- test_that("createDCA", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -358,6 +334,7 @@ test_that("pullDroneGroupsFromDCA", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -388,6 +365,7 @@ test_that("combineBeeGametes", { founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) queen <- basePop[1] @@ -404,6 +382,7 @@ test_that("combineBeeGametesHaploidDiploid", { founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) queen <- basePop[1] @@ -414,3 +393,4 @@ test_that("combineBeeGametesHaploidDiploid", { expect_equal(nInd(drones), 5) expect_equal(drones@ploidy, 1) }) + diff --git a/tests/testthat/test-L2_colony_functions.R b/tests/testthat/test-L2_colony_functions.R index bb1a729e..7154f39e 100644 --- a/tests/testthat/test-L2_colony_functions.R +++ b/tests/testthat/test-L2_colony_functions.R @@ -4,6 +4,7 @@ test_that("createColony", { founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 15, simParamBee = SP) matedQueen <- cross(basePop[2], drones = drones, simParamBee = SP) @@ -24,6 +25,7 @@ test_that("reQueen", { founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 30, simParamBee = SP) virginQueen <- basePop[2] @@ -67,6 +69,7 @@ test_that("Add functions", { founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 100, simParamBee = SP) @@ -109,9 +112,7 @@ test_that("Add functions", { expect_equal(nDrones(addDrones(colony, nInd = 5, new = TRUE, simParamBee = SP), simParamBee = SP), 5) # If input is an apiary # Empty apiary - you can add, but nothing happens - returns an empty apiary - expect_s4_class(addVirginQueens(emptyApiary, nInd = 5, simParamBee = SP), "MultiColony") - expect_s4_class(addWorkers(emptyApiary, nInd = 5, simParamBee = SP), "MultiColony") - expect_s4_class(addDrones(emptyApiary, nInd = 5, simParamBee = SP), "MultiColony") + expect_error(addVirginQueens(emptyApiary, nInd = 5, simParamBee = SP)) # Non-empty apiary expect_s4_class(addVirginQueens(apiary, nInd = 5, simParamBee = SP), "MultiColony") expect_s4_class(addWorkers(apiary, nInd = 5, simParamBee = SP), "MultiColony") @@ -124,6 +125,7 @@ test_that("BuildUpDownsize", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP) @@ -154,7 +156,7 @@ test_that("BuildUpDownsize", { # Build Up an apiary # Empty apiary - expect_s4_class(buildUp(emptyApiary, simParamBee = SP), "MultiColony") + expect_error(buildUp(emptyApiary, simParamBee = SP)) # Non-empty apiary expect_equal(nColonies(buildUp(apiary, simParamBee = SP)), 2) @@ -169,7 +171,7 @@ test_that("BuildUpDownsize", { expect_length(intersect(getId(getWorkers(downsize(colony, p = 0.1, new = TRUE, simParamBee = SP), simParamBee = SP)), workersIDs), 0) # Empty apiary - expect_s4_class(downsize(emptyApiary, simParamBee = SP), "MultiColony") + expect_error(downsize(emptyApiary, simParamBee = SP)) # Non-empty apiary downsize(apiary, simParamBee = SP) }) @@ -180,6 +182,7 @@ test_that("replaceFunctions", { founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 100, simParamBee = SP) @@ -201,9 +204,9 @@ test_that("replaceFunctions", { expect_error(replaceVirginQueens(emptyColony, p = 0.5, simParamBee = SP)) expect_error(replaceWorkers(emptyColony, p = 0, simParamBee = SP)) expect_error(replaceDrones(emptyColony, p = 1, simParamBee = SP)) - expect_s4_class(replaceVirginQueens(emptyApiary, p = 0.5, simParamBee = SP), "MultiColony") - expect_s4_class(replaceWorkers(emptyApiary, p = 0, simParamBee = SP), "MultiColony") - expect_s4_class(replaceDrones(emptyApiary, p = 1, simParamBee = SP), "MultiColony") + expect_error(replaceVirginQueens(emptyApiary, p = 0.5, simParamBee = SP)) + expect_error(replaceWorkers(emptyApiary, p = 0, simParamBee = SP)) + expect_error(replaceDrones(emptyApiary, p = 1, simParamBee = SP)) # Replace individuals in the non-empty colony/apiary expect_s4_class(replaceVirginQueens(colony, simParamBee = SP), "Colony") @@ -211,7 +214,7 @@ test_that("replaceFunctions", { expect_s4_class(replaceDrones(colony, simParamBee = SP), "Colony") expect_equal(nVirginQueens(replaceVirginQueens(colony, p = 1, simParamBee = SP), simParamBee = SP), nVirginQueens(colony, simParam = SP)) expect_equal(nWorkers(replaceWorkers(colony, p = 0.5, simParamBee = SP), simParamBee = SP), nWorkers(colony, simParamBee = SP)) - expect_equal(nDrones(replaceDrones(colony, p = 0, simParamBee = SP), simParamBee = SP), nDrones(colony, simParamBee = SP)) + expect_warning(nDrones(replaceDrones(colony, p = 0, simParamBee = SP), simParamBee = SP)) virginQueensIDs <- getId(colony@virginQueens) workerIDs <- getId(colony@workers) droneIDs <- getId(colony@drones) @@ -219,14 +222,14 @@ test_that("replaceFunctions", { virginQueensIDs), 0) expect_length(intersect(getId(replaceWorkers(colony, p = 0.5, simParamBee = SP)@workers), workerIDs), nWorkers(colony, simParamBee = SP)/2) - expect_length(intersect(getId(replaceDrones(colony, p = 0, simParamBee = SP)@drones), - droneIDs), nDrones(colony, simParamBee = SP)) + expect_warning(intersect(getId(replaceDrones(colony, p = 0, simParamBee = SP)@drones), + droneIDs)) expect_s4_class(replaceVirginQueens(apiary, simParamBee = SP), "MultiColony") expect_s4_class(replaceWorkers(apiary, simParamBee = SP), "MultiColony") expect_s4_class(replaceDrones(apiary, simParamBee = SP), "MultiColony") expect_equal(nColonies(replaceVirginQueens(apiary, p = 1, simParamBee = SP)), nColonies(apiary)) expect_equal(nColonies(replaceWorkers(apiary, p = 0.5, simParamBee = SP)), nColonies(apiary)) - expect_equal(nColonies(replaceDrones(apiary, p = 0, simParamBee = SP)), nColonies(apiary)) + expect_error(nColonies(replaceDrones(apiary, p = 0, simParamBee = SP))) }) # ---- Remove functions ---- @@ -235,6 +238,7 @@ test_that("removeFunctions", { founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 100, simParamBee = SP) @@ -256,9 +260,9 @@ test_that("removeFunctions", { expect_s4_class(removeVirginQueens(emptyColony, p = 0.5, simParamBee = SP), "Colony") expect_s4_class(removeWorkers(emptyColony, p = 0, simParamBee = SP), "Colony") expect_s4_class(removeDrones(emptyColony, p = 1, simParamBee = SP), "Colony") - expect_s4_class(removeVirginQueens(emptyApiary, p = 0.5, simParamBee = SP), "MultiColony") - expect_s4_class(removeWorkers(emptyApiary, p = 0, simParamBee = SP), "MultiColony") - expect_s4_class(removeDrones(emptyApiary, p = 1, simParamBee = SP), "MultiColony") + expect_error(removeVirginQueens(emptyApiary, p = 0.5, simParamBee = SP)) + expect_error(removeWorkers(emptyApiary, p = 0, simParamBee = SP)) + expect_error(removeDrones(emptyApiary, p = 1, simParamBee = SP)) # Remove individuals in the non-empty colony/apiary expect_s4_class(removeVirginQueens(colony, simParamBee = SP), "Colony") @@ -281,6 +285,7 @@ test_that("setLocation", { founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 50) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 1000, simParamBee = SP) droneGroups <- pullDroneGroupsFromDCA(drones, n = 4, nDrones = 10, simParamBee = SP) @@ -290,29 +295,28 @@ test_that("setLocation", { apiary <- cross(apiary, drones = droneGroups[2:4], simParamBee = SP) loc <- c(1, 1) - expect_equal(getLocation(setLocation(colony, location = loc)), loc) + expect_equal(getLocation(setLocation(colony, location = loc, simParamBee = SP)), loc) - expect_equal(getLocation(setLocation(apiary, location = loc)), + expect_equal(getLocation(setLocation(apiary, location = loc, simParamBee = SP)), list("2" = loc, "3" = loc, "4" = loc)) locList <- list("2" = c(0, 0), "3" = c(1, 1), "4" = c(2, 2)) - expect_equal(getLocation(setLocation(apiary, location = locList)), locList) + expect_equal(getLocation(setLocation(apiary, location = locList, simParamBee = SP)), locList) locDF <- data.frame(x = c(0, 1, 2), y = c(0, 1, 2)) - expect_equal(getLocation(setLocation(apiary, location = locDF)), locList) + expect_equal(getLocation(setLocation(apiary, location = locDF, simParamBee = SP)), locList) emptyColony <- createColony(simParamBee = SP) - expect_s4_class(setLocation(emptyColony, location = c(1,1)), "Colony") - expect_equal(setLocation(emptyColony, location = c(1,1))@location, c(1,1)) + expect_s4_class(setLocation(emptyColony, location = c(1,1), simParamBee = SP), "Colony") + expect_equal(setLocation(emptyColony, location = c(1,1), simParamBee = SP)@location, c(1,1)) emptyApiary <- createMultiColony(n = 3, simParamBee = SP) apiary <- createMultiColony(basePop[1:3], simParamBee = SP) - expect_s4_class(setLocation(emptyApiary, location = c(1,2)), "MultiColony") - expect_error(setLocation(emptyApiary, location = list(1,2))) # Lengths do not match - expect_s4_class(setLocation(emptyApiary, location = list(1:2, 3:4, 4:5)), "MultiColony") #Not setting anything, if all are NULL!!!! - expect_s4_class(setLocation(apiary, location = c(1,2)), "MultiColony") - expect_s4_class(setLocation(apiary, location = list(1:2, 3:4, 4:5)), "MultiColony") + expect_error(setLocation(emptyApiary, location = c(1,2), simParamBee = SP)) + expect_error(setLocation(emptyApiary, location = list(1,2), simParamBee = SP)) # Lengths do not match + expect_s4_class(setLocation(apiary, location = c(1,2), simParamBee = SP), "MultiColony") + expect_s4_class(setLocation(apiary, location = list(1:2, 3:4, 4:5), simParamBee = SP), "MultiColony") }) # ---- Supersede ---- @@ -321,6 +325,7 @@ test_that("supersede", { founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 1000, simParamBee = SP) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) @@ -357,6 +362,7 @@ test_that("split", { founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 1000, simParamBee = SP) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) @@ -398,6 +404,7 @@ test_that("resetEvents", { founderGenomes <- quickHaplo(nInd = 5, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1], nInd = 100, simParamBee = SP) @@ -421,8 +428,8 @@ test_that("resetEvents", { expect_true(hasSuperseded(colony)) expect_true(all(hasSuperseded(apiary))) - colony <- resetEvents(colony) - apiary <- resetEvents(apiary) + colony <- resetEvents(colony, simParamBee = SP) + apiary <- resetEvents(apiary, simParamBee = SP) expect_false(isProductive(colony)) expect_false(all(isProductive(apiary))) @@ -436,6 +443,7 @@ test_that("Combine", { founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 1000, simParamBee = SP) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) @@ -456,14 +464,14 @@ test_that("Combine", { apiary1 <- buildUp(x = apiary1, nWorkers = 100, nDrones = 20, simParamBee = SP) apiary2 <- buildUp(x = apiary2, nWorkers = 20, nDrones = 5, simParamBee = SP) - colony3 <- combine(strong = colony1, weak = colony2) - apiary3 <- combine(strong = apiary1, weak = apiary2) + colony3 <- combine(strong = colony1, weak = colony2, simParamBee = SP) + apiary3 <- combine(strong = apiary1, weak = apiary2, simParamBee = SP) expect_equal(nWorkers(colony3, simParamBee = SP),sum(nWorkers(colony1, simParamBee = SP), nWorkers(colony2, simParamBee = SP))) expect_equal(colony1@queen@id, colony3@queen@id) expect_equal(nWorkers(apiary3[[2]], simParamBee = SP),sum(nWorkers(apiary1[[2]], simParamBee = SP), nWorkers(apiary2[[2]], simParamBee = SP))) colony1 <- NULL colony2 <- NULL - expect_error(combine(strong = colony1, weak = colony2)) # discus the output + expect_error(combine(strong = colony1, weak = colony2, simParamBee = SP)) # discus the output }) # ---- Swarm ---- @@ -472,6 +480,7 @@ test_that("swarm", { founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 1000, simParamBee = SP) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) @@ -514,6 +523,7 @@ test_that("collapse", { founderGenomes <- quickHaplo(nInd = 10, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 1000, simParamBee = SP) fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP) @@ -526,12 +536,12 @@ test_that("collapse", { # Collapse expect_false(hasCollapsed(colony)) - colony <- collapse(colony) + colony <- collapse(colony, simParamBee = SP) expect_true(hasCollapsed(colony)) expect_false(all(hasCollapsed(apiary))) tmp <- pullColonies(apiary, n = 2, simParamBee = SP) - apiaryLost <- collapse(tmp$pulled) + apiaryLost <- collapse(tmp$pulled, simParamBee = SP) expect_true(all(hasCollapsed(apiaryLost))) apiaryLeft <- tmp$remnant expect_false(all(hasCollapsed(apiaryLeft))) diff --git a/tests/testthat/test-L3_Colonies_functions.R b/tests/testthat/test-L3_Colonies_functions.R index e73a4b8f..2598dc2b 100644 --- a/tests/testthat/test-L3_Colonies_functions.R +++ b/tests/testthat/test-L3_Colonies_functions.R @@ -1,11 +1,11 @@ # Level 3 MultiColony Functions - # ---- createMultiColony ---- test_that("createMultiColony", { founderGenomes <- quickHaplo(nInd = 6, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(basePop[1], n = 100, simParamBee = SP) # Error if individuals x are not vq or q @@ -37,12 +37,13 @@ test_that("createMultiColony", { expect_s4_class(createMultiColony(x = basePop[4:5], n = 2, simParamBee = SP), "MultiColony") }) -# ---- selectColonies ---- +# ---- selectColonies --- test_that("selectColonies", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) drones <- createDrones(x = basePop[1:4], nInd = 100, simParamBee = SP) @@ -87,6 +88,7 @@ test_that("pullColonies", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) # Error if argument multicolony isn't a multicolony class expect_error(pullColonies(basePop, simParamBee = SP)) @@ -128,6 +130,7 @@ test_that("removeColonies", { founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100) SP <- SimParamBee$new(founderGenomes, csdChr = NULL) SP$nThreads = 1L + basePop <- createVirginQueens(founderGenomes, simParamBee = SP) # Error if argument multicolony isn't a multicolony class expect_error(removeColonies(basePop, simParamBee = SP)) diff --git a/vignettes/Colony_locations.csv b/vignettes/Colony_locations.csv index 0ac2a90d..6678dcde 100644 --- a/vignettes/Colony_locations.csv +++ b/vignettes/Colony_locations.csv @@ -1,36 +1,36 @@ ColonyID,X,Y -1,0.662431162288274,5.97033145745812 -2,0.889869211813095,0.882040306233467 -3,3.20406617225118,4.2012594475023 -4,3.37467634975274,1.85276144978548 -5,6.26877271726022,4.12336288238627 -6,1.49762073729787,0.854711175179748 -7,6.27428327657999,5.28537368571472 -8,0.377119552748809,1.26003243402567 -9,0.884600259265786,2.55843304434275 -10,4.85341262281402,4.34423864421118 -11,0.439273950460384,5.78768883580839 -12,4.85267013629791,5.24037990077726 -13,6.27814888107222,1.67684867115787 -14,5.91398658831959,2.21947012261649 +1,0.211675140867833,1.58970616827141 +2,0.377119552748809,1.26003243402567 +3,0.439273950460384,5.78768883580839 +4,0.662431162288274,5.97033145745812 +5,0.884600259265786,2.55843304434275 +6,0.889869211813095,0.882040306233467 +7,1.49762073729787,0.854711175179748 +8,1.54126706057672,0.265964466739025 +9,1.59318554922445,3.95724726174676 +10,1.69897541097397,2.4435374157815 +11,1.76640287495849,1.81689439235 +12,1.94181870616625,1.04070091389605 +13,1.96273450676472,2.98552319129881 +14,2.15001173188715,5.30559199476844 15,2.2845571049277,2.76273156562477 -16,2.15001173188715,5.30559199476844 -17,3.30277055998226,3.88408253149063 -18,1.59318554922445,3.95724726174676 -19,5.14489315015939,3.48380219722517 -20,4.89542592867685,4.87175443368121 -21,4.98504294579104,4.63186113766538 -22,1.96273450676472,2.98552319129881 -23,1.94181870616625,1.04070091389605 -24,3.71355474699821,3.98892629339701 -25,1.76640287495849,1.81689439235 -26,3.49162610986539,2.007127614613 -27,4.70110619836582,1.98065883153337 -28,2.93773502070683,2.79053982429322 -29,1.69897541097397,2.4435374157815 -30,1.54126706057672,0.265964466739025 -31,0.211675140867833,1.58970616827141 -32,4.38863010920245,4.35616019770602 -33,4.3632705003701,0.955920230806015 -34,5.94574863325625,5.50420647366442 -35,2.86914251070775,0.176914088999066 +16,2.86914251070775,0.176914088999066 +17,2.93773502070683,2.79053982429322 +18,3.20406617225118,4.2012594475023 +19,3.30277055998226,3.88408253149063 +20,3.37467634975274,1.85276144978548 +21,3.49162610986539,2.007127614613 +22,3.71355474699821,3.98892629339701 +23,4.3632705003701,0.955920230806015 +24,4.38863010920245,4.35616019770602 +25,4.70110619836582,1.98065883153337 +26,4.85267013629791,5.24037990077726 +27,4.85341262281402,4.34423864421118 +28,4.89542592867685,4.87175443368121 +29,4.98504294579104,4.63186113766538 +30,5.14489315015939,3.48380219722517 +31,5.91398658831959,2.21947012261649 +32,5.94574863325625,5.50420647366442 +33,6.26877271726022,4.12336288238627 +34,6.27428327657999,5.28537368571472 +35,6.27814888107222,1.67684867115787 diff --git a/vignettes/H_Parallelisation.Rmd b/vignettes/H_Parallelisation.Rmd new file mode 100644 index 00000000..d967704b --- /dev/null +++ b/vignettes/H_Parallelisation.Rmd @@ -0,0 +1,162 @@ +--- +title: "Parallelisation and high-performing cluster setup" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Multiple colonies} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +editor_options: + markdown: + wrap: 80 + canonical: true +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + include = TRUE +) +``` + +# Quick set-up instructions + +Here, we show how you should set-up the parallel back-end on different +environments. We do recommend reading the remaining of this vignette. We +recommend running these lines straight after setting the `SimParamBee`. + +```{r quick_setup, eval=F, echo=T} +library(SIMplyBee) +library(parallel) +library(doParallel) + +founderGenomes <- quickHaplo(nInd = 30, nChr = 1, segSites = 100) +SP <- SimParamBee$new(founderGenomes) +SP$nThreads <- NCORES #Where NCORES is a specified number or all available cores (detectCores(), see below) + +# If using Linux/MACOS +registerDoParallel(cores = SP$nThreads) + +# If using Windows machine / running the simulation on HPC +cl <- makeCluster(SP$nThreads, type="PSOCK") +registerDoParallel(cl) +#Do the simulation +# At the end of everything you run +stopImplicitCluster() +``` + +# Introduction + +Honeybee simulations consist of simulating individuals as `Pop` classes, and +then on top of this also `Colony` and `MultiColony` classes, all of them with +their meta-data. This makes the simulation computationally demanding and slow. +With the aim to speed up the simulation, we parallelised the major functions +with `foreach` and `doParallel` R packages. Nothing changed in terms of running +the functions, but do functions now have the ability to run on multiple cores at +the same time. They would all search for the number of available cores in the +`SimParamBee` object, under `nThreads`. You can set that to a desired number or +make the simulation use all available cores. + +```{r nThread_setup} +library(package = "SIMplyBee") +library(package = "parallel") + +founderGenomes <- quickHaplo(nInd = 30, nChr = 1, segSites = 100) +SP <- SimParamBee$new(founderGenomes) + +# Set the number of cores to use +SP$nThreads <- 8 +# Or use all available cores +SP$nThreads <- detectCores() +``` + +In R, there are two possible options for parallelisation, `FORK` and `PSOCK`. +The forking (`FORK`) creates subprocesses that share memory and objects with the +parent process, which means it's very efficient! However, such parallelisation +is not supported on Windows machines and is not allowed on most high-performing +clusters! + +The alternative, PSOCK, can be used on all system, but it works by creating a +separate R process for each subprocess (that communicate through sockets), +meaning that the whole environment needs to be copied to each subprocess, +creating a larger memory overhead. (adapted from +). + +We profiled the following piece of code using different parallelisation options. + +```{r defining_testing_function} +create_bee_colonies <- function() { + founderGenomes <- quickHaplo(nInd = 200, nChr = 1, segSites = 50) + SP <- SimParamBee$new(founderGenomes) + + basePop <- createVirginQueens(founderGenomes) + drones = createDrones(basePop, nInd = 100) + baseColonies <- createMultiColony(cross(basePop, drones = drones, crossPlan = "create")) + baseColonies <- buildUp(baseColonies, nWorkers = 1000, nDrones = 1000) + baseColonies <- supersede(baseColonies) + baseColonies <- cross(baseColonies, drones = drones, crossPlan = "create") + tmp = split(baseColonies) +} +``` + +We set up different parallelisation back-ends with the following code, where +`ncores` was either 1 or 8. + +```{r parallelisation_options, eval=F, echo=T} + +# First one - ??? +SP$nThreads = ncores +registerDoParallel(cores = SP$nThreads) +create_bee_colonies() + +# Second one - create a FORK cluster +SP$nThreads = ncores +cl <- makeCluster(SP$nThreads, type="FORK") +registerDoParallel(cl) +create_bee_colonies() +stopImplicitCluster() + +# Third one - create a PSOCK cluster +SP$nThreads = ncores +cl <- makeCluster(SP$nThreads, type="PSOCK") +registerDoParallel(cl) +create_bee_colonies() +stopImplicitCluster() +``` + +Here are the results of running these different options. You can see that the +time can be significantly improved when running on multiple cores. When running +on a single core, setting up the parallelisation cluster via `FORK` or `PSOCK` +actually adds some overhead time. + +```{r meanTime_figure, echo=FALSE, out.width='100%'} +knitr::include_graphics("Time_mean.png") +``` + +```{r meanPCPU_figure, echo=FALSE, out.width='100%'} +knitr::include_graphics("PCPU_mean.png") +``` + +```{r meanRSS_figure, echo=FALSE, out.width='100%'} +knitr::include_graphics("RSS_mean.png") +``` + +The functions that are currently explicitily parallelised: + +- L1: cross, pullCasteP, createCastePop + +- L2: supersede, split, swarm, setEvents, combine, setLocation, reQueen, + addCastePop, removeCastePop, resetEvents, collapse + +- L3: createMultiColony + +The following figure shows the benefit of parallelised (p, in blue) vs +sequential/non-parallelised functions (np, in red) when run on a Linux machine +with 16 cores. The figure shows the mean and the standard deviation across 10 +replicates. + +```{r functions_time, echo=FALSE, out.width='100%'} +knitr::include_graphics("Profiling_parallelised_functions_Unix.png") +``` + diff --git a/vignettes/PCPU_mean.png b/vignettes/PCPU_mean.png new file mode 100644 index 00000000..5b3949aa Binary files /dev/null and b/vignettes/PCPU_mean.png differ diff --git a/vignettes/RSS_mean.png b/vignettes/RSS_mean.png new file mode 100644 index 00000000..2d2f5975 Binary files /dev/null and b/vignettes/RSS_mean.png differ diff --git a/vignettes/Time_mean.png b/vignettes/Time_mean.png new file mode 100644 index 00000000..7bad43bc Binary files /dev/null and b/vignettes/Time_mean.png differ