c++ - 用C++和Rcpp重写慢R函数

我有这行R代码:

croppedDNA <- completeDNA[,apply(completeDNA,2,function(x) any(c(FALSE,x[-length(x)]!=x[-1])))]

它的作用是识别不是普遍的(信息性的)DNA序列矩阵(1行= 1个seq)中的位点(cols)，并将其从矩阵中子集化，以形成新的“裁剪矩阵”，即摆脱所有值相同的列。对于大数据集，这大约需要6秒钟。我不知道我是否可以在C++中更快地做到这一点(仍然是C++的初学者)，但是对我来说尝试尝试将是一件好事。我的想法是使用Rcpp，遍历CharacterMatrix的列，将列(站点)拉出作为CharacterVector，以检查它们是否相同。如果它们相同，则记录该列号/索引，并继续所有列。然后最后创建一个仅包含这些列的新CharacterMatrix。重要的是，应将行名和列名保持在矩阵的“R版本”中，例如，如果有列，则列名也应如此。

我已经写了大约两分钟，到目前为止，我所拥有的(还没有完成):

#include <Rcpp.h>
#include <vector>
using namespace Rcpp;
// [[Rcpp::export]]
CharacterMatrix reduce_sequences(CharacterMatrix completeDNA)
{
  std::vector<bool> informativeSites;
  for(int i = 0; i < completeDNA.ncol(); i++)
  {
    CharacterVector bpsite = completeDNA(,i);
    if(all(bpsite == bpsite[1])
    {
      informativeSites.push_back(i);
    }
  }
CharacterMatrix cutDNA = completeDNA(,informativeSites);
return cutDNA;
}

我在这方面走对了吗？有没有更简单的方法。我的理解是我需要std::vector，因为很容易种植它们(因为我事先不知道我要保留多少个cols)。使用索引时，我是否需要在最后将+1到INFORMATIONSITES vector (因为R索引从1开始，C++从0开始)？

谢谢，
本·W。

最佳答案

样本数据:

set.seed(123)
z <- matrix(sample(c("a", "t", "c", "g", "N", "-"), 3*398508, TRUE), 3, 398508)

OP的解决方案:

system.time(y1 <- z[,apply(z,2,function(x) any(c(FALSE,x[-length(x)]!=x[-1])))])
#    user  system elapsed
#   4.929   0.043   4.976

使用基本R的更快版本:

system.time(y2 <- (z[, colSums(z[-1,] != z[-nrow(z), ]) > 0]))
#    user  system elapsed
#   0.087   0.011   0.098

结果是相同的:

identical(y1, y2)
# [1] TRUE

c++很有可能会击败它，但这真的有必要吗？

关于c++ - 用C++和Rcpp重写慢R函数，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/16555752/

RCPP

c++ - 用C++和Rcpp重写慢R函数