URL: https://linuxfr.org/users/serge_ss_paille/journaux/pythran-chatouille-cython Title: Pythran chatouille Cython Authors: serge_sans_paille Date: 2017年01月04日T21:49:31+01:00 License: CC By-SA Tags: pythran, cython et python Score: 32 Bonsoir'nal Peux-être te souviens-tu de [Pythran](http://pythonhosted.org/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](https://linuxfr.org/tags/pythran/public) - [le blog technique du projet](http://serge-sans-paille.github.io/pythran-stories/) 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](http://stackoverflow.com/questions/41314801/can-this-cython-code-be-optimized) 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 : ```python 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](http://www.numpy.org/) et des accès directs de tableau dans une boucle ``while``. les premiers sont hors d'atteinte de [cython](http://cython.org/), 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 : ```python 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 = malloc(n*sizeof(int)) double *d = malloc(n*sizeof(double)) double *v = 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 fonction``sum_axis`` par exemple) mais ça reste proche du code d'origine. Une fois compilé avec cython et transformé en un module natif: ```shell> 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](https://docs.python.org/2/library/timeit.html), j'ai joué avec [perf](http://perf.readthedocs.io/en/latest/) développé par le sympatique haypo. ```shell> 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 : ```python #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](http://serge-sans-paille.github.io/pythran-stories/from-pythran-import-typing.html) qui lui permet de détecter certaines (mais pas toutes) erreurs. Par exemple si j'avais écrit : ```python #pythran export npperm(float) ``` J'aurais eu le message suivant : ```shell 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](https://www.python.org/dev/peps/pep-0484/) 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 : ```shell> pythran perm.py ``` Puis on lance le bench, nos mains moites d'anxiété : ```shell> 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 ```shell> pythran perm.py -DUSE_BOOST_SIMD -march=native ``` Et là ```shell> 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](https://github.com/NumScale/boost.simd) et pour [perf](http://perf.readthedocs.io/en/latest/), c'est tout bon !

AltStyle によって変換されたページ (->オリジナル) /