【问题标题】:Finding duplicates using Rcpp使用 Rcpp 查找重复项
【发布时间】:2017-06-07 13:37:23
【问题描述】:

我正在尝试找到一种更快的替代方法来查找 R 中的重复项。代码的目的是将矩阵与该矩阵中的行号一起传递给 Rcpp,然后循环遍历整个矩阵以寻找匹配项排。有问题的矩阵是一个具有 1000 行和 250 列的逻辑矩阵。

听起来很简单,但下面的代码没有检测等效向量行。我不确定这是否与 equal() 函数或矩阵或向量的定义方式有关。

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::plugins]]
#include <cstddef>   // std:size_t
#include <iterator>  // std:begin, std::end
#include <vector>    // std::vector
#include <iostream>
#include <string>

// [[Rcpp::export]]
    bool dupCheckRcpp (int nVector, 
                        LogicalMatrix bigMatrix) {
    // initialize
      int i, j, nrow, ncol;
      nrow = bigMatrix.nrow();
      ncol = bigMatrix.ncol();
      LogicalVector vec(ncol);  // holds vector of interest
      LogicalVector vecMatrix(ncol); // temp vector for loop through bigMatrix
      nVector = nVector - 1;

    // copy bigMatrix data into vec based on nVector row
      for ( j = 0; j < ncol; ++j ) {
        vec(j) = bigMatrix(nVector,j);
      }

    // check loop: check vecTeam against each row in allMatrix
      for (i = 0; i < nrow; ++i) {  
        // copy bigMatrix data into vecMatrix
          for ( j = 0; j < ncol; ++j ) {
            vecMatrix(j) = bigMatrix(i,j);
          }
        // check for equality
          if (i != nVector) {  // skip if nVector row
            // compare vecTeam to vecMatrix
              if (std::equal(vec.begin(),vec.end(),vecMatrix.begin())) {
              return true;
            }
          }
      } // close check loop
      return false;
    }

【问题讨论】:

  • 如果你查看在矩阵上调用的duplicated 的调用堆栈,它几乎所有的时间都花在paste 中,而paste 已经在C 中运行。让它更快可能并不容易。
  • which(colSums(t(mat) == mat[20,]) == ncol(mat)) 在我的系统上大约需要 1.5 毫秒。真的太慢了​​吗?
  • 此外,如果您可以在前面计算 all,则类似 tmp = tcrossprod(mat); which((tmp == rowSums(mat)) &amp; lower.tri(tmp), TRUE) 的内容会返回 all 匹配行。并且中间矩阵的大小为1e3 * 1e3,这似乎是可管理的

标签: c++ r matrix duplicates rcpp


【解决方案1】:

我不确定代码中的错误在哪里,但请注意,您真的不需要像这样在 Rcpp 类型之间手动复制元素:

// copy bigMatrix data into vec based on nVector row
for (j = 0; j < ncol; ++j) {
    vec(j) = bigMatrix(nVector, j);
}

几乎总会有合适的类和/或合适的赋值运算符等,让您更简洁、更安全地完成此任务(即不易出现编程错误)。这是一个更简单的实现:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
bool is_duplicate_row(R_xlen_t r, LogicalMatrix x) {
    R_xlen_t i = 0, nr = x.nrow();
    const LogicalMatrix::Row& y = x.row(r);

    for (; i < r; i++) {
        if (is_true(all(y == x.row(i)))) {
            return true;
        }
    }
    for (i = r + 1; i < nr; i++) {
        if (is_true(all(y == x.row(i)))) {
            return true;
        }
    }

    return false;
}

本着我上述建议的精神,

  • const LogicalMatrix::Row&amp; y = x.row(r); 为我们提供了对矩阵第 rth 行的常量引用
  • x.row(i) 指的是xith 行

这两个表达式都避免了通过 for 循环进行元素复制,并且在 IMO 中更具可读性。此外,虽然使用std::equal 或任何其他标准算法肯定没有错,但使用诸如is_true(all(y == x.row(i))) 之类的Rcpp 糖表达式通常可以进一步简化您的代码。


set.seed(123)
m <- matrix(rbinom(1000 * 250, 1, 0.25) > 0, 1000)
m[600,] <- m[2,]

which(sapply(1:nrow(m) - 1, is_duplicate_row, m))
# [1]   2 600

c(which(duplicated(m, fromLast = TRUE)), which(duplicated(m)))
# [1]   2 600

【讨论】:

  • 很棒的答案,解释很清楚。我的代码从每循环 193.2 秒到 4.7 秒。速度提高了 97%。太感谢了。 PS。我在运行您的代码时确实遇到了一个小编译错误:“从 SEXP 到 long long int 的转换无效。我将 R_xlen_t 术语更改为 int,它运行良好。
  • 使用// [[Rcpp::plugins(cpp11)]] 可以使用R_xlen_t 类型,因为这是long int 的重铸,需要C++11。
猜你喜欢
  • 2022-01-02
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2016-09-29
  • 1970-01-01
相关资源
最近更新 更多