Python – vectorizando una ventana deslizante

Estoy tratando de vectorizar una operación de ventana deslizante. Para el caso 1-d, un ejemplo útil podría ir en la línea de:

x= vstack((np.array([range(10)]),np.array([range(10)]))) x[1,:]=np.where((x[0,:]0),x[1,x[0,:]+1],x[1,:]) 

El valor n + 1 para cada valor actual para los índices <5. Pero me sale este error:

 x[1,:]=np.where((x[0,:]0),x[1,x[0,:]+1],x[1,:]) IndexError: index (10) out of range (0<=index<9) in dimension 1 

Curiosamente no obtendría este error por el valor n-1, lo que significaría índices más pequeños que 0. No parece importarle:

 x[1,:]=np.where((x[0,:]0),x[1,x[0,:]-1],x[1,:]) print(x) [[0 1 2 3 4 5 6 7 8 9] [0 0 1 2 3 5 6 7 8 9]] 

¿Hay alguna manera alrededor de esto? ¿Mi enfoque es totalmente incorrecto? Cualquier comentario será bienvenido.

EDITAR:

Esto es lo que me gustaría lograr, aplanar una matriz a una matriz numpy en la que quiero calcular la media de la vecindad 6×6 de cada celda:

 matriz = np.array([[1,2,3,4,5], [6,5,4,3,2], [1,1,2,2,3], [3,3,2,2,1], [3,2,1,3,2], [1,2,3,1,2]]) # matrix to vector vector2 = ndarray.flatten(matriz) ncols = int(shape(matriz)[1]) nrows = int(shape(matriz)[0]) vector = np.zeros(nrows*ncols,dtype='float64') # Interior pixels if ( (i % ncols) != 0 and (i+1) % ncols != 0 and i>ncols and i<ncols*(nrows-1)): vector[i] = np.mean(np.array([vector2[i-ncols-1],vector2[i-ncols],vector2[i-ncols+1],vector2[i-1],vector2[i+1],vector2[i+ncols-1],vector2[i+ncols],vector2[i+ncols+1]])) 

Si entiendo el problema correctamente, le gustaría tomar la media de todos los números 1 paso alrededor del índice, descuidando el índice.

He parcheado tu función para trabajar, creo que estabas yendo por algo como esto:

 def original(matriz): vector2 = np.ndarray.flatten(matriz) nrows, ncols= matriz.shape vector = np.zeros(nrows*ncols,dtype='float64') # Interior pixels for i in range(vector.shape[0]): if ( (i % ncols) != 0 and (i+1) % ncols != 0 and i>ncols and i 

Reescribí esto usando cortes y vistas:

 def mean_around(arr): arr=arr.astype(np.float64) out= np.copy(arr[:-2,:-2]) #Top left corner out+= arr[:-2,2:] #Top right corner out+= arr[:-2,1:-1] #Top center out+= arr[2:,:-2] #etc out+= arr[2:,2:] out+= arr[2:,1:-1] out+= arr[1:-1,2:] out+= arr[1:-1,:-2] out/=8.0 #Divide by # of elements to obtain mean cout=np.empty_like(arr) #Create output array cout[1:-1,1:-1]=out #Fill with out values cout[0,:]=0;cout[-1,:]=0;cout[:,0]=0;cout[:,-1]=0 #Set edges equal to zero return cout 

Usar np.empty_like y luego rellenar los bordes parecía un poco más rápido que np.zeros_like . Primero, verifiquemos que dan lo mismo usando su matriz matriz.

 print np.allclose(mean_around(matriz),original(matriz)) True print mean_around(matriz) [[ 0. 0. 0. 0. 0. ] [ 0. 2.5 2.75 3.125 0. ] [ 0. 3.25 2.75 2.375 0. ] [ 0. 1.875 2. 2. 0. ] [ 0. 2.25 2.25 1.75 0. ] [ 0. 0. 0. 0. 0. ]] 

Algunos tiempos:

 a=np.random.rand(500,500) print np.allclose(original(a),mean_around(a)) True %timeit mean_around(a) 100 loops, best of 3: 4.4 ms per loop %timeit original(a) 1 loops, best of 3: 6.6 s per loop 

Aproximadamente ~ 1500x de aceleración.

