在 C++ 和 Rcpp 中重写慢 R 函数

Rewriting slow R function in C++ & Rcpp

本文关键字：函数 重写 C++ Rcpp      更新时间：2023-10-16

我有这一行R代码:

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

它所做的是识别DNA序列矩阵(1行=一个序列)中的位点(cols)，这些位点(cols)不是通用的(信息丰富的)，并将它们从矩阵中子集，以形成一个新的"裁剪矩阵"，即去掉所有值相同的列。对于一个大数据集，这大约需要6秒。我不知道我是否可以在c++中做得更快(仍然是c++的初学者)，但这将是对我有好处的尝试。我的想法是使用Rcpp，循环通过CharacterMatrix的列，拉出列(网站)作为CharacterVector检查它们是否相同。如果它们是相同的，记录该列号/索引，继续所有列。然后在最后创建一个只包含这些列的新CharacterMatrix。重要的是，我保留了行名和列名，因为它们在矩阵的"R版本"中，即，如果一列消失了，那么colname也应该消失。

我已经写了大约两分钟了，到目前为止我写的是(未完成):

#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)。有了索引，我是否需要在末尾对informativeSites向量加1(因为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

使用base 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++很有可能打败它，但真的有必要吗?

相关文章：