BUG DDrppi buggy for cross-corr in periodic spaces #210
Description
General information
- Corrfunc version: 2.3.2
- platform: osx
- installation method (pip/source/other?): conda
Issue description
I am computing the pairs by hand and not getting what corrfunc reports.
Expected behavior
It should match the brute force computation.
Actual behavior
It doesn't match it.
What have you tried so far?
A bunch of stuff to try and see if it will miss one particle. It does not, so the bug is more subtle. It does not appear to be a bin edge thing.
Minimal failing example
import numpy as np
import Corrfunc
lbox = 120
zmax = 60
n_pi = int(zmax)
rmin = 40
rmax = 45
n_r = 2
rbins = np.array([0, rmin, rmax])
lbox_2 = lbox / 2
npts = 11
rng = np.random.RandomState(seed=42)
x1 = rng.uniform(low=0, high=lbox, size=npts)
y1 = rng.uniform(low=0, high=lbox, size=npts)
z1 = rng.uniform(low=0, high=lbox, size=npts)
w1 = rng.uniform(size=npts) * 0 + 1
nparts = 100
px1 = rng.uniform(low=0, high=lbox, size=nparts)
py1 = rng.uniform(low=0, high=lbox, size=nparts)
pz1 = rng.uniform(low=0, high=lbox, size=nparts)
print("\n")
dd_cf = Corrfunc.theory.DDrppi(
False,
1,
zmax,
rbins,
x1,
y1,
z1,
weights1=w1,
X2=px1,
Y2=py1,
Z2=pz1,
weights2=np.ones_like(px1),
periodic=True,
boxsize=lbox,
weight_type="pair_product",
verbose=True,
isa="fallback",
)
dd_cf = dd_cf["weightavg"].reshape((n_r, n_pi)) * dd_cf["npairs"].reshape(
(n_r, n_pi)
)
print("\n")
print("testing particles one by one...")
dd = np.zeros(n_r)
for hind, (hx, hy, hz) in enumerate(zip(x1, y1, z1)):
for px, py, pz in zip(px1, py1, pz1):
dx = np.abs(hx - px)
dy = np.abs(hy - py)
dz = np.abs(hz - pz)
if dx > lbox_2:
dx = lbox - dx
if dy > lbox_2:
dy = lbox - dy
if dz > lbox_2:
dz = lbox - dz
if dz <= zmax:
r = np.sqrt(dx * dx + dy * dy)
if r < rmin:
dd[0] += w1[hind]
if dd_cf[0, int(dz)] == 0:
print(
" def missed part:",
r,
dz,
int(dz),
dd_cf[1, int(dz)],
)
elif r >= rmin and r < rmax:
dd[1] += w1[hind]
if dd_cf[1, int(dz)] == 0:
print(
" def missed part:",
r,
dz,
int(dz),
dd_cf[1, int(dz)],
)
print("dd corrfunc:", np.sum(dd_cf, axis=1))
print("dd brute :", dd)
dd_cf_sum = np.sum(dd_cf, axis=1)
assert np.array_equal(dd_cf_sum, dd)This script outputs
(anl) clarence:Desktop beckermr$ python corrfunc_test.py
[Warning] The CPU supports AVX2 but the compiler does not. Can you try another compiler?
ND1 = 11 [xmin,ymin,zmin] = [2.470139,16.739263,5.574050], [xmax,ymax,zmax] = [114.085717,116.389182,94.221115]
ND2 = 100 [xmin,ymin,zmin] = [0.662654,0.607390,1.727219], [xmax,ymax,zmax] = [118.426432,118.278054,118.806462]
Running with points in [xmin,xmax] = 0.662654,118.426432 with periodic wrapping = 120.000000
Running with points in [ymin,ymax] = 0.607390,118.278054 with periodic wrapping = 120.000000
Running with points in [zmin,zmax] = 1.727219,118.806462 with periodic wrapping = 120.000000
In gridlink_double> Running with [nmesh_x, nmesh_y, nmesh_z] = 5,5,2. Time taken = 0.000 sec
countpairs_rp_pi_double> gridlink seems inefficient. nmesh = (5, 5, 2); avg_np = 0.22. Boosting bin refine factor - should lead to better performance
xmin = 0.662654 xmax=118.426432 rpmax = 45.000000
In gridlink_double> Running with [nmesh_x, nmesh_y, nmesh_z] = 8,8,2. Time taken = 0.000 sec
In gridlink_double> Running with [nmesh_x, nmesh_y, nmesh_z] = 8,8,2. Time taken = 0.000 sec
Using fallback kernel
0%.........10%.........20%.........30%.........40%.........50%.........60%.........70%.........80%.........90%.........100% done. Time taken = 0.001 secs
testing particles one by one...
def missed part: 41.087356477185 48.94941175529149 48 0.0
dd corrfunc: [362. 93.]
dd brute : [387. 96.]
Traceback (most recent call last):
File "corrfunc_test.py", line 99, in <module>
assert np.array_equal(dd_cf_sum, dd)
AssertionError