RCPP中的（x，arr.ind = T）

Question

我無法在Rcpp或RcppArmadillo中找到很酷的which(x,arr.ind=T)功能。所以我決定快速編寫代碼。

// [[Rcpp::export]] 
arma::umat whicha(arma::mat matrix, int what){ 
    arma::uvec outp1; 
    int n = matrix.n_rows; 
    outp1 = find(matrix==what); 
    int nf = outp1.n_elem; 
    arma::mat out(nf,2); 
    arma::vec foo; 
    arma::uvec foo2; 
    foo = arma::conv_to<arma::colvec>::from(outp1) +1; 
    foo2 = arma::conv_to<arma::uvec>::from(foo); 
    for(int i=0; i<nf; i++){ 
    out(i,0) = (foo2(i) %n); 
    out(i,1) = ceil(foo(i)/n); 
    if(out(i,0)==0) { 
     out(i,0)=n; 
    } 
    } 
    return(arma::conv_to<arma::umat>::from(out)); 
} 

代碼似乎非常低效的，但microbenchmark表明，它可以比的r which功能更快。

問題：我可否進一步更改此功能以實際重現R的which功能，即將MATRIX == something傳遞給它？現在我需要第二個參數。我只是爲了方便有這個。

---

更新：修正了一個錯誤 - 需要的小區，而不是地板

如何檢查：

ma=floor(abs(rnorm(100,0,6))) 
testf=function(k) {all(which(ma==k,arr.ind=T) == whicha(ma,k))} ; sapply(1:10,testf) 

基準：

> microbenchmark(which(ma==k,arr.ind=T) , whicha(ma,k)) 
Unit: microseconds 
         expr min  lq median  uq max neval 
which(ma == k, arr.ind = T) 10.264 11.170 11.774 12.377 51.317 100 
       whicha(ma, k) 3.623 4.227 4.830 5.133 36.224 100 

Answer 1

我想通過生成的包裝做到這一點R功能和做一些醜陋的工作來處理呼叫。舉例來說，使用代碼：

whicha.cpp 
---------- 

#include <RcppArmadillo.h> 
// [[Rcpp::depends("RcppArmadillo")]] 

// [[Rcpp::export]] 
arma::umat whicha(arma::mat matrix, int what){ 
    arma::uvec outp1; 
    int n = matrix.n_rows; 
    outp1 = find(matrix==what); 
    int nf = outp1.n_elem; 
    arma::mat out(nf,2); 
    arma::vec foo; 
    arma::uvec foo2; 

    foo = arma::conv_to<arma::vec>::from(outp1) +1; 
    out.col(1) = floor( foo/n) +1; 
    foo2 = arma::conv_to<arma::uvec >::from(foo); 
    for(int i=0; i<nf; i++){ 
    out(i,0) = foo2(i) % n; 
    } 

    return(arma::conv_to<arma::umat >::from(out)); 
} 

/*** R 
whichRcpp <- function(x) { 
    call <- match.call()$x 
    xx <- eval.parent(call[[2]]) 
    what <- eval.parent(call[[3]]) 
    return(whicha(xx, what)) 
} 
x <- matrix(1:1E4, nrow=1E2) 
identical(whichRcpp(x == 100L), whicha(x, 100L)) ## TRUE 
microbenchmark::microbenchmark(whichRcpp(x == 100L), whicha(x, 100L)) 
*/ 

不幸的是，microbenchmark表明我解析調用是有點慢：

Unit: microseconds 
       expr min  lq median  uq max neval 
whichRcpp(x == 100L) 43.542 44.143 44.443 45.0440 73.271 100 
     whicha(x, 100L) 30.029 30.630 30.930 31.2305 78.075 100 

這可能是值得你花時間在C層次來解析呼叫，但我會留給你的。

Answer 2

這裏只使用RCPP我的代碼：

src <- ' 
    using namespace std; 

    NumericMatrix X(X_); 
    double what = as<double>(what_); 
    int n_rows = X.nrow(); 

    NumericVector rows(0); 
    NumericVector cols(0); 

    for(int ii = 0; ii < n_rows * n_rows; ii++) 
    { 
     if(X[ii] == what) 
     { 
      rows.push_back(ii % n_rows + 1); 
      cols.push_back(floor(ii/n_rows) + 1); 
     } 
    } 

    return List::create(rows, cols); 
' 

fun <- inline:::cxxfunction(signature(X_ = 'numeric', what_ = 'numeric'), src, 'Rcpp') 

X <- matrix(1:1E4, nrow=1E2) 

rbenchmark:::benchmark(fun(X, 100), which(X == 100L, TRUE), columns = c('test', 'replications', 'elapsed', 'relative'), replications = 1000) 

        test replications elapsed relative 
1   fun(X, 100)   1000 0.077 1.000 
2 which(X == 100, TRUE)   1000 0.100 1.299 

microbenchmark:::microbenchmark(fun(X, 100), which(X == 100L, TRUE), times = 1000L) 

        expr min  lq median  uq  max neval 
      fun(X, 100) 37.372 41.3780 43.6530 48.4825 1650.154 1000 
which(X == 100L, TRUE) 63.366 64.0745 64.3345 64.8240 1911.858 1000 

相比，從以前的海報解決

慢不了多少。有趣的是，返回數據幀而不是列表會顯着降低性能。