Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar <ruta a mi directorio> por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker.
docker run --rm -v <ruta a mi directorio>:/datos --name jupyterlab_numerical -p 8888:8888 -d palmoreck/jupyterlab_numerical:1.1.0password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_numericalDocumentación de la imagen de docker palmoreck/jupyterlab_numerical:1.1.0 en liga.
Esta nota utiliza métodos vistos en 1.5.Integracion_numerica
Documentación de cython:
La siguiente celda muestra el modo de utilizar el comando magic de %pip para instalar paquetes desde jupyterlab. Ver liga para magic commands.
Instalamos cython:
In [1]:
%pip install -q --user cython
La siguiente celda reiniciará el kernel de IPython para cargar los paquetes instalados en la celda anterior. Dar Ok en el mensaje que salga y continuar con el contenido del notebook.
In [2]:
import IPython
app = IPython.Application.instance()
app.kernel.do_shutdown(True)
Out[2]:
De las opciones más sencillas que tenemos a nuestra disposición para resolver bottlenecks en nuestro programa es hacer que nuestro código haga menos trabajo. ¿Cómo podemos hacer esto? <- compilando nuestro código a código de máquina para que el código en Python ejecute menos instrucciones.
¿Por qué puede ser lenta la ejecución de un bloque de código en Python (o en algún otro lenguaje tipo intérprete*)? La verificación de tipo de datos (si son int, double o string), los objetos temporales que se crean por tipo de dato (un objeto tipo int en Python tiene asociado un objeto de alto nivel con el que interactuamos pero que causa overhead) y las llamadas a funciones de alto nivel (por ejemplo las que ayudan a almacenar al objeto en memoria) son tres de las fuentes que hacen a Python (y a otro lenguaje tipo intérprete como R o Matlab) lento. También otras fuentes que responden la pregunta son:
Desde el punto de vista de la memoria de la máquina, el número de referencias a un objeto y las copias entre objetos.
No es posible vectorizar un cálculo sin el uso de extensiones (por ejemplo paquetes como numpy).
*Ver liga) para revisar lo que es un lenguaje tipo intérprete.
Un paquete para resolver lo anterior es Cython que requiere que escribamos en un lenguaje híbrido entre Python y C. Si bien a l@s integrantes de un equipo de desarrollo que no saben C éste cambio reducirá la velocidad de desarrollo, en la práctica si se tiene un bottleneck que no ha podido resolverse con herramientas como el cómputo en paralelo o vectorización, se recomienda utilizar Cython para regiones pequeñas del código y resolver el bottleneck del programa.
Cython es un compilador que convierte type-annotated Python y C like (instrucciones escritas en Python pero en una forma tipo lenguaje C) en un módulo compilado que funciona como una extensión de Python. Este módulo puede ser importado como un módulo regular de Python utilizando import.
Además Cython tiene ya un buen tiempo en la comunidad (2007 aprox.), es altamente usado y es de las herramientas preferidas para código tipo CPU-bound (un gran porcentaje del código es uso de CPU vs uso de memoria o I/O). También soporta a la API OpenMP que veremos en el capítulo de cómputo en paralelo para aprovechar los múltiples cores o CPU's de una máquina.
Para ver más historia de Cython ir a la referencia 1. de esta nota o a la liga.
for) en los que se realizan operaciones matemáticas típicamente no vectorizadas o que no pueden vectorizarse*. Esto es, códigos en los que las instrucciones son básicamente sólo Python sin utilizar paquetes externos. Además, si en el ciclo las variables no cambian de su tipo (por ejemplo de int a float) entonces es una blanco perfecto para ganancia en velocidad al compilar a código de máquina.*Si tu código de Python llama a operaciones vectorizadas vía numpy podría ser que no se ejecute más rápido tu código después de compilarlo.
strings o a una base de datos). Programas que tengan alta carga de I/O también es poco probable que muestren ganancias significativas.En general es poco probable que tu código compilado se ejecute más rápido que un código en C "bien aceitado" y también es poco probable que se ejecute más lento. Es muy posible que el código C generado desde tu Python pueda alcanzar las velocidades de un código escrito en C, a menos que la persona que programó en C tenga un gran conocimiento de formas de hacer que el código de C se ajuste a la arquitectura de la máquina sobre la que se ejecutan los códigos.
No olvidar: es importante fijar de forma aproximada el tiempo objetivo que se desea alcanzar para un código que escribamos. Si bien el perfilamiento y la compilación son herramientas para resolver los bottlenecks, debemos tomar en cuenta el tiempo objetivo fijado y ser práctic@s en el desarrollo, no podemos por siempre estar optimizando nuestro código.
Cython puede utilizarse vía un script setup.py que compila un módulo pero también puede utilizarse en IPython vía un comando magic.
In [1]:
import math
import time
from scipy.integrate import quad
Para este caso requerimos tres archivos:
1) El código escrito en Python que será compilado en un archivo con extensión .pyx.
2) Un archivo setup.py que contiene las instrucciones para llamar a Cython y cree el módulo compilado.
3) El código escrito en Python que importará el módulo compilado (puede pensarse como nuestro main).
1) Función a compilar en un archivo .pyx:
In [3]:
%%file Rcf_cython.pyx
def Rcf(f,a,b,n): #Rcf: rectángulo compuesto para f
"""
Compute numerical approximation using rectangle or mid-point method in
an interval.
Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
Args:
f (lambda expression): lambda expression of integrand
a (int): left point of interval
b (int): right point of interval
n (int): number of subintervals
Returns:
Rcf (float)
"""
h_hat=(b-a)/n
nodes=[a+(i+1/2)*h_hat for i in range(0,n)]
sum_res=0
for node in nodes:
sum_res=sum_res+f(node)
return h_hat*sum_res
2) Archivo setup.py que contiene las instrucciones para el build:
In [4]:
%%file setup.py
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
setup(cmdclass = {'build_ext': build_ext},
ext_modules = [Extension("Rcf_compiled",
["Rcf_cython.pyx"])]
)
Compilar desde la línea de comandos:
In [5]:
%%bash
python3 setup.py build_ext --inplace #inplace para compilar el módulo en el directorio
#actual
Notas:
La compilación debe hacerse cada vez que se modifica el código de la función Rcf del archivo Rcf_cython.pyx o cambia el setup.py.
Obsérvese en el directorio donde se encuentra la nota que se generó un archivo Rcf_cython.c.
3) Importar módulo compilado y ejecutarlo:
In [6]:
f=lambda x: math.exp(-x**2) #using math library
In [7]:
import Rcf_compiled
In [8]:
n=10**6
start_time = time.time()
aprox=Rcf_compiled.Rcf(f,0,1,n)
end_time = time.time()
In [9]:
aprox
Out[9]:
In [10]:
secs = end_time-start_time
print("Rcf tomó",secs,"segundos" )
Recuérdese revisar el error relativo:
In [11]:
def err_relativo(aprox, obj):
return math.fabs(aprox-obj)/math.fabs(obj) #obsérvese el uso de la librería math
In [12]:
obj, err = quad(f, 0, 1)
err_relativo(aprox,obj)
Out[12]:
Ejercicio: investigar por qué se tiene un error relativo del orden de $10^{-7}$ y no de mayor precisión como se verá más abajo con el archivo Rcf_cython2.pyx.
In [13]:
%timeit -n 5 -r 10 Rcf_compiled.Rcf(f,0,1,n)
Obs: El error relativo anterior es más grande del que se obtenía anteriomente, por ello se utilizará el módulo cythonize (haciendo pruebas e investigando un poco se obtuvo la precisión de antes).
In [14]:
%%file Rcf_cython2.pyx
def Rcf(f,a,b,n): #Rcf: rectángulo compuesto para f
"""
Compute numerical approximation using rectangle or mid-point method in
an interval.
Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
Args:
f (lambda expression): lambda expression of integrand
a (int): left point of interval
b (int): right point of interval
n (int): number of subintervals
Returns:
Rcf (float)
"""
h_hat=(b-a)/n
nodes=[a+(i+1/2)*h_hat for i in range(0,n)]
sum_res=0
for node in nodes:
sum_res=sum_res+f(node)
return h_hat*sum_res
In [15]:
%%file setup2.py
from distutils.core import setup
from Cython.Build import cythonize
setup(ext_modules = cythonize("Rcf_cython2.pyx",
compiler_directives={'language_level' : 3})
)
#es posible que la solución del ejercicio anterior tenga que ver con el Warning
#y uso de la directiva language_level
In [16]:
%%bash
python3 setup2.py build_ext --inplace
In [17]:
import Rcf_cython2
In [18]:
n=10**6
start_time = time.time()
aprox=Rcf_cython2.Rcf(f,0,1,n)
end_time = time.time()
In [19]:
aprox
Out[19]:
In [20]:
secs = end_time-start_time
print("Rcf tomó",secs,"segundos" )
Revisar error relativo:
In [21]:
err_relativo(aprox,obj)
Out[21]:
Este comando se ocupa de escribir el archivo .pyx, compilarlo con setup.py e importarlo en el notebook. Ver cythonmagic.
In [22]:
%load_ext Cython
In [23]:
%%cython
def Rcf(f,a,b,n): #Rcf: rectángulo compuesto para f
"""
Compute numerical approximation using rectangle or mid-point method in
an interval.
Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
Args:
f (lambda expression): lambda expression of integrand
a (int): left point of interval
b (int): right point of interval
n (int): number of subintervals
Returns:
Rcf (float)
"""
h_hat=(b-a)/n
nodes=[a+(i+1/2)*h_hat for i in range(0,n)]
sum_res=0
for node in nodes:
sum_res=sum_res+f(node)
return h_hat*sum_res
In [24]:
n=10**6
start_time = time.time()
aprox=Rcf(f,0,1,n)
end_time = time.time()
In [25]:
secs = end_time-start_time
print("Rcf tomó",secs,"segundos" )
In [26]:
err_relativo(aprox,obj)
Out[26]:
In [28]:
%timeit -n 5 -r 10 Rcf(f,0,1,n)
Cython tiene una opción de annotation que generará un archivo con extensión .html que se puede visualizar en jupyterlab o en un browser.
Cada línea puede ser expandida haciendo un doble click que mostrará el código C generado. Líneas más amarillas refieren a más llamadas en la máquina virtual de Python (por máquina virtual de Python entiéndase la maquinaria que utiliza Python para traducir el lenguaje bytecode ), mientras que líneas más blancas significan "más código en C y no Python". El objetivo es remover la mayor cantidad de líneas amarillas posibles (pues típicamente son costosas en tiempo y si las líneas están dentro de loops son todavía más costosas) y terminar con códigos cuyas annotations sean lo más blancas posibles. Concentra tu atención en las líneas que son amarillas y están dentro de los loops, no pierdas tu tiempo en líneas amarillas que están fuera de loops y que no causan una ejecución lenta (para identificar esto perfila tu código).
In [29]:
%%bash
$HOME/.local/bin/cython -a Rcf_cython.pyx
Ver archivo creado: Rcf_cython.html
In [30]:
%%cython?
In [ ]:
%%cython -a
def Rcf(f,a,b,n): #Rcf: rectángulo compuesto para f
"""
Compute numerical approximation using rectangle or mid-point method in
an interval.
Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
Args:
f (lambda expression): lambda expression of integrand
a (int): left point of interval
b (int): right point of interval
n (int): number of subintervals
Returns:
Rcf (float)
"""
h_hat=(b-a)/n
nodes=[a+(i+1/2)*h_hat for i in range(0,n)]
sum_res=0
for node in nodes:
sum_res=sum_res+f(node)
return h_hat*sum_res
Nota: la imagen anterior es un screenshot que generé ejecutando la celda anterior.
Obs: para este ejemplo la línea $18$ es muy amarilla y está dentro del loop. Recuérdese de la nota 1.6.Perfilamiento_Python.ipynb que es una línea en la que se gasta parte del tiempo total de ejecución del código.
Una primera opción que tenemos es crear los nodos para el método de integración dentro del loop y separar el llamado a la list comprehension nodes=[a+(i+1/2)*h_hat for i in range(0,n)]:
In [ ]:
%%cython -a
def Rcf2(f,a,b,n): #Rcf: rectángulo compuesto para f
"""
Compute numerical approximation using rectangle or mid-point method in
an interval.
Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
Args:
f (lambda expression): lambda expression of integrand
a (int): left point of interval
b (int): right point of interval
n (int): number of subintervals
Returns:
Rcf (float)
"""
h_hat=(b-a)/n
sum_res=0
for i in range(0,n):
x=a+(i+1/2.0)*h_hat
sum_res+=f(x)
return h_hat*sum_res
Nota: la imagen anterior es un screenshot que generé ejecutando la celda anterior.
In [41]:
n=10**6
start_time = time.time()
aprox=Rcf2(f,0,1,n)
end_time = time.time()
In [42]:
secs = end_time-start_time
print("Rcf2 tomó",secs,"segundos" )
In [43]:
err_relativo(aprox,obj)
Out[43]:
In [44]:
%timeit -n 5 -r 10 Rcf2(f,0,1,n)
Obs: para este ejemplo las líneas $17$ y $18$ son muy amarillas y están dentro del loop. Además son líneas que involucran tipos de datos que no cambiarán en la ejecución de cada loop. Nos enfocamos a hacerlas más blancas... Una primera opción es declarar los tipos de objetos que están involucrados en el loop utilizando la sintaxis cdef:
In [ ]:
%%cython -a
def Rcf3(f,double a,double b,unsigned int n): #obsérvese la declaración de los tipos
"""
Compute numerical approximation using rectangle or mid-point method in
an interval.
Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
Args:
f (lambda expression): lambda expression of integrand
a (int): left point of interval
b (int): right point of interval
n (int): number of subintervals
Returns:
Rcf (float)
"""
cdef unsigned int i #obsérvese la declaración de los tipos
cdef double x,sum_res, h_hat #obsérvese la declaración de los tipos
h_hat=(b-a)/n
sum_res=0
for i in range(0,n):
x=a+(i+1/2.0)*h_hat
sum_res+=f(x)
return h_hat*sum_res
Nota: la imagen anterior es un screenshot que generé ejecutando la celda anterior.
In [46]:
n=10**6
start_time = time.time()
aprox=Rcf3(f,0,1,n)
end_time = time.time()
In [48]:
secs = end_time-start_time
print("Rcf3 tomó",secs,"segundos" )
In [49]:
err_relativo(aprox,obj)
Out[49]:
In [50]:
%timeit -n 5 -r 10 Rcf3(f,0,1,n)
Obs: al definir tipos, éstos sólo serán entendidos por Cython y no por Python. Cython utiliza estos tipos para convertir el código de Python a objetos de C, éstos objetos no tienen que convertirse de vuelta a objetos de Python. Entonces perdemos flexibilidad pero ganamos velocidad.
Y podemos bajar más el tiempo al definir la función que será utilizada:
In [ ]:
%%cython -a
import math
def Rcf4(double a,double b,unsigned int n): #Rcf: rectángulo compuesto para f
"""
Compute numerical approximation using rectangle or mid-point method in
an interval.
Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
Args:
f (lambda expression): lambda expression of integrand
a (int): left point of interval
b (int): right point of interval
n (int): number of subintervals
Returns:
Rcf (float)
"""
cdef unsigned int i
cdef double x,sum_res, h_hat
h_hat=(b-a)/n
sum_res=0
for i in range(0,n):
x=a+(i+1/2.0)*h_hat
sum_res+=math.exp(-x**2)
return h_hat*sum_res
Nota: la imagen anterior es un screenshot que generé ejecutando la celda anterior.
In [52]:
n=10**6
start_time = time.time()
aprox=Rcf4(0,1,n)
end_time = time.time()
In [53]:
secs = end_time-start_time
print("Rcf4 tomó",secs,"segundos" )
In [54]:
err_relativo(aprox,obj)
Out[54]:
In [55]:
%timeit -n 5 -r 10 Rcf4(0,1,n)
Obs: estamos ganando velocidad pues el compilador de C gcc puede optimizar funciones de bajo nivel para operar en los bytes que están asociados a las variables y no realiza llamadas a la máquina virtual de Python.
Ejercicios
Referencias
Otras referencias: