Qual è il modo più veloce per confrontare le patch di un array?

Question

Desidero confrontare le diverse aree di un array bidimensionale $ A $ con un array specificato $ b $ di dimensioni inferiori. Dato che devo farlo molte volte, è fondamentale che questo sia eseguito molto velocemente. Ho una soluzione che funziona bene, ma speravo che potesse essere fatto più bello e più veloce.

In dettaglio:

Diciamo che abbiamo un grande array e un piccolo array. I loop su tutte le possibili "patch" all'interno del grande array che hanno le stesse dimensioni del piccolo array e confrontare queste patch con il piccolo array dato.

def get_best_fit(big_array, small_array): 

    # we assume the small array is square 
    patch_size = small_array.shape[0] 
    min_value = np.inf 
    for x in range(patch_size, big_array.shape[0] - patch_size): 
     for y in range(patch_size, big_array.shape[1] - patch_size): 
      p = get_patch_term(x, y, patch_size, big_array) 
      tmp = some_metric(p, small_array) 
      if min_value > tmp: 
       min_value = tmp 
       min_patch = p 

    return min_patch, min_value 

Al fine di ottenere le patch ho avuto questa implementazione di accesso agli array diretta:

def get_patch_term(x, y, patch_size, data): 
    """ 
    a patch has the size (patch_size)^^2 
    """ 
    patch = data[(x - (patch_size-1)/2): (x + (patch_size-1)/2 + 1), 
       (y - (patch_size-1)/2): (y + (patch_size-1)/2 + 1)] 
    return patch 

immagino che questo è il compito più importante e può essere eseguita più veloce, ma non sono sicuro su di esso.

Ho dato un'occhiata a Cython ma forse ho sbagliato, non ho molta familiarità con esso.

Il mio primo tentativo è stato una traduzione diretta in Cython:

def get_patch_term_fast(Py_ssize_t x_i, Py_ssize_t y_i, 
         Py_ssize_t patch_size, 
         np.ndarray[DTYPE_t, ndim=2] big_array): 

    assert big_array.dtype == DTYPE 
    patch_size = (patch_size - 1)/2 

    cdef np.ndarray[DTYPE_t, ndim=2] patch = <np.ndarray[DTYPE_t, ndim=2]>big_array[(x_i - patch_size):(x_i + patch_size + 1), (y_i - patch_size): (y_i + patch_size + 1)] 
    return patch 

E questo sembra essere più veloce (vedi sotto), ma ho pensato che un approccio parallelo dovrebbe essere migliore, così mi si avvicinò con questo

def get_patch_term_fast_parallel(Py_ssize_t x_i, Py_ssize_t y_i, 
           Py_ssize_t patch_size, 
           np.ndarray[DTYPE_t, ndim=2] big_array): 

    assert big_array.dtype == DTYPE 
    patch_size = (patch_size - 1)/2 

    assert big_array.dtype == DTYPE 
    cdef Py_ssize_t x 
    cdef Py_ssize_t y 


    cdef np.ndarray[DTYPE_t, ndim=1] patch = np.empty(np.power((2 * patch_size) + 1, 2)) 
    with nogil, parallel(): 
     for x in prange(x_i - patch_size, x_i + patch_size + 1): 
      for y in prange(y_i - patch_size, y_i + patch_size + 1): 
       patch[((x - (x_i - patch_size)) * (2 * patch_size + 1)) + (y - (y_i - patch_size))] = big_array[x, y] 
    #cdef np.ndarray[DTYPE_t, ndim=2] patch = <np.ndarray[DTYPE_t, ndim=2]>big_array[(x_i - patch_size):(x_i + patch_size + 1), (y_i - patch_size): (y_i + patch_size + 1)] 
    return patch 

Che, sfortunatamente, non è più veloce. Per la prova ho usato:

A = np.array(range(1200), dtype=np.float).reshape(30, 40) 
b = np.array([41, 42, 81, 84]).reshape(2, 2) 

x = 7 
y = 7 
print(timeit.timeit(lambda:get_patch_term_fast(x, y, 3, A), number=300)) 
print(timeit.timeit(lambda:get_patch_term_fast_parallel(x, y, 3, A).reshape(3,3), number=300)) 
print(timeit.timeit(lambda:get_patch_term(x, y, 3, A), number=300)) 

che dà

0.0008792859734967351 
0.0029909340664744377 
0.0029337930027395487 

Quindi, la mia prima domanda è, è possibile farlo più velocemente? La seconda domanda sarebbe: perché l'approccio parallelo non è più veloce rispetto all'implementazione numpy originale?

Edit:

ho cercato di parallelizzare ulteriormente la funzione get_best_fit ma purtroppo ho un sacco di errori affermando che non posso assegnare un oggetto Python, senza gil.

Ecco il codice:

def get_best_fit_fast(np.ndarray[DTYPE_t, ndim=2] big_array, 
         np.ndarray[DTYPE_t, ndim=2] small_array): 

    # we assume the small array is square 
    cdef Py_ssize_t patch_size = small_array.shape[0] 

    cdef Py_ssize_t x 
    cdef Py_ssize_t y 

    cdef Py_ssize_t x_range = big_array.shape[0] - patch_size 
    cdef Py_ssize_t y_range = big_array.shape[1] - patch_size 

    cdef np.ndarray[DTYPE_t, ndim=2] p 
    cdef np.ndarray[DTYPE_t, ndim=2] weights = np.empty((x_range - patch_size)*(y_range - patch_size)).reshape((x_range - patch_size), (y_range - patch_size)) 

    with nogil, parallel(): 
     for x in prange(patch_size, x_range): 
      for y in prange(patch_size, y_range): 
       p = get_patch_term_fast(x, y, patch_size, big_array) 
       weights[x - patch_size, y - patch_size] = np.linalg.norm(np.abs(p - small_array)) 

    return np.min(weights) 

PS: ho omesso la parte di restituire il più piccolo cerotto ...

Answer 1

initialy inviato un'altra risposta sulla base di pattern matching (portato via dal titolo), inviare un'altra risposta

uso un ciclo per ciclo su entrambi gli array (grandi e piccole) e applicare metriche di correlazione parziale (o altro) senza liste affettare tutto il tempo:

come sidenote, osservate il fatto (ovviamente a seconda della metrica) che una volta che una patch è stata completata, la patch successiva (una riga in basso o una colonna a destra) condivide molto con la patch precedente, solo le righe iniziale e finale (o colonne) cambiano in tal modo la nuova correlazione può essere calcolata più velocemente dalla correlazione precedente sottraendo la riga precedente e aggiungendo una nuova riga (o viceversa). Poiché la funzione metrica non è data, viene aggiunta come osservazione.

def get_best_fit(big_array, small_array): 

    # we assume the small array is square 
    patch_size = small_array.shape[0] 
    x = 0 
    y = 0 
    x2 = 0 
    y2 = 0 
    W = big_array.shape[0] 
    H = big_array.shape[1] 
    W2 = patch_size 
    H2 = patch_size 
    min_value = np.inf 
    tmp = 0 
    min_patch = None 
    start = 0 
    end = H*W 
    start2 = 0 
    end2 = W2*H2 
    while start < end: 

     tmp = 0 
     start2 = 0 
     x2 = 0 
     y2 = 0 
     valid = True 
     while start2 < end2: 
      if x+x2 >= W or y+y2 >= H: 
       valid = False 
       break 
      # !!compute partial metric!! 
      # no need to slice lists all the time 
      tmp += some_metric_partial(x, y, x2, y2, big_array, small_array) 
      x2 += 1 
      if x2>=W2: 
       x2 = 0 
       y2 += 1 
      start2 += 1 

     # one patch covered 
     # try next patch 
     if valid and min_value > tmp: 
      min_value = tmp 
      min_patch = [x, y, W2, H2] 

     x += 1 
     if x>=W: 
      x = 0 
      y += 1 

     start += 1 

    return min_patch, min_value 

Answer 2

penso dipende da ciò che la vostra funzione some_metric fa, è possibile che ci sia già un'implementazione altamente ottimizzato là fuori. Ad esempio, se si tratta solo di una convoluzione, dai un'occhiata allo Theano che ti permetterà anche di sfruttare la GPU abbastanza facilmente. Anche se la tua funzione non è tanto semplice quanto una convoluzione diretta, è probabile che all'interno di Theano ci siano blocchi predefiniti che puoi utilizzare piuttosto che provare ad andare veramente a basso livello con Cython.

Answer 3

L'altro problema con la misurazione del tempo della funzione parallela è che si chiama reshape sull'oggetto array dopo aver eseguito la funzione parallela. Potrebbe essere il caso che la funzione parallela sia più veloce, ma poi la risagoma sta aggiungendo tempo extra (anche se reshape dovrebbe essere abbastanza veloce, ma non sono sicuro di quanto velocemente).

L'altro problema è che non sappiamo quale sia il termine some_metric e non sembra che lo stiate usando in parallelo. Il modo in cui vedo il tuo codice parallelo funzionante è che stai ricevendo le patch in parallelo e poi mettendole in coda per essere analizzate dalla funzione some_metric, quindi sconfiggendo lo scopo della parallelizzazione del tuo codice.

Come ha detto John Greenhall, l'utilizzo di FFT e le convoluzioni può essere abbastanza veloce. Tuttavia, per trarre vantaggio dall'elaborazione parallela, sarà comunque necessario eseguire l'analisi della patch e del piccolo array in parallelo.

Quanto sono grandi gli array? Anche qui la memoria è un problema?