Parece un buen lugar para usar numba:

 def mean_numba(arr): out=np.zeros_like(arr) col,rows=arr.shape for x in xrange(1,col-1): for y in xrange(1,rows-1): out[x,y]=(arr[x-1,y+1]+arr[x-1,y]+arr[x-1,y-1]+arr[x,y+1]+\ arr[x,y-1]+arr[x+1,y+1]+arr[x+1,y]+arr[x+1,y-1])/8. return out nmean= autojit(mean_numba) 

Ahora vamos a comparar con todos los métodos presentados.

 a=np.random.rand(5000,5000) %timeit mean_around(a) 1 loops, best of 3: 729 ms per loop %timeit nmean(a) 10 loops, best of 3: 169 ms per loop #CT Zhu's answer %timeit it_mean(a) 1 loops, best of 3: 36.7 s per loop #Ali_m's answer %timeit fast_local_mean(a,(3,3)) 1 loops, best of 3: 4.7 s per loop #lmjohns3's answer %timeit scipy_conv(a) 1 loops, best of 3: 3.72 s per loop 

Una velocidad de 4x con numba arriba es bastante nominal, lo que indica que el código de numpy es tan bueno como se va a obtener. Saqué los otros códigos tal como se presentan, aunque tuve que cambiar la respuesta de @ CTZhu para incluir diferentes tamaños de matriz.

Parece que estás tratando de calcular una convolución 2D. Si eres capaz de usar scipy , te sugiero que intentes scipy.signal.convolve2d :

 matriz = np.random.randn(10, 10) # to average a 3x3 neighborhood kernel = np.ones((3, 3), float) # to compute the mean, divide by size of neighborhood kernel /= kernel.sum() average = scipy.signal.convolve2d(matriz, kernel) 

La razón por la que esto calcula la media de todos los barrios 3×3 se puede ver si “desenrolla” convolve2d en sus bucles constituyentes. Efectivamente (e ignorando lo que sucede en los bordes de las matrices de origen y de núcleo), se está calculando:

 X, Y = kernel.shape for i in range(matriz.shape[0]): for j in range(matriz.shape[1]): for ii in range(X): for jj in range(Y): average[i, j] += kernel[ii, jj] * matriz[i+ii, j+jj] 

Entonces, si cada valor en su kernel es 1 / (1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1) == 1/9, puede volver a escribir el código anterior como:

 for i in range(matriz.shape[0]): for j in range(matriz.shape[1]): average[i, j] = 1./9 * matriz[i:i+X, j:j+Y].sum() 

Que es exactamente lo mismo que calcular el promedio de los valores en matriz, sobre un área de 3×3, comenzando en i, j .

Una de las ventajas de hacer las cosas de esta manera es que puede cambiar fácilmente los pesos asociados con su vecindario al establecer los valores en su kernel de manera apropiada. Entonces, por ejemplo, si quisieras dar el valor central en cada vecindario el doble de peso que los demás, podrías construir tu kernel así:

 kernel = np.ones((3, 3), float) kernel[1, 1] = 2. kernel /= kernel.sum() 

y el código de convolución seguiría siendo el mismo, pero el cálculo arrojaría un tipo diferente de promedio (uno de “ponderación central”). Hay muchas posibilidades aquí; Esperemos que esto proporcione una buena abstracción para la tarea que está haciendo.

Da la casualidad de que hay una función en la biblioteca estándar de Scipy que calcula la media a través de ventanas deslizantes extremadamente rápido. Se llama uniform_filter . Puede usarlo para implementar su función de medio de vecindario de la siguiente manera:

 from scipy.ndimage.filters import uniform_filter def neighbourhood_average(arr, win=3): sums = uniform_filter(arr, win, mode='constant') * (win*win) return ((sums - arr) / (win*win - 1)) 

Esto devuelve una matriz X donde X[i,j] es el promedio de todos los vecinos de i,j en arr excluyendo i,j sí. Tenga en cuenta que la primera y la última columna y la primera y la última fila están sujetas a condiciones de límite, por lo que pueden no ser válidas para su aplicación (puede usar mode= para controlar la regla de límite si es necesario).

Debido a que uniform_filter utiliza un algoritmo de tiempo lineal altamente eficiente implementado en C directo (lineal solo en el tamaño de arr ), debería superar fácilmente cualquier otra solución, especialmente cuando la win es grande.

El problema radica en x[1,x[0,:]+1] , el índice para el segundo eje: x[0,:]+1 es [1 2 3 4 5 6 7 8 9 10] , en el que el índice 10 es más grande que la dimensión de x.

