我应该如何计算'二进制'矩阵中唯一行的数量？
假设我有一个矩阵,其条目仅为0和1,例如
set.seed(123)
m <- matrix( sample(0:1, 10, TRUE), nrow=5 )
带样本输出:
[,1] [,2]
[1,] 0 0
[2,] 1 1
[3,] 0 1
[4,] 1 1
[5,] 1 0
矩阵最多有20列,并且会有很多行.
我想要一个函数,让我们调用它rowCounts,返回:
- 特定行出现在矩阵中的次数,和
- 该行第一次出现的索引.
我怎么能解决这个问题？
Rom*_*ois 11
基于Kevin的回答,这是一个使用稍微不同的方法的C++ 11版本:
List rowCounts_2(IntegerMatrix x) {
int n = x.nrow() ;
int nc = x.ncol() ;
std::vector<int> hashes(n) ;
for( int k=0, pow=1; k<nc; k++, pow*=2){
IntegerMatrix::Column column = x.column(k) ;
std::transform( column.begin(), column.end(), hashes.begin(), hashes.begin(), [=]( int v, int h ){
return h + pow*v ;
}) ;
}
using Pair = std::pair<int,int> ;
std::unordered_map<int, Pair> map_counts ;
for( int i=0; i<n; i++){
Pair& p = map_counts[ hashes[i] ] ;
if( p.first == 0){
p.first = i+1 ; // using directly 1-based index
}
p.second++ ;
}
int nres = map_counts.size() ;
IntegerVector idx(nres), counts(nres) ;
auto it=map_counts.begin() ;
for( int i=0; i<nres; i++, ++it){
idx[i] = it->second.first ;
counts[i] = it->second.second ;
}
return List::create( _["counts"] = counts, _["idx"] = idx );
}
这个想法是为了速度交换记忆.第一个变化是我正在分配并填充一个std::vector<int>来托管哈希.这样做允许我逐列遍历输入矩阵,这更有效.
完成后,我正在训练对的哈希映射(索引,计数)std::unordered_map<int, std::pair<int,int>>.映射的关键是散列,值是一对(索引,计数).
然后我只需要遍历哈希映射并收集结果.结果不会按升序出现idx(如果我们真的想要的话很容易做到).
我得到这些结果n=1e5和n=1e7.
> m <- matrix(sample(0:1, 1e+05, TRUE), ncol = 10)
> microbenchmark(rowCounts(m), rowCountsR(m), rowCounts_2(m))
Unit: microseconds
expr min lq median uq max neval
rowCounts(m) 1194.536 1201.273 1213.1450 1231.7295 1286.458 100
rowCountsR(m) 575.004 933.637 962.8720 981.6015 23678.451 100
rowCounts_2(m) 421.744 429.118 442.5095 455.2510 530.261 100
> m <- matrix(sample(0:1, 1e+07, TRUE), ncol = 10)
> microbenchmark(rowCounts(m), rowCountsR(m), rowCounts_2(m))
Unit: milliseconds
expr min lq median uq max neval
rowCounts(m) 97.22727 98.02716 98.56641 100.42262 102.07661 100
rowCountsR(m) 57.44635 59.46188 69.34481 73.89541 100.43032 100
rowCounts_2(m) 22.95741 23.38186 23.78068 24.16814 27.44125 100
利用线程有助于进一步发挥作用.下面是我的机器上4个线程之间的时间分配方式.请参阅此要点中的代码.
以下是最新版本的基准测试:
> microbenchmark(rowCountsR(m), rowCounts_1(m), rowCounts_2(m), rowCounts_3(m,4))
Unit: milliseconds
expr min lq median uq max neval
rowCountsR(m) 93.67895 127.58762 127.81847 128.03472 151.54455 100
rowCounts_1(m) 120.47675 120.89169 121.31227 122.86422 137.86543 100
rowCounts_2(m) 28.88102 29.68101 29.83790 29.97112 38.14453 100
rowCounts_3(m, 4) 12.50059 12.68981 12.87712 13.10425 17.21966 100
- 我在这里更新了文本并将代码放在这个要点中:https://gist.github.com/romainfrancois/10016972. (2认同)
我们可以利用矩阵的结构以一种很好的方式计算唯一行的数量.因为值都是0和1,我们可以定义一个'hash'函数,它将每一行映射到一个唯一的整数值,然后计算这些哈希值.
我们将实现的哈希函数与以下R代码相同:
hash <- function(x) sum(x * 2^(0:(length(x)-1)))
其中x是0s和1s 的整数向量,表示矩阵的一行.
在我的解决方案中,因为我使用C++并且没有关联容器来维护插入顺序(在标准库中),所以我使用a std::map<int, int>来计算每行的哈希值,并使用a std::vector<int>来跟踪插入哈希值的顺序.
由于列数<= 20的限制,我们可以计算散列值并存储在一个中int,但是对于较大的矩阵是安全的,应该将散列存储在a中double(因为溢出会发生n > 31)
考虑到这一点,我们可以编写一个解决方案:
#include <Rcpp.h>
using namespace Rcpp;
inline int hash(IntegerMatrix::Row x) {
int n = x.size();
int hash = 0;
for (int j=0; j < n; ++j) {
hash += x[j] << j;
}
return hash;
}
// [[Rcpp::export]]
List rowCounts(IntegerMatrix x) {
int nrow = x.nrow();
typedef std::map<int, int> map_t;
map_t counts;
// keep track of insertion order with a separate vector
std::vector<int> ordered_hashes;
std::vector<int> insertion_order;
ordered_hashes.reserve(nrow);
insertion_order.reserve(nrow);
for (int i=0; i < nrow; ++i) {
IntegerMatrix::Row row = x(i, _);
int hashed_row = hash(row);
if (!counts[hashed_row]) {
ordered_hashes.push_back(hashed_row);
insertion_order.push_back(i);
}
++counts[hashed_row];
}
// fill the 'counts' portion of the output
int n = counts.size();
IntegerVector output = no_init(n);
for (int i=0; i < n; ++i) {
output[i] = counts[ ordered_hashes[i] ];
}
// fill the 'idx' portion of the output
IntegerVector idx = no_init(n);
for (int i=0; i < n; ++i) {
idx[i] = insertion_order[i] + 1; // 0 to 1-based indexing
}
return List::create(
_["counts"] = output,
_["idx"] = idx
);
}
/*** R
set.seed(123)
m <- matrix( sample(0:1, 10, TRUE), nrow=5 )
rowCounts(m)
m <- matrix( sample(0:1, 1E5, TRUE), ncol=5 )
str(rowCounts(m))
## Compare it to a close-ish R solution
microbenchmark( times=5,
rowCounts(m),
table(do.call(paste, as.data.frame(m)))
)
*/
打电话sourceCpp给我:
> Rcpp::sourceCpp('rowCounts.cpp')
> set.seed(123)
> m <- matrix( sample(0:1, 10, TRUE), nrow=5 )
> m
[,1] [,2]
[1,] 0 0
[2,] 1 1
[3,] 0 1
[4,] 1 1
[5,] 1 0
> rowCounts(m)
$counts
[1] 1 2 1 1
$idx
[1] 1 2 3 5
> m <- matrix( sample(0:1, 1E5, TRUE), ncol=5 )
> str(rowCounts(m))
List of 2
$ counts: int [1:32] 602 640 635 624 638 621 622 615 633 592 ...
$ idx : int [1:32] 1 2 3 4 5 6 7 8 9 10 ...
> microbenchmark( times=5,
+ rowCounts(m),
+ table(do.call(paste, as.data.frame(m)))
+ )
Unit: milliseconds
expr min lq median uq max neval
rowCounts(m) 1.14732 1.150512 1.172886 1.183854 1.184235 5
table(do.call(paste, as.data.frame(m))) 22.95222 23.146423 23.607649 24.455728 24.953177 5
我很好奇纯R解决方案将如何执行:
set.seed(123)
m <- matrix( sample(0:1, 1E5, TRUE), ncol=5 )
rowCountsR <- function(x) {
## calculate hash
h <- m %*% matrix(2^(0:(ncol(x)-1)), ncol=1)
i <- which(!duplicated(h))
counts <- tabulate(h+1)
counts[order(h[i])] <- counts
list(counts=counts, idx=i)
}
library("rbenchmark")
benchmark(rowCounts(m), rowCountsR(m))
# test replications elapsed relative user.self sys.self user.child sys.child
# 1 rowCounts(m) 100 0.189 1.000 0.188 0 0 0
# 2 rowCountsR(m) 100 0.258 1.365 0.256 0 0 0
编辑:更多列,感谢@Arun指出这一点.
set.seed(123)
m <- matrix( sample(0:1, 1e7, TRUE), ncol=10)
benchmark(rowCounts(m), rowCountsR(m), replications=100)
# test replications elapsed relative user.self sys.self user.child sys.child
#1 rowCounts(m) 100 20.659 1.077 20.533 0.024 0 0
#2 rowCountsR(m) 100 19.183 1.000 15.641 3.408 0 0