From 2c88c59a814b3eb54fcb05f7f681e35eb9abed35 Mon Sep 17 00:00:00 2001 From: Nick Koukoufilippas Date: Mon, 6 Mar 2023 19:00:15 +0200 Subject: [PATCH] `HaloModelCalculator` check cosmo in `get_ingredients`. (#983) * first commit * consistent * back compat * remove some lingering pytest warnings --- pyccl/halos/halo_model.py | 78 ++++++++++++++++--------------- pyccl/tests/test_halofit_highz.py | 12 +++-- 2 files changed, 47 insertions(+), 43 deletions(-) diff --git a/pyccl/halos/halo_model.py b/pyccl/halos/halo_model.py index bc539b16e..7a88c6b5b 100644 --- a/pyccl/halos/halo_model.py +++ b/pyccl/halos/halo_model.py @@ -12,7 +12,7 @@ from ..pyutils import _spline_integrate from .. import background from ..errors import CCLWarning -from ..parameters import physical_constants +from ..parameters import physical_constants as const import numpy as np @@ -106,46 +106,48 @@ def __init__(self, cosmo, massfunc, hbias, mass_def, else: self._integrator = self._integ_spline - self._a_current_mf = -1 - self._a_current_bf = -1 + # Cache last results for mass function and halo bias. + self._cosmo_mf = self._cosmo_bf = None + self._a_mf = self._a_bf = -1 def _integ_spline(self, fM, lM): # Spline integrator return _spline_integrate(lM, fM, lM[0], lM[-1]) - def _get_ingredients(self, a, cosmo, get_bf): - # Compute mass function and bias (if needed) at a new - # value of the scale factor. - rho0 = None - if a != self._a_current_mf: - rho0 = cosmo.rho_x(1., "matter", is_comoving=True) - self.mf = self._massfunc.get_mass_function(cosmo, self._mass, a, - mdef_other=self._mdef) - self.mf0 = (rho0 - - self._integrator(self.mf * self._mass, - self._lmass)) / self._m0 - self._a_current_mf = a - + def _get_mass_function(self, cosmo, a, rho0): + # Compute the mass function at this cosmo and a. + if a != self._a_mf or cosmo != self._cosmo_mf: + massfunc = self._massfunc.get_mass_function + self._mf = massfunc(cosmo, self._mass, a) + integ = self._integrator(self._mf*self._mass, self._lmass) + self._mf0 = (rho0 - integ) / self._m0 + self._cosmo_mf, self._a_mf = cosmo, a # cache + + def _get_halo_bias(self, cosmo, a, rho0): + # Compute the halo bias at this cosmo and a. + if cosmo != self._cosmo_bf or a != self._a_bf: + hbias = self._hbias.get_halo_bias + self._bf = hbias(cosmo, self._mass, a) + integ = self._integrator(self._mf*self._bf*self._mass, self._lmass) + self._mbf0 = (rho0 - integ) / self._m0 + self._cosmo_bf, self._a_bf = cosmo, a # cache + + def _get_ingredients(self, cosmo, a, get_bf): + """Compute mass function and halo bias at some scale factor.""" + rho0 = const.RHO_CRITICAL * cosmo["Omega_m"] * cosmo["h"]**2 + self._get_mass_function(cosmo, a, rho0) if get_bf: - if a != self._a_current_bf: - if rho0 is None: - rho0 = cosmo.rho_x(1., "matter", is_comoving=True) - self.bf = self._hbias.get_halo_bias(cosmo, self._mass, a, - mdef_other=self._mdef) - self.mbf0 = (rho0 - - self._integrator(self.mf * self.bf * self._mass, - self._lmass)) / self._m0 - self._a_current_bf = a + self._get_halo_bias(cosmo, a, rho0) def _integrate_over_mf(self, array_2): - i1 = self._integrator(self.mf[..., :] * array_2, + i1 = self._integrator(self._mf[..., :] * array_2, self._lmass) - return i1 + self.mf0 * array_2[..., 0] + return i1 + self._mf0 * array_2[..., 0] def _integrate_over_mbf(self, array_2): - i1 = self._integrator((self.mf * self.bf)[..., :] * array_2, + i1 = self._integrator((self._mf * self._bf)[..., :] * array_2, self._lmass) - return i1 + self.mbf0 * array_2[..., 0] + return i1 + self._mbf0 * array_2[..., 0] def profile_norm(self, cosmo, a, prof): """ Returns :math:`I^0_1(k\\rightarrow0,a|u)` @@ -161,7 +163,7 @@ def profile_norm(self, cosmo, a, prof): float or array_like: integral value. """ # Compute mass function - self._get_ingredients(a, cosmo, False) + self._get_ingredients(cosmo, a, False) uk0 = prof.fourier(cosmo, self._prec['k_min'], self._mass, a, mass_def=self._mdef).T norm = 1. / self._integrate_over_mf(uk0) @@ -209,17 +211,17 @@ def number_counts(self, cosmo, sel, na=128, amin=None, amax=1.0): abs_dzda = 1 / a / a dc = background.comoving_angular_distance(cosmo, a) ez = background.h_over_h0(cosmo, a) - dh = physical_constants.CLIGHT_HMPC / cosmo['h'] + dh = const.CLIGHT_HMPC / cosmo['h'] dvdz = dh * dc**2 / ez dvda = dvdz * abs_dzda # now do m intergrals in a loop mint = np.zeros_like(a) for i, _a in enumerate(a): - self._get_ingredients(_a, cosmo, False) + self._get_ingredients(cosmo, _a, False) _selm = np.atleast_2d(sel(self._mass, _a)).T mint[i] = self._integrator( - dvda[i] * self.mf[..., :] * _selm[..., :], + dvda[i] * self._mf[..., :] * _selm[..., :], self._lmass ) @@ -250,7 +252,7 @@ def I_0_1(self, cosmo, k, a, prof): value of `k`. """ # Compute mass function - self._get_ingredients(a, cosmo, False) + self._get_ingredients(cosmo, a, False) uk = prof.fourier(cosmo, k, self._mass, a, mass_def=self._mdef).T i01 = self._integrate_over_mf(uk) @@ -280,7 +282,7 @@ def I_1_1(self, cosmo, k, a, prof): value of `k`. """ # Compute mass function and halo bias - self._get_ingredients(a, cosmo, True) + self._get_ingredients(cosmo, a, True) uk = prof.fourier(cosmo, k, self._mass, a, mass_def=self._mdef).T i11 = self._integrate_over_mbf(uk) @@ -316,7 +318,7 @@ def I_0_2(self, cosmo, k, a, prof1, prof_2pt, prof2=None): value of `k`. """ # Compute mass function - self._get_ingredients(a, cosmo, False) + self._get_ingredients(cosmo, a, False) uk = prof_2pt.fourier_2pt(prof1, cosmo, k, self._mass, a, prof2=prof2, mass_def=self._mdef).T @@ -354,7 +356,7 @@ def I_1_2(self, cosmo, k, a, prof1, prof_2pt, prof2=None): value of `k`. """ # Compute mass function - self._get_ingredients(a, cosmo, True) + self._get_ingredients(cosmo, a, True) uk = prof_2pt.fourier_2pt(prof1, cosmo, k, self._mass, a, prof2=prof2, mass_def=self._mdef).T @@ -407,7 +409,7 @@ def I_0_22(self, cosmo, k, a, if prof34_2pt is None: prof34_2pt = prof12_2pt - self._get_ingredients(a, cosmo, False) + self._get_ingredients(cosmo, a, False) uk12 = prof12_2pt.fourier_2pt(prof1, cosmo, k, self._mass, a, prof2=prof2, mass_def=self._mdef).T uk34 = prof34_2pt.fourier_2pt(prof3, cosmo, k, self._mass, a, diff --git a/pyccl/tests/test_halofit_highz.py b/pyccl/tests/test_halofit_highz.py index 17c285d91..f7a26b4d7 100644 --- a/pyccl/tests/test_halofit_highz.py +++ b/pyccl/tests/test_halofit_highz.py @@ -34,6 +34,11 @@ def test_halofit_highz(cosmo): def test_halofit_background_check(): + ccl.spline_params.A_SPLINE_MIN = 0.4 + ccl.spline_params.A_SPLINE_MINLOG = 0.3 + ccl.spline_params.A_SPLINE_MIN_PK = 0.4 + ccl.spline_params.A_SPLINE_MINLOG_PK = 0.3 + cosmo = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=0.7, n_s=0.97, sigma8=0.8, @@ -41,12 +46,9 @@ def test_halofit_background_check(): matter_power_spectrum="halofit", transfer_function="eisenstein_hu") - cosmo.cosmo.spline_params.A_SPLINE_MIN = 0.4 - cosmo.cosmo.spline_params.A_SPLINE_MINLOG = 0.3 - cosmo.cosmo.spline_params.A_SPLINE_MIN_PK = 0.4 - cosmo.cosmo.spline_params.A_SPLINE_MINLOG_PK = 0.3 - k = np.geomspace(1e-3, 1, 10) with pytest.raises(CCLError): ccl.nonlin_matter_power(cosmo, k, a=0.5) + + ccl.spline_params.reload()