我是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，有没有简单的组合方法？
最新版本的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++中的off-by-one模块，因为它们是从零开始的）:
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中也有些新。
请尝试以下操作:
#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元素长的矢量重复这一步，犰狳解比其他解来得快得多。 这里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
以上所有的解决方案都是串行的。 虽然不是很简单，但是利用线程实现哪个是很可能的。 更多细节，请参阅本文。 虽然对于这么小的尺寸来说，也不会是弊大于利。 就像乘飞机飞一小段距离一样，你在机场安检处浪费了太多时间。
R实现了，通过为与输入一样大的逻辑向量分配内存，执行一次传递以将索引存储在该内存中，然后将其最终复制到正确的逻辑向量中。
直观地看，这似乎比双通循环效率低，但也不一定，因为复制一个数据范围很便宜。 请在此查看更多详细信息。
