Краткий вопрос

У меня есть большое изображение размером 10000x10000 элементов, которое я делю на несколько сотен различных секторов/корзин. Затем мне нужно выполнить некоторые итерационные вычисления для значений, содержащихся в каждом бине.

Как извлечь индексы из каждой ячейки, чтобы эффективно выполнить мой расчет с использованием значений ячейки?

То, что я ищу, - это решение, которое позволяет избежать узкого места необходимости выбирать каждый раз ind == j из моего большого массива. Есть ли способ получить сразу, за один раз, индексы элементов, принадлежащих каждой корзине?

Подробное объяснение

1. Простое решение

Одним из способов достижения того, что мне нужно, является использование кода, подобного следующему (см., Например, ЭТОТ связанный ответ), где я оцифровываю свои значения, а затем выполняю j-цикл, выбирая оцифрованные индексы, равные j, как показано ниже

import numpy as np

# This function func() is just a placemark for a much more complicated function.
# I am aware that my problem could be easily sped up in the specific case of
# of the sum() function, but I am looking for a general solution to the problem.
def func(x):
    y = np.sum(x)
    return y

vals = np.random.random(1e8)
nbins = 100
bins = np.linspace(0, 1, nbins+1)
ind = np.digitize(vals, bins)

result = [func(vals[ind == j]) for j in range(1, nbins)]

Это не то, что я хочу, так как он выбирает каждый раз ind == j из моего большого массива. Это делает это решение очень неэффективным и медленным.

2. Использование binned_statistics

Вышеописанный подход оказывается тем же, реализованным в scipy.stats.binned_statistic, для общего случая пользовательской функции. Используя Scipy напрямую, можно получить идентичный вывод с помощью следующего

import numpy as np
from scipy.stats import binned_statistics

vals = np.random.random(1e8)
results = binned_statistic(vals, vals, statistic=func, bins=100, range=[0, 1])[0]

3. Использование labeleled_comprehension

Другой альтернативой Scipy является использование scipy.ndimage.measurements.labeled_comprehension. Используя эту функцию, приведенный выше пример станет

import numpy as np
from scipy.ndimage import labeled_comprehension

vals = np.random.random(1e8)
nbins = 100
bins = np.linspace(0, 1, nbins+1)
ind = np.digitize(vals, bins)

result = labeled_comprehension(vals, ind, np.arange(1, nbins), func, float, 0)

К сожалению, эта форма также неэффективна и, в частности, она не имеет преимущества в скорости перед моим исходным примером.

4. Сравнение с языком IDL

Чтобы уточнить, я ищу функциональность, эквивалентную ключевому слову REVERSE_INDICES в функции HISTOGRAM языка IDL ЗДЕСЬ. Может ли эта очень полезная функциональность эффективно реплицироваться в Python?

В частности, с использованием языка IDL приведенный выше пример может быть записан как

vals = randomu(s, 1e8)
nbins = 100
bins = [0:1:1./nbins]
h = histogram(vals, MIN=bins[0], MAX=bins[-2], NBINS=nbins, REVERSE_INDICES=r)
result = dblarr(nbins)

for j=0, nbins-1 do begin
    jbins = r[r[j]:r[j+1]-1]  ; Selects indices of bin j
    result[j] = func(vals[jbins])
endfor

Вышеупомянутая реализация IDL примерно в 10 раз быстрее, чем Numpy, из-за того, что индексы бинов не нужно выбирать для каждого бина. И разница в скорости в пользу реализации IDL увеличивается с количеством бинов.