Python - vectorizing a sliding window
numpy, python, scipy
Solution
If I understand the problem correctly you would like to take the mean of all numbers 1 step around the index, neglecting the index.
I have patched your function to work, I believe you were going for something like this:
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<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]]))
I rewrote this using using slicing and views:
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
Using `np.empty_like` and then filling the edges seemed slightly faster then `np.zeros_like`. First lets double check they give the same thing using your `matriz` array.
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. ]]
Some timings:
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
Roughly ~1500x speedup.
Looks like a good place to use 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)
Now lets compare against all presented methods.
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
A 4x speed with numba up is pretty nominal indicating that the numpy code is about as good as its going to get. I pulled the other codes as presented, although I did have to change @CTZhu's answer to include different array sizes.
Problem
I'm trying to vectorize a sliding window operation. For the 1-d case a helpful example could go along the lines of: ``` x= vstack((np.array([range(10)]),np.array([range(10)]))) x[1,:]=np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]+1],x[1,:]) ``` The n+1 value for each current value for indices <5. But I get this error: ``` x[1,:]=np.where((x[0,:]<2)&(x[0,:]>0),x[1,x[0,:]+1],x[1,:]) IndexError: index (10) out of range (0<=index<9) in dimension 1 ``` Curiously I wouldn't get this error for the n-1 value which would mean indices smaller than 0. It doesn't seem to mind: ``` x[1,:]=np.where((x[0,:]<5)&(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]] ``` Is there anyway around this? is my approach totally wrong? any comments would be appreciated. EDIT : This is what I would like to achieve, I flatten a matrix to an numpy array on which I want to calculate the mean of the 6x6 neighborhood of each cell: ``` 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]])) ```