Simulación

Cómputo paralelo en Python

Cualquier operación donde no todos los pasos dependen de los resultados de los pasos anteriores es un candidato para ser paralelizado.

Un pool es un grupo de dos o más procesos que ejecutan sus tareas de manera simultánea.

Aplicación de funciones a iterables de forma paralela

>>> from multiprocessing import Pool
>>> def f(x):
...     return 2**x
...
>>> if __name__ == '__main__':
...     with Pool(3) as p:
...         print(p.map(f, range(1, 11)))
[2, 4, 8, 16, 32, 64, 128, 256, 512, 1024]

El paquete multiprocessing facilita la ejecución paralela en Python; forma parte de la instalación básica por lo cual no es necesario instalarlo — se carga para su uso con la instrucción import. Es importante incluir su uso dentro de main por la forma en que accede las variables.

Se puede determinar el número de núcleos disponibles en un sistema con la llamada

>>> from multiprocessing import cpu_count
>>> cpu_count()
8

>>> import psutil # instalar con pip
>>> psutil.cpu_count(logical = False)
4
>>> psutil.cpu_count(logical = True)
8

No conviene utilizar todos los núcleos en las tareas de Python, ya que el sistema operativo y su interfaz de usuario suelen requerir mantenerse responsivos durante los cálculos realizados, por lo cual es recomendable crear un cluster con un número menor de núcleos, por ejemplo with Pool(cpu_count() - 1) as p:.

El pool únicamente existe dentro del bloque de with en el cual se ha creado y se libera automáticamente al concluir el bloque.

Lo valores de las variables presentes en el ambiente de trabajo de cargan al Pool al momento de crearlo; si cambian y ese cambio ocupa reflejarse en los procesos en paralelo, hay que volver a crear el pool.

>>> desde = 3
>>> hasta = 8
>>> base = 2
>>> def f(x):
...     return base**x
...
>>> from multiprocessing import Pool
>>> if __name__ == '__main__':
...     with Pool(3) as p:
...         print(p.map(f, range(desde, hasta + 1)))
...         base = 3
...         print(p.map(f, range(desde, hasta + 1)))
...
[8, 16, 32, 64, 128, 256]
[8, 16, 32, 64, 128, 256]
>>> base = 3
>>> if __name__ == '__main__':
...     with Pool(3) as p:
...         print(p.map(f, range(desde, hasta + 1)))
[27, 81, 243, 729, 2187, 6561]

Otra opción es crear funciones que toman más de un parámetro y usar la rutina starmap:

>>> desde = 3
>>> hasta = 8
>>> def f(x, b):
...     return b**x
...
>>> from multiprocessing import Pool
>>> if __name__ == '__main__':
...     with Pool(3) as p:
...         print(p.starmap(f, zip([2] * 10, range(desde, hasta + 1))))
...         print(p.starmap(f, zip([3] * 10, range(desde, hasta + 1))))
...         print(p.starmap(f, zip([4] * 10, range(desde, hasta + 1))))
...
[9, 16, 25, 36, 49, 64]
[27, 64, 125, 216, 343, 512]
[81, 256, 625, 1296, 2401, 4096]
>>> def f(x, b, c):
...     return b**x - c
...
>>> c = [7 * i for i in range(10)]
>>> if __name__ == '__main__':
...     with Pool(3) as p:
...         for b in range(2, 5):
...             base = [b] * 10
...             vars = [j for j in range(desde, hasta + 1)]
...             params = zip(base, vars, c)
...             print(p.starmap(f, params))
[9, 9, 11, 15, 21, 29]
[27, 57, 111, 195, 315, 477]
[81, 249, 611, 1275, 2373, 4061]

Cómputo con GPU

Tarjetas gráficas contemporaneas contienen procesdaores gráficos con bastante núcleos. El paquete pyopencl permite aprovechar esos núcleos para cálculos. Su instalación requiere la herramienta pybind11 con pip3 y la presencia a las librerías de OpenCL (no todas las computadoras tienen chips o drivers compatibles). Finalizando la instalación del paquete, igual como al cargarlo al uso, se detectan GPUs presentes en la computadora en cuestión.

