Эквивалент функции 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.

Другие вопросы по тегам