From e5695c39de3b5eb8018b6367cdae3c504024a2ca Mon Sep 17 00:00:00 2001 From: David Alonso Date: Thu, 27 Jul 2023 09:31:50 +0100 Subject: [PATCH 1/7] implemented, not tested --- .github/environment.yml | 1 + pyccl/halos/hmfunc/__init__.py | 1 + pyccl/halos/hmfunc/bocquet20.py | 78 +++++++++++++++++++++++++++++++++ pyccl/tests/test_hmfunc.py | 7 +-- 4 files changed, 84 insertions(+), 3 deletions(-) create mode 100644 pyccl/halos/hmfunc/bocquet20.py diff --git a/.github/environment.yml b/.github/environment.yml index eef2065d9..c2f88f712 100644 --- a/.github/environment.yml +++ b/.github/environment.yml @@ -17,3 +17,4 @@ dependencies: - isitgr - velocileptors @ git+https://github.com/sfschen/velocileptors - baccoemu @ git+https://bitbucket.org/rangulo/baccoemu.git@master + - MiraTitanHMFemulator \ No newline at end of file diff --git a/pyccl/halos/hmfunc/__init__.py b/pyccl/halos/hmfunc/__init__.py index 950c34bd0..5276d77cb 100644 --- a/pyccl/halos/hmfunc/__init__.py +++ b/pyccl/halos/hmfunc/__init__.py @@ -8,3 +8,4 @@ from .tinker08 import * from .tinker10 import * from .watson13 import * +from .bocquet20 import * diff --git a/pyccl/halos/hmfunc/bocquet20.py b/pyccl/halos/hmfunc/bocquet20.py new file mode 100644 index 000000000..23204a689 --- /dev/null +++ b/pyccl/halos/hmfunc/bocquet20.py @@ -0,0 +1,78 @@ +__all__ = ("MassFuncBocquet20",) + +import numpy as np + +from . import MassFunc + + +class MassFuncBocquet20(MassFunc): + """Implements the mass function emulator of `Bocquet et al. 2020 + `_. This parametrization is + only valid for '200c' 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 = 'Bocquet20' + + def __init__(self, *, + mass_def="200c", + mass_def_strict=True): + import MiraTitanHMFemulator + self.emu = MiraTitanHMFemulator.Emulator() + super().__init__(mass_def=mass_def, mass_def_strict=mass_def_strict) + + def _check_mass_def_strict(self, mass_def): + return mass_def.name != '200c' + + def __call__(self, cosmo, M, a): + # Set up cosmology + h = cosmo['h'] + om = cosmo['Omega_m']*h**2 + ob = cosmo['Omega_b']*h**2 + sigma8 = cosmo['sigma8'] + if np.isnan(sigma8): + raise ValueError("sigma8 must be provided to use the MiraTitan " + "mass function emulator.") + onu = cosmo.omega_x(1.0, 'neutrinos_massive')*h**2 + hmcosmo = {'Ommh2': om, + 'Ombh2': ob, + 'Omnuh2': onu, + 'n_s': cosmo['n_s'], + 'h': h, + 'sigma_8': sigma8, + 'w_0': cosmo['w0'], + 'w_a': cosmo['wa']} + + # Filter out masses beyond emulator range + M_use = np.atleast_1d(M) + # Add h-inverse + Mh = M_use * h + m_hi = Mh > 1E16 + m_lo = Mh < 1E13 + m_good = ~(m_hi | m_lo) + + mfh = np.zeros_like(Mh) + # Populate low-halo masses through extrapolation if needed + if np.any(m_lo): + # Evaluate slope at low masses + m0 = 10**np.array([13.0, 13.1]) + mfp = self.emu.predict(hmcosmo, 1/a-1, m0, + get_errors=False)[0].flatten() + slope = np.log(mfp[1]/mfp[0])/np.log(m0[1]/m0[0]) + mfh[m_lo] = mfp[0]*(Mh[m_lo]/m0[0])**slope + # Predict in good range of masses + if np.any(m_good): + mfp = self.emu.predict(hmcosmo, 1/a-1, Mh[m_good], + get_errors=False)[0].flatten() + mfh[m_good] = mfp + # For masses above emulator range, n(M) will be set to zero + # Remove h-inverse + mf = mfh*h**3 + + if np.ndim(M) == 0: + return mf[0] + return mf diff --git a/pyccl/tests/test_hmfunc.py b/pyccl/tests/test_hmfunc.py index 26b271a7e..b93e6e700 100644 --- a/pyccl/tests/test_hmfunc.py +++ b/pyccl/tests/test_hmfunc.py @@ -4,7 +4,7 @@ COSMO = ccl.Cosmology( - Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.8, n_s=0.96, + Omega_c=0.27, Omega_b=0.05, h=0.67, sigma8=0.8, n_s=0.96, transfer_function='bbks', matter_power_spectrum='linear') HMFS = [ccl.halos.MassFuncPress74, ccl.halos.MassFuncSheth99, @@ -14,7 +14,8 @@ ccl.halos.MassFuncTinker10, ccl.halos.MassFuncWatson13, ccl.halos.MassFuncDespali16, - ccl.halos.MassFuncBocquet16] + ccl.halos.MassFuncBocquet16, + ccl.halos.MassFuncBocquet20] MS = [1E13, [1E12, 1E15], np.array([1E12, 1E15])] MFOF = ccl.halos.MassDef('fof', 'matter') MVIR = ccl.halos.MassDef('vir', 'critical') @@ -24,7 +25,7 @@ M500c = ccl.halos.MassDef(500, 'critical') M500m = ccl.halos.MassDef(500, 'matter') MDFS = [MVIR, MVIR, MVIR, MVIR, - MFOF, MFOF, MVIR, MFOF, MFOF] + MFOF, MFOF, MVIR, MFOF, MFOF, MFOF] @pytest.mark.parametrize('nM_class', HMFS) From 6dd0a50a414edb7f1a52315a1cc2fc4fe0a3d0b0 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Thu, 27 Jul 2023 15:26:39 +0100 Subject: [PATCH 2/7] docs compiled --- readthedocs/api/pyccl.emulators.cosmicemu_MTIV_pk.rst | 7 ------- readthedocs/api/pyccl.halos.hmfunc.bocquet20.rst | 7 +++++++ readthedocs/api/pyccl.halos.hmfunc.rst | 1 + 3 files changed, 8 insertions(+), 7 deletions(-) delete mode 100644 readthedocs/api/pyccl.emulators.cosmicemu_MTIV_pk.rst create mode 100644 readthedocs/api/pyccl.halos.hmfunc.bocquet20.rst diff --git a/readthedocs/api/pyccl.emulators.cosmicemu_MTIV_pk.rst b/readthedocs/api/pyccl.emulators.cosmicemu_MTIV_pk.rst deleted file mode 100644 index 0a4a7b74a..000000000 --- a/readthedocs/api/pyccl.emulators.cosmicemu_MTIV_pk.rst +++ /dev/null @@ -1,7 +0,0 @@ -pyccl.emulators.cosmicemu\_MTIV\_pk module -========================================== - -.. automodule:: pyccl.emulators.cosmicemu_MTIV_pk - :members: - :undoc-members: - :show-inheritance: diff --git a/readthedocs/api/pyccl.halos.hmfunc.bocquet20.rst b/readthedocs/api/pyccl.halos.hmfunc.bocquet20.rst new file mode 100644 index 000000000..99925cb1b --- /dev/null +++ b/readthedocs/api/pyccl.halos.hmfunc.bocquet20.rst @@ -0,0 +1,7 @@ +pyccl.halos.hmfunc.bocquet20 module +=================================== + +.. automodule:: pyccl.halos.hmfunc.bocquet20 + :members: + :undoc-members: + :show-inheritance: diff --git a/readthedocs/api/pyccl.halos.hmfunc.rst b/readthedocs/api/pyccl.halos.hmfunc.rst index c537199ff..f62249226 100644 --- a/readthedocs/api/pyccl.halos.hmfunc.rst +++ b/readthedocs/api/pyccl.halos.hmfunc.rst @@ -9,6 +9,7 @@ Submodules pyccl.halos.hmfunc.angulo12 pyccl.halos.hmfunc.bocquet16 + pyccl.halos.hmfunc.bocquet20 pyccl.halos.hmfunc.despali16 pyccl.halos.hmfunc.jenkins01 pyccl.halos.hmfunc.press74 From 181ca26307fb3f76558f01aae158fbd1df661aba Mon Sep 17 00:00:00 2001 From: David Alonso Date: Thu, 27 Jul 2023 15:37:39 +0100 Subject: [PATCH 3/7] streamlined tests --- pyccl/halos/hmfunc/bocquet20.py | 2 +- pyccl/tests/test_hmfunc.py | 21 +++++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/pyccl/halos/hmfunc/bocquet20.py b/pyccl/halos/hmfunc/bocquet20.py index 23204a689..c8b140dd4 100644 --- a/pyccl/halos/hmfunc/bocquet20.py +++ b/pyccl/halos/hmfunc/bocquet20.py @@ -21,9 +21,9 @@ class MassFuncBocquet20(MassFunc): def __init__(self, *, mass_def="200c", mass_def_strict=True): + super().__init__(mass_def=mass_def, mass_def_strict=mass_def_strict) import MiraTitanHMFemulator self.emu = MiraTitanHMFemulator.Emulator() - super().__init__(mass_def=mass_def, mass_def_strict=mass_def_strict) def _check_mass_def_strict(self, mass_def): return mass_def.name != '200c' diff --git a/pyccl/tests/test_hmfunc.py b/pyccl/tests/test_hmfunc.py index b93e6e700..901f43a47 100644 --- a/pyccl/tests/test_hmfunc.py +++ b/pyccl/tests/test_hmfunc.py @@ -154,3 +154,24 @@ def test_mass_function_mass_def_strict_always_raises(): mdef = ccl.halos.MassDef(400, "critical") with pytest.raises(ValueError): ccl.halos.MassFuncBocquet16(mass_def=mdef, mass_def_strict=False) + + +def test_nM_bocquet20_raises(): + mf = ccl.halos.MassFuncBocquet20(mass_def='200c') + Ms = np.geomspace(1E12, 1E17, 128) + + # Need A_s + cosmo = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=0.67, + A_s=2E-9, n_s=0.96) + with pytest.raises(ValueError): + mf(cosmo, Ms, 1.0) + + # Cosmo parameters out of bounds + cosmo = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=0.67, + sigma8=0.8, n_s=2.0) + with pytest.raises(ValueError): + mf(cosmo, Ms, 1.0) + + # Redshift out of bounds + with pytest.raises(ValueError): + mf(cosmo, Ms, 0.3) From 33357f01a76eee343e9b769422f3e3cd6b12a3dc Mon Sep 17 00:00:00 2001 From: David Alonso Date: Thu, 27 Jul 2023 15:50:59 +0100 Subject: [PATCH 4/7] quantitative test --- pyccl/halos/hmfunc/bocquet20.py | 5 +++-- pyccl/tests/test_hmfunc.py | 21 +++++++++++++++++---- 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/pyccl/halos/hmfunc/bocquet20.py b/pyccl/halos/hmfunc/bocquet20.py index c8b140dd4..f4be19e86 100644 --- a/pyccl/halos/hmfunc/bocquet20.py +++ b/pyccl/halos/hmfunc/bocquet20.py @@ -70,8 +70,9 @@ def __call__(self, cosmo, M, a): get_errors=False)[0].flatten() mfh[m_good] = mfp # For masses above emulator range, n(M) will be set to zero - # Remove h-inverse - mf = mfh*h**3 + # Remove h-inverse and correct for dn/dln(M) -> dn/dlog10(M) + # log10 = 2.30258509299 + mf = mfh*h**3*2.30258509299 if np.ndim(M) == 0: return mf[0] diff --git a/pyccl/tests/test_hmfunc.py b/pyccl/tests/test_hmfunc.py index 901f43a47..76d9146ab 100644 --- a/pyccl/tests/test_hmfunc.py +++ b/pyccl/tests/test_hmfunc.py @@ -26,6 +26,8 @@ M500m = ccl.halos.MassDef(500, 'matter') MDFS = [MVIR, MVIR, MVIR, MVIR, MFOF, MFOF, MVIR, MFOF, MFOF, MFOF] +# This is kind of slow to initialize, so let's do it only once +MF_emu = ccl.halos.MassFuncBocquet20(mass_def='200c') @pytest.mark.parametrize('nM_class', HMFS) @@ -156,22 +158,33 @@ def test_mass_function_mass_def_strict_always_raises(): ccl.halos.MassFuncBocquet16(mass_def=mdef, mass_def_strict=False) +def test_nM_bocquet20_compare(): + # Check that the values are sensible (they don't depart from other + # parametrisations by more than ~10% + Ms = np.geomspace(1E14, 1E15, 128) + mf1 = MF_emu + mf2 = ccl.halos.MassFuncTinker08(mass_def='200c') + + nM1 = mf1(COSMO, Ms, 1.0) + nM2 = mf2(COSMO, Ms, 1.0) + assert np.allclose(nM1, nM2, atol=0, rtol=0.1) + + def test_nM_bocquet20_raises(): - mf = ccl.halos.MassFuncBocquet20(mass_def='200c') Ms = np.geomspace(1E12, 1E17, 128) # Need A_s cosmo = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=0.67, A_s=2E-9, n_s=0.96) with pytest.raises(ValueError): - mf(cosmo, Ms, 1.0) + MF_emu(cosmo, Ms, 1.0) # Cosmo parameters out of bounds cosmo = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=0.67, sigma8=0.8, n_s=2.0) with pytest.raises(ValueError): - mf(cosmo, Ms, 1.0) + MF_emu(cosmo, Ms, 1.0) # Redshift out of bounds with pytest.raises(ValueError): - mf(cosmo, Ms, 0.3) + MF_emu(cosmo, Ms, 0.3) From 93d80f92835e418dd503950d6628a85327527800 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Thu, 27 Jul 2023 16:33:04 +0100 Subject: [PATCH 5/7] benchmarked --- benchmarks/data/codes/bocquet20mf.py | 21 +++++++++++++++++++++ benchmarks/data/mf_bocquet20mf.py.npz | Bin 0 -> 3742 bytes benchmarks/test_bocquet20_hmf.py | 24 ++++++++++++++++++++++++ 3 files changed, 45 insertions(+) create mode 100644 benchmarks/data/codes/bocquet20mf.py create mode 100644 benchmarks/data/mf_bocquet20mf.py.npz create mode 100644 benchmarks/test_bocquet20_hmf.py diff --git a/benchmarks/data/codes/bocquet20mf.py b/benchmarks/data/codes/bocquet20mf.py new file mode 100644 index 000000000..aa8402a73 --- /dev/null +++ b/benchmarks/data/codes/bocquet20mf.py @@ -0,0 +1,21 @@ +import numpy as np +import MiraTitanHMFemulator + + +emu = MiraTitanHMFemulator.Emulator() + +hmcosmo = {'Ommh2': 0.147, + 'Ombh2': 0.02205, + 'Omnuh2': 0.005, + 'n_s': 0.96, + 'h': 0.7, + 'sigma_8': 0.81, + 'w_0': -1.05, + 'w_a': 0.1} + +Ms = np.geomspace(1E13, 1E15, 32) +zs = np.array([0.0, 1.0]) +mfp = emu.predict(hmcosmo, zs, Ms, get_errors=False)[0] +np.savez("../mf_bocquet20mf.py", m=Ms, z=zs, mf=mfp, + **hmcosmo) + diff --git a/benchmarks/data/mf_bocquet20mf.py.npz b/benchmarks/data/mf_bocquet20mf.py.npz new file mode 100644 index 0000000000000000000000000000000000000000..f108aad5679e99b2ee074f3e7412678c01b3a0e1 GIT binary patch literal 3742 zcmWIWW@Zs#fB;2?Yp#168W#cOQUsC!fysWMz5$Vp3}p<} z>M5zk$wlf`3hFj#7V0_*>S_5!B}IvO@%cq5sUUH;#GK+(pm=dcVnHg9uVHMYqp71% zt3W>BVqkE1I>%sf@9K>7HBL$;Q!Ia1Zg9H4W=(#x z>}DsI1;su#y4##q2|t=#xofA>)?fC;p*QzBUH|_|d1Ampryvtk|8lV-PK#0kwE6vx zJH2$QRqNwC?PRIdY$FnP-f6nlw8Mf%mz-W&XUCqGP8&DH2kgB4*6GuuIf7pTKRW5R zJrWU?_~NwuVb9|?{@dO?X=?VXxnin76-)79@Y*F;q=$I_4vz4OC z!p?dQX9flnwXfMc&H*UV_AWG|YXUG}5J_Nd;N zD*CpofeDz%fS3)G`qE(j0EvJAjHb{Z3dTkhB|xqfx6JtBDjJM29x`j6XZJLawL{WctbXXh^zB>Fh`h27QcQw!8NAKNu| zOn&jC_O4wjLs$Q=u4{IVu5P{}-!9lq(wN@fEq~HZ^r%dnum2&tdCs0+`Pp{aeZM>Z zfy9@McFC7NJc{I5Vdt@a0zgR1@9#_e8lYZUF=7rrg z-W74pE|2U!ZRB&aNxx<1BKRir%=SxmOHcgReT3tboqNXhRL89c?fgDH?)g@^-L5)W zXw8d|b#^P|Kb7T4EVkQq`?J>Ggz0vSmXA(xoanZ@Gx_%Hr>qTjr3v=Dhd&hB>GbCB zT3efJcc{QejjcPxE^0E{js!+myUc{-Lg#BHb_VK(ukH9$>{5^5Z@nin!xmplMDAB4HUv5;3fc>Pc=7CQy9!0}E0IG>H2CMVaz~C90^Ow7Q*$;C*H_-XTnd!NS@fL&C6K?`R zvG5Yj6Vq+QPr*C^N<`)H27}cDXU~AZhyAGKDZ?H2cCB&TF;;U9_Rr^CJ|;_JpoV!1_6)`0JvX(t_8If3ey1UGi(R;8o+H*bS5w(;+)_4@C5n5)TYeLPL s$eP-KZTs?ga6U!Xgqj Date: Mon, 31 Jul 2023 16:57:13 +0100 Subject: [PATCH 6/7] comment typo --- pyccl/tests/test_hmfunc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyccl/tests/test_hmfunc.py b/pyccl/tests/test_hmfunc.py index 76d9146ab..c394f4512 100644 --- a/pyccl/tests/test_hmfunc.py +++ b/pyccl/tests/test_hmfunc.py @@ -173,7 +173,7 @@ def test_nM_bocquet20_compare(): def test_nM_bocquet20_raises(): Ms = np.geomspace(1E12, 1E17, 128) - # Need A_s + # Need sigma8 cosmo = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=0.67, A_s=2E-9, n_s=0.96) with pytest.raises(ValueError): From 982378c6ecb63e1d55ae2c177a0234f4299a2bb5 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Mon, 31 Jul 2023 17:39:51 +0100 Subject: [PATCH 7/7] flake8 changed overnight --- pyccl/_core/parameters/fftlog_params.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyccl/_core/parameters/fftlog_params.py b/pyccl/_core/parameters/fftlog_params.py index a3a32576e..5a0b4e782 100644 --- a/pyccl/_core/parameters/fftlog_params.py +++ b/pyccl/_core/parameters/fftlog_params.py @@ -52,7 +52,7 @@ def __repr__(self): def __eq__(self, other): if self is other: True - if type(self) != type(other): + if type(self) is not type(other): return False return self.to_dict() == other.to_dict()