Skip to content

Commit

Permalink
Merge pull request #877 from nikfilippas/halomod_correction
Browse files Browse the repository at this point in the history 
Enhancements of the Halo Model
  • Loading branch information
@nikfilippas
nikfilippas committed Jun 6, 2021
2 parents 4a9281e + bfa7242 commit 61a13fe
Show file tree
Hide file tree
Showing 2 changed files with 126 additions and 7 deletions.
63 changes: 56 additions & 7 deletions pyccl/halos/halo_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -463,7 +463,8 @@ def halomod_bias_1pt(cosmo, hmc, k, a, prof, normprof=False):
def halomod_power_spectrum(cosmo, hmc, k, a, prof,
prof_2pt=None, prof2=None, p_of_k_a=None,
normprof1=False, normprof2=False,
get_1h=True, get_2h=True):
get_1h=True, get_2h=True,
smooth_transition=None, supress_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
Expand Down Expand Up @@ -509,6 +510,18 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof,
term in the first equation above) won't be computed.
get_2h (bool): if `False`, the 2-halo term (i.e. the second
term in the first equation above) won't be computed.
smooth_transition (function or None):
Modify the halo model 1-halo/2-halo transition region
via a time-dependent function :math:`\\alpha(a)`,
defined as in HMCODE-2020 (``arXiv:2009.01858``): :math:`P(k,a)=
(P_{1h}^{\\alpha(a)}(k)+P_{2h}^{\\alpha(a)}(k))^{1/\\alpha}`.
If `None` the extra factor is just 1.
supress_1h (function or None):
Supress the 1-halo large scale contribution by a
time- and scale-dependent function :math:`k_*(a)`,
defined as in HMCODE-2020 (``arXiv:2009.01858``):
:math:`\\frac{(k/k_*(a))^4}{1+(k/k_*(a))^4}`.
If `None` the standard 1-halo term is returned with no damping.
Returns:
float or array_like: integral values evaluated at each
Expand All @@ -530,6 +543,21 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof,
elif not isinstance(prof_2pt, Profile2pt):
raise TypeError("prof_2pt must be of type "
"`Profile2pt` or `None`")
if smooth_transition is not None:
if not (get_1h and get_2h):
raise ValueError("transition region can only be modified "
"when both 1-halo and 2-halo terms are queried")
if not hasattr(smooth_transition, "__call__"):
raise TypeError("smooth_transition must be "
"a function of `a` or None")
if supress_1h is not None:
if not get_1h:
raise ValueError("can't supress the 1-halo term "
"when get_1h is False")
if not hasattr(supress_1h, "__call__"):
raise TypeError("supress_1h must be "
"a function of `a` or None")

# Power spectrum
if isinstance(p_of_k_a, Pk2D):
def pkf(sf):
Expand Down Expand Up @@ -579,13 +607,19 @@ def pkf(sf):
pk_2h = 0.

if get_1h:
pk_1h = hmc.I_0_2(cosmo, k_use, aa, prof,
prof_2pt, prof2)
pk_1h = hmc.I_0_2(cosmo, k_use, aa, prof, prof_2pt, prof2)
if supress_1h is not None:
ks = supress_1h(aa)
pk_1h *= (k_use / ks)**4 / (1 + (k_use / ks)**4)
else:
pk_1h = 0.

# Total power spectrum
out[ia, :] = (pk_1h + pk_2h) * norm
# Transition region
if smooth_transition is None:
out[ia, :] = (pk_1h + pk_2h) * norm
else:
alpha = smooth_transition(aa)
out[ia, :] = (pk_1h**alpha + pk_2h**alpha)**(1/alpha) * norm

if np.ndim(a) == 0:
out = np.squeeze(out, axis=0)
Expand All @@ -599,7 +633,8 @@ def halomod_Pk2D(cosmo, hmc, prof,
normprof1=False, normprof2=False,
get_1h=True, get_2h=True,
lk_arr=None, a_arr=None,
extrap_order_lok=1, extrap_order_hik=2):
extrap_order_lok=1, extrap_order_hik=2,
smooth_transition=None, supress_1h=None):
""" Returns a :class:`~pyccl.pk2d.Pk2D` object containing
the halo-model power spectrum for two quantities defined by
their respective halo profiles. See :meth:`halomod_power_spectrum`
Expand Down Expand Up @@ -649,6 +684,18 @@ def halomod_Pk2D(cosmo, hmc, prof,
extrap_order_hik (int): extrapolation order to be used on
k-values above the maximum of the splines. See
:class:`~pyccl.pk2d.Pk2D`.
smooth_transition (function or None):
Modify the halo model 1-halo/2-halo transition region
via a time-dependent function :math:`\\alpha(a)`,
defined as in HMCODE-2020 (``arXiv:2009.01858``): :math:`P(k,a)=
(P_{1h}^{\\alpha(a)}(k)+P_{2h}^{\\alpha(a)}(k))^{1/\\alpha}`.
If `None` the extra factor is just 1.
supress_1h (function or None):
Supress the 1-halo large scale contribution by a
time- and scale-dependent function :math:`k_*(a)`,
defined as in HMCODE-2020 (``arXiv:2009.01858``):
:math:`\\frac{(k/k_*(a))^4}{1+(k/k_*(a))^4}`.
If `None` the standard 1-halo term is returned with no damping.
Returns:
:class:`~pyccl.pk2d.Pk2D`: halo model power spectrum.
Expand All @@ -668,7 +715,9 @@ def halomod_Pk2D(cosmo, hmc, prof,
prof, prof_2pt=prof_2pt,
prof2=prof2, p_of_k_a=p_of_k_a,
normprof1=normprof1, normprof2=normprof2,
get_1h=get_1h, get_2h=get_2h)
get_1h=get_1h, get_2h=get_2h,
smooth_transition=smooth_transition,
supress_1h=supress_1h)

pk2d = Pk2D(a_arr=a_arr, lk_arr=lk_arr, pk_arr=pk_arr,
extrap_order_lok=extrap_order_lok,
Expand Down
70 changes: 70 additions & 0 deletions pyccl/tests/test_pkhm.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,10 @@ def f(k, a):
'h2': True, 'itg': 'spline',
'p2': None},
{'cv': None, 'norm': True,
'pk': 'linear', 'h1': True,
'h2': True, 'itg': 'simpson',
'p2': P2},
{'cv': None, 'norm': False,
'pk': 'linear', 'h1': True,
'h2': True, 'itg': 'simpson',
'p2': P2}])
Expand Down Expand Up @@ -182,6 +186,52 @@ def test_pkhm_pk2d():
assert np.all(np.fabs((pk_arr / pk_arr_2 - 1)).flatten()
< 1E-4)

# 1h/2h transition
def alpha0(a): # no smoothing
return 1.

def alpha1(a):
return 0.7

pk0 = ccl.halos.halomod_power_spectrum(COSMO, hmc, k_arr, a_arr, P1,
normprof1=True, normprof2=True,
smooth_transition=None)
pk1 = ccl.halos.halomod_power_spectrum(COSMO, hmc, k_arr, a_arr, P1,
normprof1=True, normprof2=True,
smooth_transition=alpha0)
assert np.allclose(pk0, pk1, rtol=0)
pk2 = ccl.halos.halomod_power_spectrum(COSMO, hmc, k_arr, a_arr, P1,
normprof1=True, normprof2=True,
smooth_transition=alpha1)
assert np.all(pk2/pk0 > 1.)

# 1-halo damping
def ks0(a): # no damping
return 1e-16

def ks1(a): # fully supressed
return 1e16

def ks2(a): # reasonable
return 0.04

pk0 = ccl.halos.halomod_power_spectrum(COSMO, hmc, k_arr, a_arr, P1,
normprof1=True, normprof2=True,
supress_1h=None, get_2h=False)
pk1 = ccl.halos.halomod_power_spectrum(COSMO, hmc, k_arr, a_arr, P1,
normprof1=True, normprof2=True,
supress_1h=ks0, get_2h=False)
assert np.allclose(pk0, pk1, rtol=0)
pk2 = ccl.halos.halomod_power_spectrum(COSMO, hmc, k_arr, a_arr, P1,
normprof1=True, normprof2=True,
supress_1h=ks1, get_2h=False)
assert np.allclose(pk2, 0, rtol=0)
pk3 = ccl.halos.halomod_power_spectrum(COSMO, hmc, k_arr, a_arr, P1,
normprof1=True, normprof2=True,
supress_1h=ks2, get_2h=False)
fact = (k_arr/0.04)**4 / (1 + (k_arr/0.04)**4)
assert np.allclose(pk3, pk0*fact, rtol=0)


def test_pkhm_errors():
# Wrong integration
Expand Down Expand Up @@ -221,3 +271,23 @@ def test_pkhm_errors():
with pytest.raises(TypeError):
ccl.halos.halomod_power_spectrum(COSMO, hmc, KK, AA, P1,
p_of_k_a=KK)

def func():
pass

# Wrong 1h/2h smoothing
with pytest.raises(TypeError):
ccl.halos.halomod_power_spectrum(COSMO, hmc, KK, AA, P1,
smooth_transition=True)
with pytest.raises(ValueError):
ccl.halos.halomod_power_spectrum(COSMO, hmc, KK, AA, P1,
smooth_transition=func, get_1h=False)

# Wrong 1h damping
with pytest.raises(TypeError):
ccl.halos.halomod_power_spectrum(COSMO, hmc, KK, AA, P1,
supress_1h=True)

with pytest.raises(ValueError):
ccl.halos.halomod_power_spectrum(COSMO, hmc, KK, AA, P1,
supress_1h=func, get_1h=False)

0 comments on commit 61a13fe

Please sign in to comment.