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组合起来？

- uday

---

请参见 https://github.com/romainfrancois/Rcpp-book/blob/master/chapters/threading/which.md 以获取 which 的线程版本。 - Romain Francois

@RomainFrancois，感谢分享链接。 - uday

很高兴您觉得这很有用。我知道它需要重写，但基本思路已经存在。 - Romain Francois

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中，逻辑子集也是比较新的东西。

- Dirk Eddelbuettel

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

- G. Grothendieck

---

当然可以，但为什么不使用我发布的非循环两行代码呢？ 如果您认为有必要，可以随意向序列添加+1偏移量。 - Dirk Eddelbuettel

如果在创建向量时可以指定空间的预留，那么就不需要进行复制。 - G. Grothendieck

我正在寻找一些可以分配空间但不会使向量变得过大的东西，以便将向量扩展到该大小时不会产生显着影响。data.table有类似的功能。R中的data.table具有长度（列数）和truelength（可以扩展到的列数，而不会产生显着影响）。C++中的std::vector使用reserve实现此功能。 - G. Grothendieck

实际上，我相信R向量除了实际长度之外还有一个真实长度，尽管无法从R中访问真实长度。 - G. Grothendieck

@G.Grothendieck，“真实长度”这个东西并不像名字所暗示的那样。它主要用于R中的哈希表，但与名称所暗示的最大容量无关。 - Romain Francois

Answer 3

以上所有解决方案都是串行的。虽然不是微不足道的，但是完全可以利用线程来实现 which。有关更多详细信息，请参见此处的文档。尽管对于这样的小规模数据，使用线程不会带来更多好处，就像为了短距离乘坐飞机一样，您会在机场安检处浪费太多时间。

R通过为与输入大小相同的逻辑向量分配内存，进行单次遍历以将索引存储在此内存中，然后最终将其复制到适当的逻辑向量中来实现which。 直觉上，这似乎不如双重循环高效，但并非必然如此，因为复制数据范围很便宜。请参见此处的更多详情。

- Romain Francois

---

这里有一些代码：https://github.com/romainfrancois/Rcpp-book/tree/master/chapters/threading/code 。由于代码已经4个月没有更新了，可能需要进行更新。请注意，代码需要使用 <thread> 标头文件，而据我所知，当前的 Rtools windows 工具集中并未提供该文件。 - Romain Francois

这一切都非常好，你可能会喜欢Luke在我们的会议上发表的主题演讲，因为在R中也会有基于OpenMP的类似方法“在某个时候”，但这完全与OP所要求的解决方案无关：一个可以在今天使用的Rcpp（而不是Rcpp11）的解决方案。 - Dirk Eddelbuettel

你今天可以使用Rcpp11。 - Romain Francois

你现在也可以在Rcpp中使用C++11。但是你仍然没有给出一个OP可以使用的答案。 - Dirk Eddelbuettel

@DirkEddelbuettel，Rcpp和Rcpp11之间的主要区别是什么？ - uday

@uday 主要的区别在于谁掌握主导权，以及因此每个项目的演进方式。Rcpp11 是一次完整的重写，只保留最好的部分，并专注于现代 C++ 和创新。另一方面，Rcpp 是稳定的，并兼容旧版 C++ 标准。 - Romain Francois

Answer 4

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

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

然后将此函数传递到rcpp中。

- guagua

---

欢迎来到Stack Overflow :) - www139