title: Efficiently get indices of histogram bins in Python
body: <h1>Short Question</h1>
  
  <p>I have a large 10000x10000 elements image, which I bin into a few hundred different sectors/bins. I then need to perform some iterative calculation on the values contained within each bin.</p>
  
  <p>How do I extract the indices of each bin to efficiently perform my calculation using the bins values?</p>
  
  <p>What I am looking for is a solution which avoids the bottleneck of having to select every time <code>ind == j</code> from my large array. Is there a way to obtain directly, in one go, the indices of the elements belonging to every bin?</p>
  
  <h1>Detailed Explanation</h1>
  
  <h2>1. Straightforward Solution</h2>
  
  <p>One way to achieve what I need is to use code like the following (see e.g. <a href="https://stackoverflow.com/a/6163403/3753826">THIS</a> related answer), where I digitize my values and then have a j-loop selecting digitized indices equal to j like below</p>
  
  <pre><code>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)]
  </code></pre>
  
  <p>This is not what I want as it selects every time <code>ind == j</code> from my large array. This makes this solution very inefficient and slow.</p>
  
  <h2>2. Using binned_statistics</h2>
  
  <p>The above approach turns out to be the same implemented in <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binned_statistic.html" rel="noreferrer">scipy.stats.binned_statistic</a>, for the general case of a user-defined function. Using Scipy directly an identical output can be obtained with the following</p>
  
  <pre><code>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]
  </code></pre>
  
  <h2>3. Using labeled_comprehension</h2>
  
  <p>Another Scipy alternative is to use <a href="https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.ndimage.measurements.labeled_comprehension.html" rel="noreferrer">scipy.ndimage.measurements.labeled_comprehension</a>. Using that function, the above example would become</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>Unfortunately also this form is inefficient and in particular, it has no speed advantage over my original example.</p>
  
  <h2>4. Comparison with IDL language</h2>
  
  <p>To further clarify, what I am looking for is a functionality equivalent to the <code>REVERSE_INDICES</code> keyword in the <code>HISTOGRAM</code> function of the IDL language <a href="http://www.exelisvis.com/docs/HISTOGRAM.html" rel="noreferrer">HERE</a>. Can this very useful functionality be efficiently replicated in Python?</p>
  
  <p>Specifically, using the IDL language the above example could be written as</p>
  
  <pre><code>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
  </code></pre>
  
  <p>The above IDL implementation is about 10 times faster than the Numpy one, due to the fact that the indices of the bins do not have to be selected for every bin. And the speed difference in favour of the IDL implementation increases with the number of bins.</p>