Skip to content

Commit

Permalink
Merge pull request #941 from LSSTDESC/ssc_issue938
Browse files Browse the repository at this point in the history 
Fixing missing terms in Super Sample Covariance
  • Loading branch information
@carlosggarcia
carlosggarcia committed Jun 8, 2022
2 parents 1e1c835 + 0efb7ab commit a1edba2
Show file tree
Hide file tree
Showing 3 changed files with 204 additions and 5 deletions.
56 changes: 56 additions & 0 deletions pyccl/halos/halo_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1093,6 +1093,23 @@ def halomod_Tk3D_SSC(cosmo, hmc,
if (prof34_2pt is not None) and (not isinstance(prof34_2pt, Profile2pt)):
raise TypeError("prof34_2pt must be of type `Profile2pt` or `None`")

# number counts profiles must be normalized
if (prof1.is_number_counts) and (normprof1 is False):
raise ValueError('normprof1 must be True if prof1 is of number ' +
'counts type')
if (prof2 is not None) and (prof2.is_number_counts) and \
(normprof2 is False):
raise ValueError('normprof2 must be True if prof2 is of number ' +
'counts type')
if (prof3 is not None) and (prof3.is_number_counts) and \
(normprof3 is False):
raise ValueError('normprof3 must be True if prof3 is of number ' +
'counts type')
if (prof4 is not None) and (prof4.is_number_counts) and \
(normprof4 is False):
raise ValueError('normprof4 must be True if prof3 is of number ' +
'counts type')

if prof3 is None:
prof3_bak = prof1
else:
Expand Down Expand Up @@ -1163,6 +1180,45 @@ def get_norm(normprof, prof, sf):
dpk12[ia, :] = norm12*((2.2380952381-dpk/3)*i11_1*i11_2*pk+i12_12)
dpk34[ia, :] = norm34*((2.2380952381-dpk/3)*i11_3*i11_4*pk+i12_34)

# Counter terms for clustering (i.e. - (bA + bB) * PAB
if prof1.is_number_counts or (prof2 is None or prof2.is_number_counts):
b1 = b2 = np.zeros_like(k_use)
i02_12 = hmc.I_0_2(cosmo, k_use, aa, prof1, prof12_2pt, prof2)
P_12 = norm12 * (pk * i11_1 * i11_2 + i02_12)

if prof1.is_number_counts:
b1 = i11_1 * norm1

if prof2 is None:
b2 = b1
elif prof2.is_number_counts:
b2 = i11_2 * norm2

dpk12[ia, :] -= (b1 + b2) * P_12

if prof3_bak.is_number_counts or \
((prof3_bak.is_number_counts and prof4 is None) or
(prof4 is not None) and prof4.is_number_counts):
b3 = b4 = np.zeros_like(k_use)
if (prof3 is None) and (prof4 is None) and (prof34_2pt is None):
i02_34 = i02_12
else:
i02_34 = hmc.I_0_2(cosmo, k_use, aa, prof3_bak, prof34_2pt_bak,
prof4)
P_34 = norm34 * (pk * i11_3 * i11_4 + i02_34)

if prof3 is None:
b3 = b1
elif prof3.is_number_counts:
b3 = i11_3 * norm3

if prof4 is None:
b4 = b3
elif prof4.is_number_counts:
b4 = i11_4 * norm4

dpk34[ia, :] -= (b3 + b4) * P_34

if use_log:
if np.any(dpk12 <= 0) or np.any(dpk34 <= 0):
warnings.warn(
Expand Down
2 changes: 2 additions & 0 deletions pyccl/halos/profiles.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ class HaloProfile(object):
calculation.
"""
name = 'default'
is_number_counts = False

def __init__(self):
self.precision_fftlog = {'padding_lo_fftlog': 0.1,
Expand Down Expand Up @@ -1283,6 +1284,7 @@ class HaloProfileHOD(HaloProfile):
satellites when centrals are present.
"""
name = 'HOD'
is_number_counts = True

def __init__(self, c_M_relation,
lMmin_0=12., lMmin_p=0., siglM_0=0.4,
Expand Down
151 changes: 146 additions & 5 deletions pyccl/tests/test_tkkssc.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
import numpy as np
import pytest
import pyccl as ccl
from pyccl.halos.halo_model import halomod_bias_1pt


COSMO = ccl.Cosmology(
Expand All @@ -12,7 +13,10 @@
HBF = ccl.halos.HaloBiasTinker10(COSMO, mass_def=M200)
P1 = ccl.halos.HaloProfileNFW(ccl.halos.ConcentrationDuffy08(M200),
fourier_analytic=True)
# P2 will have is_number_counts = True
P2 = ccl.halos.HaloProfileHOD(ccl.halos.ConcentrationDuffy08(M200))
P2_nogc = ccl.halos.HaloProfileHOD(ccl.halos.ConcentrationDuffy08(M200))
P2_nogc.is_number_counts = False
P3 = ccl.halos.HaloProfilePressureGNFW()
P4 = P1
Pneg = ccl.halos.HaloProfilePressureGNFW(P0=-1)
Expand All @@ -23,6 +27,34 @@
AA = 1.0


def get_ssc_counterterm_gc(k, a, hmc, prof1, prof2, prof12_2pt,
normalize=False):

P_12 = b1 = b2 = np.zeros_like(k)
if prof1.is_number_counts or prof2.is_number_counts:
norm1 = hmc.profile_norm(COSMO, a, prof1)
norm2 = hmc.profile_norm(COSMO, a, prof2)
norm12 = 1
if prof1.is_number_counts:
norm12 *= norm1
if prof2.is_number_counts:
norm12 *= norm2

i11_1 = hmc.I_1_1(COSMO, k, a, prof1)
i11_2 = hmc.I_1_1(COSMO, k, a, prof2)
i02_12 = hmc.I_0_2(COSMO, k, a, prof1, prof12_2pt, prof2)

pk = ccl.linear_matter_power(COSMO, k, a)
P_12 = norm12 * (pk * i11_1 * i11_2 + i02_12)

if prof1.is_number_counts:
b1 = halomod_bias_1pt(COSMO, hmc, k, a, prof1) * norm1
if prof2.is_number_counts:
b2 = halomod_bias_1pt(COSMO, hmc, k, a, prof2) * norm2

return (b1 + b2) * P_12


@pytest.mark.parametrize('pars',
[{'p1': P1, 'p2': None, 'cv12': None,
'p3': None, 'p4': None, 'cv34': None,
Expand All @@ -32,7 +64,7 @@
'norm': True, 'pk': None},
{'p1': P1, 'p2': P2, 'cv12': None,
'p3': None, 'p4': None, 'cv34': None,
'norm': False, 'pk': None},
'norm': True, 'pk': None},
{'p1': P1, 'p2': None, 'cv12': None,
'p3': P3, 'p4': None, 'cv34': None,
'norm': False, 'pk': None},
Expand All @@ -41,10 +73,10 @@
'norm': False, 'pk': None},
{'p1': P1, 'p2': P2, 'cv12': None,
'p3': P3, 'p4': P4, 'cv34': None,
'norm': False, 'pk': None},
'norm': True, 'pk': None},
{'p1': P2, 'p2': P2, 'cv12': PKCH,
'p3': P2, 'p4': P2, 'cv34': None,
'norm': False, 'pk': None},
'norm': True, 'pk': None},
{'p1': P1, 'p2': P2, 'cv12': PKC,
'p3': P3, 'p4': P4, 'cv34': PKC,
'norm': True, 'pk': None},
Expand All @@ -56,7 +88,8 @@
'norm': False, 'pk': 'nonlinear'},
{'p1': P1, 'p2': None, 'cv12': None,
'p3': None, 'p4': None, 'cv34': None,
'norm': False, 'pk': COSMO.get_nonlin_power()}],)
'norm': False, 'pk': COSMO.get_nonlin_power()},
],)
def test_tkkssc_smoke(pars):
hmc = ccl.halos.HMCalculator(COSMO, HMF, HBF, mass_def=M200,
nlog10M=2)
Expand All @@ -75,7 +108,8 @@ def test_tkkssc_smoke(pars):
normprof3=pars['norm'],
normprof4=pars['norm'],
p_of_k_a=pars['pk'],
lk_arr=np.log(k_arr), a_arr=a_arr)
lk_arr=np.log(k_arr), a_arr=a_arr,
)
tk = tkk.eval(0.1, 0.5)
assert np.all(np.isfinite(tk))

