Pythran 0.6 - compilation de noyaux scientifiques écrits en Python

Pythran est un compilateur pour les noyaux de calcul scientifique écrit en Python. Il permet d'écrire des modules dans un large sous-ensemble de Python + Numpy, d'ajouter quelques lignes de commentaire pour spécifier les types des fonctions exportées, enfin de compiler l'ensemble pour obtenir un module natif capable (parfois !) d'utiliser efficacement multi-cœurs et unités vectorielles. Le reste de la dépêche décrit le fonctionnement du compilateur, les évolutions récentes et propose une comparaison avec les alternatives : Cython, numba et parakeet.

Sommaire

• Exemple d'utilisation de Pythran
• Fonctionnement interne de Pythran
• Performance du code généré
• Nouveauté dans la version 0.6
• La suite ?

Exemple d'utilisation de Pythran

Prenons un petit code de calcul classique, extrait de la suite de validation de parakeet

import numpy as np
def harris(I):
  m,n = I.shape
  dx = (I[1:, :] - I[:m-1, :])[:, 1:]
  dy = (I[:, 1:] - I[:, :n-1])[1:, :]

  #
  #   At each point we build a matrix
  #   of derivative products
  #   M =
  #   | A = dx^2     C = dx * dy |
  #   | C = dy * dx  B = dy * dy |
  #
  #   and the score at that point is:
  #      det(M) - k*trace(M)^2
  #
  A = dx * dx
  B = dy * dy
  C = dx * dy
  tr = A + B
  det = A * B - C * C
  k = 0.05
  return det - k * tr * tr

C'est un bon exemple de code scientifique écrit en Python : c'est du haut niveau, on travaille pas mal sur les tableaux (à la FORTRAN ;-)). Le code est un peu polymorphique (la fonction marche sur des tableaux de float ou de double par exemple), mais il y a un contrat implicite (sur la dimension de I). On attend généralement d'un bon développeur Python + Numpy qu'il écrive - quand c'est possible et que cela reste lisible - du code de haut niveau qui aurait cette forme.

Pour compiler ce code avec Pythran, trois choses à faire :

1. le mettre dans un fichier / module à part, disons speedy.py, éventuellement avec d'autres fonctions compatibles avec Pythran. Gageons que cela permette de séparer la partie calcul du reste de l'application (<< oui, ceci est une tentative maladroite de justifier une contrainte imposée par Pythran !) ;

2. ajouter une ou plusieurs lignes pour informer Pythran des versions de cette fonction à générer, par exemple :

#pythran export harris(float64[][])
#pythran export harris(float32[][])

3. compiler le module en module natif :

$ pythran speedy.py

cette ligne aura pour effet de générer un module natif speedy.so que l'on peut importer comme un module « normal »

from speedy import harris

Fonctionnement interne de Pythran

Pythran est un compilateur statique, par opposition à un compilateur faisant de la compilation à la volée, ce qu'on trouve majoritairement dans le monde des langages interprétés (voir par exemple PyPy, Pyston, mais numba et parakeet comptent aussi). Cython est lui aussi un compilateur statique. Le prix à payer est une petite ligne pour déclarer le type de la fonction exportée (mais pas des fonctions intermédiaires!), on peut espérer y gagner en temps d'exécution (on a plus de temps pour compiler, donc on peut se permettre des analyses / optimisations plus coûteuses), même si on dispose de moins d'informations (aucune information de contexte sur la valeur des paramètres).

