25  Reescribiendo código de R en C++

25.1 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 (Eddelbuettel y François 2011) (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.

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

• La Sección 25.2 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 Sección 25.2.5 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 Sección 25.3 analiza cómo modificar los atributos de Rcpp y menciona algunas de las otras clases importantes.

• La Sección 25.4 le enseña cómo trabajar con los valores faltantes de R en C++.

• La Sección 25.5 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 Sección 25.6 muestra dos estudios de casos reales en los que se utilizó Rcpp para obtener mejoras de rendimiento considerables.

• La Sección 25.7 le enseña cómo agregar código C++ a un paquete.

• La Sección 25.8 concluye el capítulo con una lista de más recursos para ayudarlo a aprender Rcpp y C++.

Requisitos previos

Usaremos Rcpp para llamar a C++ desde R:

library(Rcpp)

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

• En Windows, instale Rtools.
• En Mac, instala Xcode desde la tienda de aplicaciones.
• En Linux, sudo apt-get install r-base-dev o similar.

25.2 Empezar con C++

cppFunction() le permite escribir funciones de C++ en R:

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
#> function (x, y, z) 
#> .Call(<pointer: 0x7f777d8fd4d0>, x, y, z)
add(1, 2, 3)
#> [1] 6

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

25.2.1 Sin entradas, salida escalar

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

one <- function() 1L

The equivalent C++ function is:

int one() {
  return 1;
}

La podemos compilar y usar desde R con cppFunction()

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

25.2.2 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:

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.

25.2.3 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!

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 Sección 25.5, 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.

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.

x <- runif(1e3)
bench::mark(
  sum(x),
  sumC(x),
  sumR(x)
)[1:6]
#> # A tibble: 3 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 sum(x)       2.06µs    2.1µs   467352.        0B     0   
#> 2 sumC(x)      1.67µs   1.76µs   537590.        0B     0   
#> 3 sumR(x)     20.14µs  20.52µs    48130.      22KB     4.81

25.2.4 Entrada vectorial, salida vectorial

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

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:

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.

y <- runif(1e6)
bench::mark(
  pdistR(0.5, y),
  pdistC(0.5, y)
)[1:6]
#> # A tibble: 2 × 6
#>   expression          min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>     <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 pdistR(0.5, y)   4.97ms   5.05ms      197.    7.63MB     103.
#> 2 pdistC(0.5, y)   4.09ms   4.14ms      240.    7.63MB     116.

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.

25.2.5 Usar 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.

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

#include <Rcpp.h>
using namespace Rcpp;

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

// [[Rcpp::export]]

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:

/*** 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:

#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".

25.2.6 Ejercicios

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.

   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. Siempre que se implemente un algoritmo numérico, siempre es bueno verificar lo que ya se sabe sobre el problema.

25.3 Otras clases

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.

25.3.1 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.

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

mod <- lm(mpg ~ wt, data = mtcars)
mpe(mod)
#> [1] -0.0154

25.3.2 Funciones

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.

#include <Rcpp.h>
using namespace Rcpp;

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

callWithOne(function(x) x + 1)
#> [1] 2
callWithOne(paste)
#> [1] "1"

Llamar funciones R con argumentos posicionales es obvio:

f("y", 1);

Pero necesita una sintaxis especial para argumentos con nombre:

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

25.3.3 Atributos

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++:

#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().

25.4 Valores faltantes

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).

25.4.1 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.

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

str(scalar_missings())
#> List of 4
#>  $ : int NA
#>  $ : chr NA
#>  $ : logi TRUE
#>  $ : num NA

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.

25.4.1.1 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.

25.4.1.2 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:

evalCpp("NAN == 1")
#> [1] FALSE
evalCpp("NAN < 1")
#> [1] FALSE
evalCpp("NAN > 1")
#> [1] FALSE
evalCpp("NAN == NAN")
#> [1] FALSE

(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:

evalCpp("NAN && TRUE")
#> [1] TRUE
evalCpp("NAN || FALSE")
#> [1] TRUE

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

evalCpp("NAN + 1")
#> [1] NaN
evalCpp("NAN - 1")
#> [1] NaN
evalCpp("NAN / 1")
#> [1] NaN
evalCpp("NAN * 1")
#> [1] NaN

25.4.2 Caracteres

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

25.4.3 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.

25.4.4 Vectores

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

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

str(missing_sampler())
#> List of 4
#>  $ : num NA
#>  $ : int NA
#>  $ : logi NA
#>  $ : chr NA

25.4.5 Ejercicios

1. Vuelva a escribir cualquiera de las funciones del primer ejercicio de la Sección 25.2.6 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.

25.5 Biblioteca de plantillas estándar

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.

Si necesita un algoritmo o una estructura de datos que no esté implementado en STL, un buen lugar para buscar es boost. 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>.

25.5.1 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:

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:

#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)]].

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

#include <numeric>
#include <Rcpp.h>
using namespace Rcpp;

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

25.5.2 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.

#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.

25.5.3 Estructuras de datos

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.

25.5.4 Vectores

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.

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.

#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.

25.5.5 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.

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.

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

25.5.6 Map

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:

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

25.5.7 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() o 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.

25.6 Caso de estudio

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

25.6.1 Muestreador de Gibbs

El siguiente estudio de caso actualiza un ejemplo sobre el que se escribió en un blog 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.

El código R es el siguiente:

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.

#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:

bench::mark(
  gibbs_r(100, 10),
  gibbs_cpp(100, 10),
  check = FALSE
)
#> # A tibble: 2 × 6
#>   expression              min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>         <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 gibbs_r(100, 10)     2.42ms   2.51ms      387.   107.5KB     20.9
#> 2 gibbs_cpp(100, 10) 244.33µs 258.25µs     3802.    1.61KB     16.8

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

Este ejemplo está adaptado de “Rcpp fuma rápido para modelos basados en agentes en 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í:

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.

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:

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++:

#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:

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:

bench::mark(
  vacc1 = vacc1(age, female, ily),
  vacc2 = vacc2(age, female, ily),
  vacc3 = vacc3(age, female, ily)
)
#> # A tibble: 3 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 vacc1        1.51ms   1.55ms      638.    7.86KB    19.5 
#> 2 vacc2        42.7µs  45.08µs    21105.  148.56KB    42.0 
#> 3 vacc3       11.41µs  11.73µs    84445.   14.48KB     8.45

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.

25.7 Using Rcpp in a 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:

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").

25.8 Aprendiendo más

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 y registrarse en la lista de correo de Rcpp.

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

• Effective C++ (Meyers 2005) y Effective STL (Meyers 2001).

• C++ Annotations, 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, 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 (Skiena 1998), Introducción a los algoritmos del MIT, Algorithms de Robert Sedgewick y Kevin Wayne, que tiene un libro de texto en línea gratuito y un [curso de Coursera] correspondiente (https://www.coursera.org/learn/algorithms-part1).

25.9 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.

Eddelbuettel, Dirk, y Romain François. 2011. «Rcpp: Seamless R and C++ Integration». Journal of Statistical Software 40 (8): 1-18. https://doi.org/10.18637/jss.v040.i08.

Meyers, Scott. 2001. Effective STL: 50 specific ways to improve your use of the standard template library. Pearson Education. http://amzn.com/0201749629.

———. 2005. Effective C++: 55 specific ways to improve your programs and designs. Pearson Education. http://amzn.com/0321334876.

Skiena, Steven S. 1998. The Algorithm Design Manual. Springer. http://amzn.com/0387948600.