Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

No normprof #1072

Merged
merged 9 commits into from
Apr 19, 2023
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 7 additions & 5 deletions benchmarks/test_covariances.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
import os

import numpy as np

import pytest
import pyccl as ccl


Expand Down Expand Up @@ -31,9 +30,12 @@ def test_ssc_WL():
a = np.linspace(1/(1+6), 1, n_z)
k = np.geomspace(k_min, k_max, n_k)

tk3D = ccl.halos.halomod_Tk3D_SSC(cosmo=cosmo, hmc=hmc,
prof=nfw, prof2=nfw, prof12_2pt=None,
lk_arr=np.log(k), a_arr=a, use_log=True)
with pytest.warns():
nikfilippas marked this conversation as resolved.
Show resolved Hide resolved
tk3D = ccl.halos.halomod_Tk3D_SSC(cosmo=cosmo, hmc=hmc,
prof=nfw, prof2=nfw, prof12_2pt=None,
lk_arr=np.log(k), a_arr=a,
use_log=True,
normprof1=True, normprof2=True)
nikfilippas marked this conversation as resolved.
Show resolved Hide resolved

z, nofz = np.loadtxt(os.path.join(data_dir, "ssc_WL_nofz.txt"),
unpack=True)
Expand Down
8 changes: 6 additions & 2 deletions benchmarks/test_halomod.py
Original file line number Diff line number Diff line change
Expand Up @@ -43,15 +43,19 @@ def test_halomod(model):
z = 0.
k = data_z0[:, 0] * cosmo['h']
pk = data_z0[:, -1] / (cosmo['h']**3)
pk_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k, 1./(1+z), prf)
with pytest.warns():
nikfilippas marked this conversation as resolved.
Show resolved Hide resolved
pk_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k, 1./(1+z), prf,
normprof1=True)
nikfilippas marked this conversation as resolved.
Show resolved Hide resolved
tol = pk * HALOMOD_TOLERANCE
err = np.abs(pk_ccl - pk)
assert np.all(err <= tol)

z = 1.
k = data_z1[:, 0] * cosmo['h']
pk = data_z1[:, -1] / (cosmo['h']**3)
pk_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k, 1./(1+z), prf)
with pytest.warns():
nikfilippas marked this conversation as resolved.
Show resolved Hide resolved
pk_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k, 1./(1+z), prf,
normprof1=True)
nikfilippas marked this conversation as resolved.
Show resolved Hide resolved
tol = pk * HALOMOD_TOLERANCE
err = np.abs(pk_ccl - pk)
assert np.all(err <= tol)
3 changes: 2 additions & 1 deletion benchmarks/test_hod.py
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,8 @@ def _nz_2mrs(z):
# P(k)
a_arr, lk_arr, _ = cosmo.get_linear_power().get_spline_arrays()
pk_hod = ccl.halos.halomod_Pk2D(cosmo, hmc, prf, prof_2pt=prf2pt,
lk_arr=lk_arr, a_arr=a_arr)
lk_arr=lk_arr, a_arr=a_arr,
normprof1=True)
nikfilippas marked this conversation as resolved.
Show resolved Hide resolved
# C_ell
tr = ccl.NumberCountsTracer(cosmo, has_rsd=False, dndz=(z_arr, dndz),
bias=(z_arr, np.ones(len(dndz))))
Expand Down
5 changes: 2 additions & 3 deletions pyccl/base/deprecations.py
Original file line number Diff line number Diff line change
Expand Up @@ -146,9 +146,8 @@ def wrapper(*args, **kwargs):
if any([par.startswith("normprof") and kwargs.get(par) is not None
for par in kwargs]):
warnings.warn(
"Argument `normprof` has been deprecated. Change the default "
"value only by subclassing. More comprehensive profile "
"normalization options will be provided with CCLv3.0.0.",
"Argument `normprof` will be deprecated in CCL v3. All "
"profiles will carry their own normalisation.",
nikfilippas marked this conversation as resolved.
Show resolved Hide resolved
CCLDeprecationWarning)

