Функция быстрого статистического режима Rcpp с векторным вводом

Чтобы заставить функцию работать для любого векторного ввода, вы можете реализовать алгоритм @ JosephWood для любого типа данных, который вы хотите поддерживать, и вызвать его из switch(TYPEOF(x)). Но это было бы много дублирования кода. Вместо этого лучше создать универсальную функцию, которая может работать с любым аргументом Vector<RTYPE>. Если мы следуем парадигме R, что все является вектором, и пусть функция также возвращает a Vector<RTYPE>, то мы можем использовать RCPP_RETURN_VECTOR . Обратите внимание, что нам нужен C ++ 11, чтобы иметь возможность передавать дополнительные аргументы в функцию, вызываемую RCPP_RETURN_VECTOR. Одна хитрость в том, что вам нужен тип хранилища для Vector<RTYPE>, чтобы создать подходящий std::unordered_map. Здесь Rcpp::traits::storage_type<RTYPE>::type приходит на помощь. Однако std::unordered_map не знает, как обращаться с комплексными числами из R. Для простоты я отключаю этот особый случай.

Собираем все вместе:

#include <Rcpp.h>
using namespace Rcpp ;

// [[Rcpp::plugins(cpp11)]]
#include <unordered_map>

template <int RTYPE>
Vector<RTYPE> fastModeImpl(Vector<RTYPE> x, bool narm){
  if (narm) x = x[!is_na(x)];
  int myMax = 1;
  Vector<RTYPE> myMode(1);
  // special case for factors == INTSXP with "class" and "levels" attribute
  if (x.hasAttribute("levels")){
    myMode.attr("class") = x.attr("class");
    myMode.attr("levels") = x.attr("levels");
  }
  std::unordered_map<typename Rcpp::traits::storage_type<RTYPE>::type, int> modeMap;
  modeMap.reserve(x.size());

  for (std::size_t i = 0, len = x.size(); i < len; ++i) {
    auto it = modeMap.find(x[i]);

    if (it != modeMap.end()) {
      ++(it->second);
      if (it->second > myMax) {
        myMax = it->second;
        myMode[0] = x[i];
      }
    } else {
      modeMap.insert({x[i], 1});
    }
  }

  return myMode;
}

template <>
Vector<CPLXSXP> fastModeImpl(Vector<CPLXSXP> x, bool narm) {
  stop("Not supported SEXP type!");
}

// [[Rcpp::export]]
SEXP fastMode( SEXP x, bool narm = false ){
  RCPP_RETURN_VECTOR(fastModeImpl, x, narm);
}

/*** R
set.seed(1234)
s <- sample(1e5, replace = TRUE)
fastMode(s)
fastMode(s + 0.1)
l <- sample(c(TRUE, FALSE), 11, replace = TRUE) 
fastMode(l)
c <- sample(letters, 1e5, replace = TRUE)
fastMode(c)
f <- as.factor(c)
fastMode(f) 
*/

Выход:

> set.seed(1234)

> s <- sample(1e5, replace = TRUE)

> fastMode(s)
[1] 85433

> fastMode(s + 0.1)
[1] 85433.1

> l <- sample(c(TRUE, FALSE), 11, replace = TRUE) 

> fastMode(l)
[1] TRUE

> c <- sample(letters, 1e5, replace = TRUE)

> fastMode(c)
[1] "z"

> f <- as.factor(c)

> fastMode(f) 
[1] z
Levels: a b c d e f g h i j k l m n o p q r s t u v w x y z

В вашей функции Mode, поскольку вы в основном вызываете функции-оболочки, вы не увидите такого большого улучшения по сравнению с базой R. На самом деле, просто написав точный базовый перевод R, мы имеем:

baseMode <- function(x, narm = FALSE) {
    if (narm) x <- x[!is.na(x)]
    ux <- unique(x)
    ux[which.max(table(match(x, ux)))]
}

И для сравнительного анализа имеем:

set.seed(1234)
s <- sample(1e5, replace = TRUE)

library(microbenchmark)
microbenchmark(Mode(s), baseMode(s), times = 10, unit = "relative")
Unit: relative
       expr      min       lq     mean   median       uq      max neval
    Mode(s) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    10
baseMode(s) 1.490765 1.645367 1.571132 1.616061 1.637181 1.448306    10

Как правило, когда мы предпринимаем усилия по написанию собственного скомпилированного код, мы ожидали бы большие выгоды. Простое объединение этих и без того эффективных скомпилированных функций в Rcpp не сможет волшебным образом принести вам ожидаемые результаты. На самом деле, на больших примерах базовое решение быстрее. Заметьте:

set.seed(1234)
sBig <- sample(1e6, replace = TRUE)

system.time(Mode(sBig))
 user  system elapsed 
1.410   0.036   1.450 

system.time(baseMode(sBig))
 user  system elapsed 
0.915   0.025   0.943 

Чтобы ответить на ваш вопрос о написании функции более быстрого режима, мы можем использовать std::unordered_map, который очень похож на table под капотом (т.е. они обе являются хеш-таблицами в их сердце). Кроме того, поскольку вы возвращаете одно целое число, мы можем с уверенностью предположить, что мы можем заменить NumericVector на IntegerVector, а также то, что вам не нужно возвращать каждое значение , которое встречается чаще всего.

Алгоритм, приведенный ниже, можно изменить, чтобы он возвращал истинный режим , но я оставлю это в качестве упражнения (подсказка: вам потребуется std::vector вместе с выполнением какого-либо действия, когда it->second == myMax) , Нотабене вам также нужно добавить // [[Rcpp::plugins(cpp11)]] вверху вашего файла cpp для std::unordered_map и auto.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::plugins(cpp11)]]
#include <unordered_map>

// [[Rcpp::export]]
int fastIntMode(IntegerVector x, bool narm = false) {
    if (narm) x = x[!is_na(x)];
    int myMax = 1;
    int myMode = 0;
    std::unordered_map<int, int> modeMap;
    modeMap.reserve(x.size());

    for (std::size_t i = 0, len = x.size(); i < len; ++i) {
        auto it = modeMap.find(x[i]);

        if (it != modeMap.end()) {
            ++(it->second);
            if (it->second > myMax) {
                myMax = it->second;
                myMode = x[i];
            }
        } else {
            modeMap.insert({x[i], 1});
        }
    }

    return myMode;
}

И тесты:

microbenchmark(Mode(s), baseMode(s), fastIntMode(s), times = 15, unit = "relative")
Unit: relative
          expr      min       lq     mean   median       uq      max neval
       Mode(s) 6.428343 6.268131 6.622914 6.134388 6.881746  7.78522    15
   baseMode(s) 9.757491 9.404101 9.454857 9.169315 9.018938 10.16640    15
fastIntMode(s) 1.000000 1.000000 1.000000 1.000000 1.000000  1.00000    15

Теперь мы говорим ... примерно в 6 раз быстрее, чем оригинал, и в 9 раз быстрее, чем база. Все они возвращают одно и то же значение:

fastIntMode(s)
##[1] 85433

baseMode(s)
##[1] 85433

Mode(s)
##[1] 85433

И для нашего более крупного примера:

## base R returned in 0.943s
system.time(fastIntMode(s))
 user  system elapsed 
0.217   0.006   0.224