Une séance de compilation (c'est un peu la séance de massage du .py ;-)) se déroule ainsi :

1. Transformation du .py en un Arbre syntaxique abstrait (AST) en utilisant le module standard ast ;
2. Vérification des contraintes de Pythran (pas d'eval !) et simplification de l'AST en une Représentation Interne (IR) qui reste proche de l'AST;
3. Analyses de l'IR (par exemple fonctions pures, use-def chains, cfg, etc.) et optimisations (indépendantes du type) de l'IR, par exemple propagation de constantes interprocédurales, élimination de code mort, forward substitution, passage en évaluation paresseuse
4. Génération de code, soit du Python, dans ce cas Pythran n'a pas besoin des annotations d'export de fonction et on s'arrête à cette étape, soit du C++, dans ce cas un Méta-programme (comprendre : une soupe de template) est généré et on passe à l'étape suivante.
5. Instanciation du méta-programme à l'aide des annotations d'export de fonction par le compilateur C++ de votre choix. Cette étape repose sur la présence d'une bibliothèque, pythonic basée sur NT2, qui fournit des conteneurs et des opérations optimisées pour la calcul scientifique. Certaines optimisations ont lieu à ce moment (Expression templates mes amours).

Performance du code généré

Il est toujours difficile (et bien souvent biaisé) de sortir des courbes de performance et de se comparer aux autres. Prenons le cas de la fonction harris. Elle vient de la suite de validation de parakeet, donc au moins, elle n'est pas biaisée en notre faveur ;-)

Pour mesurer les perfs, j'utilise le très standard module timeit. Bon en fait je l'ai monkey patch pour qu'il me donne aussi l'écart type, mais ça a peu d'importance. Donc :

~~&xx000A;$~~ 
python -m timeit -s 'from harris import harris; import numpy as np ; a = np.random.randn(512, 512)' 'harris(a)'

Pour du code Python (2.7, désolé, numpy 1.8.2), j'obtiens ça (les temps sont en nanosecondes, c'est ce que sort timeit):

bench  engine min  average dev
-----  ------ ---  ------- ---
harris Python 5441 5471    25

En compilant avec Pythran, sans options particulières et en utilisant GCC (4.9) comme backend, on a grosso modo un facteur d'accélération de deux, qui correspond principalement à la suppression des tableaux intermédiaires :

bench  engine  min  average dev
-----  ------  ---  ------- ---
harris pythran 2512 2644    150

Cela correspond à peu près aux temps obtenus par parakeet (0.24):

bench  engine   min  average dev
-----  ------   ---  ------- ---
harris parakeet 2919 2939    18

numba (0.15.1) est hors course, mais le code de haut niveau, ce n'est pas trop son truc, il préfère les boucles :

bench  engine min  average dev
-----  ------ ---  ------- ---
harris numba  6420 6435    10

On peut s'amuser à utiliser Clang (3.5) au lieu de GCC (4.9) comme backend, dans ce cas on obtient:

bench  engine        min  average dev
-----  ------        ---  ------- ---
harris pythran-clang 2563 2695    144

c'est-à-dire pas de grosse différence.

Pas de benchmark Cython, parce qu'utiliser Cython, ça veut dire réécrire son code, que c'est berk et que je n'ai pas envie de benchmarker ma (mé)connaissance de leur sous-langage.

Pour faciliter nos développement, j'ai regroupé des benchmarks glanés ici et là sur un dépôt https://github.com/serge-sans-paille/numpy-benchmarks avec deux trois outils pour collecter les temps et les formater. Je ne vais pas me lancer dans l'exercice périlleux de leur interprétation, mais voilà une sortie qui donnera peut-être lieu à des commentaires intéressants !

Les temps sont ceux sortis par timeit en nanosecondes (c'est des moyennes hein), 0 veut juste dire que le compilateur a échoué à compiler le code. Ni vectorisation ni parallélisation là dedans ;-)

                          Python   numba parakeet pythran pythran-clang
      allpairs_distances   37601       0     1752  *1718*          1763
allpairs_distances_loops   50650   55759     2611    1757        *1671*
            arc_distance    1377    1373     1444     855         *816*
                    conv 1947702 2596320   *1805*    1868          1886
             create_grid  *3788*    3895        0    3836          4059
                cronbach    1702    1642        0  *1479*          1782
               diffusion   22597   23201    14684    4913        *4153*
                  evolve    5958       0        0    3887        *3408*
                    fdtd 1571130 2028434        0  *1184*          2802
                     fft   24023       0        0   *813*          1138
                grouping    2117       0        0     836         *785*
                 growcut 1799480    4123   *2114*    3861          5545
                  harris    5471    6435     2939  *2644*          2695
                 hasting       9      10       61     *1*           *1*
                 hyantes  258354  294639     1863  *1707*          1920
                   julia 2678910    2801     *49*    2465          2498
                  l2norm    5807    5753    12221   *876*           881
            local_maxima   57142       0        0    5198        *1299*
                  lstsqr    7544    6851        0    2305        *2282*
                  mandel  443242  *4581*        0    5153          5758
            multiple_sum    2999    3132     7640    1673        *1254*
                pairwise 3937063    5665   *3397*    3471          3424
           periodic_dist    1779    1870        0  *1002*          1136
               repeating    1245    1282        0     659         *560*
          reverse_cumsum    2636    2642        0  *2370*          2372
                   rosen   14192   12891     1522    1446        *1232*
               slowparts    5606       0    *983*    2470          2507
               smoothing  982609    4878     6266    4866        *3146*
         specialconvolve    8171    7394     7917    2197        *1373*
             vibr_energy    2584    2557     2282  *1192*          1494
                    wave   52355   51844        0  *1078*          1239
                   wdist   82199   91824     2456    1536        *1364*

Nouveauté dans la version 0.6

Pour les courageux, d'après le Changelog, que je vous traduis pour l'occasion :

• Support complet de la vectorisation. Toutes les expressions numpy sont vectorisées (AVX/SSE) grâce à Boost.SIMD !
• Meilleure gestion de la mémoire à la frontière entre Python et C++
• Support des appels par paramètres nommés (e.g. np.zeros(10, dtype=np.int32))
• Meilleur support des nombres complexes
• Beaucoup de nettoyage de code
• Génération de code plus efficace pour les boucles explicites
• Paquet ArchLinux, guide d'installation sur MacOS
• Validation par Travis, qui valide en mode séquentiel, vectoriel et parallèle !
• Améliorations des performances pour les expressions Numpy
• Utilisation de memcpy pour les copies de tableaux / sous-tableaux dès que possible
• Propagation de constante bien plus efficace
• Possibilité de compiler du code Pythran à travers distutils
• Plus de fonctions numpy supportées, et amélioration de l'existant
• Forward substitution de meilleur qualité
• Passage à la nouvelle version de NT2
• Dépendance sur libgomp maintenant optionnelle

La suite ?

Il commence à y avoir une petite communauté d'utilisateurs, c'est très agréable d'avoir des retours (même si ce sont souvent des retours de bugs), voir des demandes d'améliorations ! Passez faire un coucou sur l'IRC (FreeNode, #pythran) si ça vous tente de contribuer ;-)

Aller plus loin

• Site du projet (206 clics)
• Code source (49 clics)
• Chez le fromager (67 clics)
• Présentation récente (48 clics)
• Le tag pythran sur LinuxFR (53 clics)