Эквивалент функции which в Rcpp
Я новичок в 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) есть ли эквивалент which
в рспп?
(ii) есть ли способ сделать вышеуказанную функцию короче / четче? В частности, есть ли простой способ объединить вышеприведенные шаги 1, 2, 3?
4 ответа
Последние версии 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)
*/
дает желаемый ответ (модуль off-by-one в 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.
Попробуй это:
#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
Все вышеперечисленные решения являются серийными. Хотя это и не тривиально, но вполне возможно использовать преимущества многопоточности для реализации which
, Смотрите это написать для более подробной информации. Хотя для таких небольших размеров это принесло бы больше вреда, чем пользы. Как если бы вы взяли самолет на небольшое расстояние, вы слишком много времени теряли в безопасности аэропорта.
R реализует which
выделяя память для логического вектора размером с вход, делает один проход для сохранения индексов в этой памяти, а затем копирует его в соответствующий логический вектор.
Интуитивно это кажется менее эффективным, чем двухпроходный цикл, но не обязательно, так как копирование диапазона данных обходится дешево. Смотрите подробности здесь.
Просто напишите для себя функцию, например:
which_1<-function(a,b){
return(which(a>b))
}
Затем передайте эту функцию в rcpp.