Tutorial Rcpp: Cómo usar C++ en R con sólo 3 líneas!

George Vega
23 Abr 2014

Grupo de Usuarios de R en Chile

Usua

Agenda

En esta presentación…

  • Introducción
    • Requerimientos
    • Verdades de Rcpp
    • Patada inicial!
  • Cómo utilizar Rcpp
    • Cómo 'compilar' con Rcpp
    • Definición del lenguaje
    • Ejemplos

Requerimientos

Para poder utilizar Rcpp es necesario contar con lo siguiente:

  • R (>= 3.0.0)
  • gcc (Incluido en Rtools)

Adicionalmente se sugiere

Verdades de Rcpp

Algunos motivos para utilizar Rcpp?

Calentando motores...

Nos saltaremos la parte del típico 'Hola mundo' para pasar de lleno a Rcpp

Para eso partiremos con un ejemplo propuesto por uno de los creadores de Rcpp, Dirk Eddelbuettel

Considere la siguiente función

\[ f(n) = \begin{cases} n & \text{si $n<2$} \\ f(n-1) + f(n - 2) & \text{si $n\geq2$} \end{cases} \]

Su implementación en R sería

# Definiendo la funcion
f <- function(n) {
  if (n<2) return(n)

  return(f(n-1) + f(n-2))
}

# Aplicando la funcion a un vector de 0 a 9
sapply(0:9, f)

 [1]  0  1  1  2  3  5  8 13 21 34

Y su equivalente en Rcpp

# Cargando Rcpp
library(Rcpp)

# Definiendo la funcion
cppFunction('
int g(int n) {
  if (n < 2) return(n);
  return(g(n-1) + g(n-2));
}'
)

# Aplicando la funcion a un vector de 0 a 9
sapply(0:9, g)

 [1]  0  1  1  2  3  5  8 13 21 34

Y por qué la molestia?

Misma función… diferente desempeño!

library(rbenchmark)
benchmark(f(20), g(20))[,1:4]

   test replications elapsed relative
1 f(20)          100   2.673    445.5
2 g(20)          100   0.006      1.0
  • Rcpp es al menos x400 más rápido!

Rcpp vs R API

Lo que en Rcpp implica acceder a una función del paquete stats …

Environment stats("package:stats");
Function rnorm = stats["rnorm"];
return rnorm(10, Named("sd", 100.0));

Rcpp vs R API

… en la API de R equivale a un desastre!

SEXP stats = PROTECT(
  R_FindNamespace(
    mkString("stats")));
SEXP rnorm = PROTECT(
  findVarInFrame(stats,
    install("rnorm")));
SEXP call = PROTECT(
  LCONS( rnorm,
    CONS(ScalarInteger(10),
      CONS(ScalarReal(100.0),
        R_NilValue))));
SET_TAG(CDDR(call),install("sd"));
SEXP res = PROTECT(eval(call,
  R_GlobalEnv));
UNPROTECT(4);
return res;

Maneras de utilizar `Rcpp`

El paquete Rcpp permite compilar código a través de tres funciones clave:

  • cppFunction Define funciones de Rcpp de manera dinámica.
  • sourceCpp Ejecuta archivos .cpp, evalúa e importa funciones (además permite incorporar bloques de R).
  • evalCpp Evalúa expresiones escritas en C++

(Nota: En Rstudio presionar [Ctrl] + [Enter] equivale a sourceCpp)

Ejemplos evalCpp:

Como cppFunction ya lo vimos, partiremos con evalCpp

# Obteniendo constantes
evalCpp("PI")

[1] 3.142

evalCpp("std::numeric_limits<double>::max()")

[1] 1.798e+308

evalCpp("__cplusplus")

[1] 199711

Ejemplos con `sourceCpp`

Utilizando la etiqueta // [[Rcpp::export]], sourceCpp importa funciones escritas en C++ a R de manera dinámica. Más aún, permite incluir bloques de código de R.

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
int g(int n) {
  if (n < 2) return(n);
  return(g(n-1) + g(n-2));
}

/*** R
g(10)
***/

Tipos de Objetos en `Rcpp`

En lo que sigue, revisaremos las clases de objetos en Rcpp y como manipularlos.

  • Vectores IntegerVector, NumericVector, LogicalVector, CharacterVector (Se accede con [i]).
  • Matrices IntegerMatrix, NumericMatrix, … (Se accede con (i,j)).
  • Listas List (Se accede con _["objeto"]).
  • Funciones Function

Ejemplo con NumericVector (R 1)

x <- -5:5
x
log(abs(x))

 [1] -5 -4 -3 -2 -1  0  1  2  3  4  5

 [1] 1.6094 1.3863 1.0986 0.6931 0.0000   -Inf 0.0000 0.6931 1.0986 1.3863
[11] 1.6094

Ejemplo con NumericVector (R 2)

