Эквивалент Rcpp для rowsum

Вот моя попытка сделать это с помощью Rcpp (впервые использую пакет, поэтому укажите на мою неэффективность):

library(inline)
library(Rcpp)

rowsum_helper = cxxfunction(signature(x = "numeric", y = "integer"), '
  NumericVector var(x);
  IntegerVector factor(y);

  std::vector<double> sum(*std::max_element(factor.begin(), factor.end()) + 1,
                          std::numeric_limits<double>::quiet_NaN());
  for (int i = 0, size = var.size(); i < size; ++i) {
    if (sum[factor[i]] != sum[factor[i]]) sum[factor[i]] = var[i];
    else sum[factor[i]] += var[i];
  }

  return NumericVector(sum.begin(), sum.end());
', plugin = "Rcpp")

rowsum_fast = function(x, y) {
  res = rowsum_helper(x, y)
  elements = which(!is.nan(res))
  list(elements - 1, res[elements])
}

Это довольно быстро для данных примера Мартина, но будет работать только в том случае, если фактор состоит из неотрицательных целых чисел и будет потреблять память порядка наибольшего целого числа в векторе факторов (одно очевидное улучшение вышеизложенного состоит в том, чтобы вычесть min из max до уменьшить использование памяти - что можно сделать либо в функции R, либо в C++).

n = 1e7; x = runif(n); f = sample(n/2, n, T)

system.time(rowsum(x,f))
#    user  system elapsed 
#   14.241  0.170  14.412

system.time({tabulate(f); sum(x)})
#    user  system elapsed 
#   0.216   0.027   0.252

system.time(rowsum_fast(x,f))
#    user  system elapsed 
#   0.313   0.045   0.358

Также обратите внимание, что значительное замедление (по сравнению с tabulate) происходит в коде R, поэтому, если вместо этого вы перенесете его на C++, вы должны увидеть больше улучшений:

system.time(rowsum_helper(x,f))
#    user  system elapsed 
#   0.210   0.018   0.228

Вот обобщение, которое будет обрабатывать почти любой y, но будет немного медленнее (на самом деле я бы предпочел делать это в Rcpp, но не знаю, как там обрабатывать произвольные типы R):

rowsum_fast = function(x, y) {
  if (is.numeric(y)) {
    y.min = min(y)
    y = y - y.min
    res = rowsum_helper(x, y)
  } else {
    y = as.factor(y)
    res = rowsum_helper(x, as.numeric(y))
  }

  elements = which(!is.nan(res))

  if (is.factor(y)) {
    list(levels(y)[elements-1], res[elements])
  } else {
    list(elements - 1 + y.min, res[elements])
  }
}

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

> n = 100000; x = runif(n); f = sample(n/2, n, TRUE)
> system.time(rowsum(x, f))
   user  system elapsed 
  0.228   0.000   0.229 
> n = 1000000; x = runif(n); f = sample(n/2, n, TRUE)
> system.time(rowsum(x, f)) 
   user  system elapsed 
  1.468   0.040   1.514 
> n = 10000000; x = runif(n); f = sample(n/2, n, TRUE)
> system.time(rowsum(x, f))
   user  system elapsed 
 17.369   0.748  18.166 

Кажется, есть два коротких пути, позволяющих избежать повторного заказа

> n = 10000000; x = runif(n); f = sample(n/2, n, TRUE)
> system.time(rowsum(x, f, reorder=FALSE))
   user  system elapsed 
 16.501   0.476  17.025 

и избегая внутреннего принуждения к характеру

> n = 10000000; x = runif(n); f = as.character(sample(n/2, n, TRUE)); 
> system.time(rowsum(x, f, reorder=FALSE))
   user  system elapsed 
  8.652   0.268   8.949 

И затем основные операции, которые, казалось бы, должны быть задействованы — определение уникальных значений фактора группировки (для предварительного выделения результирующего вектора) и выполнение суммы

> n = 10000000; x = runif(n); f = sample(n/2, n, TRUE)
> system.time({ t = tabulate(f); sum(x) })
   user  system elapsed 
  0.640   0.000   0.643 

так что да, похоже, что есть много возможностей для более быстрой одноцелевой реализации. Это кажется естественным для data.table, и его не так уж сложно реализовать на C. Вот смешанное решение, использующее R для табуляции и «классический» интерфейс C для суммирования.

library(inline)

rowsum1.1 <- function(x, f) {
    t <- tabulate(f)
    crowsum1(x, f, t)
}

crowsum1 = cfunction(c(x_in="numeric", f_in="integer", t_in = "integer"), "
    SEXP res_out;
    double *x = REAL(x_in), *res;
    int len = Rf_length(x_in), *f = INTEGER(f_in);

    res_out = PROTECT(Rf_allocVector(REALSXP, Rf_length(t_in)));
    res = REAL(res_out);
    memset(res, 0, Rf_length(t_in) * sizeof(double));
    for (int i = 0; i < len; ++i)
        res[f[i] - 1] += x[i];
    UNPROTECT(1);
    return res_out;
")

с

> system.time(r1.1 <- rowsum1.1(x, f))
   user  system elapsed 
  1.276   0.092   1.373 

Чтобы на самом деле вернуть результат, идентичный rowsum, он должен быть оформлен в виде матрицы с соответствующими тусклыми именами.

rowsum1 <- function(x, f) {
    t <- tabulate(f)
    r <- crowsum1(x, f, t)
    keep <- which(t != 0)
    matrix(r[keep], ncol=1, dimnames=list(keep, NULL))
}

> system.time(r1 <- rowsum1(x, f))
   user  system elapsed 
  9.312   0.300   9.641

так что для всей этой работы мы только в 2 раза быстрее (и гораздо менее общие - x должен быть числовым, f должен быть целым; без значений NA). Да, есть неэффективность, например, выделение уровней пространства, которые не имеют счетчиков (хотя это позволяет избежать дорогостоящего приведения к вектору символов для имен).

Чтобы дополнить код Мартина, вот некоторая версия на основе Rcpp.

int increment_maybe(int value, double vec_i){
    return vec_i == 0 ? value : ( value +1 ) ;  
}

// [[Rcpp::export]]
NumericVector cpprowsum2(NumericVector x, IntegerVector f){
    std::vector<double> vec(10) ;
    vec.reserve(1000); 
    int n=x.size(); 
    for( int i=0; i<n; i++){
        int index=f[i]; 
        while( index >= vec.size() ){
            vec.resize( vec.size() * 2 ) ;    
        }
        vec[ index ] += x[i] ;
    }
    // count the number of non zeros
    int s = std::accumulate( vec.begin(), vec.end(), 0, increment_maybe) ; 
    NumericVector result(s) ;
    CharacterVector names(s) ;

    std::vector<double>::iterator it = vec.begin() ;
    for( int i=0, j=0 ; j<s; j++ ,++it, ++i ){
        // move until the next non zero value
        while( ! *it ){ i++ ; ++it ;}
        result[j] = *it ;
        names[j]  = i ;
    }
    result.attr( "dim" ) = IntegerVector::create(s, 1) ;
    result.attr( "dimnames" ) = List::create(names, R_NilValue) ; 
    return result ;
}

Код C++ имеет дело со всем, включая форматирование в матричный формат, заданный rowsum< /a> и показывает (немного) лучшую производительность (по крайней мере, в примере).

# from Martin's answer
> system.time(r1 <- rowsum1(x, f))
   user  system elapsed
  0.014   0.001   0.015

> system.time(r3 <- cpprowsum2(x, f))
   user  system elapsed
  0.011   0.001   0.013

> identical(r1, r3)
[1] TRUE

В комментарии и «ответе», которые @Ben удалил, оказывается, что f упорядочено и увеличивается.

n = 1e7; x = runif(n);
f <- cumsum(c(1L, sample(c(TRUE, FALSE), n - 1, TRUE)))

So

rowsum3 <- function(x, f)
{
    y <- cumsum(x)
    end <- c(f[-length(f)] != f[-1], TRUE)
    diff(c(0, y[end]))
}

является общим решением R (если кто-то не слишком беспокоится о точности), и

crowsum3 <- cfunction(c(x_in="numeric", f_in="integer"), "
    int j = 0, *f = INTEGER(f_in), len = Rf_length(f_in), 
        len_out = len == 0 ? 0 : f[len - 1];
    SEXP res = Rf_allocVector(REALSXP, len_out);
    double *x = REAL(x_in), *r = REAL(res);
    memset(r, 0, len_out * sizeof(double));
    for (int i = 0; i < len; ++i) {
        if (i != 0 && f[i] != f[i-1]) ++j;
        r[j] += x[i];
    }
    return res;
")

может быть решением C. У них есть тайминги

> system.time(r3 <- rowsum3(x, f))
   user  system elapsed 
  1.116   0.120   1.238 
> system.time(c3 <- crowsum3(x, f))
   user  system elapsed 
  0.080   0.000   0.081 

и потеря точности в реализации R очевидна

> all.equal(r3, c3)
[1] TRUE
> identical(r3, c3)
[1] FALSE

rowsum_helper есть

> system.time(r2 <- rowsum_helper(x, f))
   user  system elapsed 
  0.464   0.004   0.470 

но также предполагает индексацию на основе 0, поэтому

> head(rowsum_helper(x, f))
[1]       NaN 0.9166577 0.4380485 0.7777094 2.0866507 0.7300764
> head(crowsum3(x, f))
[1] 0.9166577 0.4380485 0.7777094 2.0866507 0.7300764 0.7195091