From ab5b4e2e4b2858d129037e2f87a0886f39925803 Mon Sep 17 00:00:00 2001
From: David Alonso <dam.phys@gmail.com>
Date: Wed, 19 Apr 2023 10:27:43 +0100
Subject: [PATCH 1/9] extricated normprof

---
 pyccl/base/deprecations.py            |  5 +-
 pyccl/halos/halo_model.py             | 10 ----
 pyccl/halos/pk_1pt.py                 | 12 ++---
 pyccl/halos/pk_2pt.py                 | 20 +++-----
 pyccl/halos/pk_4pt.py                 | 68 ++++++++++++---------------
 pyccl/halos/profiles/cib_shang12.py   |  2 +-
 pyccl/halos/profiles/einasto.py       |  2 +-
 pyccl/halos/profiles/gaussian.py      |  4 +-
 pyccl/halos/profiles/hernquist.py     |  2 +-
 pyccl/halos/profiles/hod.py           | 16 ++++++-
 pyccl/halos/profiles/nfw.py           |  2 +-
 pyccl/halos/profiles/powerlaw.py      |  4 +-
 pyccl/halos/profiles/pressure_gnfw.py |  2 +-
 pyccl/halos/profiles/profile_base.py  | 54 +++++++++------------
 pyccl/tests/test_halomodel.py         |  3 +-
 pyccl/tests/test_tkkssc.py            |  3 +-
 16 files changed, 92 insertions(+), 117 deletions(-)

diff --git a/pyccl/base/deprecations.py b/pyccl/base/deprecations.py
index 9f6f0d6af..fb4ec7ba9 100644
--- a/pyccl/base/deprecations.py
+++ b/pyccl/base/deprecations.py
@@ -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.",
                 CCLDeprecationWarning)
 
         # API compatibility for deprecated HMCalculator argument k_min.
diff --git a/pyccl/halos/halo_model.py b/pyccl/halos/halo_model.py
index 4a64efe37..1d1f79131 100644
--- a/pyccl/halos/halo_model.py
+++ b/pyccl/halos/halo_model.py
@@ -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")],
diff --git a/pyccl/halos/pk_1pt.py b/pyccl/halos/pk_1pt.py
index b0d528607..2a7a0b3d7 100644
--- a/pyccl/halos/pk_1pt.py
+++ b/pyccl/halos/pk_1pt.py
@@ -1,4 +1,4 @@
-from ..base import UnlockInstance, warn_api
+from ..base import warn_api
 import numpy as np
 
 
@@ -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)
@@ -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
 
     if np.ndim(a) == 0:
         out = np.squeeze(out, axis=0)
diff --git a/pyccl/halos/pk_2pt.py b/pyccl/halos/pk_2pt.py
index 8c1bcee64..e9d70f107 100644
--- a/pyccl/halos/pk_2pt.py
+++ b/pyccl/halos/pk_2pt.py
@@ -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
@@ -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
 
@@ -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
@@ -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)
diff --git a/pyccl/halos/pk_4pt.py b/pyccl/halos/pk_4pt.py
index 0daaaef8b..85296553d 100644
--- a/pyccl/halos/pk_4pt.py
+++ b/pyccl/halos/pk_4pt.py
@@ -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
@@ -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,
@@ -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)
@@ -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)
@@ -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]
@@ -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)
@@ -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
@@ -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):
@@ -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
@@ -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
 
 
diff --git a/pyccl/halos/profiles/cib_shang12.py b/pyccl/halos/profiles/cib_shang12.py
index e2d2d9909..b73ec4350 100644
--- a/pyccl/halos/profiles/cib_shang12.py
+++ b/pyccl/halos/profiles/cib_shang12.py
@@ -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__)
diff --git a/pyccl/halos/profiles/einasto.py b/pyccl/halos/profiles/einasto.py
index 88e723a52..08cc9c8cb 100644
--- a/pyccl/halos/profiles/einasto.py
+++ b/pyccl/halos/profiles/einasto.py
@@ -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'):
diff --git a/pyccl/halos/profiles/gaussian.py b/pyccl/halos/profiles/gaussian.py
index 623708f18..6b050fcb7 100644
--- a/pyccl/halos/profiles/gaussian.py
+++ b/pyccl/halos/profiles/gaussian.py
@@ -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
diff --git a/pyccl/halos/profiles/hernquist.py b/pyccl/halos/profiles/hernquist.py
index 01f2180fe..c46f375dd 100644
--- a/pyccl/halos/profiles/hernquist.py
+++ b/pyccl/halos/profiles/hernquist.py
@@ -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,
diff --git a/pyccl/halos/profiles/hod.py b/pyccl/halos/profiles/hod.py
index 3efa94d0c..d71a618b2 100644
--- a/pyccl/halos/profiles/hod.py
+++ b/pyccl/halos/profiles/hod.py
@@ -113,7 +113,7 @@ class HaloProfileHOD(HaloProfileNumberCounts):
         "log10Mmin_0", "log10Mmin_p", "siglnM_0", "siglnM_p", "log10M0_0",
         "log10M0_p", "log10M1_0", "log10M1_p", "alpha_0", "alpha_p", "fc_0",
         "fc_p", "bg_0", "bg_p", "bmax_0", "bmax_p", "a_pivot",