Expand Down Expand Up @@ -103,8 +137,115 @@ def test_tkkssc_errors():
ccl.halos.halomod_Tk3D_SSC(COSMO, hmc, P1,
prof34_2pt=P2)

# No normalization for number counts profile
with pytest.raises(ValueError):
ccl.halos.halomod_Tk3D_SSC(COSMO, hmc, P2, normprof1=False)
with pytest.raises(ValueError):
ccl.halos.halomod_Tk3D_SSC(COSMO, hmc, P1, prof2=P2, normprof2=False)
with pytest.raises(ValueError):
ccl.halos.halomod_Tk3D_SSC(COSMO, hmc, P1, prof3=P2, normprof3=False)
with pytest.raises(ValueError):
ccl.halos.halomod_Tk3D_SSC(COSMO, hmc, P1, prof4=P2, normprof4=False)

# Negative profile in logspace
assert_warns(ccl.CCLWarning, ccl.halos.halomod_Tk3D_1h,
COSMO, hmc, P3, prof2=Pneg,
lk_arr=np.log(k_arr), a_arr=a_arr,
use_log=True)


@pytest.mark.parametrize('kwargs', [
# All is_number_counts = False
{'prof1': P1, 'prof2': P1,
'prof3': P1, 'prof4': P1},
#
{'prof1': P2, 'prof2': P1,
'prof3': P1, 'prof4': P1},
#
{'prof1': P2, 'prof2': P2,
'prof3': P1, 'prof4': P1},
#
{'prof1': P2, 'prof2': P2,
'prof3': P2, 'prof4': P1},
#
{'prof1': P1, 'prof2': P1,
'prof3': P2, 'prof4': P2},
#
{'prof1': P2, 'prof2': P1,
'prof3': P2, 'prof4': P2},
#
{'prof1': P2, 'prof2': P2,
'prof3': P1, 'prof4': P2},
#
{'prof1': P2, 'prof2': None,
'prof3': P1, 'prof4': None},
#
{'prof1': P2, 'prof2': None,
'prof3': None, 'prof4': None},
#
{'prof1': P2, 'prof2': P2,
'prof3': None, 'prof4': None},
# As in benchmarks/test_covariances.py
{'prof1': P1, 'prof2': P1,
'prof3': None, 'prof4': None},
# Setting prof34_2pt
{'prof1': P1, 'prof2': P1,
'prof3': None, 'prof4': None,
'prof34_2pt': PKC},
# All is_number_counts = True
{'prof1': P2, 'prof2': P2,
'prof3': P2, 'prof4': P2},
]
)
def test_tkkssc_counterterms_gc(kwargs):
hmc = ccl.halos.HMCalculator(COSMO, HMF, HBF, mass_def=M200,
nlog10M=2)
k_arr = KK
a_arr = np.array([0.3, 0.5, 0.7, 1.0])

