Merge pull request #958 from LSSTDESC/analytic_hernquist_2D_densities

Add analytic _projected, _cumul2d and _fourier to Hernquist profile

pyccl/halos/profiles.py

@@ -639,7 +639,7 @@ class HaloProfileNFW(HaloProfile):
             relation to use with this profile.
         fourier_analytic (bool): set to `True` if you want to compute
             the Fourier profile analytically (and not through FFTLog).
-            Default: `False`.
+            Default: `True`.
         projected_analytic (bool): set to `True` if you want to
             compute the 2D projected profile analytically (and not
             through FFTLog). Default: `False`.
@@ -919,21 +919,48 @@ class HaloProfileHernquist(HaloProfile):
     Args:
         c_M_relation (:obj:`Concentration`): concentration-mass
             relation to use with this profile.
+        fourier_analytic (bool): set to `True` if you want to compute
+            the Fourier profile analytically (and not through FFTLog).
+            Default: `False`.
+        projected_analytic (bool): set to `True` if you want to
+            compute the 2D projected profile analytically (and not
+            through FFTLog). Default: `False`.
+        cumul2d_analytic (bool): set to `True` if you want to
+            compute the 2D cumulative surface density analytically
+            (and not through FFTLog). Default: `False`.
         truncated (bool): set to `True` if the profile should be
             truncated at :math:`r = R_\\Delta` (i.e. zero at larger
             radii.
     """
     name = 'Hernquist'
 
-    def __init__(self, c_M_relation, truncated=True):
+    def __init__(self, c_M_relation,
+                 truncated=True,
+                 fourier_analytic=False,
+                 projected_analytic=False,
+                 cumul2d_analytic=False):
         if not isinstance(c_M_relation, Concentration):
             raise TypeError("c_M_relation must be of type `Concentration`)")
 
         self.cM = c_M_relation
         self.truncated = truncated
+        if fourier_analytic:
+            self._fourier = self._fourier_analytic
+        if projected_analytic:
+            if truncated:
+                raise ValueError("Analytic projected profile not supported "
+                                 "for truncated Hernquist. Set `truncated` or "
+                                 "`projected_analytic` to `False`.")
+            self._projected = self._projected_analytic
+        if cumul2d_analytic:
+            if truncated:
+                raise ValueError("Analytic cumuative 2d profile not supported "
+                                 "for truncated Hernquist. Set `truncated` or "
+                                 "`cumul2d_analytic` to `False`.")
+            self._cumul2d = self._cumul2d_analytic
         super(HaloProfileHernquist, self).__init__()
         self.update_precision_fftlog(padding_hi_fftlog=1E2,
-                                     padding_lo_fftlog=1E-2,
+                                     padding_lo_fftlog=1E-4,
                                      n_per_decade=1000,
                                      plaw_fourier=-2.)
 
@@ -966,6 +993,112 @@ def _real(self, cosmo, r, M, a, mass_def):
             prof = np.squeeze(prof, axis=0)
         return prof
 
+    def _fx_projected(self, x):
+
+        def f1(xx):
+            x2m1 = xx * xx - 1
+            return (-3 / 2 / x2m1**2
+                    + (x2m1+3) * np.arccosh(1 / xx) / 2 / np.fabs(x2m1)**2.5)
+
+        def f2(xx):
+            x2m1 = xx * xx - 1
+            return (-3 / 2 / x2m1**2
+                    + (x2m1+3) * np.arccos(1 / xx) / 2 / np.fabs(x2m1)**2.5)
+
+        xf = x.flatten()
+        return np.piecewise(xf,
+                            [xf < 1, xf > 1],
+                            [f1, f2, 2./15.]).reshape(x.shape)
+
+    def _projected_analytic(self, cosmo, r, M, a, mass_def):
+        r_use = np.atleast_1d(r)
+        M_use = np.atleast_1d(M)
+
+        # Comoving virial radius
+        R_M = mass_def.get_radius(cosmo, M_use, a) / a
+        c_M = self._get_cM(cosmo, M_use, a, mdef=mass_def)
+        R_s = R_M / c_M
+
+        x = r_use[None, :] / R_s[:, None]
+        prof = self._fx_projected(x)
+        norm = 2 * R_s * self._norm(M_use, R_s, c_M)
+        prof = prof[:, :] * norm[:, None]
+
+        if np.ndim(r) == 0:
+            prof = np.squeeze(prof, axis=-1)
+        if np.ndim(M) == 0:
+            prof = np.squeeze(prof, axis=0)
+        return prof
+
+    def _fx_cumul2d(self, x):
+
+        def f1(xx):
+            x2m1 = xx * xx - 1
+            return (1 + 1 / x2m1
+                    + (x2m1 + 1) * np.arccosh(1 / xx) / np.fabs(x2m1)**1.5)
+
+        def f2(xx):
+            x2m1 = xx * xx - 1
+            return (1 + 1 / x2m1
+                    - (x2m1 + 1) * np.arccos(1 / xx) / np.fabs(x2m1)**1.5)
+
+        xf = x.flatten()
+        f = np.piecewise(xf,
+                         [xf < 1, xf > 1],
+                         [f1, f2, 1./3.]).reshape(x.shape)
+
+        return f / x**2
+
+    def _cumul2d_analytic(self, cosmo, r, M, a, mass_def):
+        r_use = np.atleast_1d(r)
+        M_use = np.atleast_1d(M)
+
+        # Comoving virial radius
+        R_M = mass_def.get_radius(cosmo, M_use, a) / a
+        c_M = self._get_cM(cosmo, M_use, a, mdef=mass_def)
+        R_s = R_M / c_M
+
+        x = r_use[None, :] / R_s[:, None]
+        prof = self._fx_cumul2d(x)
+        norm = 2 * R_s * self._norm(M_use, R_s, c_M)
+        prof = prof[:, :] * norm[:, None]
+
+        if np.ndim(r) == 0:
+            prof = np.squeeze(prof, axis=-1)
+        if np.ndim(M) == 0:
+            prof = np.squeeze(prof, axis=0)
+        return prof
+
+    def _fourier_analytic(self, cosmo, k, M, a, mass_def):
+        M_use = np.atleast_1d(M)
+        k_use = np.atleast_1d(k)
+
+        # Comoving virial radius
+        R_M = mass_def.get_radius(cosmo, M_use, a) / a
+        c_M = self._get_cM(cosmo, M_use, a, mdef=mass_def)
+        R_s = R_M / c_M
+
+        x = k_use[None, :] * R_s[:, None]
+        Si2, Ci2 = sici(x)
+        c_Mp1 = c_M + 1
+        P1 = M / ((c_M / c_Mp1)**2 / 2)
+        if self.truncated:
+            Si1, Ci1 = sici(c_Mp1 * x)
+            P2 = x * np.sin(x) * (Ci1 - Ci2) - x * np.cos(x) * (Si1 - Si2)
+            P3 = (-1 + np.sin(c_M[:, None] * x) / (c_Mp1**2 * x)
+                  + c_Mp1 * np.cos(c_M[:, None] * x) / (c_Mp1**2))
+            prof = P1[:, None] * (P2 - P3) / 2
+        else:
+            P2 = (-x * (2 * np.sin(x) * Ci2 + np.pi * np.cos(x))
+                  + 2 * x * np.cos(x) * Si2 + 2) / 4
+            prof = P1[:, None] * P2
+
+        if np.ndim(k) == 0:
+            prof = np.squeeze(prof, axis=-1)
+        if np.ndim(M) == 0:
+            prof = np.squeeze(prof, axis=0)
+        return prof
+
 
 class HaloProfilePressureGNFW(HaloProfile):
     """ Generalized NFW pressure profile by Arnaud et al.

pyccl/tests/test_profiles.py

@@ -61,11 +61,16 @@ def test_empirical_smoke(prof_class):
     with pytest.raises(TypeError):
         p = prof_class(None)
 
-    if prof_class == ccl.halos.HaloProfileNFW:
+    if prof_class in [ccl.halos.HaloProfileNFW,
+                      ccl.halos.HaloProfileHernquist]:
         with pytest.raises(ValueError):
             p = prof_class(c,
                            projected_analytic=True,
                            truncated=True)
+        with pytest.raises(ValueError):
+            p = prof_class(c,
+                           cumul2d_analytic=True,
+                           truncated=True)
 
     p = prof_class(c)
     smoke_assert_prof_real(p)
@@ -328,7 +333,7 @@ def fk(k):
     assert np.all(res < tol)
 
 
-def test_f2r():
+def test_nfw_f2r():
     cM = ccl.halos.ConcentrationDuffy08(M200)
     p1 = ccl.halos.HaloProfileNFW(cM)
     # We force p2 to compute the real-space profile
@@ -388,3 +393,96 @@ def test_nfw_cumul2d_accuracy(fourier_analytic):
 
     res2 = np.fabs(srt2/srt1-1)
     assert np.all(res2 < 5E-3)
+
+
+@pytest.mark.parametrize("use_analytic", [True, False])
+def test_hernquist_accuracy(use_analytic):
+    from scipy.special import sici
+
+    tol = 1E-10 if use_analytic else 5E-3
+    cM = ccl.halos.ConcentrationDuffy08(M200)
+    p = ccl.halos.HaloProfileHernquist(cM, fourier_analytic=use_analytic)
+    M = 1E14
+    a = 0.5
+    c = cM.get_concentration(COSMO, M, a)
+    r_Delta = M200.get_radius(COSMO, M, a) / a
+    r_s = r_Delta / c
+
+    def fk(k):
+        x = k * r_s
+        cp1 = c + 1
+        Si1, Ci1 = sici(cp1 * x)
+        Si2, Ci2 = sici(x)
+        P1 = 1 / ((c / cp1)**2 / 2)
+        P2 = x * np.sin(x) * (Ci1 - Ci2) - x * np.cos(x) * (Si1 - Si2)
+        P3 = (-1 + np.sin(c * x) / (cp1**2 * x)
+              + cp1 * np.cos(c * x) / (cp1**2))
+        return M * P1 * (P2 - P3) / 2
+
+    k_arr = np.logspace(-2, 2, 256) / r_Delta
+    fk_arr = p.fourier(COSMO, k_arr, M, a, M200)
+    fk_arr_pred = fk(k_arr)
+    # Normalize to  1 at low k
+    res = np.fabs((fk_arr - fk_arr_pred) / M)
+    assert np.all(res < tol)
+
+
+def test_hernquist_f2r():
+    cM = ccl.halos.ConcentrationDuffy08(M200)
+    p1 = ccl.halos.HaloProfileHernquist(cM)
+    # We force p2 to compute the real-space profile
+    # by FFT-ing the Fourier-space one.
+    p2 = ccl.halos.HaloProfileHernquist(cM, fourier_analytic=True)
+    p2._real = None
+    p2.update_precision_fftlog(padding_hi_fftlog=1E3)
+
+    M = 1E14
+    a = 0.5
+    r_Delta = M200.get_radius(COSMO, M, a) / a
+
+    r_arr = np.logspace(-2, 1, ) * r_Delta
+    pr_1 = p1.real(COSMO, r_arr, M, a, M200)
+    pr_2 = p2.real(COSMO, r_arr, M, a, M200)
+
+    id_good = r_arr < r_Delta  # Profile is 0 otherwise
+    res = np.fabs(pr_2[id_good] / pr_1[id_good] - 1)
+    assert np.all(res < 0.01)
+
+
+@pytest.mark.parametrize('fourier_analytic', [True, False])
+def test_hernquist_projected_accuracy(fourier_analytic):
+    cM = ccl.halos.ConcentrationDuffy08(M200)
+    # Analytic projected profile
+    p1 = ccl.halos.HaloProfileHernquist(cM, truncated=False,
+                                        projected_analytic=True)
+    # FFTLog
+    p2 = ccl.halos.HaloProfileHernquist(cM, truncated=False,
+                                        fourier_analytic=fourier_analytic)
+
+    M = 1E14
+    a = 0.5
+    rt = np.logspace(-3, 2, 1024)
+    srt1 = p1.projected(COSMO, rt, M, a, M200)
+    srt2 = p2.projected(COSMO, rt, M, a, M200)
+
+    res2 = np.fabs(srt2/srt1-1)
+    assert np.all(res2 < 5E-3)
+
+
+@pytest.mark.parametrize('fourier_analytic', [True, False])
+def test_hernquist_cumul2d_accuracy(fourier_analytic):
+    cM = ccl.halos.ConcentrationDuffy08(M200)
+    # Analytic cumul2d profile
+    p1 = ccl.halos.HaloProfileHernquist(cM, truncated=False,
+                                        cumul2d_analytic=True)
+    # FFTLog
+    p2 = ccl.halos.HaloProfileHernquist(cM, truncated=False,
+                                        fourier_analytic=fourier_analytic)
+    M = 1E14
+    a = 1.0
+    rt = np.logspace(-3, 2, 1024)
+    srt1 = p1.cumul2d(COSMO, rt, M, a, M200)
+    srt2 = p2.cumul2d(COSMO, rt, M, a, M200)
+
+    res2 = np.fabs(srt2/srt1-1)
+    assert np.all(res2 < 5E-3)