Каков самый быстрый способ сравнить патчи массива?

Question

Я хочу сравнить различные области 2-мерного массива $ A $ с заданным массивом $ b $ меньшего размера. Поскольку я должен делать это много раз, очень важно, чтобы это выполнялось очень быстро. У меня есть решение, которое отлично работает, но я надеялся, что это можно сделать лучше и быстрее.

Подробно:

Допустим, у нас есть большой массив и небольшой массив. Я перебираю все возможные «патчи» в большом массиве, которые имеют тот же размер, что и малый массив, и сравнивают эти патчи с данным небольшим массивом.

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 

Для того, чтобы получить патчи я получил это непосредственное осуществление доступа к массиву:

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 

Я думаю, что это наиболее важная задача и может быть выполнена быстрее, но я не уверен в этом.

Я посмотрел на Китона, но, возможно, я сделал это неправильно, я не очень-то знаком с ним.

Моя первая попытка была прямой перевод в 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 

И это, кажется, быстрее (см ниже), но я подумал, что параллельный подход должен быть лучше, так что я пришел с этим

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 

Это, к сожалению, не быстрее. Для тестирования я использовал:

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

Что дает

0.0008792859734967351 
0.0029909340664744377 
0.0029337930027395487 

Итак, мой первый вопрос, можно ли сделать это быстрее? Второй вопрос: почему параллельный подход не быстрее первоначальной реализации numpy?

Edit:

Я пытался дальше распараллелить функцию get_best_fit, но, к сожалению, я получаю много ошибок, заявив, что я не могу назначить объект Python без Gil.

Вот код:

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: Я опустил часть возвращения наименьшую патч ...

Answer 1

initialy опубликовал еще один ответ на основе сопоставления с образцом (уносится по названию), опубликовать еще один ответ

использовать один цикл для перебора обоих массивах (большой и маленький) и применить частичную корреляционную метрику (или любой другой) без нарезки списков все время:

в качестве побочного элемента, обратите внимание на факт (в зависимости от метрики, конечно), что после завершения одного патча следующий патч (либо с одной строкой вниз, либо с одним столбцом справа) имеет много общего с предыдущим патчем, только его начальные и конечные строки (или столбцы), таким образом, новая корреляция может быть вычислена быстрее из предыдущей корреляции путем вычитания предыдущей строки и добавления новой строки (или наоборот). Поскольку метрическая функция не задана, добавляется как наблюдение.

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

Я думаю, что зависит от того, что делает ваша some_metric функция, возможно уже есть высоко оптимизированная реализация там. Например, если это просто свертка, взгляните на Theano, что даже позволит вам легко использовать GPU. Даже если ваша функция не такая простая, как простая свертка, скорее всего, в Theanano вы можете использовать строительные блоки, а не пытаться идти на очень низкий уровень с Cython.

Answer 3

Другая проблема с измерением времени параллельной функции заключается в том, что вы вызываете reshape на объект массива после запуска вашей параллельной функции. Это может быть так, что параллельная функция выполняется быстрее, но затем изменение формы добавляет дополнительное время (хотя reshape должно быть довольно быстрым, но я не уверен, насколько быстро).

Другая проблема заключается в том, что мы не знаем, каков ваш термин some_metric, и он не кажется, что вы используете это параллельно. То, как я вижу, что ваш параллельный код работает, заключается в том, что вы получаете патчи параллельно, а затем помещаете их в очередь, которую нужно проанализировать с помощью функции some_metric, следовательно, побеждая цель распараллеливания вашего кода.

Как сказал Джон Гринхолл, использование БПФ и сверток может быть довольно быстрым. Тем не менее, чтобы воспользоваться преимуществами параллельной обработки, вам все равно придется анализировать патч и малый массив параллельно.

Насколько велики массивы? Проблема памяти здесь тоже?