Merge pull request #105 from gaynorr/devel

Version 1.3.4

.Rbuildignore

@@ -1,6 +1,4 @@
-revdep
-vignettes/web_only/
-docs
+^vignettes/web_only$
 ^.*\.Rproj$
 ^\.Rproj\.user$
 ^Notes\.txt$

.github/workflows/R-CMD-check.yaml

@@ -33,7 +33,7 @@ jobs:
 
       - uses: r-lib/actions/setup-pandoc@v2
 
-      - uses: r-lib/actions/setup-tinytex@master
+      - uses: r-lib/actions/setup-tinytex@v2
 
       - uses: r-lib/actions/setup-r@v2
         with:

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: AlphaSimR
 Type: Package
 Title: Breeding Program Simulations
-Version: 1.3.2
-Date: 2022-11-2
+Version: 1.3.4
+Date: 2022-12-8
 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", 
@@ -35,7 +35,7 @@ Depends: R (>= 4.0.0), methods, R6
 Imports: Rcpp (>= 0.12.7), Rdpack
 RdMacros: Rdpack
 LinkingTo: Rcpp, RcppArmadillo (>= 0.7.500.0.0), BH
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.2
 Suggests: knitr, rmarkdown, testthat
 VignetteBuilder: knitr
 NeedsCompilation: true

NEWS.md

@@ -1,3 +1,15 @@
+# AlphaSimR 1.3.4
+
+*changed C++ using `sprintf` to use `snprintf`
+
+# AlphaSimR 1.3.3
+
+*fixed bug in calculation of genic variance
+
+*fixed `importHaplo` not passing ploidy to `newMapPop`
+
+*fixed bug with correlated error variances
+
 # AlphaSimR 1.3.2
 
 *fixed column name bug with multiple traits in `setEBV`

R/importData.R

