Эффективно получить индексы гистограмм в Python

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

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

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

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

Детальное объяснение

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

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

 import numpy as np # This function func() is just a place mark for a much more complicated function. # I am aware that my problem could be easily speed 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. Использование labeled_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 функции 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 увеличивается с количеством ящиков.

4 Solutions collect form web for “Эффективно получить индексы гистограмм в Python”

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

 def binned_statistic(x, values, func, nbins, range): '''The usage is approximately the same as the scipy one''' from scipy.sparse import csr_matrix N = len(values) r0, r1 = range digitized = (float(nbins) / (r1-r0) * (x-r0)).astype(int) S = csr_matrix((values, [digitized, np.arange(N)]), shape=(nbins, N)) return [func(group) for group in np.split(S.data, S.indptr[1:-1])] 

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

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

 vals = np.random.random(1e4) nbins = 100 bins = np.linspace(0, 1, nbins+1) ind = np.digitize(vals, bins) new_order = argsort(ind) ind = ind[new_order] ordered_vals = vals[new_order] # slower way of calculating first_hit (first version of this post) # _,first_hit = unique(ind,return_index=True) # faster way: first_hit = searchsorted(ind,arange(1,nbins-1)) first_hit.sort() #example of using the data: for j in range(nbins-1): #I am using a plotting function for your f, to show that they cluster plot(ordered_vals[first_hit[j]:first_hit[j+1]],'o') 

На рисунке показано, что ящики на самом деле являются кластерами, как ожидалось: введите описание изображения здесь

Вы можете вдвое сократить время вычисления, сначала отсортировав массив, затем используйте np.searchsorted .

 vals = np.random.random(1e8) vals.sort() nbins = 100 bins = np.linspace(0, 1, nbins+1) ind = np.digitize(vals, bins) results = [func(vals[np.searchsorted(ind,j,side='left'): np.searchsorted(ind,j,side='right')]) for j in range(1,nbins)] 

Используя 1e8 качестве тестового примера, я 1e8 с 34 секунд вычисления до примерно 17.

Одним из эффективных решений является использование пакета numpy_indexed (отказ от ответственности: я являюсь его автором):

 import numpy_indexed as npi npi.group_by(ind).split(vals) 
  • Разбор и вычисление определений логических множеств
  • коэффициенты интерполяции сплайна в scipy
  • База данных или табличное решение для временных массивов Numpy
  • matplotlib.mlab.griddata очень медленный и возвращает массив nan при вводе действительных данных
  • эффективно находить интервал с не-нулями в scipy / numpy в Python?
  • Java Scientific Packages похож на SciPy?
  • Scipy Отрицательное расстояние? Какие?
  • Оценка хроматической аберрации в питоне
  • Python - лучший язык программирования в мире.