在 C++ 和 Rcpp 中使用 times 参数重现 R rep
Reproducing R rep with the times argument in C++ and Rcpp
我正在学习使用 Rcpp。我想使用 C++ 复制 R 中的 rep 函数。Rcpp 包含几个与 R 中的 rep 相对应的糖函数。（请参阅第 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 糖函数的任何方法。
我创建了以下有效的 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
它比 R 中的 rep 慢一点。我想知道是否有任何方法可以加快速度，或者是否有人有更好的主意。据我了解，rep 正在调用 C 代码，因此可能几乎不可能改进 rep。我的目标是加速 MCMC 循环（使用 rep 函数），它在 R 中需要很多时间才能达到 运行，因此任何加速都会很有用。 MCMC 循环的其他部分是慢速部分，不是 rep，但我的循环中需要相同的功能。
一种加快速度的方法是使用 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
这是两个初始版本的快速重复。它还添加了 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>
我们可以通过 no_init:
实现 R base rep 的性能
// [[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;
}
基准：
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
我正在学习使用 Rcpp。我想使用 C++ 复制 R 中的 rep 函数。Rcpp 包含几个与 R 中的 rep 相对应的糖函数。（请参阅第 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 糖函数的任何方法。
我创建了以下有效的 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
它比 R 中的 rep 慢一点。我想知道是否有任何方法可以加快速度，或者是否有人有更好的主意。据我了解，rep 正在调用 C 代码，因此可能几乎不可能改进 rep。我的目标是加速 MCMC 循环（使用 rep 函数），它在 R 中需要很多时间才能达到 运行，因此任何加速都会很有用。 MCMC 循环的其他部分是慢速部分，不是 rep，但我的循环中需要相同的功能。
一种加快速度的方法是使用 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
这是两个初始版本的快速重复。它还添加了 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>
我们可以通过 no_init:
rep 的性能
// [[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;
}
基准：
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