在 C++ 和 Rcpp 中使用时间参数再现 R rep

Question

我正在学习使用 Rcpp。我想用 C++ 来复制 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 一样多值(value)。据我了解，我看不到为此使用 Rcpp 糖函数的任何方法。

我创建了以下有效的 C++ 函数:

// [[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 .我的目标是加速需要大量时间在 R 中运行的 MCMC 循环(使用 rep 函数)，因此任何加速都会有用。 MCMC 循环的其他部分是慢速部分，而不是 rep ，但我的循环中需要相同的功能。

Answer 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>