Implement Issue rqtl#184, calc_het can use multi-core

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: qtl2
-Version: 0.23-8
-Date: 2020-11-18
+Version: 0.23-9
+Date: 2020-12-02
 Title: Quantitative Trait Locus Mapping in Experimental Crosses
 Description: Provides a set of tools to perform quantitative
     trait locus (QTL) analysis in experimental crosses. It is a

NEWS.md

@@ -1,4 +1,4 @@
-## qtl2 0.23-8 (2020-11-18)
+## qtl2 0.23-9 (2020-12-02)
 
 ### Major changes
 
@@ -16,6 +16,8 @@
 - Updated mouse gene database with 2020-09-07 data from
   [MGI](http://www.informatics.jax.org/downloads/mgigff3/archive/monthly/).
 
+- Implemented Issue #184, to make `calc_het()` multi-core.
+
 ### Bug fixes
 
 - Fixed [Issue #181](https://github.com/rqtl/qtl2/issues/181), where

R/calc_het.R

@@ -8,6 +8,11 @@
 #'
 #' @param omit_x If TRUE, omit the X chromosome.
 #'
+#' @param cores Number of CPU cores to use, for parallel calculations.
+#' (If `0`, use [parallel::detectCores()].)
+#' Alternatively, this can be links to a set of cluster sockets, as
+#' produced by [parallel::makeCluster()].
+#'
 #' @return
 #' The result is a vector of estimated heterozygosities
 #'
@@ -23,13 +28,16 @@
 #'
 #' @export
 calc_het <-
-    function(probs, by=c("individual", "marker"), omit_x=TRUE)
+    function(probs, by=c("individual", "marker"), omit_x=TRUE, cores=1)
 {
     by <- match.arg(by)
     if(is.null(probs)) stop("probs is NULL")
     if(is.cross2(probs))
         stop('Input probs is a "cross2" object but should be genotype probabilities, as from calc_genoprob')
 
+    # set up cluster
+    cores <- setup_cluster(cores, quiet=TRUE)
+
     is_x_chr <- attr(probs, "is_x_chr")
     if(any(is_x_chr) && !all(is_x_chr) && omit_x) {
         probs <- probs[,!is_x_chr]
@@ -54,7 +62,7 @@ calc_het <-
         total_mar <- sum(dim(probs)[3,])
 
         # summarize each chromosome
-        result <- lapply(seq_len(n_chr), function(chr) apply(probs[[chr]][,het_col[[chr]],,drop=FALSE], 1, sum))
+        result <- cluster_lapply(cores, seq_len(n_chr), function(chr) apply(probs[[chr]][,het_col[[chr]],,drop=FALSE], 1, sum))
 
         if(length(result)>1) {
             for(i in seq_along(result)[-1])
@@ -64,6 +72,6 @@ calc_het <-
     }
 
     # else: by marker
-    unlist(lapply(seq_len(n_chr), function(chr) apply(probs[[chr]][,het_col[[chr]],,drop=FALSE], 3, sum)/
+    unlist(cluster_lapply(cores, seq_len(n_chr), function(chr) apply(probs[[chr]][,het_col[[chr]],,drop=FALSE], 3, sum)/
                                                           nrow(probs[[chr]])))
 }

tests/testthat/test-calc_het.R

@@ -25,3 +25,18 @@ test_that("calc_het works for an intercross", {
     expect_equal(calc_het(pr[,c(19,"X")], "marker", omit_x=FALSE), expected)
 
 })
+
+test_that("calc_het works with multi-core", {
+
+    skip_if(isnt_karl(), "This test only run locally")
+
+    iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
+    pr <- calc_genoprob(iron, err=0.002)
+
+    # test multicore
+    expect_equal(calc_het(pr), calc_het(pr, cores=4))
+    expect_equal(calc_het(pr, "marker"), calc_het(pr, "marker", cores=4))
+    expect_equal(calc_het(pr, omit_x=FALSE), calc_het(pr, omit_x=FALSE, cores=4))
+    expect_equal(calc_het(pr, "marker", omit_x=FALSE), calc_het(pr, "marker", omit_x=FALSE, cores=4))
+
+})