Skip to content

Commit

Permalink
HaloModelCalculator check cosmo in get_ingredients. (#983)
Browse files Browse the repository at this point in the history 
* first commit

* consistent

* back compat

* remove some lingering pytest warnings
  • Loading branch information
@nikfilippas
nikfilippas committed Mar 6, 2023
1 parent 757c70c commit 2c88c59
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 43 deletions.
78 changes: 40 additions & 38 deletions pyccl/halos/halo_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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)`
Expand All @@ -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)
Expand Down Expand Up @@ -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
)

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down
12 changes: 7 additions & 5 deletions pyccl/tests/test_halofit_highz.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,19 +34,21 @@ 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,
w0=-1.04, wa=-0.1,
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()

0 comments on commit 2c88c59

Please sign in to comment.