Ecrire des belles fonctions de valorisation c’est bien… mais il faut aussi penser “performance”

(je vous vois vous frotter les mains, désolé je parle de performance technique/informatique pas financière)

Dans le dernier article Il est l’hooore, mon seignooor (d’écrire une fonction de valorisation) nous avons exploré plusieurs façons d’écrire notre fonction de valorisation.

Et vous avez peut-être remarqué que je me suis arrêté un peu abruptement. C’est pour pouvoir vous expliquer un détail CAPITAL, qui vaut bien un article à lui tout-seul, lorsqu’on traite des fichiers reprenant des dizaines de milliers d’enregistrements tels des fichiers PMSI.

Dans le précédent article, le principe de calcul paraissait simple avec diverses implémentations plus ou moins élégantes, certes. Cependant il ne faut pas perdre de vue qu’au regard de la taille potentielle de nos fichiers d’export, le temps d’exécution de cette évaluation – appelée à être relancée des dizaines de milliers de fois par traitement – n’est pas anodin.

Résumé de l’épisode précédent

Nous avons défini plusieurs variations de la fonction de valorisation :

• valoriseV() : la Vectorize-ation d’une fonction initialement non vectorisée nécessitant une jointure préalable
• valorise2() : une fonction vectorisée nativement nécessitant une jointure préalable
• valorise3() : une version contractée en code de valorise2 (utilisant min/max à la place de ifelse)
• valorise5() : une version vectorisée avec un traitement de liste sans jointure préalable mais passage de la table en paramètre
• valorise.ghs_df() : un version appelable directement depuis la table de GHS/DS sans passer par un mutate(), qui dans l’article précédent utilisait valorise5()

(valorise4 n’a pas été écrite)

Tester la performance

Pour tester la vitesse d’exécution de fonctions, il existe une librairie très pratique : microbenchmark important la fonction du même nom.

Pour tester le temps d’exécution d’un code, on fait :

microbenchmark(fonction_a_tester1,fonction_a_tester2,...)

Il est possible de compléter par le nombre de fois que l’on veut faire tourner les fonctions via le paramètre times= qui vaut 100 par défaut. Attention :

• les fonctions sont réellement évaluées et si elles modifient des données, celles-ci seront modifiées times fois.
• fonction_a_tester(…) doit être une expression
• Enfin, si les traitements sont longs, le benchmark va durer :
  <nombre de fonction_a_tester> * <times> * <durée unitaire du traitement>
  Ca peut lasser…

Pour illustrer, lançons le test sur des expressions triviales :

> microbenchmark(1 + 2, 1 - 2, 1 / 2, 1 * 2)
Unit: nanoseconds
  expr min lq  mean median uq  max neval
 1 + 2   0  1 53.98      1  2 2300   100
 1 - 2   0  1 18.06      1  2  401   100
   1/2   0  1 16.95      1  2  400   100
 1 * 2   0  1 18.00      1  2  301   100

La colonne qui nous intéresse le plus est bien sûr “mean” qui donne la durée d’exécution moyenne (les autres : “min” = durée minimale, “max” = durée maximale, “median” = durée médiane, “lq/uq” = lower/upper quartile), sur 100 évaluations (“neval”, valeur par défaut) de chacune des expressions passées en paramètre. L’unité est précisée en haut du tableau, ici la nanoseconde, elle est réglable via unit = "seconds" par exemple, mais s’adapte automatiquement à l’ordre de grandeur des résultats.

Passons au cas réel

Nous allons vouloir mettre en compétition valorise[V, 2, 3 et 5] ainsi que la version “5_df”. Il nous faut donc écrire :

microbenchmark(
     RSA %>% left_join(GHSs, by = c("GHS")) %>%
             mutate(valo = valoriseV(GHS, DS, BB, BH, GHSPRIX, FORFAIT_EXB, EXBPRIX, EXHPRIX)),
     RSA %>% left_join(GHSs, by = c("GHS")) %>%
             mutate(valo = valorise2(GHS, DS, BB, BH, GHSPRIX, FORFAIT_EXB, EXBPRIX, EXHPRIX)),
     RSA %>% mutate(valo = valorise5(GHS,DS,GHSs)),
     RSA %>% valorise5_df(GHSs),
times=2)

Je n’ai lancé le benchmark que pour 2 tours car le résultat est flagrant (je vous conseille d’en faire autant lors des premiers essais au risque de devoir arrêter brutalement le traitement par épuisement de patience). Aussi beau et astucieux que soit le code de valorise5(), son exécution est une catastrophe (je n’ai gardé que les moyennes pour la lisibilité) :

Nous passons de la méthode “jointure” à la méthode “apply” par un rapport de plus de x1700 ! Un RSA de 25.000 lignes durerait de 0.125 secondes à 3 minutes (l’échantillon du test était de 1000 lignes).

On note aussi que la version “auto”-vectorisée met tout de même 2 fois plus de temps que la version “à la main” (c’est normal, Vectorize utilise un apply() en interne) et que la version à base de min()/max() est plus longue que celle à base d’ifelse().

La conclusion est sans appel, il ne faut PAS utiliser un apply() (ou une boucle, c’est encore pire) si une jointure fonctionne, aussi “belle” soit la solution, et on écrit soi-même ses fonctions vectorisées si les volumes de traitements sont importants.

Victoire donc à :

valorise2 <- function(ghs, ds, bb, bh, ghsprix, forfait_exb = 0, exb, exh){
    nexb = bb - ds
    nexh = ds - bh
    vexb = ifelse(nexb>0, ifelse(forfait_exb > 0,
                                 forfait_exb,
                                 nexb * exb
    ), 0)
    vexh = ifelse(nexh > 0, nexh * exh, 0)
    
    ghsprix - vexb + vexh
}

RSA %>% left_join(GHSs, by = c("GHS")) %>% mutate(valo = valorise2(GHS, DS, BB, BH, GHSPRIX, FORFAIT_EXB, EXBPRIX, EXHPRIX))

Mais pour le coup, nous n’avons plus l’écriture simplifiée RSA %>% valorise.ghs_df(...). Bon, et bien…

Retour à la planche à dessin

Nous avons donc vu que valorise2() est la fonction la plus rapide, nous allons donc essayer de la réemployer pour disposer d’une fonction agissant directement sur une table. Seul la ligne produisant la nouvelle table de données va changer :

valorise2_df <- function (x, ghss, ghs_col = "GHS", ds_col = "DS", out_col="valo"){

  if (out_col %in% names(x)) stop("la colonne de sortie `", out_col, "`` existe déjà.")
  if (!(ghs_col %in% names(x))) stop("la colonne du GHS `", ghs_col, "` est absente du jeu source.")
  if (!(ds_col %in% names(x))) stop("la colonne de durée de séjours `", ds_col, "` est absente du jeu source.")
  suppressWarnings(if (is.na(ghss) | !("data.frame" %in% class(ghss))) error("ghss doit être une data.frame au format attendu."))

# inchangé jusque ici

  bind_cols(x,
            x %>% select(all_of(c(ghs_col, ds_col))) %>%
              left_join(ghss, by = setNames(ghs_col, "GHS")) %>%
    mutate("{out_col}" := valorise2(.data[[ghs_col]],.data[[ds_col]],BB, BH, GHSPRIX, FORFAIT_EXB, EXBPRIX, EXHPRIX)) %>% select(out_col)
    )
}

Décortiquons un peu cette ligne :

• bind_cols(...) colle les n dataframes passées en paramètres l’une à l’autre en ajoutant, sans jointure, les colonnes les unes à la suite des autres. Dans la situation actuelle, il va coller au jeu de départ x la colonne de valorisation calculée dans la seconde table. (Les data.frames étant généralement sur le même niveau sémantique, on n’utilise pas la notation x %>% bind_cols(...), mais on pourrait surtout lorsqu’on ne rajoute comme ici qu’une colonne)
• x : la table source, inchangée, que l’on a reçue en paramètre
• x %>% select(all_of(c(ghs_col, ds_col))) : isole du jeu source les 2 colonnes qui nous intéressent. Le but est d’éviter les doublons de colonnes si par exemple dans le RSA, on a déjà collé `BB`, il risque d’y avoir ambiguïté et comportement inattendu…
• %>% left_join(...) : rien d’exceptionnel. On utilise setNames() en son sein pour obtenir le c("GHS" = "GHS") (si vous avez suivi l’article sur les types de base, c’est un vecteur de chaine de caractères auquel on a ajouté un nom à l’unique valeur) à partir de la donnée ghs_colpour la partie à gauche.
• %>% mutate(...): C’est grosso-modo le même principe que pour valorise.ghs_df() de l’article précédent.
• %>% select(out_col) : on ne garde que la colonne de résultat pour l’ajouter au jeu source via le bind_cols()

Nous passons alors de 4.8 millisecondes (valorise2) à 13 millisecondes (valorise2_df), c’est certes plus long mais reste raisonnable, cela nous fait une vitesse de valorisation 77.000 RSA par seconde.

Monsieur Plus

Pour la petite histoire, on pourrait se passer partiellement ou complètement de dplyr mais au prix d’une lisibilité moindre et gagner encore quelques millisecondes à tel point qu’on devient équivalent à la version “jointure/mutate” :

# Il est nécéssaire de faire une adaptation subtile de valorise2 
# en rajoutant des unlist() englobant les ifelse() pour retransformer
# les lists produites par ifelse() en vectors unitaires numériques.

