diff --git a/.github/environment.yml b/.github/environment.yml index 3a60efe78..e79ec9d4f 100644 --- a/.github/environment.yml +++ b/.github/environment.yml @@ -30,4 +30,4 @@ dependencies: - MiraTitanHMFemulator - dark_emulator==1.1.2 - colossus - + \ No newline at end of file diff --git a/benchmarks/data/pkresponse.npz b/benchmarks/data/pkresponse.npz new file mode 100644 index 000000000..16028f0db Binary files /dev/null and b/benchmarks/data/pkresponse.npz differ diff --git a/benchmarks/test_pk_response.py b/benchmarks/test_pk_response.py new file mode 100644 index 000000000..72e394119 --- /dev/null +++ b/benchmarks/test_pk_response.py @@ -0,0 +1,181 @@ +import numpy as np +from pyccl.pk_response import ( + resp_Pmm_hresponse, + resp_Pgm_darkemu, + resp_Pgg_darkemu, +) +import pyccl as ccl + +# Set cosmology +Om = 0.3156 +ob = 0.02225 +oc = 0.1198 +OL = 0.6844 +As = 2.2065e-9 +ns = 0.9645 +h = 0.6727 + +cosmo = ccl.Cosmology( + Omega_c=oc / (h**2), + Omega_b=ob / (h**2), + h=h, + n_s=ns, + A_s=As, + m_nu=0.06, + transfer_function="boltzmann_camb", + extra_parameters={"camb": {"halofit_version": "takahashi"}}, +) + + +# Load data +pkresponse_data = np.load("benchmarks/data/pkresponse.npz") +k_data = pkresponse_data["k_h"] +k_data_mm = pkresponse_data["k_h_mm"] +Pmm_resp_data = pkresponse_data["Pmm_resp"] / h**3 +Pgm_resp_data = pkresponse_data["Pgm_resp"] / h**3 +Pgg_resp_data = pkresponse_data["Pgg_resp"] / h**3 + + +# HOD parameters +logMhmin = 13.94 +logMh1 = 14.46 +alpha = 1.192 +kappa = 0.60 +sigma_logM = 0.5 +sigma_lM = sigma_logM * np.log(10) +logMh0 = logMhmin + np.log10(kappa) + +logMmin = np.log10(10**logMhmin / h) +logM0 = np.log10(10**logMh0 / h) +logM1 = np.log10(10**logMh1 / h) + +# Construct HOD +mass_def = ccl.halos.MassDef(200, "matter") +cm = ccl.halos.ConcentrationDuffy08(mass_def=mass_def) +prof_hod = ccl.halos.HaloProfileHOD( + mass_def=mass_def, + concentration=cm, + log10Mmin_0=logMmin, + siglnM_0=sigma_lM, + log10M0_0=logM0, + log10M1_0=logM1, + alpha_0=alpha, +) + +# Define input parameters for pkresponse functions +log10Mh_min = np.log10(2.6e12) +log10Mh_max = 15.9 +a_arr = np.array([1.0, 0.64516129, 0.49382716, 0.40387722]) +indx = (k_data > 1e-2) & (k_data < 3) +lk_arr = np.log(k_data[indx] * h) # Using loaded k_data +indx_mm = (k_data_mm > 1e-2) & (k_data_mm < 3) +lk_arr_mm = np.log(k_data_mm[indx_mm] * h) # Using loaded k_data + +# Generate power spectrum responses using pkresponse.py functions +generated_Pmm_resp = resp_Pmm_hresponse( + cosmo, deltah=0.02, lk_arr=lk_arr_mm, a_arr=a_arr +) + +generated_Pgm_resp = resp_Pgm_darkemu( + cosmo, + prof_hod, + deltah=0.02, + log10Mh_min=log10Mh_min, + log10Mh_max=log10Mh_max, + lk_arr=lk_arr, + a_arr=a_arr, +) + +generated_Pgg_resp = resp_Pgg_darkemu( + cosmo, + prof_hod, + deltalnAs=0.03, + log10Mh_min=log10Mh_min, + log10Mh_max=log10Mh_max, + lk_arr=lk_arr, + a_arr=a_arr, +) + + +# Compare the generated responses with simulation data. +# The error bars from the simulations are too strict, +# probably because we use the pair and fixed initial density method. +# Hence instead of "atol" (absolute difference compared to the error bars), +# we use "rtol" (relative difference compared to the simulated response). +# The difference between the model and the simulations are +# typically few ten percents. +# Higher redshifts of Pgm, Pgg (especially Pgg) have worse accuracy due to the +# accuracy of darkemulator itself, +# but they are still a factor of difference at most. + + +def test_pmm_resp(): + assert np.allclose( + generated_Pmm_resp, + Pmm_resp_data[:, indx_mm], + rtol=0.15, + ) + + +def test_pgm_resp(): + assert np.allclose( + generated_Pgm_resp[0], + Pgm_resp_data[0, indx], + rtol=0.3, + ) + + +def test_pgm_resp1(): + assert np.allclose( + generated_Pgm_resp[1], + Pgm_resp_data[1, indx], + rtol=0.36, + ) + + +def test_pgm_resp2(): + assert np.allclose( + generated_Pgm_resp[2], + Pgm_resp_data[2, indx], + rtol=0.78, + ) + + +def test_pgm_resp3(): + assert np.allclose( + generated_Pgm_resp[3], + Pgm_resp_data[3, indx], + rtol=0.85, + ) + + +def test_pgg_resp(): + assert np.allclose( + generated_Pgg_resp[0], + Pgg_resp_data[0, indx], + rtol=0.65, + ) + + +def test_pgg_resp1(): + assert np.allclose( + generated_Pgg_resp[1], + Pgg_resp_data[1, indx], + rtol=1.7, + ) + + +def test_pgg_resp2(): + assert np.allclose( + generated_Pgg_resp[2], + Pgg_resp_data[2, indx], + rtol=3.1, + ) + + +def test_pgg_resp3(): + assert np.allclose( + generated_Pgg_resp[3], + Pgg_resp_data[3, indx], + rtol=2.6, + ) diff --git a/pyccl/__init__.py b/pyccl/__init__.py index b16d66e2f..d95a579a0 100644 --- a/pyccl/__init__.py +++ b/pyccl/__init__.py @@ -45,3 +45,5 @@ from . import nl_pt from .cosmology import * + +from .pk_response import * \ No newline at end of file diff --git a/pyccl/pk_response.py b/pyccl/pk_response.py new file mode 100755 index 000000000..7a03489e7 --- /dev/null +++ b/pyccl/pk_response.py @@ -0,0 +1,829 @@ +import numpy as np + +from . import cosmology + +from dark_emulator import darkemu +from scipy import integrate +from . import halos, lib, check + +# use the perturbation theory below khmin +khmin = 1e-2 # [h/Mpc] + + +def resp_Pmm_hresponse( + cosmo, + deltah=0.02, + extra_parameters={ + "camb": { + "halofit_version": "original", + } + }, + lk_arr=None, + a_arr=None, + khmin=khmin, +): + """Implements the response of matter power spectrum to the long wavelength + modes developed in Terasawa et al. 2023 (arXiv:2310.13330) as: + + .. math:: + \\frac{\\partial P_{mm}(k)}{\\partial\\delta_b} = + \\left(1 + \\frac{26}{21}T_{h}^{mm}(k) - \\frac{1}{3} + \\frac{d\\log P_{mm}(k)}{d\\log k}\\right)P_{mm}(k), + + where the :math:`T_{h}^{mm}(k)` is the normalized growth response to + the Hubble parameter defined as + :math:`T_{h}^{mm}(k) + = \\frac{d\\log P_{mm}(k)}{dh}/(2\\frac{d\\log D}{dh})`. + + Args: + cosmo (:class:`~pyccl.core.Cosmology`): a Cosmology object. + deltah (float): the variation of :math:`h` to compute :math:`T_{h}(k)` + by the two-sided numerical derivative method. + extra_parameters (:obj:`dict`): Dictionary holding extra + parameters. Currently supports extra parameters for CAMB. + lk_arr (array): an array holding values of the natural + logarithm of the wavenumber (in units of :math:`Mpc^{-1}`) at + which the response is calculated. + a_arr (array): an array holding values of the scale factor + at which the response is calculated. + khmin (float): the wavenumber (in units of :math:`h Mpc^{-1}`) below + which the response is calculated using the perturbation theory. + Returns: + Response of the matter power spectrum. + """ + + # Set k and a sampling from CCL parameters + if lk_arr is None: + status = 0 + lk_arr, status = lib.get_pk_spline_lk( + cosmo.cosmo, lib.get_pk_spline_nk(cosmo.cosmo), status + ) + check(status, cosmo=cosmo) + if a_arr is None: + status = 0 + a_arr, status = lib.get_pk_spline_a( + cosmo.cosmo, lib.get_pk_spline_na(cosmo.cosmo), status + ) + check(status, cosmo=cosmo) + + k_use = np.exp(lk_arr) + + # set h-modified cosmology to take finite differencing + cosmo_hp, cosmo_hm = _set_hmodified_cosmology( + cosmo, deltah, extra_parameters + ) + + # Growth factor + Dp = cosmo_hp.growth_factor_unnorm(a_arr) + Dm = cosmo_hm.growth_factor_unnorm(a_arr) + + # Linear power spectrum + cosmo.compute_linear_power() + cosmo_hp.compute_linear_power() + cosmo_hm.compute_linear_power() + + pk2dlin = cosmo.get_linear_power("delta_matter:delta_matter") + + pk2d = cosmo.get_nonlin_power("delta_matter:delta_matter") + pk2d_hp = cosmo_hp.get_nonlin_power("delta_matter:delta_matter") + pk2d_hm = cosmo_hm.get_nonlin_power("delta_matter:delta_matter") + + na = len(a_arr) + nk = len(k_use) + dpk12 = np.zeros([na, nk]) + dpk = np.zeros(nk) + T_h = np.zeros(nk) + + # use the perturbation theory below kmin + kmin = khmin * cosmo["h"] + T_h[k_use <= kmin] = 1 + + for ia, aa in enumerate(a_arr): + pk = pk2d(k_use, aa, cosmo) + pk_hp = pk2d_hp(k_use, aa, cosmo_hp) + pk_hm = pk2d_hm(k_use, aa, cosmo_hm) + + dpknl = pk2d(k_use, aa, cosmo, derivative=True) + dpklin = pk2dlin(k_use, aa, cosmo, derivative=True) + + # Eq. 11 ((hp-hm) term is cancelled out) + T_h[k_use > kmin] = ( + np.log(pk_hp[k_use > kmin]) - np.log(pk_hm[k_use > kmin]) + ) / (2 * (np.log(Dp[ia]) - np.log(Dm[ia]))) + + dpk[k_use <= kmin] = dpklin[k_use <= kmin] + dpk[k_use > kmin] = dpknl[k_use > kmin] + + # Eq. 23 + dpk12[ia, :] = pk * (1.0 + (26.0 / 21.0) * T_h - dpk / 3.0) + + return dpk12 + + +def resp_Pgm_darkemu( + cosmo, + prof_hod, + deltah=0.02, + log10Mh_min=12.0, + log10Mh_max=15.9, + lk_arr=None, + a_arr=None, + khmin=khmin, +): + """Implements the response of galaxy-matter power spectrum + to the long wavelength modes + developed in Terasawa et al. 2023 (arXiv:2310.13330) as: + + .. math:: + \\frac{\\partial P_{gm}(k)}{\\partial\\delta_b} = + \\left(\\partial P_{gm}(k)}{\\partial\\delta_b}|_{G} - \\frac{1}{3} + \\frac{d\\log P_{gm}(k)}{d\\log k}\\right)P_{gm}(k), + + where the :math:`\\partial P_{gm}(k)}{\\partial\\delta_b}|_{G}` + is the growth response to the long wavelength modes. + + Args: + cosmo (:class:`~pyccl.core.Cosmology`): a Cosmology object. + prof_hod (:class:`~pyccl.halos.profiles.hod.HaloProfileHOD`): + HOD profile. + deltah (float): the variation of :math:`h` to compute + the growth reponse to :math:`h` + by the two-sided numerical derivative method. + log10Mh_min (float): the minimum halo mass + (in units of :math:`M_\\odot/h`) for integration. + log10Mh_max (float): the maximum halo mass + (in units of :math:`M_\\odot/h`) for integration. + lk_arr (array): an array holding values of the natural + logarithm of the wavenumber (in units of :math:`Mpc^{-1}`) at + which the response is calculated. + a_arr (array): an array holding values of the scale factor + at which the response is calculated. + khmin (float): the wavenumber (in units of :math:`h Mpc^{-1}`) below + which the response is calculated using the perturbation theory. + Returns: + Response of the galaxy-matter power spectrum. + """ + + # Set k and a sampling from CCL parameters + if lk_arr is None: + status = 0 + lk_arr, status = lib.get_pk_spline_lk( + cosmo.cosmo, lib.get_pk_spline_nk(cosmo.cosmo), status + ) + check(status, cosmo=cosmo) + if a_arr is None: + status = 0 + a_arr, status = lib.get_pk_spline_a( + cosmo.cosmo, lib.get_pk_spline_na(cosmo.cosmo), status + ) + check(status, cosmo=cosmo) + + k_use = np.exp(lk_arr) + + # Check inputs + if not isinstance(prof_hod, halos.profiles.HaloProfile): + raise TypeError("prof_hod must be of type `HaloProfile`") + + # dark emulator is valid for 0 =< z <= 1.48 + if np.any((1.0 / a_arr - 1) > 1.48): + raise ValueError("dark emulator is valid for z<=1.48") + + # dark emulator support range is 10^12 <= M200m <= 10^16 Msun/h + if log10Mh_min < 12.0 or log10Mh_max > 16.0: + raise ValueError( + "Input mass range is not supported." + "The supported range is from 10^12 to 10^16 Msun/h." + ) + + h = cosmo["h"] + k_emu = k_use / h # [h/Mpc] + cosmo.compute_linear_power() + pk2dlin = cosmo.get_linear_power("delta_matter:delta_matter") + + # initialize dark emulator class + emu = darkemu.de_interface.base_class() + + # set h-modified cosmology to take finite differencing + hp = h + deltah + hm = h - deltah + cosmo_hp, cosmo_hm = _set_hmodified_cosmology(cosmo, deltah) + + # Growth factor + Dp = cosmo_hp.growth_factor_unnorm(a_arr) + Dm = cosmo_hm.growth_factor_unnorm(a_arr) + + na = len(a_arr) + nk = len(k_use) + dpk12 = np.zeros([na, nk]) + nM = 2**5 + 1 + log10M_min = np.log10(10**log10Mh_min / h) + log10M_max = np.log10(10**log10Mh_max / h) + + b1_th_tink = np.zeros(nM) + Pth = np.zeros((nM, nk)) + Pnth_hp = np.zeros((nM, nk)) + Pnth_hm = np.zeros((nM, nk)) + Pbin = np.zeros((nM, nk)) + nths = np.zeros(nM) + + mass_def = halos.MassDef200m + hmf = halos.MassFuncNishimichi19(mass_def=mass_def, extrapolate=True) + hbf = halos.HaloBiasTinker10(mass_def=mass_def) + hmc = halos.HMCalculator( + mass_function=hmf, + halo_bias=hbf, + mass_def=mass_def, + log10M_min=log10M_min, + log10M_max=log10M_max, + nM=nM, + ) + + logM = hmc._lmass + M = hmc._mass + Mh = M * h + dlogM = logM[1] - logM[0] + + for ia, aa in enumerate(a_arr): + z = 1.0 / aa - 1 + hmc._get_ingredients(cosmo, aa, get_bf=True) + + _darkemu_set_cosmology(emu, cosmo) + for m in range(nM): + hmc_m = halos.HMCalculator( + mass_function=hmf, + halo_bias=hbf, + mass_def=mass_def, + log10M_min=logM[m], + log10M_max=17.0, + ) + hmc_m._get_ingredients(cosmo, aa, get_bf=True) + + Pth[m] = emu.get_phm_massthreshold(k_emu, Mh[m], z) * (1 / h) ** 3 + Pbin[m] = emu.get_phm_mass(k_emu, Mh[m], z) * (1 / h) ** 3 + nths[m] = emu.mass_to_dens(Mh[m], z) * h**3 + + array_2 = np.ones(len(hmc_m._mass)) + array_2[..., 0] = 0 + b1_th_tink[m] = hmc_m._integrate_over_mbf( + array_2 + ) / hmc_m._integrate_over_mf(array_2) + + _darkemu_set_cosmology(emu, cosmo_hp) + for m in range(nM): + Pnth_hp[m] = ( + emu.get_phm( + k_emu * (h / hp), np.log10(nths[m] * (1 / hp) ** 3), z + ) + * (1 / hp) ** 3 + ) + + _darkemu_set_cosmology(emu, cosmo_hm) + for m in range(nM): + Pnth_hm[m] = ( + emu.get_phm( + k_emu * (h / hm), np.log10(nths[m] * (1 / hm) ** 3), z + ) + * (1 / hm) ** 3 + ) + + Nc = prof_hod._Nc(M, aa) + Ns = prof_hod._Ns(M, aa) + fc = prof_hod._fc(aa) + Ng = Nc * (fc + Ns) + logMps = logM + dlogM + logMms = logM - dlogM + + prof_Mp = prof_hod.fourier(cosmo, k_use, (10**logMps), aa) + prof_Mm = prof_hod.fourier(cosmo, k_use, (10**logMms), aa) + + dprof_dlogM = (prof_Mp - prof_Mm) / (2 * dlogM) + + nth_mat = np.tile(nths, (len(k_use), 1)).transpose() + + # Eq. 18 + ng = hmc._integrate_over_mf(Ng) + + # Eq. 17 + bgE = hmc._integrate_over_mbf(Ng) / ng + + # Eq. 19 + bgE2 = hmc._integrate_over_mf(Ng * _b2H17(hmc._bf)) / ng + + bgL = bgE - 1 + + b1L_th_mat = np.tile(b1_th_tink - 1, (len(k_emu), 1)).transpose() + + dPhm_db_nfix = ( + (26.0 / 21.0) + * np.log(np.array(Pnth_hp) / np.array(Pnth_hm)) + * np.array(Pth) + / (2 * (np.log(Dp[ia]) - np.log(Dm[ia]))) + ) # Mpc^3 + + dnP_hm_db_emu = nth_mat * (dPhm_db_nfix + b1L_th_mat * np.array(Pbin)) + + # Eq. A2 + nP = nth_mat * np.array(Pth) + + # Eq. A7 + Pgm = integrate.romb(dprof_dlogM * nP, dx=dlogM, axis=0) / ng + + # The first term of Eq. A8 + dnP_gm_db = integrate.romb( + dprof_dlogM * (dnP_hm_db_emu), dx=dlogM, axis=0 + ) + + # Here we assume the galaxies' profile around host halo + # is fixed at physical corrdinate. + dprof_dlogM_dlogk = np.zeros((nM, nk)) + for i in range(nM): + dprof_dlogM_dlogk[i] = np.gradient((dprof_dlogM[i])) / np.gradient( + np.log(k_use) + ) + + # The second term of Eq. A8 + G_prof = ( + +1.0 + / 3.0 + * integrate.romb(dprof_dlogM_dlogk * nP, dx=dlogM, axis=0) + ) + + # Eq. 25 + Pgm_growth = (dnP_gm_db + G_prof) / ng - bgL * Pgm + + Pgm_d = ( + -1.0 + / 3.0 + * np.gradient(np.log(Pgm)) + / np.gradient(np.log(k_use)) + * Pgm + ) + + # Eq. 22 + dPgm_db_emu = Pgm_growth + Pgm_d + + dpklin = pk2dlin(k_use, aa, cosmo, derivative=True) + + # Eq. 16 + dPgm_db_lin = ( + (47 / 21 + bgE2 / bgE - bgE - 1 / 3 * dpklin) + * bgE + * pk2dlin(k_use, aa, cosmo) + ) + + # stitching + k_switch = 0.08 # [h/Mpc] + + # Eq. 27 + dPgm_db = dPgm_db_lin * np.exp(-k_emu / k_switch) + dPgm_db_emu * ( + 1 - np.exp(-k_emu / k_switch) + ) + + # use the perturbation theory below khmin + dPgm_db[k_emu < khmin] = dPgm_db_lin[k_emu < khmin] + dpk12[ia, :] = dPgm_db + + return dpk12 + + +def resp_Pgg_darkemu( + cosmo, + prof_hod, + deltalnAs=0.03, + log10Mh_min=12.0, + log10Mh_max=15.9, + lk_arr=None, + a_arr=None, + khmin=khmin, +): + """Implements the response of galaxy power spectrum to the long wavelength + modes developed in Terasawa et al. 2023 (arXiv:2310.13330) as: + + .. math:: + \\frac{\\partial P_{gg}(k)}{\\partial\\delta_b} = + \\left(-1 + \\partial P_{gg}(k)}{\\partial\\delta_b}|_{G} + - \\frac{1}{3} + \\frac{d\\log P_{gg}(k)}{d\\log k}\\right)P_{gg}(k), + + where the :math:`\\partial P_{gg}(k)}{\\partial\\delta_b}|_{G}` + is the growth response to the long wavelength modes. + + Args: + cosmo (:class:`~pyccl.core.Cosmology`): a Cosmology object. + prof_hod (:class:`~pyccl.halos.profiles.hod.HaloProfileHOD`): + HOD profile. + deltalnAs (float): the variation of :math:`\\ln A_{s}` + to compute the growth reponse to :math:`\\ln A_{s}` by + the two-sided numerical derivative method. + log10Mh_min (float): the minimum halo mass + (in units of :math:`M_\\odot/h`) for integration. + log10Mh_max (float): the maximum halo mass + (in units of :math:`M_\\odot/h`) for integration. + lk_arr (array): an array holding values of the natural + logarithm of the wavenumber (in units of :math:`Mpc^{-1}`) at + which the response is calculated. + a_arr (array): an array holding values of the scale factor + at which the response is calculated. + khmin (float): the wavenumber (in units of :math:`h Mpc^{-1}`) below + which the response is calculated using the perturbation theory. + Returns: + Response of the galaxy power spectrum. + """ + + # Set k and a sampling from CCL parameters + if lk_arr is None: + status = 0 + lk_arr, status = lib.get_pk_spline_lk( + cosmo.cosmo, lib.get_pk_spline_nk(cosmo.cosmo), status + ) + check(status, cosmo=cosmo) + if a_arr is None: + status = 0 + a_arr, status = lib.get_pk_spline_a( + cosmo.cosmo, lib.get_pk_spline_na(cosmo.cosmo), status + ) + check(status, cosmo=cosmo) + + k_use = np.exp(lk_arr) + + # Check inputs + if not isinstance(prof_hod, halos.profiles.HaloProfile): + raise TypeError("prof_hod must be of type `HaloProfile`") + + # dark emulator is valid for 0 =< z <= 1.48 + if np.any((1.0 / a_arr - 1) > 1.48): + raise ValueError("dark emulator is valid for z<=1.48") + + # dark emulator support range is 10^12 <= M200m <= 10^16 Msun/h + if log10Mh_min < 12.0 or log10Mh_max > 16.0: + raise ValueError( + "Input mass range is not supported." + "The supported range is from 10^12 to 10^16 Msun/h." + ) + + h = cosmo["h"] + k_emu = k_use / h # [h/Mpc] + cosmo.compute_linear_power() + pk2dlin = cosmo.get_linear_power("delta_matter:delta_matter") + + # initialize dark emulator class + emu = darkemu.de_interface.base_class() + + na = len(a_arr) + nk = len(k_use) + dpk12 = np.zeros([na, nk]) + nM = 2**5 + 1 + log10M_min = np.log10(10**log10Mh_min / h) + log10M_max = np.log10(10**log10Mh_max / h) + + b1_th_tink = np.zeros(nM) + Pth = np.zeros((nM, nM, nk)) + Pth_Ap = np.zeros((nM, nM, nk)) + Pth_Am = np.zeros((nM, nM, nk)) + Pth_bin = np.zeros((nM, nM, nk)) + nths = np.zeros(nM) + + mass_def = halos.MassDef200m + hmf = halos.MassFuncNishimichi19(mass_def=mass_def, extrapolate=True) + hbf = halos.HaloBiasTinker10(mass_def=mass_def) + prof_2pt = halos.profiles_2pt.Profile2ptHOD() + + hmc = halos.HMCalculator( + mass_function=hmf, + halo_bias=hbf, + mass_def=mass_def, + log10M_min=log10M_min, + log10M_max=log10M_max, + nM=nM, + ) + + logM = hmc._lmass + M = hmc._mass + Mh = M * h + dlogM = logM[1] - logM[0] + + for ia, aa in enumerate(a_arr): + z = 1.0 / aa - 1 + hmc._get_ingredients(cosmo, aa, get_bf=True) + + for m in range(nM): + hmc_m = halos.HMCalculator( + mass_function=hmf, + halo_bias=hbf, + mass_def=mass_def, + log10M_min=logM[m], + log10M_max=17.0, + ) + hmc_m._get_ingredients(cosmo, aa, get_bf=True) + + nths[m] = emu.mass_to_dens(Mh[m], z) * h**3 + + array_2 = np.ones(len(hmc_m._mass)) + array_2[..., 0] = 0 + b1_th_tink[m] = hmc_m._integrate_over_mbf( + array_2 + ) / hmc_m._integrate_over_mf(array_2) + + # set cosmology for dark emulator + _darkemu_set_cosmology(emu, cosmo) + for m in range(nM): + for n in range(nM): + Pth[m, n] = ( + emu.get_phh( + k_emu, + np.log10(nths[m] / (h**3)), + np.log10(nths[n] / (h**3)), + z, + ) + * (1 / h) ** 3 + ) + + Pth_bin[m, n] = ( + _get_phh_massthreshold_mass( + emu, k_emu, nths[m] / (h**3), Mh[n], z + ) + * (1 / h) ** 3 + ) + + _darkemu_set_cosmology_forAsresp(emu, cosmo, deltalnAs) + for m in range(nM): + for n in range(nM): + Pth_Ap[m, n] = ( + emu.get_phh( + k_emu, + np.log10(nths[m] / (h**3)), + np.log10(nths[n] / (h**3)), + z, + ) + * (1 / h) ** 3 + ) + + _darkemu_set_cosmology_forAsresp(emu, cosmo, -deltalnAs) + for m in range(nM): + for n in range(nM): + Pth_Am[m, n] = ( + emu.get_phh( + k_emu, + np.log10(nths[m] / (h**3)), + np.log10(nths[n] / (h**3)), + z, + ) + * (1 / h) ** 3 + ) + + Nc = prof_hod._Nc(M, aa) + Ns = prof_hod._Ns(M, aa) + fc = prof_hod._fc(aa) + Ng = Nc * (fc + Ns) + logMps = logM + dlogM + logMms = logM - dlogM + + prof_Mp = prof_hod.fourier(cosmo, k_use, (10**logMps), aa) + prof_Mm = prof_hod.fourier(cosmo, k_use, (10**logMms), aa) + prof_1h = prof_2pt.fourier_2pt(cosmo, k_use, M, aa, prof_hod) + + dprof_dlogM = (prof_Mp - prof_Mm) / (2 * dlogM) + dprof_dlogM_dlogk = np.zeros((nM, nk)) + dprof_1h_dlogk = np.zeros((nM, nk)) + for i in range(nM): + dprof_dlogM_dlogk[i] = np.gradient((dprof_dlogM[i])) / np.gradient( + np.log(k_use) + ) + dprof_1h_dlogk[i] = np.gradient((prof_1h[i])) / np.gradient( + np.log(k_use) + ) + + nth_mat = np.tile(nths, (len(k_use), 1)).transpose() + + # Eq. 18 + ng = hmc._integrate_over_mf(Ng) + b1 = hbf(cosmo, M, aa) + + # Eq. 17 + bgE = hmc._integrate_over_mbf(Ng) / ng + + # Eq. 19 + bgE2 = hmc._integrate_over_mf(Ng * _b2H17(b1)) / ng + + bgL = bgE - 1 + b1L_mat = np.tile(b1 - 1, (len(k_emu), 1)).transpose() + b1L_th_mat = np.tile(b1_th_tink - 1, (len(k_emu), 1)).transpose() + + # P_gg(k) + _Pgg_1h = hmc._integrate_over_mf(prof_1h.T) / (ng**2) + + Pgg_2h_int = list() + for m in range(nM): + Pgg_2h_int.append( + integrate.romb( + Pth[m] * nth_mat * dprof_dlogM, axis=0, dx=dlogM + ) + ) + Pgg_2h_int = np.array(Pgg_2h_int) + + # Eq. A12 + _Pgg_2h = integrate.romb( + Pgg_2h_int * nth_mat * dprof_dlogM, axis=0, dx=dlogM + ) / (ng**2) + + Pgg = _Pgg_2h + _Pgg_1h + + # 2-halo response + dPhh_db_nfix = (26.0 / 21.0) * (Pth_Ap - Pth_Am) / (2 * deltalnAs) + + resp_2h_int = list() + for m in range(nM): + dP_hh_db_tot = dPhh_db_nfix[m] + 2 * b1L_th_mat * Pth_bin[m] + resp_2h_int.append( + integrate.romb( + dP_hh_db_tot * nth_mat * dprof_dlogM, axis=0, dx=dlogM + ) + ) + resp_2h_int = np.array(resp_2h_int) + resp_2h = integrate.romb( + resp_2h_int * nth_mat * dprof_dlogM, axis=0, dx=dlogM + ) / (ng**2) + + # Here we assume the galaxies' profile around host halo + # is fixed at physical corrdinate. + G_prof = ( + +2.0 + / 3.0 + * integrate.romb( + resp_2h_int * nth_mat * dprof_dlogM_dlogk, axis=0, dx=dlogM + ) + / (ng**2) + ) + + # The first term of Eq. A14 + resp_2h = resp_2h + G_prof + + # 1-halo response + resp_1h = hmc._integrate_over_mf(b1L_mat.T * prof_1h.T) / (ng**2) + + # Here we assume the galaxies' profile around host halo + # is fixed at physical corrdinate. + G_prof = hmc._integrate_over_mf(dprof_1h_dlogk.T) / (ng**2) / 3.0 + + # The first term of Eq. A15 + resp_1h = resp_1h + G_prof + + Pgg_growth = (resp_1h + resp_2h) - 2 * bgL * Pgg + + Pgg_d = ( + -1.0 + / 3.0 + * np.gradient(np.log(Pgg)) + / np.gradient(np.log(k_use)) + * Pgg + ) + + # Eq. 22 + dPgg_db_emu = Pgg_growth + Pgg_d - Pgg + + dpklin = pk2dlin(k_use, aa, cosmo, derivative=True) + + Pgg_lin = bgE**2 * pk2dlin(k_use, aa, cosmo) + + # Eq. 16 + dPgg_db_lin = ( + 47 / 21 + 2 * bgE2 / bgE - 2 * bgE - 1 / 3 * dpklin + ) * Pgg_lin + # stitching + k_switch = 0.08 # [h/Mpc] + + # Eq. 27 + dPgg_db = dPgg_db_lin * np.exp(-k_emu / k_switch) + dPgg_db_emu * ( + 1 - np.exp(-k_emu / k_switch) + ) + + Pgg = Pgg_lin * np.exp(-k_emu / k_switch) + Pgg * ( + 1 - np.exp(-k_emu / k_switch) + ) + + # use the perturbation theory below khmin + dPgg_db[k_emu < khmin] = dPgg_db_lin[k_emu < khmin] + dpk12[ia, :] = dPgg_db + + return dpk12 + + +# Utility functions #################### +def _get_phh_massthreshold_mass(emu, k_emu, dens1, Mbin, redshift): + """Compute the halo-halo power spectrum between + mass bin halo sample and mass threshold halo sample + specified by the corresponding cumulative number density. + """ + M2p = Mbin * 1.01 + M2m = Mbin * 0.99 + dens2p = emu.mass_to_dens(M2p, redshift) + dens2m = emu.mass_to_dens(M2m, redshift) + logdens1, logdens2p, logdens2m = ( + np.log10(dens1), + np.log10(dens2p), + np.log10(dens2m), + ) + + p_p = emu.get_phh(k_emu, logdens1, logdens2p, redshift) + p_m = emu.get_phh(k_emu, logdens1, logdens2m, redshift) + + numer = dens2m * p_m - dens2p * p_p + denom = dens2m - dens2p + + return numer / denom + + +def _b2H17(b1): + """Implements fitting formula for secondary halo bias, b_2, described in + arXiv:1607.01024. + """ + b2 = 0.77 - (2.43 * b1) + (b1 * b1) + return b2 + + +def _darkemu_set_cosmology(emu, cosmo): + """Input cosmology and initiallize the base class of DarkEmulator.""" + h = cosmo["h"] + n_s = cosmo["n_s"] + A_s = cosmo["A_s"] + if np.isnan(A_s): + raise ValueError("A_s must be provided to use the Dark Emulator") + + omega_c = cosmo["Omega_c"] * h**2 + omega_b = cosmo["Omega_b"] * h**2 + omega_nu = 0.00064 # we fix this value (Nishimichi et al. 2019) + Omega_L = 1 - ((omega_c + omega_b + omega_nu) / h**2) + + # Parameters cparam (numpy array) : Cosmological parameters + # (omega_b,omega_c,Omega_de,ln(10^10As),ns,w) + cparam = np.array( + [omega_b, omega_c, Omega_L, np.log(10**10 * A_s), n_s, -1.0] + ) + if darkemu.cosmo_util.test_cosm_range(cparam): + raise ValueError( + ("cosmological parameter out of supported range of DarkEmulator") + ) + + emu.set_cosmology(cparam) + + +def _darkemu_set_cosmology_forAsresp(emu, cosmo, deltalnAs): + """Input cosmology and initiallize the base class of DarkEmulator + for cosmology with modified A_s. + """ + h = cosmo["h"] + n_s = cosmo["n_s"] + A_s = cosmo["A_s"] + if np.isnan(A_s): + raise ValueError("A_s must be provided to use the Dark Emulator") + + omega_c = cosmo["Omega_c"] * h**2 + omega_b = cosmo["Omega_b"] * h**2 + omega_nu = 0.00064 # we fix this value (Nishimichi et al. 2019) + Omega_L = 1 - ((omega_c + omega_b + omega_nu) / h**2) + + # Parameters cparam (numpy array) : Cosmological parameters + # (omega_b,omega_c,Omega_de,ln(10^10As),ns,w) + cparam = np.array( + [ + omega_b, + omega_c, + Omega_L, + np.log(10**10 * A_s) + deltalnAs, + n_s, + -1.0, + ] + ) + if darkemu.cosmo_util.test_cosm_range(cparam): + raise ValueError( + ("cosmological parameter out of supported range of DarkEmulator") + ) + + emu.set_cosmology(cparam) + + return emu + + +def _set_hmodified_cosmology(cosmo, deltah, extra_parameters=None): + """Create the Cosmology objects with modified Hubble parameter h.""" + Omega_c = cosmo["Omega_c"] + Omega_b = cosmo["Omega_b"] + h = cosmo["h"] + + cosmo_hmodified = [] + for i in [+1, -1]: + hp = h + i * deltah + + # \Omega_c h^2, \Omega_b h^2 is fixed + Omega_c_p = np.power((h / hp), 2) * Omega_c + Omega_b_p = np.power((h / hp), 2) * Omega_b + + cosmo_hp_dict = cosmo.to_dict() + cosmo_hp_dict["h"] = hp + cosmo_hp_dict["Omega_c"] = Omega_c_p + cosmo_hp_dict["Omega_b"] = Omega_b_p + cosmo_hp_dict["extra_parameters"] = extra_parameters + cosmo_hp = cosmology.Cosmology(**cosmo_hp_dict) + cosmo_hmodified.append(cosmo_hp) + + return cosmo_hmodified[0], cosmo_hmodified[1] diff --git a/pyccl/tests/test_pk_response.py b/pyccl/tests/test_pk_response.py new file mode 100644 index 000000000..39b9fc9d8 --- /dev/null +++ b/pyccl/tests/test_pk_response.py @@ -0,0 +1,232 @@ +import numpy as np +import pyccl as ccl +import pytest +from dark_emulator import darkemu +from pyccl.pk_response import ( + resp_Pmm_hresponse, + resp_Pgm_darkemu, + resp_Pgg_darkemu, + _get_phh_massthreshold_mass, + _b2H17, + _darkemu_set_cosmology, + _darkemu_set_cosmology_forAsresp, + _set_hmodified_cosmology, +) +from .test_cclobject import check_eq_repr_hash + +# Set cosmology +Om = 0.3156 +ob = 0.02225 +oc = 0.1198 +OL = 0.6844 +As = 2.2065e-9 +ns = 0.9645 +h = 0.6727 + +cosmo = ccl.Cosmology( + Omega_c=oc / (h**2), + Omega_b=ob / (h**2), + h=h, + n_s=ns, + A_s=As, + m_nu=0.06, + transfer_function="boltzmann_camb", + extra_parameters={"camb": {"halofit_version": "takahashi"}}, +) + +deltah = 0.02 +deltalnAs = 0.03 +lk_arr = np.log(np.geomspace(1e-3, 1e1, 100)) +k_use = np.exp(lk_arr) +k_emu = k_use / h # [h/Mpc] +a_arr = np.array([1.0]) + + +# HOD parameters +logMhmin = 13.94 +logMh1 = 14.46 +alpha = 1.192 +kappa = 0.60 +sigma_logM = 0.5 +sigma_lM = sigma_logM * np.log(10) +logMh0 = logMhmin + np.log10(kappa) + +logMmin = np.log10(10**logMhmin / h) +logM0 = np.log10(10**logMh0 / h) +logM1 = np.log10(10**logMh1 / h) + +# Construct HOD +mass_def = ccl.halos.MassDef(200, "matter") +cm = ccl.halos.ConcentrationDuffy08(mass_def=mass_def) +prof_hod = ccl.halos.HaloProfileHOD( + mass_def=mass_def, + concentration=cm, + log10Mmin_0=logMmin, + siglnM_0=sigma_lM, + log10M0_0=logM0, + log10M1_0=logM1, + alpha_0=alpha, +) + +# initialize dark emulator class +emu = darkemu.de_interface.base_class() + +cosmo.compute_linear_power() +pk2dlin = cosmo.get_linear_power("delta_matter:delta_matter") + + +def test_resp_Pgm_darkemu_init(): + with pytest.raises(TypeError): + resp_Pgm_darkemu( + cosmo, None, deltah=deltah, lk_arr=lk_arr, a_arr=a_arr + ) + + # dark emulator is valid for z<=1.48 + with pytest.raises(ValueError): + resp_Pgm_darkemu( + cosmo, + prof_hod, + deltah=deltah, + lk_arr=lk_arr, + a_arr=np.array([0.3]), + ) + + # dark emulator support range is 10^12 <= M200m <= 10^16 Msun/h + with pytest.raises(ValueError): + resp_Pgm_darkemu( + cosmo, + prof_hod, + deltah=deltah, + lk_arr=lk_arr, + a_arr=a_arr, + log10Mh_min=11.9, + log10Mh_max=16.1, + ) + + +def test_resp_Pgg_darkemu_init(): + with pytest.raises(TypeError): + resp_Pgg_darkemu( + cosmo, None, deltalnAs=deltalnAs, lk_arr=lk_arr, a_arr=a_arr + ) + + # dark emulator is valid for z<=1.48 + with pytest.raises(ValueError): + resp_Pgg_darkemu( + cosmo, + prof_hod, + deltalnAs=deltalnAs, + lk_arr=lk_arr, + a_arr=np.array([0.3]), + ) + + # dark emulator support range is 10^12 <= M200m <= 10^16 Msun/h + with pytest.raises(ValueError): + resp_Pgg_darkemu( + cosmo, + prof_hod, + deltalnAs=deltalnAs, + lk_arr=lk_arr, + a_arr=a_arr, + log10Mh_min=11.9, + log10Mh_max=16.1, + ) + + +def test_resp_Pmm_hresponse(): + response = resp_Pmm_hresponse( + cosmo, deltah=deltah, lk_arr=lk_arr, a_arr=a_arr + ) + valid = (k_emu > 1e-2) & (k_emu < 4) + + assert np.all(response[0][valid] > 0) + + +def test_resp_Pgm_darkemu(): + response = resp_Pgm_darkemu( + cosmo, prof_hod, deltah=deltah, lk_arr=lk_arr, a_arr=a_arr + ) + valid = (k_emu > 1e-2) & (k_emu < 4) + + assert np.all(response[0][valid] > 0) + + +def test_resp_Pgg_darkemu(): + response = resp_Pgg_darkemu( + cosmo, prof_hod, deltalnAs=deltalnAs, lk_arr=lk_arr, a_arr=a_arr + ) + valid = (k_emu > 1e-2) & (k_emu < 4) + + assert np.all(response[0][valid] < 0) + + +# Tests for the utility functions +def test_get_phh_massthreshold_mass(): + # set cosmology for dark emulator + _darkemu_set_cosmology(emu, cosmo) + dens1 = 1e-3 + Mbin = 1e13 + redshift = 0.0 + phh = _get_phh_massthreshold_mass(emu, k_emu, dens1, Mbin, redshift) + pklin = pk2dlin(k_use, 1.0, cosmo) * h**3 + + valid = (k_emu > 1e-3) & (k_emu < 1e-2) + + # phh is power spectrum of biased tracer + assert np.all(phh[valid] > pklin[valid]) + + +def test_b2H17(): + b2 = _b2H17(0.0) + + assert b2 == 0.77 + + +def test_darkemu_set_cosmology(): + # Cosmo parameters out of bounds + cosmo_wr = ccl.Cosmology( + Omega_c=0.25, Omega_b=0.05, h=0.67, A_s=2.2e-9, n_s=2.0 + ) + with pytest.raises(ValueError): + _darkemu_set_cosmology(emu, cosmo_wr) + + # A_s must be provided + cosmo_noAs = ccl.Cosmology( + Omega_c=0.25, Omega_b=0.05, h=0.67, sigma8=0.8, n_s=0.96 + ) + with pytest.raises(ValueError): + _darkemu_set_cosmology(emu, cosmo_noAs) + + +def test_darkemu_set_cosmology_forAsresp(): + # Cosmo parameters out of bounds + cosmo_wr = ccl.Cosmology( + Omega_c=0.25, Omega_b=0.05, h=0.67, A_s=2.2e-9, n_s=2.0 + ) + with pytest.raises(ValueError): + _darkemu_set_cosmology_forAsresp(emu, cosmo_wr, deltalnAs=100.0) + + # A_s must be provided + cosmo_noAs = ccl.Cosmology( + Omega_c=0.25, Omega_b=0.05, h=0.67, sigma8=0.8, n_s=0.96 + ) + with pytest.raises(ValueError): + _darkemu_set_cosmology(emu, cosmo_noAs) + + +def test_set_hmodified_cosmology(): + cosmo_hp, cosmo_hm = _set_hmodified_cosmology( + cosmo, deltah, extra_parameters=None + ) + + for cosmo_test in [cosmo_hp, cosmo_hm]: + cosmo_dict = cosmo_test.to_dict() + cosmo_dict["h"] = cosmo["h"] + cosmo_dict["Omega_c"] = cosmo["Omega_c"] + cosmo_dict["Omega_b"] = cosmo["Omega_b"] + cosmo_dict["extra_parameters"] = cosmo["extra_parameters"] + cosmo_new = ccl.Cosmology(**cosmo_dict) + + # make sure the output cosmologies are exactly the same as + # the input one, except for the modified parameters. + assert check_eq_repr_hash(cosmo, cosmo_new)