Python - vettorizzazione di una finestra scorrevole

Question

Sto provando a vettorizzare un'operazione a finestra scorrevole. Per il caso 1-d un esempio utile potrebbe andare lungo le linee di:

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,:]) 

Il n + 1 valore per ogni valore di corrente per gli indici < 5. Ma ottengo questo errore:

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 

curiosamente non vorrei ottenere questo errore per il valore n-1 che significherebbe indici inferiori a 0. non sembra in mente:

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]] 

c'è comunque intorno a questo? il mio approccio è totalmente sbagliato? eventuali commenti sarebbero apprezzati.

EDIT:

Questo è quello che vorrei realizzare, ho appiattire una matrice a una matrice NumPy su cui voglio calcolare la media del quartiere 6x6 di ogni cella:

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]])) 

Answer 1

Se capisco correttamente il problema, si desidera prendere la media di tutti i numeri di 1 livello attorno all'indice, trascurando l'indice.

ho patchato la funzione di lavorare, credo che stavano andando per qualcosa di simile:

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]])) 

ho riscritto questo usando usando affettare e viste:

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 

Utilizzando np.empty_like e poi riempire i bordi sembravano leggermente più veloci dello np.zeros_like. Per prima cosa ricontrollare danno la stessa cosa usando l'array 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. ]] 

Alcune timing:

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 

Circa ~ 1500x aumento di velocità.

sembra un buon posto per usare 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) 

consente ora di confrontare contro tutti i metodi presentati.

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 

velocità Un 4x con numba up è abbastanza nominale che indica che il codice NumPy è circa buono come la sua intenzione di ottenere. Ho tirato gli altri codici come presentati, anche se ho dovuto cambiare la risposta di @ CTZhu per includere diverse dimensioni di array.

Answer 2

Il problema si trova in x[1,x[0,:]+1], l'indice per il 2o asse: x[0,:]+1 è [1 2 3 4 5 6 7 8 9 10], in cui l'indice 10 è più grande della dimensione di x.

Nel caso della x[1,x[0,:]-1], l'indice del 2 ° asse è [-1 0 1 2 3 4 5 6 7 8 9], si finisce per ottenere [9 0 1 2 3 4 5 6 7 8], come 9 è l'ultimo elemento e ha un indice di -1. L'indice del secondo elemento dalla fine è -2 e così via.

Con np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]-1],x[1,:]) e x[0,:]=[0 1 2 3 4 5 6 7 8 9], ciò essenzialmente sta succedendo è che la prima cella viene acquistata forma x[1,:] perché x[0,0] è 0 e x[0,:]<5)&(x[0,:]>0 è False. I successivi quattro elementi sono presi da x[1,x[0,:]-1]. Il resto è da x[1,:]. Infine, il risultato è [0 0 1 2 3 4 5 6 7 8]

Può sembrare essere OK per finestra scorrevole di appena 1 cellulare, ma è gonna a sorpresa 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]) 

Quando si tenta di spostarlo da una finestra di due celle .

Per questo specifico problema, se vogliamo mantenere ogni cosa in una linea, questa, farà:

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

[0 0 1 2 3 5 6 7 8 9] 
[0 0 0 1 2 5 6 7 8 9] 
[0 0 0 0 1 5 6 7 8 9] 
[0 0 0 0 0 5 6 7 8 9] 
[0 0 0 0 0 5 6 7 8 9] 
[0 0 0 0 0 5 6 7 8 9] 

Edit: Ora ho capito la tua domanda iniziale meglio, in pratica si vuole prendere un 2D array e calcolare la media delle celle N * N attorno a ciascuna cella. Questo è abbastanza comune. Per prima cosa probabilmente si vuole limitare N a numeri dispari, altrimenti una cosa come la media 2 * 2 attorno a una cella è difficile da definire. Supponiamo di voler 3 * 3 media:

#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]]) 

penso che non c'è bisogno di appiattire voi 2D-array, che provoca confusione. Inoltre, se si desidera gestire gli elementi di bordo in modo diverso da quelli appena tagliati, considerare la possibilità di creare array mascherati utilizzando np.ma in Passaggio "Sposta l'array originale".

Answer 3

Sembra che tu stia cercando di calcolare una convoluzione 2D. Se siete in grado di utilizzare scipy, vorrei suggerire cercando 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 ragione per questo calcola la media di tutti i quartieri 3x3 si può vedere se si "srotolare" convolve2d nei suoi cicli costituenti. Effettivamente (e ignorando quello che succede ai bordi della fonte e gli array del kernel), è calcolo:

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] 

Quindi, se tutti i valori nel vostro kernel è 1/(1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1) == 1/9, si può riscrivere il codice di cui sopra come:

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() 

che è esattamente la stessa di calcolare la media dei valori di Matriz, su un'area di 3x3, a partire da i, j.

Un vantaggio di fare le cose in questo modo è che è possibile modificare facilmente i pesi associati al proprio vicinato impostando opportunamente i valori nel proprio kernel. Così, per esempio, se si voleva dare al valore centrale in ogni quartiere due volte tanto peso come gli altri, si potrebbe costruire il kernel come questo:

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

e il codice di convoluzione dovrebbe rimanere lo stesso, ma il il calcolo darebbe un diverso tipo di media (uno "ponderato al centro"). Ci sono molte possibilità qui; speriamo che questo fornisca una bella astrazione per il compito che stai facendo.

Answer 4

Esiste solo una funzione nella libreria standard di Scipy che calcola la media delle finestre scorrevoli estremamente veloce. Si chiama uniform_filter. Si può usare per implementare la funzione media-di-quartiere come segue:

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)) 

Ciò restituisce un array X dove X[i,j] è la media di tutti i vicini di i,j in arr escluso i,j stessa. Si noti che la prima e l'ultima colonna e la prima e l'ultima riga sono soggette a condizioni al contorno e pertanto potrebbero non essere valide per l'applicazione (è possibile utilizzare mode= per controllare la regola dei limiti, se necessario).

Poiché uniform_filter utilizza un algoritmo lineare tempo molto efficace implementato in rettilineo C (lineare solo nella dimensione arr), occorre facilmente superare qualsiasi altro soluzioni, specialmente quando win è grande.