Skip to content

Commit

Permalink
Tests created and pass locally
Browse files Browse the repository at this point in the history
  • Loading branch information
@paulrogozenski
paulrogozenski committed Dec 13, 2024
1 parent 7d4619a commit ecd7f0b
Show file tree
Hide file tree
Showing 2 changed files with 119 additions and 34 deletions.
108 changes: 100 additions & 8 deletions tests/test_covariance_fourier_cNG.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,37 @@ def get_config():
return config


def get_hod_model():
obj = FouriercNGHaloModel(INPUT_YML_cNG)
mass_def = ccl.halos.MassDef200m
cM = ccl.halos.ConcentrationDuffy08(mass_def=mass_def)
hod = ccl.halos.HaloProfileHOD(
mass_def=mass_def,
concentration=cM,
log10Mmin_0=obj.HOD_dict["log10Mmin_0"],
log10Mmin_p=obj.HOD_dict["log10Mmin_p"],
siglnM_0=obj.HOD_dict["siglnM_0"],
siglnM_p=obj.HOD_dict["siglnM_p"],
log10M0_0=obj.HOD_dict["log10M0_0"],
log10M0_p=obj.HOD_dict["log10M0_p"],
log10M1_0=obj.HOD_dict["log10M1_0"],
log10M1_p=obj.HOD_dict["log10M1_p"],
alpha_0=obj.HOD_dict["alpha_0"],
alpha_p=obj.HOD_dict["alpha_p"],
fc_0=obj.HOD_dict["fc_0"],
fc_p=obj.HOD_dict["fc_p"],
bg_0=obj.HOD_dict["bg_0"],
bg_p=obj.HOD_dict["bg_p"],
bmax_0=obj.HOD_dict["bmax_0"],
bmax_p=obj.HOD_dict["bmax_p"],
a_pivot=obj.HOD_dict["a_pivot"],
ns_independent=obj.HOD_dict["ns_independent"],
is_number_counts=obj.HOD_dict["is_number_counts"],
)

return hod


def get_halo_model(cosmo):
md = ccl.halos.MassDef200m
mf = ccl.halos.MassFuncTinker08(mass_def=md)
Expand Down Expand Up @@ -115,15 +146,29 @@ def test_get_covariance_block(cov_fcNG, tracer_comb1, tracer_comb2):
np.max(np.abs((covf["cov_nob"] + 1e-100) / (cov_cNG + 1e-100) - 1))
< 1e-10
)

# CCL covariance
na = ccl.ccllib.get_pk_spline_na(cosmo.cosmo)
a_arr, _ = ccl.ccllib.get_pk_spline_a(cosmo.cosmo, na, 0)

# TODO: Need to make 1h TK3D object with HOD
# & weight non-HOD TK3D object with gbias factors
# & combine together before call to
# angular_cl_cov_cNG for proper comparison
tr = {}
tr[1], tr[2] = tracer_comb1
tr[3], tr[4] = tracer_comb2
z_max = []
for i in range(4):
tr_sacc = s.tracers[tr[i + 1]]
z = tr_sacc.z
z_max.append(z.max())
# Divide by zero errors happen when default a_arr used for 1h term
z_max = np.min(z_max)

# Array of a.
# Use the a's in the pk spline
na = ccl.ccllib.get_pk_spline_na(cosmo.cosmo)
a_arr, _ = ccl.ccllib.get_pk_spline_a(cosmo.cosmo, na, 0)
# Cut the array for efficiency
sel = 1 / a_arr < z_max + 1
# Include the next node so that z_max is in the range
sel[np.sum(~sel) - 1] = True
a_arr = a_arr[sel]
bias1 = bias2 = bias3 = bias4 = 1
if "gc" in tracer_comb1[0]:
bias1 = cov_fcNG.bias_lens[tracer_comb1[0]]
Expand All @@ -136,13 +181,34 @@ def test_get_covariance_block(cov_fcNG, tracer_comb1, tracer_comb2):

if "gc" in tracer_comb2[0]:
bias4 = cov_fcNG.bias_lens[tracer_comb2[1]]


biases = bias1 * bias2 * bias3 * bias4

hmc = get_halo_model(cosmo)
nfw_profile = get_NFW_profile()
hod = get_hod_model()
prof_2pt = ccl.halos.profiles_2pt.Profile2ptHOD()

tkk_cNG = ccl.halos.halomod_Tk3D_cNG(
cosmo,
hmc,
prof=nfw_profile,
separable_growth=True,
a_arr=a_arr,
)
tkk_1h_nfw = ccl.halos.halomod_Tk3D_1h(
cosmo,
hmc,
prof=nfw_profile,
a_arr=a_arr,
)
tkk_1h_hod = ccl.halos.halomod_Tk3D_1h(
cosmo,
hmc,
prof=hod,
prof12_2pt=prof_2pt,
prof34_2pt=prof_2pt,
a_arr=a_arr,
)

ccl_tracers, _ = cov_fcNG.get_tracer_info()
Expand All @@ -151,7 +217,7 @@ def test_get_covariance_block(cov_fcNG, tracer_comb1, tracer_comb2):
tr3 = ccl_tracers[tracer_comb2[0]]
tr4 = ccl_tracers[tracer_comb2[1]]

fsky = get_fsky(tr1, tr2, tr3, tr4)
fsky = get_fsky(*tracer_comb1, *tracer_comb2)

cov_ccl = ccl.covariances.angular_cl_cov_cNG(
cosmo,
Expand All @@ -164,6 +230,32 @@ def test_get_covariance_block(cov_fcNG, tracer_comb1, tracer_comb2):
fsky=fsky,
)

