From 103d3a525521edc864e03f035ed0b484763836d7 Mon Sep 17 00:00:00 2001 From: Keitaro Ishikawa Date: Tue, 8 Jul 2025 23:29:15 +0100 Subject: [PATCH 1/7] added dark emulator pkgg/gm based on pk_2pt --- pyccl/halos/__init__.py | 1 + pyccl/halos/emulators/__init__.py | 1 + pyccl/halos/emulators/darkemu1_pk.py | 241 +++++++++++++++++++++++++++ 3 files changed, 243 insertions(+) create mode 100644 pyccl/halos/emulators/__init__.py create mode 100644 pyccl/halos/emulators/darkemu1_pk.py diff --git a/pyccl/halos/__init__.py b/pyccl/halos/__init__.py index 7d8a75b65..2be15d97a 100644 --- a/pyccl/halos/__init__.py +++ b/pyccl/halos/__init__.py @@ -9,3 +9,4 @@ from .pk_2pt import * from .pk_4pt import * from .halo_model import * +from .emulators import * diff --git a/pyccl/halos/emulators/__init__.py b/pyccl/halos/emulators/__init__.py new file mode 100644 index 000000000..d2e45e13d --- /dev/null +++ b/pyccl/halos/emulators/__init__.py @@ -0,0 +1 @@ +from .darkemu1_pk import * diff --git a/pyccl/halos/emulators/darkemu1_pk.py b/pyccl/halos/emulators/darkemu1_pk.py new file mode 100644 index 000000000..ee959a30e --- /dev/null +++ b/pyccl/halos/emulators/darkemu1_pk.py @@ -0,0 +1,241 @@ +__all__ = ("darkemu_power_spectrum",) + +import warnings +import numpy as np +from scipy import integrate +from scipy.interpolate import RectBivariateSpline as rbs +from scipy.interpolate import InterpolatedUnivariateSpline as ius + +from ... import Pk2D +from .. import Profile2pt + +from dark_emulator import darkemu +from dark_emulator import model_hod +hod = model_hod.darkemu_x_hod() + +def _compute_logdens(Mh, redshift): + Mh = np.atleast_1d(Mh) + logdens = np.log10([hod._convert_mass_to_dens( + Mh[i], redshift, integration=hod.config["hmf_int_algorithm"]) for i in range(len(Mh))]) + return logdens + +def _compute_p_hh(ks, Mh, redshift): # referred from dark_emulator_public/dark_emulator/model_hod/hod_interface.py + hod._compute_p_hh_spl(redshift) + + logdens_de = hod.g1.logdens_list + logdens = _compute_logdens(Mh, redshift) + + p_hh_tmp = np.zeros( + (len(logdens_de), len(Mh), len(ks))) + for i in range(len(logdens_de)): + p_hh_tmp[i] = rbs(-logdens_de, hod.fftlog_2h.k, + hod.p_hh_base[i])(-logdens, ks) + + p_hh = np.zeros((len(Mh), len(Mh), len(ks))) + for i in range(len(Mh)): + p_hh[:, i] = rbs(-logdens_de, hod.fftlog_2h.k, + p_hh_tmp[:, i, :])(-logdens, ks) + return p_hh.transpose(2, 0, 1) + + +def _compute_p_hm(ks, Mh, redshift): # referred from dark_emulator_public/dark_emulator/model_hod/hod_interface.py + if hod.logdens_computed == False: + hod._compute_logdens(redshift) + if hod.xi_hm_computed == False: + hod._compute_xi_hm(redshift) + + logdens = _compute_logdens(Mh, redshift) + + if hod.do_linear_correction: + pm_lin = hod.get_pklin(ks) + + p_hm = np.zeros((len(logdens), len(ks))) + for i in range(len(logdens)): + p_hm[i] = hod.fftlog_1h.xi2pk(hod.xi_hm[i], 1.01, N_extrap_high=1024)[1] + if hod.do_linear_correction: + p_hm[i] *= (hod.pm_lin_k_1h_out_of_range/pm_lin) + + return p_hm.T + + +def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, + prof2=None, prof_2pt=None, + p_of_k_a=None, + get_1h=True, get_2h=True, + suppress_1h=None, + extrap_pk=False): # KI add demu + """ Computes the halo model power spectrum for two + quantities defined by their respective halo profiles. + The halo model power spectrum for two profiles + :math:`u` and :math:`v` is: + + .. math:: + P_{u,v}(k,a) = I^0_2(k,a|u,v) + + I^1_1(k,a|u)\\,I^1_1(k,a|v)\\,P_{\\rm lin}(k,a) + + where :math:`P_{\\rm lin}(k,a)` is the linear matter + power spectrum, :math:`I^1_1` is defined in the documentation + of :meth:`~pyccl.halos.halo_model.HMCalculator.I_1_1`, and :math:`I^0_2` + is defined in the documentation of + :meth:`~pyccl.halos.halo_model.HMCalculator.I_0_2`. + + Args: + demu (:class:`dark_emulator.darkemu.base_class`): the dark emulator object. + cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object. + hmc (:class:`~pyccl.halos.halo_model.HMCalculator`): a halo model calculator. + k (:obj:`float` or `array`): comoving wavenumber in Mpc^-1. + a (:obj:`float` or `array`): scale factor. + prof (:class:`~pyccl.halos.profiles.profile_base.HaloProfile`): halo + profile. + prof2 (:class:`~pyccl.halos.profiles.profile_base.HaloProfile`): a + second halo profile. If ``None``, ``prof`` will be used as + ``prof2``. + prof_2pt (:class:`~pyccl.halos.profiles_2pt.Profile2pt`): + a profile covariance object + returning the the two-point moment of the two profiles + being correlated. If ``None``, the default second moment + will be used, corresponding to the products of the means + of both profiles. + p_of_k_a (:class:`~pyccl.pk2d.Pk2D`): a `Pk2D` object to + be used as the linear matter power spectrum. If ``None``, + the power spectrum stored within `cosmo` will be used. + get_1h (:obj:`bool`): if ``False``, the 1-halo term (i.e. the first + term in the first equation above) won't be computed. + get_2h (:obj:`bool`): if ``False``, the 2-halo term (i.e. the second + term in the first equation above) won't be computed. + suppress_1h (:obj:`callable` or :obj:`None`): + Suppress the 1-halo large scale contribution by a + time- and scale-dependent function :math:`k_*(a)`, + defined as in `HMCODE-2020 `_: + :math:`1/[1+(k_*(a)/k)^4]`. + If ``None`` the standard 1-halo term is returned with no damping. + extrap_pk (:obj:`bool`): + Whether to extrapolate ``p_of_k_a`` in case ``a`` is out of its + support. If ```False```, and the queried values are out of bounds, + an error is raised. + + Returns: + (:obj:`float` or `array`): integral values evaluated at each + combination of ``k`` and ``a``. The shape of the output will + be ``(N_a, N_k)`` where ``N_k`` and ``N_a`` are the sizes of + ``k`` and ``a`` respectively. If ``k`` or ``a`` are scalars, the + corresponding dimension will be squeezed out on output. + """ # noqa + a_use = np.atleast_1d(a).astype(float) + k_use = np.atleast_1d(k).astype(float) + + good_a = (1./a_use - 1. <= 1.48) & (1./a_use - 1. >= 0.) + if np.any(good_a == False): + warnings.warn('Dark Emulator I supports 0 <= z <= 1.48. The redshifts rather than 1.48 will be truncated.') + a_use = a_use[good_a] + + # KI add: cosmology set for the dark emulator + h_ = cosmo['h'] + if np.isnan(cosmo['A_s']): + raise ValueError("Dark Emulator needs A_s as input parameters.") + # [omegab0.,omega_c0.,Omega_Lambda0,np.log(As*10**10),ns,w0] + demu.set_cosmology(np.array([cosmo['Omega_b']*h_**2., + cosmo['Omega_c']*h_**2., + 1-(cosmo['Omega_c']+cosmo['Omega_b']), + np.log(cosmo['A_s']*10**10), + cosmo['n_s'], + cosmo['w0']])) + + if suppress_1h is not None: + if not get_1h: + raise ValueError("Can't suppress the 1-halo term " + "when get_1h is False.") + + if prof2 is None: + prof2 = prof + if prof_2pt is None: + prof_2pt = Profile2pt() + demu_phh_mass = demu.get_phh_mass + demu_phm_mass = demu.get_phm_mass + + # pk2d = cosmo.parse_pk(p_of_k_a) + extrap = cosmo if extrap_pk else None # extrapolation rule for pk2d + + na = len(a_use) + nk = len(k_use) + out = np.zeros([na, nk]) + M_array = hmc._mass.copy() + good_m = (M_array <= 1E16 / h_) & (M_array >= 1E12 / h_) + M_array = M_array[good_m] + if np.any(good_m == False): + warnings.warn('Dark Emulator I supports 1E12 <= Mh [Msun/h] <= 1E16. Others will be truncated.') + nM = len(M_array) + print('nM', nM) + lmass = np.log10(M_array) + for ia, aa in enumerate(a_use): + + # normalizations + norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) + if prof2 == prof: + norm2 = norm1 + else: + norm2 = prof2.get_normalization(cosmo, aa, hmc=hmc) + + # call components: hmf, uk, + # hmc._mf: dn/dlogM [Mpc^-3] (comoving) not dn/dM + hmc._get_ingredients(cosmo, aa, get_bf=False) + hmf = ius(hmc._lmass, hmc._mf)(lmass) + u1k = prof.fourier(cosmo, k_use, M_array, aa).T + + # calc pkhx + if prof2 == prof: # pkgg + u2k = u1k + # else: + # u2k = prof2.fourier(cosmo, k_use, M_array, aa).T + + ######### + # pkhh + # Dark Emulator accepts k [h/Mpc] and returns Pk [Mpc/h]^3 as units + pkhh = np.zeros([nk, nM, nM]) + for i in range(nM): + for j in range(i+1): + vals = demu_phh_mass(ks=k_use/h_, + M1=M_array[i]*h_,M2=M_array[j]*h_, + redshift=1./aa-1.) / h_**3. + pkhh[:,i,j] = vals + pkhh[:,j,i] = vals + # pkhh = _compute_p_hh(k_use, M_array, redshift=1./aa-1) # (Nk,NM,NM) + + # integration + pk_2h_M2_int = list() + for i in range(nM): + pk_2h_M2_int.append(hmc._integrator( + pkhh[:,i] * u1k * hmf[None,:], lmass)) + pk_2h = hmc._integrator(np.array(pk_2h_M2_int).T * u2k * hmf[None,:], lmass) + ######### + + else: + pkhm = np.zeros([nk, nM]) + for i in range(nM): + pkhm[:,i] = demu_phm_mass(ks=k_use/h_, M=M_array[i]*h_, redshift=1./aa-1.) / h_**3. + # pkhm = _compute_p_hm(k_use, M_array, redshift=1./aa-1) # (Nk,NM) + + pk_2h = hmc._integrator(pkhm * u1k * hmf[None,:], lmass) + + if get_1h and (prof2 == prof): + pk_1h = hmc.I_0_2(cosmo, k_use, aa, prof, + prof2=prof2, prof_2pt=prof_2pt) # 1h term + + if suppress_1h is not None: + # large-scale damping of 1-halo term + ks = suppress_1h(aa) + pk_1h *= (k_use / ks)**4 / (1 + (k_use / ks)**4) + else: + pk_1h = 0 + + # smooth 1h/2h transition region + if prof2 == prof: + out[ia] = (pk_1h + pk_2h) / (norm1 * norm2) + else: + out[ia] = pk_2h / norm1 + + if np.ndim(a) == 0: + out = np.squeeze(out, axis=0) + if np.ndim(k) == 0: + out = np.squeeze(out, axis=-1) + return out \ No newline at end of file From 065e05c605d36487f5c2d9b229ecf54e33db643c Mon Sep 17 00:00:00 2001 From: Keitaro Ishikawa Date: Wed, 9 Jul 2025 13:35:57 +0100 Subject: [PATCH 2/7] Fixed flake8 error etc. --- pyccl/halos/emulators/darkemu1_pk.py | 99 ++++++++++++++-------------- 1 file changed, 50 insertions(+), 49 deletions(-) diff --git a/pyccl/halos/emulators/darkemu1_pk.py b/pyccl/halos/emulators/darkemu1_pk.py index ee959a30e..b046ccc33 100644 --- a/pyccl/halos/emulators/darkemu1_pk.py +++ b/pyccl/halos/emulators/darkemu1_pk.py @@ -2,24 +2,26 @@ import warnings import numpy as np -from scipy import integrate + from scipy.interpolate import RectBivariateSpline as rbs from scipy.interpolate import InterpolatedUnivariateSpline as ius -from ... import Pk2D from .. import Profile2pt -from dark_emulator import darkemu from dark_emulator import model_hod hod = model_hod.darkemu_x_hod() + def _compute_logdens(Mh, redshift): Mh = np.atleast_1d(Mh) logdens = np.log10([hod._convert_mass_to_dens( - Mh[i], redshift, integration=hod.config["hmf_int_algorithm"]) for i in range(len(Mh))]) + Mh[i], redshift, integration=hod.config["hmf_int_algorithm"]) + for i in range(len(Mh))]) return logdens -def _compute_p_hh(ks, Mh, redshift): # referred from dark_emulator_public/dark_emulator/model_hod/hod_interface.py + +# referred from dark_emulator_public/dark_emulator/model_hod/hod_interface.py +def _compute_p_hh(ks, Mh, redshift): hod._compute_p_hh_spl(redshift) logdens_de = hod.g1.logdens_list @@ -29,19 +31,20 @@ def _compute_p_hh(ks, Mh, redshift): # referred from dark_emulator_public/dark_e (len(logdens_de), len(Mh), len(ks))) for i in range(len(logdens_de)): p_hh_tmp[i] = rbs(-logdens_de, hod.fftlog_2h.k, - hod.p_hh_base[i])(-logdens, ks) + hod.p_hh_base[i])(-logdens, ks) p_hh = np.zeros((len(Mh), len(Mh), len(ks))) for i in range(len(Mh)): p_hh[:, i] = rbs(-logdens_de, hod.fftlog_2h.k, - p_hh_tmp[:, i, :])(-logdens, ks) + p_hh_tmp[:, i, :])(-logdens, ks) return p_hh.transpose(2, 0, 1) -def _compute_p_hm(ks, Mh, redshift): # referred from dark_emulator_public/dark_emulator/model_hod/hod_interface.py - if hod.logdens_computed == False: +# referred from dark_emulator_public/dark_emulator/model_hod/hod_interface.py +def _compute_p_hm(ks, Mh, redshift): + if not hod.logdens_computed: hod._compute_logdens(redshift) - if hod.xi_hm_computed == False: + if not hod.xi_hm_computed: hod._compute_xi_hm(redshift) logdens = _compute_logdens(Mh, redshift) @@ -63,7 +66,7 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, p_of_k_a=None, get_1h=True, get_2h=True, suppress_1h=None, - extrap_pk=False): # KI add demu + extrap_pk=False): """ Computes the halo model power spectrum for two quantities defined by their respective halo profiles. The halo model power spectrum for two profiles @@ -109,10 +112,6 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, defined as in `HMCODE-2020 `_: :math:`1/[1+(k_*(a)/k)^4]`. If ``None`` the standard 1-halo term is returned with no damping. - extrap_pk (:obj:`bool`): - Whether to extrapolate ``p_of_k_a`` in case ``a`` is out of its - support. If ```False```, and the queried values are out of bounds, - an error is raised. Returns: (:obj:`float` or `array`): integral values evaluated at each @@ -125,8 +124,9 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, k_use = np.atleast_1d(k).astype(float) good_a = (1./a_use - 1. <= 1.48) & (1./a_use - 1. >= 0.) - if np.any(good_a == False): - warnings.warn('Dark Emulator I supports 0 <= z <= 1.48. The redshifts rather than 1.48 will be truncated.') + if np.any(good_a is False): + warnings.warn('Dark Emulator I supports 0 <= z <= 1.48. ' + 'The redshifts rather than 1.48 will be truncated.') a_use = a_use[good_a] # KI add: cosmology set for the dark emulator @@ -140,35 +140,32 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, np.log(cosmo['A_s']*10**10), cosmo['n_s'], cosmo['w0']])) - + if suppress_1h is not None: if not get_1h: raise ValueError("Can't suppress the 1-halo term " "when get_1h is False.") - + if prof2 is None: prof2 = prof + demu_phh_mass = demu.get_phh_mass if prof_2pt is None: prof_2pt = Profile2pt() - demu_phh_mass = demu.get_phh_mass - demu_phm_mass = demu.get_phm_mass - - # pk2d = cosmo.parse_pk(p_of_k_a) - extrap = cosmo if extrap_pk else None # extrapolation rule for pk2d - + demu_phm_mass = demu.get_phm_mass + na = len(a_use) nk = len(k_use) out = np.zeros([na, nk]) M_array = hmc._mass.copy() good_m = (M_array <= 1E16 / h_) & (M_array >= 1E12 / h_) M_array = M_array[good_m] - if np.any(good_m == False): - warnings.warn('Dark Emulator I supports 1E12 <= Mh [Msun/h] <= 1E16. Others will be truncated.') + if np.any(good_m is False): + warnings.warn('Dark Emulator I supports 1E12 <= Mh [Msun/h] <= 1E16. ' + 'Others will be truncated.') nM = len(M_array) print('nM', nM) lmass = np.log10(M_array) for ia, aa in enumerate(a_use): - # normalizations norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) if prof2 == prof: @@ -176,18 +173,18 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, else: norm2 = prof2.get_normalization(cosmo, aa, hmc=hmc) - # call components: hmf, uk, - # hmc._mf: dn/dlogM [Mpc^-3] (comoving) not dn/dM + # call components: hmf, uk + # hmc._mf: dn/dlogM [Mpc^-3] (comoving) not dn/dM hmc._get_ingredients(cosmo, aa, get_bf=False) hmf = ius(hmc._lmass, hmc._mf)(lmass) u1k = prof.fourier(cosmo, k_use, M_array, aa).T - + # calc pkhx - if prof2 == prof: # pkgg + if prof2 == prof: # pkgg u2k = u1k # else: # u2k = prof2.fourier(cosmo, k_use, M_array, aa).T - + ######### # pkhh # Dark Emulator accepts k [h/Mpc] and returns Pk [Mpc/h]^3 as units @@ -195,47 +192,51 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, for i in range(nM): for j in range(i+1): vals = demu_phh_mass(ks=k_use/h_, - M1=M_array[i]*h_,M2=M_array[j]*h_, + M1=M_array[i]*h_, M2=M_array[j]*h_, redshift=1./aa-1.) / h_**3. - pkhh[:,i,j] = vals - pkhh[:,j,i] = vals - # pkhh = _compute_p_hh(k_use, M_array, redshift=1./aa-1) # (Nk,NM,NM) - + pkhh[:, i, j] = vals + pkhh[:, j, i] = vals + # pkhh = _compute_p_hh(k_use, M_array, redshift=1./aa-1) + # (Nk,NM,NM) + # integration pk_2h_M2_int = list() for i in range(nM): pk_2h_M2_int.append(hmc._integrator( - pkhh[:,i] * u1k * hmf[None,:], lmass)) - pk_2h = hmc._integrator(np.array(pk_2h_M2_int).T * u2k * hmf[None,:], lmass) + pkhh[:, i] * u1k * hmf[None, :], lmass)) + pk_2h = hmc._integrator(np.array(pk_2h_M2_int).T * + u2k * hmf[None, :], lmass) ######### else: pkhm = np.zeros([nk, nM]) for i in range(nM): - pkhm[:,i] = demu_phm_mass(ks=k_use/h_, M=M_array[i]*h_, redshift=1./aa-1.) / h_**3. - # pkhm = _compute_p_hm(k_use, M_array, redshift=1./aa-1) # (Nk,NM) - - pk_2h = hmc._integrator(pkhm * u1k * hmf[None,:], lmass) - + pkhm[:, i] = demu_phm_mass(ks=k_use/h_, M=M_array[i]*h_, + redshift=1./aa-1.) / h_**3. + # pkhm = _compute_p_hm(k_use, M_array, redshift=1./aa-1) + # (Nk,NM) + + pk_2h = hmc._integrator(pkhm * u1k * hmf[None, :], lmass) + if get_1h and (prof2 == prof): pk_1h = hmc.I_0_2(cosmo, k_use, aa, prof, prof2=prof2, prof_2pt=prof_2pt) # 1h term - + if suppress_1h is not None: # large-scale damping of 1-halo term ks = suppress_1h(aa) pk_1h *= (k_use / ks)**4 / (1 + (k_use / ks)**4) else: pk_1h = 0 - + # smooth 1h/2h transition region if prof2 == prof: out[ia] = (pk_1h + pk_2h) / (norm1 * norm2) else: out[ia] = pk_2h / norm1 - + if np.ndim(a) == 0: out = np.squeeze(out, axis=0) if np.ndim(k) == 0: out = np.squeeze(out, axis=-1) - return out \ No newline at end of file + return out From fa8a4139639edb43f4a5684a592f3fa8db566c6d Mon Sep 17 00:00:00 2001 From: Keitaro Ishikawa Date: Wed, 9 Jul 2025 13:39:05 +0100 Subject: [PATCH 3/7] Fixed flake8 error etc. --- pyccl/halos/emulators/darkemu1_pk.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pyccl/halos/emulators/darkemu1_pk.py b/pyccl/halos/emulators/darkemu1_pk.py index b046ccc33..a0411d89f 100644 --- a/pyccl/halos/emulators/darkemu1_pk.py +++ b/pyccl/halos/emulators/darkemu1_pk.py @@ -15,7 +15,7 @@ def _compute_logdens(Mh, redshift): Mh = np.atleast_1d(Mh) logdens = np.log10([hod._convert_mass_to_dens( - Mh[i], redshift, integration=hod.config["hmf_int_algorithm"]) + Mh[i], redshift, integration=hod.config["hmf_int_algorithm"]) for i in range(len(Mh))]) return logdens @@ -54,7 +54,8 @@ def _compute_p_hm(ks, Mh, redshift): p_hm = np.zeros((len(logdens), len(ks))) for i in range(len(logdens)): - p_hm[i] = hod.fftlog_1h.xi2pk(hod.xi_hm[i], 1.01, N_extrap_high=1024)[1] + p_hm[i] = hod.fftlog_1h.xi2pk(hod.xi_hm[i], + 1.01, N_extrap_high=1024)[1] if hod.do_linear_correction: p_hm[i] *= (hod.pm_lin_k_1h_out_of_range/pm_lin) @@ -204,14 +205,14 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, for i in range(nM): pk_2h_M2_int.append(hmc._integrator( pkhh[:, i] * u1k * hmf[None, :], lmass)) - pk_2h = hmc._integrator(np.array(pk_2h_M2_int).T * + pk_2h = hmc._integrator(np.array(pk_2h_M2_int).T * u2k * hmf[None, :], lmass) ######### else: pkhm = np.zeros([nk, nM]) for i in range(nM): - pkhm[:, i] = demu_phm_mass(ks=k_use/h_, M=M_array[i]*h_, + pkhm[:, i] = demu_phm_mass(ks=k_use/h_, M=M_array[i]*h_, redshift=1./aa-1.) / h_**3. # pkhm = _compute_p_hm(k_use, M_array, redshift=1./aa-1) # (Nk,NM) From 4682fab3252e4274d371eb8761f01c94153371b8 Mon Sep 17 00:00:00 2001 From: Keitaro Ishikawa Date: Wed, 9 Jul 2025 13:39:57 +0100 Subject: [PATCH 4/7] Fixed flake8 error etc. --- pyccl/halos/emulators/darkemu1_pk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyccl/halos/emulators/darkemu1_pk.py b/pyccl/halos/emulators/darkemu1_pk.py index a0411d89f..0a32f9a79 100644 --- a/pyccl/halos/emulators/darkemu1_pk.py +++ b/pyccl/halos/emulators/darkemu1_pk.py @@ -54,7 +54,7 @@ def _compute_p_hm(ks, Mh, redshift): p_hm = np.zeros((len(logdens), len(ks))) for i in range(len(logdens)): - p_hm[i] = hod.fftlog_1h.xi2pk(hod.xi_hm[i], + p_hm[i] = hod.fftlog_1h.xi2pk(hod.xi_hm[i], 1.01, N_extrap_high=1024)[1] if hod.do_linear_correction: p_hm[i] *= (hod.pm_lin_k_1h_out_of_range/pm_lin) From fc004aa58757bf046fb7aafc867d5f2f5b782e9e Mon Sep 17 00:00:00 2001 From: Keitaro Ishikawa Date: Mon, 4 Aug 2025 18:04:02 +0900 Subject: [PATCH 5/7] Confirmed this is consistent with Dark Emulator I output. --- pyccl/halos/emulators/darkemu1_pk.py | 252 +++++++++++++++++---------- 1 file changed, 159 insertions(+), 93 deletions(-) diff --git a/pyccl/halos/emulators/darkemu1_pk.py b/pyccl/halos/emulators/darkemu1_pk.py index 0a32f9a79..0a6f244b4 100644 --- a/pyccl/halos/emulators/darkemu1_pk.py +++ b/pyccl/halos/emulators/darkemu1_pk.py @@ -1,73 +1,145 @@ -__all__ = ("darkemu_power_spectrum",) +__all__ = ("darkemu_power_spectrum", "galaxy_bias", "galaxy_bias_DEmuxHOD", ) import warnings import numpy as np -from scipy.interpolate import RectBivariateSpline as rbs from scipy.interpolate import InterpolatedUnivariateSpline as ius +from scipy.interpolate import RectBivariateSpline as rbs +from scipy.interpolate import RegularGridInterpolator as rgi from .. import Profile2pt from dark_emulator import model_hod -hod = model_hod.darkemu_x_hod() - +demuhod = model_hod.darkemu_x_hod() + + +def galaxy_bias(cosmo, hmc, a, prof): + """ Computes the galaxy bias for a given halo profile.""" + Mh = np.logspace(12.2,15.0,2**5+1) / cosmo['h'] # [Msun] + lmass = np.log10(Mh) + NgalM = _Ngal(Mh, 1., profile=prof) + # hmc._mf: dn/dlogM [Mpc^-3] (comoving) not dn/dM + hmc._get_ingredients(cosmo, a, get_bf=False) + hmf = ius(hmc._lmass, hmc._mf)(lmass) + norm1 = hmc._integrator(hmf * NgalM, lmass) # integral of dn/dM + hbias = np.array([demuhod.get_bias_mass(Mi, redshift=1./a - 1)[0,0] + for Mi in Mh*cosmo['h']]) + galbias = hmc._integrator(hmf * hbias * NgalM, lmass) / norm1 + # integral of dn/dM + return galbias + + +def galaxy_bias_DEmuxHOD(cosmo, a, prof): + a_use = np.atleast_1d(a).astype(float) + good_a = (1./a_use - 1. <= 1.48) & (1./a_use - 1. >= 0.) + if np.any(good_a is False): + warnings.warn('Dark Emulator I supports 0 <= z <= 1.48. ' + 'The redshifts larger than 1.48 will be truncated.') + a_use = a_use[good_a] -def _compute_logdens(Mh, redshift): + # cosmological params for the dark emulator + h_ = cosmo['h'] + if np.isnan(cosmo['A_s']): + raise ValueError("Dark Emulator needs A_s as an input parameter.") + # [omegab0.,omega_c0.,Omega_Lambda0,np.log(As*10**10),ns,w0] + demuhod.set_cosmology(np.array([cosmo['Omega_b']*h_**2., + cosmo['Omega_c']*h_**2., + 1-(cosmo['Omega_c']+cosmo['Omega_b']), + np.log(cosmo['A_s']*10**10), + cosmo['n_s'], + cosmo['w0']])) + # HOD parameters for the dark emulator + gparam_input = {"logMmin": prof.log10Mmin_0 + np.log10(h_), + "sigma_sq": (prof.siglnM_0/np.log(10))**2, + "logM1": prof.log10M1_0 + np.log10(h_), + "alpha": prof.alpha_0, + "kappa": 10**(prof.log10M0_0 - prof.log10Mmin_0), + "poff": 0., "Roff": 0., + "sat_dist_type": "NFW", + "alpha_inc": 0., "logM_inc": 0.} + demuhod.set_galaxy(gparam_input) + + galbias = np.array([demuhod._get_effective_bias(redshift=z_) + for z_ in (1./a_use - 1.)]) + if len(galbias) == 1: + return galbias[0] + else: + return galbias + + +def _compute_logdens(dehod, Mh, redshift): Mh = np.atleast_1d(Mh) - logdens = np.log10([hod._convert_mass_to_dens( - Mh[i], redshift, integration=hod.config["hmf_int_algorithm"]) + logdens = np.log10([dehod._convert_mass_to_dens( + Mh[i], redshift, integration=dehod.config["hmf_int_algorithm"]) for i in range(len(Mh))]) return logdens -# referred from dark_emulator_public/dark_emulator/model_hod/hod_interface.py -def _compute_p_hh(ks, Mh, redshift): - hod._compute_p_hh_spl(redshift) - - logdens_de = hod.g1.logdens_list - logdens = _compute_logdens(Mh, redshift) - - p_hh_tmp = np.zeros( - (len(logdens_de), len(Mh), len(ks))) - for i in range(len(logdens_de)): - p_hh_tmp[i] = rbs(-logdens_de, hod.fftlog_2h.k, - hod.p_hh_base[i])(-logdens, ks) - - p_hh = np.zeros((len(Mh), len(Mh), len(ks))) - for i in range(len(Mh)): - p_hh[:, i] = rbs(-logdens_de, hod.fftlog_2h.k, - p_hh_tmp[:, i, :])(-logdens, ks) +def _compute_p_hh(dehod, ks, Mh, redshift): + if not dehod.logdens_computed: + dehod._compute_logdens(redshift) + dehod._compute_p_hh_spl(redshift) + + points_known = (-dehod.g1.logdens_list, -dehod.g1.logdens_list, dehod.fftlog_2h.k) + logdens = _compute_logdens(dehod=dehod, Mh=Mh, redshift=redshift) + grid_d1, grid_d2, grid_k = np.meshgrid(-logdens, -logdens, + ks, indexing='ij') + points_target = np.stack([grid_d1.ravel(), grid_d2.ravel(), + grid_k.ravel()], axis=-1) + + interpolator = rgi(points_known, dehod.p_hh_base, method='cubic', + bounds_error=False, fill_value=None) + + p_hh = interpolator(points_target).reshape(len(logdens), len(logdens), len(ks)) return p_hh.transpose(2, 0, 1) + # referred from dark_emulator_public/dark_emulator/model_hod/hod_interface.py -def _compute_p_hm(ks, Mh, redshift): - if not hod.logdens_computed: - hod._compute_logdens(redshift) - if not hod.xi_hm_computed: - hod._compute_xi_hm(redshift) +def _compute_p_hm(dehod, ks, Mh, redshift): + if not dehod.logdens_computed: + dehod._compute_logdens(redshift) + if not dehod.xi_hm_computed: + dehod._compute_xi_hm(redshift) - logdens = _compute_logdens(Mh, redshift) + logdens = _compute_logdens(dehod=dehod, Mh=Mh, redshift=redshift) - if hod.do_linear_correction: - pm_lin = hod.get_pklin(ks) + if dehod.do_linear_correction: + pm_lin = dehod.get_pklin(ks) - p_hm = np.zeros((len(logdens), len(ks))) - for i in range(len(logdens)): - p_hm[i] = hod.fftlog_1h.xi2pk(hod.xi_hm[i], - 1.01, N_extrap_high=1024)[1] - if hod.do_linear_correction: - p_hm[i] *= (hod.pm_lin_k_1h_out_of_range/pm_lin) + p_hm = np.zeros((len(dehod.logdens), len(dehod.fftlog_1h.k))) + for i in range(len(dehod.logdens)): + p_hm[i] = dehod.fftlog_1h.xi2pk(dehod.xi_hm[i], + 1.01, N_extrap_high=1024)[1] + if dehod.do_linear_correction: + p_hm[i] *= (dehod.pm_lin_k_1h_out_of_range/pm_lin) - return p_hm.T + return rbs(dehod.fftlog_1h.k, -dehod.logdens, p_hm.T)(ks, -logdens) -def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, +def _Ngal(M, a, profile): + Nc = profile._Nc(M, a) + Ns = profile._Ns(M, a) + fc = profile._fc(a) + if profile.ns_independent: + return Nc*fc + Ns + return Nc*(fc + Ns) + + +def _I_0_2(hmc, cosmo, k, a, prof, *, prof2=None, prof_2pt, lmass, hmf): + if prof2 is None: + prof2 = prof + + M_array = 10**lmass + hmc._check_mass_def(prof, prof2) + uk = prof_2pt.fourier_2pt(cosmo, k, M_array, a, prof, prof2=prof2).T + return hmc._integrator(hmf[None, :] * uk, lmass) + + +def darkemu_power_spectrum(cosmo, hmc, k, a, prof, *, prof2=None, prof_2pt=None, - p_of_k_a=None, get_1h=True, get_2h=True, - suppress_1h=None, - extrap_pk=False): + suppress_1h=None): """ Computes the halo model power spectrum for two quantities defined by their respective halo profiles. The halo model power spectrum for two profiles @@ -84,7 +156,6 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, :meth:`~pyccl.halos.halo_model.HMCalculator.I_0_2`. Args: - demu (:class:`dark_emulator.darkemu.base_class`): the dark emulator object. cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object. hmc (:class:`~pyccl.halos.halo_model.HMCalculator`): a halo model calculator. k (:obj:`float` or `array`): comoving wavenumber in Mpc^-1. @@ -100,9 +171,6 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, being correlated. If ``None``, the default second moment will be used, corresponding to the products of the means of both profiles. - p_of_k_a (:class:`~pyccl.pk2d.Pk2D`): a `Pk2D` object to - be used as the linear matter power spectrum. If ``None``, - the power spectrum stored within `cosmo` will be used. get_1h (:obj:`bool`): if ``False``, the 1-halo term (i.e. the first term in the first equation above) won't be computed. get_2h (:obj:`bool`): if ``False``, the 2-halo term (i.e. the second @@ -127,20 +195,30 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, good_a = (1./a_use - 1. <= 1.48) & (1./a_use - 1. >= 0.) if np.any(good_a is False): warnings.warn('Dark Emulator I supports 0 <= z <= 1.48. ' - 'The redshifts rather than 1.48 will be truncated.') + 'The redshifts larger than 1.48 will be truncated.') a_use = a_use[good_a] - # KI add: cosmology set for the dark emulator + # cosmological params for the dark emulator h_ = cosmo['h'] if np.isnan(cosmo['A_s']): - raise ValueError("Dark Emulator needs A_s as input parameters.") + raise ValueError("Dark Emulator needs A_s as an input parameter.") # [omegab0.,omega_c0.,Omega_Lambda0,np.log(As*10**10),ns,w0] - demu.set_cosmology(np.array([cosmo['Omega_b']*h_**2., + demuhod.set_cosmology(np.array([cosmo['Omega_b']*h_**2., cosmo['Omega_c']*h_**2., 1-(cosmo['Omega_c']+cosmo['Omega_b']), np.log(cosmo['A_s']*10**10), cosmo['n_s'], cosmo['w0']])) + # HOD parameters for the dark emulator + gparam_input = {"logMmin": prof.log10Mmin_0 + np.log10(h_), + "sigma_sq": (prof.siglnM_0/np.log(10))**2, + "logM1": prof.log10M1_0 + np.log10(h_), + "alpha": prof.alpha_0, + "kappa": 10**(prof.log10M0_0 - prof.log10Mmin_0), + "poff": 0., "Roff": 0., + "sat_dist_type": "NFW", + "alpha_inc": 0., "logM_inc": 0.} + demuhod.set_galaxy(gparam_input) if suppress_1h is not None: if not get_1h: @@ -149,80 +227,67 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, if prof2 is None: prof2 = prof - demu_phh_mass = demu.get_phh_mass if prof_2pt is None: prof_2pt = Profile2pt() - demu_phm_mass = demu.get_phm_mass na = len(a_use) nk = len(k_use) out = np.zeros([na, nk]) M_array = hmc._mass.copy() good_m = (M_array <= 1E16 / h_) & (M_array >= 1E12 / h_) - M_array = M_array[good_m] if np.any(good_m is False): warnings.warn('Dark Emulator I supports 1E12 <= Mh [Msun/h] <= 1E16. ' 'Others will be truncated.') - nM = len(M_array) - print('nM', nM) + M_array = M_array[good_m] lmass = np.log10(M_array) + nM = len(M_array) + + for ia, aa in enumerate(a_use): # normalizations - norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) - if prof2 == prof: - norm2 = norm1 - else: - norm2 = prof2.get_normalization(cosmo, aa, hmc=hmc) - - # call components: hmf, uk + # norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) + # redshift dependence is not consistent with the original darkemulator. + NgalM = _Ngal(M_array, 1., profile=prof) # hmc._mf: dn/dlogM [Mpc^-3] (comoving) not dn/dM hmc._get_ingredients(cosmo, aa, get_bf=False) hmf = ius(hmc._lmass, hmc._mf)(lmass) - u1k = prof.fourier(cosmo, k_use, M_array, aa).T + norm1 = hmc._integrator(hmf * NgalM, lmass) # integral of dn/dM + if prof2 == prof: + norm2 = norm1 + - # calc pkhx + # calc profile in fourier space + u1k = prof.fourier(cosmo, k_use, M_array, aa).T # (Nk,NM) if prof2 == prof: # pkgg u2k = u1k - # else: - # u2k = prof2.fourier(cosmo, k_use, M_array, aa).T - - ######### - # pkhh - # Dark Emulator accepts k [h/Mpc] and returns Pk [Mpc/h]^3 as units - pkhh = np.zeros([nk, nM, nM]) - for i in range(nM): - for j in range(i+1): - vals = demu_phh_mass(ks=k_use/h_, - M1=M_array[i]*h_, M2=M_array[j]*h_, - redshift=1./aa-1.) / h_**3. - pkhh[:, i, j] = vals - pkhh[:, j, i] = vals - # pkhh = _compute_p_hh(k_use, M_array, redshift=1./aa-1) + # (Nk,NM,NM) + pkhh = _compute_p_hh(dehod=demuhod, ks=k_use/h_, + Mh=M_array*h_, redshift=1./aa-1) / h_**3. # integration pk_2h_M2_int = list() for i in range(nM): + res = pkhh[:, i] * u1k pk_2h_M2_int.append(hmc._integrator( - pkhh[:, i] * u1k * hmf[None, :], lmass)) - pk_2h = hmc._integrator(np.array(pk_2h_M2_int).T * - u2k * hmf[None, :], lmass) - ######### + hmf[None, :] * res, lmass)) + res = np.array(pk_2h_M2_int).T * u2k + pk_2h = hmc._integrator(hmf[None, :] * res, lmass) else: - pkhm = np.zeros([nk, nM]) - for i in range(nM): - pkhm[:, i] = demu_phm_mass(ks=k_use/h_, M=M_array[i]*h_, - redshift=1./aa-1.) / h_**3. - # pkhm = _compute_p_hm(k_use, M_array, redshift=1./aa-1) # (Nk,NM) + pkhm = _compute_p_hm(demuhod, k_use/h_, M_array*h_, + redshift=1./aa-1) / h_**3. - pk_2h = hmc._integrator(pkhm * u1k * hmf[None, :], lmass) + # integration + res = pkhm * u1k + pk_2h = hmc._integrator(hmf[None, :] * res, lmass) if get_1h and (prof2 == prof): - pk_1h = hmc.I_0_2(cosmo, k_use, aa, prof, - prof2=prof2, prof_2pt=prof_2pt) # 1h term - + pk_1h = _I_0_2(hmc, cosmo, k_use, aa, prof, + prof2=prof2, prof_2pt=prof_2pt, + lmass=lmass, hmf=hmf) # 1h term + if suppress_1h is not None: # large-scale damping of 1-halo term ks = suppress_1h(aa) @@ -241,3 +306,4 @@ def darkemu_power_spectrum(demu, cosmo, hmc, k, a, prof, *, if np.ndim(k) == 0: out = np.squeeze(out, axis=-1) return out + From 9c93ee4e87f1271ad70128d41f4a6eea71793951 Mon Sep 17 00:00:00 2001 From: Keitaro Ishikawa Date: Mon, 4 Aug 2025 18:27:28 +0900 Subject: [PATCH 6/7] fixed flake8 errors --- pyccl/halos/emulators/darkemu1_pk.py | 46 +++++++++++++--------------- 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/pyccl/halos/emulators/darkemu1_pk.py b/pyccl/halos/emulators/darkemu1_pk.py index 0a6f244b4..603b69b1e 100644 --- a/pyccl/halos/emulators/darkemu1_pk.py +++ b/pyccl/halos/emulators/darkemu1_pk.py @@ -15,14 +15,14 @@ def galaxy_bias(cosmo, hmc, a, prof): """ Computes the galaxy bias for a given halo profile.""" - Mh = np.logspace(12.2,15.0,2**5+1) / cosmo['h'] # [Msun] + Mh = np.logspace(12.2, 15.0, 2**5+1) / cosmo['h'] # [Msun] lmass = np.log10(Mh) NgalM = _Ngal(Mh, 1., profile=prof) # hmc._mf: dn/dlogM [Mpc^-3] (comoving) not dn/dM hmc._get_ingredients(cosmo, a, get_bf=False) hmf = ius(hmc._lmass, hmc._mf)(lmass) norm1 = hmc._integrator(hmf * NgalM, lmass) # integral of dn/dM - hbias = np.array([demuhod.get_bias_mass(Mi, redshift=1./a - 1)[0,0] + hbias = np.array([demuhod.get_bias_mass(Mi, redshift=1./a - 1)[0, 0] for Mi in Mh*cosmo['h']]) galbias = hmc._integrator(hmf * hbias * NgalM, lmass) / norm1 # integral of dn/dM @@ -43,11 +43,11 @@ def galaxy_bias_DEmuxHOD(cosmo, a, prof): raise ValueError("Dark Emulator needs A_s as an input parameter.") # [omegab0.,omega_c0.,Omega_Lambda0,np.log(As*10**10),ns,w0] demuhod.set_cosmology(np.array([cosmo['Omega_b']*h_**2., - cosmo['Omega_c']*h_**2., - 1-(cosmo['Omega_c']+cosmo['Omega_b']), - np.log(cosmo['A_s']*10**10), - cosmo['n_s'], - cosmo['w0']])) + cosmo['Omega_c']*h_**2., + 1-(cosmo['Omega_c']+cosmo['Omega_b']), + np.log(cosmo['A_s']*10**10), + cosmo['n_s'], + cosmo['w0']])) # HOD parameters for the dark emulator gparam_input = {"logMmin": prof.log10Mmin_0 + np.log10(h_), "sigma_sq": (prof.siglnM_0/np.log(10))**2, @@ -59,7 +59,7 @@ def galaxy_bias_DEmuxHOD(cosmo, a, prof): "alpha_inc": 0., "logM_inc": 0.} demuhod.set_galaxy(gparam_input) - galbias = np.array([demuhod._get_effective_bias(redshift=z_) + galbias = np.array([demuhod._get_effective_bias(redshift=z_) for z_ in (1./a_use - 1.)]) if len(galbias) == 1: return galbias[0] @@ -80,21 +80,22 @@ def _compute_p_hh(dehod, ks, Mh, redshift): dehod._compute_logdens(redshift) dehod._compute_p_hh_spl(redshift) - points_known = (-dehod.g1.logdens_list, -dehod.g1.logdens_list, dehod.fftlog_2h.k) + points_known = (-dehod.g1.logdens_list, -dehod.g1.logdens_list, + dehod.fftlog_2h.k) logdens = _compute_logdens(dehod=dehod, Mh=Mh, redshift=redshift) grid_d1, grid_d2, grid_k = np.meshgrid(-logdens, -logdens, ks, indexing='ij') points_target = np.stack([grid_d1.ravel(), grid_d2.ravel(), grid_k.ravel()], axis=-1) - interpolator = rgi(points_known, dehod.p_hh_base, method='cubic', + interpolator = rgi(points_known, dehod.p_hh_base, method='cubic', bounds_error=False, fill_value=None) - p_hh = interpolator(points_target).reshape(len(logdens), len(logdens), len(ks)) + p_hh = interpolator(points_target).reshape(len(logdens), + len(logdens), len(ks)) return p_hh.transpose(2, 0, 1) - # referred from dark_emulator_public/dark_emulator/model_hod/hod_interface.py def _compute_p_hm(dehod, ks, Mh, redshift): if not dehod.logdens_computed: @@ -204,11 +205,11 @@ def darkemu_power_spectrum(cosmo, hmc, k, a, prof, *, raise ValueError("Dark Emulator needs A_s as an input parameter.") # [omegab0.,omega_c0.,Omega_Lambda0,np.log(As*10**10),ns,w0] demuhod.set_cosmology(np.array([cosmo['Omega_b']*h_**2., - cosmo['Omega_c']*h_**2., - 1-(cosmo['Omega_c']+cosmo['Omega_b']), - np.log(cosmo['A_s']*10**10), - cosmo['n_s'], - cosmo['w0']])) + cosmo['Omega_c']*h_**2., + 1-(cosmo['Omega_c']+cosmo['Omega_b']), + np.log(cosmo['A_s']*10**10), + cosmo['n_s'], + cosmo['w0']])) # HOD parameters for the dark emulator gparam_input = {"logMmin": prof.log10Mmin_0 + np.log10(h_), "sigma_sq": (prof.siglnM_0/np.log(10))**2, @@ -242,7 +243,6 @@ def darkemu_power_spectrum(cosmo, hmc, k, a, prof, *, lmass = np.log10(M_array) nM = len(M_array) - for ia, aa in enumerate(a_use): # normalizations # norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) @@ -255,15 +255,14 @@ def darkemu_power_spectrum(cosmo, hmc, k, a, prof, *, if prof2 == prof: norm2 = norm1 - # calc profile in fourier space u1k = prof.fourier(cosmo, k_use, M_array, aa).T # (Nk,NM) if prof2 == prof: # pkgg u2k = u1k - + # (Nk,NM,NM) pkhh = _compute_p_hh(dehod=demuhod, ks=k_use/h_, - Mh=M_array*h_, redshift=1./aa-1) / h_**3. + Mh=M_array*h_, redshift=1./aa-1) / h_**3. # integration pk_2h_M2_int = list() @@ -285,9 +284,9 @@ def darkemu_power_spectrum(cosmo, hmc, k, a, prof, *, if get_1h and (prof2 == prof): pk_1h = _I_0_2(hmc, cosmo, k_use, aa, prof, - prof2=prof2, prof_2pt=prof_2pt, + prof2=prof2, prof_2pt=prof_2pt, lmass=lmass, hmf=hmf) # 1h term - + if suppress_1h is not None: # large-scale damping of 1-halo term ks = suppress_1h(aa) @@ -306,4 +305,3 @@ def darkemu_power_spectrum(cosmo, hmc, k, a, prof, *, if np.ndim(k) == 0: out = np.squeeze(out, axis=-1) return out - From a0e26397d226059590ceb7d0adf7948c4aa53d25 Mon Sep 17 00:00:00 2001 From: Keitaro Ishikawa Date: Fri, 28 Nov 2025 16:11:59 +0900 Subject: [PATCH 7/7] add "hybrid" option for Pgg at low k --- pyccl/halos/emulators/darkemu1_pk.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/pyccl/halos/emulators/darkemu1_pk.py b/pyccl/halos/emulators/darkemu1_pk.py index 603b69b1e..84b854eb5 100644 --- a/pyccl/halos/emulators/darkemu1_pk.py +++ b/pyccl/halos/emulators/darkemu1_pk.py @@ -8,6 +8,9 @@ from scipy.interpolate import RegularGridInterpolator as rgi from .. import Profile2pt +from .. import halomod_power_spectrum + +import pyccl as ccl from dark_emulator import model_hod demuhod = model_hod.darkemu_x_hod() @@ -140,7 +143,8 @@ def _I_0_2(hmc, cosmo, k, a, prof, *, prof2=None, prof_2pt, lmass, hmf): def darkemu_power_spectrum(cosmo, hmc, k, a, prof, *, prof2=None, prof_2pt=None, get_1h=True, get_2h=True, - suppress_1h=None): + suppress_1h=None, + hybrid=False): """ Computes the halo model power spectrum for two quantities defined by their respective halo profiles. The halo model power spectrum for two profiles @@ -272,6 +276,14 @@ def darkemu_power_spectrum(cosmo, hmc, k, a, prof, *, hmf[None, :] * res, lmass)) res = np.array(pk_2h_M2_int).T * u2k pk_2h = hmc._integrator(hmf[None, :] * res, lmass) + if hybrid: + # hybrid approach + pk_gg_halomodel = halomod_power_spectrum(cosmo, hmc, k_use, aa, prof, prof_2pt=prof_2pt, get_1h=False) * (norm1 * norm2) + # smooth transition between the two approaches + k_transition = 0.01 # [Mpc^-1] + weight = np.exp(-(k_use / k_transition)**4.) + pk_2h = (1 - weight) * pk_2h + weight * pk_gg_halomodel + else: # (Nk,NM)