Rcpp.qmd

# Reescribiendo código de R en C++ {#sec-rcpp}

```{r, include = FALSE}
source("common.R")
```

## Introducción

A veces, el código R simplemente no es lo suficientemente rápido. Ha utilizado la creación de perfiles para descubrir dónde están los cuellos de botella y ha hecho todo lo posible en R, pero su código aún no es lo suficientemente rápido. En este capítulo, aprenderá a mejorar el rendimiento reescribiendo funciones clave en C++. Esta magia viene a través del paquete [Rcpp](http://www.rcpp.org/) [@Rcpp] (con contribuciones clave de Doug Bates, John Chambers y JJ Allaire).

Rcpp hace que sea muy simple conectar C++ a R. Si bien es *posible* escribir código C o Fortran para usar en R, será doloroso en comparación. Rcpp proporciona una API limpia y accesible que le permite escribir código de alto rendimiento, aislado de la compleja API C de R. \index{Rcpp} \index{C++}

Los cuellos de botella típicos que C++ puede abordar incluyen:

-   Bucles que no se pueden vectorizar fácilmente porque las iteraciones posteriores dependen de las anteriores.

-   Funciones recursivas, o problemas que implican llamar funciones millones de veces. La sobrecarga de llamar a una función en C++ es mucho menor que en R.

-   Problemas que requieren estructuras de datos y algoritmos avanzados que R no proporciona. A través de la biblioteca de plantillas estándar (STL), C++ tiene implementaciones eficientes de muchas estructuras de datos importantes, desde mapas ordenados hasta colas de dos extremos.

El objetivo de este capítulo es discutir solo aquellos aspectos de C++ y Rcpp que son absolutamente necesarios para ayudarlo a eliminar los cuellos de botella en su código. No dedicaremos mucho tiempo a funciones avanzadas como la programación orientada a objetos o las plantillas porque el enfoque está en escribir funciones pequeñas e independientes, no grandes programas. Un conocimiento práctico de C++ es útil, pero no esencial. Muchos buenos tutoriales y referencias están disponibles gratuitamente, incluidos <http://www.learncpp.com/> y <https://en.cppreference.com/w/cpp>. Para temas más avanzados, la serie *Effective C++* de Scott Meyers es una opción popular.

### Estructura {.unnumbered}

-   La @sec-rcpp-intro teaches you how to write C++ by converting simple R functions to their C++ equivalents. You'll learn how C++ differs from R, and what the key scalar, vector, and matrix classes are called.

-   La @sec-sourceCpp le muestra cómo usar `sourceCpp()` para cargar un archivo C++ desde el disco de la misma manera que usa `source()` para cargar un archivo de código R.

-   La @sec-rcpp-classes analiza cómo modificar los atributos de Rcpp y menciona algunas de las otras clases importantes.

-   La @sec-rcpp-na le enseña cómo trabajar con los valores faltantes de R en C++.

-   La @sec-stl le muestra cómo usar algunas de las estructuras de datos y algoritmos más importantes de la biblioteca de plantillas estándar, o STL, integrada en C++.

-   La @sec-rcpp-case-studies muestra dos estudios de casos reales en los que se utilizó Rcpp para obtener mejoras de rendimiento considerables.

-   La @sec-rcpp-package le enseña cómo agregar código C++ a un paquete.

-   La @sec-rcpp-more concluye el capítulo con una lista de más recursos para ayudarlo a aprender Rcpp y C++.

### Requisitos previos {.unnumbered}

Usaremos [Rcpp](http://www.rcpp.org/) para llamar a C++ desde R:

```{r setup}
library(Rcpp)
```

También necesitará un compilador de C++ que funcione. Para conseguirlo:

-   En Windows, instale [Rtools](http://cran.r-project.org/bin/windows/Rtools/).
-   En Mac, instala Xcode desde la tienda de aplicaciones.
-   En Linux, `sudo apt-get install r-base-dev` o similar.

## Empezar con C++ {#sec-rcpp-intro}

`cppFunction()` le permite escribir funciones de C++ en R: \index{cppFunction()}

```{r add}
cppFunction('int add(int x, int y, int z) {
  int sum = x + y + z;
  return sum;
}')
# add funciona como una función R normal
add
add(1, 2, 3)
```

Cuando ejecute este código, Rcpp compilará el código C++ y construirá una función R que se conecta a la función C++ compilada. Suceden muchas cosas debajo del capó, pero Rcpp se encarga de todos los detalles para que no tengas que preocuparte por ellos.

Las siguientes secciones le enseñarán los conceptos básicos mediante la traducción de funciones simples de R a sus equivalentes de C++. Comenzaremos de manera simple con una función que no tiene entradas y una salida escalar, y luego la complicaremos progresivamente:

-   Entrada escalar y salida escalar
-   Entrada vectorial y salida escalar
-   Entrada vectorial y salida vectorial
-   Entrada matricial y salida vectorial

### Sin entradas, salida escalar

Comencemos con una función muy simple. No tiene argumentos y siempre devuelve el entero 1:

```{r one-r}
one <- function() 1L
```

The equivalent C++ function is:

``` cpp
int one() {
  return 1;
}
```

La podemos compilar y usar desde R con `cppFunction()`

```{r one-cpp}
cppFunction('int one() {
  return 1;
}')
```

Esta pequeña función ilustra una serie de diferencias importantes entre R y C++:

-   La sintaxis para crear una función se parece a la sintaxis para llamar a una función; no usa la asignación para crear funciones como lo hace en R.

-   Debe declarar el tipo de salida que devuelve la función. Esta función devuelve un `int` (un entero escalar). Las clases para los tipos más comunes de vectores R son: `NumericVector`, `IntegerVector`, `CharacterVector` y `LogicalVector`.

-   Los escalares y los vectores son diferentes. Los equivalentes escalares de vectores numéricos, enteros, de caracteres y lógicos son: `doble`, `int`, `String` y `bool`.

-   Debe usar una declaración `return` explícita para devolver un valor de una función.

-   Cada instrucción termina con un `;`.

### Entrada escalar, salida escalar

La siguiente función de ejemplo implementa una versión escalar de la función `sign()` que devuelve 1 si la entrada es positiva y -1 si es negativa:

```{r sign}
signR <- function(x) {
  if (x > 0) {
    1
  } else if (x == 0) {
    0
  } else {
    -1
  }
}

cppFunction('int signC(int x) {
  if (x > 0) {
    return 1;
  } else if (x == 0) {
    return 0;
  } else {
    return -1;
  }
}')
```

En la versión C++:

-   Declaramos el tipo de cada entrada de la misma forma que declaramos el tipo de la salida. Si bien esto hace que el código sea un poco más detallado, también aclara el tipo de entrada que necesita la función.

-   La sintaxis `if` es idéntica --- si bien existen grandes diferencias entre R y C++, ¡también hay muchas similitudes! C++ también tiene una instrucción `while` que funciona de la misma manera que la de R. Como en R, puede usar `break` para salir del ciclo, pero para omitir una iteración necesita usar `continue` en lugar de `next`.

### Entrada vectorial, salida escalar

Una gran diferencia entre R y C++ es que el costo de los bucles es mucho menor en C++. Por ejemplo, podríamos implementar la función `sum` en R usando un bucle. Si ha estado programando en R por un tiempo, ¡probablemente tendrá una reacción visceral a esta función!

```{r sum-r}
sumR <- function(x) {
  total <- 0
  for (i in seq_along(x)) {
    total <- total + x[i]
  }
  total
}
```

En C++, los bucles tienen muy poca sobrecarga, por lo que está bien usarlos. En la @sec-stl, verá alternativas a los bucles `for` que expresan más claramente su intención; no son más rápidos, pero pueden hacer que su código sea más fácil de entender.

```{r sum-cpp}
cppFunction('double sumC(NumericVector x) {
  int n = x.size();
  double total = 0;
  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return total;
}')
```

La versión C++ es similar, pero:

-   Para encontrar la longitud del vector, usamos el método `.size()`, que devuelve un número entero. Los métodos de C++ se llaman con `.` (es decir, un punto).

-   La sentencia `for` tiene una sintaxis diferente: `for(init; check; increment)`. Este bucle se inicializa creando una nueva variable llamada `i` con valor 0. Antes de cada iteración comprobamos que `i < n`, y terminamos el bucle si no es así. Después de cada iteración, incrementamos el valor de `i` en uno, utilizando el operador de prefijo especial `++` que aumenta el valor de `i` en 1.

-   En C++, los índices vectoriales comienzan en 0, lo que significa que el último elemento está en la posición `n - 1`. Repetiré esto porque es muy importante: **¡EN C++, LOS ÍNDICES VECTORIALES EMPIEZAN EN 0**! Esta es una fuente muy común de errores al convertir funciones de R a C++.

-   Utilice `=` para la asignación, no `<-`.

-   C++ proporciona operadores que modifican en el lugar: `total += x[i]` es equivalente a `total = total + x[i]`. Operadores in situ similares son `-=`, `*=` y `/=`.

Este es un buen ejemplo de dónde C++ es mucho más eficiente que R. Como se muestra en el siguiente micropunto de referencia, `sumC()` es competitivo con el integrado (y altamente optimizado) `sum()`, mientras que `sumR()` es varios órdenes de magnitud más lento.

```{r sum-bench}
x <- runif(1e3)
bench::mark(
  sum(x),
  sumC(x),
  sumR(x)
)[1:6]
```

### Entrada vectorial, salida vectorial

<!-- FIXME: Proponer un mejor ejemplo. También arreglar en otros dos lugares se produce -->

A continuación, crearemos una función que calcule la distancia euclidiana entre un valor y un vector de valores:

```{r pdist-r}
pdistR <- function(x, ys) {
  sqrt((x - ys) ^ 2)
}
```

En R, no es obvio que queremos que `x` sea un escalar de la definición de la función, y debemos dejarlo claro en la documentación. Eso no es un problema en la versión de C++ porque tenemos que ser explícitos con los tipos:

```{r pdist-cpp}
cppFunction('NumericVector pdistC(double x, NumericVector ys) {
  int n = ys.size();
  NumericVector out(n);

  for(int i = 0; i < n; ++i) {
    out[i] = sqrt(pow(ys[i] - x, 2.0));
  }
  return out;
}')
```

Esta función introduce solo algunos conceptos nuevos:

-   Creamos un nuevo vector numérico de longitud `n` con un constructor: `NumericVector out(n)`. Otra forma útil de hacer un vector es copiar uno existente: `NumericVector zs = clone(ys)`.

-   C++ usa `pow()`, no `^`, para la exponenciación.

Tenga en cuenta que debido a que la versión R está completamente vectorizada, ya será rápida.

```{r}
y <- runif(1e6)
bench::mark(
  pdistR(0.5, y),
  pdistC(0.5, y)
)[1:6]
```

En mi computadora, toma alrededor de 5 ms con un vector `y` de 1 millón de elementos. La función de C++ es unas 2,5 veces más rápida, \~2 ms, pero suponiendo que le haya llevado 10 minutos escribir la función de C++, necesitará ejecutarla \~200.000 veces para que valga la pena reescribirla. La razón por la que la función de C++ es más rápida es sutil y se relaciona con la administración de la memoria. La versión R necesita crear un vector intermedio de la misma longitud que y (`x - ys`), y asignar memoria es una operación costosa. La función de C++ evita esta sobrecarga porque usa un escalar intermedio.

```{r, include = FALSE}
# 5e-3 * x == 2e-3 * x + 10 * 60
600 / (5e-3 - 2e-3) 
```

### Usar sourceCpp {#sec-sourceCpp}

Hasta ahora, hemos usado C++ en línea con `cppFunction()`. Esto hace que la presentación sea más simple, pero para problemas reales, generalmente es más fácil usar archivos independientes de C++ y luego enviarlos a R usando `sourceCpp()`. Esto le permite aprovechar la compatibilidad con el editor de texto para archivos C++ (p. ej., resaltado de sintaxis), además de facilitar la identificación de los números de línea en los errores de compilación. \index{sourceCpp()}

Su archivo independiente de C++ debe tener la extensión `.cpp` y debe comenzar con:

``` cpp
#include <Rcpp.h>
using namespace Rcpp;
```

Y para cada función que desee que esté disponible dentro de R, debe prefijarla con:

``` cpp
// [[Rcpp::export]]
```

::: sidebar
Si está familiarizado con roxygen2, puede preguntarse cómo se relaciona esto con `@export`. `Rcpp::export` controla si una función se exporta de C++ a R; `@export` controla si una función se exporta desde un paquete y se pone a disposición del usuario.
:::

Puede incrustar código R en bloques de comentarios especiales de C++. Esto es realmente conveniente si desea ejecutar algún código de prueba:

``` cpp
/*** R
# This is R code
*/
```

El código R se ejecuta con `source(echo = TRUE)`, por lo que no es necesario que imprima la salida explícitamente.

Para compilar el código C++, usa `sourceCpp("ruta/al/archivo.cpp")`. Esto creará las funciones R coincidentes y las agregará a su sesión actual. Tenga en cuenta que estas funciones no se pueden guardar en un archivo `.Rdata` y volver a cargar en una sesión posterior; deben volver a crearse cada vez que reinicie R.

Por ejemplo, ejecutar `sourceCpp()` en el siguiente archivo implementa mean en C++ y luego lo compara con `mean()` integrado:

```{r, engine = "Rcpp", eval = FALSE}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double meanC(NumericVector x) {
  int n = x.size();
  double total = 0;

  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return total / n;
}

/*** R
x <- runif(1e5)
bench::mark(
  mean(x),
  meanC(x)
)
*/
```

NB: Si ejecuta este código, notará que `meanC()` es mucho más rápido que el `mean()` integrado. Esto se debe a que cambia la precisión numérica por la velocidad.

En el resto de este capítulo, el código C++ se presentará de forma independiente en lugar de envuelto en una llamada a `cppFunction`. Si desea intentar compilar y/o modificar los ejemplos, debe pegarlos en un archivo fuente de C++ que incluya los elementos descritos anteriormente. Esto es fácil de hacer en RMarkdown: todo lo que necesita hacer es especificar `engine = "Rcpp"`.

### Ejercicios {#sec-exercise-started}

1.  Con los conceptos básicos de C++ en la mano, ahora es un buen momento para practicar leyendo y escribiendo algunas funciones simples de C++. Para cada una de las siguientes funciones, lea el código y averigüe cuál es la función base R correspondiente. Es posible que aún no comprenda cada parte del código, pero debería poder descubrir los conceptos básicos de lo que hace la función.

    ``` cpp
    double f1(NumericVector x) {
      int n = x.size();
      double y = 0;

      for(int i = 0; i < n; ++i) {
        y += x[i] / n;
      }
      return y;
    }

    NumericVector f2(NumericVector x) {
      int n = x.size();
      NumericVector out(n);

      out[0] = x[0];
      for(int i = 1; i < n; ++i) {
        out[i] = out[i - 1] + x[i];
      }
      return out;
    }

    bool f3(LogicalVector x) {
      int n = x.size();

      for(int i = 0; i < n; ++i) {
        if (x[i]) return true;
      }
      return false;
    }

    int f4(Function pred, List x) {
      int n = x.size();

      for(int i = 0; i < n; ++i) {
        LogicalVector res = pred(x[i]);
        if (res[0]) return i + 1;
      }
      return 0;
    }

    NumericVector f5(NumericVector x, NumericVector y) {
      int n = std::max(x.size(), y.size());
      NumericVector x1 = rep_len(x, n);
      NumericVector y1 = rep_len(y, n);

      NumericVector out(n);

      for (int i = 0; i < n; ++i) {
        out[i] = std::min(x1[i], y1[i]);
      }

      return out;
    }
    ```

2.  Para practicar sus habilidades de escritura de funciones, convierta las siguientes funciones a C++. Por ahora, suponga que las entradas no tienen valores faltantes.

    1.  `all()`.

    2.  `cumprod()`, `cummin()`, `cummax()`.

    3.  `diff()`. Start by assuming lag 1, and then generalise for lag `n`.

    4.  `range()`.

    5.  `var()`. Lea acerca de los enfoques que puede adoptar en [Wikipedia](http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance). Siempre que se implemente un algoritmo numérico, siempre es bueno verificar lo que ya se sabe sobre el problema.

## Otras clases {#sec-rcpp-classes}

Ya has visto las clases vectoriales básicas (`IntegerVector`, `NumericVector`, `LogicalVector`, `CharacterVector`) y sus equivalentes escalares (`int`, `double`, `bool`, `String`). Rcpp también proporciona contenedores para todos los demás tipos de datos básicos. Los más importantes son para listas y marcos de datos, funciones y atributos, como se describe a continuación. Rcpp también proporciona clases para más tipos como `Environment`, `DottedPair`, `Language`, `Symbol`, etc., pero estos están más allá del alcance de este capítulo.

### Listas y data frames

Rcpp también proporciona las clases `List` y `DataFrame`, pero son más útiles para la salida que para la entrada. Esto se debe a que las listas y los marcos de datos pueden contener clases arbitrarias, pero C++ necesita conocer sus clases de antemano. Si la lista tiene una estructura conocida (p. ej., es un objeto de S3), puede extraer los componentes y convertirlos manualmente a sus equivalentes de C++ con `as()`. Por ejemplo, el objeto creado por `lm()`, la función que ajusta un modelo lineal, es una lista cuyos componentes son siempre del mismo tipo. El siguiente código ilustra cómo puede extraer el error porcentual medio (`mpe()`) de un modelo lineal. Este no es un buen ejemplo de cuándo usar C++, porque se implementa muy fácilmente en R, pero muestra cómo trabajar con una clase S3 importante. Tenga en cuenta el uso de `.inherits()` y `stop()` para verificar que el objeto realmente es un modelo lineal. \index{lists!in C++} \index{data frames!in C++}

<!-- FIXME: necesita una mejor motivación -->

```{r, engine = "Rcpp"}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double mpe(List mod) {
  if (!mod.inherits("lm")) stop("Input must be a linear model");

  NumericVector resid = as<NumericVector>(mod["residuals"]);
  NumericVector fitted = as<NumericVector>(mod["fitted.values"]);

  int n = resid.size();
  double err = 0;
  for(int i = 0; i < n; ++i) {
    err += resid[i] / (fitted[i] + resid[i]);
  }
  return err / n;
}
```

```{r}
mod <- lm(mpg ~ wt, data = mtcars)
mpe(mod)
```

### Funciones {#sec-functions-rcpp}

\index{functions!in C++}

Puede poner funciones R en un objeto de tipo `Function`. Esto hace que llamar a una función R desde C++ sea sencillo. El único desafío es que no sabemos qué tipo de salida devolverá la función, por lo que usamos el tipo general `RObject`.

```{r, engine = "Rcpp"}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
RObject callWithOne(Function f) {
  return f(1);
}
```

```{r}
callWithOne(function(x) x + 1)
callWithOne(paste)
```

Llamar funciones R con argumentos posicionales es obvio:

``` cpp
f("y", 1);
```

Pero necesita una sintaxis especial para argumentos con nombre:

``` cpp
f(_["x"] = "y", _["value"] = 1);
```

### Atributos

\index{attributes!in C++}

Todos los objetos R tienen atributos, que se pueden consultar y modificar con `.attr()`. Rcpp también proporciona `.names()` como un alias para el atributo de nombre. El siguiente fragmento de código ilustra estos métodos. Tenga en cuenta el uso de `::create()`, un método de *clase*. Esto le permite crear un vector R a partir de valores escalares de C++:

```{r attribs, engine = "Rcpp"}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector attribs() {
  NumericVector out = NumericVector::create(1, 2, 3);

  out.names() = CharacterVector::create("a", "b", "c");
  out.attr("my-attr") = "my-value";
  out.attr("class") = "my-class";

  return out;
}
```

Para los objetos S4, `.slot()` juega un papel similar a `.attr()`.

## Valores faltantes {#sec-rcpp-na}

\index{NA}

Si está trabajando con valores perdidos, necesita saber dos cosas:

-   Cómo se comportan los valores faltantes de R en los escalares de C++ (por ejemplo, `doble`).
-   Cómo obtener y establecer valores faltantes en vectores (por ejemplo, `NumericVector`).

### Escalares

El siguiente código explora lo que sucede cuando tomas uno de los valores faltantes de R, lo conviertes en un escalar y luego vuelves a convertirlo en un vector R. Tenga en cuenta que este tipo de experimentación es una forma útil de descubrir qué hace cualquier operación.

```{r missings, engine = "Rcpp"}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List scalar_missings() {
  int int_s = NA_INTEGER;
  String chr_s = NA_STRING;
  bool lgl_s = NA_LOGICAL;
  double num_s = NA_REAL;

  return List::create(int_s, chr_s, lgl_s, num_s);
}
```

```{r}
str(scalar_missings())
```

Con la excepción de `bool`, las cosas se ven bastante bien aquí: todos los valores que faltan se han conservado. Sin embargo, como veremos en las siguientes secciones, las cosas no son tan sencillas como parecen.

#### Enteros

Con números enteros, los valores faltantes se almacenan como el número entero más pequeño. Si no les haces nada, se conservarán. Pero, dado que C++ no sabe que el entero más pequeño tiene este comportamiento especial, si hace algo al respecto, es probable que obtenga un valor incorrecto: por ejemplo, `evalCpp('NA_INTEGER + 1')` da -2147483647.

Entonces, si desea trabajar con valores faltantes en números enteros, use un `IntegerVector` de longitud 1 o tenga mucho cuidado con su código.

#### Dobles

Con los dobles, es posible que pueda ignorar los valores faltantes y trabajar con NaN (no un número). Esto se debe a que el NA de R es un tipo especial de número NaN de punto flotante IEEE 754. Entonces, cualquier expresión lógica que involucre un NaN (o en C++, NAN) siempre se evalúa como FALSE:

```{r, echo = FALSE, message = FALSE}
library(Rcpp)
```

```{r}
evalCpp("NAN == 1")
evalCpp("NAN < 1")
evalCpp("NAN > 1")
evalCpp("NAN == NAN")
```

(Aquí estoy usando `evalCpp()` que le permite ver el resultado de ejecutar una sola expresión de C++, lo que la hace excelente para este tipo de experimentación interactiva.)

Pero tenga cuidado al combinarlos con valores booleanos:

```{r}
evalCpp("NAN && TRUE")
evalCpp("NAN || FALSE")
```

Sin embargo, en contextos numéricos, los NaN propagarán los NA:

```{r}
evalCpp("NAN + 1")
evalCpp("NAN - 1")
evalCpp("NAN / 1")
evalCpp("NAN * 1")
```

### Caracteres

`String` es una clase de cadena de caracteres escalar introducida por Rcpp, por lo que sabe cómo lidiar con los valores faltantes.

### Booleano

Mientras que `bool` de C++ tiene dos valores posibles (`true` o `false`), un vector lógico en R tiene tres (`TRUE`, `FALSE` y `NA`). Si fuerza un vector lógico de longitud 1, asegúrese de que no contenga ningún valor faltante; de lo contrario, se convertirán en VERDADERO. Una solución fácil es usar `int` en su lugar, ya que esto puede representar `TRUE`, `FALSE` y `NA`.

### Vectores {#sec-vectors-rcpp}

Con los vectores, debe usar un valor perdido específico para el tipo de vector, `NA_REAL`, `NA_INTEGER`, `NA_LOGICAL`, `NA_STRING`:

```{r, engine = "Rcpp"}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List missing_sampler() {
  return List::create(
    NumericVector::create(NA_REAL),
    IntegerVector::create(NA_INTEGER),
    LogicalVector::create(NA_LOGICAL),
    CharacterVector::create(NA_STRING)
  );
}
```

```{r}
str(missing_sampler())
```

### Ejercicios

1.  Vuelva a escribir cualquiera de las funciones del primer ejercicio de la @sec-exercise-started para tratar con los valores faltantes. Si `na.rm` es verdadero, ignore los valores faltantes. Si `na.rm` es falso, devuelve un valor faltante si la entrada contiene valores faltantes. Algunas buenas funciones para practicar son `min()`, `max()`, `range()`, `mean()` y `var()`.

2.  Vuelva a escribir `cumsum()` y `diff()` para que puedan manejar los valores faltantes. Tenga en cuenta que estas funciones tienen un comportamiento un poco más complicado.

## Biblioteca de plantillas estándar {#sec-stl}

La verdadera fortaleza de C++ se revela cuando necesita implementar algoritmos más complejos. La biblioteca de plantillas estándar (STL) proporciona un conjunto de estructuras de datos y algoritmos extremadamente útiles. Esta sección explicará algunos de los algoritmos y estructuras de datos más importantes y le indicará la dirección correcta para obtener más información. No puedo enseñarte todo lo que necesitas saber sobre STL, pero espero que los ejemplos te muestren el poder de STL y te convenzan de que es útil aprender más. \index{standard template library}

Si necesita un algoritmo o una estructura de datos que no esté implementado en STL, un buen lugar para buscar es [boost](http://www.boost.org/doc/). La instalación de boost en su computadora está más allá del alcance de este capítulo, pero una vez que lo haya instalado, puede usar las estructuras de datos y los algoritmos de boost al incluir el archivo de encabezado apropiado con (p. ej.) `#include <boost/array.hpp>`.

### Usar iteradores

Los iteradores se utilizan mucho en STL: muchas funciones aceptan o devuelven iteradores. Son el siguiente paso de los bucles básicos, abstrayendo los detalles de la estructura de datos subyacente. Los iteradores tienen tres operadores principales: \index{iterators}

1.  Avanza con `++`.
2.  Obtener el valor al que se refieren, o **desreferenciar**, con `*`.
3.  Comparar con `==`.

Por ejemplo, podríamos reescribir nuestra función de suma usando iteradores:

```{r, engine = "Rcpp"}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double sum3(NumericVector x) {
  double total = 0;
  
  NumericVector::iterator it;
  for(it = x.begin(); it != x.end(); ++it) {
    total += *it;
  }
  return total;
}
```

Los principales cambios están en el bucle for:

-   Comenzamos en `x.begin()` y repetimos hasta llegar a `x.end()`. Una pequeña optimización es almacenar el valor del iterador final para que no tengamos que buscarlo cada vez. Esto solo ahorra alrededor de 2 ns por iteración, por lo que solo es importante cuando los cálculos en el ciclo son muy simples.

-   En lugar de indexar en x, usamos el operador de desreferencia para obtener su valor actual: `*it`.

-   Observe el tipo de iterador: `NumericVector::iterator`. Cada tipo de vector tiene su propio tipo de iterador: `LogicalVector::iterator`, `CharacterVector::iterator`, etc.

Este código se puede simplificar aún más mediante el uso de una característica de C++11: bucles for basados en rangos. C ++ 11 está ampliamente disponible y se puede activar fácilmente para usar con Rcpp agregando `[[Rcpp::plugins(cpp11)]]`.

```{r, engine = "Rcpp"}
// [[Rcpp::plugins(cpp11)]]
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double sum4(NumericVector xs) {
  double total = 0;
  
  for(const auto &x : xs) {
    total += x;
  }
  return total;
}
```

Los iteradores también nos permiten usar los equivalentes de C++ de la familia de funciones apply. Por ejemplo, podríamos volver a escribir `sum()` para usar la función `accumulate()`, que toma un iterador inicial y final, y suma todos los valores en el vector. El tercer argumento para `acumular` da el valor inicial: es particularmente importante porque esto también determina el tipo de datos que usa `acumular` (así que usamos `0.0` y no `0` para que `acumule` use un `doble`, no un `int`.). Para usar `accumulate()` necesitamos incluir el encabezado `<numeric>`.

```{r, engine = "Rcpp"}
#include <numeric>
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double sum5(NumericVector x) {
  return std::accumulate(x.begin(), x.end(), 0.0);
}
```

### Algoritmos

El encabezado `<algorithm>` proporciona una gran cantidad de algoritmos que funcionan con iteradores. Una buena referencia está disponible en <https://en.cppreference.com/w/cpp/algorithm>. Por ejemplo, podríamos escribir una versión básica de Rcpp de `findInterval()` que toma dos argumentos, un vector de valores y un vector de rupturas, y ubica el contenedor en el que cae cada x. Esto muestra algunas características más avanzadas del iterador. Lea el código a continuación y vea si puede descubrir cómo funciona. \index{findInterval()}

```{r, engine = "Rcpp"}
#include <algorithm>
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerVector findInterval2(NumericVector x, NumericVector breaks) {
  IntegerVector out(x.size());

  NumericVector::iterator it, pos;
  IntegerVector::iterator out_it;

  for(it = x.begin(), out_it = out.begin(); it != x.end(); 
      ++it, ++out_it) {
    pos = std::upper_bound(breaks.begin(), breaks.end(), *it);
    *out_it = std::distance(breaks.begin(), pos);
  }

  return out;
}
```

Los puntos clave son:

-   Pasamos por dos iteradores (entrada y salida) simultáneamente.

-   Podemos asignar a un iterador desreferenciado (`out_it`) para cambiar los valores en `out`.

-   `upper_bound()` devuelve un iterador. Si quisiéramos el valor de `upper_bound()`, podríamos quitarle la referencia; para averiguar su ubicación, usamos la función `distance()`.

-   Pequeña nota: si queremos que esta función sea tan rápida como `findInterval()` en R (que usa código C escrito a mano), necesitamos calcular las llamadas a `.begin()` y `.end()` una vez y guardar los resultados. Esto es fácil, pero distrae la atención de este ejemplo, por lo que se ha omitido. Hacer este cambio produce una función que es un poco más rápida que la función `findInterval()` de R, pero es aproximadamente 1/10 del código.

En general, es mejor usar algoritmos de STL que bucles enrollados a mano. En *Effective STL*, Scott Meyers da tres razones: eficiencia, corrección y mantenibilidad. Los algoritmos de STL están escritos por expertos en C++ para que sean extremadamente eficientes y han existido durante mucho tiempo, por lo que están bien probados. El uso de algoritmos estándar también hace que la intención de su código sea más clara, lo que ayuda a que sea más legible y fácil de mantener.

### Estructuras de datos {#sec-data-structures-rcpp}

El STL proporciona un gran conjunto de estructuras de datos: `array`, `bitset`, `list`, `forward_list`, `map`, `multimap`, `multiset`, `priority_queue`, `queue`, `deque`, `set`, `stack`, `unordered_map`, `unordered_set`, `unordered_multimap`, `unordered_multiset`, y `vector`. Las más importantes de estas estructuras de datos son el `vector`, el `unordered_set` y el `unordered_map`. Nos centraremos en estos tres en esta sección, pero el uso de los otros es similar: solo tienen diferentes compensaciones de rendimiento. Por ejemplo, `deque` (pronunciado "deck") tiene una interfaz muy similar a los vectores pero una implementación subyacente diferente que tiene diferentes compensaciones de rendimiento. Es posible que desee probarlo para su problema. Una buena referencia para las estructuras de datos STL es <https://en.cppreference.com/w/cpp/container> --- Le recomiendo que lo mantenga abierto mientras trabaja con STL.

Rcpp sabe cómo convertir muchas estructuras de datos STL a sus equivalentes R, por lo que puede devolverlas desde sus funciones sin convertirlas explícitamente en estructuras de datos R.

### Vectores {#sec-vectors-stl}

Un vector STL es muy similar a un vector R, excepto que crece de manera eficiente. Esto hace que los vectores sean apropiados para usar cuando no se sabe de antemano qué tan grande será la salida. Los vectores tienen una plantilla, lo que significa que debe especificar el tipo de objeto que contendrá el vector cuando lo cree: `vector<int>`, `vector<bool>`, `vector<double>`, `vector<String>`. Puede acceder a elementos individuales de un vector usando la notación estándar `[]`, y puede agregar un nuevo elemento al final del vector usando `.push_back()`. Si tiene una idea de antemano de qué tan grande será el vector, puede usar `.reserve()` para asignar suficiente almacenamiento. \index{vectors!in C++}

El siguiente código implementa la codificación de longitud de ejecución (`rle()`). Produce dos vectores de salida: un vector de valores y un vector de "longitudes" que indica cuántas veces se repite cada elemento. Funciona recorriendo el vector de entrada `x` comparando cada valor con el anterior: si es el mismo, incrementa el último valor en `lengths`; si es diferente, agrega el valor al final de `values` y establece la longitud correspondiente en 1.

```{r, engine = "Rcpp"}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List rleC(NumericVector x) {
  std::vector<int> lengths;
  std::vector<double> values;

  // Initialise first value
  int i = 0;
  double prev = x[0];
  values.push_back(prev);
  lengths.push_back(1);

  NumericVector::iterator it;
  for(it = x.begin() + 1; it != x.end(); ++it) {
    if (prev == *it) {
      lengths[i]++;
    } else {
      values.push_back(*it);
      lengths.push_back(1);

      i++;
      prev = *it;
    }
  }

  return List::create(
    _["lengths"] = lengths, 
    _["values"] = values
  );
}
```

(Una implementación alternativa sería reemplazar `i` con el iterador `lengths.rbegin()` que siempre apunta al último elemento del vector. Es posible que desee intentar implementar eso.)

Otros métodos de un vector se describen en <https://en.cppreference.com/w/cpp/container/vector>.

### Conjuntos

Los conjuntos mantienen un conjunto único de valores y pueden indicar de manera eficiente si ha visto un valor antes. Son útiles para problemas que involucran duplicados o valores únicos (como `unique`, `duplicated` o `in`). C++ proporciona tanto conjuntos ordenados (`std::set`) como desordenados (`std::unordered_set`), dependiendo de si el orden es importante para usted o no. Los conjuntos desordenados tienden a ser mucho más rápidos (porque usan una tabla hash internamente en lugar de un árbol), por lo que incluso si necesita un conjunto ordenado, debe considerar usar un conjunto desordenado y luego ordenar la salida. Al igual que los vectores, los conjuntos tienen plantillas, por lo que debe solicitar el tipo de conjunto adecuado para su propósito: `unordered_set<int>`, `unordered_set<bool>`, etc. Hay más detalles disponibles en <https://en.cppreference.com/w/cpp/container/set> y <https://en.cppreference.com/w/cpp/container/unordered_set>. \index{sets}

La siguiente función usa un conjunto desordenado para implementar un equivalente a `duplicated()` para vectores enteros. Tenga en cuenta el uso de `seen.insert(x[i]).second`. `insert()` devuelve un par, el valor `.first` es un iterador que apunta al elemento y el valor `.second` es un valor booleano que es verdadero si el valor fue una nueva adición al conjunto.

```{r, engine = "Rcpp"}
// [[Rcpp::plugins(cpp11)]]
#include <Rcpp.h>
#include <unordered_set>
using namespace Rcpp;

// [[Rcpp::export]]
LogicalVector duplicatedC(IntegerVector x) {
  std::unordered_set<int> seen;
  int n = x.size();
  LogicalVector out(n);

  for (int i = 0; i < n; ++i) {
    out[i] = !seen.insert(x[i]).second;
  }

  return out;
}
```

### Map

\index{hashmaps}

Un mapa es similar a un conjunto, pero en lugar de almacenar presencia o ausencia, puede almacenar datos adicionales. Es útil para funciones como `table()` o `match()` que necesitan buscar un valor. Al igual que con los conjuntos, existen versiones ordenadas (`std::map`) y desordenadas (`std::unordered_map`). Dado que los mapas tienen un valor y una clave, debe especificar ambos tipos al inicializar un mapa: `mapa<doble, int>`, `mapa_desordenado<int, doble>`, etc. El siguiente ejemplo muestra cómo podrías usar un `map` para implementar `table()` para vectores numéricos:

```{r, engine = "Rcpp"}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
std::map<double, int> tableC(NumericVector x) {
  std::map<double, int> counts;

  int n = x.size();
  for (int i = 0; i < n; i++) {
    counts[x[i]]++;
  }

  return counts;
}
```

### Ejercicios

Para practicar el uso de estructuras de datos y algoritmos STL, implemente lo siguiente mediante funciones R en C++, utilizando las sugerencias proporcionadas:

1.  `median.default()` usando `partial_sort`.

2.  `%in%` usando `unordered_set` y los [`find()`](https://en.cppreference.com/w/cpp/container/unordered_set/find) o [`count()`](https://en.cppreference.com/w/cpp/container/unordered_set/count) métodos.

3.  `unique()` usando un `unordered_set` (desafío: ¡hazlo en una línea!).

4.  `min()` usando `std::min()`, o `max()` usando `std::max()`.

5.  `which.min()` usando `min_element`, o `which.max()` usando `max_element`.

6.  `setdiff()`, `union()`, y `intersect()` para números enteros usando rangos ordenados y `set_union`, `set_intersection` y `set_difference`.

## Caso de estudio {#sec-rcpp-case-studies}

Los siguientes estudios de casos ilustran algunos usos reales de C++ para reemplazar el código R lento.

### Muestreador de Gibbs

<!-- FIXME: necesita más contexto? -->

El siguiente estudio de caso actualiza un ejemplo [sobre el que se escribió en un blog](http://dirk.eddelbuettel.com/blog/2011/07/14/) de Dirk Eddelbuettel, que ilustra la conversión de una muestra de Gibbs en R a C++. El código R y C++ que se muestra a continuación es muy similar (solo tomó unos minutos convertir la versión R a la versión C++), pero se ejecuta unas 20 veces más rápido en mi computadora. La publicación del blog de Dirk también muestra otra forma de hacerlo aún más rápido: usar las funciones de generador de números aleatorios más rápidas en GSL (fácilmente accesible desde R a través del paquete RcppGSL) puede hacerlo otras dos o tres veces más rápido. \index{Gibbs sampler}

El código R es el siguiente:

```{r}
gibbs_r <- function(N, thin) {
  mat <- matrix(nrow = N, ncol = 2)
  x <- y <- 0

  for (i in 1:N) {
    for (j in 1:thin) {
      x <- rgamma(1, 3, y * y + 4)
      y <- rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))
    }
    mat[i, ] <- c(x, y)
  }
  mat
}
```

Esto es fácil de convertir a C++. Nosotros:

-   Agregar declaraciones de tipo a todas las variables.

-   Utilice `(` en lugar de `[` para indexar en la matriz.

-   Subíndice los resultados de `rgamma` y `rnorm` para convertir de un vector a un escalar.

```{r, engine = "Rcpp"}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix gibbs_cpp(int N, int thin) {
  NumericMatrix mat(N, 2);
  double x = 0, y = 0;

  for(int i = 0; i < N; i++) {
    for(int j = 0; j < thin; j++) {
      x = rgamma(1, 3, 1 / (y * y + 4))[0];
      y = rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))[0];
    }
    mat(i, 0) = x;
    mat(i, 1) = y;
  }

  return(mat);
}
```

La evaluación comparativa de los dos rendimientos de implementación:

```{r}
bench::mark(
  gibbs_r(100, 10),
  gibbs_cpp(100, 10),
  check = FALSE
)
```

### R vectorización frente a vectorización C++

<!-- FIXME: necesita más contexto? -->

Este ejemplo está adaptado de ["Rcpp fuma rápido para modelos basados en agentes en marcos de datos"](https://gweissman.github.io/post/rcpp-is-Smoking-fast-for-agent-based-models-%20en-marcos-de-datos/). El desafío es predecir la respuesta de un modelo a partir de tres entradas. La versión básica de R del predictor se ve así:

```{r}
vacc1a <- function(age, female, ily) {
  p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
  p <- p * if (female) 1.25 else 0.75
  p <- max(0, p)
  p <- min(1, p)
  p
}
```

Queremos poder aplicar esta función a muchas entradas, por lo que podríamos escribir una versión de entrada vectorial usando un bucle for.

```{r}
vacc1 <- function(age, female, ily) {
  n <- length(age)
  out <- numeric(n)
  for (i in seq_len(n)) {
    out[i] <- vacc1a(age[i], female[i], ily[i])
  }
  out
}
```

Si está familiarizado con R, tendrá el presentimiento de que esto será lento, y de hecho lo es. Hay dos formas en las que podemos atacar este problema. Si tiene un buen vocabulario de R, puede ver inmediatamente cómo vectorizar la función (usando `ifelse()`, `pmin()` y `pmax()`). Alternativamente, podríamos reescribir `vacc1a()` y `vacc1()` en C++, utilizando nuestro conocimiento de que los bucles y las llamadas a funciones tienen una sobrecarga mucho menor en C++.

Cualquiera de los dos enfoques es bastante sencillo. En R:

```{r}
vacc2 <- function(age, female, ily) {
  p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
  p <- p * ifelse(female, 1.25, 0.75)
  p <- pmax(0, p)
  p <- pmin(1, p)
  p
}
```

(Si ha trabajado mucho con R, es posible que reconozca algunos cuellos de botella potenciales en este código: se sabe que `ifelse`, `pmin` y `pmax` son lentos y podrían reemplazarse con `p * 0.75 + p * 0.5 * mujer`, `p[p < 0] <- 0`, `p[p > 1] <- 1`. Es posible que desee intentar cronometrar esas variaciones.)

O en C++:

```{r engine = "Rcpp"}
#include <Rcpp.h>
using namespace Rcpp;

double vacc3a(double age, bool female, bool ily){
  double p = 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily;
  p = p * (female ? 1.25 : 0.75);
  p = std::max(p, 0.0);
  p = std::min(p, 1.0);
  return p;
}

// [[Rcpp::export]]
NumericVector vacc3(NumericVector age, LogicalVector female, 
                    LogicalVector ily) {
  int n = age.size();
  NumericVector out(n);

  for(int i = 0; i < n; ++i) {
    out[i] = vacc3a(age[i], female[i], ily[i]);
  }

  return out;
}
```

A continuación, generamos algunos datos de muestra y verificamos que las tres versiones devuelvan los mismos valores:

```{r}
n <- 1000
age <- rnorm(n, mean = 50, sd = 10)
female <- sample(c(T, F), n, rep = TRUE)
ily <- sample(c(T, F), n, prob = c(0.8, 0.2), rep = TRUE)

stopifnot(
  all.equal(vacc1(age, female, ily), vacc2(age, female, ily)),
  all.equal(vacc1(age, female, ily), vacc3(age, female, ily))
)
```

La publicación del blog original olvidó hacer esto e introdujo un error en la versión C++: usaba `0.004` en lugar de `0.04`. Finalmente, podemos comparar nuestros tres enfoques:

```{r}
bench::mark(
  vacc1 = vacc1(age, female, ily),
  vacc2 = vacc2(age, female, ily),
  vacc3 = vacc3(age, female, ily)
)
```

No es sorprendente que nuestro enfoque original con bucles sea muy lento. La vectorización en R proporciona una gran aceleración, y podemos obtener aún más rendimiento (unas diez veces) con el ciclo de C++. Me sorprendió un poco que C++ fuera mucho más rápido, pero se debe a que la versión R tiene que crear 11 vectores para almacenar resultados intermedios, mientras que el código C++ solo necesita crear 1.

## Using Rcpp in a package {#sec-rcpp-package}

El mismo código C++ que se usa con `sourceCpp()` también se puede agrupar en un paquete. Hay varios beneficios de mover el código de un archivo fuente de C++ independiente a un paquete: \index{Rcpp!in a package}

1.  Su código puede estar disponible para los usuarios sin las herramientas de desarrollo de C++.

2.  El sistema de compilación de paquetes R gestiona automáticamente varios archivos de origen y sus dependencias.

3.  Los paquetes brindan infraestructura adicional para pruebas, documentación y consistencia.

Para agregar `Rcpp` a un paquete existente, coloque sus archivos C++ en el directorio `src/` y cree o modifique los siguientes archivos de configuración:

-   En `DESCRIPTION` añade

    ```         
    LinkingTo: Rcpp
    Imports: Rcpp
    ```

-   Asegúrese de que su `NAMESPACE` incluye:

    ```         
    useDynLib(mypackage)
    importFrom(Rcpp, sourceCpp)
    ```

    Necesitamos importar algo (cualquier cosa) de Rcpp para que el código interno de Rcpp se cargue correctamente. Este es un error en R y, con suerte, se solucionará en el futuro.

La forma más fácil de configurar esto automáticamente es llamar a `usethis::use_rcpp()`.

Antes de compilar el paquete, deberá ejecutar `Rcpp::compileAttributes()`. Esta función escanea los archivos C++ en busca de atributos `Rcpp::export` y genera el código necesario para que las funciones estén disponibles en R. Vuelva a ejecutar `compileAttributes()` cada vez que se agreguen, eliminen o cambien sus firmas. Esto lo hace automáticamente el paquete devtools y Rstudio.

Para obtener más detalles, consulte la viñeta del paquete Rcpp, `vignette("Rcpp-package")`.

## Aprendiendo más {#sec-rcpp-more}

Este capítulo solo ha tocado una pequeña parte de Rcpp, brindándole las herramientas básicas para reescribir código R de bajo rendimiento en C++. Como se señaló, Rcpp tiene muchas otras capacidades que facilitan la interfaz de R con el código C++ existente, que incluyen:

-   Características adicionales de los atributos, incluida la especificación de argumentos predeterminados, la vinculación en dependencias externas de C++ y la exportación de interfaces de C++ desde paquetes. Estas características y más están cubiertas en la viñeta de atributos de Rcpp, `vignette("Rcpp-attributes")`.

-   Creación automática de contenedores entre estructuras de datos de C++ y estructuras de datos de R, incluida la asignación de clases de C++ a clases de referencia. Una buena introducción a este tema es la viñeta de módulos Rcpp, `vignette("Rcpp-modules")`.

-   La guía de referencia rápida de Rcpp, `vignette("Rcpp-quickref")`, contiene un resumen útil de las clases de Rcpp y lenguajes de programación comunes.

Recomiendo encarecidamente estar atento a la [página de inicio de Rcpp](http://www.rcpp.org) y registrarse en la [lista de correo de Rcpp](http://lists.r-forge.r-project.org/%20cgi-bin/mailman/listinfo/rcpp-devel).

Otros recursos que he encontrado útiles para aprender C++ son:

-   *Effective C++* [@effective-cpp] y *Effective STL* [@effective-stl].

-   [*C++ Annotations*](http://www.icce.rug.nl/documents/cplusplus/cplusplus.html), dirigido a usuarios expertos en C (o cualquier otro lenguaje que use una gramática similar a C, como Perl o Java) que deseen saber más sobre C++ o hacer la transición a C++.

-   [*Algorithm Libraries*](http://www.cs.helsinki.fi/u/tpkarkka/alglib/k06/), que proporciona una descripción más técnica, pero aún concisa, de conceptos importantes de STL. (Siga los enlaces debajo de las notas).

Escribir código de rendimiento también puede requerir que reconsidere su enfoque básico: una sólida comprensión de las estructuras de datos y los algoritmos básicos es muy útil aquí. Eso está más allá del alcance de este libro, pero sugeriría el *Manual de diseño de algoritmos* [@alg-design-man], [*Introducción a los algoritmos* del MIT](http://ocw.mit.edu/courses/electrical%20-engineering-and-computer-science/6-046j-introduction-to-algorithms-sma-5503-fall-2005/), *Algorithms* de Robert Sedgewick y Kevin Wayne, que tiene un [libro de texto en línea gratuito](http:/%20/algs4.cs.princeton.edu/home/) y un \[curso de Coursera\] correspondiente (https://www.coursera.org/learn/algorithms-part1).

## Reconocimientos

Me gustaría agradecer a la lista de correo de Rcpp por muchas conversaciones útiles, en particular a Romain Francois y Dirk Eddelbuettel, quienes no solo han brindado respuestas detalladas a muchas de mis preguntas, sino que también han sido increíblemente receptivos para mejorar Rcpp. Este capítulo no hubiera sido posible sin JJ Allaire; me animó a aprender C++ y luego respondió muchas de mis preguntas tontas en el camino.