Fast 1D and 2D histogram functions in Python

Hi all,

I happen to have a (scientific) application where I need fast histograms, and I found this package very promising to be exactly what I need! If it would only support histograms in arbitrary dimensions, it would be perfect... and I think that I did in fact manage to pull it off. :)

This PR adds the function histogramdd to the package.

It is a straight generalization of the 1D and 2D code, only that the input sample is given as a tuple. The code discovers automatically how many arrays there are in the tuple and generates the output dimensions accordingly. If the dimension is D and the number of samples N, then the tuple must contain D arrays with N elements. Bins and ranges are passed as arrays as well with shapes (D,) and (D, 2), respectively to pass the number of bins as well as the ranges for each dimension.

The output of histogram2d(x, y, range, bins) should be the same as histogramdd((x, y), range, bins). The speed of the generalized function is a bit slower, so it is still beneficial to use the dedicated 1D and 2D functions where appropriate.

I'm also working on the weighted version but I thought I'd make a PR to get some feedback. This is still missing tests and documentation, and maybe someone with more C and numpy API experience could have a look if there is any possibility to make things a little bit faster still to match the performance of the 1D and 2D code.

Nominally in C code that takes some time one can release the GIL with:

Py_BEGIN_ALLOW_THREADS
<Non-Python calling code>
Py_END_ALLOW_THREADS

It looks like it would be pretty trivial to enclose your for loops with these macros, is this something you're interested in?

