1A.soft - Calcul numérique et Cython

Links: notebook, html, PDF, python, slides, GitHub

Python est très lent. Il est possible d’écrire certains parties en C mais le dialogue entre les deux langages est fastidieux. Cython propose un mélange de C et Python qui accélère la conception.

from jyquickhelper import add_notebook_menu
add_notebook_menu()

Calcul numérique

On peut mesurer le temps que met en programme comme ceci (qui ne marche qu’avec IPython…timeit) :

def racine_carree1(x) :
    return x**0.5

%timeit -r 10 [ racine_carree1(x) for x in range(0,1000) ]
537 µs ± 86.5 µs per loop (mean ± std. dev. of 10 runs, 1000 loops each)
import math
def racine_carree2(x) :
    return math.sqrt(x)

%timeit -r 10 [ racine_carree2(x) for x in range(0,1000) ]
442 µs ± 41.7 µs per loop (mean ± std. dev. of 10 runs, 1000 loops each)

La seconde fonction est plus rapide. Seconde vérification :

%timeit -r 10 [ x**0.5 for x in range(0,1000) ]
%timeit -r 10 [ math.sqrt(x) for x in range(0,1000) ]
379 µs ± 36.7 µs per loop (mean ± std. dev. of 10 runs, 1000 loops each)
359 µs ± 88.2 µs per loop (mean ± std. dev. of 10 runs, 1000 loops each)

On remarque également que l’appel à une fonction pour ensuite effectuer le calcul a coûté environ 100 \mu s pour 1000 appels. L’instruction timeit effectue 10 boucles qui calcule 1000 fois une racine carrée.

Cython