rla <- function(x) {
  for (i in 1:length(x)) {
    x[i] <- log(abs(x[i]))
  }
  return(x)
}
rla(x)

 [1] 1.6094 1.3863 1.0986 0.6931 0.0000   -Inf 0.0000 0.6931 1.0986 1.3863
[11] 1.6094

Ejemplo con NumericVector (Rcpp 1)

cppFunction("
NumericVector cppla1(NumericVector x) {
  NumericVector out(x.size());
  for(int i=0;i<x.size();i++)
    out[i] = log(fabs(x[i]));
  return out;
}")

cppla1(x)

 [1] 1.6094 1.3863 1.0986 0.6931 0.0000   -Inf 0.0000 0.6931 1.0986 1.3863
[11] 1.6094

Ejemplo con NumericVector (Rcpp 2)

Utilizando sugar podemos escribirlo de una forma más sencilla!

cppFunction(" 
NumericVector cppla2(NumericVector x) {
  return log(abs(x)); 
}")

cppla1(x)

 [1] 1.6094 1.3863 1.0986 0.6931 0.0000   -Inf 0.0000 0.6931 1.0986 1.3863
[11] 1.6094

Algunas de las funciones incorporadas son abs(x), exp(x), floor(x), ceil(x), pow(x, z). (más en breve)

Ejemplo con NumericVector (R vs Rcpp)

x <- 1:10000
rbenchmark::benchmark(log(abs(x)), rla(x),cppla1(x),cppla2(x))[,1:4]

         test replications elapsed relative
3   cppla1(x)          100   0.061    1.271
4   cppla2(x)          100   0.048    1.000
1 log(abs(x))          100   0.056    1.167
2      rla(x)          100   2.057   42.854

Ejemplo con NumericMatrix (R)

mimat <- matrix(1:9, nrow=3)
mimat

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9


diag(mimat)

[1] 1 5 9

Ejemplo con NumericMatrix (Rcpp)

cppFunction("
NumericVector cppdiag(NumericMatrix x) {
  NumericVector out(x.nrow());
  for(int i=0;i<x.nrow();i++)
    out[i] = x(i,i);
  return out;
}")

cppdiag(mimat)

[1] 1 5 9

Ejemplo con List (R)

Crearemos una lista con números aleatorios un vector char y un vector 'al cuadrado'

rlist <- function(x,y,z) {
    return(list(
      x=rnorm(x),
      y=y,
      z=z^2
      ))
}

Ejemplo con List (Rcpp)

Su equivalente en Rcpp (una vez mas incluyendo sugar)

cppFunction('
List cpplist(
  int x, CharacterVector y, NumericVector z) {
  return List::create(
    _["x"]=rnorm(x,0,1),
    _["y"]=y,
    _["z"]=pow(z,2) /* sugar */
    );
}')

Ejemplo con List (R vs Rcpp)

set.seed(123)
rlist(3, c("a","hola","Si"),
      seq(-3,3))

$x
[1] -0.5605 -0.2302  1.5587

$y
[1] "a"    "hola" "Si"  

$z
[1] 9 4 1 0 1 4 9

set.seed(123)
cpplist(3, c("a","hola","Si"),
        seq(-3,3))

$x
[1] -0.5605 -0.2302  1.5587

$y
[1] "a"    "hola" "Si"  

$z
[1] 9 4 1 0 1 4 9

Ejemplo de funciones (R)

# Funcion que llama funciones
rf <- function(x,f) {
  for(i in 1:100)
    out <- f(x) 
  return(out)
}

# Ejemplos
rf(c(1,4,9), sqrt)
rf(c(1,4,9), function(x) 1/x)

[1] 1 2 3

[1] 1.0000 0.2500 0.1111

Ejemplo de funciones (Rcpp)

# Funcion que llama funciones
cppFunction('
NumericVector cppf(NumericVector x, Function f) {
  NumericVector out(x.length());
  for(int i=0;i<100;i++)
    out=f(x);
  return out; 
}')

# Ejemplos
cppf(c(1,4,9), sqrt)
cppf(c(1,4,9), function(x) 1/x)

[1] 1 2 3

[1] 1.0000 0.2500 0.1111

Ejemplo de funciones (R vs Rcpp)

Con Rcpp es posible acelerar los tiempos de cómputo de funciones escritas por el usuario en bucles.

x <- 1:3
f <- sqrt
benchmark(rf(x,f), cppf(x, f))[,1:4]

        test replications elapsed relative
2 cppf(x, f)          100   0.148       37
1   rf(x, f)          100   0.004        1

Genial! Ahora... donde buscar más info?

Muchas Gracias!

George Vega (g.vegayon@gmail.com)

http://ggvega.com

Grupo de Usuarios de R en Chile, 23 Abr 2014

(presentación creada con Rstudio + knitr)

Usua

Referencias