Skip to content

Correctly convert counts to cf when doing an autocorr #201

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 26 additions & 9 deletions Corrfunc/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@

def convert_3d_counts_to_cf(ND1, ND2, NR1, NR2,
D1D2, D1R2, D2R1, R1R2,
estimator='LS'):
estimator='LS', autocorr=False):
"""
Converts raw pair counts to a correlation function.

Expand Down Expand Up @@ -61,6 +61,9 @@ def convert_3d_counts_to_cf(ND1, ND2, NR1, NR2,
The kind of estimator to use for computing the correlation
function. Currently, only supports Landy-Szalay

autocorr: bool, default=False
Whether the counts are from an autocorrelation

Returns
---------

Expand Down Expand Up @@ -138,14 +141,24 @@ def convert_3d_counts_to_cf(ND1, ND2, NR1, NR2,

nonzero = pair_counts['R1R2'] > 0
if 'LS' in estimator or 'Landy' in estimator:
fN1 = np.float(NR1) / np.float(ND1)
fN2 = np.float(NR2) / np.float(ND2)
cf = np.zeros(nbins)
cf[:] = np.nan
cf[nonzero] = (fN1 * fN2 * pair_counts['D1D2'][nonzero] -
fN1 * pair_counts['D1R2'][nonzero] -
fN2 * pair_counts['D2R1'][nonzero] +
pair_counts['R1R2'][nonzero]) / pair_counts['R1R2'][nonzero]
if autocorr:
fN1 = np.float(NR1) / np.float(ND1)
fN2 = (np.float(NR1) - 1) / (np.float(ND1) - 1)
fN3 = (np.float(NR1) - 1) / np.float(ND1)
cf[nonzero] = (fN1 * fN2 * pair_counts['D1D2'][nonzero] -
2 * fN3 * pair_counts['D1R2'][nonzero] +
pair_counts['R1R2'][nonzero]
) / pair_counts['R1R2'][nonzero]
else:
fN1 = np.float(NR1) / np.float(ND1)
fN2 = np.float(NR2) / np.float(ND2)
cf[nonzero] = (fN1 * fN2 * pair_counts['D1D2'][nonzero] -
fN1 * pair_counts['D1R2'][nonzero] -
fN2 * pair_counts['D2R1'][nonzero] +
pair_counts['R1R2'][nonzero]
) / pair_counts['R1R2'][nonzero]
if len(cf) != nbins:
msg = 'Bug in code. Calculated correlation function does not '\
'have the same number of bins as input arrays. Input bins '\
Expand All @@ -163,7 +176,7 @@ def convert_3d_counts_to_cf(ND1, ND2, NR1, NR2,
def convert_rp_pi_counts_to_wp(ND1, ND2, NR1, NR2,
D1D2, D1R2, D2R1, R1R2,
nrpbins, pimax, dpi=1.0,
estimator='LS'):
estimator='LS', autocorr=False):
"""
Converts raw pair counts to a correlation function.

Expand Down Expand Up @@ -210,6 +223,10 @@ def convert_rp_pi_counts_to_wp(ND1, ND2, NR1, NR2,
The kind of estimator to use for computing the correlation
function. Currently, only supports Landy-Szalay

autocorr: bool, default=False
Whether the counts are from an autocorrelation


Returns
---------

Expand Down Expand Up @@ -289,7 +306,7 @@ def convert_rp_pi_counts_to_wp(ND1, ND2, NR1, NR2,

xirppi = convert_3d_counts_to_cf(ND1, ND2, NR1, NR2,
D1D2, D1R2, D2R1, R1R2,
estimator=estimator)
estimator=estimator, autocorr=autocorr)
wp = np.empty(nrpbins)
npibins = len(xirppi) // nrpbins
if ((npibins * nrpbins) != len(xirppi)):
Expand Down