search

m3g / CellListMap.jl

Flexible implementation of cell lists to map the calculations of particle-pair dependent functions, such as forces, energies, neighbor lists, etc.
https://m3g.github.io/CellListMap.jl/
MIT License
87 stars 4 forks source link

Pair counts changing depending on number of threads #47

Closed epaillas closed 2 years ago

epaillas commented 2 years ago

Hi! First of all, thanks for putting up together this code. It's been pretty useful for my calculations of galaxy two-point clustering.

I have run a few tests where I compute the PDF of galaxy number densities (which is pretty similar to a histogram of particle distances), and I've found that the number of pairs that are returned by map_pairwise change depending on how many threads are used. I only get sensible results when using a single thread. I have setup up a test repository to exemplify the problem:

https://github.com/epaillas/julia-sandbox

I call Julia from Python by running call_julia.py, which will generate a plot of the galaxy density PDF. I get different results if I change the number of threads in Line 6.

I am using an M1 Macbook Pro 2020 with Big Sur, Python 3.8.8, Julia 1.7.0, and the latest version of CellListMap.

Is there any obvious problem in my code that would cause this discrepancy, or could it be something else?

lmiq commented 2 years ago

Here:

(x, y, i, j, d2, output) -> count_pairs!(i, DD),

output and DD (which name doesn't matter)

must be the same variable. This is just a function definition that for the input variables on the left returns the result of the function on the right.

On the right a variable that is not on the left only appears if it is some data that has to be used but is not provided by the standard interface (like masses, for example)

lmiq commented 2 years ago

(otherwise you are updating the DD array, which is not threaded, and you have concurrency problems)

lmiq commented 2 years ago

Also note that within the count_pairs function you may want to update both count[i] and count[j], because the function only runs once per pair. But this is not related to threading.

epaillas commented 2 years ago

Thanks a lot for the prompt reply - this indeed solved the issue, and I'm glad to know it was just a silly mistake from my side. You're also correct about updating counts[j], as I would be underestimating the pair counts otherwise.