Power spectrum responses for SSC

.github/environment.yml

@@ -1,7 +1,7 @@
 name: test  # default testing environment name from conda-incubator
 dependencies:
   - gfortran=12.2.0  # for isitgr
-  - python=3.8
+  - python=3.11
   - pip
   - setuptools_scm
   - cmake
@@ -17,9 +17,11 @@ dependencies:
   - pytest-cov
   - pytables  # for baccoemu
   - pyfftw  # for velocileptors
+  - george  # for dark_emu
   - pip:
       #- classy<3
       - isitgr
       - velocileptors @ git+https://github.com/sfschen/velocileptors
       - baccoemu @ git+https://bitbucket.org/rangulo/baccoemu.git@master
       - MiraTitanHMFemulator
+      - dark_emulator

benchmarks/data/SSC-Terasawa24/README.md

@@ -0,0 +1,5 @@
+k_h.npy: k [h/Mpc]
+
+PgX_resp_zYY.npy: response of P_gX [(Mpc/h)^3] (X= m or g) at redshift YY: dP_gX/d\delta_b
+
+PgX_resp_err_zYY.npy: error bars on PgX_resp_zYY

benchmarks/test_pkresponse.py

@@ -0,0 +1,143 @@
+import numpy as np
+import os
+from pyccl.pkresponse import Pmm_resp, darkemu_Pgm_resp, darkemu_Pgg_resp
+import pyccl as ccl
+
+data_directory_path = os.path.expanduser("benchmarks/data/SSC-Terasawa24/")
+
+# 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"}},
+)
+
+
+# Construct the full path for each data file (z=0)
+k_data_path = os.path.join(data_directory_path, "k_h.npy")
+k_data_mm_path = os.path.join(data_directory_path, "k_h_mm.npy")
+Pmm_resp_data_path = os.path.join(data_directory_path, "Pmm_resp_z0.npy")
+Pmm_resp_err_data_path = os.path.join(
+    data_directory_path, "Pmm_resp_err_z0.npy"
+)
+Pgm_resp_data_path = os.path.join(data_directory_path, "Pgm_resp_z0.npy")
+Pgm_resp_err_data_path = os.path.join(
+    data_directory_path, "Pgm_resp_err_z0.npy"
+)
+Pgg_resp_data_path = os.path.join(data_directory_path, "Pgg_resp_z0.npy")
+Pgg_resp_err_data_path = os.path.join(
+    data_directory_path, "Pgg_resp_err_z0.npy"
+)
+
+# Load data
+k_data = np.load(k_data_path)
+k_data_mm = np.load(k_data_mm_path)
+Pmm_resp_data = np.load(Pmm_resp_data_path) / h**3
+Pmm_resp_err_data = np.load(Pmm_resp_err_data_path) / h**3
+Pgm_resp_data = np.load(Pgm_resp_data_path) / h**3
+Pgm_resp_err_data = np.load(Pgm_resp_err_data_path) / h**3
+Pgg_resp_data = np.load(Pgg_resp_data_path) / h**3
+Pgg_resp_err_data = np.load(Pgg_resp_err_data_path) / 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])
+indx = (k_data > 1e-2) & (k_data < 4)
+lk_arr = np.log(k_data[indx] * h)  # Using loaded k_data
+indx_mm = (k_data_mm > 1e-2) & (k_data_mm < 4)
+lk_arr_mm = np.log(k_data_mm[indx_mm] * h)  # Using loaded k_data
+
+# Generate power spectrum responses using pkresponse.py functions
+use_log = False
+
+generated_Pmm_resp = Pmm_resp(
+    cosmo, deltah=0.02, lk_arr=lk_arr_mm, a_arr=a_arr, use_log=use_log
+)
+
+generated_Pgm_resp = darkemu_Pgm_resp(
+    cosmo,
+    prof_hod,
+    deltah=0.02,
+    log10Mh_min=log10Mh_min,
+    log10Mh_max=log10Mh_max,
+    lk_arr=lk_arr,
+    a_arr=a_arr,
+    use_log=use_log,
+)
+
+generated_Pgg_resp = darkemu_Pgg_resp(
+    cosmo,
+    prof_hod,
+    deltalnAs=0.03,
+    log10Mh_min=log10Mh_min,
+    log10Mh_max=log10Mh_max,
+    lk_arr=lk_arr,
+    a_arr=a_arr,
+    use_log=use_log,
+)
+
+
+# Compare the generated responses with simulation data
+def test_pmm_resp():
+    assert np.allclose(
+        Pmm_resp_data[indx_mm],
+        generated_Pmm_resp,
+        atol=6 * Pmm_resp_err_data[indx_mm],
+    )
+
+
+def test_pgm_resp():
+    assert np.allclose(
+        Pgm_resp_data[indx],
+        generated_Pgm_resp,
+        atol=2 * Pgm_resp_err_data[indx],
+    )
+
+
+def test_pgg_resp():
+    assert np.allclose(
+        Pgg_resp_data[indx],
+        generated_Pgg_resp,
+        atol=3 * Pgg_resp_err_data[indx],
+    )

