Journal Pythran chatouille Cython

Posté par  (site web personnel) . Licence CC By‑SA.
Étiquettes :
32
4
jan.
2017

Bonsoir'nal

Peux-être te souviens-tu de Pythran, ce compilateur pour Python spécialisé pour le calcul scientifique ? Si ce n'est pas le cas Alzheimer te guette, mais une médicamentation est toujours possible :
- le tag pythran sur LinuxFr
- le blog technique du projet
Quoiqu'il en soit, j'ai eu l'envie soudaine, inopinée et irrésistible de me re-balader sur |StackOverflow](http://stackoverflow.com/) à le recherche d'un cas test à optimiser. Et là je tombe (ouille) sur ce post qui explore l'optimisation en Cython. Et c'est parti pour une comparaison que j'espère pas trop biaisée.

Le code en question ressemble à ça :

import numpy as np
def npperm(M):
    n = M.shape[0]
    d = np.ones(n)
    j =  0
    s = 1
    f = np.arange(n)
    v = M.sum(axis=0)
    p = np.prod(v)
    while (j < n-1):
        v -= 2*d[j]*M[j]
        d[j] = -d[j]
        s = -s
        prod = np.prod(v)
        p += s*prod
        f[0] = 0
        f[j] = f[j+1]
        f[j+1] = j+1
        j = f[0]
    return p/2**(n-1)

Y a des appels au package de calcul scientifique numpy et des accès directs de tableau dans une boucle while. les premiers sont hors d'atteinte de cython, mais les deuxièmes --- si on annote le code efficacement --- ouvres de belles perspectives d'optimisation : pour faire simple l'utilisateur va annoter toutes les déclarations de variables pour permettre à cython de transformer le code Python en code C faisant le moins d'appels possible à l'API C de CPython. Et le meilleur code trouvé par la communauté donne ça :

import numpy as np
cimport numpy as np
cimport cython
from libc.stdlib cimport malloc, free
from libc.math cimport pow

cdef inline double sum_axis(double *v, double *M, int n):
    cdef:
        int i, j
    for i in range(n):
        for j in range(n):
            v[i] += M[j*n+i]


@cython.boundscheck(False) 
@cython.wraparound(False)
def permfunc_modified(np.ndarray [double, ndim =2, mode='c'] M):
    cdef:
        int n = M.shape[0], j=0, s=1, i
        int *f = <int*>malloc(n*sizeof(int))
        double *d = <double*>malloc(n*sizeof(double))
        double *v = <double*>malloc(n*sizeof(double))
        double p = 1, prod

    sum_axis(v,&M[0,0],n)

    for i in range(n):
        p *= v[i]
        f[i] = i
        d[i] = 1

    while (j < n-1):
        for i in range(n):
            v[i] -= 2.*d[j]*M[j, i]
        d[j] = -d[j]
        s = -s
        prod = 1
        for i in range(n):
            prod *= v[i]
        p += s*prod
        f[0] = 0
        f[j] = f[j+1]
        f[j+1] = j+1
        j = f[0]

    free(d)
    free(f)
    free(v)
    return p/pow(2.,(n-1))

C'est (un peu) verbeux et certains bouts de code ont été recodés à la main (la fonctionsum_axis par exemple) mais ça reste proche du code d'origine.

Une fois compilé avec cython et transformé en un module natif:

> cython perm.pyx
> gcc -O3 -march=native perm.c -shared -fPIC `python-config --cflags --libs` -o perm.so

On peut lancer quelques benchmarks. Et au lieu d'utiliser le vénérable module timeit, j'ai joué avec perf développé par le sympatique haypo.

> python -m perf timeit -s 'from scipy.stats import ortho_group; from perm import permfunc_modified as npperm ; import numpy as np; np.random.seed(7);M = ortho_group.rvs(23)' 'npperm(M)'
.....................
Median +- std dev: 106.1 ms +- 1 ms

Sympa ! Par rapport à la sortie de timeit j'ai des infos sur la déviation standard (et là ça va mon bench a l'air stable) et ça me sort la médiane au lieu du minimum.

On a donc notre score à battre : 106.1ms.

Passons donc à Pythran. Là c'est facile, on a juste à ajouter un commentaire sur le fichier initial :

#pythran export npperm(float[:,:])
import numpy as np
def npperm(M):
    n = M.shape[0]
    d = np.ones(n)
    j =  0
    s = 1
    f = np.arange(n)
    v = M.sum(axis=0)
    p = np.prod(v)
    while (j < n-1):
        v -= 2*d[j]*M[j]
        d[j] = -d[j]
        s = -s
        prod = np.prod(v)
        p += s*prod
        f[0] = 0
        f[j] = f[j+1]
        f[j+1] = j+1
        j = f[0]
    return p/2**(n-1)

