diff --git a/tests/test_covariance_fourier_cNG.py b/tests/test_covariance_fourier_cNG.py index 9b6b3d89..60de56dd 100644 --- a/tests/test_covariance_fourier_cNG.py +++ b/tests/test_covariance_fourier_cNG.py @@ -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) @@ -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]] @@ -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() @@ -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, @@ -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 diff --git a/tjpcov/covariance_fourier_cNG.py b/tjpcov/covariance_fourier_cNG.py index 4b77f380..267ae035 100644 --- a/tjpcov/covariance_fourier_cNG.py +++ b/tjpcov/covariance_fourier_cNG.py @@ -3,7 +3,6 @@ # import healpy as hp import numpy as np import pyccl as ccl -import warnings from .covariance_builder import CovarianceFourier @@ -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( @@ -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) @@ -179,10 +178,10 @@ 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( @@ -190,33 +189,28 @@ def get_covariance_block( ) 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, @@ -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: