Jaki jest najszybszy sposób porównywania poprawek z tablicy?

Question

Chcę porównać różne obszary dwuwymiarowej tablicy $ A $ z podaną tablicą $ b $ o mniejszym rozmiarze. Ponieważ muszę to robić wiele razy, ważne jest, aby było to wykonywane bardzo szybko. Mam rozwiązanie, które działa dobrze, ale miałem nadzieję, że da się to zrobić ładniej i szybciej.

W szczegółach:

Powiedzmy mamy duży wachlarz i małą tablicę. Pętlę na wszystkie możliwe "łaty" w dużej tablicy, które mają ten sam rozmiar co mała tablica i porównuję te poprawki z podaną małą tablicą.

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 

W celu uzyskania poprawki Mam ten bezpośrednią realizację dostępu do tablicy:

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 

myślę, że jest to najważniejsze zadanie i może być wykonywane szybciej, ale nie jestem pewien o tym.

Zajrzałem do Cythona, ale może zrobiłem to źle, nie jestem do końca obeznany z tym.

Moja pierwsza próba była bezpośrednim tłumaczeniem na 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 

I to wydaje się być szybszy (patrz poniżej), ale myślałem, że równoległe podejście powinno być lepiej, więc wymyśliłem tego

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 

Co niestety nie jest szybsze. Do testów użyłem:

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

Który daje

0.0008792859734967351 
0.0029909340664744377 
0.0029337930027395487 

Tak, moje pierwsze pytanie brzmi, czy to możliwe, aby zrobić to szybciej? Drugie pytanie brzmi: dlaczego równoległe podejście nie jest szybsze niż oryginalna implementacja numpy?

Edit:

próbowałem dalej parallelize funkcję get_best_fit ale niestety mam dużo błędów, stwierdzające, że nie można przypisać obiekt Pythona bez gil.

Oto kod:

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: I pominął część zwracając najmniejszej łatę ...

Answer 1

initialy pisał kolejną odpowiedź na podstawie dopasowania wzorca (ponieść tytule), zakładać inna odpowiedź

użyć jednej pętli do pętli na obu tablicach (duże i małe) i zastosowanie częściowej korelacji metrykę (lub cokolwiek) bez listach krojenie cały czas:

jako sidenote, obserwuj fakt (w zależności od metryki oczywiście), że po ukończeniu jednej łatki następna łatka (jeden wiersz w dół lub jedna kolumna w prawo) ma wiele wspólnego z poprzednią poprawką, a jedynie początkowe i końcowe wiersze (lub kolumny) zmieniają się, więc nowa korelacja może zostać obliczona szybciej od poprzedniej korelacji poprzez odjęcie poprzedniego wiersza i dodanie nowego wiersza (lub odwrotnie). Ponieważ funkcja metryczna nie jest podana, jest dodawana jako obserwacja.

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

myślę zależy od tego, co robi twoja funkcja some_metric, możliwe jest już wysoce zoptymalizowany realizacji tam. Na przykład, jeśli jest to tylko splot, należy spojrzeć na Theano, który nawet pozwoli Ci w łatwy sposób wykorzystać GPU. Nawet jeśli twoja funkcja nie jest tak prosta, jak prosta splot, najprawdopodobniej będą to klocki w Theano, których możesz użyć, zamiast próbować przejść naprawdę niski poziom z Cythonem.

Answer 3

Innym problemem związanym z pomiarem czasu funkcji równoległej jest wywołanie reshape na obiekcie macierzy po uruchomieniu funkcji równoległej. Może się zdarzyć, że funkcja równoległa jest szybsza, ale następnie zmiana kształtu dodaje dodatkowy czas (chociaż reshape powinien być dość szybki, ale nie jestem pewien, jak szybko).

Inną kwestią jest to, że nie wiemy, jaki jest twój termin some_metric i nie wydaje Ci się, że używasz go równolegle. Sposób, w jaki widzę twój równoległy kod działa, polega na tym, że otrzymujesz poprawki równolegle, a następnie umieszczasz je w kolejce do analizy za pomocą funkcji some_metric, tym samym pokonując cel zrównoleglania twojego kodu.

Jak powiedział John Greenhall, używanie FFT i zwojów może być dość szybkie. Jednakże, aby skorzystać z przetwarzania równoległego, trzeba by równolegle przeprowadzić analizę zarówno łaty, jak i małej macierzy.

Jak duże są tablice? Czy pamięć jest tu również problemem?