How to use the treecorr.KGCorrelation function in TreeCorr

To help you get started, we’ve selected a few TreeCorr examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github rmjarvis / TreeCorr / tests / test_patch.py View on Github external
print('NK2: ',s1, t1-t0, s1-s0)
    gk_cat2.unload()
    np.testing.assert_allclose(nk1.xi, nk2.xi)
    np.testing.assert_allclose(nk1.weight, nk2.weight)

    # KG without r_col, and not from a file.
    gk_cat1 = treecorr.Catalog(ra=ra[:ngal//100], dec=dec[:ngal//100],
                               g1=g1[:ngal//100], g2=g2[:ngal//100], k=k[:ngal//100],
                               ra_units='rad', dec_units='rad',
                               patch_centers=patch_centers)
    gk_cat2 = treecorr.Catalog(ra=ra[:ngal//100], dec=dec[:ngal//100],
                               g1=g1[:ngal//100], g2=g2[:ngal//100], k=k[:ngal//100],
                               ra_units='rad', dec_units='rad',
                               patch_centers=patch_centers, save_patch_dir='output')

    kg1 = treecorr.KGCorrelation(bin_size=0.5, min_sep=1., max_sep=30., sep_units='arcmin')
    t0 = time.time()
    s0 = hp.heap().size if hp else 0
    kg1.process(gk_cat1, gk_cat1)
    t1 = time.time()
    s1 = hp.heap().size if hp else 0
    print('KG1: ',s1, t1-t0, s1-s0)
    gk_cat1.unload()
    kg2 = treecorr.KGCorrelation(bin_size=0.5, min_sep=1., max_sep=30., sep_units='arcmin')
    t0 = time.time()
    s0 = hp.heap().size if hp else 0
    kg2.process(gk_cat2, gk_cat2, low_mem=True)
    t1 = time.time()
    s1 = hp.heap().size if hp else 0
    print('KG2: ',s1, t1-t0, s1-s0)
    gk_cat2.unload()
    np.testing.assert_allclose(kg1.xi, kg2.xi)
github rmjarvis / TreeCorr / tests / test_kg.py View on Github external
z2 = rng.normal(0,s, (ngal,) )
    w2 = rng.random_sample(ngal)
    g12 = rng.normal(0,0.2, (ngal,) )
    g22 = rng.normal(0,0.2, (ngal,) )

    ra1, dec1 = coord.CelestialCoord.xyz_to_radec(x1,y1,z1)
    ra2, dec2 = coord.CelestialCoord.xyz_to_radec(x2,y2,z2)

    cat1 = treecorr.Catalog(ra=ra1, dec=dec1, ra_units='rad', dec_units='rad', w=w1, k=k1)
    cat2 = treecorr.Catalog(ra=ra2, dec=dec2, ra_units='rad', dec_units='rad', w=w2, g1=g12, g2=g22)

    min_sep = 1.
    max_sep = 10.
    nbins = 50
    bin_size = np.log(max_sep/min_sep) / nbins
    kg = treecorr.KGCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                sep_units='deg', brute=True)
    kg.process(cat1, cat2)

    r1 = np.sqrt(x1**2 + y1**2 + z1**2)
    r2 = np.sqrt(x2**2 + y2**2 + z2**2)
    x1 /= r1;  y1 /= r1;  z1 /= r1
    x2 /= r2;  y2 /= r2;  z2 /= r2

    north_pole = coord.CelestialCoord(0*coord.radians, 90*coord.degrees)

    true_npairs = np.zeros(nbins, dtype=int)
    true_weight = np.zeros(nbins, dtype=float)
    true_xi = np.zeros(nbins, dtype=complex)

    c1 = [coord.CelestialCoord(r*coord.radians, d*coord.radians) for (r,d) in zip(ra1, dec1)]
    c2 = [coord.CelestialCoord(r*coord.radians, d*coord.radians) for (r,d) in zip(ra2, dec2)]
github rmjarvis / TreeCorr / tests / test_kg.py View on Github external
treecorr.corr2(config)
    data = fitsio.read(config['kg_file_name'])
    np.testing.assert_allclose(data['r_nom'], kg.rnom)
    np.testing.assert_allclose(data['npairs'], kg.npairs)
    np.testing.assert_allclose(data['weight'], kg.weight)
    np.testing.assert_allclose(data['kgamT'], kg.xi, rtol=1.e-3)
    np.testing.assert_allclose(data['kgamX'], kg.xi_im, rtol=1.e-3)

    # Invalid with only one file_name
    del config['file_name2']
    with assert_raises(TypeError):
        treecorr.corr2(config)

    # Repeat with binslop = 0, since code is different for bin_slop=0 and brute=True.
    # And don't do any top-level recursion so we actually test not going to the leaves.
    kg = treecorr.KGCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins, bin_slop=0,
                                max_top=0)
    kg.process(cat1, cat2)
    np.testing.assert_array_equal(kg.npairs, true_npairs)
    np.testing.assert_allclose(kg.weight, true_weight, rtol=1.e-5, atol=1.e-8)
    np.testing.assert_allclose(kg.xi, true_xi.real, rtol=1.e-3, atol=1.e-3)
    np.testing.assert_allclose(kg.xi_im, true_xi.imag, rtol=1.e-3, atol=1.e-3)

    # Check a few basic operations with a KGCorrelation object.
    do_pickle(kg)

    kg2 = kg.copy()
    kg2 += kg
    np.testing.assert_allclose(kg2.npairs, 2*kg.npairs)
    np.testing.assert_allclose(kg2.weight, 2*kg.weight)
    np.testing.assert_allclose(kg2.meanr, 2*kg.meanr)
    np.testing.assert_allclose(kg2.meanlogr, 2*kg.meanlogr)
github rmjarvis / TreeCorr / tests / test_patch.py View on Github external
ra_units='rad', dec_units='rad',
                               patch_centers=patch_centers)
    gk_cat2 = treecorr.Catalog(ra=ra[:ngal//100], dec=dec[:ngal//100],
                               g1=g1[:ngal//100], g2=g2[:ngal//100], k=k[:ngal//100],
                               ra_units='rad', dec_units='rad',
                               patch_centers=patch_centers, save_patch_dir='output')

    kg1 = treecorr.KGCorrelation(bin_size=0.5, min_sep=1., max_sep=30., sep_units='arcmin')
    t0 = time.time()
    s0 = hp.heap().size if hp else 0
    kg1.process(gk_cat1, gk_cat1)
    t1 = time.time()
    s1 = hp.heap().size if hp else 0
    print('KG1: ',s1, t1-t0, s1-s0)
    gk_cat1.unload()
    kg2 = treecorr.KGCorrelation(bin_size=0.5, min_sep=1., max_sep=30., sep_units='arcmin')
    t0 = time.time()
    s0 = hp.heap().size if hp else 0
    kg2.process(gk_cat2, gk_cat2, low_mem=True)
    t1 = time.time()
    s1 = hp.heap().size if hp else 0
    print('KG2: ',s1, t1-t0, s1-s0)
    gk_cat2.unload()
    np.testing.assert_allclose(kg1.xi, kg2.xi)
    np.testing.assert_allclose(kg1.weight, kg2.weight)
github rmjarvis / TreeCorr / tests / test_kg.py View on Github external
# In addition to the shape noise below, there is shot noise from the random x,y positions.
        x = (rng.random_sample(ngal)-0.5) * L
        y = (rng.random_sample(ngal)-0.5) * L
        # Varied weights are hard, but at least check that non-unit weights work correctly.
        w = np.ones_like(x) * 5
        r2 = (x**2 + y**2)/r0**2
        g1 = -gamma0 * np.exp(-r2/2.) * (x**2-y**2)/r0**2
        g2 = -gamma0 * np.exp(-r2/2.) * (2.*x*y)/r0**2
        k = kappa0 * np.exp(-r2/2.)
        # This time, add some shape noise (different each run).
        g1 += rng.normal(0, 0.3, size=ngal)
        g2 += rng.normal(0, 0.3, size=ngal)
        k += rng.normal(0, 0.1, size=ngal)

        cat = treecorr.Catalog(x=x, y=y, w=w, g1=g1, g2=g2, k=k)
        kg = treecorr.KGCorrelation(bin_size=0.1, min_sep=10., max_sep=100.)
        kg.process(cat, cat)
        all_kgs.append(kg)

    mean_xi = np.mean([kg.xi for kg in all_kgs], axis=0)
    var_xi = np.var([kg.xi for kg in all_kgs], axis=0)
    mean_varxi = np.mean([kg.varxi for kg in all_kgs], axis=0)

    print('mean_xi = ',mean_xi)
    print('mean_varxi = ',mean_varxi)
    print('var_xi = ',var_xi)
    print('ratio = ',var_xi / mean_varxi)
    print('max relerr for xi = ',np.max(np.abs((var_xi - mean_varxi)/var_xi)))
    np.testing.assert_allclose(mean_varxi, var_xi, rtol=0.02 * tol_factor)
github rmjarvis / TreeCorr / tests / test_twod.py View on Github external
print('max abs diff ng.xi = ',np.max(np.abs(ng1.xi - ng0.xi)))
    np.testing.assert_array_equal(ng1.npairs, ng0.npairs)
    np.testing.assert_allclose(ng1.xi, ng0.xi, atol=2.e-4)

    ng2 = treecorr.NGCorrelation(max_sep=max_sep, nbins=nbins, bin_type='TwoD', bin_slop=0.1)
    t0 = time.time()
    ng2.process(cat, cat)
    t1 = time.time()
    print('t for bs=0.1 = ',t1-t0)
    print('max abs diff npairs = ',np.max(np.abs(ng2.npairs - ng0.npairs)))
    print('max rel diff npairs = ',np.max(np.abs(ng2.npairs - ng0.npairs)/np.abs(ng0.npairs)))
    print('max abs diff ng.xi = ',np.max(np.abs(ng2.xi - ng0.xi)))
    np.testing.assert_allclose(ng2.npairs, ng0.npairs, rtol=1.e-2)
    np.testing.assert_allclose(ng2.xi, ng0.xi, atol=5.e-4)

    kg1 = treecorr.KGCorrelation(max_sep=max_sep, nbins=nbins, bin_type='TwoD', bin_slop=0)
    t0 = time.time()
    kg1.process(cat, cat)
    t1 = time.time()
    print('t for bs=1.e-10 = ',t1-t0)
    print('ng0.xi = ',ng0.xi)
    print('ng1.xi = ',ng1.xi)
    print('kg1.xi = ',kg1.xi)
    print('max abs diff kg.xi = ',np.max(np.abs(kg1.xi - 10.*ng0.xi)))
    np.testing.assert_array_equal(kg1.npairs, ng0.npairs)
    np.testing.assert_allclose(kg1.xi, 10.*ng0.xi, atol=2.e-3)
    np.testing.assert_array_equal(kg1.npairs, ng1.npairs)
    np.testing.assert_allclose(kg1.xi, 10.*ng1.xi, atol=1.e-10)

    kg2 = treecorr.KGCorrelation(max_sep=max_sep, nbins=nbins, bin_type='TwoD', bin_slop=0.1)
    t0 = time.time()
    kg2.process(cat, cat)
github rmjarvis / TreeCorr / tests / test_kg.py View on Github external
np.testing.assert_allclose(kg2.meanlogr, 2*kg.meanlogr)
    np.testing.assert_allclose(kg2.xi, 2*kg.xi)
    np.testing.assert_allclose(kg2.xi_im, 2*kg.xi_im)

    kg2.clear()
    kg2 += kg
    np.testing.assert_allclose(kg2.npairs, kg.npairs)
    np.testing.assert_allclose(kg2.weight, kg.weight)
    np.testing.assert_allclose(kg2.meanr, kg.meanr)
    np.testing.assert_allclose(kg2.meanlogr, kg.meanlogr)
    np.testing.assert_allclose(kg2.xi, kg.xi)
    np.testing.assert_allclose(kg2.xi_im, kg.xi_im)

    ascii_name = 'output/kg_ascii.txt'
    kg.write(ascii_name, precision=16)
    kg3 = treecorr.KGCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins)
    kg3.read(ascii_name)
    np.testing.assert_allclose(kg3.npairs, kg.npairs)
    np.testing.assert_allclose(kg3.weight, kg.weight)
    np.testing.assert_allclose(kg3.meanr, kg.meanr)
    np.testing.assert_allclose(kg3.meanlogr, kg.meanlogr)
    np.testing.assert_allclose(kg3.xi, kg.xi)
    np.testing.assert_allclose(kg3.xi_im, kg.xi_im)

    fits_name = 'output/kg_fits.fits'
    kg.write(fits_name)
    kg4 = treecorr.KGCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins)
    kg4.read(fits_name)
    np.testing.assert_allclose(kg4.npairs, kg.npairs)
    np.testing.assert_allclose(kg4.weight, kg.weight)
    np.testing.assert_allclose(kg4.meanr, kg.meanr)
    np.testing.assert_allclose(kg4.meanlogr, kg.meanlogr)
github rmjarvis / TreeCorr / treecorr / corr2.py View on Github external
if config.get('nk_statistic',None) == 'compensated':
                raise TypeError("rand_files is required for nk_statistic = compensated")
        elif config.get('nk_statistic','compensated') == 'compensated':
            rk = treecorr.NKCorrelation(config,logger)
            rk.process(rand1,cat2)
            logger.info("Done RK calculation.")

        nk.write(config['nk_file_name'], rk)
        logger.warning("Wrote NK correlation to %s",config['nk_file_name'])

    # Do KG correlation function if necessary
    if 'kg_file_name' in config:
        if cat2 is None:
            raise TypeError("file_name2 is required for kg correlation")
        logger.warning("Performing KG calculations...")
        kg = treecorr.KGCorrelation(config,logger)
        kg.process(cat1,cat2)
        logger.info("Done KG calculation.")
        kg.write(config['kg_file_name'])
        logger.warning("Wrote KG correlation to %s",config['kg_file_name'])