相当于 Rcpp 中的 'which' 函数

Question

我是 C++ 和 Rcpp 的新手。假设，我有一个向量

t1<-c(1,2,NA,NA,3,4,1,NA,5)

我想获得 t1 元素的索引 NA。我会写：

NumericVector retIdxNA(NumericVector x) {

    // Step 1: get the positions of NA in the vector
    LogicalVector y=is_na(x);

    // Step 2: count the number of NA
    int Cnt=0;
    for (int i=0;i<x.size();i++) {
       if (y[i]) {
         Cnt++;
       }
    }

    // Step 3: create an output matrix whose size is same as that of NA
    // and return the answer
    NumericVector retIdx(Cnt);
    int Cnt1=0;
    for (int i=0;i<x.size();i++) {
       if (y[i]) {
          retIdx[Cnt1]=i+1;
          Cnt1++;
       }
    }
    return retIdx;
}

然后我得到

retIdxNA(t1)
[1] 3 4 8

我想知道：

(i) Rcpp 中是否有任何与which 等效的东西？

(ii) 有什么方法可以使上述功能更短/更清晰？特别是，有没有什么简单的方法可以把上面的第1、2、3步结合起来？

Answer 1

RcppArmadillo 的最新版本具有识别有限和非有限值索引的功能。

所以这段代码

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

// [[Rcpp::export]]
arma::uvec whichNA(arma::vec x) {
  return arma::find_nonfinite(x);
}

/*** R
t1 <- c(1,2,NA,NA,3,4,1,NA,5)
whichNA(t1)
*/

产生您想要的答案（在 C/C++ 中将它们逐一模块化，因为它们是从零开始的）：

R> sourceCpp("/tmp/uday.cpp")

R> t1 <- c(1,2,NA,NA,3,4,1,NA,5)

R> whichNA(t1)
     [,1]
[1,]    2
[2,]    3
[3,]    7
R> 

如果您首先创建要子集化的序列，Rcpp 也可以这样做：

// [[Rcpp::export]]
Rcpp::IntegerVector which2(Rcpp::NumericVector x) {
  Rcpp::IntegerVector v = Rcpp::seq(0, x.size()-1);
  return v[Rcpp::is_na(x)];
}

添加到上面的代码会产生：

R> which2(t1)
[1] 2 3 7
R> 

逻辑子集在 Rcpp 中也有些新。

Answer 2

试试这个：

#include <Rcpp.h> 
using namespace Rcpp; 

// [[Rcpp::export]]
IntegerVector which4( NumericVector x) {

    int nx = x.size();
    std::vector<int> y;
    y.reserve(nx);

    for(int i = 0; i < nx; i++) {
        if (R_IsNA(x[i])) y.push_back(i+1);
    }

    return wrap(y);
}

我们可以在 R 中这样运行：

> which4(t1)
[1] 3 4 8

性能

请注意，我们更改了上述解决方案，为输出向量保留空间。这将替换 which3，即：

// [[Rcpp::export]]
IntegerVector which3( NumericVector x) {
    int nx = x.size();
    IntegerVector y;
    for(int i = 0; i < nx; i++) {
        // if (internal::Rcpp_IsNA(x[i])) y.push_back(i+1);
        if (R_IsNA(x[i])) y.push_back(i+1);
    }
    return y;
}

那么在 9 个元素长的向量上的性能如下，which4 最快：

> library(rbenchmark)
> benchmark(retIdxNA(t1), whichNA(t1), which2(t1), which3(t1), which4(t1), 
+    replications = 10000, order = "relative")[1:4]
          test replications elapsed relative
5   which4(t1)        10000    0.14    1.000
4   which3(t1)        10000    0.16    1.143
1 retIdxNA(t1)        10000    0.17    1.214
2  whichNA(t1)        10000    0.17    1.214
3   which2(t1)        10000    0.25    1.786

对长度为 9000 个元素的向量重复此操作，Armadillo 解决方案的速度比其他解决方案快很多。这里which3（与which4 相同，只是它不为输出向量保留空间）最差，而which4 排在第二位。

> tt <- rep(t1, 1000)
> benchmark(retIdxNA(tt), whichNA(tt), which2(tt), which3(tt), which4(tt), 
+   replications = 1000, order = "relative")[1:4]
          test replications elapsed relative
2  whichNA(tt)         1000    0.09    1.000
5   which4(tt)         1000    0.79    8.778
3   which2(tt)         1000    1.03   11.444
1 retIdxNA(tt)         1000    1.19   13.222
4   which3(tt)         1000   23.58  262.000

Answer 3

以上所有解决方案都是串行的。虽然不是微不足道的，但很有可能利用线程来实现which。有关详细信息，请参阅this write up。虽然对于这么小的尺寸，它不会弊大于利。就像坐飞机走一小段路一样，你在机场安检上浪费了太多时间..

R 通过为与输入一样大的逻辑向量分配内存来实现which，执行单次传递以将索引存储在此内存中，然后最终将其复制到适当的逻辑向量中。

直观地说，这似乎比双循环循环效率低，但不一定，因为复制数据范围很便宜。查看更多详情here。

Answer 4

只需为自己编写一个函数，例如：

which_1<-function(a,b){
return(which(a>b))
}

然后把这个函数传入rcpp。