This PR contains one critical fix and one fix for consistency [issue #10] .

Critically, when a data point was exactly at the upper edge of all binning dimensions, then the bin_idx in histogramdd, L636 could become equal to the array size and thus would cause it to write outside of the allocated memory, causing a segfault.

I added checks whether a data point is at the uppermost border in each binning, and if so, the count will be added to the last bin, consistently with numpy's behavior. This does cause a little drop in performance and I tried various different ways of implementing this (if statements, triple operators, subtraction of ix/nx) and found that adding simple if statements cost the least in the end. The speedup is ~80 % of what it used to be on my machine vs. numpy.

This is the benchmark result on my machine without the checks:

And this is what I am getting now with the check:

In my mind this would still be an acceptable performance, but we could also choose to break with consistency with numpy here. In that case, I would fix the issue with histogramdd by changing one > to an >= when I check whether the point is in range.

First, thanks for the great library. I suspect I am hitting a memory issue. Unless I explictly copy the arrays before passing them to histogram2D, the results differ from numpy. So far, I only saw this happen when I am using the weighted version of the 2Dhistogram. Code and output follows

#! /usr/bin/env python
import numpy as np
import fast_histogram as ft

print(f"Numpy version {np.__version__}")
print(f"Fast-histogram version {ft.__version__}")

# Initial data
#######################################################
nbins=100
max_val=17.64
min_val=0.0
set_range=[[min_val, max_val],[min_val, max_val]]
with np.load("data.npz") as d:
    pot = d["pot"]
    displ_new = d["displ_new"]
    bond_ind = d["bond_ind"]


# Without copy
print("NOT EXPLICITLY COPYING THE DATA")
pot_array = pot[:, bond_ind[0,0], bond_ind[0,1]]
arr1 = displ_new[:,bond_ind[0,0]]
arr2 = displ_new[:,bond_ind[0,1]]

#Run the histograms
h_np = np.histogram2d(arr1, arr2, nbins, range=set_range, weights=pot_array)[0]
h_ft = ft.histogram2d(arr1, arr2, nbins, range=set_range, weights=pot_array)

print("===FAST-HIST===")
print(h_ft)
print("===NUMPY===")
print(h_np)


# With copy
print("EXPLICITLY COPYING THE DATA")
pot_array = pot[:, bond_ind[0,0], bond_ind[0,1]].copy()
arr1 = displ_new[:,bond_ind[0,0]].copy()
arr2 = displ_new[:,bond_ind[0,1]].copy()

#Run the histograms
h_np = np.histogram2d(arr1, arr2, nbins, range=set_range, weights=pot_array)[0]
h_ft = ft.histogram2d(arr1, arr2, nbins, range=set_range, weights=pot_array)

print("===FAST-HIST===")
print(h_ft)
print("===NUMPY===")
print(h_np)

with output

Numpy version 1.18.1
Fast-histogram version 0.8
NOT EXPLICITLY COPYING THE DATA
===FAST-HIST===
[[0.         0.         0.         ... 0.         0.         0.        ]
 [0.0425809  0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.01108671 ... 0.         0.         0.        ]
 ...
 [0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.         ... 0.         0.         0.        ]]
===NUMPY===
[[0.         0.         0.         ... 0.         0.         0.        ]
 [0.00524697 0.11627543 0.00218829 ... 0.         0.         0.        ]
 [0.         0.0233833  0.08906353 ... 0.         0.         0.        ]
 ...
 [0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.         ... 0.         0.         0.        ]]
EXPLICITLY COPYING THE DATA
===FAST-HIST===
[[0.         0.         0.         ... 0.         0.         0.        ]
 [0.00524697 0.11627543 0.00218829 ... 0.         0.         0.        ]
 [0.         0.0233833  0.08906353 ... 0.         0.         0.        ]
 ...
 [0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.         ... 0.         0.         0.        ]]
===NUMPY===
[[0.         0.         0.         ... 0.         0.         0.        ]
 [0.00524697 0.11627543 0.00218829 ... 0.         0.         0.        ]
 [0.         0.0233833  0.08906353 ... 0.         0.         0.        ]
 ...
 [0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.         ... 0.         0.         0.        ]]

The data file necessary to run the code is attached. Due to github restrictions, I had to zip it first. data.npz.zip

Hi,

I'd be interested in returning the sum of weights squared for histograms. I took a stab at implementing the 1D case. If this is something you're not interested in, no hard feelings; if you think the project has a place for this, I'd be happy to implement unit tests and a 2D version.

Thanks, Doug

Data file : https://ufile.io/cnj9l

Python Code:

import numpy as np import matplotlib.pyplot as plt from fast_histogram import histogram1d

data = np.load('x.npz') h_np, _ = np.histogram(data['x'], bins=1100, range=[0, 1100] ) h_fast = histogram1d(data['x'], bins=1100, range=[0, 1100] )

plt.plot(h_np[:-1], 'r--', label='numpy') plt.plot(h_fast[:-1], label='fast') plt.legend()

The result is as :

Also in histogram2d, the spikes are there.

• Adds npy_float32 functions. Casting determination done by a casting dictionary, fast_histogram._cast_to which mimics numpy.histogram.
• Releases the GIL in pure-C operation, although this does not appear to beneficially impact performance.

Sorry for the autistic changing of double -> npy_float64. I can revert that if you want.

Do you have any intention of upstreaming this? I feel like it would be really useful (together with having a bincount2d which is what I'm usually missing).

In setup.py the build-time dependency on NumPy is defined via a setup_requires without an install_requires. This is poorly supported, and leads to the following problem on a clean python install

pip3 install fast-histogram
python3 -c 'import fast_histogram'
ImportError: No module named 'numpy.core._multiarray_umath'

While this might be an issue with Numpy (see numpy/numpy#6599), SciPy solved the problem via scipy/scipy#3566.

Hi, I would like to try your package, but my code keeps throwing errors. First, I do not use range parameter at all (I have no need). That is OK, I can add a few numbers. However, when I try to run it again, I get The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() since I use bins=[bin_x,bin_y].

I can in general rewrite my code but my sets of bins are quite sophisticated and it would be a bit painful. Is it possible to update your code somehow?

Sincerely, V.

The current version of setup.py imports numpy to get its include path. I'm finding this to be a problem when you add fast-histogram as a dependency in your requirements.txt. Even if you also add numpy as a requirement, there is no way to enforce that numpy will have been installed by the time the setup.py in fast-histogram gets run.

Hi, Thomas

First of all, thank you for the package, it's really handy. However, I have a particular problem that might require an additional modification of the code. I have to calculate a lot of histograms (starting from the trivial case with one bin only and up to probably 10^6 bins) using the same grid but passing different sets of data. That means constant allocations of memory whereas it is possible to do this once and just keep clearing the array.

Is it potentially possible to implement this?

V.

At the moment, the test_1d_compare_with_numpy and test_2d_compare_with_numpy tests are slightly hacky in how they generate test cases. Here's the code for the 1-d case:

https://github.com/astrofrog/fast-histogram/blob/4dbb898d1ad221587262cfb609ef5ae93bda7715/fast_histogram/tests/test_histogram.py#L17-L39

What I'm trying to do is generate two arrays x and w which have the same dtype (either 32-bit or 64-bit floats, big or little endian), have the same 1-d size (sampled between 0 and 200), and have values in the range -1000, 1000. So at the moment I generate a 64-bit array, cast it inside the test, then split it into two. It would be cleaner to be able to directly generate the two arrays directly with the correct dtype, but I can't figure out how to do this.

@Zac-HD - do you have any suggestions how I might be able to achieve this in a cleaner way?

I'm experimenting with various optimizations in the pulsar search code in HENDRICS (https://github.com/StingraySoftware/HENDRICS/pull/65). One of my bottlenecks was a 2d histogram operation, and I gave a shot to yours. Tested with simulated data, everything worked and histogram2d was not a bottleneck anymore. Yay! However... When I used it on real data, Jupyter notebook was giving "dead kernel" messages with no information whatsoever. After spending one day banging my head, I tried to cast the longdouble arrays to double, and everything worked again.

It would help to add a type check, raising a ValueError or something similar if the number type is unsupported.

Hey. I was just playing around with this and was trying to see if there's a way to implement this efficiently with standard libs.

My usual way to do things like this is using the scipy.sparse.coo_matrix construct.

import scipy.sparse as sp

def bincount2d(x, y, bins):
	return sp.coo_matrix((np.ones(x.shape[0]), (x, y)), shape=(bins, bins), dtype=np.int)

If the data was scaled so that making it ints would put it in the right bins, this would work.

import numpy as np
x = np.random.random(10_000_000)
y = np.random.random(10_000_000)

from fast_histogram import histogram2d
%timeit _ = histogram2d(x, y, range=[[0, 1], [0, 1]], bins=30)

36.8 ms ± 4.14 ms per loop

xx = (x * 30)
yy = (y * 30)

%timeit bincount2d(xx, yy, bins=30)

153 ms ± 4.04 ms per loop

So your code would "only" be 5x faster, so it's about a 4x speedup over numpy. Unfortunately I cheated and didn't include shifting ``xx` so that the data aligns with the bins. I don't think it's possible to make this work without copying the data at least once, which is why I'm giving up on this route.

Thought it might be of interest, though ...

I think it should be possible to rewrite the loops to include OpenMP pragmas. The final update could be protected with #pragma omp atomic and you can use default(shared) private(tx ty ix iy) firstprivate(dataptr) or something similar. You'll need to redo the loops, however, you'll need for loops instead of the inner while (or maybe the outer do while), and the stride part might need to be precalculated in some way.