論理文に一致するRcpp行列のサブセット

Question

Rでは、データ行列、たとえば100×10行列X、および100要素のベクトルtが可能な値（0,1,2,3）を持つ場合、問題は、私はRcppのNumericMatrixでそれを行うことができる方法、である、

y = X[t == 1, ] 

しかし：、我々は簡単に単純な構文を使用してXの部分行列yを見つけることができますか？
（または、より一般的には、私が行うことができますどのようにC++の任意のコンテナを？に）ディルクのヒントに

おかげで、それは

NumericMatrix X(dataX); 
IntegerVector T(dataT); 
mat Xmat(X.begin(), X.nrow(), X.ncol(), false); 
vec tIdx(T.begin(), T.size(), false); 
mat y = X.rows(find(tIdx == 1)); 

は私がやりたいことができているようですが、それはあまりにも長いようです。

Answer 1

私の知っている最も近いRcppArmadilloを介してアクセスArmadilloでsubmat()機能と組み合わせるfind()機能を組み合わせたものです。

編集：これは、パッチで追加できるコースです。誰かがこれを試すのに十分な動機があれば、rcpp-develメーリングリストに来てください。

Answer 2

私はこれを砂糖と見なしたいと思います。残念ながら、私はそれを実装する資格はありません。ここには、私がやったさまざまなソリューションがまだあります。

まず、私は（これは仕事を得るためにcolvec代わりにvectIdxためとXmat.rows(...代わりのX.rows(...功・李遼コードにいくつかの変更をしなければならなかった：

mat Xmat(X.begin(), X.nrow(), X.ncol(), false); 
colvec tIdx(T.begin(), T.size(), false); 
mat y = Xmat.rows(find(tIdx == 1)); 

第二に、ここでは3つの機能付きですarmaまたはrcppの引数と戻り値をとる関数2つはGong-Yi Liaoの解を基にしており、1つは単純なループベースの解であり、もう1つは単純なループベースの解です。

N（行）= 100、P（T == 1）= 0.3

   expr min  lq median  uq max 
1 submat_arma(X, T) 5.009 5.3955 5.8250 6.2250 28.320 
2 submat_arma2(X, T) 4.859 5.2995 5.6895 6.1685 45.122 
3 submat_rcpp(X, T) 5.831 6.3690 6.7465 7.3825 20.876 
4  X[T == 1, ] 3.411 3.9380 4.1475 4.5345 27.981 

N（行）= 10000、P（T == 1）= 0.3

   expr  min  lq median  uq  max 
1 submat_arma(X, T) 107.070 113.4000 125.5455 141.3700 1468.539 
2 submat_arma2(X, T) 76.179 80.4295 88.2890 100.7525 1153.810 
3 submat_rcpp(X, T) 244.242 247.3120 276.6385 309.2710 1934.126 
4  X[T == 1, ] 229.884 236.1445 263.5240 289.2370 1876.980 

submat.cpp

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

using namespace Rcpp; 
using namespace arma; 

// arma in; arma out 
// [[Rcpp::export]] 
mat submat_arma(arma::mat X, arma::colvec T) { 
    mat y = X.rows(find(T == 1)); 
    return y; 
} 

// rcpp in; arma out 
// [[Rcpp::export]] 
mat submat_arma2(NumericMatrix X, NumericVector T) { 
    mat Xmat(X.begin(), X.nrow(), X.ncol(), false); 
    colvec tIdx(T.begin(), T.size(), false); 
    mat y = Xmat.rows(find(tIdx == 1)); 
    return y; 
} 

// rcpp in; rcpp out 
// [[Rcpp::export]] 
NumericMatrix submat_rcpp(NumericMatrix X, LogicalVector condition) { 
    int n=X.nrow(), k=X.ncol(); 
    NumericMatrix out(sum(condition),k); 
    for (int i = 0, j = 0; i < n; i++) { 
     if(condition[i]) { 
      out(j,_) = X(i,_); 
      j = j+1; 
     } 
    } 
    return(out); 
} 


/*** R 
library("microbenchmark") 

# simulate data 
n=100 
p=0.3 
T=rbinom(n,1,p) 
X=as.matrix(cbind(rnorm(n),rnorm(n))) 

# compare output 
identical(X[T==1,],submat_arma(X,T)) 
identical(X[T==1,],submat_arma2(X,T)) 
identical(X[T==1,],submat_rcpp(X,T)) 

# benchmark 
microbenchmark(X[T==1,],submat_arma(X,T),submat_arma2(X,T),submat_rcpp(X,T),times=500) 

# increase n 
n=10000 
p=0.3 
T=rbinom(n,1,p) 
X=as.matrix(cbind(rnorm(n),rnorm(n))) 
# benchmark 
microbenchmark(X[T==1,],submat_arma(X,T),submat_arma2(X,T),submat_rcpp(X,T),times=500) 

*/