Векторизация рулет

Question

Сейчас я делаю выборку из массива и использую логическую маску для подсчета количества значений, которые лежат в пределах определенного радиуса от заданного центроида. Мне было интересно, может ли кто-нибудь показать мне, как векторизовать мой подход, вместо того, чтобы использовать цикл python для зацикливания элементов списка.

У меня есть массив 1d numpy, называемый центрами, где каждый элемент центров представляет собой центроид, представленный в виде массива 1d numpy. Я также создал логическую маску, которая центрируется по центру массива ('centre_Y', 'centre_X' в приведенном ниже коде), которую я перемещаю из середины массива в позицию центроида. Затем я делаю маску разреженной, потому что мой фактический массив, из которого я выбираю ('item' в коде ниже), является разреженным. Затем я подсчитываю количество ненулевых элементов под маской. Ниже мой код

for item in data:

    ### Do stuff

    for radius in radii:

        ### Do stuff

        # roll mask to centroid and count the number of elements within radius
        for centroid in centres:

            # roll in the vertical direction to centroid y coordinate
            mask_roll_y = np.roll(mask,centroid[0]-centre_Y,axis=0)

            # roll in the horizontal direction to centroid x coordinate and make sparse
            roll_mask = sparse.csr_matrix(np.roll(mask_roll_y,centroid[1]-centre_X,axis=1))

            # apply sparse mask to sparse data and count number of nonzero elements below
            number_count.append(np.count_nonzero(item[roll_mask]))

Теперь приведенный выше код прекрасно работает и дает желаемые результаты. Моя проблема заключается в том, что после некоторого времени цикл «для центроида в центрах» вычисляется примерно за 0,4 секунды (для одного массива данных при использовании 50 радиусов и взятии 100 выборок для каждого радиуса), что, несомненно, является самым трудоемкая часть моего кода. Мне нужно сделать это для примерно 100 наборов данных с примерно 1000 массивами в каждом наборе данных, и я хотел бы взять гораздо больше радиусов и выборок, чем я делал ранее. Итак, мне было интересно, как я могу устранить петлю «для центроида в центрах»? Я попробовал следующий ужасный кусок кода

number_count.append(np.count_nonzero(item[sparse.csr_matrix(np.roll(np.roll(mask,centres[:][0]-centre_Y,axis=0),centres[:][1]-centre_X,axis=1))]))

, где я пытался векторизовать центры, но это просто дало мне список нулей в числе чисел длины len (центры). Может ли кто-нибудь помочь, пожалуйста?

EDIT:

Я забыл указать, что мне нужно повернуть маску от центра к определенным позициям, чтобы поддерживать периодичность маски. Применение маски сразу к краям логического массива не оборачивает маску параллельным краем, и я не хочу, чтобы маска имела обрезание. Это потому, что данные, к которым я применяю маску, получены в результате моделирования с периодическими граничными условиями.

UPDATE:

Так что мне удалось написать кусок кода с использованием Numba, который избавил от необходимости использовать маски и довольно быстро. Он принимает двумерный массив, список позиций (центроидов) для выборки и радиус для выборки. Код перебирает позиции и ищет ненулевые элементы на расстоянии радиуса от каждой позиции. Границы обрабатываются периодически. Некоторые аспекты кода являются избыточными, и я уверен, что его можно сделать лучше и более общим, но он работает для того, что мне нужно. Спасибо тем, кто внес свой вклад в проблему.

@nb.njit(fastmath=True)
def nonzero_count(array,positions,radius):

    stored_values = list()
    y,x = array.shape

    for k in range(len(positions)):

        pos_y,pos_x = positions[k][0],positions[k][1]
        nonzero = 0

        # this section handles the periodic boundary conditions
        if pos_y+radius+1>y:
            # recenter pos_y so it is given a negative value
            aa = y-(pos_y+radius+1)
            # iterate around new pos_y, from -'ve to +'ve
            yy = (aa-radius,aa+radius+1)
        else:
            aa = pos_y
            yy = (pos_y-radius,pos_y+radius+1)

        if pos_x+radius+1>x:
            # recenter pos_x so it is given a negative value
            bb = x-(pos_x+radius+1)
            # iterate around new pos_x, from -'ve to +'ve
            xx = (bb-radius,bb+radius+1)
        else:
            bb = pos_x
            xx = (pos_x-radius,pos_x+radius+1)

        # this section handles the count
        for i in range(yy[0],yy[1]):
            for j in range(xx[0],xx[1]):
                # check nonzero elements lie within radius
                if array[i,j] != 0 and (bb-j)**2+(aa-i)**2<=radius**2:
                    nonzero += 1

        stored_values.append(nonzero)

    return stored_values

Answer 1

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

import numba as nb

@nb.njit(nopython=True)
def slow_code():
    # roll in the vertical direction to centroid y coordinate
    mask_roll_y = np.roll(mask,centroid[0]-centre_Y,axis=0)

    # roll in the horizontal direction to centroid x coordinate and make sparse
    roll_mask = sparse.csr_matrix(np.roll(mask_roll_y,centroid[1]-centre_X,axis=1))

    # apply sparse mask to sparse data and count number of nonzero elements below
    number_count.append(np.count_nonzero(item[roll_mask]))

for item in data:

    ### Do stuff

    for radius in radii:

        ### Do stuff

        # roll mask to centroid and count the number of elements within radius

        for centroid in centres:
            slow_code()

Answer 2

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

Наблюдать:

import numpy as np
import scipy.sparse as sparse

zero_matrix = np.zeros((1024, 1024))
%timeit np.roll(zero_matrix, (10, 10))
>>1000 loops, best of 3: 1.36 ms per loop

sparse_matrix = sparse.random(1024, 1024, density = 0.001)
%timeit np.roll(sparse_matrix, (10, 10))
>>100000 loops, best of 3: 15.4 µs per loop