Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
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)
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)]
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)
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)
# 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)
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)
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)
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'])