Merge pull request #207 from gaynorr/devel

CRAN release 1.6.1
gaynorr · Nov 1, 2024 · a171cf6 · a171cf6
2 parents 832c1e4 + fd654bd
commit a171cf6
Show file tree

Hide file tree

Showing 10 changed files with 111 additions and 37 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: AlphaSimR
 Type: Package
 Title: Breeding Program Simulations
-Version: 1.6.0
-Date: 2024-08-15
+Version: 1.6.1
+Date: 2024-11-01
 Authors@R: c(person("Chris", "Gaynor", email = "[email protected]",
   role = c("aut", "cre"), comment = c(ORCID = "0000-0003-0558-6656")),
   person("Gregor", "Gorjanc", role = "ctb", 

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,7 @@
+# AlphaSimR 1.6.1
+
+*fixed bug in `mergePops` and `[` (subset) methods - they were failing for populations that had a misc slot with a matrix - we now check if a misc slot element is a matrix and rbind them for `mergePops` and subset rows for `[` (assuming the first dimension represents individuals)
+
 # AlphaSimR 1.6.0
 
 *exported `meanEBV` and added `varEBV` to complement `meanP`/`varP` and `meanG`/`varG`

diff --git a/R/Class-Pop.R b/R/Class-Pop.R
@@ -377,7 +377,8 @@ isNamedMapPop = function(x) {
 #' Used by genomic selection models but otherwise ignored.
 #' @slot misc a list whose elements correspond to additional miscellaneous
 #' nodes with the items for individuals in the population (see example in
-#' \code{\link{newPop}}).
+#' \code{\link{newPop}}) - we support vectors and matrices or objects that
+#' have a generic length and subset method.
 #' This list is normally empty and exists solely as an
 #' open slot available for uses to store extra information about
 #' individuals.
@@ -462,8 +463,16 @@ setValidity("Pop",function(object){
   if(object@nInd!=length(object@fixEff)){
     errors = c(errors,"nInd!=length(fixEff)")
   }
-  if(any(object@nInd!=sapply(object@misc, length))){
-    errors = c(errors,"any(nInd!=sapply(misc, length))")
+  length2 = function(x){
+    if(is.matrix(x)){
+      ret = dim(x)[1]
+    }else{
+      ret = length(x)
+    }
+    return(ret)
+  }
+  if(any(object@nInd!=sapply(object@misc, length2))){
+    errors = c(errors,"any(nInd!=sapply(misc, length(x) or dim(x)[1]))")
   }
   if(length(errors)==0){
     return(TRUE)
@@ -489,12 +498,19 @@ setMethod("[",
                 stop("Trying to select invalid individuals")
               }
             }
+            subset2 = function(z){
+              if(is.matrix(z)){
+                return(z[i,,drop=FALSE])
+              }else{
+                return(z[i])
+              }
+            }
             x@id = x@id[i]
             x@iid = x@iid[i]
             x@mother = x@mother[i]
             x@father = x@father[i]
             x@fixEff = x@fixEff[i]
-            x@misc = lapply(x@misc, FUN = function(z) z[i])
+            x@misc = lapply(x@misc, FUN = subset2)
             x@miscPop = list()
             x@gv = x@gv[i,,drop=FALSE]
             x@pheno = x@pheno[i,,drop=FALSE]
@@ -796,7 +812,7 @@ resetPop = function(pop,simParam=NULL){
   pop@gv = matrix(NA_real_,nrow=pop@nInd,
                   ncol=simParam$nTraits)
   pop@fixEff = rep(1L,pop@nInd)
-  
+
   # Calculate genetic values
   for(i in seq_len(simParam$nTraits)){
     tmp = getGv(simParam$traits[[i]],pop,simParam$nThreads)
@@ -805,7 +821,7 @@ resetPop = function(pop,simParam=NULL){
       pop@gxe[[i]] = tmp[[2]]
     }
   }
-  
+
   # Add back trait names
   colnames(pop@pheno) = colnames(pop@gv) = traitNames
 
@@ -875,12 +891,12 @@ newEmptyPop = function(ploidy=2L, simParam=NULL){
                     ncol = simParam$nTraits)
 
   traitNames = character(simParam$nTraits)
-  
+
   # Get trait names
   for(i in seq_len(simParam$nTraits)){
     traitNames[i] = simParam$traits[[i]]@name
   }
-  
+
   colnames(traitMat) = traitNames
 
   # Create empty geno list

diff --git a/R/mergePops.R b/R/mergePops.R
@@ -43,52 +43,52 @@ mergePops = function(popList){
     classes = classes[-remove]
   }
   stopifnot(all(classes=="Pop"))
-  
+
   #nChr
   nChr = do.call("c",lapply(popList,
                             function(x) x@nChr))
   stopifnot(all(nChr==nChr[1]))
   nChr = nChr[1]
-  
+
   #ploidy
   ploidy = do.call("c",lapply(popList,
                               function(x) x@ploidy))
   stopifnot(all(ploidy==ploidy[1]))
   ploidy = ploidy[1]
-  
+
   #nLoci
   nLoci = do.call("c",lapply(popList,
                              function(x){
                                all(x@nLoci==popList[[1]]@nLoci)
                              }))
   stopifnot(all(nLoci))
   nLoci = popList[[1]]@nLoci
-  
+
   #id
   id = do.call("c",
                lapply(popList,
                       function(x) x@id))
-  
+
   #iid
   iid = do.call("c",
                 lapply(popList,
                        function(x) x@iid))
-  
+
   #mother
   mother = do.call("c",
                    lapply(popList,
                           function(x) x@mother))
-  
+
   #father
   father= do.call("c",
                   lapply(popList,
                          function(x) x@father))
-  
+
   #fixEff
   fixEff= do.call("c",
                   lapply(popList,
                          function(x) x@fixEff))
-  
+
   #misc
   tmp = sapply(popList, function(x) length(x@misc))
   if(all(tmp == tmp[1]) & tmp[1]>0) {
@@ -106,7 +106,11 @@ mergePops = function(popList){
       misc = vector("list", length=length(tmp[[1]]))
       for(i in seq_len(length(tmp[[1]]))){
         miscTmp = lapply(popList, function(x) x@misc[[i]])
-        misc[[i]] = do.call("c", miscTmp)
+        if (is.matrix(miscTmp[[1]])) {
+          misc[[i]] = do.call("rbind", miscTmp)
+        } else {
+          misc[[i]] = do.call("c", miscTmp)
+        }
       }
       names(misc) = tmp[[1]]
     }else{
@@ -115,30 +119,30 @@ mergePops = function(popList){
   } else {
     misc = list()
   }
-  
+
   #sex
   sex = do.call("c",
                    lapply(popList,
                           function(x) x@sex))
-  
+
   #nTraits
   nTraits = do.call("c",lapply(popList,
                                function(x) x@nTraits))
   stopifnot(all(nTraits==nTraits[1]))
   nTraits = nTraits[1]
-  
+
   #nInd
   nInd = do.call("c",lapply(popList,
                             function(x) x@nInd))
-  
+
   #gv
   gv = do.call("rbind",lapply(popList,
                               function(x) x@gv))
-  
+
   #pheno
   pheno = do.call("rbind",lapply(popList,
                                  function(x) x@pheno))
-  
+
   #ebv
   ebv = do.call("c",lapply(popList,
                            function(x) ncol(x@ebv)))
@@ -148,7 +152,7 @@ mergePops = function(popList){
   }else{
     ebv = matrix(NA_real_,nrow=sum(nInd),ncol=0)
   }
-  
+
   #gxe
   if(nTraits>=1){
     gxe = vector("list",length=nTraits)
@@ -162,12 +166,12 @@ mergePops = function(popList){
   }else{
     gxe = list()
   }
-  
+
   #geno
   nBin = as.integer(nLoci%/%8L + (nLoci%%8L > 0L))
   geno = mergeMultGeno(popList,nInd=nInd,nBin=nBin,ploidy=ploidy)
   dim(geno) = NULL # Account for matrix bug in RcppArmadillo
-  
+
   #wrap it all up into a Pop
   nInd = sum(nInd)
   return(new("Pop",

diff --git a/R/pullGeno.R b/R/pullGeno.R
@@ -43,7 +43,20 @@ getLociNames = function(lociPerChr, lociLoc, genMap){
   return(lociNames)
 }
 
-# Finds loci on a genetic map and return a list of positions
+#' Finds positions of loci by marker name
+#' 
+#' @description Used to generate lociPerChr and lociLoc 
+#' objects for a set of markers. These objects can be passed 
+#' other functions for pulling genotypes or haplotypes.
+#' 
+#' @param markers a vector of marker names
+#' @param genMap a genetic map in AlphaSimR's internal 
+#' genetic map format
+#'
+#' @return A list containing lociPerChr and lociLoc 
+#' that can be
+#' 
+#' @keywords internal
 mapLoci = function(markers, genMap){
   # Check that the markers are present on the map
   genMapMarkerNames = unlist(lapply(genMap, names))

diff --git a/man/Pop-class.Rd b/man/Pop-class.Rd
diff --git a/man/mapLoci.Rd b/man/mapLoci.Rd
diff --git a/src/algorithm.cpp b/src/algorithm.cpp
@@ -321,7 +321,6 @@ void GraphBuilder::markCurrentTree(){
     pTreeEdgesToCoalesceArray[i]=node->getTopEdgeByIndex(0);
   }
   unsigned int iIndex = 0;
-  int iIterations = 0;
   bool bFirstSample = true;
   while(iIndex<iTotalSamples){
     EdgePtr & curEdge = pTreeEdgesToCoalesceArray[iIndex];
@@ -332,7 +331,6 @@ void GraphBuilder::markCurrentTree(){
       iIndex = 0;
     }
     bFirstSample = false;
-    ++iIterations;
   }
 }
 

diff --git a/src/simulator.cpp b/src/simulator.cpp
@@ -169,7 +169,6 @@ void Simulator::readInputParameters(CommandArguments arguments){
         if (inFile.is_open()) {
           pConfig->pAlleleFreqBinPtrSet = new AlleleFreqBinPtrSet;
           string line;
-          int total=0;
           double lastStart=0.;
           double cumFreq=0.;
           double maxFreq = flipAllele?0.5:1.;
@@ -183,7 +182,6 @@ void Simulator::readInputParameters(CommandArguments arguments){
             AlleleFreqBinPtr bin = AlleleFreqBinPtr(new AlleleFreqBin(start,end,freq));
             pConfig->pAlleleFreqBinPtrSet->insert(bin);
             lastStart = end;
-            ++total;
           }
           inFile.close();
           pConfig->bSNPAscertainment = true;
@@ -230,14 +228,12 @@ void Simulator::readInputParameters(CommandArguments arguments){
         if (inFile.is_open()) {
           pConfig->pHotSpotBinPtrList = new HotSpotBinPtrList;
           string line;
-          int total=0;
           while(getline(inFile,line)){
             istringstream inStr(line);
             double start,end,ratio;
             inStr>>start>>end>>ratio;
             HotSpotBinPtr bin(new HotSpotBin(start,end,ratio));
             pConfig->pHotSpotBinPtrList->push_back(bin);
-            ++total;
           }
           inFile.close();
           pConfig->bVariableRecomb = true;

diff --git a/tests/testthat/test-misc.R b/tests/testthat/test-misc.R
@@ -7,36 +7,54 @@ test_that("misc_and_miscPop",{
   popOrig = newPop(founderPop, simParam=SP)
   multiPop = newMultiPop(popOrig, popOrig)
 
+  expect_equal(popOrig@misc, list())
+  expect_equal(popOrig@miscPop, list())
+
   popSub = popOrig[1]
   expect_equal(popSub@misc, list())
   expect_equal(popSub@miscPop, list())
 
   pop = popOrig
   pop@misc$vec = rnorm(n=2)
+  pop@misc$mat = matrix(1:4, nrow=2)
   pop@misc$mtP = popOrig # hmm, should this actually be multiple pop objects or one with multiple individuals?
   pop@misc$mtLP = list(popOrig, popOrig)
   pop@misc$mtMP = multiPop
+  # setting these miscPop elements just as an example - we are not testing them below,
+  # because they get dropped in most/all operations on a pop
   pop@miscPop$vec = sum(pop@misc$vec)
   pop@miscPop$af = colMeans(pullSegSiteGeno(pop, simParam=SP))
+  pop@miscPop$mat = matrix(1:4, nrow=2)
 
   popSub = pop[1]
   expect_equal(popSub@misc$vec, pop@misc$vec[1])
+  expect_equal(popSub@misc$mat, pop@misc$mat[1, , drop=FALSE])
   expect_equal(popSub@misc$mtP, pop@misc$mtP[1])
   expect_equal(popSub@misc$mtLP, pop@misc$mtLP[1])
   expect_equal(popSub@misc$mtMP, pop@misc$mtMP[1])
   expect_equal(popSub@miscPop, list())
 
   popSub = pop[0]
   expect_equal(popSub@misc$vec, numeric(0))
+  expect_equal(popSub@misc$mat, pop@misc$mat[0, , drop=FALSE])
   expect_equal(popSub@misc$mtP, newEmptyPop(simParam=SP))
   expect_equal(popSub@misc$mtLP, list())
   expect_equal(popSub@misc$mtMP, new("MultiPop", pops=list()))
   expect_equal(popSub@miscPop, list())
 
   popC = c(pop, pop)
   expect_equal(popC@misc$vec, c(pop@misc$vec, pop@misc$vec))
+  expect_equal(popC@misc$mat, rbind(pop@misc$mat, pop@misc$mat))
   expect_equal(popC@misc$mtP, c(pop@misc$mtP, pop@misc$mtP))
   expect_equal(popC@misc$mtLP, c(pop@misc$mtLP, pop@misc$mtLP))
   expect_equal(popC@misc$mtMP, c(pop@misc$mtMP, pop@misc$mtMP))
   expect_equal(popC@miscPop, list())
+
+  popC = c(pop, pop[1])
+  expect_equal(popC@misc$vec, c(pop@misc$vec, pop@misc$vec[1]))
+  expect_equal(popC@misc$mat, rbind(pop@misc$mat, pop@misc$mat[1, ]))
+  expect_equal(popC@misc$mtP, c(pop@misc$mtP, pop@misc$mtP[1]))
+  expect_equal(popC@misc$mtLP, c(pop@misc$mtLP, pop@misc$mtLP[1]))
+  expect_equal(popC@misc$mtMP, c(pop@misc$mtMP, pop@misc$mtMP[1]))
+  expect_equal(popC@miscPop, list())
 })