programacion - funciones r en r ejemplos
¿Cómo usar Rcpp para acelerar un bucle for? (4)
Esto puede ser fácilmente vectorizado en R.
Por ejemplo, con diff y findInterval :
TR <- sign(diff(price))
TR[TR==0] <- TR[which(TR != 0)][findInterval(which(TR == 0), which(TR != 0))]
He creado un ciclo for y me gustaría acelerarlo usando la biblioteca Rcpp. No estoy muy familiarizado con C ++. ¿Puedes ayudarme a hacer mi función más rápido? ¡Gracias por tu ayuda!
He incluido mi algoritmo, el código, junto con la entrada y la salida, con sessionInfo.
Aquí está mi algoritmo:
si el precio actual está por encima del precio anterior, marque (+1) en la columna llamada TR
si el precio actual está por debajo del precio anterior, marque (-1) en la columna llamada TR
si el precio actual es el mismo que el precio anterior, marque lo mismo que en el precio anterior en la columna llamada TR
Aquí está mi código:
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}
}
Y aquí está mi resultado: rendimientos 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)
y aquí está mi 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
Puede traducir directamente el bucle 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;
}")
Como es habitual, puede obtener una aceleración bastante grande con Rcpp en comparación con el bucle for en la base 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
Vemos una aceleración de más de 100x en comparación con un bucle for y 10x en comparación con las alternativas vectorizadas para el vector proporcionado.
La aceleración es aún más significativa con un vector más grande (longitud 1 millón probado aquí):
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
Ahora tenemos una aceleración de 1000x en comparación con el bucle for y una velocidad de ~ 70x en comparación con las funciones R vectorizadas. Incluso a este tamaño, no está claro que haya mucha ventaja de Rcpp sobre las soluciones R vectorizadas si la función solo se llama una vez, ya que ciertamente lleva al menos 100 ms compilar el código Rcpp. La aceleración es bastante atractiva si se trata de una pieza de código que se llama repetidamente en su análisis.
Puedes probar la compilación de bytes. También es útil observar un bucle R que usa la misma lógica if-else-if-else que el código Rcpp. Con R 3.1.2 consigo
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, que se lanzará como R 3.2.0 en abril, tiene una serie de mejoras en el motor de código de bytes para cálculos escalares como este; ahí me pongo
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
Esto te lleva al mismo estadio general que las soluciones vectorizadas en este ejemplo. Todavía hay lugar para nuevas mejoras en el motor de código de bytes que debería convertirse en versiones futuras.
Todas las soluciones difieren en la forma en que manejan los valores de NA/NaN , que pueden o no ser importantes para usted.
Quizás intente la vectorización primero. Aunque probablemente no sea tan rápido como Rcpp, es más sencillo.
f2 <- function(price) {
change <- sign(diff(price))
good <- change != 0
goodval <- change[good]
c(NA, goodval[cumsum(good)])
}
Y sigue siendo una velocidad sustancial sobre un R para bucle.
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