pyccl/__init__.py

@@ -45,3 +45,5 @@
 from . import nl_pt
 
 from .cosmology import *
+
+from .pkresponse import *

pyccl/halos/hmfunc/__init__.py

@@ -9,3 +9,4 @@
 from .tinker10 import *
 from .watson13 import *
 from .bocquet20 import *
+from .darkemulator import *

pyccl/halos/hmfunc/darkemulator.py

@@ -0,0 +1,65 @@
+__all__ = ("MassFuncDarkEmulator",)
+
+import numpy as np
+
+from dark_emulator import darkemu
+from . import MassFunc
+
+
+class MassFuncDarkEmulator(MassFunc):
+    """Implements mass function described in 2019ApJ...884..29P.
+    This parametrization is only valid for '200m' masses.
+
+    Args:
+        mass_def (:class:`~pyccl.halos.massdef.MassDef` or :obj:`str`):
+            a mass definition object, or a name string.
+        mass_def_strict (:obj:`bool`): if ``False``, consistency of the mass
+            definition will be ignored.
+    """
+
+    name = "DarkEmulator"
+
+    def __init__(self, *, mass_def="200m", mass_def_strict=True):
+        super().__init__(mass_def=mass_def, mass_def_strict=mass_def_strict)
+        self.emu = darkemu.de_interface.base_class()
+
+    def _check_mass_def_strict(self, mass_def):
+        if isinstance(mass_def.Delta, str):
+            return True
+        elif int(mass_def.Delta) == 200:
+            if mass_def.rho_type != "matter":
+                return True
+        return False
+
+    def _setup(self):
+        pass
+
+    def _get_fsigma(self, cosmo, sigM, a, lnM):
+        z = 1.0 / a - 1
+
+        Omega_c = cosmo["Omega_c"]
+        Omega_b = cosmo["Omega_b"]
+        h = cosmo["h"]
+        n_s = cosmo["n_s"]
+        A_s = cosmo["A_s"]
+
+        omega_c = Omega_c * h**2
+        omega_b = Omega_b * h**2
+        omega_nu = 0.00064
+        Omega_L = 1 - ((omega_c + omega_b + omega_nu) / h**2)
+
+        # Parameters cparam (numpy array) : Cosmological parameters
+        # (𝜔𝑏, 𝜔𝑐, Ω𝑑𝑒, ln(10^10 𝐴𝑠), 𝑛𝑠, 𝑤)
+        cparam = np.array(
+            [omega_b, omega_c, Omega_L, np.log(10**10 * A_s), n_s, -1.0]
+        )
+        self.emu.set_cosmology(cparam)
+
+        alpha = 10 ** (-((0.75 / (np.log10(200 / 75.0))) ** 1.2))
+
+        pA = self.emu.massfunc.coeff_Anorm_spl(-z)
+        pa = self.emu.massfunc.coeff_a_spl(-z)
+        pb = 2.57 * pa**alpha
+        pc = 1.19
+
+        return pA * ((pb / sigM) ** pa + 1.0) * np.exp(-pc / sigM**2)