# API compatibility for deprecated HMCalculator argument k_min.
Expand Down
10 changes: 0 additions & 10 deletions pyccl/halos/halo_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -156,16 +156,6 @@ def profile_norm(self, cosmo, a, prof):
self._mass, a, mass_def=self.mass_def).T
return 1. / self._integrate_over_mf(uk0)

def get_profile_norm(self, cosmo, a, prof):
"""Compute the normalization of a profile."""
if not prof.normprof: # TODO: Remove for CCLv3.
return 1
uk0 = prof._normalization(self)(cosmo=cosmo, a=a)
if isinstance(uk0, (int, float)):
return 1 / uk0
self._get_ingredients(cosmo, a, get_bf=False)
return 1 / self._integrate_over_mf(uk0)

@warn_api(pairs=[("sel", "selection"),
("amin", "a_min"),
("amax", "a_max")],
Expand Down
12 changes: 4 additions & 8 deletions pyccl/halos/pk_1pt.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
from ..base import UnlockInstance, warn_api
from ..base import warn_api
import numpy as np


Expand Down Expand Up @@ -28,11 +28,6 @@ def _Ix1(func, cosmo, hmc, k, a, prof, normprof):
`k` and `a` respectively. If `k` or `a` are scalars, the
corresponding dimension will be squeezed out on output.
"""
# TODO: Remove for CCLv3.
if normprof is not None:
with UnlockInstance(prof):
prof.normprof = normprof

func = getattr(hmc, func)

a_use = np.atleast_1d(a).astype(float)
Expand All @@ -43,8 +38,9 @@ def _Ix1(func, cosmo, hmc, k, a, prof, normprof):
out = np.zeros([na, nk])
for ia, aa in enumerate(a_use):
i11 = func(cosmo, k_use, aa, prof)
norm = hmc.get_profile_norm(cosmo, aa, prof)
out[ia] = i11 * norm
norm = prof.normalization(hmc, cosmo, aa) if normprof else 1
# TODO: CCLv3 remove if
out[ia] = i11 / norm
nikfilippas marked this conversation as resolved.
Show resolved Hide resolved

if np.ndim(a) == 0:
out = np.squeeze(out, axis=0)
Expand Down
20 changes: 7 additions & 13 deletions pyccl/halos/pk_2pt.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
from ..base import UnlockInstance, warn_api
from ..base import warn_api
from ..pk2d import Pk2D, parse_pk
from .profiles_2pt import Profile2pt
import numpy as np
Expand Down Expand Up @@ -105,14 +105,6 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof, *,
if prof_2pt is None:
prof_2pt = Profile2pt()

# TODO: Remove for CCLv3.
if normprof1 is not None:
with UnlockInstance(prof):
prof.normprof = normprof1
if normprof2 is not None:
with UnlockInstance(prof2):
prof2.normprof = normprof2

pk2d = parse_pk(cosmo, p_of_k_a)
extrap = cosmo if extrap_pk else None # extrapolation rule for pk2d

Expand All @@ -121,12 +113,14 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof, *,
out = np.zeros([na, nk])
for ia, aa in enumerate(a_use):
# normalizations
norm1 = hmc.get_profile_norm(cosmo, aa, prof)
norm1 = prof.normalization(hmc, cosmo, aa) if normprof1 else 1
# TODO: CCLv3, remove if

if prof2 == prof:
norm2 = norm1
else:
norm2 = hmc.get_profile_norm(cosmo, aa, prof2)
norm2 = prof2.normalization(hmc, cosmo, aa) if normprof2 else 1
# TODO: CCLv3, remove if

if get_2h:
# bias factors
Expand Down Expand Up @@ -154,10 +148,10 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof, *,

# smooth 1h/2h transition region
if smooth_transition is None:
out[ia] = (pk_1h + pk_2h) * norm1 * norm2
out[ia] = (pk_1h + pk_2h) / (norm1 * norm2)
else:
alpha = smooth_transition(aa)
out[ia] = (pk_1h**alpha + pk_2h**alpha)**(1/alpha) * norm1 * norm2
out[ia] = (pk_1h**alpha + pk_2h**alpha)**(1/alpha) / (norm1*norm2)

if np.ndim(a) == 0:
out = np.squeeze(out, axis=0)
Expand Down
68 changes: 30 additions & 38 deletions pyccl/halos/pk_4pt.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
from ..base import UnlockInstance, warn_api
from ..base import warn_api
from ..pk2d import parse_pk
from ..tk3d import Tk3D
from ..errors import CCLWarning
Expand Down Expand Up @@ -82,30 +82,33 @@ def halomod_trispectrum_1h(cosmo, hmc, k, a, prof, *,

# define all the profiles
prof, prof2, prof3, prof4, prof12_2pt, prof34_2pt = \
_allocate_profiles(prof, prof2, prof3, prof4, prof12_2pt, prof34_2pt,
normprof1, normprof2, normprof3, normprof4)
_allocate_profiles(prof, prof2, prof3, prof4, prof12_2pt, prof34_2pt)

na = len(a_use)
nk = len(k_use)
out = np.zeros([na, nk, nk])
for ia, aa in enumerate(a_use):
# normalizations
norm1 = hmc.get_profile_norm(cosmo, aa, prof)
norm1 = prof.normalization(hmc, cosmo, aa) if normprof1 else 1
# TODO: CCLv3 remove if

if prof2 == prof:
norm2 = norm1
else:
norm2 = hmc.get_profile_norm(cosmo, aa, prof2)
norm2 = prof2.normalization(hmc, cosmo, aa) if normprof2 else 1
# TODO: CCLv3 remove if

if prof3 == prof:
norm3 = norm1
else:
norm3 = hmc.get_profile_norm(cosmo, aa, prof3)
norm3 = prof3.normalization(hmc, cosmo, aa) if normprof3 else 1
# TODO: CCLv3 remove if

if prof4 == prof2:
norm4 = norm2
else:
norm4 = hmc.get_profile_norm(cosmo, aa, prof4)
norm4 = prof4.normalization(hmc, cosmo, aa) if normprof4 else 1
# TODO: CCLv3 remove if

# trispectrum
tk_1h = hmc.I_0_22(cosmo, k_use, aa,
Expand All @@ -114,7 +117,7 @@ def halomod_trispectrum_1h(cosmo, hmc, k, a, prof, *,
prof12_2pt=prof12_2pt,
prof34_2pt=prof34_2pt)

out[ia] = tk_1h * norm1 * norm2 * norm3 * norm4 # assign
out[ia] = tk_1h / (norm1 * norm2 * norm3 * norm4) # assign

if np.ndim(a) == 0:
out = np.squeeze(out, axis=0)
Expand Down Expand Up @@ -306,14 +309,15 @@ def halomod_Tk3D_SSC_linear_bias(cosmo, hmc, *, prof,
nk = len(k_use)
dpk12, dpk34 = [np.zeros([na, nk]) for _ in range(2)]
for ia, aa in enumerate(a_arr):
norm = hmc.get_profile_norm(cosmo, aa, prof)**2
norm = prof.normalization(hmc, cosmo, aa)**2
# TODO: CCLv3 remove if
i12 = hmc.I_1_2(cosmo, k_use, aa, prof, prof2=prof, prof_2pt=prof_2pt)

pk = pk2d(k_use, aa, cosmo=extrap)
dpk = pk2d(k_use, aa, derivative=True, cosmo=extrap)

# ~ (47/21 - 1/3 dlogPk/dlogk) * Pk + I12
dpk12[ia] = (47/21 - dpk/3)*pk + i12 * norm
dpk12[ia] = (47/21 - dpk/3)*pk + i12 / norm
dpk34[ia] = dpk12[ia].copy()

# Counter terms for clustering (i.e. - (bA + bB) * PAB)
Expand All @@ -323,7 +327,7 @@ def halomod_Tk3D_SSC_linear_bias(cosmo, hmc, *, prof,
i02 = hmc.I_0_2(cosmo, k_use, aa, prof,
prof2=prof, prof_2pt=prof_2pt)

P_12 = P_34 = pk + i02 * norm
P_12 = P_34 = pk + i02 / norm

if is_number_counts1:
b1 = bias1[ia]
Expand Down Expand Up @@ -439,8 +443,7 @@ def halomod_Tk3D_SSC(

# define all the profiles
prof, prof2, prof3, prof4, prof12_2pt, prof34_2pt = \
_allocate_profiles(prof, prof2, prof3, prof4, prof12_2pt, prof34_2pt,
normprof1, normprof2, normprof3, normprof4)
_allocate_profiles(prof, prof2, prof3, prof4, prof12_2pt, prof34_2pt)

k_use = np.exp(lk_arr)
pk2d = parse_pk(cosmo, p_of_k_a)
Expand All @@ -449,28 +452,32 @@ def halomod_Tk3D_SSC(
dpk12, dpk34 = [np.zeros((len(a_arr), len(k_use))) for _ in range(2)]
for ia, aa in enumerate(a_arr):
# normalizations & I11 integral
norm1 = hmc.get_profile_norm(cosmo, aa, prof)
norm1 = prof.normalization(hmc, cosmo, aa) if normprof1 else 1
# TODO: CCLv3 remove if
i11_1 = hmc.I_1_1(cosmo, k_use, aa, prof)

if prof2 == prof:
norm2 = norm1
i11_2 = i11_1
else:
norm2 = hmc.get_profile_norm(cosmo, aa, prof2)
norm2 = prof2.normalization(hmc, cosmo, aa) if normprof2 else 1
# TODO: CCLv3 remove if
i11_2 = hmc.I_1_1(cosmo, k_use, aa, prof2)

if prof3 == prof:
norm3 = norm1
i11_3 = i11_1
else:
norm3 = hmc.get_profile_norm(cosmo, aa, prof3)
norm3 = prof3.normalization(hmc, cosmo, aa) if normprof3 else 1
# TODO: CCLv3 remove if
i11_3 = hmc.I_1_1(cosmo, k_use, aa, prof3)

if prof4 == prof2:
norm4 = norm2
i11_4 = i11_2
else:
norm4 = hmc.get_profile_norm(cosmo, aa, prof4)
norm4 = prof4.normalization(hmc, cosmo, aa) if normprof4 else 1
# TODO: CCLv3 remove if
i11_4 = hmc.I_1_1(cosmo, k_use, aa, prof4)

# I12 integral
Expand All @@ -487,17 +494,17 @@ def halomod_Tk3D_SSC(
dpk = pk2d(k_use, aa, derivative=True, cosmo=extrap)

# (47/21 - 1/3 dlogPk/dlogk) * I11 * I11 * Pk + I12
dpk12[ia] = norm1 * norm2 * ((47/21 - dpk/3)*i11_1*i11_2*pk + i12_12)
dpk34[ia] = norm3 * norm4 * ((47/21 - dpk/3)*i11_3*i11_4*pk + i12_34)
dpk12[ia] = ((47/21 - dpk/3)*i11_1*i11_2*pk + i12_12) / (norm1 * norm2)
dpk34[ia] = ((47/21 - dpk/3)*i11_3*i11_4*pk + i12_34) / (norm3 * norm4)

# Counter terms for clustering (i.e. - (bA + bB) * PAB)
def _get_counterterm(pA, pB, p2pt, nA, nB, i11_A, i11_B):
"""Helper to compute counter-terms."""
# p : profiles | p2pt : 2-point | n : norms | i11 : I_1_1 integral
bA = i11_A * nA if isinstance(pA, ProfNC) else np.zeros_like(k_use)
bB = i11_B * nB if isinstance(pB, ProfNC) else np.zeros_like(k_use)
bA = i11_A / nA if isinstance(pA, ProfNC) else np.zeros_like(k_use)
bB = i11_B / nB if isinstance(pB, ProfNC) else np.zeros_like(k_use)
i02 = hmc.I_0_2(cosmo, k_use, aa, pA, prof2=pB, prof_2pt=p2pt)
P = nA * nB * (pk * i11_A * i11_B + i02)
P = (pk * i11_A * i11_B + i02) / (nA * nB)
return (bA + bB) * P

if isinstance(prof, ProfNC) or isinstance(prof2, ProfNC):
Expand All @@ -519,8 +526,7 @@ def _get_counterterm(pA, pB, p2pt, nA, nB, i11_A, i11_B):
extrap_order_hik=extrap_order_hik, is_logt=use_log)


def _allocate_profiles(prof, prof2, prof3, prof4, prof12_2pt, prof34_2pt,
normprof1, normprof2, normprof3, normprof4):
def _allocate_profiles(prof, prof2, prof3, prof4, prof12_2pt, prof34_2pt):
"""Helper that controls how the undefined profiles are allocated."""
if prof2 is None:
prof2 = prof
Expand All @@ -533,20 +539,6 @@ def _allocate_profiles(prof, prof2, prof3, prof4, prof12_2pt, prof34_2pt,
if prof34_2pt is None:
prof34_2pt = prof12_2pt

# TODO: Remove for CCLv3.
if normprof1 is not None:
with UnlockInstance(prof):
prof.normprof = normprof1
if normprof2 is not None:
with UnlockInstance(prof2):
prof2.normprof = normprof2
if normprof3 is not None:
with UnlockInstance(prof3):
prof3.normprof = normprof3
if normprof4 is not None:
with UnlockInstance(prof4):
prof4.normprof = normprof4

return prof, prof2, prof3, prof4, prof12_2pt, prof34_2pt


Expand Down
2 changes: 1 addition & 1 deletion pyccl/halos/profiles/cib_shang12.py
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ class HaloProfileCIBShang12(HaloProfileCIB):
"""
__repr_attrs__ = __eq_attrs__ = (
"nu", "alpha", "T0", "beta", "gamma", "s_z", "log10Meff", "siglog10M",
"Mmin", "L0", "concentration", "precision_fftlog", "normprof")
"Mmin", "L0", "concentration", "precision_fftlog",)
__getattr__ = deprecate_attr(pairs=[('l10meff', 'log10Meff'),
('sigLM', 'siglog10M')]
)(super.__getattribute__)
Expand Down
2 changes: 1 addition & 1 deletion pyccl/halos/profiles/einasto.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ class HaloProfileEinasto(HaloProfileMatter):
'cosmo' to calculate the value from cosmology. Default: 'cosmo'
"""
__repr_attrs__ = __eq_attrs__ = (
"truncated", "alpha", "precision_fftlog", "concentration", "normprof",)
"truncated", "alpha", "precision_fftlog", "concentration",)

@warn_api(pairs=[("c_M_relation", "concentration")])
def __init__(self, *, concentration, truncated=True, alpha='cosmo'):
Expand Down
4 changes: 1 addition & 3 deletions pyccl/halos/profiles/gaussian.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,9 +22,7 @@ class HaloProfileGaussian(HaloProfile):
rho0 (:obj:`function`): the amplitude of the profile.
It should have the same signature as `r_scale`.
"""
__repr_attrs__ = __eq_attrs__ = ("r_scale", "rho_0", "precision_fftlog",
"normprof",)
normprof = False
__repr_attrs__ = __eq_attrs__ = ("r_scale", "rho_0", "precision_fftlog",)

@deprecated()
@warn_api
Expand Down
2 changes: 1 addition & 1 deletion pyccl/halos/profiles/hernquist.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ class HaloProfileHernquist(HaloProfileMatter):
"""
__repr_attrs__ = __eq_attrs__ = (
"fourier_analytic", "projected_analytic", "cumul2d_analytic",
"truncated", "concentration", "precision_fftlog", "normprof",)
"truncated", "concentration", "precision_fftlog",)

@warn_api(pairs=[("c_M_relation", "concentration")])
def __init__(self, *, concentration,
Expand Down
Loading