0

我一直在为我的主管研究一个简单的函数集合，这些函数将执行一些简单的初始基因组规模统计，这很容易让我的团队快速了解可能需要更多时间的未来分析——例如 RDP4 或BioC（只是为了解释为什么我没有直接去 BioConductor）。我想加快一些速度以允许更大的 contig 大小，所以我决定使用 doParallel 和 foreach 来编辑一些 for 循环以允许这样做。下面是一个简单的函数，它识别某些序列（存储为矩阵）中相同的碱基并将它们删除。

strip.invar <- function(x) {
  cat("
          Now removing invariant sites from DNA data matrix, this may take some time...
      ")
  prog <- txtProgressBar(min=0, max=ncol(x), style=3)
  removals<-c()
  for(i in 1:ncol(x)){
    setTxtProgressBar(prog, i)
    if(length(unique(x[,i])) == 1) { 
      removals <- append(removals, i)
    }
  }
  newDnaMatrix <- x[,-removals]
  return(newDnaMatrix)
}

在阅读了 doParallel 和 foreach 的介绍之后，我尝试制作一个版本以容纳更多内核 - 在我的 Mac 上，这是 8 个 - 每个内核有两个线程的四核 - 8 个虚拟内核：

strip.invar <- function(x, coresnum=detectCores()){
  cat("
          Now removing invariant sites from DNA data matrix, this may take some time...
          ")
  prog <- txtProgressBar(min=0, max=ncol(x), style=3)
  removals<-c()
  if(coresnum > 1) {
    cl <- makeCluster(coresnum)
    registerDoParallel(cl)
    foreach(i=1:ncol(x)) %dopar% {
      setTxtProgressBar(prog, i)
      if(all(x[,i] == x[[1,i]])){
        removals <- append(removals, i)
      }
    }
  } else {
    for(i in 1:ncol(x)){
      setTxtProgressBar(prog, i)
      if(length(unique(x[,i])) == 1) { 
        removals <- append(removals, i)
      }
    }
  }
  newDnaMatrix <- x[,-removals]
  return(newDnaMatrix)
}

但是，如果我运行它并将核心数设置为 8，我不完全确定它是否有效 - 我看不到进度条做任何事情，但我听说打印到屏幕和涉及图形设备的东西很棘手在 R 中使用并行计算。但它似乎仍然需要一些时间，而且我的笔记本电脑变得“非常”热，所以我不确定我是否正确地做到了这一点，我在看到几个例子后尝试过（我成功运行了小插图中很好的引导示例），但我一定会遇到学习障碍。顺便说一句，我想我也会问人们的意见，对于涉及循环或应用的 R 代码瓶颈的最佳加速是什么 - 并行化它，还是 Rcpp？

谢谢。

4

2 回答 2

0

首先，尝试运行cl <- makeCluster( coresnum-1 ). 主 R 进程已经在使用您的一个核心，并用于分派和接收从属作业的结果，因此您有 7 个空闲核心用于从属作业。我认为您将有效地排队您的 foreach 循环之一，以等待之前的循环之一完成，因此这项工作将需要更长的时间才能完成。

其次，您通常会在非并行环境中运行此功能的控制台上看到的内容仍会打印到控制台，只是每个作业输出都会打印到从属进程控制台，因此您不会看到它。但是，您可以将不同 foreach 循环的输出保存到文本文件中以检查它们。这是如何保存控制台输出的示例。将代码粘贴在您的foreach语句中。

您的笔记本电脑会变得非常热，因为在您运行此作业时，您的所有内核都以 100% 的容量工作。

我发现该foreach包提供了一组出色的功能来提供简单的并行处理。Rcpp可能（会？！）给您带来更好的性能，但是您在编写 C++ 代码方面的表现如何？这个函数的运行时间是多少？它的使用频率如何？我总是首先考虑这些事情。

于 2013-03-02T13:32:20.563 回答
0

我的另一个答案不正确，因为 colmean 等于第一个值不足以作为唯一值数量的测试。所以这是另一个答案：

您可以使用 优化循环apply。

set.seed(42)
dat <- matrix(sample(1:5,1e5,replace=TRUE),nrow=2)
res1 <- strip.invar(dat)


strip.invar2 <- function(dat) {
  ix <- apply(dat,2,function(x) length(unique(x))>1)
  dat[,ix]}

res2 <- strip.invar2(dat)

all.equal(res1,res2)
#TRUE
library(microbenchmark)
microbenchmark(strip.invar(dat),strip.invar2(dat),times=10)
#Unit: milliseconds
#             expr       min        lq    median       uq      max neval
#strip.invar(dat)  2514.7995 2529.2827 2547.6751 2678.464 2792.405    10
#strip.invar2(dat)  933.3701  945.5689  974.7564 1008.589 1018.400    10

这大大提高了性能，尽管如果矢量化是可能的，那么你可以达到的效果并不多。

并行化在这里不会提供更好的性能，因为每次迭代本身不需要太多性能，因此并行化开销实际上会增加所需的时间。但是，您可以并行拆分数据和处理块。

于 2013-03-02T15:10:50.023 回答