valorise2 <- function(ghs, ds, bb, bh, ghsprix, forfait_exb = 0, exb, exh){
  
  
  nexb = bb - ds
  nexh = ds - bh
  vexb = unlist(ifelse(nexb>0, ifelse(forfait_exb > 0,
                               forfait_exb,
                               nexb * exb
  ), 0))
  vexh = unlist(ifelse(nexh > 0, nexh * exh, 0))
  
  ghsprix - vexb + vexh
}

valorise7_df <- function (x, ghss, ghs_col = "GHS", ds_col = "DS", out_col="valo"){
  if (out_col %in% names(x))
     stop("la colonne de sortie `", out_col, "`` existe déjà.")
  if (!(ghs_col %in% names(x)))
     stop("la colonne du GHS `", ghs_col, "` est absente du jeu source.")
  if (!(ds_col %in% names(x)))
     stop("la colonne de durée de séjours `", ds_col, "` est absente du jeu source.")
  suppressWarnings(
     if (is.na(ghss) | !("data.frame" %in% class(ghss))) 
        error("ghss doit être une data.frame au format attendu.")
     )
  
  y <- x[, c(ghs_col, ds_col)] %>%
       merge(ghss, by.x = ghs_col, by.y = "GHS", all.x = TRUE)
  x[,out_col] <-  valorise2(NA,        # Le GHS est inutile
                            y[, ds_col],
                            y[, "BB"],
                            y[, "BH"],
                            y[, "GHSPRIX"],
                            y[, "FORFAIT_EXB"],
                            y[, "EXBPRIX"],
                            y[, "EXHPRIX"])
   x
  }

Le benchmark de RSA %>% valorise7_df(GHSs) donne alors 5.18 millisecondes contre 5.32 pour left_join() %>% mutate(valorise2()) (on a très légèrement dégradé valorise2 en rajoutant les unlist()) .

Attention cependant, merge() est moins efficace que *_join() sur les grands sets tels qu’un RSA annuel, il vaut donc mieux garder la solution de jointure par *_join()car nous y gagnerons sur des vrais RSA.

Capillotractage

Pour la beauté du geste, gagnons encore quelques 1000èmes de seconde en vous présentant en remplacement de valorise2(), une approche exotique mais plutôt séduisante évitant les if_else() et les unlist() résultants🤪:

valorise8 <- function(ghs, ds, bb, bh, ghsprix, forfait_exb = 0, exb, exh){
  dna = !is.na(ghsprix)
  dexb=  dna & (ds < bb) & (forfait_exb == 0)
  dexh = dna & (ds > bh) 
  dforfait = dna & (ds < bb) & (forfait_exb >= 0) 
  dghs = dna & (ds >= bb) & (ds <= bh)
  
  x<-double()
  x[!dna]      <- NA
  x[dghs]     <- ghsprix[dghs]
  x[dforfait] <- ghsprix[dforfait] - forfait_exb[dforfait]
  x[dexb]     <- bb[dexb] - ds[dexb]
  x[dexh]     <- ds[dexh] - bh[dexh]
  x
}

(nous créons des vecteurs logiques afin d’appliquer les expressions définies uniquement aux subsets correspondants à la situation de chaque séjour). Sur un ordinateur puissant cela me fait gagner encore 15% sur un set d’environ 21K lignes de RSA ; sur mon professionnel -moins puissant et moins doté en mémoire- la différence est > 50%.
[sur mon ordinateur personnel, je monte à un potentiel de plus de 2 millions de RSA valorisés par seconde sur 1 cœur/1 thread, autant dire instantané… et la charge annexe de l’ordinateur (bouger la souris par exemple fait perdre… 4300 RSA par seconde) aura plus d’impact que l’optimisation]

Là, je commence à ne pas avoir plus rapide sans passer par l’utilisation de code C++ ! Sinon, il ne resterait que la parallélisation mais avec des réserves sur les “petits” volumes et le traitement lui-même qui reste rapide, le tri/préchargement des données dans chaque process parallèle rendant la procédure probablement peu intéressante eut égard aux gains anticipables…

Conclusion

Cet article nous a permis de toucher du doigt le problème de la performance liée aux traitements de masse, à leur optimisation pour R et proposé quelques pistes de prise en charge.

C’ était un clin d’oeil à la plus désagréable des situations quand on génère ou retraite des fichiers PMSI : Rester devant son écran et… attendre… que… les… calculs… se… terminent… (on te voit module de “commande de traitement” d’e-pmsi le 31 janvier !). Alors que comme nous venons de le voir, une adaptation du code permet de faire des gains monstrueux en terme de vitesse de traitement.

Nous avons enfin notre fonction de valorisation finale d’un séjour, qui de plus s’exécute dans un temps raisonnable, cependant il reste encore un détail à régler : l’application de l’arrêté tarifaire. Ce sera le sujet du prochain article. Seulement alors pourrons-nous reprendre le cours de notre évaluation de la RAAC.