Le code source de R

Je commence à connaître un peu le code source de R, en raison d’améliorations que j’y ai apportées : possibilité d’interrompre les gros calculs, accélération de la génération de nombre aléatoires, suppression de la limite à 128 connections, accélération de la fonction mean() sur les ALTREPs, refactoring de la fonction sum().

Voilà mon opinion : R contient des algorithmes conçus pour une précision numérique maximale, souvent au détriment des performances, ce qui se comprend étant donné les objectifs du projet. Par exemple, dès que possible, des nombres à virgule flottante 80 bits sont utilisés dans les calculs intermédiaires, quand bien même le résultat final est seulement stocké sur 64 bits. Cependant, la structure du code source est vraiment mauvaise. Cela rend le code moins maintenable et plus obscur. Voilà, les éléments que j’ai constatés en lisant le code source des quelques fichiers dont j’ai eu besoin pour mes patches:

Un peu de Fortran 77

Pour commencer par un des aspects les plus acceptables du code source, celui-ci contient du code Fotran 77 avec le style de l’époque, notamment le respect de la limite à 6 caractères pour la longueur des identifiants, y compris des variables locales. Ce n’est pas très agréable à lire, mais heureusement, cela ne correspond qu’à une minorité du code source, et le code est basé sur des interfaces déjà connues (telles que BLAS/LAPACK, LINPACK), même si cela a engendré des problèmes d’interopérabilité Fortran/C.

Il contient aussi un peu de Fortran 90 qui le rend impossible à compiler, pour l’instant, sur Apple M1; d’ici quelques mois il est probable que GCC soit porté sur cette nouvelle plateforme.

Des gros switch() partout

Dans le fichier src/main/RNG.c contenant les 8 algorithmes de génération de nombres pseudo-aléatoires (dont l’algorithme « user-supplied »), les algorithmes sont implémentés de manière éparse dans diverses fonctions.

C’est-à-dire qu’un algorithme devrait implémenter:

  1. Une fonction d’initialisation de l’état du générateur à partir d’une graine pseudo-aléatoire 32 bits (RNG_Init)
  2. Une fonction de remise en cohérence de l’état du générateur au cas où l’utilisateur aurait modifié l’état interne par un accès direct à l’objet .Random.seed (FixupSeeds).
  3. Une fonction de génération de nombre pseudo-aléatoire proprement dite (unif_rand)
  4. Des informations telles que le nombre de bits de précision effective du nombre pseudo-aléatoire

Plutôt que de définir une interface, sous forme d’une structure avec pointeurs vers fonction, l’ensemble des fonctionnalités est implémentée sous forme de gros switch() pour chacune des méthodes.