Algunas versiones de pip requieren un ligero hack en /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pip/req.py (ajustar según su sistema) para que no marquen error al agregar pybind11:

class Hack: # agregado para compatibilidad
    def __init__(self, text):
        self.req = text

def parse_requirements(filename, session=None): # otro hack para que haya session
    return [Hack(line) for line in (line.strip() for line in open(filename)) if line and not line.startswith("#")]
  

En una iMac con una GPU de Radeon dice

>>> import pyopencl
>>> from pyopencl.tools import get_test_platforms_and_devices
>>> get_test_platforms_and_devices()
[(<pyopencl.Platform 'Apple' at 0x7fff0000>, [<pyopencl.Device 'Intel(R) Core(TM) i7-6700K CPU @ 4.00GHz' on 'Apple' at 0xffffffff>, <pyopencl.Device 'ATI Radeon R9 M395X Compute Engine' on 'Apple' at 0x1021c00>])]

En una MacBook con una GPU de Intel dice

>>> import pyopencl
>>> from pyopencl.tools import get_test_platforms_and_devices
>>> get_test_platforms_and_devices()
[(<pyopencl.Platform 'Apple' at 0x7fff0000>, [<pyopencl.Device 'Intel(R) Core(TM) M-5Y71 CPU @ 1.20GHz' on 'Apple' at 0xffffffff>, <pyopencl.Device 'Intel(R) HD Graphics 5300' on 'Apple' at 0x1024500>])]

Por ejemplo, el producto de matrices que suele ser la prueba estándar para ver si la GPU está haciendo lo que debe (véase las notas de Jan Verschelde hace lo siguiente con el GPU de Radeon. Nota que el código para el GPU viene escrito en el lenguaje C.

# adaptado de http://homepages.math.uic.edu/~jan/mcs572/mcs572notes/lec29.html
import pyopencl as cl
import numpy as np
import os
os.environ['PYOPENCL_COMPILER_OUTPUT'] = '1'
os.environ['PYOPENCL_CTX'] = '1'
n = 7000
M1 = np.random.rand(n, n).astype(np.float32)
M2 = np.random.rand(n, n).astype(np.float32)
M3 = np.zeros((n, n), dtype=np.float32)
platforms = cl.get_platforms()
ctx = cl.Context(dev_type=cl.device_type.ALL, 
                 properties=[(cl.context_properties.PLATFORM, platforms[0])])
queue = cl.CommandQueue(ctx)
mf = cl.mem_flags
# memoria para procesar en GPU
b1 = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=M1)
b2 = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=M2)
b3 = cl.Buffer(ctx, mf.WRITE_ONLY, M3.nbytes)
# lo que se hace en GPU
cCode = '''PONER AQUI LO QUE VIENE EN C ABAJO'''
prg = cl.Program(ctx, cCode).build()
from time import time
antes = time()
prg.multiply(queue, M3.shape, None, np.uint16(n), b1, b2, b3)
cl.enqueue_copy(queue, np.empty_like(M3), b3)
print(time() - antes)
antes = time()
np.matmul(M1, M2)
print(time() - antes)

__kernel void multiply(ushort n, __global float *a, __global float *b, __global float *c) {
  int gid = get_global_id(0);
  c[gid] = 0.0f;
  int rowC = gid / n;
  int colC = gid % n;
  __global float *pA = &a[rowC * n];
  __global float *pB = &b[colC];
  for (int k = 0; k < n; k++) {
    pB = &b[colC + k * n];
    c[gid] += (*(pA++)) * (*pB); 
  }
}

Lo chistoso es que las rutinas de numpy son tan eficientes que no se logra ganarlos en este caso particular. Abajo están las salidas con la iMac y la MacBook, respectivamente. En la laptop que tiene un procesador más modesto, es casi lo mismo, pero en la iMac que tiene un procesador potente, el GPU pierde por el setup overhead.

$ python3 gpu.py
10.52643609046936
1.9072740077972412

$ python3 gpu.py
8.462943077087402
8.137696027755737