0

which(x,arr.ind=T)Rcpp や 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::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));
}

このコードは非常に効率が悪いように見えますが、R の関数microbenchmarkよりも高速であることがわかります。which

質問: この関数をさらに変更して、R のwhich関数を実際に正確に再現する、つまり渡すMATRIX == somethingことはできますか? 今のところ、そのための2番目の引数が必要です。利便性のためにこれを持っているだけです。

更新: バグを修正しました - 床の代わりに天井が必要でした

確認方法:

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
4

2 に答える 2

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 レベルで呼び出しを解析するのに時間をかける価値があるかもしれませんが、それはあなたに任せます。

于 2013-08-26T05:56:56.570 に答える