cov_ccl_1h_nfw = ccl.covariances.angular_cl_cov_cNG(
cosmo,
tracer1=tr1,
tracer2=tr2,
tracer3=tr3,
tracer4=tr4,
ell=ell,
t_of_kk_a=tkk_1h_nfw,
fsky=fsky,
)

cov_ccl_1h_hod = ccl.covariances.angular_cl_cov_cNG(
cosmo,
tracer1=tr1,
tracer2=tr2,
tracer3=tr3,
tracer4=tr4,
ell=ell,
t_of_kk_a=tkk_1h_hod,
fsky=fsky,
)
# An unfortunately messy way to to calculate the 234h terms
# with an NFW Profile and only the 1h term with an HOD
# using current CCL infrastructure.
cov_ccl = biases * (cov_ccl - cov_ccl_1h_nfw) + cov_ccl_1h_hod

assert np.max(np.fabs(np.diag(cov_cNG / cov_ccl - 1))) < 1e-5
assert np.max(np.fabs(cov_cNG / cov_ccl - 1)) < 1e-3

Expand Down
45 changes: 19 additions & 26 deletions tjpcov/covariance_fourier_cNG.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
# import healpy as hp
import numpy as np
import pyccl as ccl
import warnings

from .covariance_builder import CovarianceFourier

Expand Down Expand Up @@ -104,17 +103,19 @@ def get_covariance_block(
tr = {}
tr[1], tr[2] = tracer_comb1
tr[3], tr[4] = tracer_comb2

cosmo = self.get_cosmology()
mass_def = ccl.halos.MassDef200m
hmf = ccl.halos.MassFuncTinker08(mass_def=mass_def)
hbf = ccl.halos.HaloBiasTinker10(mass_def=mass_def)
cM = ccl.halos.ConcentrationDuffy08(mass_def=mass_def)
prof_2pt = ccl.halos.profiles_2pt.Profile2ptHOD()
nfw = ccl.halos.HaloProfileNFW(
mass_def=mass_def, concentration=cM, fourier_analytic=True
)
hmc = ccl.halos.HMCalculator(
mass_function=hmf, halo_bias=hbf, mass_def=mass_def
mass_function=hmf,
halo_bias=hbf,
mass_def=mass_def,
)

hod = ccl.halos.HaloProfileHOD(
Expand Down Expand Up @@ -148,10 +149,8 @@ def get_covariance_block(
for i in range(4):
tr_sacc = sacc_file.tracers[tr[i + 1]]
z = tr_sacc.z
# z, nz = tr_sacc.z, tr_sacc.nz
# z_min.append(z[np.where(nz > 0)[0][0]])
# z_max.append(z[np.where(np.cumsum(nz)/np.sum(nz) > 0.999)[0][0]])
z_max.append(z.max())
# Divide by zero errors happen when default a_arr used for 1h term

z_max = np.min(z_max)

Expand Down Expand Up @@ -179,44 +178,39 @@ def get_covariance_block(
# TODO: This should be unified with the other classes in
# CovarianceBuilder.
fsky = np.mean(masks[1] * masks[2] * masks[3] * masks[4])

# Tk3D = b1*b2*b3*b4 * T_234h (NFW) + T_1h (HOD)

tkk = ccl.halos.pk_4pt.halomod_trispectrum_2h_22(
cosmo, hmc, np.exp(lk_arr), a_arr, prof=nfw
cosmo, hmc, np.exp(lk_arr), a_arr, prof=nfw, separable_growth=True
)

tkk += ccl.halos.halomod_trispectrum_2h_13(
cosmo, hmc, np.exp(lk_arr), a_arr, prof=nfw
)

tkk += ccl.halos.halomod_trispectrum_3h(
cosmo, hmc, np.exp(lk_arr), a_arr, prof=nfw
cosmo, hmc, np.exp(lk_arr), a_arr, prof=nfw, separable_growth=True
)

tkk += ccl.halos.halomod_trispectrum_4h(
cosmo, hmc, np.exp(lk_arr), a_arr, prof=nfw
cosmo, hmc, np.exp(lk_arr), a_arr, prof=nfw, separable_growth=True
)

tkk *= bias1 * bias2 * bias3 * bias4

tkk += ccl.halos.halomod_trispectrum_1h(
cosmo, hmc, np.exp(lk_arr), a_arr, prof=hod
cosmo,
hmc,
np.exp(lk_arr),
a_arr,
prof=hod,
prof12_2pt=prof_2pt,
prof34_2pt=prof_2pt,
)

s = self.io.get_sacc_file()
isnc = []
for i in range(1, 5):
isnc[i] = (s.tracers[tr[i]].quantity == "galaxy_density") or (
"lens" in tr[i]
)
if any(isnc):
warnings.warn(
"Using linear galaxy bias with 1h term. This should "
"be checked. HOD version need implementation."
)

tk3D = ccl.tk3d.Tk3D(a_arr=a_arr, lk_arr=lk_arr, tkk_arr=tkk)

tk3D = ccl.tk3d.Tk3D(
a_arr=a_arr, lk_arr=lk_arr, tkk_arr=tkk, is_logt=False
)
ell = self.get_ell_eff()
cov_cng = ccl.covariances.angular_cl_cov_cNG(
cosmo,
Expand All @@ -236,7 +230,6 @@ def get_covariance_block(
cov_full = np.zeros((nbpw, ncell1, nbpw, ncell2))
cov_full[:, 0, :, 0] = cov_cng
cov_full = cov_full.reshape((nbpw * ncell1, nbpw * ncell2))

np.savez_compressed(fname, cov=cov_full, cov_nob=cov_cng)

if not include_b_modes:
Expand Down

0 comments on commit ecd7f0b

Please sign in to comment.