# Tk's without clustering terms. Set is_number_counts=False for HOD
# profiles
# Ensure HOD profiles are normalized
kwargs_nogc = kwargs.copy()
keys = list(kwargs.keys())
for k in keys:
v = kwargs[k]
if isinstance(v, ccl.halos.HaloProfileHOD):
kwargs_nogc[k] = P2_nogc
kwargs_nogc['norm' + k] = True
kwargs['norm' + k] = True
tkk_nogc = ccl.halos.halomod_Tk3D_SSC(COSMO, hmc,
lk_arr=np.log(k_arr), a_arr=a_arr,
**kwargs_nogc)
_, _, _, tkk_nogc_arrs = tkk_nogc.get_spline_arrays()
tk_nogc_12, tk_nogc_34 = tkk_nogc_arrs

# Tk's with clustering terms
tkk_gc = ccl.halos.halomod_Tk3D_SSC(COSMO, hmc,
lk_arr=np.log(k_arr), a_arr=a_arr,
**kwargs)
_, _, _, tkk_gc_arrs = tkk_gc.get_spline_arrays()
tk_gc_12, tk_gc_34 = tkk_gc_arrs

# Update the None's to their corresponding values
if kwargs['prof2'] is None:
kwargs['prof2'] = kwargs['prof1']
if kwargs['prof3'] is None:
kwargs['prof3'] = kwargs['prof1']
if kwargs['prof4'] is None:
kwargs['prof4'] = kwargs['prof3']

# Tk's of the clustering terms
tkc12 = []
tkc34 = []
for aa in a_arr:
tkc12.append(get_ssc_counterterm_gc(k_arr, aa, hmc, kwargs['prof1'],
kwargs['prof2'], PKC))
tkc34.append(get_ssc_counterterm_gc(k_arr, aa, hmc, kwargs['prof3'],
kwargs['prof4'], PKC))
tkc12 = np.array(tkc12)
tkc34 = np.array(tkc34)

assert np.abs((tk_nogc_12 - tkc12) / tk_gc_12 - 1).max() < 1e-5
assert np.abs((tk_nogc_34 - tkc34) / tk_gc_34 - 1).max() < 1e-5

0 comments on commit a1edba2

Please sign in to comment.