DDrppi function seems not to return expected pair counts

[flepri95]

General information

• Corrfunc version: 2.3.4
• platform: MacOS
• installation method (pip/source/other?): pip

Issue description

DDrppi function doesn't return expected pair counts, calculated through two-loop pairs counting. I can't figure out if there is a mistake in my thought process (and code) or something more in Corrfunc. If there is a mistake on my part, I apologize for wasting your time.

Minimal failing example

import numpy as np
import math as mt
import Corrfunc
from Corrfunc.theory.DDrppi import DDrppi

def pi(x1, y1, z1, x2, y2, z2):
        return (mt.fabs((x1*x1-x2*x2)+(y1*y1-y2*y2)+(z1*z1-z2*z2))/mt.sqrt((x1+x2)*(x1+x2)+(y1+y2)*(y1+y2)+(z1+z2)*(z1+z2)))

def rp (x1, y1, z1, x2, y2, z2, pi):
        return(mt.sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)-pi*pi))

rpmin=0.5
rpmax=40.0
s_rp=0.5
nbin_rp=int(np.floor((rpmax-rpmin)/s_rp))

pimin=0.0
pimax=40.0
s_pi=1.0
nbin_pi=int(np.floor((pimax-pimin)/s_pi))

counts_rppi=np.zeros(shape=(nbin_rp,nbin_pi), dtype='float')

rpbins=np.linspace(rpmin, rpmax, nbin_rp+1)

N = 10000

boxsize = 420.0
nthreads = 4
autocorr = 1
seed = 42
np.random.seed(seed)
X = np.random.uniform(0, boxsize, N)
Y = np.random.uniform(0, boxsize, N)
Z = np.random.uniform(0, boxsize, N)
weights = np.ones_like(X)

results = DDrppi(autocorr, nthreads, pimax, rpbins, X, Y, Z)

for i in range(0, N):
    for j in range(0, i):
        r_pi=pi(X[i], Y[i], Z[i], X[j], Y[j], Z[j])
        r_rp=rp(X[i], Y[i], Z[i], X[j], Y[j], Z[j], r_pi)
        if(r_rp>rpmin and r_rp<rpmax and r_pi>pimin and r_pi<pimax):
            bin_rp=int(np.floor((r_rp-rpmin)/s_rp))
            bin_pi=int(np.floor((r_pi-pimin)/s_pi))
            counts_rppi[bin_rp][bin_pi]+=2.

Thanks!

[lgarrison]

Hi, at a glance, the pi function looks suspicious, as the line-of-sight distance (in the theory module) is just abs(z1-z2). Could that be the cause of the difference?

[flepri95]

It was indeed that. I'm so sorry for wasting your time! Fixing rp and pi definition, the pair counts are the same.

[lgarrison]

No worries, glad it's fixed!