<jupyter_start><jupyter_text>Acelerando Python con Cython La [implementación](https://wiki.python.org/moin/PythonImplementations) más utilizada del interprete/compilador de Python está escrita en C y se llama [CPython](https://github.com/python/cpython). Recordemos que la implementación es aquello que interpreta y corre el código escrito en Python. :::{note}También existen otras implementaciones de Python menos usadas como [PyPy](https://www.pypy.org/), [Jython](https://www.jython.org/) o IronPython. :::CPython compila el código escrito en Python en un código de máquina (binario) de forma transparente al usuario. Luego CPython interpreta el binario. :::{note}Si es de nuestro interés, podemos estudiar dicho código binario usando el módulo [dis](https://docs.python.org/3.7/library/dis.html) de Python::: Por lo tanto, si sabemos C, podemos usar la [API de Python/C](https://docs.python.org/3/c-api/index.html) para 1. Escribir módulos de C que puedan llamarse desde Python1. Hacer interfaces entre código C y PythonSin embargo la API es un poco complicada, por lo que han surgido alternativas menos tediosas para lograr estos objetivos, como por ejemplo [ctypes](https://docs.python.org/2/library/ctypes.html), [CFFI](https://cffi.readthedocs.io/en/latest/), [SWIG](http://swig.org/) y [Cython](https://cython.org/)En este capítulo veremos en detalle esta última<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Cython: C Extensions for PythonRecordemos que **Python** es un lenguaje interpretado con tipos dinámicos. Esto hace que cada operación tenga un overhead. Por ejemplo```python    z = x + y     overhead: Inferir el tipo de x     overhead: Inferir el tipo de y     Hacer la operación suma     overhead: Darle el tipo adecuado a z```en cambio, **Cython** es un lenguaje de programación que le agrega a Python algunas propiedades de C y C++, una de ellas son los **tipos estáticos**```python    int x = 1    int y = 2    int z = x + y     No hay que inferir el tipo de x, y, z```Esto hace que Cython sea menos flexible pero decenas de veces más rápido que Python. En términos de sintaxis Cython es casi tan simple como Python. Sin embargo a diferencia de Python, el lenguaje Cython debe compilarse- El compilador de Cython convierte el código fuente en código C - Luego el código C se compila como un módulo de Python con la implementación **CPython**:::{important}Una vez compilado el código escrito en Cython este puede llamarse desde Python.:::¿Por qué considerar Cython?- Cython es casi tan simple como Python y casi tan rápido como C- Con Cython se pueden llamar funciones y librerías de C- Cython se integra de buena manera con NumPyPor ende Cython es muy atractivo para proyectos que usen Python y tengan requisitos de alto rendimiento. Estudiemos la sintaxis de Cython mediante algunos ejemplos **Ejemplo** A lo largo de esta lección utilizaremos como ejemplo el cálculo de la **Distancia euclidiana todos-con-todos**Sea un conjunto de con $N$ datos bidimensionales ($D=2$) donde queremos calcular la distancia euclidiana de cada dato con todos los demás, es decir una matriz donde el elemento $ij$ es $$\text{dist}_{ij} = \sqrt {\sum_{k=1}^D (x_{ik} - x_{jk})^2}$$A continuación se muestran dos códigos que cumplen este proposito y obtienen un resultado equivalente- El primero usa "Python puro" y calcula las distancias de forma secuencial- El segundo usa operaciones vectoriales de NumPy y calcula las distancias "al mismo tiempo"Estudie ambos códigos y asegúrese de comprenderlos antes de continuar<jupyter_code>data = np.random.randn(1000, 2)

def distancia_pares(data):    
    N, D = data.shape
    dist = np.zeros(shape=(N, N))
    for i in range(N):
        for j in range(i+1, N):
            for k in range(D):
                dist[i, j] += (data[i, k] - data[j, k])**2
            dist[i, j] = np.sqrt(dist[i, j])
            dist[j, i] = dist[i, j]
    return dist
            
def distancia_pares_numpy(data):
    return np.sqrt(np.sum((data.reshape(-1, 1, 2) - data.reshape(1, -1, 2))**2, axis=-1))

np.allclose(distancia_pares(data), distancia_pares_numpy(data))
time_pure = %timeit -r3 -n1 -q -o distancia_pares(data)
time_pure.average
time_numpy = %timeit -r3 -n1 -q -o distancia_pares_numpy(data)
time_numpy.average<jupyter_output><empty_output><jupyter_text>NumPy es<jupyter_code>time_pure.average/time_numpy.average<jupyter_output><empty_output><jupyter_text>veces más rápido que Python puro  Cython desde IPython/JupyterInstale cython en su ambiente de conda utilizando    conda install cythonLuego en IPython podemos cargar la extensión `cython` como se muestra a continuación<jupyter_code>%load_ext cython<jupyter_output><empty_output><jupyter_text>con lo que tendremos disponible la magia de bloque `%%cython`:::{important}Un bloque donde se use esta magia acepta lenguaje cython y se compila al ejecutarlo. Luego podremos llamar las funciones de ese bloque desde otros bloques regulares de Python:::Si surgen errores de compulación estos aparecen en la salida del bloque. :::{note}El bloque que empieza con `%%cython` está "desconectado" del resto del cuadernillo, por lo que debe tener sus propios `import`:::La magia tiene las siguientes opciones- `-a` (annotate) retorna un profile linea a linea indicando con amarillo las llamadas a CPython (mientras más llamadas más lento es nuestro código)- `-+` Usar C++ en lugar de C- `-c` Argumentos de compilación- `-l` librerías para linkear a nuestro código- `-L` directorio con librerías- `-I` directorio con cabeceras (include)   Veamos que ocurre al agregarle la magia al ejemplo anterior<jupyter_code>%%cython

import numpy as np

def distancia_pares_cython_inocente(data):    
    N, D = data.shape
    dist = np.zeros(shape=(N, N))
    for i in range(N):
        for j in range(i+1, N):
            for k in range(D):
                dist[i, j] += (data[i, k] - data[j, k])**2
            dist[i, j] = np.sqrt(dist[i, j])
            dist[j, i] = dist[i, j]
    return dist
np.allclose(distancia_pares_numpy(data), distancia_pares_cython_inocente(data))
time_cython_naive = %timeit -r3 -n1 -q -o distancia_pares_cython_inocente(data)
time_cython_naive.average
time_pure.average/time_cython_naive.average
time_numpy.average/time_cython_naive.average<jupyter_output><empty_output><jupyter_text>El resultado es idéntico pero como no hemos hecho ningún cambio el tiempo de ejecución es sólo levemente mejor que la versión en Python puroA continuación veremos como "cythonizar" nuestro código para ganar rendimiento  Mejora 1: Definiendo tipos en CythonEn Cython se definen tipos estáticos con el keyword `cdef` seguido del tipo y luego el nombre. Por ejemplo una variable de tipo `double` llamada `mi_variable`:```cythoncdef double mi_variable ```    Para los arreglos (ndarray) se usan ["memory-views"](https://cython.readthedocs.io/en/latest/src/userguide/memoryviews.htmlmemoryviews). Por ejemplo un *memory-view* para una arreglo de tres dimensiones:```cythoncdef double [:, :, :] mi_arreglo```    Y se puede ganar un poco más de rendimiento usando el operador `::1` para especificar si el arreglo es *row-major* (estilo C)```cythoncdef double [:, :, ::1] mi_arreglo```    o *column-major* (estilo Fortran)```cythoncdef double [::1, :, :] mi_arreglo```Veamos como queda el ejemplo introduciendo tipos estáticos<jupyter_code>%%cython 

import numpy as np

def distancia_pares_cython_estatico(double [:, ::1] data):
    # Definimos el tipo de N, D, dist y data
    cdef int N = data.shape[0]
    cdef int D = data.shape[1]
    dist = np.zeros(shape=(N, N), dtype=np.double)
    cdef double [:, ::1] dist_view = dist
    # También definimos los índices, se puede usar int o Py_ssize_t 
    cdef int i, j, k
    for i in range(N):
        for j in range(i+1, N):
            for k in range(D):
                dist_view[i, j] += (data[i, k] - data[j, k])**2
            dist_view[i, j] = np.sqrt(dist_view[i, j])
            dist_view[j, i] = dist_view[i, j]
    return dist
np.allclose(distancia_pares_numpy(data), distancia_pares_cython_estatico(data))
time_cython_static = %timeit -r3 -n1 -o -q distancia_pares_cython_estatico(data)
time_cython_static.average<jupyter_output><empty_output><jupyter_text>Speed-up con respecto a Python puro:<jupyter_code>time_pure.average/time_cython_static.average<jupyter_output><empty_output><jupyter_text>Speed-up con respecto a NumPy:<jupyter_code>time_numpy.average/time_cython_static.average<jupyter_output><empty_output><jupyter_text>Con solo definir el tipo de data, dist, $N$, $D$ y los índices hemos obtenido un *speed-up* importante con respecto a "Python puro" aunque sigue siendo más lento que NumPySi observaramos el código anotado por cython y revisemos la cantidad de llamadas a CPython de cada linea notaríamos que la linea 17 es particularmente conflictiva - involucra una gran cantidad de instrucciones- se llama NxN vecesesto se debe a que estamos usando la función de NumPy `np.sqrt`. Podemos obtener un rendimiento mucho mejor si usamos la implementación de raíz cuadrada de C  Mejora 2: Llamando funciones de C desde CythonEs posible llamar funciones de C desde Cython de forma sencilla. Consideremos como ejemplo la función `sqrt` de la [librería matemática estándar de C](http://pubs.opengroup.org/onlinepubs/9699919799/basedefs/math.h.html). Necesitamos dos para usar `sqrt` desde cython:- Importar la función `sqrt` desde la cabecera `math.h`: Esto se hace con el keyword `cdef extern from`- Compilar contra `libm`Veamos como queda nuestro ejemplo con esta modificación<jupyter_code>%%cython -l m

import numpy as np

cdef extern from "math.h":
    double sqrt(double)

def distancia_pares_cython_sqrtC(double [:, ::1] data):
    # Definimos el tipo de N, D, dist y data
    cdef int N = data.shape[0]
    cdef int D = data.shape[1]
    dist = np.zeros(shape=(N, N), dtype=np.double)
    cdef double [:, ::1] dist_view = dist
    # También definimos los índices, se puede usar int o Py_ssize_t 
    cdef Py_ssize_t i, j, k
    for i in range(N):
        for j in range(i+1, N):
            for k in range(D):
                dist_view[i, j] += (data[i, k] - data[j, k])**2
            dist_view[i, j] = sqrt(dist_view[i, j])
            dist_view[j, i] = dist_view[i, j]
    return dist
np.allclose(distancia_pares_numpy(data), distancia_pares_cython_sqrtC(data))
time_static_cfun = %timeit -r3 -n1 -q -o distancia_pares_cython_sqrtC(data)
time_static_cfun.average
time_numpy.average/time_static_cfun.average<jupyter_output><empty_output><jupyter_text>Con esta simple modificación hemos obtenido un tiempo de ejecución incluso menor que la implementación en NumPy, excelente!Es **sumamente importante** notar que usando `extern` podemos hacer interfaces [entre Python y casi cualquier código escrito en C](https://cython.readthedocs.io/en/latest/src/tutorial/clibraries.html)  Mejora 3: Deshabilitando comprobaciones para ir aun más rápidoPodemos hacer nuestro código más rápido (y más inseguro) deshabilitando las verificaciones que Python realiza por defectoEn Cython esto se logra usando decoradores que funcionan como directivas de compilación. Las opciones disponibles se encuentrán [aquí](https://cython.readthedocs.io/en/latest/src/userguide/source_files_and_compilation.htmlcompiler-directives). En este caso particular deshabilitaremos dos verificaciones: *boundscheck* y *wraparound*Al deshabilitarlas el código no comprobará si escribimos fuera del arreglo y tampoco traducirá índices negativosTambién aprovecharemos de hacer un cambio menor que nos será de utilidad más adelante:- Definiremos el tipo de data y dist de forma más conveniente usando `ctypedef`- Para esto incluiremos un modulo de Cython llamado numpy usando `cimport` que contiene definiciones tipo C<jupyter_code>%%cython -l m
import cython
cimport numpy as npc
import numpy as np

# Por conveniencia podemos definir el tipo de data y dist con 
ctypedef npc.float64_t TIPO_t
TIPO = np.float64

cdef extern from "math.h":
    TIPO_t sqrt(TIPO_t)

# Deshabilitamos las comprobaciones de Python:
@cython.boundscheck(False)
@cython.wraparound(False)
def distancia_pares_cython(TIPO_t [:, ::1] data):
    cdef int N = data.shape[0]
    cdef int M = data.shape[1]
    dist = np.zeros(shape=(N, N), dtype=TIPO)
    cdef TIPO_t [:, ::1] dist_view = dist
    cdef Py_ssize_t i, j, k
    for i in range(N):
        for j in range(i+1, N):
            for k in range(M):
                dist_view[i, j] += (data[i, k] - data[j, k])**2
            dist_view[i, j] = sqrt(dist_view[i, j])
            dist_view[j, i] = dist_view[i, j]
    return dist
np.allclose(distancia_pares_numpy(data), distancia_pares_cython(data))
time_cython_static_cfun_nochecks = %timeit -r3 -n1 -q -o distancia_pares_cython_sqrtC(data)
time_cython_static_cfun_nochecks.average
time_numpy.average/time_cython_static_cfun_nochecks.average<jupyter_output><empty_output><jupyter_text>Una modificación marginal en el código que nos da un leve speed-up con respecto al caso anterior    Mejora 4: Mayor flexibilidad con tipos de Cython fusionados Lo que hemos perdido al usar Cython con respecto a Python: **Flexibilidad**Si utilizamos un argumento que no sea `double` nuestra función en Cython retornará un error. Es nuestra responsabilidad evitar que esto ocurra<jupyter_code>data_float32 = data.astype(np.float32)<jupyter_output><empty_output><jupyter_text>Si ejecutaramos la siguiente instrucción: ```pythondistancia_pares_cython(data_float32)```Se nos retornaría una excepción ```pythonValueError: Buffer dtype mismatch, expected 'TIPO_t' but got 'float'``` Podemos resolver este tipo de problemas con los "tipos fusionados" de Cython. Si queremos fusionar dos tipos de datos "tipo1" y "tipo2" bajo el nombre "tipo3", lo escribimos usando```cythonctypedef fused tipo3:    tipo1    tipo2```Con esto Cython creará, de forma transparente, dos funciones en lugar de unaModifiquemos la función `distancia_pares_cython` para que acepte los tipos double e float como un tipo fusionado. En este caso no debemos olvidar importar la definición de `sqrtf` de `math.h`<jupyter_code>%%cython -l m
import cython
cimport numpy as npc
import numpy as np

# Tipo fusionado
ctypedef fused TIPO_t:
    npc.float32_t
    npc.float64_t # double
    
cdef extern from "math.h":
    npc.float32_t sqrtf(npc.float32_t) #Definición para float32
    npc.float64_t sqrt(npc.float64_t) # Definición para float64

@cython.boundscheck(False)
@cython.wraparound(False)
def distancia_pares_cython_multitipo(TIPO_t [:, ::1] data):
    cdef int N = data.shape[0]
    cdef int M = data.shape[1]
    # Comprobamos el tipo antes de crear el arreglo de numpy
    if TIPO_t is npc.float32_t:
        TIPO = np.float32
    else:
        TIPO = np.float64
    dist = np.zeros(shape=(N, N), dtype=TIPO)
    cdef TIPO_t [:, ::1] dist_view = dist
    cdef Py_ssize_t i, j, k
    for i in range(N):
        for j in range(i+1, N):
            for k in range(M):
                dist_view[i, j] += (data[i, k] - data[j, k])**2
            if TIPO_t is npc.float32_t:
                dist_view[i, j] = sqrtf(dist_view[i, j])
            else:
                dist_view[i, j] = sqrt(dist_view[i, j])
            dist_view[j, i] = dist_view[i, j]
    return dist
np.allclose(distancia_pares_numpy(data),
            distancia_pares_cython_multitipo(data))
np.allclose(distancia_pares_numpy(data_float32),
            distancia_pares_cython_multitipo(data_float32))
time_numpy_f32 = %timeit -r3 -n1 -q -o distancia_pares_numpy(data_float32)
time_numpy_f32.average
time_cython_f32 = %timeit -r3 -n1 -q -o distancia_pares_cython_multitipo(data_float32)
time_cython_f32.average
time_numpy_f32.average/time_cython_f32.average<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09024081115335868, 'fraction_numerical': 0.009315589353612166, 'mean_word_length': 4.376831345826235, 'pattern_counts': {'":': 3, '<': 49, '<?xml version=': 0, '>': 49, 'https://': 12, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '647179', 'license': 'Unlicense', 'n_tokens_mistral': 5687, 'n_tokens_neox': 5565, 'n_words': 1786, 'path': 'contents/hpc/cython.ipynb', 'repo_name': 'phuijse/PythonBook', 'size': '15780'}