One generic comment: by using bincount you are forcing your weighted histogram to be convertible to double. Not sure if there is a use case for complex, or object weights, but it may be a good idea to route those to the generic, slower algorithm.

@jaimefrio - thanks for the feedback! I agree about needing to be careful about the weights type. I'm not sure if the code was originally intended to be used with non-float weights since there was the comment # Histogram is an integer or a float array depending on the weights. but I'll try and find some examples of non-float weights and will add some unit tests (and will make sure we use the slower algorithm there).

@jaimefrio - I've now made sure that if the weights are complex, we do the right thing with bincount, and if they are not complex, float, or int, we revert to the original algorithm (and I added a test for all these).

Note that there is still an issue with bincount mentioned in #6100 (comment) but I think that can be dealt with separately. It's not specific to this PR, and it's a definite corner case.

Quick question: In one of my use cases, instead of providing the number of bins, I provide an array of equally spaced bins. The obvious fix is on my end but I believe you are technically speaking not using the optimized code for the case where an array of equally spaced bins is provided. Is that correct?

Would there be problems associated doing that check as well? I.e. something along the lines of

bin_widths = np.diff(bins)
if not iterable(bins) or np.all(bin_widths == bin_widths[0])

There's some overhead for that check of course and there are probably rounding issues so the simplest way out may just be to use the function as intended.

@ctk3b, you might want the range keyword.

The issue with checking the bin widths is that it can fail due to rounding error, since a set of bins that are equally spaced in floating point doesn't have to pass that check. e.g.:

In [136]: a = np.arange(10) * (1./7)

In [137]: bins = np.arange(10) * (1./7)

In [138]: bin_widths = np.diff(bins)

In [139]: bin_widths == bin_widths[0]
Out[139]: array([ True,  True,  True,  True, False, False, False, False, False], dtype=bool)

edit: doh I see you are aware that rounding issues exist.

Yeah that's the tricky part I guess. One thing to maybe consider, if the cost of the check is acceptable, is to do a np.allclose() check and suggest to the user that they might want to consider using range and specify number of bins instead.

*but they have to be absolutely sure that's what they want

I personally don't think that we should check for equal bin spacing in the bins array (if it's not an int), but happy to add it if there is a consensus. An alternative is to simply make the docstring clearer and add a note about how to get optimal performance.

@jaimefrio - apart from that question, I think this is pretty much ready for a final review.

It mostly LGTM. Have quite a few questions/comments, but aside from the int-that-should-be-an-intp, there is nothing fundamentally wrong that I see.

@jaimefrio - I'll reply here about the blocks/chunks to make sure the comment stays here in the long term.

I copied the use of the blocks from the older/slower method. It turns out that doing this results in a factor of 2 improvement in performance. Let's take the case of a 10^8 array. The times to do the histogram (with the new method in this PR) are:

BLOCK         TIME (s)
2**16         1.15103316307 (current implementation)
2**17         1.10001802444
2**18         1.08962392807
2**19         1.10686206818
2**20         1.23000407219
2**21         1.69827389717
2**22         1.80364918709
2**23         1.91272115707
...
2**28         2.28366589546

So with BLOCK = 65536 we are basically 2x faster than with no chunking (it turns out we could do even better with BLOCK=2**18). I presume this must have something to do with larger arrays needing disproportionately more time to be allocated? In addition, the memory requirements are very different - for the 10^8 array, the peak RAM with BLOCK = 65536 is a little under 900MB, while if we don't use blocks, the memory requirement goes up to 2.7Gb.

While slightly inelegant, I do think a 2x improvement in speed and 3x improvement in memory usage is worth it. In fact, I think that having a 3x worsening in memory usage would be considered a performance regression by many, even if the code ran 5x faster. Do you agree that it makes sense to leave it as-is? I can add a comment about why this is done.

I'll buy that 2x. Thanks for takinig the time to substantiate with measurements all my existential doubts,

Correct me if I am wrong, but it seems we both agree that, eventually, the right thing to do would be to rewrite at least this particular path in C, right?

If you could fix that int-that-should-be-an-intp, and get rid of some of the extraneous white space, I'll go ahead and merge it.

@jaimefrio - yes, I think in the long term we could benefit from switching over to C and I'd be happy to look into it (but that will need to wait a few weeks). I'll finish up the current pull request now so that you can go ahead and merge it.

@jaimefrio - I've made the fixes you requested.

I'd be interested in trying to tackle the switching over to C - should I open an issue to keep track of it, or should I just open a pull request once I have something working?

@jaimefrio - I can also look into applying the same kind of optimization to the histogramdd (and by extension histogram2d).

6 commits seem a few too many: if you can squash them into a single one, we can put it in.

@jaimefrio - done!

Thanks a lot Thomas.

@jaimefrio - thanks for all your comments! Once I get a chance (reasonably soon hopefully) I'll investigate some of the further improvements discussed here, including for example switching to C for part of the code.