Воспроизведение R rep с аргументом times в С++ и Rcpp

Я изучаю использование Rcpp. Я хотел бы использовать С++ для репликации функции rep в R. Rcpp включает в себя несколько функций сахара, которые соответствуют rep в R. (см. Нижнюю часть страницы 3 по адресу: http://cran.r-project.org/web/packages/Rcpp/vignettes/Rcpp-quickref.pdf. Насколько я понимаю документацию, функции сахара rep, rep_each и rep_len принимают два аргумента - вектор и целое число. Однако, что бы я хотел чтобы повторить функцию rep в R, когда я использую аргумент times. В этом случае вы можете предоставить два вектора. Быстрый пример в R:

x <- c(10, 5, 12)
y <- c(2, 6, 3)

rep(x, times = y)
[1] 10 10  5  5  5  5  5  5 12 12 12

Таким образом, rep с аргументом times реплицирует каждый элемент x столько раз, сколько соответствующее значение y. Насколько я понимаю, я не вижу возможности использовать для этого функции сахара Rcpp.

Я создал следующую С++-функцию, которая работает:

// [[Rcpp::export]]
NumericVector reptest(NumericVector x, NumericVector y) {
    int n = y.size();
    NumericVector myvector(sum(y));
    int ind = 0;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < y(i); ++j) {
            myvector(ind) = x[i];
            ind = ind + 1;
            }
        }
    return myvector;
}

x <- c(10, 5, 12)
y <- c(2, 6, 3)

reptest(x, y)
[1] 10 10  5  5  5  5  5  5 12 12 12

Это немного медленнее, чем rep в R. Мне интересно, есть ли способ ускорить это, или если у кого-то есть лучшая идея. Насколько я понимаю, rep вызывает код C, поэтому, возможно, почти невозможно улучшить rep. Моя цель - ускорить цикл MCMC (который использует функцию rep), которая требует много времени для запуска в R, поэтому любое ускорение было бы полезно. Другие части цикла MCMC - это медленные части, а не rep, но мне нужна такая же функциональность в моем цикле.

Ответ 1

Вот быстрый рифф в двух основных версиях. Он также добавляет rep.int():

#include <algorithm>
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector reptest(NumericVector x, NumericVector y) {
    int n = y.size();
    NumericVector myvector(sum(y));
    int ind = 0;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < y[i]; ++j) {
            myvector[ind] = x[i];
            ind = ind + 1;
            }
        }
    return myvector;
}

// [[Rcpp::export]]
NumericVector reptest2(NumericVector x, NumericVector y) {
    int n = y.size();
    NumericVector myvector(sum(y));
    int ind=0;
    for (int i=0; i < n; ++i) {
        int p = y[i];
        std::fill(myvector.begin()+ind, myvector.begin()+ind+p, x[i]);
        ind += p;
    }
    return myvector;
}


/*** R
x <- rep(c(10, 5, 12), 10000)
y <- rep(c(20, 60, 30), 10000)
all.equal(reptest(x, y), reptest2(x, y), rep(x, times=y))

library(microbenchmark)
microbenchmark(reptest(x, y), reptest2(x, y), rep(x, times=y), rep.int(x, y))
***/

С этим мы немного приближаемся, но R все еще выигрывает:

R> Rcpp::sourceCpp("/tmp/rep.cpp")

R> x <- rep(c(10, 5, 12), 10000)

R> y <- rep(c(20, 60, 30), 10000)

R> all.equal(reptest(x, y), reptest2(x, y), rep(x, times=y))
[1] TRUE

R> library(microbenchmark)

R> microbenchmark(reptest(x, y), reptest2(x, y), rep(x, times=y), rep.int(x, y))
Unit: milliseconds
              expr     min      lq    mean  median      uq       max neval
     reptest(x, y) 4.61604 4.74203 5.47543 4.78120 6.78039   7.01879   100
    reptest2(x, y) 3.14788 3.27507 5.25515 3.33166 5.24583 140.64080   100
 rep(x, times = y) 2.45876 2.56025 3.26857 2.60669 4.60116   6.76278   100
     rep.int(x, y) 2.42390 2.50241 3.38362 2.53987 4.56338   6.44241   100
R>

Ответ 2

Одним из способов ускорить это было бы использование std::fill вместо повторения через каждый элемент, который нужно заполнить:

#include <algorithm>
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector reptest2(NumericVector x, NumericVector y) {
  int n = y.size();
  std::vector<double> myvector(sum(y));
  int ind=0;
  for (int i=0; i < n; ++i) {
    std::fill(myvector.begin()+ind, myvector.begin()+ind+y[i], x[i]);
    ind += y[i];
  }
  return Rcpp::wrap(myvector);
}

