Implementation of the lensing power spectrum computation with CCL (#166)

* Implementing the lensing power spectrum computation with CCL. Added a flag to choose whether to use CAMB or CCL to compute this power spectrum.

* Adding test to check CCL and CAMB agree to predict the CMB lensing power spectrum

---------

Co-authored-by: Raphaël Kou <raphael.kou@apc.in2p3.fr>

soliket/lensing/lensing.py

@@ -52,6 +52,8 @@ class LensingLikelihood(BinnedPSLikelihood, InstallableLikelihood):
     sim_number = 0
     lmax = 3000
     theory_lmax = 10000
+    # flag about whether CCL should be used to compute the cmb lensing power spectrum
+    pp_ccl = False
 
     fiducial_params = {
         "ombh2": 0.02219218,
@@ -151,15 +153,33 @@ def get_requirements(self):
 
         :return: Dictionary ``Cl`` of lmax for each spectrum type.
         """
-        return {
-            "Cl": {
-                "pp": self.theory_lmax,
-                "tt": self.theory_lmax,
-                "te": self.theory_lmax,
-                "ee": self.theory_lmax,
-                "bb": self.theory_lmax,
+        if self.pp_ccl == False:
+            return {
+                "Cl": {
+                    "pp": self.theory_lmax,
+                    "tt": self.theory_lmax,
+                    "te": self.theory_lmax,
+                    "ee": self.theory_lmax,
+                    "bb": self.theory_lmax,
+                }
             }
-        }
+        else:
+            return {
+                "Cl": {
+                    "pp": self.theory_lmax,
+                    "tt": self.theory_lmax,
+                    "te": self.theory_lmax,
+                    "ee": self.theory_lmax,
+                    "bb": self.theory_lmax,
+                },
+                "CCL": {"kmax": 10,
+                        "nonlinear": True},
+                "zstar": None
+            }
+
+    def _get_CCL_results(self):
+        cosmo_dict = self.provider.get_CCL()
+        return cosmo_dict["ccl"], cosmo_dict["cosmo"]
 
     def _get_data(self):
         bin_centers, bandpowers, cov = \
@@ -193,15 +213,21 @@ def _get_theory(self, **params_values):
         """
         cl = self.provider.get_Cl(ell_factor=False)
 
-        Cl_theo = cl["pp"][0: self.lmax]
+        if self.pp_ccl == False:
+            Cl_theo = cl["pp"][0: self.lmax]
+            ls = self.ls
+            Clkk_theo = (ls * (ls + 1)) ** 2 * Cl_theo * 0.25
+        else:
+            ccl, cosmo = self._get_CCL_results()
+            zstar = self.provider.get_param("zstar")
+            cmbk = ccl.CMBLensingTracer(cosmo, z_source=zstar)
+            Clkk_theo = ccl.angular_cl(cosmo, cmbk, cmbk, self.ls)
+
         Cl_tt = cl["tt"][0: self.lmax]
         Cl_ee = cl["ee"][0: self.lmax]
         Cl_te = cl["te"][0: self.lmax]
         Cl_bb = cl["bb"][0: self.lmax]
 
-        ls = self.ls
-        Clkk_theo = (ls * (ls + 1)) ** 2 * Cl_theo * 0.25
-
         Clkk_binned = self.binning_matrix.dot(Clkk_theo)
 
         correction = (

soliket/tests/test_lensing.py

@@ -41,3 +41,35 @@ def test_lensing_like(request):
     loglikes, derived = model.loglikes()
 
     assert np.isclose(loglikes[0], 335.8560097798468, atol=0.2, rtol=0.0)
+
+
+def test_lensing_ccl_limber(request):
+    """
+    Test whether the CMB lensing power spectrum predicted by CCL is the same as with CAMB
+    """
+
+    from cobaya.install import install
+    install({"likelihood": {"soliket.lensing.LensingLikelihood": None}},
+            path=packages_path, skip_global=False, force=True, debug=True)
+
+    from soliket.lensing import LensingLikelihood
+    from copy import deepcopy
+
+    info_dict = deepcopy(info)
+    # Neutrino mass put to 0 as far as it is not included in the ccl wrapper
+    info_dict['params']["mnu"] = 0
+    info_dict['params']["omnuh2"] = 0
+    info_dict['likelihood'] = {"LensingLikelihood": {"external": LensingLikelihood}}
+    model = get_model(info_dict)
+    model.loglikes({})
+    cl_camb = model.likelihood['LensingLikelihood']._get_theory()
+
+    info_dict["likelihood"] = {"LensingLikelihood": {"external": LensingLikelihood, 
+                                                     "pp_ccl": True}}
+    info_dict["theory"]["soliket.CCL"] = {"kmax": 10, "nonlinear": True}
+    model = get_model(info_dict)
+    model.loglikes({})
+    cl_ccl = model.likelihood['LensingLikelihood']._get_theory()
+
+    assert np.any(np.not_equal(cl_ccl, cl_camb))
+    assert np.allclose(cl_ccl, cl_camb, rtol=1e-2, atol=0)