@@ -70,10 +70,17 @@ importGenMap = function(genMap){
 #' map position (Morgans). Marker name and chromosome are 
 #' coerced using as.character. See \link{importGenMap}
 #' @param ped an optional pedigree for the supplied 
-#' genotypes. The first three columns must be: id, 
-#' mother, and father. All values are coerced using 
-#' as.character.
+#' genotypes. See details. 
 #'
+#' @details 
+#' The optional pedigree can be a data.frame, matrix or a vector. 
+#' If the object is a data.frame or matrix, the first three 
+#' columns must include information in the following order: id, 
+#' mother, and father. All values are coerced using 
+#' as.character. If the object is a vector, it is assumed to only 
+#' include the id. In this case, the mother and father will be set 
+#' to "0" for all individuals.
+#' 
 #' @return a \code{\link{MapPop-class}} if ped is NULL,
 #' otherwise a \code{\link{NamedMapPop-class}}
 #' 
@@ -98,11 +105,18 @@ importGenMap = function(genMap){
 importInbredGeno = function(geno, genMap, ped=NULL){
   # Extract pedigree, if supplied
   if(!is.null(ped)){
-    id = as.character(ped[,1])
-    stopifnot(length(id)==nrow(geno),
-              !any(duplicated(id)))
-    mother = as.character(ped[,2])
-    father = as.character(ped[,3])
+    if(is.vector(ped)){
+      id = as.character(ped[,1])
+      stopifnot(length(id)==nrow(geno),
+                !any(duplicated(id)))
+      mother = father = rep("0", length(id))
+    }else{
+      id = as.character(ped[,1])
+      stopifnot(length(id)==nrow(geno),
+                !any(duplicated(id)))
+      mother = as.character(ped[,2])
+      father = as.character(ped[,3])
+    }
   }
 
   genMap = importGenMap(genMap)
@@ -185,9 +199,16 @@ importInbredGeno = function(geno, genMap, ped=NULL){
 #' coerced using as.character. See \code{\link{importGenMap}}
 #' @param ploidy ploidy level of the organism
 #' @param ped an optional pedigree for the supplied 
-#' genotypes. The first three columns must be: id, 
+#' genotypes. See details. 
+#' 
+#' @details 
+#' The optional pedigree can be a data.frame, matrix or a vector. 
+#' If the object is a data.frame or matrix, the first three 
+#' columns must include information in the following order: id, 
 #' mother, and father. All values are coerced using 
-#' as.character.
+#' as.character. If the object is a vector, it is assumed to only 
+#' include the id. In this case, the mother and father will be set 
+#' to "0" for all individuals.
 #'
 #' @return a \code{\link{MapPop-class}} if ped is NULL,
 #' otherwise a \code{\link{NamedMapPop-class}}
@@ -216,11 +237,18 @@ importInbredGeno = function(geno, genMap, ped=NULL){
 importHaplo = function(haplo, genMap, ploidy=2L, ped=NULL){
   # Extract pedigree, if supplied
   if(!is.null(ped)){
-    id = as.character(ped[,1])
-    stopifnot(length(id)==(nrow(haplo)/ploidy),
-              !any(duplicated(id)))
-    mother = as.character(ped[,2])
-    father = as.character(ped[,3])
+    if(is.vector(ped)){
+      id = as.character(ped[,1])
+      stopifnot(length(id)==(nrow(haplo)/ploidy),
+                !any(duplicated(id)))
+      mother = father = rep("0", length(id))
+    }else{
+      id = as.character(ped[,1])
+      stopifnot(length(id)==(nrow(haplo)/ploidy),
+                !any(duplicated(id)))
+      mother = as.character(ped[,2])
+      father = as.character(ped[,3])
+    }
   }
 
   genMap = importGenMap(genMap)
@@ -252,7 +280,8 @@ importHaplo = function(haplo, genMap, ploidy=2L, ped=NULL){
   }
 
   founderPop = newMapPop(genMap=genMap,
-                         haplotypes=haplotypes)
+                         haplotypes=haplotypes,
+                         ploidy=ploidy)
 
   if(!is.null(ped)){
     founderPop = new("NamedMapPop",

R/misc.R

@@ -558,7 +558,7 @@ usefulness = function(pop,trait=1,use="gv",p=0.1,
 #' Creates an m by m linear transformation matrix that 
 #' can be applied to n by m uncorrelated deviates 
 #' sampled from a standard normal distribution to produce
-#' create correlated deviates with an arbitrary correlation 
+#' correlated deviates with an arbitrary correlation 
 #' of R. If R is not positive semi-definite, the function 
 #' returns smoothing and returns a warning (see details).
 #' 
@@ -569,11 +569,12 @@ usefulness = function(pop,trait=1,use="gv",p=0.1,
 #' matrix and used to test if it is positive semi-definite. 
 #' If the matrix is not positive semi-definite, it is not a 
 #' valid correlation matrix. In this case, smoothing is 
-#' applied to the matrix (as described in the 'psych' library) 
-#' to obtain a valid correlation matrix. The resulting 
-#' deviates will thus not exactly match the desired correlation, 
-#' but will hopefully be close if the the input matrix wasn't 
-#' too far removed from a valid correlation matrix.
+#' applied to the matrix (as described in the 'cor.smooth' of 
+#' the 'psych' library) to obtain a valid correlation matrix. 
+#' The resulting deviates will thus not exactly match the 
+#' desired correlation, but will hopefully be close if the 
+#' input matrix wasn't too far removed from a valid 
+#' correlation matrix.
 #' 
 #' @examples 
 #' # Create an 2x2 correlation matrix

R/phenotypes.R

@@ -209,7 +209,11 @@ setPheno = function(pop, h2=NULL, H2=NULL, varE=NULL, corE=NULL,
       stopifnot(length(varE)==nTraits)
     }
   }else{
-    varE = simParam$varE[traits]
+    if(is.matrix(simParam$varE)){
+      varE = simParam$varE[traits, traits]
+    }else{
+      varE = simParam$varE[traits]
+    }
   }
 
   # Set error correlations

src/alphaSuite.cpp

@@ -32,7 +32,7 @@ void writeASGenotypes(const arma::Cube<unsigned char> & g,
 
     arma::Col<unsigned char> selectedg = selected0 + selected1;
 
-    std::sprintf(name,"%s",names[i].c_str());
+    std::snprintf(name,0,"%s",names[i].c_str());
     ASout << name;
 
     if (snpchips(i) == 0) {
@@ -77,7 +77,7 @@ void writeASHaplotypes(const arma::Cube<unsigned char> & g,
     for (arma::uword j = 0; j < 2; j ++){
       arma::Col<unsigned char> all = g.slice(i).col(j);
 
-      std::sprintf(name,"%s",names[i].c_str());
+      std::snprintf(name,0,"%s",names[i].c_str());
       ASout << name;
 
       if (snpchips(i) == 0) {

src/calcGenParam.cpp

@@ -338,6 +338,7 @@ Rcpp::List calcGenParam(const Rcpp::S4& trait,
     arma::vec bvE(ploidy+1), ddE(ploidy+1); //Expected for random mating
     double gvMu, gvEMu, genoMu, p, q, dK, alpha, alphaE;
 
+    // Compute genotype frequencies
     for(arma::uword j=0; j<nInd; ++j){
       freq(genoMat(j,i)) += 1;
     }
@@ -346,33 +347,38 @@ Rcpp::List calcGenParam(const Rcpp::S4& trait,
     p = genoMu/dP;
     q = 1-p;
 
-    // Set effects, means and expected frequencies
+    // Expected genotype frequencies
+    freqE.zeros();
+    for(arma::uword k=0; k<(ploidy+1); ++k){
+      dK = double(k);
+      freqE(k) = choose(dP,dK)*std::pow(p,dK)*std::pow(q,dP-dK);
+    }
+
+    // Set genetic values
     aEff = xa*a(i);
     if(hasD){
       dEff = xd*d(i);
       gv = aEff+dEff;
     }else{
       gv = aEff;
     }
+
+    // Mean genetic values
     gvMu = accu(freq%gv);
+    gvEMu =  accu(freqE%gv);
     mu(tid) += gvMu;
-    freqE.zeros();
-    for(arma::uword k=0; k<(ploidy+1); ++k){
-      dK = double(k);
-      freqE(k) = choose(dP,dK)*std::pow(p,dK)*std::pow(q,dP-dK);
-    }
+    eMu(tid) += gvEMu;
+
+    // Average effect
     alpha = accu(freq%(gv-gvMu)%(x-genoMu))/
       accu(freq%(x-genoMu)%(x-genoMu));
-    alphaE = accu(freqE%(gv-gvMu)%(x-genoMu))/
+    alphaE = accu(freqE%(gv-gvEMu)%(x-genoMu))/
       accu(freqE%(x-genoMu)%(x-genoMu)); 
 
-    //Check for divide by zero
+    // Check for divide by zero
     if(!std::isfinite(alpha)) alpha=0;
     if(!std::isfinite(alphaE)) alphaE=0;
 
-    gvEMu =  accu(freqE%gv);
-    eMu(tid) += gvEMu;
-
     // Set additive genic variances
     bv = (x-genoMu)*alpha; //Breeding values
     bvE = (x-genoMu)*alphaE; //Random mating breeding value