Chapitre 5 Rcpp ou comment intégrer facilement du code C++dans un package R

Rcpp (R-C-Plus-Plus) est un package qui facilite l’interface entre C++ et R. R est un langage interprété, ce qui facilite un certain nombre de choses (notamment nous donne accès à la console dans laquelle on peut évaluer du code à la volée). Néanmoins, cette facilité d’utilisation se compense entre autre par des temps de calcul supérieurs à ceux de langages de plus bas niveau, tels que C, Fortran et C++ (mais qui nécessitent eux une compilation).

On dirigera le lecteur curieux vers le livre en ligne Rcpp for everyone de Masaki E. Tsuda, qui constitue une ressource très complète pour comprendre l’utilisation de Rcpp en plus de l’introduction que l’on peut trouver dans le livre Advanced R d’Hadley Wickham.

5.1 Première fonction en Rcpp

A vous de jouer !

  1. Afin de rendre votre package prêt pour l’utilisation avec Rcpp, commencez par executer la commande suivante :

  1. Constatez les changements apportés

  2. il faut également ajouter les 2 commentaires roxygen suivants dans la page d’aide du package dans son ensemble :

Nous allons maintenant créer une première fonction en Rcpp permettant d’inverser une matrice. Pour cela, nous allons nous appuyer sur la library C++ Armadillo. Il s’agit d’une library d’algèbre linéaire moderne et simple, hautement optimisée, et interfacée avec R via le package RcppArmadillo.

C++ n’est pas un langage très différent de R. Les principales différences qui nous concernent :

  • C++est très efficaces pour le boucles for (y compris les boucles for emboîtées). Attention : il y a souvent un sens qui est plus rapide que l’autre (ceci est dû à la manière dont C++ attribue et parcours la mémoire).

  • Chaque commande doit se terminer par un point virgule ‘;’

  • C++est un langage typé : il faut déclarer le type de chaque variable avant de pouvoir l’utiliser.

A vous de jouer !

  1. Créez un nouveau fichier C++ depuis RStudio (via le menu File > New File > C++ File), et enregistrez le dans le dossier src. Prenez le temps de le lire et essayez de comprendre chaque ligne.

  2. Compilez et chargez votre package (via le bouton “Install and Restart”) et essayez d’utiliser la fonction timesTwo() depuis la console.

  3. Installez le package RcppArmadillo, et n’oubliez pas de faire les ajouts nécessaires dans DESCRIPTION (cf. Rcpp précédement - vous pouvez expérimentez avec la fonction RcppArmadillo::RcppArmadillo.package.skeleton() qui a le désavantage de créer beaucoup de fichiers inutiles)

  4. À l’aide de la documentation des packages Rcpp et RcppArmadillo de celle de la library Armadillo, tentez d’écrire une courte fonction invC en C++ calculant l’inverse d’une matrice.

  5. Lorsque vous avez réussi à compiler votre fonction invC et qu’elle est accèssible depuis R créer une fonction mvnpdf_invC() à partir de l’implémentation de mvnpdfsmart en remplaçant uniquement les calculs d’inverse matriciel par un appel à invC.

  6. Evaluer le gain en performance de cette nouvelle implémentation mvnpdf_invC

## Unit: microseconds
##                                                            expr      min
##              mvtnorm::dmvnorm(matrix(1.96, nrow = n, ncol = 2))  143.906
##       mvnpdf(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 5199.792
##  mvnpdfsmart(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 3874.640
##  mvnpdfoptim(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 3030.025
##  mvnpdf_invC(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 3851.251
##        lq      mean   median       uq       max neval
##   209.250  256.7828  262.985  292.773   411.976   100
##  5403.868 6266.2112 5742.189 6390.564 11586.222   100
##  3986.232 4615.3506 4144.007 4545.059  8907.410   100
##  3222.184 3715.9829 3339.993 3623.265 10141.433   100
##  4020.351 4467.4800 4134.856 4469.074  8875.418   100

5.2 Optimisation grâce à C++

En règle générale, on ne gagne pas beaucoup en temps de calcul en remplaçant une fonction R optimisée par une fonction en C++. En effet, la plupart des fonctions de base de R s’appuie en réalité déjà sur des routines C ou Fortran bien optimisée. Le gain se limite alors simplement à la suppression des vérifications des arguments et de la gestion des différents types.

A vous de jouer !

  1. À partir de mvnpdfsmart, proposez une implémentation completement en C++ du calcul de densité de la loi Normale multivariée mvnpdfC().

  2. Evaluer le gain en performance de cette nouvelle implémentation mvnpdf_invC

Vous pouvez télécharger notre proposition de mvnpdfC.cpp ici.

## Unit: microseconds
##                                                                                                  expr
##                                                    mvtnorm::dmvnorm(matrix(1.96, nrow = n, ncol = 2))
##                                             mvnpdf(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE)
##                                        mvnpdfsmart(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE)
##                                        mvnpdfoptim(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE)
##                                        mvnpdf_invC(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE)
##  mvnpdfC(x = matrix(1.96, nrow = 2, ncol = n), mean = rep(0, 2),      varcovM = diag(2), Log = FALSE)
##       min       lq       mean    median        uq       max neval
##   153.970  215.692  248.40744  251.3340  274.9605   403.513   100
##  5159.574 5303.610 5942.66552 5624.4390 6267.0645 10993.952   100
##  3799.345 3912.713 5002.47948 3997.8350 4401.0525 65655.801   100
##  2965.691 3091.124 3576.33037 3209.7125 3433.1545 10726.958   100
##  3774.164 3884.534 4576.82329 3985.9375 4211.6280 14368.164   100
##    33.316   44.611   52.14602   50.1765   58.6415    98.589   100

À noter que vous pouvez utiliser des fonctions Rcpp en dehors de l’architecture d’un package grâce à la fonction Rcpp::sourceCpp(). Mais comme nous avons qu’il est préférable de gérer tous ces code sous la forme de package, il est peu probable que vous en ayez besoin !

5.3 Annexe 5.1 : l’optimisation prématurée n’est pas une bonne idée

Chambers, Software for Data Analysis: Programming with R, Springer, 2008 :

Including additional C code is a serious step, with some added dangers and often a substantial amount of programming and debugging required. You should have a good reason.