重复加入大型数据集

Question

需要组合多个具有2列的矩阵,如下所示

matrix1
1,3
1,5
3,6

matrix2
1,4
1,5
3,6
3,7

output
1,3,1
1,4,1
1,5,2
3,6,2
3,7,1

输出中的第三列是在所有矩阵中看到一对的次数.我写了一些代码来做到这一点

require(data.table)

set.seed(1000)
data.lst <- lapply(1:200, function(n) { x <- matrix(sample(1:1000,2000,replace=T), ncol=2); x[!duplicated(x),] })

#method 1
pair1.dt <- data.table(i=integer(0), j=integer(0), cnt=integer(0))
for(mat in data.lst) {
    pair1.dt <- rbind(pair1.dt, data.table(i=mat[,1],j=mat[,2],cnt=1))[, .(cnt=sum(cnt)), .(i,j)]
}

#method 2
pair2.dt <- data.table(i=integer(0), j=integer(0), cnt=integer(0))
for(mat in data.lst) {
    pair2.dt <- merge(pair2.dt, data.table(i=mat[,1],j=mat[,2],cnt=1), by=c("i","j"), all=T)[, 
        cnt:=rowSums(.SD,na.rm=T), .SDcols=c("cnt.x","cnt.y")][, c("cnt.x","cnt.y"):=NULL]
}

cat(sprintf("num.rows  =>  pair1: %d,  pair2: %d", pair1.dt[,.N], pair2.dt[,.N]), "\n")

在实际问题中,每个矩阵具有数百万行,并且可能有30-40%的重叠.我试图找出最快的方法来做到这一点.我尝试使用Matrix :: sparseMatrix.虽然这要快得多,但我遇到了一个错误"尚未支持的长向量".我在这里有几种不同的基于data.table的方法.我正在寻找加快此代码和/或建议替代方法的建议.

Answer 1

首先,制作data.tables:

dt.lst = lapply(data.lst, as.data.table)

堆叠.为了比较,这是涉及堆叠的快速方式:

res0 = rbindlist(dt.lst)[, .(n = .N), by=V1:V2]

OP表示这是不可行的,因为由此产生的中间结果rbindlist会太大.

首先列举.对于一小部分值,我建议先预先列出它们:

res1 = CJ(V1 = 1:1000, V2 = 1:1000)[, n := 0L]
for (k in seq_along(dt.lst)) res1[ dt.lst[[k]], n := n + .N, by=.EACHI ] 

fsetequal(res0, res1[n>0]) # TRUE

OP表示有1e12个可能的值,所以这似乎不是一个好主意.相反,我们可以使用

res2 = dt.lst[[1L]][0L]
for (k in seq_along(dt.lst)) res2 = funion(res2, dt.lst[[k]])
res2[, n := 0L]
setkey(res2, V1, V2)
for (k in seq_along(dt.lst)) res2[ dt.lst[[k]], n := n + .N, by=.EACHI ]     

fsetequal(res0, res2) # TRUE

这是给出的例子中三个选项中最慢的一个,但鉴于OP的担忧,对我来说似乎是最好的.

在循环中成长.最后...

res3 = dt.lst[[1L]][0L][, n := NA_integer_][]
for (k in seq_along(dt.lst)){
  setkey(res3, V1, V2)
  res3[dt.lst[[k]], n := n + .N, by=.EACHI ]
  res3 = rbind(
    res3, 
    fsetdiff(dt.lst[[k]], res3[, !"n", with=FALSE], all=TRUE)[, .(n = .N), by=V1:V2]
  )
} 

fsetequal(res0, res3) # TRUE

在循环中生长对象在R中是非常不鼓励和低效的,但是这允许您在单个循环而不是两个循环中执行它.

其他选项和说明.我怀疑你最好使用哈希.这些可以在散列包中找到,也可能通过Rcpp包.

fsetequal,fsetdiff并且funion是该软件包的开发版本的最新补充.在data.table项目的官方网站上查找详细信息.

顺便说一句,如果每个矩阵内的条目是不同的,你可以替换.N与1L无处不在上面和下降by=.EACHI和all=TRUE.

Answer 2

使用 Rcpp。此方法将利用 std::unordered_map 的哈希属性。

#include "Rcpp.h"
#include <stdint.h>
#include <unordered_map>

using namespace std;
using namespace Rcpp;

//[[Rcpp::export]]
Rcpp::XPtr<int> CreateMap(){
  std::unordered_map<int64_t, int>* myMap = new std::unordered_map<int64_t, int>();
  Rcpp::XPtr<int> p((int*)myMap,false);
  return p;
}

//[[Rcpp::export]]
void FreeMap(Rcpp::XPtr<int> map_ptr){
  std::unordered_map<int64_t, int>* myMap =  (std::unordered_map<int64_t, int>*)(int*)map_ptr;
  delete myMap;
}

//[[Rcpp::export]]
void AccumulateValues(Rcpp::XPtr<int> map_ptr, SEXP mat){

  NumericMatrix m(mat);

  std::unordered_map<int64_t, int>* myMap =  (std::unordered_map<int64_t, int>*)(int*)map_ptr;
  for(int i = 0; i<m.nrow(); i++){
    int c1 = m(i, 0);
    int c2 = m(i, 1);
    int64_t key = ((int64_t)c1 << 32) + c2;
    (*myMap)[key] ++;

  }
}
//[[Rcpp::export]]
SEXP AsMatrix(Rcpp::XPtr<int> map_ptr){
  std::unordered_map<int64_t, int>* myMap =  (std::unordered_map<int64_t, int>*)(int*)map_ptr;
  NumericMatrix m(myMap->size(),3);
  int index = 0;
  for ( auto it = myMap->begin(); it != myMap->end(); ++it ){
    int64_t key = it->first;
    m(index, 0) = (int)(key >> 32);
    m(index, 1) = (int)key;
    m(index, 2) = it->second;
    index++;
  }
  return m;
}

那么R代码是：

myMap<-CreateMap()
AccumulateValues(myMap, matrix1)
AccumulateValues(myMap, matrix2)
result<-AsMatrix(myMap)
FreeMap(myMap)

还要求

PKG_CXXFLAGS = "-std=c++0x"

在 makevars 包中