HOD profile: relax requirement to only form satellites when centrals …

…are present (#875)

* implemented r_corr in 1h 2-pt correlator

* HOD profile: decouple satellite formation from centrals (ns_independent)

* HOD profile: decouple satellite formation from centrals (ns_independent)

* *ignore*: reverted r_corr changes

* *ignore*: reverted r_corr changes

* test of HOD ns_independent

* HOD ns_indep r_scale smaller than r_truncated

* more coverage

* flaked

pyccl/halos/profiles.py

@@ -1165,7 +1165,9 @@ class HaloProfileHOD(HaloProfile):
         bmax_p (float): tilt parameter for
             :math:`\\beta_{\\rm max}`.
         a_pivot (float): pivot scale factor :math:`a_*`.
-    """
+        ns_independent (bool): drop requirement to only form
+            satellites when centrals are present.
+        """
     name = 'HOD'
 
     def __init__(self, c_M_relation,
@@ -1174,7 +1176,7 @@ def __init__(self, c_M_relation,
                  lM1_0=13.3, lM1_p=0., alpha_0=1.,
                  alpha_p=0., fc_0=1., fc_p=0.,
                  bg_0=1., bg_p=0., bmax_0=1., bmax_p=0.,
-                 a_pivot=1.):
+                 a_pivot=1., ns_independent=False):
         if not isinstance(c_M_relation, Concentration):
             raise TypeError("c_M_relation must be of type `Concentration`)")
 
@@ -1196,6 +1198,7 @@ def __init__(self, c_M_relation,
         self.bmax_0 = bmax_0
         self.bmax_p = bmax_p
         self.a_pivot = a_pivot
+        self.ns_independent = ns_independent
         super(HaloProfileHOD, self).__init__()
 
     def _get_cM(self, cosmo, M, a, mdef=None):
@@ -1209,7 +1212,8 @@ def update_parameters(self, lMmin_0=None, lMmin_p=None,
                           fc_0=None, fc_p=None,
                           bg_0=None, bg_p=None,
                           bmax_0=None, bmax_p=None,
-                          a_pivot=None):
+                          a_pivot=None,
+                          ns_independent=None):
         """ Update any of the parameters associated with
         this profile. Any parameter set to `None` won't be updated.
 
@@ -1247,6 +1251,8 @@ def update_parameters(self, lMmin_0=None, lMmin_p=None,
             bmax_p (float): tilt parameter for
                 :math:`\\beta_{\\rm max}`.
             a_pivot (float): pivot scale factor :math:`a_*`.
+            ns_independent (bool): drop requirement to only form
+                satellites when centrals are present
         """
         if lMmin_0 is not None:
             self.lMmin_0 = lMmin_0
@@ -1282,6 +1288,8 @@ def update_parameters(self, lMmin_0=None, lMmin_p=None,
             self.bmax_p = bmax_p
         if a_pivot is not None:
             self.a_pivot = a_pivot
+        if ns_independent is not None:
+            self.ns_independent = ns_independent
 
     def _usat_real(self, cosmo, r, M, a, mass_def):
         r_use = np.atleast_1d(r)
@@ -1346,7 +1354,10 @@ def _real(self, cosmo, r, M, a, mass_def):
         # NFW profile
         ur = self._usat_real(cosmo, r_use, M_use, a, mass_def)
 
-        prof = Nc[:, None] * (fc + Ns[:, None] * ur)
+        if self.ns_independent:
+            prof = Nc[:, None] * fc + Ns[:, None] * ur
+        else:
+            prof = Nc[:, None] * (fc + Ns[:, None] * ur)
 
         if np.ndim(r) == 0:
             prof = np.squeeze(prof, axis=-1)
@@ -1364,7 +1375,10 @@ def _fourier(self, cosmo, k, M, a, mass_def):
         # NFW profile
         uk = self._usat_fourier(cosmo, k_use, M_use, a, mass_def)
 
-        prof = Nc[:, None] * (fc + Ns[:, None] * uk)
+        if self.ns_independent:
+            prof = Nc[:, None] * fc + Ns[:, None] * uk
+        else:
+            prof = Nc[:, None] * (fc + Ns[:, None] * uk)
 
         if np.ndim(k) == 0:
             prof = np.squeeze(prof, axis=-1)
@@ -1384,7 +1398,10 @@ def _fourier_variance(self, cosmo, k, M, a, mass_def):
         uk = self._usat_fourier(cosmo, k_use, M_use, a, mass_def)
 
         prof = Ns[:, None] * uk
-        prof = Nc[:, None] * (2 * fc * prof + prof**2)
+        if self.ns_independent:
+            prof = 2 * Nc[:, None] * fc * prof + prof**2
+        else:
+            prof = Nc[:, None] * (2 * fc * prof + prof**2)
 
         if np.ndim(k) == 0:
             prof = np.squeeze(prof, axis=-1)

pyccl/tests/test_profiles.py

@@ -120,6 +120,42 @@ def test_hod_smoke():
         assert getattr(p, n) == 1234.
 
 
+@pytest.mark.parametrize('real_prof', [True, False])
+def test_hod_ns_independent(real_prof):
+    def func(prof):
+        return prof._real if real_prof else prof._fourier
+
+    c = ccl.halos.ConcentrationDuffy08(M200)
+    hmd = c.mdef
+    p1 = ccl.halos.HaloProfileHOD(c_M_relation=c,
+                                  lMmin_0=12.,
+                                  ns_independent=False)
+    p2 = ccl.halos.HaloProfileHOD(c_M_relation=c,
+                                  lMmin_0=12.,
+                                  ns_independent=True)
+    # M < Mmin
+    f1 = func(p1)(COSMO, 0.01, 1e10, 1., hmd)
+    assert np.all(f1 == 0)
+    f2 = func(p2)(COSMO, 0.01, 1e10, 1., hmd)
+    assert np.all(f2 > 0)
+    # M > Mmin
+    f1 = func(p1)(COSMO, 0.01, 1e14, 1., hmd)
+    f2 = func(p2)(COSMO, 0.01, 1e14, 1., hmd)
+    assert np.allclose(f1, f2, rtol=0)
+    # M == Mmin
+    f1 = func(p1)(COSMO, 0.01, 1e12, 1., hmd)
+    f2 = func(p2)(COSMO, 0.01, 1e12, 1., hmd)
+    assert np.allclose(2*f1, f2+0.5, rtol=0)
+
+    if not real_prof:
+        f1 = p1._fourier_variance(COSMO, 0.01, 1e10, 1., hmd)
+        f2 = p2._fourier_variance(COSMO, 0.01, 1e10, 1., hmd)
+        assert f2 > f1 == 0
+
+    p1.update_parameters(ns_independent=True)
+    assert p1.ns_independent is True
+
+
 def test_hod_2pt_raises():
     pbad = ccl.halos.HaloProfilePressureGNFW()
     c = ccl.halos.ConcentrationDuffy08(M200)