C++ ve Rcpp

Question

'daki zaman argümanı ile R rep reprodüksiyonu Rcpp kullanmayı öğreniyorum. . Ben R. Rcpp içinde rep fonksiyon R. içinde rep (en sayfa 3 alt kısmını görmek için gelen birkaç şeker işlevleri içerir çoğaltmak için C++ kullanmak istiyorum: http://cran.r-project.org/web/packages/Rcpp/vignettes/Rcpp-quickref.pdf Ben belgeleri anlamak gibi şeker fonksiyonları rep, rep_each, ve rep_len iki argüman alır - bir vektör ve bir tamsayı.Ama yapmak istediğim, times argümanını kullandığımda R içinde rep'un işlevselliğini çoğaltmaktır.Bu durumda, iki vektör sağlayabilirsiniz. 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 

Böylece

times bağımsız değişken ile repher öğe çoğaltır, karşılık gelen y değeri kadar çok kez. Anladığım kadarıyla, bunun için Rcpp şeker fonksiyonlarını kullanmanın bir yolunu göremiyorum.

// [[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 

Ben bu hızlandırmak için zaten var olup olmadığını merak ya da herhangi bir iyi fikri varsa am R. içinde rep biraz daha yavaştır:

ben çalışır aşağıdaki C++ işlevi oluşturduk. Anladığım kadarıyla, rep C kodunu çağırıyor, bu yüzden rep'u iyileştirmek neredeyse imkansız olacaktır. Amacım, bir MCMC döngüsünü hızlandırarak (rep işlevini kullanıyor) R'de çalışmak için çok zaman harcıyor, böylece herhangi bir hızlanma yararlı olacaktır. MCMC döngüsünün diğer bölümleri rep değil, yavaş kısımlardır, ancak döngümde aynı işlevselliğe ihtiyacım var. Bunu hızlandırmak için

Answer 1

İki intial versiyonda hızlı bir riff. Ayrıca rep.int() ekler: Bununla

#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)) 
***/ 

, biz biraz daha yakın olsun ama R hala kazanır:

// [[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; 
} 

:

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> 

Answer 2

bir yolu yerine doldurulacak her bir unsur ile yinelemek std::fill kullanmak olacaktır: Daha büyük bir örnek üzerinde

#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); 
} 

, bu rep yaklaşmak için görünü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)) 
# [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 

Answer 3

Biz no_init R tabanı rep performans elde edebilirsiniz Karşılaştırma:

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 

Ayrıca RcppParallel ile bu kodu iyileştirebilir:

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; 
} 

Karşılaştırma:

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