Как использовать Rcpp для ускорения цикла for?

Question

Я создал цикл for, и я хотел бы ускорить его с помощью библиотеки Rcpp. Я не очень хорошо знаком с C++. Не могли бы вы помочь мне быстрее выполнить мою функцию? Благодарим за помощь!

Я включил свой алгоритм, код вместе с вводом и выходом, с sessionInfo.

Вот мой алгоритм:

, если текущая цена выше предыдущей цены, знак (+1) в колонке под названием TR

, если текущая цена ниже предыдущей цены, марки (-1) в колонке под названием TR

, если текущая цена такая же, как и предыдущая цена, отметить то же самое, как и в предыдущей цены в колонке с именем TR

Вот мой код:

price <- c(71.91, 71.82, 71.81, 71.81, 71.81, 71.82, 71.81, 71.81, 71.81, 
      71.82, 71.81, 71.81, 71.8, 71.81, 71.8, 71.81, 71.8, 71.8, 71.8, 
      71.8, 71.81, 71.81, 71.81, 71.81, 71.81, 71.81, 71.81, 71.81, 
      71.81, 71.82, 71.81, 71.81, 71.81, 71.81, 71.81, 71.81, 71.8, 
      71.8, 71.81, 71.81, 71.81, 71.81, 71.82, 71.82, 71.81, 71.81, 
      71.81, 71.81, 71.81, 71.81, 71.81, 71.82, 71.82, 71.82, 71.82, 
      71.82, 71.82, 71.82, 71.81, 71.82, 71.82, 71.82, 71.82, 71.82, 
      71.82, 71.82, 71.82, 71.82, 71.81, 71.81, 71.81, 71.82, 71.82, 
      71.81, 71.82, 71.82, 71.82, 71.81, 71.82, 71.82, 71.82, 71.81, 
      71.81, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 
      71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 
      71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 
      71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 
      71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 
      71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 71.81, 
      71.82, 71.82, 71.82, 71.82, 71.83, 71.82, 71.82, 71.82, 71.81, 
      71.81, 71.81, 71.81, 71.81, 71.81, 71.81, 71.82, 71.82, 71.82, 
      71.81, 71.81, 71.81, 71.82, 71.82, 71.82, 71.82, 71.82, 71.82, 
      71.82, 71.82, 71.82, 71.82, 71.82, 71.83, 71.83, 71.83, 71.83, 
      71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 
      71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 
      71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 71.83, 
      71.83) 

TR <- numeric(length(price)-1) 
TR <- c(NA,TR) 

for (i in 1: (length(price)-1)){ 

    if (price[i] == price[i+1]) {TR[i+1] = TR[i]} 

    if (price[i] < price[i+1]) {TR[i+1] = 1} 

    if (price[i] > price[i+1]) {TR[i+1] = -1} 

} 

А вот мой выход: dput (TR) дает

c(NA, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 
    -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 
    -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, 
    1, 1, 1, 1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 
    1, 1, -1, 1, 1, 1, -1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 
    -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1) 

и вот мой sessionInfo:

> sessionInfo() 
R version 3.1.2 (2014-10-31) 
Platform: x86_64-w64-mingw32/x64 (64-bit) 

locale: 
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C       
[5] LC_TIME=English_United States.1252  

attached base packages: 
[1] stats  graphics grDevices utils  datasets methods base  

other attached packages: 
[1] data.table_1.9.4 

loaded via a namespace (and not attached): 
[1] chron_2.3-45 plyr_1.8.1 Rcpp_0.11.1 reshape2_1.4 stringr_0.6.2 tools_3.1.2 

Answer 1

Вы можете довольно просто перевести цикл for:

library(Rcpp) 
cppFunction(
"IntegerVector proc(NumericVector x) { 
    const int n = x.size(); 
    IntegerVector y(n); 
    y[0] = NA_INTEGER; 
    for (int i=1; i < n; ++i) { 
    if (x[i] == x[i-1]) y[i] = y[i-1]; 
    else if (x[i] > x[i-1]) y[i] = 1; 
    else y[i] = -1; 
    } 
    return y; 
}") 

Как обычно, вы можете получить довольно большой прирост скорости с Rcpp по сравнению с для цикла в базовой R:

proc.for <- function(price) { 
    TR <- numeric(length(price)-1) 
    TR <- c(NA,TR) 
    for (i in 1: (length(price)-1)){ 
    if (price[i] == price[i+1]) {TR[i+1] = TR[i]} 
    if (price[i] < price[i+1]) {TR[i+1] = 1} 
    if (price[i] > price[i+1]) {TR[i+1] = -1} 
    } 
    return(TR) 
} 
proc.aaron <- function(price) { 
    change <- sign(diff(price)) 
    good <- change != 0 
    goodval <- change[good] 
    c(NA, goodval[cumsum(good)]) 
} 
proc.jbaums <- function(price) { 
    TR <- sign(diff(price)) 
    TR[TR==0] <- TR[which(TR != 0)][findInterval(which(TR == 0), which(TR != 0))] 
    TR 
} 

all.equal(proc(price), proc.for(price), proc.aaron(price), proc.jbaums(price)) 
# [1] TRUE 
library(microbenchmark) 
microbenchmark(proc(price), proc.for(price), proc.aaron(price), proc.jbaums(price)) 
# Unit: microseconds 
#    expr  min  lq  mean median  uq  max neval 
#   proc(price) 1.871 2.5380 3.92111 3.1110 4.5880 15.318 100 
#  proc.for(price) 408.200 448.2830 542.19766 484.1265 546.3255 1821.104 100 
# proc.aaron(price) 23.916 25.5770 33.53259 31.5420 35.8575 190.372 100 
# proc.jbaums(price) 33.536 38.8995 46.80109 43.4510 49.3555 112.306 100 