À noter pour les habitués que depuis peu, pythran s'est vu doté d'un vérificateur de type qui lui permet de détecter certaines (mais pas toutes) erreurs. Par exemple si j'avais écrit :

#pythran export npperm(float)

J'aurais eu le message suivant :

CRITICAL I am in trouble. Your input file does not seem to match Pythran's constraints...
E: Specification for `npperm` does not match inferred type:
expected `Callable[[Array[1 d+, T0]], Array[0 d+, T1]]`
got `Callable[[float], ...]`

Comprenne qui voudra, il vaut mieux avoir lu la PEP484 pour piger les notations. Mais en gros Pythran a estimé que le type de npperm c'était fonction qui prend un tableau au moins 1D et renvoie un tableau au moins 1D et il est pas content car on lui file un nombre flottant.

Donc on compile :

> pythran perm.py

Puis on lance le bench, nos mains moites d'anxiété :

> python -m perf timeit -s 'from scipy.stats import ortho_group; from permpy import npperm ; import numpy as np; np.random.seed(7);M = ortho_group.rvs(23)' 'npperm(M)'
.....................
Median +- std dev: 173 ms +- 6 ms

Moarf, c'est moisi :-/ Ah mais on a oublié d'utiliser la vectorisation, tsss. On recommence

> pythran perm.py -DUSE_BOOST_SIMD -march=native

Et là

> python -m perf timeit -s 'from scipy.stats import ortho_group; from permpy import npperm ; import numpy as np; np.random.seed(7);M = ortho_group.rvs(23)' 'npperm(M)'
.....................
Median +- std dev: 93.7 ms +- 3.9 ms