En el caso de x[1,x[0,:]-1] , el índice del 2do eje es [-1 0 1 2 3 4 5 6 7 8 9] , terminas obteniendo [9 0 1 2 3 4 5 6 7 8] , ya que 9 es el último elemento y tiene un índice de -1 . El índice del segundo elemento desde el final es -2 y así sucesivamente.

Con np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]-1],x[1,:]) y x[0,:]=[0 1 2 3 4 5 6 7 8 9] x[1,:] x[0,:]=[0 1 2 3 4 5 6 7 8 9] , lo que esencialmente está sucediendo es que la primera celda se toma de x[1,:] porque x[0,0] es 0 y x[0,:]<5)&(x[0,:]>0 x[1,:] x[0,:]<5)&(x[0,:]>0 es False . Los siguientes cuatro elementos se toman de x[1,x[0,:]-1] x[1,:] x[1,x[0,:]-1] . El rest es de x[1,:] Finalmente, el resultado es [0 0 1 2 3 4 5 6 7 8]

Puede parecer que está bien para una ventana deslizante de solo 1 celda, pero te sorprenderá con:

 >>> np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]-2],x[1,:]) array([0, 9, 0, 1, 2, 5, 6, 7, 8, 9]) 

Cuando intentas moverlo por una ventana de dos celdas.

Para este problema específico, si queremos mantener todo en una línea, esto hará:

 >>> for i in [1, 2, 3, 4, 5, 6]: print hstack((np.where(x[1,x[0,:]-i] 

Edición: ahora entiendo mejor su pregunta original, básicamente desea tomar una matriz 2D y calcular el promedio de las celdas N * N en cada celda. Eso es bastante común. En primer lugar, es probable que desee limitar N a números impares, de lo contrario es difícil definir un promedio de 2 * 2 alrededor de una celda. Supongamos que queremos un promedio de 3 * 3:

 #In this example, the shape is (10,10) >>> a1=\ array([[3, 7, 0, 9, 0, 8, 1, 4, 3, 3], [5, 6, 5, 2, 9, 2, 3, 5, 2, 9], [0, 9, 8, 5, 3, 1, 8, 1, 9, 4], [7, 4, 0, 0, 9, 3, 3, 3, 5, 4], [3, 1, 2, 4, 8, 8, 2, 1, 9, 6], [0, 0, 3, 9, 3, 0, 9, 1, 3, 3], [1, 2, 7, 4, 6, 6, 2, 6, 2, 1], [3, 9, 8, 5, 0, 3, 1, 4, 0, 5], [0, 3, 1, 4, 9, 9, 7, 5, 4, 5], [4, 3, 8, 7, 8, 6, 8, 1, 1, 8]]) #move your original array 'a1' around, use range(-2,2) for 5*5 average and so on >>> movea1=[a1[np.clip(np.arange(10)+i, 0, 9)][:,np.clip(np.arange(10)+j, 0, 9)] for i, j in itertools.product(*[range(-1,2),]*2)] #then just take the average >>> averagea1=np.mean(np.array(movea1), axis=0) #trim the result array, because the cells among the edges do not have 3*3 average >>> averagea1[1:10-1, 1:10-1] array([[ 4.77777778, 5.66666667, 4.55555556, 4.33333333, 3.88888889, 3.66666667, 4. , 4.44444444], [ 4.88888889, 4.33333333, 4.55555556, 3.77777778, 4.55555556, 3.22222222, 4.33333333, 4.66666667], [ 3.77777778, 3.66666667, 4.33333333, 4.55555556, 5. , 3.33333333, 4.55555556, 4.66666667], [ 2.22222222, 2.55555556, 4.22222222, 4.88888889, 5. , 3.33333333, 4. , 3.88888889], [ 2.11111111, 3.55555556, 5.11111111, 5.33333333, 4.88888889, 3.88888889, 3.88888889, 3.55555556], [ 3.66666667, 5.22222222, 5. , 4. , 3.33333333, 3.55555556, 3.11111111, 2.77777778], [ 3.77777778, 4.77777778, 4.88888889, 5.11111111, 4.77777778, 4.77777778, 3.44444444, 3.55555556], [ 4.33333333, 5.33333333, 5.55555556, 5.66666667, 5.66666667, 4.88888889, 3.44444444, 3.66666667]]) 

Creo que no necesitas aplanar tu matriz 2D, eso causa confusión. Además, si desea manejar los elementos de borde de manera diferente a no solo recortarlos, considere hacer matrices enmascaradas usando np.ma en el paso 'Mover su matriz original'.