API restoration (#1083)

* up to test_halomodel

* bug in deprecated halo model. better error message

* mass def fixing

* SSC bug

* add ct test

* broken test

pyccl/cosmology.py

@@ -694,7 +694,7 @@ def _get_halo_model_nonlin_power(self):
         prf = hal.HaloProfileNFW(concentration=c)
         hmc = hal.HMCalculator(mass_function=hmf, halo_bias=hbf,
                                mass_def=mdef)
-        return hal.halomod_Pk2D(self, hmc, prf)
+        return hal.halomod_Pk2D(self, hmc, prf, normprof1=True)
 
     @cache(maxsize=3)
     def _compute_nonlin_power(self):

pyccl/halos/halo_model.py

@@ -57,8 +57,12 @@ def __init__(self, *, mass_function, halo_bias, mass_def=None,
                  log10M_min=8., log10M_max=16., nM=128,
                  integration_method_M='simpson', k_min=1E-5):
         # Initialize halo model ingredients.
-        self.mass_def, self.mass_function, self.halo_bias = MassDef.from_specs(
-            mass_def, mass_function=mass_function, halo_bias=halo_bias)
+        out = MassDef.from_specs(mass_def, mass_function=mass_function,
+                                 halo_bias=halo_bias)
+        if len(out) != 3:
+            raise ValueError("A valid mass function and halo bias is "
+                             "needed")
+        self.mass_def, self.mass_function, self.halo_bias = out
 
         self.precision = {
             'log10M_min': log10M_min, 'log10M_max': log10M_max, 'nM': nM,
@@ -79,6 +83,14 @@ def __init__(self, *, mass_function, halo_bias, mass_def=None,
         self._cosmo_mf = self._cosmo_bf = None
         self._a_mf = self._a_bf = -1
 
+    def _fix_profile_mass_def(self, prof):
+        # TODO v3: remove this (in v3 all profiles have a mass_def).
+        # If profile has no mass definition assigned, assign one.
+        if prof.mass_def is None:
+            warnings.warn("In v3 all profiles will need an associated "
+                          "mass definition.", CCLDeprecationWarning)
+            prof.mass_def = self.mass_def
+
     def _integ_spline(self, fM, log10M):
         # Spline integrator
         return _spline_integrate(log10M, fM, log10M[0], log10M[-1])
@@ -156,6 +168,7 @@ def profile_norm(self, cosmo, a, prof):
         Returns:
             float or array_like: integral value.
         """
+        self._fix_profile_mass_def(prof)
         self._check_mass_def(prof)
         self._get_ingredients(cosmo, a, get_bf=False)
         uk0 = prof.fourier(cosmo, self.precision['k_min'], self._mass, a).T
@@ -246,6 +259,7 @@ def I_0_1(self, cosmo, k, a, prof):
             float or array_like: integral values evaluated at each
             value of `k`.
         """
+        self._fix_profile_mass_def(prof)
         self._check_mass_def(prof)
         self._get_ingredients(cosmo, a, get_bf=False)
         uk = prof.fourier(cosmo, k, self._mass, a).T
@@ -274,6 +288,7 @@ def I_1_1(self, cosmo, k, a, prof):
             float or array_like: integral values evaluated at each
             value of `k`.
         """
+        self._fix_profile_mass_def(prof)
         self._check_mass_def(prof)
         self._get_ingredients(cosmo, a, get_bf=True)
         uk = prof.fourier(cosmo, k, self._mass, a).T
@@ -313,6 +328,8 @@ def I_0_2(self, cosmo, k, a, prof, *, prof2=None, prof_2pt):
         """
         if prof2 is None:
             prof2 = prof
+        self._fix_profile_mass_def(prof)
+        self._fix_profile_mass_def(prof2)
 
         self._check_mass_def(prof, prof2)
         self._get_ingredients(cosmo, a, get_bf=False)
@@ -352,6 +369,8 @@ def I_1_2(self, cosmo, k, a, prof, *, prof2=None, prof_2pt):
         """
         if prof2 is None:
             prof2 = prof
+        self._fix_profile_mass_def(prof)
+        self._fix_profile_mass_def(prof2)
 
         self._check_mass_def(prof, prof2)
         self._get_ingredients(cosmo, a, get_bf=True)
@@ -409,6 +428,10 @@ def I_0_22(self, cosmo, k, a, prof, *,
         if prof34_2pt is None:
             prof34_2pt = prof12_2pt
 
+        self._fix_profile_mass_def(prof)
+        self._fix_profile_mass_def(prof2)
+        self._fix_profile_mass_def(prof3)
+        self._fix_profile_mass_def(prof4)
         self._check_mass_def(prof, prof2, prof3, prof4)
         self._get_ingredients(cosmo, a, get_bf=False)
         uk12 = prof12_2pt.fourier_2pt(

pyccl/halos/pk_1pt.py

@@ -60,6 +60,7 @@ def halomod_mean_profile_1pt(cosmo, hmc, k, a, prof, *, normprof=None):
     :math:`\\langle u(k,a|M)\\rangle` is the halo profile as a
     function of scale, scale factor and halo mass.
     """
+    hmc._fix_profile_mass_def(prof)
     return _Ix1("I_0_1", cosmo, hmc, k, a, prof, normprof)
 
 
@@ -76,6 +77,7 @@ def halomod_bias_1pt(cosmo, hmc, k, a, prof, *, normprof=None):
     :math:`\\langle u(k,a|M)\\rangle` is the halo profile as a
     function of scale, scale factor and halo mass.
     """
+    hmc._fix_profile_mass_def(prof)
     return _Ix1("I_1_1", cosmo, hmc, k, a, prof, normprof)
 
 

pyccl/halos/pk_2pt.py

@@ -104,6 +104,8 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof, *,
         prof2 = prof
     if prof_2pt is None:
         prof_2pt = Profile2pt()
+    hmc._fix_profile_mass_def(prof)
+    hmc._fix_profile_mass_def(prof2)
 
     pk2d = cosmo.parse_pk(p_of_k_a)
     extrap = cosmo if extrap_pk else None  # extrapolation rule for pk2d

pyccl/halos/pk_4pt.py

@@ -7,7 +7,6 @@
 
 from .. import CCLWarning, Tk3D, warn_api
 from . import HaloProfileNFW, Profile2pt
-from . import HaloProfileNumberCounts as ProfNC
 
 
 @warn_api(pairs=[("prof1", "prof")], reorder=["prof12_2pt", "prof3", "prof4"])
@@ -80,6 +79,10 @@ 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)
+    hmc._fix_profile_mass_def(prof)
+    hmc._fix_profile_mass_def(prof2)
+    hmc._fix_profile_mass_def(prof3)
+    hmc._fix_profile_mass_def(prof4)
 
     na = len(a_use)
     nk = len(k_use)
@@ -292,6 +295,7 @@ def halomod_Tk3D_SSC_linear_bias(cosmo, hmc, *, prof,
 
     if not isinstance(prof, HaloProfileNFW):
         raise TypeError("prof should be HaloProfileNFW.")
+    hmc._fix_profile_mass_def(prof)
 
     # Make sure biases are of the form number of a x number of k
     ones = np.ones_like(a_arr)
@@ -445,6 +449,10 @@ 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)
+    hmc._fix_profile_mass_def(prof)
+    hmc._fix_profile_mass_def(prof2)
+    hmc._fix_profile_mass_def(prof3)
+    hmc._fix_profile_mass_def(prof4)
 
     k_use = np.exp(lk_arr)
     pk2d = cosmo.parse_pk(p_of_k_a)
@@ -506,19 +514,19 @@ def halomod_Tk3D_SSC(
         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 pA.is_number_counts else np.zeros_like(k_use)
+            bB = i11_B / nB if pB.is_number_counts else np.zeros_like(k_use)
             i02 = hmc.I_0_2(cosmo, k_use, aa, pA, prof2=pB, prof_2pt=p2pt)
             P = (pk * i11_A * i11_B + i02) / (nA * nB)
             return (bA + bB) * P
 
-        if isinstance(prof, ProfNC) or isinstance(prof2, ProfNC):
+        if prof.is_number_counts or prof2.is_number_counts:
             dpk12[ia] -= _get_counterterm(prof, prof2, prof12_2pt,
                                           norm1, norm2, i11_1, i11_2)
 
-        if isinstance(prof3, ProfNC) or isinstance(prof4, ProfNC):
+        if prof3.is_number_counts or prof4.is_number_counts:
             if (prof, prof2, prof12_2pt) == (prof3, prof4, prof34_2pt):
-                dpk34[ia] -= dpk12[ia]
+                dpk34[ia] = dpk12[ia]
             else:
                 dpk34[ia] -= _get_counterterm(prof3, prof4, prof34_2pt,
                                               norm3, norm4, i11_3, i11_4)

pyccl/pk2d.py

@@ -601,7 +601,13 @@ def parse_pk2d(cosmo, p_of_k_a=DEFAULT_POWER_SPECTRUM, *, is_linear=False):
     if isinstance(p_of_k_a, Pk2D):
         psp = p_of_k_a.psp
     else:
-        if p_of_k_a is None or isinstance(p_of_k_a, str):
+        if p_of_k_a is None:
+            warnings.warn("The default power spectrum can is now designated "
+                          "via ``ccl.DEFAULT_POWER_SPECTRUM``. The use of "
+                          "``None`` will be deprecated in future versions.",
+                          CCLDeprecationWarning)
+            name = DEFAULT_POWER_SPECTRUM
+        elif isinstance(p_of_k_a, str):
             name = p_of_k_a
         else:
             raise ValueError("p_of_k_a must be a pyccl.Pk2D object, "

pyccl/pyutils.py

@@ -10,7 +10,7 @@
 
 import numpy as np
 
-from . import CCLError, lib, spline_params
+from . import CCLError, lib, spline_params, deprecated
 
 
 NoneArr = np.array([])
@@ -677,3 +677,22 @@ def check_openmp_threads():
     """
 
     return lib.openmp_threads()
+
+
+@deprecated
+def assert_warns(wtype, f, *args, **kwargs):
+    """Check that a function call `f(*args, **kwargs)` raises a warning of
+    type wtype.
+    Returns the output of `f(*args, **kwargs)` unless there was no warning,
+    in which case an AssertionError is raised.
+    """
+    import warnings
+    # Check that f() raises a warning, but not an error.
+    with warnings.catch_warnings(record=True) as w:
+        warnings.simplefilter("always")
+        res = f(*args, **kwargs)
+    assert len(w) >= 1, "Expected warning was not raised."
+    assert issubclass(w[0].category, wtype), \
+        "Warning raised was the wrong type (got %s, expected %s)" % (
+            w[0].category, wtype)
+    return res

pyccl/tests/test_power.py

@@ -25,7 +25,7 @@ def test_halomod_f2d_copy():
     prf = ccl.halos.HaloProfileNFW(concentration=cc)
     hmc = ccl.halos.HMCalculator(mass_function=hmf, halo_bias=hbf,
                                  mass_def=mdef)
-    pk2d = ccl.halos.halomod_Pk2D(COSMO_HM, hmc, prf)
+    pk2d = ccl.halos.halomod_Pk2D(COSMO_HM, hmc, prf, normprof1=True)
     psp_new = pk2d.psp
     # This just triggers the internal calculation
     with pytest.warns(ccl.CCLDeprecationWarning):
@@ -68,7 +68,8 @@ def test_nonlin_matter_power_halomod(k):
     prf = ccl.halos.HaloProfileNFW(concentration=cc)
     hmc = ccl.halos.HMCalculator(mass_function=hmf, halo_bias=hbf,
                                  mass_def=mdef)
-    pkb = ccl.halos.halomod_power_spectrum(COSMO_HM, hmc, k, a, prf)
+    pkb = ccl.halos.halomod_power_spectrum(COSMO_HM, hmc, k, a, prf,
+                                           normprof1=True)
 
     assert np.allclose(pk, pkb)
     assert np.all(np.isfinite(pk))

pyccl/tests/test_tkkssc.py

@@ -16,6 +16,8 @@
 
 NFW = ccl.halos.HaloProfileNFW(concentration=CON, fourier_analytic=True)
 HOD = ccl.halos.HaloProfileHOD(concentration=CON)
+HOD_nogc = ccl.halos.HaloProfileHOD(concentration=CON)
+HOD_nogc.is_number_counts = False
 GNFW = ccl.halos.HaloProfilePressureGNFW(mass_def=M200)
 
 PKC = ccl.halos.Profile2pt()
@@ -128,3 +130,119 @@ def test_tkkssc_linear_bias_raises():
     """Test that it raises if the profile is not NFW."""
     with pytest.raises(TypeError):
         ccl.halos.halomod_Tk3D_SSC_linear_bias(COSMO, HMC, prof=GNFW)
+
+
+def get_ssc_counterterm_gc(k, a, hmc, prof1, prof2, prof12_2pt,
+                           normalize=False):
+    P_12 = b1 = b2 = np.zeros_like(k)
+    if prof1.is_number_counts or prof2.is_number_counts:
+        norm1 = 1./prof1.get_normalization(COSMO, a, hmc=hmc)
+        norm2 = 1./prof2.get_normalization(COSMO, a, hmc=hmc)
+        norm12 = 1
+        if prof1.is_number_counts or normalize:
+            norm12 *= norm1
+        if prof2.is_number_counts or normalize:
+            norm12 *= norm2
+
+        i11_1 = hmc.I_1_1(COSMO, k, a, prof1)
+        i11_2 = hmc.I_1_1(COSMO, k, a, prof2)
+        i02_12 = hmc.I_0_2(COSMO, k, a, prof1,
+                           prof_2pt=prof12_2pt, prof2=prof2)
+
+        pk = ccl.linear_matter_power(COSMO, k, a)
+        P_12 = norm12 * (pk * i11_1 * i11_2 + i02_12)
+
+        if prof1.is_number_counts:
+            b1 = ccl.halos.halomod_bias_1pt(COSMO, hmc, k, a, prof1) * norm1
+        if prof2.is_number_counts:
+            b2 = ccl.halos.halomod_bias_1pt(COSMO, hmc, k, a, prof2) * norm2
+
+    return (b1 + b2) * P_12
+
+
+@pytest.mark.parametrize('kwargs', [
+                         # All is_number_counts = False
+                         {'prof': NFW, 'prof2': NFW,
+                          'prof3': NFW, 'prof4': NFW},
+                         {'prof': HOD, 'prof2': NFW,
+                          'prof3': NFW, 'prof4': NFW},
+                         {'prof': HOD, 'prof2': HOD,
+                          'prof3': NFW, 'prof4': NFW},
+                         {'prof': HOD, 'prof2': HOD,
+                          'prof3': HOD, 'prof4': NFW},
+                         {'prof': NFW, 'prof2': NFW,
+                          'prof3': HOD, 'prof4': HOD},
+                         {'prof': HOD, 'prof2': NFW,
+                          'prof3': HOD, 'prof4': HOD},
+                         {'prof': HOD, 'prof2': HOD,
+                          'prof3': NFW, 'prof4': HOD},
+                         {'prof': HOD, 'prof2': None,
+                          'prof3': NFW, 'prof4': None},
+                         {'prof': HOD, 'prof2': None,
+                          'prof3': None, 'prof4': None},
+                         {'prof': HOD, 'prof2': HOD,
+                          'prof3': None, 'prof4': None},
+                         # As in benchmarks/test_covariances.py
+                         {'prof': NFW, 'prof2': NFW,
+                          'prof3': None, 'prof4': None},
+                         # Setting prof34_2pt
+                         {'prof': NFW, 'prof2': NFW,
+                          'prof3': None, 'prof4': None,
+                          'prof34_2pt': PKC},
+                         # All is_number_counts = True
+                         {'prof': HOD, 'prof2': HOD,
+                          'prof3': HOD, 'prof4': HOD},
+                         ]
+                         )
+def test_tkkssc_counterterms_gc(kwargs):
+    k_arr = KK
+    a_arr = np.array([0.3, 0.5, 0.7, 1.0])
+
+    # Tk's without clustering terms. Set is_number_counts=False for HOD
+    # profiles
+    # Ensure HOD profiles are normalized
+    kwargs_nogc = kwargs.copy()
+    keys = list(kwargs.keys())
+    for k in keys:
+        v = kwargs[k]
+        if isinstance(v, ccl.halos.HaloProfileHOD):
+            kwargs_nogc[k] = HOD_nogc
+            normname = 'norm'+k
+            if k == 'prof':
+                normname = 'normprof1'
+            kwargs_nogc[normname] = True
+            kwargs[normname] = True
+    tkk_nogc = ccl.halos.halomod_Tk3D_SSC(COSMO, HMC,
+                                          lk_arr=np.log(k_arr), a_arr=a_arr,
+                                          **kwargs_nogc)
+    _, _, _, tkk_nogc_arrs = tkk_nogc.get_spline_arrays()
+    tk_nogc_12, tk_nogc_34 = tkk_nogc_arrs
+
+    # Tk's with clustering terms
+    tkk_gc = ccl.halos.halomod_Tk3D_SSC(COSMO, HMC,
+                                        lk_arr=np.log(k_arr), a_arr=a_arr,
+                                        **kwargs)
+    _, _, _, tkk_gc_arrs = tkk_gc.get_spline_arrays()
+    tk_gc_12, tk_gc_34 = tkk_gc_arrs
+
+    # Update the None's to their corresponding values
+    if kwargs['prof2'] is None:
+        kwargs['prof2'] = kwargs['prof']
+    if kwargs['prof3'] is None:
+        kwargs['prof3'] = kwargs['prof']
+    if kwargs['prof4'] is None:
+        kwargs['prof4'] = kwargs['prof2']
+
+    # Tk's of the clustering terms
+    tkc12 = []
+    tkc34 = []
+    for aa in a_arr:
+        tkc12.append(get_ssc_counterterm_gc(k_arr, aa, HMC, kwargs['prof'],
+                                            kwargs['prof2'], PKC))
+        tkc34.append(get_ssc_counterterm_gc(k_arr, aa, HMC, kwargs['prof3'],
+                                            kwargs['prof4'], PKC))
+    tkc12 = np.array(tkc12)
+    tkc34 = np.array(tkc34)
+
+    assert np.abs((tk_nogc_12 - tkc12) / tk_gc_12 - 1).max() < 1e-5
+    assert np.abs((tk_nogc_34 - tkc34) / tk_gc_34 - 1).max() < 1e-5