Skip to main content

Formation avancée: Outils de développement et de performance

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

Rcpp (R-C-Plus-Plus) est un package qui facilite l’interface entre C++ et . 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.

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 :
```r
usethis::use_rcpp()
```
  1. Constatez les changements apportés
  2. Comme mentionner dans la console, il faut également ajouter les 2 commentaires Roxygen suivants dans la page d’aide globale du package :
```r
#' @useDynLib mypkgr
#' @importFrom Rcpp sourceCpp, .registration = TRUE
NULL
```

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 via le package RcppArmadillo.

C++ n’est pas un langage très différent de . 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. usethis::use_rcpp_armadillo())
  4. À l’aide de l’introduction à Rcpp de Hadley Wickham dans son livre Advanced R 6, ainsi que de la documentation du package RcppArmadillo et de celle de la library C++ 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 accessible depuis 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.
n <- 1000
mb <- microbenchmark(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),
                     times=100L)
mb
## Unit: microseconds
##                                                            expr      min
##              mvtnorm::dmvnorm(matrix(1.96, nrow = n, ncol = 2))   44.198
##       mvnpdf(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 3150.358
##  mvnpdfsmart(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 2329.251
##  mvnpdfoptim(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 1763.984
##  mvnpdf_invC(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 2316.746
##        lq       mean   median       uq      max neval  cld
##    55.842   75.18457   77.203   90.897  122.795   100 a   
##  3287.462 3530.78511 3356.649 3490.576 9381.292   100  b  
##  2363.015 2554.23973 2390.361 2435.851 7216.287   100   c 
##  1810.499 1979.80062 1848.157 1889.198 7138.059   100    d
##  2353.584 2499.20133 2380.337 2433.658 6020.686   100   c
profvis::profvis(mvnpdfoptim(x=matrix(1.96, 
    nrow = 2, ncol = 1000), Log=FALSE))
## Error in parse_rprof_lines(lines, expr_source): No parsing data available. Maybe your function was too fast?
profvis::profvis(mvnpdfoptim(x=matrix(1.96, 
    nrow = 100, ncol = 1000), Log=FALSE))

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 optimisée par une fonction en C++. En effet, la plupart des fonctions de base de 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 mvnpdfC

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

Pour un gain (relativement faible small) supplémentaire de temps de calcul (au prix d’un code plus difficile à lire), vous pouvez jettez un oeil à notre implémentation optimizée utilisant C++ et Armadillo dans le fichier mvnpdfoptimC.cpp.

n <- 1000
mb <- microbenchmark(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),
                     mvnpdfsmartC(x=matrix(1.96, nrow = 2, ncol = n), mean = rep(0, 2), varcovM = diag(2), Log=FALSE),
                     mvnpdfoptimC(x=matrix(1.96, nrow = 2, ncol = n), mean = rep(0, 2), varcovM = diag(2), Log=FALSE),
                     times=100L)
mb
## 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)
##  mvnpdfsmartC(x = matrix(1.96, nrow = 2, ncol = n), mean = rep(0,      2), varcovM = diag(2), Log = FALSE)
##  mvnpdfoptimC(x = matrix(1.96, nrow = 2, ncol = n), mean = rep(0,      2), varcovM = diag(2), Log = FALSE)
##       min        lq       mean    median        uq      max neval  cld
##    44.362   57.2360   75.42319   77.8590   88.7035  113.365   100 a   
##  3173.359 3297.3840 3516.53966 3380.4295 3508.6980 8638.864   100  b  
##  2335.155 2370.7840 2451.90824 2397.0035 2427.5690 6993.452   100   c 
##  1765.337 1810.7240 2002.75488 1836.4720 1889.8950 6392.761   100    d
##  2321.871 2358.7710 2565.19452 2385.4620 2421.3780 7517.391   100   c 
##    51.496   55.4320   61.39750   60.3725   65.6615   88.806   100 a   
##    35.834   38.9295   45.58872   43.9315   50.8810  112.873   100 a

À 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 l’avons vu, il est préférable de gérer tous ses codes sous la forme de package : il est donc peu probable que vous en ayez besoin !

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