On est dans les mêmes ordres de grandeur, on a fait moins d'effort, on a honteusement fait de la pub pour boost.simd et pour perf, c'est tout bon !

  • # Temps d'exécution sans optimisation ?

    Posté par  (site web personnel) . Évalué à 10.

    J'ai trouvé les temps d'exécution pour la version Cython, pour la version Pythran sans SIMD et la version Pythran avec SIMD. Hum, mais quel est le point de départ ? Quelles sont les perfs sans Cython ni Pythran, juste avec numpy avec CPython standard ?

    • [^] # Re: Temps d'exécution sans optimisation ?

      Posté par  (site web personnel) . Évalué à 8.

      Et qu'en est-il des perfs en C, C++, ou fortran ? Est-ce que cette recherche d'optimisation en python est bien raisonnable ? Je pose honnêtement la question. Ce langage n'a-t-il pas pour objectif de privilégier la simplicité d'écriture sur la vitesse ? Du coup ces efforts ne sont-ils pas quelque peu contre nature ? Le code optimisé en Cython paraît s'être tellement compliqué que le béotien que je suis se demande bien s'il n'aurait pas été plus efficace de l'écrire directement dans un langage plus bas niveau ? Si la version Pythran paraît plus naturelle, il reste tout de même délicat de comprendre la finalité de ces efforts. Est-ce que ça ne revient pas à transformer du code interpréter en code compilé, en perdant au passage l'essentiel des avantages du Python ?

      « IRAFURORBREVISESTANIMUMREGEQUINISIPARETIMPERAT » — Odes — Horace

      • [^] # Re: Temps d'exécution sans optimisation ?

        Posté par  . Évalué à 10.

        Si je ne m'abuse, le passage du code Python naïf à la version optimisée Pythran n'a nécessité qu'une ligne, donc je trouve que c'est plus que raisonnable : on conserve l'avantage du langage Python pour un prototype vite fait et on a des performances plus honnêtes (voire bonnes?) en réfléchissant juste sur une ligne qui permet à Pythran de traiter ce code!

      • [^] # Re: Temps d'exécution sans optimisation ?

        Posté par  (site web personnel) . Évalué à 4.

        Il y a un état de fai[tl] : la demande est là. Il suffit de regarder les compilos existants :

        • cython
        • numba
        • pyston
        • pypy

        Et chacun de ces compilos a pour objectif principal ou auxiliaire mais en tout cas affiché d'optimiser du code Python utilisant numpy :-)

        • [^] # Re: Temps d'exécution sans optimisation ?

          Posté par  (site web personnel) . Évalué à 2.

          Certes. Mais en ce qui me concerne, ces questions se posent dans une autre perspective ; celle d'un choix didactique. Dans le cadre d'une initiation au calcul scientifique, vaut-il mieux enseigner aux étudiants du Python ou du C/C⁺⁺ ? D'où l'intérêt pour ces compilateurs de python qui changent certains éléments dans la balance.

          « IRAFURORBREVISESTANIMUMREGEQUINISIPARETIMPERAT » — Odes — Horace

          • [^] # Re: Temps d'exécution sans optimisation ?

            Posté par  (site web personnel) . Évalué à 10.

            L'idée c'est que tu n'as pas besoin de tout optimiser. Python est beaucoup plus pratique que C++ pour toute la glu qui entoure ton code scientifique: manipulation de fichiers, gestion des données et des tests, visualisation, etc.

            Donc en ce qui me concerne, pouvoir tout coder en Python et garder des performances acceptables me semble très intéressant. Regarde le projet scikit-learn par exemple: l'utilisateur n'a que quelques lignes de python idiomatique à écrire pour exécuter des algorithmes complexes de machine learning.

          • [^] # Re: Temps d'exécution sans optimisation ?

            Posté par  . Évalué à 4.

            Dans le cadre d'une initiation, je dirais que les performances pures importent probablement moins que le fait de pouvoir prototyper et tester rapidement des algorithmes non-triviaux.

            Cela dépend aussi du niveau de tes étudiants en programmation : s'ils sont déjà à leur aise en C ou en C++, cela vaut peut-être le coup de continuer avec ces langages.

      • [^] # Re: Temps d'exécution sans optimisation ?

        Posté par  . Évalué à 5.

        Pour te répondre, tu gardes la mise en donnée de python, tu manipules un langage que tu connais, donc même si cela à l'air plus compliqué pour celui qui maîtrise le C++, ça l'est beaucoup moins que d'apprendre un nouveau langage, maîtriser la chaîne de compilation et savoir comment transférer les différents structures/objets de C/C++ vers python. Tout cela se fait bien avec la lib python.h mais cela demande de fortes compétences. Donc si tu veux juste accélérer tes boucles et conditions car c'est principalement de cela dont il s'agit, alors tu es bien content, même s'il faut jouer les contorsionnistes d'avoir une solution simple telle que cython/pythran qui te permet de tout créer directement depuis le langage que tu connais même s'il a de nouveaux concepts à apprendre, il sont accessibles car tu peux toujours simplement utiliser les libs que tu connais, donc en terme ingénierie logicielle cela n'a rien à voir. Je trouve que ces démarches au contraire sont vraiment salutaires.

      • [^] # Re: Temps d'exécution sans optimisation ?

        Posté par  . Évalué à 3.

        Tu pars du principe que le programmeur est meilleur que la machine pour optimiser et donc qu'il est plus efficace de partir d'un langage bas niveau qui donne toute la latitude pour piloter le calcul.

        Mais c'est généralement faux. Il est souvent préférable d'avoir un langage haut niveau plus expressif, qui exprime l'intention algorithmique et de laisser un compilateur faire des optimisations difficilement appréhendable par un programmeur.

        C'est ce que tente de faire pythran, cython ou pypy avec Python : permettre, par des annotations, de décrire l'intention algorithmique (et aussi d'ajouter du typage, moyen le plus efficace de permettre des optimisations et qui manque cruellement à python, ne nous leurrons pas).

        Les langages fonctionnelles sont très efficaces pour cela. Les notions de fonctions pures, de typage fort (avec ou sans inférence), d'évaluation paresseuse, les map (qui sont une forme de SIMD), donnent beaucoup d'information à un compilateur pour faire des optimisations très efficaces.

        • [^] # Re: Temps d'exécution sans optimisation ?

          Posté par  . Évalué à 5.

          Tu pars du principe que le programmeur est meilleur que la machine pour optimiser et donc qu'il est plus efficace de partir d'un langage bas niveau qui donne toute la latitude pour piloter le calcul.

          Ça dépend de quel type d'optimisation on discute (voir plus bas). Dans le cas normal je suis d'accord avec toi, le compilateur est souvent meilleur que toi. Lorsqu'on parle de gens qui cherchent à optimiser une boucle après avoir observé que c'est celle-ci prend beaucoup de temps à s'exécuter, alors c'est un cas différent : tu as déjà compilé avec -O3, possiblement en promettant à ton compilo que les pointeurs ne pointent pas sur des zones qui se recouvrent1

          Les langages fonctionnelles sont très efficaces pour cela. Les notions de fonctions pures, de typage fort (avec ou sans inférence), d'évaluation paresseuse, les map (qui sont une forme de SIMD)

          (c'est moi qui graisse) Si le concept est similaire (une seule instruction appliquée à un jeu de données vs. une seule fonction appliquée à une liste d'éléments), les implications en termes de génération de code n'ont strictement rien à voir.

          Quand il s'agit de vectorisation, les compilateurs (y compris Intel) sont quand même pas complètement au point. Ils font souvent de la « vectorisation partielle », dans le sens où ils utilisent correctement les registres SSE/AVX (avec le bon type d'instructions du coup), mais où ils ne savent pas s'ils ont le droit de correctement charger les données 128 ou 256 bits à la fois, et donc utilisent ces registres principalement pour charger des données 64 bits dans des registres qui font le double ou le quadruple en capacité. Depuis peu, GCC fait de l'auto-vectorisation partielle (qui me semble légèrement moins au point que celle d'Intel), mais dans tous les cas, à moins d'avoir des garanties en béton sur les pointeurs, sur le code, etc., la vectorisation complète (càd : qui utilise pleinement les registres vectoriels des machines x86/x64) est rarement achevée. C'est pour ça qu'on continue d'utiliser les intrinsics à la main dans beaucoup de codes à partir du moment où on se dit que ça va pas assez vite. Après, il existe des outils pour déterminer si ça vaut le coup de se prendre la tête à vectoriser, en regardant le rapport entre le nombre d'instructions arithmétiques et le nombre de mouvements mémoire : un code peut être vectorisable sans pour autant faire gagner en performance, car le processeur passe son temps à attendre que la mémoire lui fournisse les prochains éléments à traiter.

          D'ailleurs, c'est tellement reconnu que le dernier standard en date pour OpenMP (qui est le langage de programmation parallèle le plus populaire dans la communauté du calcul intensif) fournit une clause (le mot-clef « simd ») pour dire au compilateur que le code est vectorisable (donc en gros le programmeur prend la responsabilité d'avoir vérifié que l'agencement des données permet effectivement de vectoriser le code, mais au moins il n'a plus à utiliser les intrinsics — enfin, on espère).


          1. En C/C++, ça s'exprime avec restrict dans la signature des fonctions, ou au niveau du compilateur tu peux le spécifier « par fichier » avec un flag de type -fno-alias chez Intel ou -fstrict-aliasing chez GCC je crois — et il semblerait qu'au moins avec gcc ce soit activé à partir de -O2.  

    • [^] # Re: Temps d'exécution sans optimisation ?

      Posté par  (site web personnel) . Évalué à 8.

      CPython pur :

      > rm *.so
      > python -m perf timeit -s 'from scipy.stats import ortho_group; from perm import npperm ; import numpy as np; np.random.seed(7);M = ortho_group.rvs(23)' 'npperm(M)'
      .....................
      Median +- std dev: 20.4 sec +- 1.1 sec
      

      c'est ultra lent :-)

  • # SIMD

    Posté par  . Évalué à 4.

    Bravo pour les performances, c'est impressionnant !
    La version cython utilise-t-elle le SIMD ? (pour savoir si on compare une version SIMD à une version non-SIMD, et le cas échéant, combien de coeurs sont utilisés)

    Pour ma part j'avais considéré l'utilisation de pythran mais j'ai finalement opté pour cython car j'avais vraiment besoin de classes, et de toute façon je ne partais pas d'un code python existant. Mais ça reste très intéressant pour accélérer des fonctions simples sans se prendre la tête.

    • [^] # Re: SIMD

      Posté par  (site web personnel) . Évalué à 3.

      La version cython utilise-t-elle le SIMD ?

      Elle est compilée avec des drapeaux qui permettent la vectorisation, mais après avoir inspecté le binaire, gcc a pas réussi à vectoriser, non. Du coup ça veut dire que la version pythran a encore un soucis, mais je n'ai pas encore trouvé lequel :-)

      car j'avais vraiment besoin de classes

      C'est dans ma TODO liste, tiens tu attendrais quoi comme annotation #pythran export pour une classe ? C'est pas bien fixé pour moi…

      • [^] # Re: SIMD

        Posté par  . Évalué à 3.

        tu attendrais quoi comme annotation #pythran export pour une classe ?

        J'imagine qu'il faudrait annoter chaque méthode individuellement, comme on annote une fonction… (?)
        En cython il est appréciable de pouvoir cythoniser (cdef) certaines méthodes, et pas d'autres. Par exemple dans mon cas, les constructeurs n'ont pas besoin d'être optimisés car appelés une seule fois, ils restent en python, et on peut donc implémenter librement toute la partie lecture de données et initialisation qui font appel à du code "non pythranisable". L'idée est d'optimiser uniquement les méthodes "lourdes" qui concentrent l'essentiel des calculs.

Suivre le flux des commentaires

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