ciclo - Vectorización de Python anidada para bucles

ciclo for en python

Question

Agradecería un poco de ayuda para encontrar y comprender una forma pitónica para optimizar las siguientes manipulaciones de matriz en bucles anidados:

def _func(a, b, radius):
"Return 0 if a>b, otherwise return 1"
if distance.euclidean(a, b) < radius:
return 1
else:
return 0
def _make_mask(volume, roi, radius):
mask = numpy.zeros(volume.shape)
for x in range(volume.shape[0]):
for y in range(volume.shape[1]):
for z in range(volume.shape[2]):
mask[x, y, z] = _func((x, y, z), roi, radius)
return mask

Donde volume.shape (182, 218, 200) y roi.shape (3,) son ambos tipos de ndarray ; y el radius es un int

Answer 1

Digamos que primero construyes una matriz xyzy :

import itertools
xyz = [np.array(p) for p in itertools.product(range(volume.shape[0]), range(volume.shape[1]), range(volume.shape[2]))]

Ahora, usando numpy.linalg.norm ,

np.linalg.norm(xyz - roi, axis=1) < radius

comprueba si la distancia para cada tupla desde roi es menor que el radio.

Finalmente, simplemente reshape el resultado a las dimensiones que necesita.

Answer 2

Enfoque n. ° 1

Aquí hay un enfoque vectorizado:

m,n,r = volume.shape
x,y,z = np.mgrid[0:m,0:n,0:r]
X = x - roi[0]
Y = y - roi[1]
Z = z - roi[2]
mask = X**2 + Y**2 + Z**2 < radius**2

Posible mejora: probablemente podamos acelerar el último paso con el módulo numexpr -

import numexpr as ne
mask = ne.evaluate(''X**2 + Y**2 + Z**2 < radius**2'')

Enfoque n. ° 2

También podemos construir gradualmente los tres rangos correspondientes a los parámetros de forma y realizar la resta contra los tres elementos de roi sobre la marcha sin crear las mallas como se hizo anteriormente con np.mgrid . Esto se beneficiaría con el uso de la broadcasting con fines de eficiencia. La implementación se vería así:

m,n,r = volume.shape
vals = ((np.arange(m)-roi[0])**2)[:,None,None] + /
((np.arange(n)-roi[1])**2)[:,None] + ((np.arange(r)-roi[2])**2)
mask = vals < radius**2

Versión simplificada: Gracias a @Bi Rico por sugerir una mejora aquí, ya que podemos usar np.ogrid para realizar esas operaciones de una manera un poco más concisa, así:

m,n,r = volume.shape
x,y,z = np.ogrid[0:m,0:n,0:r]-roi
mask = (x**2+y**2+z**2) < radius**2

Prueba de tiempo de ejecución

Definiciones de funciones

def vectorized_app1(volume, roi, radius):
m,n,r = volume.shape
x,y,z = np.mgrid[0:m,0:n,0:r]
X = x - roi[0]
Y = y - roi[1]
Z = z - roi[2]
return X**2 + Y**2 + Z**2 < radius**2
def vectorized_app1_improved(volume, roi, radius):
m,n,r = volume.shape
x,y,z = np.mgrid[0:m,0:n,0:r]
X = x - roi[0]
Y = y - roi[1]
Z = z - roi[2]
return ne.evaluate(''X**2 + Y**2 + Z**2 < radius**2'')
def vectorized_app2(volume, roi, radius):
m,n,r = volume.shape
vals = ((np.arange(m)-roi[0])**2)[:,None,None] + /
((np.arange(n)-roi[1])**2)[:,None] + ((np.arange(r)-roi[2])**2)
return vals < radius**2
def vectorized_app2_simplified(volume, roi, radius):
m,n,r = volume.shape
x,y,z = np.ogrid[0:m,0:n,0:r]-roi
return (x**2+y**2+z**2) < radius**2

Tiempos -

In [106]: # Setup input arrays
...: volume = np.random.rand(90,110,100) # Half of original input sizes
...: roi = np.random.rand(3)
...: radius = 3.4
...:
In [107]: %timeit _make_mask(volume, roi, radius)
1 loops, best of 3: 41.4 s per loop
In [108]: %timeit vectorized_app1(volume, roi, radius)
10 loops, best of 3: 62.3 ms per loop
In [109]: %timeit vectorized_app1_improved(volume, roi, radius)
10 loops, best of 3: 47 ms per loop
In [110]: %timeit vectorized_app2(volume, roi, radius)
100 loops, best of 3: 4.26 ms per loop
In [139]: %timeit vectorized_app2_simplified(volume, roi, radius)
100 loops, best of 3: 4.36 ms per loop

Entonces, como siempre broadcasting mostrando su magia para una aceleración loca de casi 10,000x sobre el código original y más de 10x mejor que la creación de mallas mediante el uso de operaciones transmitidas sobre la marcha.