Pythran 0.6 - compilation de noyaux scientifiques écrits en Python

36
6
nov.
2014
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

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 ;-)

  • # Le même en Fortran SVP

    Posté par (page perso) . Évalué à 3.

    Ça serait intéressant d'avoir la même routine en Fortran et la compiler avec gfortran (compilo intégré avec GCC) afin de pouvoir comparer les performances.

  • # Support de Numpy

    Posté par . Évalué à 2.

    Tout d'abord, un joli projet ! J'utilise python pour des petites simulations physiques et jusqu'ici ses performances m'ont suffit, mais vu la facilité d'utilisation de pythran, il faut vraiment que je vois ça de plus près. Surtout qu'il permet de faire du parallélisme en python, alors que j'ai toujours lu que ce n'était pas possible à faire directement parce que GIL. (Note : justement, là, ce n'est pas direct. Mais pour les programmes scientifiques, on s'en contrefiche que ce soit compilé plutôt qu’interprété.)

    Le seul point qui pourrait me bloquer est le support de numpy, random et scipy. Pour numpy, je n'est pas trouvé de liste des fonctions qui marchent (ou, au contraire, qui ne marchent pas), pour random, ça semble dans la TODO liste, mais je n'ai rien trouvé sur scipy.

    Un beau projet, qui me pousse plus encore à jouer avec python plutôt que C !

    • [^] # Re: Support de Numpy

      Posté par (page perso) . Évalué à 4.

      Content que ça te plaise \o/

      La liste des fonctions numpy est présente cough cough

      le support de numpy.random n'est pas problématique, pour scipy, on a une PR en cours pour que pythran compile les modules importés directement, ça devrait faire avancer la chose. mais c'est sur qu'avoir une bonne couverture de numpy/scipy, c'est un long chemin.

      Astuce : quand une fonction te manque, ouvre un bug, envoie un mail sur la mliste ou poke nous sur IRC, ça motive les troupes !

    • [^] # Re: Support de Numpy

      Posté par . Évalué à 1.

      Surtout qu'il permet de faire du parallélisme en python, alors que j'ai toujours lu que ce n'était pas possible à faire directement parce que GIL

      Autant que je sache, le problème du GIL peut être contourné soit en utilisant multiprocessing, soit mpi4py.

      • [^] # Re: Support de Numpy

        Posté par (page perso) . Évalué à 1.

        Certes, mais avec une granularité toute autre.

        • [^] # Re: Support de Numpy

          Posté par . Évalué à 2.

          • multiprocessing ==> fork (resp thread pour multithread) : OK variables pas partagées par défaut (mais assez simple de partager des arrays) mais c'est le même "grain", on bosse sur son multicœurs, voire même entre plusieurs grâce à l'utilisation de multiprocessing.Manager.
          • mpi4py ==> idem, tu bosses au choix sur une ou plusieurs machines multicœurs (sans parler de jouer avec openMP)

          J'ai bon ? Du coup, pas sûr de comprendre : qu'entends-tu par "granularité toute autre" ?

          Au passage, merci pour ce super projet !!

          • [^] # Re: Support de Numpy

            Posté par (page perso) . Évalué à 1.

            De rien ;-) (Faut remercier Pierrick Brunet qui fait un boulot de ouf dessus aussi !)

            Je pensais à la granularité de l'unité de calcul minimale à partir de laquelle il devient intéressant de paralléliser. En gros

            a + b ** 2
            

            où a et b sont des tableaux, ni mpi4py ni multiprocessing ne te donneront de speedup là dessus (à cause du coup de transfert mémoire) alors qu'en pure OpenMP on peut gagner (un peu sur cette exemple, mais pas beaucoup) et avec de l'AVX / SSE on gagne pas mal, suivant le type de a et b.

            Je ne sais pas si mes explications sont claires…

            • [^] # Re: Support de Numpy

              Posté par . Évalué à 2.

              à cause du coup de transfert mémoire

              Tu peux utiliser de la mémoire partagée. Pas de copie dans ce cas. Après pour une seule ligne, ça ne vaut sans doute pas souvent le coup de sortir l'artillerie multiprocessing avec l'overhead que ça comporte.

              • [^] # Re: Support de Numpy

                Posté par (page perso) . Évalué à 1.

                Je comprends bien que multiprocessing permet d'utiliser de la mémoire partagée, à condition d'utiliser Value ou Array (c'est ce que je comprends de la doc) mais faut copier tes données utilisateurs dans ces conteneurs, donc on a bien un transfert mémoire qu'il faut amortir, non?

                • [^] # Re: Support de Numpy

                  Posté par . Évalué à 3.

                  Ça dépend comment on procède. On peut facilement convertir des array numpy en RawArray et vice versa. Selon les besoins on n'a pas forcément à transférer grand chose.

                  Par exemple :

                      size = 42
                      MyRawArray1 = RawArray(ctypes.c_double, size*1000)
                      MyRawArray2 = RawArray(ctypes.c_double, size*1000)
                  
                      initargs = (MyRawArray1, MyRawArray2)
                      with mp.Pool(processes=4, initializer=initializer, initargs=initargs) as pool:
                          pool.map(worker, range(size))

                  Avec

                      def initializer(arg1, arg2):
                          global MyArray1, MyArray2
                  
                          MyArray1 = np.ctypeslib.as_array(arg1)
                          MyArray2 = np.ctypeslib.as_array(arg2)

                  et

                      def worker(arg):
                          MyArray1[arg*1000:(arg+1)*1000] = DoSomething(arg, 1)
                          MyArray2[arg*1000:(arg+1)*1000] = DoSomething(arg, 2)

                  De cette façon il n'y a pratiquement rien à transféré sauf le pointeur pour que chaque sous-processus crée son tableau numpy qui utilise la mémoire partagée.

  • # Parakeet winner?

    Posté par . Évalué à 2.

    Je ne vais pas me lancer dans l'exercice périlleux de leur interprétation

    En même temps c’est tout l’intérêt de faire des mesures. En voilà une :

    Partout où on voudrait utiliser un accélérateur tel que parakeet, numba ou pythran, c’est à dire quand les temps de calculs sont lents en CPython, parakeet semble être la solution. Si on ne considère que les temps pour lesquels CPython prend plus de 1e6, parakeet est systématiquement meilleur que pythran. Après il faudrait savoir si c’est signifiant, peut-être en regardant si la distribution est normale et en ayant les écart-types, sinon ça ne veut rien dire (déjà que les micro-benchs on sait ce que ça vaut).

    Après, je rejoint l’idée de la comparaison avec du fortran ou l’une des bibliothèques équivalentes en C++.

    • [^] # Re: Parakeet winner?

      Posté par (page perso) . Évalué à 2.

      si on ne considère que les temps pour lesquels CPython prend plus de 1e6, parakeet est systématiquement meilleur que pythran.

      Je peux même pousser l'analyse un peu plus loin : parakeet est meilleur que Pythran pour la gestion des boucles explicites. Les boucles explicites sont plutôt découragées en Numpy, car pas terrible niveau perf - c'est pour cela que tu identifies ces cas comme les cas où CPython est très lent. C'est un problème connu, sur lequel je dois m'atteler.

      En fait parakeet est un très bon outil, plus facile à déployer que Pythran et avec souvent de très bonne performances. Les seuls limitations queje lui donne, c'est qu'il supporte peu de fonctions Numpy (regarde tous les 0 dans sa colonne) et qu'il ne vectorise pas les opérations - mais j'ai fais exprès de ne pas mettre ça en avant dans mes mesures.

      Un concurrent sérieux donc ;-)

Suivre le flux des commentaires

Note : les commentaires appartiennent à ceux qui les ont postés. Nous n'en sommes pas responsables.