匹配邏輯語句的Rcpp矩陣的子集

Question

在R中，如果我們有一個數據矩陣，比如一個100×10的矩陣X和一個具有可能值（0,1,2,3）的100個元素的向量t，我們可以用一個簡單的語法很容易地找到X的子矩陣Y：

y = X[t == 1, ] 

但是，問題是，我怎麼能做到這一點與RCPP的NumericMatrix？
（或者更一般地說，我該怎麼做，在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

我所知道的最接近的是find()函數加上submat()函數的組合，Armadillo可以通過RcppArmadillo訪問。

編輯：這當然是我們可以通過補丁添加的東西。如果有人有足夠的動機去嘗試這個，請到rcpp-devel郵件列表。

Answer 2

我很想看到這個糖。不幸的是，我沒有資格實施它。這仍然是我玩過的許多不同的解決方案。

首先，我不得不做出一些修改功一撩代碼得到這個工作（colvec而不是vec爲tIdx和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)); 

其次，這裏有三個函數所有基於邏輯語句的子集矩陣的基準函數採用arma或rcpp的參數和返回值兩個基於Gong-Yi Liao的解決方案，一個是簡單的基於循環的解決方案

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) 

*/