-        "ns_independent", "concentration", "precision_fftlog", "normprof",)
+        "ns_independent", "concentration", "precision_fftlog",)
     __getattr__ = deprecate_attr(pairs=[
         ("lMmin_0", "log10Mmin_0"), ("lMmin_p", "log10Mmin_p"),
         ("siglM_0", "siglnM_0"), ("siglM_p", "siglnM_p"),
@@ -322,6 +322,20 @@ def _real(self, cosmo, r, M, a, mass_def):
             prof = np.squeeze(prof, axis=0)
         return prof
 
+    def normalization(self, hmc, cosmo, a):
+        """Returns the normalization of this profile, which is the
+        mean galaxy number density.
+        """
+        Nc = self._Nc(hmc._mass, a)
+        Ns = self._Ns(hmc._mass, a)
+        fc = self._fc(a)
+        if self.ns_independent:
+            Ngal = Nc*fc + Ns
+        else:
+            Ngal = Nc*(fc + Ns)
+        hmc._get_ingredients(cosmo, a, get_bf=False)
+        return hmc._integrate_over_mf(Ngal)
+
     def _fourier(self, cosmo, k, M, a, mass_def):
         M_use = np.atleast_1d(M)
         k_use = np.atleast_1d(k)
diff --git a/pyccl/halos/profiles/nfw.py b/pyccl/halos/profiles/nfw.py
index c072a3e98..2bbe37855 100644
--- a/pyccl/halos/profiles/nfw.py
+++ b/pyccl/halos/profiles/nfw.py
@@ -47,7 +47,7 @@ class HaloProfileNFW(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,
diff --git a/pyccl/halos/profiles/powerlaw.py b/pyccl/halos/profiles/powerlaw.py
index 624e8e1c4..83baf1385 100644
--- a/pyccl/halos/profiles/powerlaw.py
+++ b/pyccl/halos/profiles/powerlaw.py
@@ -22,9 +22,7 @@ class HaloProfilePowerLaw(HaloProfile):
             profile. The signature of this function should
             be `f(cosmo, a)`.
     """
-    __repr_attrs__ = __eq_attrs__ = ("r_scale", "tilt", "precision_fftlog",
-                                     "normprof",)
-    normprof = False
+    __repr_attrs__ = __eq_attrs__ = ("r_scale", "tilt", "precision_fftlog",)
 
     @deprecated()
     @warn_api
diff --git a/pyccl/halos/profiles/pressure_gnfw.py b/pyccl/halos/profiles/pressure_gnfw.py
index c15567596..6ade1d277 100644
--- a/pyccl/halos/profiles/pressure_gnfw.py
+++ b/pyccl/halos/profiles/pressure_gnfw.py
@@ -67,7 +67,7 @@ class HaloProfilePressureGNFW(HaloProfilePressure):
     """
     __repr_attrs__ = __eq_attrs__ = (
         "mass_bias", "P0", "c500", "alpha", "alpha_P", "beta", "gamma",
-        "P0_hexp", "qrange", "nq", "x_out", "precision_fftlog", "normprof",)
+        "P0_hexp", "qrange", "nq", "x_out", "precision_fftlog",)
 
     @warn_api
     def __init__(self, *, mass_bias=0.8, P0=6.41,
diff --git a/pyccl/halos/profiles/profile_base.py b/pyccl/halos/profiles/profile_base.py
index dcbf5c848..4cc72b491 100644
--- a/pyccl/halos/profiles/profile_base.py
+++ b/pyccl/halos/profiles/profile_base.py
@@ -1,8 +1,8 @@
 from ...pyutils import resample_array, _fftlog_transform
 from ...base import CCLAutoRepr, unlock_instance, warn_api, deprecate_attr
 from ...parameters import FFTLogParams
+from ...parameters import physical_constants as const
 import numpy as np
-from abc import abstractmethod
 import functools
 from typing import Callable
 
@@ -47,35 +47,25 @@ def __init__(self):
                             "_real or _fourier implementation.")
         self.precision_fftlog = FFTLogParams()
 
-    @property
-    @abstractmethod
-    def normprof(self) -> bool:
-        """Normalize the profile in auto- and cross-correlations by
-        :math:`I^0_1(k\\rightarrow 0, a|u)`
-        (see :meth:`~pyccl.halos.halo_model.HMCalculator.I_0_1`).
-        """
-
-    # TODO: CCLv3 - Rename & allocate _normprof_bool to the subclasses.
+    def normalization(self, hmc, cosmo, a):
+        """Returns the overall normalisation factor of this profile.
 
-    def _normprof_false(self, hmc, **settings):
-        """Option for ``normprof = False``."""
-        return lambda *args, cosmo, a, **kwargs: 1.
-
-    def _normprof_true(self, hmc, k_min=1e-5):
-        """Option for ``normprof = True``."""
-        # TODO: remove the first two lines in CCLv3.
-        k_hmc = hmc.precision["k_min"]
-        k_min = k_hmc if k_hmc != k_min else k_min
-        M, mass_def = hmc._mass, hmc.mass_def
-        return functools.partial(self.fourier, k=k_min, M=M, mass_def=mass_def)
+        Args:
+            hmc (:class:`~pyccl.halos.HMCalculator`): a halo model calculator
+                object.
+            cosmo (:class:`~pyccl.core.Cosmology`): a Cosmology object.
+            a (float): scale factor.
 
-    def _normalization(self, hmc, **settings):
-        """This is the API adapter and it decides which norm to use.
-        It returns a function of ``cosmo`` and ``a``. Optional args & kwargs.
+        Reurns:
+            float: normalisation of this profile.
         """
-        if self.normprof:
-            return self._normprof_true(hmc, **settings)
-        return self._normprof_false(hmc, **settings)
+        uk0 = self.fourier(cosmo=cosmo, k=hmc.precision['k_min'],
+                           M=hmc._mass, a=a, mass_def=hmc.mass_def)
+        hmc._get_ingredients(cosmo, a, get_bf=False)
+        return hmc._integrate_over_mf(uk0)
+        # TODO: CCLv3 replace by the below in v3 (profiles will all have a
+        # default normalisation of 1. Normalisation will always be applied).
+        # return lambda *args, cosmo, a, **kwargs: 1.
 
     @unlock_instance(mutate=True)
     @functools.wraps(FFTLogParams.update_parameters)
@@ -482,19 +472,21 @@ def _projected_fftlog_wrap(self, cosmo, r_t, M, a, mass_def,
 
 class HaloProfileNumberCounts(HaloProfile):
     """Base for number counts halo profiles."""
-    normprof = True
 
 
 class HaloProfileMatter(HaloProfile):
     """Base for matter halo profiles."""
-    normprof = True
+
+    def normalization(self, hmc, cosmo, a):
+        """Returns the normalization of all matter overdensity
+        profiles, which is simply the comoving matter density.
+        """
+        return const.RHO_CRITICAL * cosmo["Omega_m"] * cosmo["h"]**2
 
 
 class HaloProfilePressure(HaloProfile):
     """Base for pressure halo profiles."""
-    normprof = False
 
 
 class HaloProfileCIB(HaloProfile):
     """Base for CIB halo profiles."""
-    normprof = False
diff --git a/pyccl/tests/test_halomodel.py b/pyccl/tests/test_halomodel.py
index d2cb8a05f..e8a0378a7 100644
--- a/pyccl/tests/test_halomodel.py
+++ b/pyccl/tests/test_halomodel.py
@@ -69,7 +69,8 @@ def get_pk_new(mf, c, cosmo, a, k, get_1h, get_2h):
                                  mass_def=mdef)
     return ccl.halos.halomod_power_spectrum(cosmo, hmc, k, a, prf,
                                             get_1h=get_1h,
-                                            get_2h=get_2h)
+                                            get_2h=get_2h,
+                                            normprof1=True)
 
 
 @pytest.mark.parametrize('mf_c', [['shethtormen', 'bhattacharya2011'],
diff --git a/pyccl/tests/test_tkkssc.py b/pyccl/tests/test_tkkssc.py
index accc92dc2..73223ef24 100644
--- a/pyccl/tests/test_tkkssc.py
+++ b/pyccl/tests/test_tkkssc.py
@@ -85,7 +85,8 @@ def test_tkkssc_linear_bias(isNC1, isNC2, isNC3, isNC4):
 
     # Test with the full T(k1,k2,a) for an NFW profile with bias ~1.
     tkk = ccl.halos.halomod_Tk3D_SSC(
-        COSMO, HMC, prof=NFW, lk_arr=np.log(KK), a_arr=AA)
+        COSMO, HMC, prof=NFW, lk_arr=np.log(KK), a_arr=AA,
+        normprof1=True)
     *_, (tkk_12, tkk_34) = tkk.get_spline_arrays()
     assert np.allclose(tkkl_12, tkk_12, atol=0, rtol=5e-3)
     assert np.allclose(tkkl_34, tkk_34, atol=0, rtol=5e-3)

From 2ab427cf5397b02c997a6ddfbd122d99b9fa266d Mon Sep 17 00:00:00 2001
From: David Alonso <dam.phys@gmail.com>
Date: Wed, 19 Apr 2023 10:36:30 +0100
Subject: [PATCH 2/9] fixed benchmarks

---
 benchmarks/test_covariances.py | 12 +++++++-----
 benchmarks/test_halomod.py     |  8 ++++++--
 benchmarks/test_hod.py         |  3 ++-
 3 files changed, 15 insertions(+), 8 deletions(-)

diff --git a/benchmarks/test_covariances.py b/benchmarks/test_covariances.py
index 560f82757..57f08f865 100644
--- a/benchmarks/test_covariances.py
+++ b/benchmarks/test_covariances.py
@@ -1,7 +1,6 @@
 import os
-
 import numpy as np
-
+import pytest
 import pyccl as ccl
 
 
@@ -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():
+        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)
 
     z, nofz = np.loadtxt(os.path.join(data_dir, "ssc_WL_nofz.txt"),
                          unpack=True)
diff --git a/benchmarks/test_halomod.py b/benchmarks/test_halomod.py
index a82b4d022..bd1ef8e9b 100644
--- a/benchmarks/test_halomod.py
+++ b/benchmarks/test_halomod.py
@@ -43,7 +43,9 @@ 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():
+        pk_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k, 1./(1+z), prf,
+                                                  normprof1=True)
     tol = pk * HALOMOD_TOLERANCE
     err = np.abs(pk_ccl - pk)
     assert np.all(err <= tol)
@@ -51,7 +53,9 @@ def test_halomod(model):
     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():
+        pk_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k, 1./(1+z), prf,
+                                                  normprof1=True)
     tol = pk * HALOMOD_TOLERANCE
     err = np.abs(pk_ccl - pk)
     assert np.all(err <= tol)
diff --git a/benchmarks/test_hod.py b/benchmarks/test_hod.py
index cf82290c7..88d312305 100644
--- a/benchmarks/test_hod.py
+++ b/benchmarks/test_hod.py
@@ -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)
     # C_ell
     tr = ccl.NumberCountsTracer(cosmo, has_rsd=False, dndz=(z_arr, dndz),
                                 bias=(z_arr, np.ones(len(dndz))))

From c31009db1829768a9968f4012f8ee3871994d4c1 Mon Sep 17 00:00:00 2001
From: David Alonso <dam.phys@gmail.com>
Date: Wed, 19 Apr 2023 11:53:11 +0100
Subject: [PATCH 3/9] comments

---
 benchmarks/test_covariances.py       |  2 +-
 benchmarks/test_halomod.py           |  4 +-
 benchmarks/test_hod.py               |  8 ++--
 pyccl/base/deprecations.py           |  2 +-
 pyccl/halos/halo_model.py            | 60 ++++++++++++++--------------
 pyccl/halos/pk_1pt.py                |  2 +-
 pyccl/halos/pk_2pt.py                |  4 +-
 pyccl/halos/pk_4pt.py                | 18 ++++-----
 pyccl/halos/profiles/hod.py          | 10 ++---
 pyccl/halos/profiles/profile_base.py | 20 ++++++----
 pyccl/tests/test_halomodel.py        |  9 +++--
 pyccl/tests/test_tkkssc.py           |  7 ++--
 12 files changed, 78 insertions(+), 68 deletions(-)

diff --git a/benchmarks/test_covariances.py b/benchmarks/test_covariances.py
index 57f08f865..4465630d0 100644
--- a/benchmarks/test_covariances.py
+++ b/benchmarks/test_covariances.py
@@ -30,7 +30,7 @@ def test_ssc_WL():
     a = np.linspace(1/(1+6), 1, n_z)
     k = np.geomspace(k_min, k_max, n_k)
 
-    with pytest.warns():
+    with pytest.warns(ccl.CCLDeprecationWarning):
         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,
diff --git a/benchmarks/test_halomod.py b/benchmarks/test_halomod.py
index bd1ef8e9b..9ccad4612 100644
--- a/benchmarks/test_halomod.py
+++ b/benchmarks/test_halomod.py
@@ -43,7 +43,7 @@ def test_halomod(model):
     z = 0.
     k = data_z0[:, 0] * cosmo['h']
     pk = data_z0[:, -1] / (cosmo['h']**3)
-    with pytest.warns():
+    with pytest.warns(ccl.CCLDeprecationWarning):
         pk_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k, 1./(1+z), prf,
                                                   normprof1=True)
     tol = pk * HALOMOD_TOLERANCE
@@ -53,7 +53,7 @@ def test_halomod(model):
     z = 1.
     k = data_z1[:, 0] * cosmo['h']
     pk = data_z1[:, -1] / (cosmo['h']**3)
-    with pytest.warns():
+    with pytest.warns(ccl.CCLDeprecationWarning):
         pk_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k, 1./(1+z), prf,
                                                   normprof1=True)
     tol = pk * HALOMOD_TOLERANCE
diff --git a/benchmarks/test_hod.py b/benchmarks/test_hod.py
index 88d312305..623c52d2d 100644
--- a/benchmarks/test_hod.py
+++ b/benchmarks/test_hod.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pytest
 import pyccl as ccl
 
 
@@ -66,9 +67,10 @@ def _nz_2mrs(z):
     prf2pt = ccl.halos.Profile2ptHOD()
     # 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,
-                                    normprof1=True)
+    with pytest.warns(ccl.CCLDeprecationWarning):
+        pk_hod = ccl.halos.halomod_Pk2D(cosmo, hmc, prf, prof_2pt=prf2pt,
+                                        lk_arr=lk_arr, a_arr=a_arr,
+                                        normprof1=True)
     # C_ell
     tr = ccl.NumberCountsTracer(cosmo, has_rsd=False, dndz=(z_arr, dndz),
                                 bias=(z_arr, np.ones(len(dndz))))
diff --git a/pyccl/base/deprecations.py b/pyccl/base/deprecations.py
index fb4ec7ba9..e307e94b2 100644
--- a/pyccl/base/deprecations.py
+++ b/pyccl/base/deprecations.py
@@ -147,7 +147,7 @@ def wrapper(*args, **kwargs):
                 for par in kwargs]):
             warnings.warn(
                 "Argument `normprof` will be deprecated in CCL v3. All "
-                "profiles will carry their own normalisation.",
+                "profiles will carry their own normalization.",
                 CCLDeprecationWarning)
 
         # API compatibility for deprecated HMCalculator argument k_min.
diff --git a/pyccl/halos/halo_model.py b/pyccl/halos/halo_model.py
index 1d1f79131..c0bf372d4 100644
--- a/pyccl/halos/halo_model.py
+++ b/pyccl/halos/halo_model.py
@@ -82,8 +82,8 @@ def __init__(self, *, mass_function, halo_bias, mass_def=None,
             'log10M_min': log10M_min, 'log10M_max': log10M_max, 'nM': nM,
             'integration_method_M': integration_method_M, 'k_min': k_min}
         self._lmass = np.linspace(log10M_min, log10M_max, nM)
-        self._mass = 10.**self._lmass
-        self._m0 = self._mass[0]
+        self.M = 10.**self._lmass
+        self._m0 = self.M[0]
 
         if integration_method_M == "simpson":
             from scipy.integrate import simpson
@@ -106,8 +106,8 @@ def _integ_spline(self, fM, log10M):
     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:
-            self._mf = self.mass_function(cosmo, self._mass, a)
-            integ = self._integrator(self._mf*self._mass, self._lmass)
+            self._mf = self.mass_function(cosmo, self.M, a)
+            integ = self._integrator(self._mf*self.M, self._lmass)
             self._mf0 = (rho0 - integ) / self._m0
             self._cosmo_mf, self._a_mf = cosmo, a  # cache
 
@@ -115,24 +115,24 @@ def _get_mass_function(self, cosmo, a, rho0):
     def _get_halo_bias(self, cosmo, a, rho0):
         # Compute the halo bias at this cosmo and a.
         if a != self._a_bf or cosmo != self._cosmo_bf:
-            self._bf = self.halo_bias(cosmo, self._mass, a)
-            integ = self._integrator(self._mf*self._bf*self._mass, self._lmass)
+            self._bf = self.halo_bias(cosmo, self.M, a)
+            integ = self._integrator(self._mf*self._bf*self.M, 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):
+    def update_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:
             self._get_halo_bias(cosmo, a, rho0)
 
-    def _integrate_over_mf(self, array_2):
+    def integrate_massfunc(self, array_2):
         #  ∫ dM n(M) f(M)
         i1 = self._integrator(self._mf * array_2, self._lmass)
         return i1 + self._mf0 * array_2[..., 0]
 
-    def _integrate_over_mbf(self, array_2):
+    def integrate_massfunc_and_bias(self, array_2):
         #  ∫ dM n(M) b(M) f(M)
         i1 = self._integrator(self._mf * self._bf * array_2, self._lmass)
         return i1 + self._mbf0 * array_2[..., 0]
@@ -151,10 +151,10 @@ def profile_norm(self, cosmo, a, prof):
         Returns:
             float or array_like: integral value.
         """
-        self._get_ingredients(cosmo, a, get_bf=False)
+        self.update_ingredients(cosmo, a, get_bf=False)
         uk0 = prof.fourier(cosmo, self.precision['k_min'],
-                           self._mass, a, mass_def=self.mass_def).T
-        return 1. / self._integrate_over_mf(uk0)
+                           self.M, a, mass_def=self.mass_def).T
+        return 1. / self.integrate_massfunc(uk0)
 
     @warn_api(pairs=[("sel", "selection"),
                      ("amin", "a_min"),
@@ -210,8 +210,8 @@ def number_counts(self, cosmo, *, selection,
         # now do m intergrals in a loop
         mint = np.zeros_like(a)
         for i, _a in enumerate(a):
-            self._get_ingredients(cosmo, _a, get_bf=False)
-            _selm = np.atleast_2d(selection(self._mass, _a)).T
+            self.update_ingredients(cosmo, _a, get_bf=False)
+            _selm = np.atleast_2d(selection(self.M, _a)).T
             mint[i] = self._integrator(
                 dvda[i] * self._mf[..., :] * _selm[..., :],
                 self._lmass
@@ -241,9 +241,9 @@ def I_0_1(self, cosmo, k, a, prof):
             float or array_like: integral values evaluated at each
             value of `k`.
         """
-        self._get_ingredients(cosmo, a, get_bf=False)
-        uk = prof.fourier(cosmo, k, self._mass, a, mass_def=self.mass_def).T
-        return self._integrate_over_mf(uk)
+        self.update_ingredients(cosmo, a, get_bf=False)
+        uk = prof.fourier(cosmo, k, self.M, a, mass_def=self.mass_def).T
+        return self.integrate_massfunc(uk)
 
     def I_1_1(self, cosmo, k, a, prof):
         """ Solves the integral:
@@ -268,9 +268,9 @@ def I_1_1(self, cosmo, k, a, prof):
             float or array_like: integral values evaluated at each
             value of `k`.
         """
-        self._get_ingredients(cosmo, a, get_bf=True)
-        uk = prof.fourier(cosmo, k, self._mass, a, mass_def=self.mass_def).T
-        return self._integrate_over_mbf(uk)
+        self.update_ingredients(cosmo, a, get_bf=True)
+        uk = prof.fourier(cosmo, k, self.M, a, mass_def=self.mass_def).T
+        return self.integrate_massfunc_and_bias(uk)
 
     @warn_api(pairs=[("prof1", "prof")], reorder=["prof_2pt", "prof2"])
     def I_0_2(self, cosmo, k, a, prof, *, prof2=None, prof_2pt):
@@ -307,10 +307,10 @@ def I_0_2(self, cosmo, k, a, prof, *, prof2=None, prof_2pt):
         if prof2 is None:
             prof2 = prof
 
-        self._get_ingredients(cosmo, a, get_bf=False)
-        uk = prof_2pt.fourier_2pt(cosmo, k, self._mass, a, prof,
+        self.update_ingredients(cosmo, a, get_bf=False)
+        uk = prof_2pt.fourier_2pt(cosmo, k, self.M, a, prof,
                                   prof2=prof2, mass_def=self.mass_def).T
-        return self._integrate_over_mf(uk)
+        return self.integrate_massfunc(uk)
 
     @warn_api(pairs=[("prof1", "prof")], reorder=["prof_2pt", "prof2"])
     def I_1_2(self, cosmo, k, a, prof, *, prof2=None, prof_2pt):
@@ -346,10 +346,10 @@ def I_1_2(self, cosmo, k, a, prof, *, prof2=None, prof_2pt):
         if prof2 is None:
             prof2 = prof
 
-        self._get_ingredients(cosmo, a, get_bf=True)
-        uk = prof_2pt.fourier_2pt(cosmo, k, self._mass, a, prof,
+        self.update_ingredients(cosmo, a, get_bf=True)
+        uk = prof_2pt.fourier_2pt(cosmo, k, self.M, a, prof,
                                   prof2=prof2, mass_def=self.mass_def).T
-        return self._integrate_over_mbf(uk)
+        return self.integrate_massfunc_and_bias(uk)
 
     @warn_api(pairs=[("prof1", "prof")],
               reorder=["prof12_2pt", "prof2", "prof3", "prof34_2pt", "prof4"])
@@ -402,9 +402,9 @@ def I_0_22(self, cosmo, k, a, prof, *,
         if prof34_2pt is None:
             prof34_2pt = prof12_2pt
 
-        self._get_ingredients(cosmo, a, get_bf=False)
+        self.update_ingredients(cosmo, a, get_bf=False)
         uk12 = prof12_2pt.fourier_2pt(
-            cosmo, k, self._mass, a, prof,
+            cosmo, k, self.M, a, prof,
             prof2=prof2, mass_def=self.mass_def).T
 
         if (prof, prof2, prof12_2pt) == (prof3, prof4, prof34_2pt):
@@ -412,7 +412,7 @@ def I_0_22(self, cosmo, k, a, prof, *,
             uk34 = uk12
         else:
             uk34 = prof34_2pt.fourier_2pt(
-                cosmo, k, self._mass, a, prof3,
+                cosmo, k, self.M, a, prof3,
                 prof2=prof4, mass_def=self.mass_def).T
 
-        return self._integrate_over_mf(uk12[None, :, :] * uk34[:, None, :])
+        return self.integrate_massfunc(uk12[None, :, :] * uk34[:, None, :])
diff --git a/pyccl/halos/pk_1pt.py b/pyccl/halos/pk_1pt.py
index 2a7a0b3d7..a432e5c5d 100644
--- a/pyccl/halos/pk_1pt.py
+++ b/pyccl/halos/pk_1pt.py
@@ -38,7 +38,7 @@ 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 = prof.normalization(hmc, cosmo, aa) if normprof else 1
+        norm = prof.get_normalization(cosmo, aa, hmc) if normprof else 1
         # TODO: CCLv3 remove if
         out[ia] = i11 / norm
 
diff --git a/pyccl/halos/pk_2pt.py b/pyccl/halos/pk_2pt.py
index e9d70f107..caa3335cc 100644
--- a/pyccl/halos/pk_2pt.py
+++ b/pyccl/halos/pk_2pt.py
@@ -113,13 +113,13 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof, *,
     out = np.zeros([na, nk])
     for ia, aa in enumerate(a_use):
         # normalizations
-        norm1 = prof.normalization(hmc, cosmo, aa) if normprof1 else 1
+        norm1 = prof.get_normalization(cosmo, aa, hmc) if normprof1 else 1
         # TODO: CCLv3, remove if
 
         if prof2 == prof:
             norm2 = norm1
         else:
-            norm2 = prof2.normalization(hmc, cosmo, aa) if normprof2 else 1
+            norm2 = prof2.get_normalization(cosmo, aa, hmc) if normprof2 else 1
             # TODO: CCLv3, remove if
 
         if get_2h:
diff --git a/pyccl/halos/pk_4pt.py b/pyccl/halos/pk_4pt.py
index 85296553d..753e047da 100644
--- a/pyccl/halos/pk_4pt.py
+++ b/pyccl/halos/pk_4pt.py
@@ -89,25 +89,25 @@ def halomod_trispectrum_1h(cosmo, hmc, k, a, prof, *,
     out = np.zeros([na, nk, nk])
     for ia, aa in enumerate(a_use):
         # normalizations
-        norm1 = prof.normalization(hmc, cosmo, aa) if normprof1 else 1
+        norm1 = prof.get_normalization(cosmo, aa, hmc) if normprof1 else 1
         # TODO: CCLv3 remove if
 
         if prof2 == prof:
             norm2 = norm1
         else:
-            norm2 = prof2.normalization(hmc, cosmo, aa) if normprof2 else 1
+            norm2 = prof2.get_normalization(cosmo, aa, hmc) if normprof2 else 1
             # TODO: CCLv3 remove if
 
         if prof3 == prof:
             norm3 = norm1
         else:
-            norm3 = prof3.normalization(hmc, cosmo, aa) if normprof3 else 1
+            norm3 = prof3.get_normalization(cosmo, aa, hmc) if normprof3 else 1
             # TODO: CCLv3 remove if
 
         if prof4 == prof2:
             norm4 = norm2
         else:
-            norm4 = prof4.normalization(hmc, cosmo, aa) if normprof4 else 1
+            norm4 = prof4.get_normalization(cosmo, aa, hmc) if normprof4 else 1
             # TODO: CCLv3 remove if
 
         # trispectrum
@@ -309,7 +309,7 @@ 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 = prof.normalization(hmc, cosmo, aa)**2
+        norm = prof.get_normalization(cosmo, aa, hmc)**2
         # TODO: CCLv3 remove if
         i12 = hmc.I_1_2(cosmo, k_use, aa, prof, prof2=prof, prof_2pt=prof_2pt)
 
@@ -452,7 +452,7 @@ 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 = prof.normalization(hmc, cosmo, aa) if normprof1 else 1
+        norm1 = prof.get_normalization(cosmo, aa, hmc) if normprof1 else 1
         # TODO: CCLv3 remove if
         i11_1 = hmc.I_1_1(cosmo, k_use, aa, prof)
 
@@ -460,7 +460,7 @@ def halomod_Tk3D_SSC(
             norm2 = norm1
             i11_2 = i11_1
         else:
-            norm2 = prof2.normalization(hmc, cosmo, aa) if normprof2 else 1
+            norm2 = prof2.get_normalization(cosmo, aa, hmc) if normprof2 else 1
             # TODO: CCLv3 remove if
             i11_2 = hmc.I_1_1(cosmo, k_use, aa, prof2)
 
@@ -468,7 +468,7 @@ def halomod_Tk3D_SSC(
             norm3 = norm1
             i11_3 = i11_1
         else:
-            norm3 = prof3.normalization(hmc, cosmo, aa) if normprof3 else 1
+            norm3 = prof3.get_normalization(cosmo, aa, hmc) if normprof3 else 1
             # TODO: CCLv3 remove if
             i11_3 = hmc.I_1_1(cosmo, k_use, aa, prof3)
 
@@ -476,7 +476,7 @@ def halomod_Tk3D_SSC(
             norm4 = norm2
             i11_4 = i11_2
         else:
-            norm4 = prof4.normalization(hmc, cosmo, aa) if normprof4 else 1
+            norm4 = prof4.get_normalization(cosmo, aa, hmc) if normprof4 else 1
             # TODO: CCLv3 remove if
             i11_4 = hmc.I_1_1(cosmo, k_use, aa, prof4)
 
diff --git a/pyccl/halos/profiles/hod.py b/pyccl/halos/profiles/hod.py
index d71a618b2..dfbd79666 100644
--- a/pyccl/halos/profiles/hod.py
+++ b/pyccl/halos/profiles/hod.py
@@ -322,19 +322,19 @@ def _real(self, cosmo, r, M, a, mass_def):
             prof = np.squeeze(prof, axis=0)
         return prof
 
-    def normalization(self, hmc, cosmo, a):
+    def get_normalization(self, cosmo, a, hmc):
         """Returns the normalization of this profile, which is the
         mean galaxy number density.
         """
-        Nc = self._Nc(hmc._mass, a)
-        Ns = self._Ns(hmc._mass, a)
+        Nc = self._Nc(hmc.M, a)
+        Ns = self._Ns(hmc.M, a)
         fc = self._fc(a)
         if self.ns_independent:
             Ngal = Nc*fc + Ns
         else:
             Ngal = Nc*(fc + Ns)
-        hmc._get_ingredients(cosmo, a, get_bf=False)
-        return hmc._integrate_over_mf(Ngal)
+        hmc.update_ingredients(cosmo, a, get_bf=False)
+        return hmc.integrate_massfunc(Ngal)
 
     def _fourier(self, cosmo, k, M, a, mass_def):
         M_use = np.atleast_1d(M)
diff --git a/pyccl/halos/profiles/profile_base.py b/pyccl/halos/profiles/profile_base.py
index 4cc72b491..363f47121 100644
--- a/pyccl/halos/profiles/profile_base.py
+++ b/pyccl/halos/profiles/profile_base.py
@@ -47,8 +47,14 @@ def __init__(self):
                             "_real or _fourier implementation.")
         self.precision_fftlog = FFTLogParams()
 
-    def normalization(self, hmc, cosmo, a):
-        """Returns the overall normalisation factor of this profile.
+    def get_normalization(self, cosmo, a, hmc):
+        """Profiles may be normalised by an overall function of redshift
+        (or scale factor). This function may be cosmology dependent and
+        often comes from integrating certain halo properties over mass.
+        This method returns this normalising factor. For example,
+        to get the normalised profile in real space, one would call
+        the `real` method, and then **divide** the result by the value
+        returned by this method.
 
         Args:
             hmc (:class:`~pyccl.halos.HMCalculator`): a halo model calculator
@@ -57,12 +63,12 @@ def normalization(self, hmc, cosmo, a):
             a (float): scale factor.
 
         Reurns:
-            float: normalisation of this profile.
+            float: normalisation factor of this profile.
         """
         uk0 = self.fourier(cosmo=cosmo, k=hmc.precision['k_min'],
-                           M=hmc._mass, a=a, mass_def=hmc.mass_def)
-        hmc._get_ingredients(cosmo, a, get_bf=False)
-        return hmc._integrate_over_mf(uk0)
+                           M=hmc.M, a=a, mass_def=hmc.mass_def)
+        hmc.update_ingredients(cosmo, a, get_bf=False)
+        return hmc.integrate_massfunc(uk0)
         # TODO: CCLv3 replace by the below in v3 (profiles will all have a
         # default normalisation of 1. Normalisation will always be applied).
         # return lambda *args, cosmo, a, **kwargs: 1.
@@ -477,7 +483,7 @@ class HaloProfileNumberCounts(HaloProfile):
 class HaloProfileMatter(HaloProfile):
     """Base for matter halo profiles."""
 
-    def normalization(self, hmc, cosmo, a):
+    def get_normalization(self, cosmo, a, hmc):
         """Returns the normalization of all matter overdensity
         profiles, which is simply the comoving matter density.
         """
diff --git a/pyccl/tests/test_halomodel.py b/pyccl/tests/test_halomodel.py
index e8a0378a7..39d42c26f 100644
--- a/pyccl/tests/test_halomodel.py
+++ b/pyccl/tests/test_halomodel.py
@@ -67,10 +67,11 @@ def get_pk_new(mf, c, cosmo, a, k, get_1h, get_2h):
     prf = ccl.halos.HaloProfileNFW(concentration=cc)
     hmc = ccl.halos.HMCalculator(mass_function=hmf, halo_bias=hbf,
                                  mass_def=mdef)
-    return ccl.halos.halomod_power_spectrum(cosmo, hmc, k, a, prf,
-                                            get_1h=get_1h,
-                                            get_2h=get_2h,
-                                            normprof1=True)
+    with pytest.warns(ccl.CCLDeprecationWarning):
+        return ccl.halos.halomod_power_spectrum(cosmo, hmc, k, a, prf,
+                                                get_1h=get_1h,
+                                                get_2h=get_2h,
+                                                normprof1=True)
 
 
 @pytest.mark.parametrize('mf_c', [['shethtormen', 'bhattacharya2011'],
diff --git a/pyccl/tests/test_tkkssc.py b/pyccl/tests/test_tkkssc.py
index 73223ef24..05b502d3f 100644
--- a/pyccl/tests/test_tkkssc.py
+++ b/pyccl/tests/test_tkkssc.py
@@ -84,9 +84,10 @@ def test_tkkssc_linear_bias(isNC1, isNC2, isNC3, isNC4):
     assert np.allclose(tkkl_12, tkkl_34, atol=0, rtol=1e-12)
 
     # Test with the full T(k1,k2,a) for an NFW profile with bias ~1.
-    tkk = ccl.halos.halomod_Tk3D_SSC(
-        COSMO, HMC, prof=NFW, lk_arr=np.log(KK), a_arr=AA,
-        normprof1=True)
+    with pytest.warns(ccl.CCLDeprecationWarning):
+        tkk = ccl.halos.halomod_Tk3D_SSC(
+            COSMO, HMC, prof=NFW, lk_arr=np.log(KK), a_arr=AA,
+            normprof1=True)
     *_, (tkk_12, tkk_34) = tkk.get_spline_arrays()
     assert np.allclose(tkkl_12, tkk_12, atol=0, rtol=5e-3)
     assert np.allclose(tkkl_34, tkk_34, atol=0, rtol=5e-3)

From 96d3af320e217fe713652de99715468202adae6b Mon Sep 17 00:00:00 2001
From: David Alonso <dam.phys@gmail.com>
Date: Wed, 19 Apr 2023 12:00:15 +0100
Subject: [PATCH 4/9] comments

---
 benchmarks/test_covariances.py | 2 +-
 benchmarks/test_halomod.py     | 4 ++--
 benchmarks/test_hod.py         | 2 +-
 pyccl/tests/test_halomodel.py  | 2 +-
 pyccl/tests/test_tkkssc.py     | 2 +-
 5 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/benchmarks/test_covariances.py b/benchmarks/test_covariances.py
index 4465630d0..f40070a88 100644
--- a/benchmarks/test_covariances.py
+++ b/benchmarks/test_covariances.py
@@ -30,7 +30,7 @@ def test_ssc_WL():
     a = np.linspace(1/(1+6), 1, n_z)
     k = np.geomspace(k_min, k_max, n_k)
 
-    with pytest.warns(ccl.CCLDeprecationWarning):
+    with pytest.warns(ccl.CCLDeprecationWarning):  # TODO: remove normprof v3
         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,
diff --git a/benchmarks/test_halomod.py b/benchmarks/test_halomod.py
index 9ccad4612..26c025c0f 100644
--- a/benchmarks/test_halomod.py
+++ b/benchmarks/test_halomod.py
@@ -43,7 +43,7 @@ def test_halomod(model):
     z = 0.
     k = data_z0[:, 0] * cosmo['h']
     pk = data_z0[:, -1] / (cosmo['h']**3)
-    with pytest.warns(ccl.CCLDeprecationWarning):
+    with pytest.warns(ccl.CCLDeprecationWarning):  # TODO: remove normprof v3
         pk_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k, 1./(1+z), prf,
                                                   normprof1=True)
     tol = pk * HALOMOD_TOLERANCE
@@ -53,7 +53,7 @@ def test_halomod(model):
     z = 1.
     k = data_z1[:, 0] * cosmo['h']
     pk = data_z1[:, -1] / (cosmo['h']**3)
-    with pytest.warns(ccl.CCLDeprecationWarning):
+    with pytest.warns(ccl.CCLDeprecationWarning):  # TODO: remove normprof v3
         pk_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k, 1./(1+z), prf,
                                                   normprof1=True)
     tol = pk * HALOMOD_TOLERANCE
diff --git a/benchmarks/test_hod.py b/benchmarks/test_hod.py
index 623c52d2d..e9848f1d3 100644
--- a/benchmarks/test_hod.py
+++ b/benchmarks/test_hod.py
@@ -67,7 +67,7 @@ def _nz_2mrs(z):
     prf2pt = ccl.halos.Profile2ptHOD()
     # P(k)
     a_arr, lk_arr, _ = cosmo.get_linear_power().get_spline_arrays()
-    with pytest.warns(ccl.CCLDeprecationWarning):
+    with pytest.warns(ccl.CCLDeprecationWarning):  # TODO: remove normprof v3
         pk_hod = ccl.halos.halomod_Pk2D(cosmo, hmc, prf, prof_2pt=prf2pt,
                                         lk_arr=lk_arr, a_arr=a_arr,
                                         normprof1=True)
diff --git a/pyccl/tests/test_halomodel.py b/pyccl/tests/test_halomodel.py
index 39d42c26f..454f15d46 100644
--- a/pyccl/tests/test_halomodel.py
+++ b/pyccl/tests/test_halomodel.py
@@ -67,7 +67,7 @@ def get_pk_new(mf, c, cosmo, a, k, get_1h, get_2h):
     prf = ccl.halos.HaloProfileNFW(concentration=cc)
     hmc = ccl.halos.HMCalculator(mass_function=hmf, halo_bias=hbf,
                                  mass_def=mdef)
-    with pytest.warns(ccl.CCLDeprecationWarning):
+    with pytest.warns(ccl.CCLDeprecationWarning):  # TODO: remove normprof v3
         return ccl.halos.halomod_power_spectrum(cosmo, hmc, k, a, prf,
                                                 get_1h=get_1h,
                                                 get_2h=get_2h,
diff --git a/pyccl/tests/test_tkkssc.py b/pyccl/tests/test_tkkssc.py
index 05b502d3f..8d7252174 100644
--- a/pyccl/tests/test_tkkssc.py
+++ b/pyccl/tests/test_tkkssc.py
@@ -84,7 +84,7 @@ def test_tkkssc_linear_bias(isNC1, isNC2, isNC3, isNC4):
     assert np.allclose(tkkl_12, tkkl_34, atol=0, rtol=1e-12)
 
     # Test with the full T(k1,k2,a) for an NFW profile with bias ~1.
-    with pytest.warns(ccl.CCLDeprecationWarning):
+    with pytest.warns(ccl.CCLDeprecationWarning):  # TODO: remove normprof v3
         tkk = ccl.halos.halomod_Tk3D_SSC(
             COSMO, HMC, prof=NFW, lk_arr=np.log(KK), a_arr=AA,
             normprof1=True)

From ef00c6ca4f0dc85290390b37edd558a5416c96e3 Mon Sep 17 00:00:00 2001
From: David Alonso <dam.phys@gmail.com>
Date: Wed, 19 Apr 2023 15:04:34 +0100
Subject: [PATCH 5/9] comments

---
 pyccl/halos/halo_model.py            | 79 +++++++++++++++++-----------
 pyccl/halos/pk_1pt.py                |  2 +-
 pyccl/halos/pk_2pt.py                |  4 +-
 pyccl/halos/pk_4pt.py                | 18 +++----
 pyccl/halos/profiles/hod.py          | 18 +++----
 pyccl/halos/profiles/profile_base.py | 15 +++---
 6 files changed, 78 insertions(+), 58 deletions(-)

diff --git a/pyccl/halos/halo_model.py b/pyccl/halos/halo_model.py
index c0bf372d4..79782a6c5 100644
--- a/pyccl/halos/halo_model.py
+++ b/pyccl/halos/halo_model.py
@@ -82,8 +82,8 @@ def __init__(self, *, mass_function, halo_bias, mass_def=None,
             'log10M_min': log10M_min, 'log10M_max': log10M_max, 'nM': nM,
             'integration_method_M': integration_method_M, 'k_min': k_min}
         self._lmass = np.linspace(log10M_min, log10M_max, nM)
-        self.M = 10.**self._lmass
-        self._m0 = self.M[0]
+        self._mass = 10.**self._lmass
+        self._m0 = self._mass[0]
 
         if integration_method_M == "simpson":
             from scipy.integrate import simpson
@@ -106,8 +106,8 @@ def _integ_spline(self, fM, log10M):
     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:
-            self._mf = self.mass_function(cosmo, self.M, a)
-            integ = self._integrator(self._mf*self.M, self._lmass)
+            self._mf = self.mass_function(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
 
@@ -115,28 +115,47 @@ def _get_mass_function(self, cosmo, a, rho0):
     def _get_halo_bias(self, cosmo, a, rho0):
         # Compute the halo bias at this cosmo and a.
         if a != self._a_bf or cosmo != self._cosmo_bf:
-            self._bf = self.halo_bias(cosmo, self.M, a)
-            integ = self._integrator(self._mf*self._bf*self.M, self._lmass)
+            self._bf = self.halo_bias(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 update_ingredients(self, cosmo, a, *, get_bf):
+    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:
             self._get_halo_bias(cosmo, a, rho0)
 
-    def integrate_massfunc(self, array_2):
+    def _integrate_over_mf(self, array_2):
         #  ∫ dM n(M) f(M)
         i1 = self._integrator(self._mf * array_2, self._lmass)
         return i1 + self._mf0 * array_2[..., 0]
 
-    def integrate_massfunc_and_bias(self, array_2):
+    def _integrate_over_mbf(self, array_2):
         #  ∫ dM n(M) b(M) f(M)
         i1 = self._integrator(self._mf * self._bf * array_2, self._lmass)
         return i1 + self._mbf0 * array_2[..., 0]
 
+    def integrate_over_massfunc(self, func, cosmo, a):
+        """ Returns the integral over mass of a given funcion times
+        the mass function.
+
+        Args:
+            func (callable): a function accepting an array of halo masses
+                as a single argument, and returning an array of the
+                same size.
+            cosmo (:class:`~pyccl.core.Cosmology`): a Cosmology object.
+            a (float): scale factor.
+
+        Returns:
+            float or array_like: integral over mass of the function times
+                the mass function.
+        """
+        fM = func(self._mass)
+        hmc.update_ingredients(cosmo, a, get_bf=False)
+        return self._integrate_over_mf(fM)
+
     @deprecated()
     def profile_norm(self, cosmo, a, prof):
         """ Returns :math:`I^0_1(k\\rightarrow0,a|u)`
@@ -151,10 +170,10 @@ def profile_norm(self, cosmo, a, prof):
         Returns:
             float or array_like: integral value.
         """
-        self.update_ingredients(cosmo, a, get_bf=False)
+        self._get_ingredients(cosmo, a, get_bf=False)
         uk0 = prof.fourier(cosmo, self.precision['k_min'],
-                           self.M, a, mass_def=self.mass_def).T
-        return 1. / self.integrate_massfunc(uk0)
+                           self._mass, a, mass_def=self.mass_def).T
+        return 1. / self._integrate_over_mf(uk0)
 
     @warn_api(pairs=[("sel", "selection"),
                      ("amin", "a_min"),
@@ -210,8 +229,8 @@ def number_counts(self, cosmo, *, selection,
         # now do m intergrals in a loop
         mint = np.zeros_like(a)
         for i, _a in enumerate(a):
-            self.update_ingredients(cosmo, _a, get_bf=False)
-            _selm = np.atleast_2d(selection(self.M, _a)).T
+            self._get_ingredients(cosmo, _a, get_bf=False)
+            _selm = np.atleast_2d(selection(self._mass, _a)).T
             mint[i] = self._integrator(
                 dvda[i] * self._mf[..., :] * _selm[..., :],
                 self._lmass
@@ -241,9 +260,9 @@ def I_0_1(self, cosmo, k, a, prof):
             float or array_like: integral values evaluated at each
             value of `k`.
         """
-        self.update_ingredients(cosmo, a, get_bf=False)
-        uk = prof.fourier(cosmo, k, self.M, a, mass_def=self.mass_def).T
-        return self.integrate_massfunc(uk)
+        self._get_ingredients(cosmo, a, get_bf=False)
+        uk = prof.fourier(cosmo, k, self._mass, a, mass_def=self.mass_def).T
+        return self._integrate_over_mf(uk)
 
     def I_1_1(self, cosmo, k, a, prof):
         """ Solves the integral:
@@ -268,9 +287,9 @@ def I_1_1(self, cosmo, k, a, prof):
             float or array_like: integral values evaluated at each
             value of `k`.
         """
-        self.update_ingredients(cosmo, a, get_bf=True)
-        uk = prof.fourier(cosmo, k, self.M, a, mass_def=self.mass_def).T
-        return self.integrate_massfunc_and_bias(uk)
+        self._get_ingredients(cosmo, a, get_bf=True)
+        uk = prof.fourier(cosmo, k, self._mass, a, mass_def=self.mass_def).T
+        return self._integrate_over_mbf(uk)
 
     @warn_api(pairs=[("prof1", "prof")], reorder=["prof_2pt", "prof2"])
     def I_0_2(self, cosmo, k, a, prof, *, prof2=None, prof_2pt):
@@ -307,10 +326,10 @@ def I_0_2(self, cosmo, k, a, prof, *, prof2=None, prof_2pt):
         if prof2 is None:
             prof2 = prof
 
-        self.update_ingredients(cosmo, a, get_bf=False)
-        uk = prof_2pt.fourier_2pt(cosmo, k, self.M, a, prof,
+        self._get_ingredients(cosmo, a, get_bf=False)
+        uk = prof_2pt.fourier_2pt(cosmo, k, self._mass, a, prof,
                                   prof2=prof2, mass_def=self.mass_def).T
-        return self.integrate_massfunc(uk)
+        return self._integrate_over_mf(uk)
 
     @warn_api(pairs=[("prof1", "prof")], reorder=["prof_2pt", "prof2"])
     def I_1_2(self, cosmo, k, a, prof, *, prof2=None, prof_2pt):
@@ -346,10 +365,10 @@ def I_1_2(self, cosmo, k, a, prof, *, prof2=None, prof_2pt):
         if prof2 is None:
             prof2 = prof
 
-        self.update_ingredients(cosmo, a, get_bf=True)
-        uk = prof_2pt.fourier_2pt(cosmo, k, self.M, a, prof,
+        self._get_ingredients(cosmo, a, get_bf=True)
+        uk = prof_2pt.fourier_2pt(cosmo, k, self._mass, a, prof,
                                   prof2=prof2, mass_def=self.mass_def).T
-        return self.integrate_massfunc_and_bias(uk)
+        return self._integrate_over_mbf(uk)
 
     @warn_api(pairs=[("prof1", "prof")],
               reorder=["prof12_2pt", "prof2", "prof3", "prof34_2pt", "prof4"])
@@ -402,9 +421,9 @@ def I_0_22(self, cosmo, k, a, prof, *,
         if prof34_2pt is None:
             prof34_2pt = prof12_2pt
 
-        self.update_ingredients(cosmo, a, get_bf=False)
+        self._get_ingredients(cosmo, a, get_bf=False)
         uk12 = prof12_2pt.fourier_2pt(
-            cosmo, k, self.M, a, prof,
+            cosmo, k, self._mass, a, prof,
             prof2=prof2, mass_def=self.mass_def).T
 
         if (prof, prof2, prof12_2pt) == (prof3, prof4, prof34_2pt):
@@ -412,7 +431,7 @@ def I_0_22(self, cosmo, k, a, prof, *,
             uk34 = uk12
         else:
             uk34 = prof34_2pt.fourier_2pt(
-                cosmo, k, self.M, a, prof3,
+                cosmo, k, self._mass, a, prof3,
                 prof2=prof4, mass_def=self.mass_def).T
 
-        return self.integrate_massfunc(uk12[None, :, :] * uk34[:, None, :])
+        return self._integrate_over_mf(uk12[None, :, :] * uk34[:, None, :])
diff --git a/pyccl/halos/pk_1pt.py b/pyccl/halos/pk_1pt.py
index a432e5c5d..abc342ba5 100644
--- a/pyccl/halos/pk_1pt.py
+++ b/pyccl/halos/pk_1pt.py
@@ -38,7 +38,7 @@ 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 = prof.get_normalization(cosmo, aa, hmc) if normprof else 1
+        norm = prof.get_normalization(cosmo, aa, hmc=hmc) if normprof else 1
         # TODO: CCLv3 remove if
         out[ia] = i11 / norm
 
diff --git a/pyccl/halos/pk_2pt.py b/pyccl/halos/pk_2pt.py
index caa3335cc..e7a9d5768 100644
--- a/pyccl/halos/pk_2pt.py
+++ b/pyccl/halos/pk_2pt.py
@@ -113,13 +113,13 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof, *,
     out = np.zeros([na, nk])
     for ia, aa in enumerate(a_use):
         # normalizations
-        norm1 = prof.get_normalization(cosmo, aa, hmc) if normprof1 else 1
+        norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) if normprof1 else 1
         # TODO: CCLv3, remove if
 
         if prof2 == prof:
             norm2 = norm1
         else:
-            norm2 = prof2.get_normalization(cosmo, aa, hmc) if normprof2 else 1
+            norm2 = prof2.get_normalization(cosmo, aa, hmc=hmc) if normprof2 else 1
             # TODO: CCLv3, remove if
 
         if get_2h:
diff --git a/pyccl/halos/pk_4pt.py b/pyccl/halos/pk_4pt.py
index 753e047da..a39b8643a 100644
--- a/pyccl/halos/pk_4pt.py
+++ b/pyccl/halos/pk_4pt.py
@@ -89,25 +89,25 @@ def halomod_trispectrum_1h(cosmo, hmc, k, a, prof, *,
     out = np.zeros([na, nk, nk])
     for ia, aa in enumerate(a_use):
         # normalizations
-        norm1 = prof.get_normalization(cosmo, aa, hmc) if normprof1 else 1
+        norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) if normprof1 else 1
         # TODO: CCLv3 remove if
 
         if prof2 == prof:
             norm2 = norm1
         else:
-            norm2 = prof2.get_normalization(cosmo, aa, hmc) if normprof2 else 1
+            norm2 = prof2.get_normalization(cosmo, aa, hmc=hmc) if normprof2 else 1
             # TODO: CCLv3 remove if
 
         if prof3 == prof:
             norm3 = norm1
         else:
-            norm3 = prof3.get_normalization(cosmo, aa, hmc) if normprof3 else 1
+            norm3 = prof3.get_normalization(cosmo, aa, hmc=hmc) if normprof3 else 1
             # TODO: CCLv3 remove if
 
         if prof4 == prof2:
             norm4 = norm2
         else:
-            norm4 = prof4.get_normalization(cosmo, aa, hmc) if normprof4 else 1
+            norm4 = prof4.get_normalization(cosmo, aa, hmc=hmc) if normprof4 else 1
             # TODO: CCLv3 remove if
 
         # trispectrum
@@ -309,7 +309,7 @@ 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 = prof.get_normalization(cosmo, aa, hmc)**2
+        norm = prof.get_normalization(cosmo, aa, hmc=hmc)**2
         # TODO: CCLv3 remove if
         i12 = hmc.I_1_2(cosmo, k_use, aa, prof, prof2=prof, prof_2pt=prof_2pt)
 
@@ -452,7 +452,7 @@ 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 = prof.get_normalization(cosmo, aa, hmc) if normprof1 else 1
+        norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) if normprof1 else 1
         # TODO: CCLv3 remove if
         i11_1 = hmc.I_1_1(cosmo, k_use, aa, prof)
 
@@ -460,7 +460,7 @@ def halomod_Tk3D_SSC(
             norm2 = norm1
             i11_2 = i11_1
         else:
-            norm2 = prof2.get_normalization(cosmo, aa, hmc) if normprof2 else 1
+            norm2 = prof2.get_normalization(cosmo, aa, hmc=hmc) if normprof2 else 1
             # TODO: CCLv3 remove if
             i11_2 = hmc.I_1_1(cosmo, k_use, aa, prof2)
 
@@ -468,7 +468,7 @@ def halomod_Tk3D_SSC(
             norm3 = norm1
             i11_3 = i11_1
         else:
-            norm3 = prof3.get_normalization(cosmo, aa, hmc) if normprof3 else 1
+            norm3 = prof3.get_normalization(cosmo, aa, hmc=hmc) if normprof3 else 1
             # TODO: CCLv3 remove if
             i11_3 = hmc.I_1_1(cosmo, k_use, aa, prof3)
 
@@ -476,7 +476,7 @@ def halomod_Tk3D_SSC(
             norm4 = norm2
             i11_4 = i11_2
         else:
-            norm4 = prof4.get_normalization(cosmo, aa, hmc) if normprof4 else 1
+            norm4 = prof4.get_normalization(cosmo, aa, hmc=hmc) if normprof4 else 1
             # TODO: CCLv3 remove if
             i11_4 = hmc.I_1_1(cosmo, k_use, aa, prof4)
 
diff --git a/pyccl/halos/profiles/hod.py b/pyccl/halos/profiles/hod.py
index dfbd79666..75012fdf3 100644
--- a/pyccl/halos/profiles/hod.py
+++ b/pyccl/halos/profiles/hod.py
@@ -326,15 +326,15 @@ def get_normalization(self, cosmo, a, hmc):
         """Returns the normalization of this profile, which is the
         mean galaxy number density.
         """
-        Nc = self._Nc(hmc.M, a)
-        Ns = self._Ns(hmc.M, a)
-        fc = self._fc(a)
-        if self.ns_independent:
-            Ngal = Nc*fc + Ns
-        else:
-            Ngal = Nc*(fc + Ns)
-        hmc.update_ingredients(cosmo, a, get_bf=False)
-        return hmc.integrate_massfunc(Ngal)
+        def integ(M):
+            Nc = self._Nc(M, a)
+            Ns = self._Ns(M, a)
+            fc = self._fc(a)
+            if self.ns_independent:
+                return Nc*fc + Ns
+            else:
+                return Nc*(fc + Ns)
+        return hmc.integrate_over_massfunc(integ)
 
     def _fourier(self, cosmo, k, M, a, mass_def):
         M_use = np.atleast_1d(M)
diff --git a/pyccl/halos/profiles/profile_base.py b/pyccl/halos/profiles/profile_base.py
index 363f47121..70dfe6ff8 100644
--- a/pyccl/halos/profiles/profile_base.py
+++ b/pyccl/halos/profiles/profile_base.py
@@ -47,7 +47,7 @@ def __init__(self):
                             "_real or _fourier implementation.")
         self.precision_fftlog = FFTLogParams()
 
-    def get_normalization(self, cosmo, a, hmc):
+    def get_normalization(self, cosmo, a, *, hmc=None):
         """Profiles may be normalised by an overall function of redshift
         (or scale factor). This function may be cosmology dependent and
         often comes from integrating certain halo properties over mass.
@@ -65,13 +65,14 @@ def get_normalization(self, cosmo, a, hmc):
         Reurns:
             float: normalisation factor of this profile.
         """
-        uk0 = self.fourier(cosmo=cosmo, k=hmc.precision['k_min'],
-                           M=hmc.M, a=a, mass_def=hmc.mass_def)
-        hmc.update_ingredients(cosmo, a, get_bf=False)
-        return hmc.integrate_massfunc(uk0)
+        def integ(M):
+            return self.fourier(cosmo=cosmo,
+                                k=hmc.precision['k_min'],
+                                M=M, a=a, mass_def=hmc.mass_def)
+        return hmc.integrate_over_massfunc(integ)
         # TODO: CCLv3 replace by the below in v3 (profiles will all have a
         # default normalisation of 1. Normalisation will always be applied).
-        # return lambda *args, cosmo, a, **kwargs: 1.
+        # return 1.0
 
     @unlock_instance(mutate=True)
     @functools.wraps(FFTLogParams.update_parameters)
@@ -483,7 +484,7 @@ class HaloProfileNumberCounts(HaloProfile):
 class HaloProfileMatter(HaloProfile):
     """Base for matter halo profiles."""
 
-    def get_normalization(self, cosmo, a, hmc):
+    def get_normalization(self, cosmo, a, *, hmc=None):
         """Returns the normalization of all matter overdensity
         profiles, which is simply the comoving matter density.
         """

From a863b28de35c9f17377bff6b8b750a5ba665cbc1 Mon Sep 17 00:00:00 2001
From: David Alonso <dam.phys@gmail.com>
Date: Wed, 19 Apr 2023 15:10:22 +0100
Subject: [PATCH 6/9] flaked

---
 pyccl/halos/halo_model.py |  2 +-
 pyccl/halos/pk_2pt.py     |  3 ++-
 pyccl/halos/pk_4pt.py     | 24 ++++++++++++++++--------
 3 files changed, 19 insertions(+), 10 deletions(-)

diff --git a/pyccl/halos/halo_model.py b/pyccl/halos/halo_model.py
index 79782a6c5..7693bb98a 100644
--- a/pyccl/halos/halo_model.py
+++ b/pyccl/halos/halo_model.py
@@ -153,7 +153,7 @@ def integrate_over_massfunc(self, func, cosmo, a):
                 the mass function.
         """
         fM = func(self._mass)
-        hmc.update_ingredients(cosmo, a, get_bf=False)
+        self.update_ingredients(cosmo, a, get_bf=False)
         return self._integrate_over_mf(fM)
 
     @deprecated()
diff --git a/pyccl/halos/pk_2pt.py b/pyccl/halos/pk_2pt.py
index e7a9d5768..74f2659d5 100644
--- a/pyccl/halos/pk_2pt.py
+++ b/pyccl/halos/pk_2pt.py
@@ -119,7 +119,8 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof, *,
         if prof2 == prof:
             norm2 = norm1
         else:
-            norm2 = prof2.get_normalization(cosmo, aa, hmc=hmc) if normprof2 else 1
+            norm2 = prof2.get_normalization(cosmo, aa,
+                                            hmc=hmc) if normprof2 else 1
             # TODO: CCLv3, remove if
 
         if get_2h:
diff --git a/pyccl/halos/pk_4pt.py b/pyccl/halos/pk_4pt.py
index a39b8643a..1e46a4e3e 100644
--- a/pyccl/halos/pk_4pt.py
+++ b/pyccl/halos/pk_4pt.py
@@ -89,25 +89,29 @@ def halomod_trispectrum_1h(cosmo, hmc, k, a, prof, *,
     out = np.zeros([na, nk, nk])
     for ia, aa in enumerate(a_use):
         # normalizations
-        norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) if normprof1 else 1
+        norm1 = prof.get_normalization(cosmo, aa,
+                                       hmc=hmc) if normprof1 else 1
         # TODO: CCLv3 remove if
 
         if prof2 == prof:
             norm2 = norm1
         else:
-            norm2 = prof2.get_normalization(cosmo, aa, hmc=hmc) if normprof2 else 1
+            norm2 = prof2.get_normalization(cosmo, aa,
+                                            hmc=hmc) if normprof2 else 1
             # TODO: CCLv3 remove if
 
         if prof3 == prof:
             norm3 = norm1
         else:
-            norm3 = prof3.get_normalization(cosmo, aa, hmc=hmc) if normprof3 else 1
+            norm3 = prof3.get_normalization(cosmo, aa,
+                                            hmc=hmc) if normprof3 else 1
             # TODO: CCLv3 remove if
 
         if prof4 == prof2:
             norm4 = norm2
         else:
-            norm4 = prof4.get_normalization(cosmo, aa, hmc=hmc) if normprof4 else 1
+            norm4 = prof4.get_normalization(cosmo, aa,
+                                            hmc=hmc) if normprof4 else 1
             # TODO: CCLv3 remove if
 
         # trispectrum
@@ -452,7 +456,8 @@ 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 = prof.get_normalization(cosmo, aa, hmc=hmc) if normprof1 else 1
+        norm1 = prof.get_normalization(cosmo, aa,
+                                       hmc=hmc) if normprof1 else 1
         # TODO: CCLv3 remove if
         i11_1 = hmc.I_1_1(cosmo, k_use, aa, prof)
 
@@ -460,7 +465,8 @@ def halomod_Tk3D_SSC(
             norm2 = norm1
             i11_2 = i11_1
         else:
-            norm2 = prof2.get_normalization(cosmo, aa, hmc=hmc) if normprof2 else 1
+            norm2 = prof2.get_normalization(cosmo, aa,
+                                            hmc=hmc) if normprof2 else 1
             # TODO: CCLv3 remove if
             i11_2 = hmc.I_1_1(cosmo, k_use, aa, prof2)
 
@@ -468,7 +474,8 @@ def halomod_Tk3D_SSC(
             norm3 = norm1
             i11_3 = i11_1
         else:
-            norm3 = prof3.get_normalization(cosmo, aa, hmc=hmc) if normprof3 else 1
+            norm3 = prof3.get_normalization(cosmo, aa,
+                                            hmc=hmc) if normprof3 else 1
             # TODO: CCLv3 remove if
             i11_3 = hmc.I_1_1(cosmo, k_use, aa, prof3)
 
@@ -476,7 +483,8 @@ def halomod_Tk3D_SSC(
             norm4 = norm2
             i11_4 = i11_2
         else:
-            norm4 = prof4.get_normalization(cosmo, aa, hmc=hmc) if normprof4 else 1
+            norm4 = prof4.get_normalization(cosmo, aa,
+                                            hmc=hmc) if normprof4 else 1
             # TODO: CCLv3 remove if
             i11_4 = hmc.I_1_1(cosmo, k_use, aa, prof4)
 

From 36529ab17938b4ae4dc82c42225ac5f2cf8dda5b Mon Sep 17 00:00:00 2001
From: David Alonso <dam.phys@gmail.com>
Date: Wed, 19 Apr 2023 15:11:28 +0100
Subject: [PATCH 7/9] comments

---
 pyccl/halos/profiles/pressure_gnfw.py |  4 ++--
 pyccl/halos/profiles/profile_base.py  | 10 +++++-----
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/pyccl/halos/profiles/pressure_gnfw.py b/pyccl/halos/profiles/pressure_gnfw.py
index 6ade1d277..347a296af 100644
--- a/pyccl/halos/profiles/pressure_gnfw.py
+++ b/pyccl/halos/profiles/pressure_gnfw.py
@@ -182,8 +182,8 @@ def integrand(x):
         return Fq
 
     def _norm(self, cosmo, M, a, mb):
-        # Computes the normalisation factor of the GNFW profile.
-        # Normalisation factor is given in units of eV/cm^3.
+        # Computes the normalization factor of the GNFW profile.
+        # Normalization factor is given in units of eV/cm^3.
         # (Bolliet et al. 2017).
         h70 = cosmo["h"]/0.7
         C0 = 1.65*h70**2
diff --git a/pyccl/halos/profiles/profile_base.py b/pyccl/halos/profiles/profile_base.py
index 70dfe6ff8..f2928fd21 100644
--- a/pyccl/halos/profiles/profile_base.py
+++ b/pyccl/halos/profiles/profile_base.py
@@ -48,11 +48,11 @@ def __init__(self):
         self.precision_fftlog = FFTLogParams()
 
     def get_normalization(self, cosmo, a, *, hmc=None):
-        """Profiles may be normalised by an overall function of redshift
+        """Profiles may be normalized by an overall function of redshift
         (or scale factor). This function may be cosmology dependent and
         often comes from integrating certain halo properties over mass.
-        This method returns this normalising factor. For example,
-        to get the normalised profile in real space, one would call
+        This method returns this normalizing factor. For example,
+        to get the normalized profile in real space, one would call
         the `real` method, and then **divide** the result by the value
         returned by this method.
 
@@ -63,7 +63,7 @@ def get_normalization(self, cosmo, a, *, hmc=None):
             a (float): scale factor.
 
         Reurns:
-            float: normalisation factor of this profile.
+            float: normalization factor of this profile.
         """
         def integ(M):
             return self.fourier(cosmo=cosmo,
@@ -71,7 +71,7 @@ def integ(M):
                                 M=M, a=a, mass_def=hmc.mass_def)
         return hmc.integrate_over_massfunc(integ)
         # TODO: CCLv3 replace by the below in v3 (profiles will all have a
-        # default normalisation of 1. Normalisation will always be applied).
+        # default normalization of 1. Normalization will always be applied).
         # return 1.0
 
     @unlock_instance(mutate=True)

From 7e334834c99255d2782d542f8462d17b82786fe8 Mon Sep 17 00:00:00 2001
From: David Alonso <dam.phys@gmail.com>
Date: Wed, 19 Apr 2023 15:23:17 +0100
Subject: [PATCH 8/9] bad syntax

---
 pyccl/halos/profiles/hod.py          | 2 +-
 pyccl/halos/profiles/profile_base.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/pyccl/halos/profiles/hod.py b/pyccl/halos/profiles/hod.py
index 75012fdf3..569f6ee11 100644
--- a/pyccl/halos/profiles/hod.py
+++ b/pyccl/halos/profiles/hod.py
@@ -334,7 +334,7 @@ def integ(M):
                 return Nc*fc + Ns
             else:
                 return Nc*(fc + Ns)
-        return hmc.integrate_over_massfunc(integ)
+        return hmc.integrate_over_massfunc(integ, cosmo, a)
 
     def _fourier(self, cosmo, k, M, a, mass_def):
         M_use = np.atleast_1d(M)
diff --git a/pyccl/halos/profiles/profile_base.py b/pyccl/halos/profiles/profile_base.py
index f2928fd21..02a67472e 100644
--- a/pyccl/halos/profiles/profile_base.py
+++ b/pyccl/halos/profiles/profile_base.py
@@ -69,7 +69,7 @@ def integ(M):
             return self.fourier(cosmo=cosmo,
                                 k=hmc.precision['k_min'],
                                 M=M, a=a, mass_def=hmc.mass_def)
-        return hmc.integrate_over_massfunc(integ)
+        return hmc.integrate_over_massfunc(integ, cosmo, a)
         # TODO: CCLv3 replace by the below in v3 (profiles will all have a
         # default normalization of 1. Normalization will always be applied).
         # return 1.0

From 3f1485ed107730e8eded0abba7a39245bb23f0a7 Mon Sep 17 00:00:00 2001
From: David Alonso <dam.phys@gmail.com>
Date: Wed, 19 Apr 2023 15:24:35 +0100
Subject: [PATCH 9/9] old

---
 pyccl/halos/halo_model.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyccl/halos/halo_model.py b/pyccl/halos/halo_model.py
index 7693bb98a..42f4668ee 100644
--- a/pyccl/halos/halo_model.py
+++ b/pyccl/halos/halo_model.py
@@ -153,7 +153,7 @@ def integrate_over_massfunc(self, func, cosmo, a):
                 the mass function.
         """
         fM = func(self._mass)
-        self.update_ingredients(cosmo, a, get_bf=False)
+        self._get_ingredients(cosmo, a, get_bf=False)
         return self._integrate_over_mf(fM)
 
     @deprecated()