Convertir un cálculo de bucle nested a Numpy para acelerar

Parte de mi progtwig de Python contiene el siguiente fragmento de código, donde se calcula una nueva cuadrícula según los datos encontrados en la cuadrícula anterior.

La cuadrícula es una lista bidimensional de flotadores. El código utiliza tres bucles for:

for t in xrange(0, t, step): for h in xrange(1, height-1): for w in xrange(1, width-1): new_gr[h][w] = gr[h][w] + gr[h][w-1] + gr[h-1][w] + t * gr[h+1][w-1]-2 * (gr[h][w-1] + t * gr[h-1][w]) gr = new_gr return gr 

El código es extremadamente lento para una cuadrícula grande y un tiempo grande t .

He intentado usar Numpy para acelerar este código, sustituyendo el bucle interno con:

 J = np.arange(1, width-1) new_gr[h][J] = gr[h][J] + gr[h][J-1] ... 

Pero los resultados producidos (los flotantes en la matriz) son aproximadamente un 10% más pequeños que sus contrapartes de cálculo de lista.

  • ¿Qué pérdida de precisión se espera al convertir listas de flotadores en una matriz de flotadores Numpy usando np.array (lista de pilotes) y luego hacer un cálculo?

  • ¿Cómo debo hacer para convertir un triple for-loop en un bonito y rápido código Numpy? (¿O hay otras sugerencias para acelerar el código significativamente?)

Si gr es una lista de flotadores, el primer paso si está buscando vectorizar con NumPy sería convertir gr en una matriz NumPy con np.array() .

A continuación, new_gr que tienes new_gr inicializado con ceros de forma (height,width) . Los cálculos que se realizan en los dos bucles más internos representan básicamente 2D convolution . Por lo tanto, puede usar signal.convolve2d con un kernel apropiado. Para decidir sobre el kernel , debemos observar los factores de escalado, hacer un kernel de 3 x 3 y negarlos para simular los cálculos que estamos haciendo con cada iteración. Por lo tanto, tendría una solución vectorizada con los dos bucles más internos eliminados para un mejor rendimiento, como por ejemplo:

 import numpy as np from scipy import signal # Get the scaling factors and negate them to get kernel kernel = -np.array([[0,1-2*t,0],[-1,1,0,],[t,0,0]]) # Initialize output array and run 2D convolution and set values into it out = np.zeros((height,width)) out[1:-1,1:-1] = signal.convolve2d(gr, kernel, mode='same')[1:-1,:-2] 

Verifique las pruebas de salida y tiempo de ejecución

Definir funciones:

 def org_app(gr,t): new_gr = np.zeros((height,width)) for h in xrange(1, height-1): for w in xrange(1, width-1): new_gr[h][w] = gr[h][w] + gr[h][w-1] + gr[h-1][w] + t * gr[h+1][w-1]-2 * (gr[h][w-1] + t * gr[h-1][w]) return new_gr def proposed_app(gr,t): kernel = -np.array([[0,1-2*t,0],[-1,1,0,],[t,0,0]]) out = np.zeros((height,width)) out[1:-1,1:-1] = signal.convolve2d(gr, kernel, mode='same')[1:-1,:-2] return out 

Verificar

 In [244]: # Inputs ...: gr = np.random.rand(40,50) ...: height,width = gr.shape ...: t = 1 ...: In [245]: np.allclose(org_app(gr,t),proposed_app(gr,t)) Out[245]: True Timings - In [246]: # Inputs ...: gr = np.random.rand(400,500) ...: height,width = gr.shape ...: t = 1 ...: In [247]: %timeit org_app(gr,t) 1 loops, best of 3: 2.13 s per loop In [248]: %timeit proposed_app(gr,t) 10 loops, best of 3: 19.4 ms per loop 

@Divakar, probé un par de variaciones en su org_app . La versión totalmente vectorizada es:

 def org_app4(gr,t): new_gr = np.zeros((height,width)) I = np.arange(1,height-1)[:,None] J = np.arange(1,width-1) new_gr[I,J] = gr[I,J] + gr[I,J-1] + gr[I-1,J] + t * gr[I+1,J-1]-2 * (gr[I,J-1] + t * gr[I-1,J]) return new_gr 

Si bien la mitad de la velocidad de tu proposed_app , está más cercana en estilo al original. Y, por lo tanto, puede ayudar a comprender cómo se pueden vectorizar los bucles nesteds.

Un paso importante es la conversión de I en una matriz de columnas, de modo que, en conjunto I,J indexen un bloque de valores.