Corrfunc npibins bug (issue #96)

test_npibins_bug.py

from __future__ import print_function

import numpy as np
from Corrfunc.mocks import DDrppi_mocks
from Corrfunc.io import read_catalog
from Corrfunc.utils import convert_rp_pi_counts_to_wp
from Corrfunc.utils import convert_3d_counts_to_cf
import Corrfunc
from os.path import dirname, abspath, join as pjoin
from IPython.core.debugger import Tracer


def main():
    randoms_file = pjoin(dirname(abspath(Corrfunc.__file__)),
                         '../mocks/tests/data/',
                         'Mr19_randoms_northonly.rdcz.ff')
    galaxy_file = pjoin(dirname(abspath(Corrfunc.__file__)),
                        '../mocks/tests/data/',
                        'Mr19_mock_northonly.rdcz.ff')

    ra, dec, cz = read_catalog(galaxy_file)
    rand_ra, rand_dec, rand_cz = read_catalog(randoms_file)

    npts = len(ra)
    npts_randoms = 2 * npts
    seed = 42
    np.random.seed(seed)
    inds = np.random.randint(0, len(rand_ra), size=npts_randoms)
    rand_ra = rand_ra[inds]
    rand_dec = rand_dec[inds]
    rand_cz = rand_cz[inds]

    verbose = False
    nbins = 10
    
    # BUG
    bins = np.logspace(-1, 1.35, nbins)

    # CORRECT
    bins = np.logspace(-1, 1.35, nbins + 1)
    
    pimax = 40.0
    cosmology = 1
    nthreads = 4
    
    autocorr = 1
    DD_counts = DDrppi_mocks(autocorr, cosmology, nthreads, pimax, bins,
                             ra, dec, cz, verbose=verbose)
    
    autocorr = 0
    DR_counts = DDrppi_mocks(autocorr, cosmology, nthreads, pimax, bins,
                             ra, dec, cz,
                             RA2=rand_ra,
                             DEC2=rand_dec,
                             CZ2=rand_cz, verbose=verbose)

    autocorr = 1
    RR_counts = DDrppi_mocks(autocorr, cosmology, nthreads, pimax, bins,
                             rand_ra,
                             rand_dec,
                             rand_cz, verbose=verbose)

    npibins = len(DD_counts) // nbins
    dpi = pimax/npibins
    
    try:
        if dpi <= 0.0:
            msg = 'Binsize along the line of sight (dpi) = {0}'\
                   'must be positive'.format(dpi)
            raise ValueError(msg)
                        
        xirppi = convert_3d_counts_to_cf(npts, npts,
                                         npts_randoms, npts_randoms,
                                         DD_counts,
                                         DR_counts,
                                         DR_counts,
                                         RR_counts)
        wp = np.empty(nbins)
        npibins = len(xirppi) // nbins
        print("len(xirppi) = {0}".format(len(xirppi)))
        if ((npibins * nbins) != len(xirppi)):
            msg = 'Number of pi bins could not be calculated correctly.'\
                  'Expected to find that the total number of bins = {0} '\
                  'would be the product of the number of pi bins = {1} '\
                  'and the number of rp bins = {2}'.format(len(xirppi),
                                                           npibins,
                                                           nbins)
            raise ValueError(msg)
        
        for i in range(nbins):
            wp[i] = 2.0 * dpi * np.sum(xirppi[i * npibins:(i + 1) * npibins])

        print("wp size = {0}\nwp = {1}".format(wp.size, wp))
        
        wp = convert_rp_pi_counts_to_wp(npts, npts, npts_randoms, npts_randoms,
                                        DD_counts,
                                        DR_counts,
                                        DR_counts,
                                        RR_counts,
                                        nbins)
        print("wp size = {0}\nwp = {1}".format(wp.size, wp))
    except:
        Tracer()()
        

if __name__ == '__main__':
    main()

this gist will not run after the latest commit (4546ca6) that fixed issue 96. The latest commit requires a pimax argument to convert_rp_pi_counts_to_wp as a safety check.