double unif_rand(void)
{
    double value;

    switch(RNG_kind) {

    case WICHMANN_HILL:
        I1 = I1 * 171 % 30269;
        I2 = I2 * 172 % 30307;
        I3 = I3 * 170 % 30323;
        value = I1 / 30269.0 + I2 / 30307.0 + I3 / 30323.0;
        return fixup(value - (int) value);/* in [0,1) */

    case MARSAGLIA_MULTICARRY:/* 0177777(octal) == 65535(decimal)*/
        I1= 36969*(I1 & 0177777) + (I1>>16);
        I2= 18000*(I2 & 0177777) + (I2>>16);
        return fixup(((I1 << 16)^(I2 & 0177777)) * i2_32m1); /* in [0,1) */

/* ... */

Cela rend l’ajout d’algorithmes très délicat, parce qu’il faut rajouter une entrée dans chacun des switch() de src/main/RNG.c.

Par ailleurs, la liste des générateurs pseudo-aléatoires est identifié par une énumération:

typedef enum {
    WICHMANN_HILL,
    MARSAGLIA_MULTICARRY,
    SUPER_DUPER,
    MERSENNE_TWISTER,
    KNUTH_TAOCP,
    USER_UNIF,
    KNUTH_TAOCP2,
    LECUYER_CMRG
} RNGtype;

Mais dans la fonction PutRNGstate, il est vérifié que le numéro du générateur est « valide »; cette vérification est simplement la comparaison de RNG_king à la valeur LECUYER_CMRG. En bref, le fait que l’algorithme LECUYER_CMRG est le dernier de la liste est hard-codé dans le code source, obligeant à modifier cette partie du code source lors de l’ajout d’un nouveau générateur pseudo-aléatoire.

Des constantes magiques

Les fonctions R de base (primitives) mean, sum, prod, min et max sont implémentées dans la même fonction C : do_summary dans src/main/summary.c. Les cinq fonctions (mean, sum, prod, min et max) sont distinguées par un numéro d’opération (0=sum, 1=mean, 2=min, 3=max, 4=prod) et des gros switch() disséminés dans la fonction do_summary adaptant le comportement à l’opération. Plutôt que de définir une énumération permettant d’associer un numéro à un nom (p.e. MIN_IOP = 2), le code contient des constantes magiques, 0, 1, 2, 3 et 4, conduisant à ce genre de bout de code :

    if(empty && (iop == 2 || iop == 3)) {
	if(ans_type == STRSXP) {
	    warningcall(call, _("no non-missing arguments, returning NA"));
	} else {
	    if(iop == 2)
		warningcall(call, _("no non-missing arguments to min; returning Inf"));
	    else
		warningcall(call, _("no non-missing arguments to max; returning -Inf"));
	    ans_type = REALSXP;
	}
    }

Si on ne sait pas que 2=min et 3=max, on peut difficilement comprendre ce code.

Dans la même fonction do_summary, il y a une variable « empty », partagée entre toutes les opérations, valant TRUE tant qu’on n’a rencontré que des vecteurs vides (longueur zéro) ou ne présentant que des données manquantes (NA) éliminées par un na.rm=TRUE. Cela permet aux fonctions min et max d’émettre un message d’avertissement dans ce genre de cas:

> min(c(NA,NA),NULL,na.rm=TRUE)
[1] Inf
Warning message:
In min(c(NA, NA), NULL, na.rm = TRUE) :
  no non-missing arguments to min; returning Inf

Lors du traitement de chaque paramètre de la fonction, une variable « updated » est mise à jour afin de dire si des valeurs non manquantes (sauf si na.rm=FALSE) ont été intégrées dans l’opération. Eh bien, dans certaines situations la variable updated, qui est booléenne à la base, peut devenir un trooléen pour décrire une situation de dépassement de capacité (overflow) des nombres entiers 64 bits, concernant la somme de nombres entiers. La troisième valeur, qui n’est ni TRUE, ni FALSE, vaut 42, sans la moindre explication quant à ce choix, même si on devine qu’il y a un rapport direct à la grande question sur la vie, l’univers et le reste.

Du code dupliqué en fonctionnalité

Résumé

R contient deux implémentations de fonctions associées aux distributions aléatoires (rnorm, rbinom, …), différant de manière subtile et accidentelle en comportement. Une des deux implémentations n’est jamais utilisée, sauf dans le cas où l’utilisateur outrepasserait volontairement les fonctions proposées par R (p.e. stats::rnorm), afin d’accéder directement aux fonctionnalités internes (.Internal(rnorm(a,b,c)).

Détail

L’implémentation des fonctions de hasard est séparée en deux parties :

  1. D’un côté, les algorithmes pseudo-aléatoires élémentaires, implémentés dans src/main/RNG.c : ceux-ci permettent de générer des nombres selon une loi uniforme entre 0 et 1 (à quelques détails près) ;
  2. D’un autre côté, des algorithmes permettant de générer des distributions particulières (p.e. binomiale, Student, exponentielle, etc.), conçues pour être directement accessibles depuis du code R (rbinom, rt, rexp, etc.).

En apparence, la deuxième partie est implémentée dans le fichier src/main/random.c, sauf que la grande majorité du code de ce fichier ne semble servir à rien ! Il ne semble jamais être appelé, jamais être utilisé ! La seule partie utilisée concerne la fonction R sample(), mais tout ce qui concerne rchisq, rexp, rgeom, rpois, rt, rsignrank, rbeta, rbinom, rcauchy, rf, rgamma, rlnorm, rlogis, rnbinom, rnorm, runif, rweibull, rwilcox, rhyper, n’est jamais appelé par du code R (nous verrons par la suite quand ça peut encore être utilisé). J’ai pu supprimer 257 lignes de code de src/main/random.c, sans perte évidente de fonctionnalité de R. Comment cela est-il possible ? Eh bien, ces fonctionnalités ont été réimplémentées dans src/library/stats/src/random.c, compilé dans une bibliothèque dynamique stats.so ou stats.dll qui sera ensuite chargée en même temps que le packages ‘stats’ qui fait partie des package pré-chargés de R. Cette réimplémentation est quasiment identique, en fonctionnalité, à l’implémentation de src/main/random.c, avec néanmoins de très subtiles différences, telles que la fonction rbinom() qui renverra des nombres entiers ou flottants, selon qu’ils dépassent ou pas INT_MAX (2147483647) dans src/library/stats/src/random.c, alors qu’elle renverra toujours des flottants dans src/main/random.c. Même si les deux implémentations sont presque identiques en fonctionnalité, elles diffèrent nettement en code source, avec néanmoins des points de convergence assez nets.

On peut en comprendre l’histoire, en remontant l’historique des commits. En avril 2012, il n’y avait pas de doublon de fonctionnalité, puisque src/library/stats/src/random.c ne contenait que le code de rmultinom et r2dtable qui ne sont pas implémentés dans src/main/random.c. En mai 2012, le code de src/main/random.c a été dupliqué dans src/library/stats/src/random.c. Déjà, à l’époque, la duplication du code n’était pas fidèle : c’était une réimplémentation compatible mais de style différent. Il y avait déjà des divergences de fonctionnalité apparente; notamment rbinom() renvoyait toujours des nombres entiers avec l’implémentation de src/library/stats/src/random.c, avec possible overflow, alors que le code de src/main/random.c renvoyait des nombres à virgule flottante. Par la suite les deux implémentations ont divergé de plus en plus, y compris avec des patches qui modifiaient le code de src/main/random.c sans savoir qu’il n’était pas vraiment utilisé. Le grand responsable semble donc être le patch de mai 2012. En même temps que de dupliquer tous les algorithmes de génération de nombre aléatoires (rbinom, rchisq, …), ce patch a aussi dupliqué tous les algorithmes de densité de probabilité, de fonction de répartition et de fonctions des quantiles (dnorm, pnorm, qnorm, dbinom, dbinom, …) en créant le fichier src/library/stats/src/distn.c qui faisait encore 481 lignes de fonctionnalité dupliquée. Ce patch semble correspondre au CHANGELOG de R 3.0.0 qui précise (citation exacte):

CODE MIGRATION

  • The C code underlying base graphics has been migrated to the graphics package (and hence no longer uses .Internal() calls).
  • Most of the .Internal() calls used in the stats package have been migrated to C code in that package. This means that a number of .Internal() calls which have been used by packages no longer exist, including .Internal(cor) .Internal(cov), .Internal(optimhess) and .Internal(update.formula).
  • Some .External() calls to the base package (really to the R executable or shared library) have been moved to more appropriate packages. Packages should not have been using such calls, but some did (mainly those used by integrate()).

Les fonctions .Internal sont celles que l’on retrouve dans le dossier src/main, telles que les fonctions de src/main/random.c. Cela suggère donc une migration de fonctionnalité dont l’intérêt n’est pas complètement évident, mais cela n’explique pas la duplication de fonctionnalité, qui a mené à de dangereuses divergences de code par la suite. Pourquoi cette duplication ? On a un élément de réponse dans le message venant avec le commit de mai 2012:

port C code for [dpqr]xxx functions to stats
leave .Internals for now, as they are used directly as an optimization by byte-compiler

Ainsi, le code inutilisé en apparence pourrait encore être utilisé par le byte-compiler; sauf que c’est faux, bien heureusement. Le byte-compiler est une fonction introduite dans R 2.13.0, avec la précompilation de tous les packages par défaut à partir de R 2.14.0. C’est R 2.15.0 et R 3.0.0 que le patch de mai 2012 a été introduit. Cette compilation convertit le code R en code de machine virtuelle basée sur une pile (stack based virtual machine), avec une amélioration très légère en performances, notamment en raison de certaines optimisations. Ces optimisations consistent en le remplacement de certaines fonctions de bases (p.e. addition, boucle for, parenthèse ouvrante, accolade ouvrante) par une référence directe à la fonction interne (builtin ou special), outrepassant les accès aux environnements de travail, qui sont très chronophages. Dans l’idée du commentaire du commit de mai 2012, le byte-code permettrait l’accès direct aux fonctions telles que rnorm() sans passer par la version réécrite pour le package stats.

Mais cela est presque totalement faux, pour plusieurs raisons :

  1. Le code responsable de la compilation dans src/library/compiler/R/cmp.R n’outrepassera le code de la fonction stats::rnorm() et autres fonctions aléatoires seulement dans le cas où le corps de celle-ci serait limité à un appel direct à .Internal, grace à la vérification faire dans is.simpleInternal. Ainsi, le byte-code compilé par R 3.0.0 ne fera pas d’optimisation et appellera toujours la version implémentée dans src/library/stats/src/random.c.
  2. R 3.0.0 refusera de charger un package compilé pour R 2.15, évitant ainsi tout risque d’utiliser du byte-code optimisé pour la version précédente et faisant un accès direct à la version interne implémentée dans src/main/random.c.

Il reste théoriquement possible d’exécuter l’ancienne version avec la procédure suivante

  1. Appeler compiler::setCompilerOptions(optimize=3) sous R 2.15
  2. Compiler du byte-code toujours sous R 2.15 et l’exporter dans un fichier rds avec saveRDS()
  3. Charger ce byte-code sous R 3.0.0 avec readerRDS()
  4. Évaluer ce code avec eval()

Compiler du code sous une version de R et l’exécuter sous la suivante est à coup sûr une mauvaise idée de toute façon. Évidemment, il existe aussi la possibilité qu’un programme accède directement et explicitement à la fonction interne avec un appel à .Internal. Il n’est pas étonnant que ce genre de programme se casse dans une nouvelle version de R. En pratique, le byte-code compilé avec R 2.15 s’exécutera correctement sur une version récente de R (p.e. 4.1.0) parce que le byte-code ayant évolué entre temps, R ignorera le byte-code de l’expression compilée et se contentera d’utiliser l’expression brute, non compilée, évitant alors d’exécuter la version de src/main/random.c.

Dans tout les cas, la solution retenue est la pire qui soit. Dans les très rares cas où la fonction interne de src/main/random.c serait appelée, par l’exécution de code byte-compilé pour la version précédente de R, ou à cause d’un appel direct par .Internal, le comportement va diverger de celui de la fonction standard (maintenant src/library/stats/src/random.c) en exécutant un code C de plus en plus mal testé au gré des versions successives de R.

Les solutions pertinentes auraient été:

  1. Supprimer le code de src/main/random.c et les fonctions internes correspondantes afin de ne conserver que la nouvelle implémentation avec src/library/stats/src/random.c
  2. Remplacer le code de src/main/random.c par des redirections vers les implémentations de src/library/stats/src/random.c conservant ainsi la compatibilité avec les cas d’usage les plus bizarres
  3. Ne pas migrer le code dans src/library/stats/src/random.c

La même réflexion s’applique aux fonctions dupliquées dans src/library/stats/src/distn.c et src/main/arithmetic.c.

On peut aussi s’étonner de l’incohérence avec la migration d’autres parties du code, telles que les fonctions cor et cov dont l’implémentation interne dans src/main avait été totalement supprimée lors de la migration. La raison semble être due à la tentative d’optimisation du byte-code qui n’existait pas pour ces dernières fonctions.

J’ai mentionné ce problème sur la mailing list r-devel, mais ça ne semble pas intéresser les développeurs R, le message étant resté sans réponse.

Conclusion

Le code de R aurait besoin d’un refactoring important, mais ce n’est pas le style de la maison. Les recommandations de développement de R sont claires :

DO NOT fix exotic bugs that haven't bugged anyone
DO NOT fix things that are not broken
DO NOT restructure code

Mise à jour

J’ai rapporté le bug du code dupliqué sur le bugzilla de R et ai rapidement reçu la réponse des développeurs de R:

WONTFIX « I prefer to leave this as is for now.« 

Ce qui signifie que l’équipe de développement a conscience de ce problème mais préfère conserver le code dupliqué, tel qu’il est; probablement par peur de créer un problème de compatibilité, sans bénéfice immédiat évident. Cela est en cohérence avec la politique :

DO NOT fix things that are not broken
DO NOT restructure code

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée. Les champs obligatoires sont indiqués avec *