Le module Cython est une façon d’accélérer les calculs en insérant dans un programme python du code écrit dans une syntaxe proche de celle du C. Il existe différentes approches pour accélérer un programme python :

  • Cython : on insère du code [C](http://fr.wikipedia.org/wiki/C_(langage) dans le programme python, on peut gagné un facteur 10 sur des fonctions qui utilisent des boucles de façon intensives.

  • autres alternatives :

  • PyPy : on compile le programme python de façon statique au lieu de l’interpréter au fur et à mesure de l’exécution, cette solution n’est praticable que si on a déjà programmé avec un langage compilé ou plus exactement un langage où le typage est fort. Le langage python, parce qu’il autorise une variable à changer de type peut créer des problèmes d’inférence de type.

  • module implémenté en C : c’est le cas le plus fréquent et une des raisons pour lesquelles Python a été rapidement adopté. Beaucoup de librairies se sont ainsi retrouvées disponibles en Python. Néanmoins, l’API C du Python nécessite un investissement conséquent pour éviter les erreurs. Il est préférable de passer par des outils tels que

    • boost python : facile d’accès, le module sera disponible sous forme compilée,

    • SWIG : un peu plus difficile, le module sera soit compilé par la librairie soit packagé de telle sorte qu’il soit compilé lors de son l’installation.

Parmi les trois solutions, la première est la plus accessible, et en développement constant (Cython changes).

L’exemple qui suit ne peut pas fonctionner directement sous notebook car Cython compile un module (fichier *.pyd) avant de l’utiliser. Si la compilation ne fonctionne pas et fait apparaître un message avec unable for find file vcvarsall.bat, il vous faut lire l’article Build a Python 64 bit extension on Windows 8 après avoir noté la version de Visual Studio que vous utilisez. Il est préférable d’avoir programmé en C/C++ même si ce n’est pas indispensable.

Cython dans un notebook

Le module IPython propose une façon simplifiée de se servir de Cython illustrée ici : Some Linear Algebra with Cython. Vous trouverez plus bas la façon de faire sans IPython que nous n’utiliserons pas pour cette séance. On commence par les préliminaires à n’exécuter d’une fois :

%load_ext cython

Puis on décrit la fonction avec la syntaxe Cython :

%%cython --annotate
cimport cython

def cprimes(int kmax):
    cdef int n, k, i
    cdef int p[1000]
    result = []
    if kmax > 1000:
        kmax = 1000
    k = 0
    n = 2
    while k < kmax:
        i = 0
        while i < k and n % p[i] != 0:
            i = i + 1
        if i == k:
            p[k] = n
            k = k + 1
            result.append(n)
        n = n + 1
    return result
Cython: _cython_magic_8fa17b4b01e0386d07eadd55c699ff6a.pyx

Generated by Cython 0.29.21

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

 01: cimport cython
 02: 
+03: def cprimes(int kmax):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_8fa17b4b01e0386d07eadd55c699ff6a_1cprimes(PyObject *__pyx_self, PyObject *__pyx_arg_kmax); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_8fa17b4b01e0386d07eadd55c699ff6a_1cprimes = {"cprimes", (PyCFunction)__pyx_pw_46_cython_magic_8fa17b4b01e0386d07eadd55c699ff6a_1cprimes, METH_O, 0};
static PyObject *__pyx_pw_46_cython_magic_8fa17b4b01e0386d07eadd55c699ff6a_1cprimes(PyObject *__pyx_self, PyObject *__pyx_arg_kmax) {
  int __pyx_v_kmax;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("cprimes (wrapper)", 0);
  assert(__pyx_arg_kmax); {
    __pyx_v_kmax = __Pyx_PyInt_As_int(__pyx_arg_kmax); if (unlikely((__pyx_v_kmax == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_8fa17b4b01e0386d07eadd55c699ff6a.cprimes", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_8fa17b4b01e0386d07eadd55c699ff6a_cprimes(__pyx_self, ((int)__pyx_v_kmax));
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_8fa17b4b01e0386d07eadd55c699ff6a_cprimes(CYTHON_UNUSED PyObject *__pyx_self, int __pyx_v_kmax) {
  int __pyx_v_n;
  int __pyx_v_k;
  int __pyx_v_i;
  int __pyx_v_p[0x3E8];
  PyObject *__pyx_v_result = NULL;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("cprimes", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_AddTraceback("_cython_magic_8fa17b4b01e0386d07eadd55c699ff6a.cprimes", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XDECREF(__pyx_v_result);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple_ = PyTuple_Pack(7, __pyx_n_s_kmax, __pyx_n_s_kmax, __pyx_n_s_n, __pyx_n_s_k, __pyx_n_s_i, __pyx_n_s_p, __pyx_n_s_result); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple_);
  __Pyx_GIVEREF(__pyx_tuple_);
/* … */
  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_8fa17b4b01e0386d07eadd55c699ff6a_1cprimes, NULL, __pyx_n_s_cython_magic_8fa17b4b01e0386d07); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_cprimes, __pyx_t_1) < 0) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 04:     cdef int n, k, i
 05:     cdef int p[1000]
+06:     result = []
  __pyx_t_1 = PyList_New(0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 6, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_v_result = ((PyObject*)__pyx_t_1);
  __pyx_t_1 = 0;
+07:     if kmax > 1000:
  __pyx_t_2 = ((__pyx_v_kmax > 0x3E8) != 0);
  if (__pyx_t_2) {
/* … */
  }
+08:         kmax = 1000
    __pyx_v_kmax = 0x3E8;
+09:     k = 0
  __pyx_v_k = 0;
+10:     n = 2
  __pyx_v_n = 2;
+11:     while k < kmax:
  while (1) {
    __pyx_t_2 = ((__pyx_v_k < __pyx_v_kmax) != 0);
    if (!__pyx_t_2) break;
+12:         i = 0
    __pyx_v_i = 0;
+13:         while i < k and n % p[i] != 0:
    while (1) {
      __pyx_t_3 = ((__pyx_v_i < __pyx_v_k) != 0);
      if (__pyx_t_3) {
      } else {
        __pyx_t_2 = __pyx_t_3;
        goto __pyx_L8_bool_binop_done;
      }
      if (unlikely((__pyx_v_p[__pyx_v_i]) == 0)) {
        PyErr_SetString(PyExc_ZeroDivisionError, "integer division or modulo by zero");
        __PYX_ERR(0, 13, __pyx_L1_error)
      }
      __pyx_t_3 = ((__Pyx_mod_int(__pyx_v_n, (__pyx_v_p[__pyx_v_i])) != 0) != 0);
      __pyx_t_2 = __pyx_t_3;
      __pyx_L8_bool_binop_done:;
      if (!__pyx_t_2) break;
+14:             i = i + 1
      __pyx_v_i = (__pyx_v_i + 1);
    }
+15:         if i == k:
    __pyx_t_2 = ((__pyx_v_i == __pyx_v_k) != 0);
    if (__pyx_t_2) {
/* … */
    }
+16:             p[k] = n
      (__pyx_v_p[__pyx_v_k]) = __pyx_v_n;
+17:             k = k + 1
      __pyx_v_k = (__pyx_v_k + 1);
+18:             result.append(n)
      __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 18, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_1);
      __pyx_t_4 = __Pyx_PyList_Append(__pyx_v_result, __pyx_t_1); if (unlikely(__pyx_t_4 == ((int)-1))) __PYX_ERR(0, 18, __pyx_L1_error)
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+19:         n = n + 1
    __pyx_v_n = (__pyx_v_n + 1);
  }
+20:     return result
  __Pyx_XDECREF(__pyx_r);
  __Pyx_INCREF(__pyx_v_result);
  __pyx_r = __pyx_v_result;
  goto __pyx_L0;

On termine en estimant son temps d’exécution. Il faut noter aussi que ce code ne peut pas être déplacé dans la section précédente qui doit être entièrement écrite en cython.

%timeit [ cprimes (567) for i in range(10) ]
6 ms ± 201 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Exercice : python/C appliqué à une distance d’édition

La distance de Levenshtein aussi appelé distance d’édition calcule une distance entre deux séquences d’éléments. Elle s’applique en particulier à deux mots comme illustré par Distance d’édition et programmation dynamique. L’objectif est de modifier la fonction suivante de façon à utiliser Cython puis de comparer les temps d’exécution.

def distance_edition(mot1, mot2):
    dist = { (-1,-1): 0 }
    for i,c in enumerate(mot1) :
        dist[i,-1] = dist[i-1,-1] + 1
        dist[-1,i] = dist[-1,i-1] + 1
        for j,d in enumerate(mot2) :
            opt = [ ]
            if (i-1,j) in dist :
                x = dist[i-1,j] + 1
                opt.append(x)
            if (i,j-1) in dist :
                x = dist[i,j-1] + 1
                opt.append(x)
            if (i-1,j-1) in dist :
                x = dist[i-1,j-1] + (1 if c != d else 0)
                opt.append(x)
            dist[i,j] = min(opt)
    return dist[len(mot1)-1,len(mot2)-1]

%timeit distance_edition("idstzance","distances")
145 µs ± 19.9 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Auparavant, il est probablement nécessaire de suivre ces indications :

  • Si vous souhaitez remplacer le dictionnaire par un tableau à deux dimensions, comme le langage C n’autorise pas la création de tableau de longueur variables, il faut allouer un pointeur (c’est du C par du C++). Toutefois, je déconseille cette solution :

Je suggère donc de remplacer dist par un tableau cdef int dist [500][500]. La signature de la fonction est la suivante : def cdistance_edition(str mot1, str mot2). Enfin, Cython a été optimisé pour une utilisation conjointe avec numpy, à chaque fois que vous avez le choix, il vaut mieux utiliser les container numpy plutôt que d’allouer de grands tableaux sur la pile des fonctions ou d’allouer soi-même ses propres pointeurs.

Cython sans les notebooks

Cette partie n’est utile que si vous avez l’intention d’utiliser Cython sans IPython. Les lignes suivantes implémentent toujours avec Cython la fonction primes qui retourne les entiers premiers entiers compris entre 1 et N. On suit maintenant la méthode préconisée dans le tutoriel de Cython. Il faut d’abord créer deux fichiers :

  • example_cython.pyx qui contient le code de la fonction

  • setup.py qui compile le module avec le compilateur Visual Studio C++

code = """
def primes(int kmax):
    cdef int n, k, i
    cdef int p[1000]
    result = []
    if kmax > 1000:
        kmax = 1000
    k = 0
    n = 2
    while k < kmax:
        i = 0
        while i < k and n % p[i] != 0:
            i = i + 1
        if i == k:
            p[k] = n
            k = k + 1
            result.append(n)
        n = n + 1
    return result
"""

name = "example_cython"
with open(name + ".pyx","w") as f : f.write(code)

setup_code = """
from distutils.core import setup
from Cython.Build import cythonize

setup(
    ext_modules = cythonize("__NAME__.pyx",
                            compiler_directives={'language_level' : "3"})
)
""".replace("__NAME__",name)

with open("setup.py","w") as f:
    f.write(setup_code)

Puis on compile le fichier .pyx créé en exécutant le fichier setup.py avec des paramètres précis :

import os
import sys
cmd = "{0} setup.py build_ext --inplace".format(sys.executable)
from pyquickhelper.loghelper import run_cmd
out,err = run_cmd(cmd)
if err != '' and err is not None:
    raise Exception(err)

[ _ for _ in os.listdir(".") if "cython" in _ or "setup.py" in _ ]
['example_cython.c',
 'example_cython.cp38-win_amd64.pyd',
 'example_cython.pyx',
 'setup.py',
 'td1a_cython_edit.ipynb',
 'td1a_cython_edit_correction.ipynb']

Puis on importe le module :

import pyximport
pyximport.install()
import example_cython
c:python387_x64libsite-packagesCythonCompilerMain.py:369: FutureWarning: Cython directive 'language_level' not set, using 2 for now (Py2). This will change in a later release! File: C:xavierdupre__home_GitHubensae_teaching_cs_docnotebookstd1a_softexample_cython.pyx
  tree = Parsing.p_module(s, pxd, full_module_name)

Si votre dernière modification n’apparaît pas, il faut redémarrer le kernel. Lorsque Python importe le module example_cython la première fois, il charge le fichier example_cython.pyd. Lors d’une modification du module, ce fichier est bloqué en lecture et ne peut être modifié. Or cela est nécessaire car le module doit être recompilé. Pour cette raison, il est plus pratique d’implémenter sa fonction dans un éditeur de texte qui n’utilise pas IPython.

On teste le temps mis par la fonction primes :

%timeit [ example_cython.primes (567) for i in range(10) ]
7.47 ms ± 902 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Puis on compare avec la version écrites un Python :

def py_primes(kmax):
    p = [ 0 for _ in range(1000) ]
    result = []
    if kmax > 1000:
        kmax = 1000
    k = 0
    n = 2
    while k < kmax:
        i = 0
        while i < k and n % p[i] != 0:
            i = i + 1
        if i == k:
            p[k] = n
            k = k + 1
            result.append(n)
        n = n + 1
    return result

%timeit [ py_primes (567) for i in range(10) ]
201 ms ± 43 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)