Мы видим ускорение более чем 100х по сравнению с для цикла и 10x по сравнению с векторизованными альтернативами для предоставленного вектора.

Ускорение является еще более значительным с большим вектором (длина +1000000 здесь испытывались):

price.big <- rep(price, times=5000) 
all.equal(proc(price.big), proc.for(price.big), proc.aaron(price.big), proc.jbaums(price.big)) 
# [1] TRUE 
microbenchmark(proc(price.big), proc.for(price.big), proc.aaron(price.big), proc.jbaums(price.big)) 
# Unit: milliseconds 
#     expr   min   lq  mean  median   uq  max neval 
#   proc(price.big) 1.442119 1.818494 5.094274 2.020437 2.771903 56.54321 100 
#  proc.for(price.big) 2639.819536 2699.493613 2949.962241 2781.636460 3062.277930 4472.35369 100 
# proc.aaron(price.big) 91.499940 99.859418 132.519296 140.521212 147.462259 207.72813 100 
# proc.jbaums(price.big) 117.242451 138.528214 170.989065 170.606048 180.337074 487.13615 100 

Теперь у нас есть 1000x ускорение по сравнению с цикл и ~ 70x ускорение по сравнению с векторизованного R функции. Даже при таком размере неясно, есть много преимуществ Rcpp над векторизованными R-решениями, если функция вызывается только один раз, так как для компиляции кода Rcpp требуется, как минимум, 100 мс. Ускорение довольно привлекательно, если это фрагмент кода, который неоднократно вызывается в вашем анализе.

Answer 2

Это может быть легко vectorised в R.

Например с diff и findInterval:

TR <- sign(diff(price)) 
TR[TR==0] <- TR[which(TR != 0)][findInterval(which(TR == 0), which(TR != 0))] 

Answer 3

Может попробовать векторизации первый. Хотя это, вероятно, будет не так быстро, как Rcpp, это более прямолинейно.

f2 <- function(price) { 
    change <- sign(diff(price)) 
    good <- change != 0 
    goodval <- change[good] 
    c(NA, goodval[cumsum(good)]) 
} 

И это по-прежнему значительно ускоряет работу над циклом R.

f1 <- function(price) { 
    TR <- numeric(length(price)-1) 
    TR <- c(NA,TR) 
    for (i in 1: (length(price)-1)){ 
     if (price[i] == price[i+1]) {TR[i+1] = TR[i]} 
     if (price[i] < price[i+1]) {TR[i+1] = 1} 
     if (price[i] > price[i+1]) {TR[i+1] = -1} 
    } 
    TR 
} 

microbenchmark(f1(price), f2(price), times=100) 
## Unit: microseconds 
##  expr  min  lq  mean median  uq  max neval cld 
## f1(price) 550.037 592.9830 756.20095 618.7910 703.8335 3042.530 100 b 
## f2(price) 36.915 39.3285 56.45267 45.5225 60.1965 184.536 100 a 

Answer 4

Вы можете попробовать выполнить составление байт. Также полезно посмотреть на цикл R, который использует ту же логику if-else-if-else, что и код Rcpp. С R 3.1.2 я получаю

f1 <- function(price) { 
    TR <- numeric(length(price)-1) 
    TR <- c(NA,TR) 
    for (i in 1: (length(price)-1)){ 
     if (price[i] == price[i+1]) {TR[i+1] = TR[i]} 
     if (price[i] < price[i+1]) {TR[i+1] = 1} 
     if (price[i] > price[i+1]) {TR[i+1] = -1} 
    } 
    return(TR) 
} 

f2 <- function(price) { 
    TR <- numeric(length(price)-1) 
    TR <- c(NA,TR) 
    for (i in 1: (length(price)-1)){ 
     if (price[i] == price[i+1]) {TR[i+1] = TR[i]} 
     else if (price[i] < price[i+1]) {TR[i+1] = 1} 
     else {TR[i+1] = -1} 
    } 
    return(TR) 
} 

library(compiler) 
f1c <- cmpfun(f1) 
f2c <- cmpfun(f2) 

library(microbenchmark) 
microbenchmark(f1(price), f2(price), f1c(price), f2c(price), times = 1000) 
## Unit: microseconds 
##  expr  min  lq  mean median  uq  max neval cld 
## f1(price) 536.619 570.3715 667.3520 586.2465 609.9280 45046.462 1000 d 
## f2(price) 328.592 351.2070 386.5895 365.0245 381.4850 1302.497 1000 c 
## f1c(price) 167.570 182.4645 218.9537 192.4780 204.7810 7843.291 1000 b 
## f2c(price) 96.644 107.4465 124.1324 113.5470 121.5365 1019.389 1000 a 

R-Devel, который будет выпущен в качестве R 3.2.0 в апреле, имеет ряд усовершенствований в байтовый код двигателя для скалярных вычислений, как это; там я получаю

microbenchmark(f1(price), f2(price), f1c(price), f2c(price), times = 1000) 
## Unit: microseconds 
##  expr  min  lq  mean median  uq  max neval cld 
## f1(price) 490.300 520.3845 559.19539 533.2050 548.6850 1330.219 1000 d 
## f2(price) 298.375 319.7475 348.71384 330.4535 342.6405 1813.113 1000 c 
## f1c(price) 61.947 66.3255 68.01555 67.7270 69.5470 138.308 1000 b 
## f2c(price) 36.334 38.9500 40.45085 40.1830 41.8610 55.909 1000 a 

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

Все решения различаются тем, как они обрабатывают значения NA/NaN, которые могут или не имеют значения для вас.