В более крупном примере это похоже на приближение к rep:

x <- rep(c(10, 5, 12), 10000)
y <- rep(c(20, 60, 30), 10000)
all.equal(reptest(x, y), reptest2(x, y), rep(x, times=y))
# [1] TRUE

library(microbenchmark)
microbenchmark(reptest(x, y), reptest2(x, y), rep(x, times=y))
# Unit: milliseconds
#               expr      min       lq      mean   median        uq      max neval
#      reptest(x, y) 9.072083 9.297573 11.469345 9.522182 13.015692 20.47905   100
#     reptest2(x, y) 5.097358 5.270827  7.367577 5.436549  8.961004 15.68812   100
#  rep(x, times = y) 1.457933 1.499051  2.884887 1.561408  1.949750 13.21706   100

Ответ 3

Мы можем достичь производительности R base rep с помощью no_init:

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
NumericVector reptest3(const NumericVector& x, const IntegerVector& times) {
    std::size_t n = times.size();
    if (n != 1 && n != x.size())
        stop("Invalid 'times' value");
    std::size_t n_out = std::accumulate(times.begin(), times.end(), 0);
    NumericVector res = no_init(n_out);
    auto begin = res.begin();
    for (std::size_t i = 0, ind = 0; i < n; ind += times[i], ++i) {
        auto start = begin + ind;
        auto end = start + times[i];
        std::fill(start, end, x[i]);
    }
    return res;
}

Benchmark:

library(microbenchmark)
x <- rep(c(10, 5, 12), 10000)
y <- rep(c(20, 60, 30), 10000)
microbenchmark(
    reptest(x, y), reptest2(x, y), reptest3(x, y),
    rep(x, times = y), rep.int(x, y))
#> Unit: milliseconds
#>              expr       min        lq      mean    median        uq       max neval
#>     reptest(x, y) 13.209912 14.014886 15.129395 14.457418 15.123676 56.655527   100
#>    reptest2(x, y)  4.289786  4.653088  5.789094  5.105859  5.782284 46.679824   100
#>    reptest3(x, y)  1.812713  2.810637  3.860590  3.194529  3.809141 44.111422   100
#> rep(x, times = y)  2.510219  2.877324  3.576183  3.461315  3.927312  5.961317   100
#>     rep.int(x, y)  2.496481  2.901303  3.422384  3.318761  3.831794  5.283187   100

Также мы можем улучшить этот код с помощью RcppParallel:

struct Sum : Worker {
    const RVector<int> input;
    int value;
    Sum(const IntegerVector& input) : input(input), value(0) {}
    Sum(const Sum& sum, Split) : input(sum.input), value(0) {}
    void operator()(std::size_t begin, std::size_t end) {
        value += std::accumulate(input.begin() + begin, input.begin() + end, 0);
    }
    void join(const Sum& rhs) {
        value += rhs.value;
    }
};

struct Fill: Worker {
    const RVector<double> input;
    const RVector<int> times;
    RVector<double> output;
    std::size_t ind;
    Fill(const NumericVector& input, const IntegerVector& times, NumericVector& output)
        : input(input), times(times), output(output), ind(0) {}
    void operator()(std::size_t begin, std::size_t end) {
        for (std::size_t i = begin; i < end; ind += times[i], ++i)
            std::fill(output.begin() + ind, output.begin() + ind + times[i], input[i]);
    }
};

// [[Rcpp::export]]
NumericVector reptest4(const NumericVector& x, const IntegerVector& times) {
    std::size_t n = times.size();
    if (n != 1 && n != x.size())
        stop("Invalid 'times' value");
    Sum s(times);
    parallelReduce(0, n, s);
    NumericVector res = no_init(s.value);
    Fill f(x, times, res);
    parallelFor(0, n, f);
    return res;
}

Сравнение:

library(microbenchmark)
x <- rep(c(10, 5, 12), 10000)
y <- rep(c(20, 60, 30), 10000)
microbenchmark(
    reptest(x, y), reptest2(x, y), reptest3(x, y),
    rep(x, times = y), rep.int(x, y))
#> Unit: milliseconds
#>               expr      min       lq     mean   median       uq       max neval
#>     reptest3(x, y) 2.442446 3.410985 5.143627 3.893345 5.054285 57.871429   100
#>     reptest4(x, y) 1.211256 1.534428 1.979526 1.821398 2.170999  4.073395   100
#>  rep(x, times = y) 2.435122 3.173904 4.447954 3.795285 4.687695 54.000920   100
#>      rep.int(x, y) 2.444310 3.208522 4.026722 3.913618 4.798793  6.690333   100