add functions for partially delensed spectra: get_lensed_cls_with_spe…

…ctrum and get_partially_lensed_cls; update example notebook with example and correlation function example

camb/correlations.py

@@ -391,6 +391,7 @@ def lensed_cls(cls, clpp, lmax=None, lmax_lensed=None, sampling_factor=1.4, delt
     accurate results should be called with lmax high enough that input cls are effectively band limited
     (lmax >= 2500, or higher for accurate BB to small scales).
     Usually lmax truncation errors are far larger than other numerical errors for lmax<4000.
+    For a faster result use get_lensed_cls_with_spectrum.
 
     :param cls: 2D array of unlensed cls(L,ix), with L starting at zero and ix=0,1,2,3 in order TT, EE, BB, TE.
         cls should include :math:`\ell(\ell+1)/2\pi` factors.

camb/results.py

@@ -423,6 +423,10 @@ def get_cmb_power_spectra(self, params=None, lmax=None,
         are TT, EE, BB, TE, unless raw_cl is True in which case return just :math:`C_\ell`.
         For the lens_potential the power spectrum returned is that of the deflection.
 
+        Note that even if lmax is None, all spectra a returned to the same lmax, appropriate
+        for lensed spectra. Use the individual functions instead if you want to the full unlensed
+        and lensing potential power spectra to the higher lmax actually computed.
+
         :param params: optional :class:`~.model.CAMBparams` instance with parameters to use. If None, must have
           previously set parameters and called `calc_power_spectra` (e.g. if you got this instance
           using :func:`.camb.get_results`),
@@ -1275,6 +1279,46 @@ def get_lensed_gradient_cls(self, lmax=None, CMB_unit=None, raw_cl=False):
         self._scale_cls(res, CMB_unit, raw_cl)
         return res
 
+    def get_lensed_cls_with_spectrum(self, clpp, lmax=None, CMB_unit=None, raw_cl=False):
+        r"""
+        Get lensed CMB power spectra using curved-sky correlation function method, using
+        cpp as the lensing spectrum. Useful for e.g. getting partially-delensed spectra.
+
+        :param clpp: array of :math:`[L(L+1)]^2 C_L^{\phi\phi}/2\pi` lensing potential power spectrum (zero based)
+        :param lmax: lmax to output to
+        :param CMB_unit: scale results from dimensionless. Use 'muK' for :math:`\mu K^2` units for CMB :math:`C_\ell`
+        :param raw_cl: return :math:`C_\ell` rather than :math:`\ell(\ell+1)C_\ell/2\pi`
+        :return: numpy array CL[0:lmax+1,0:4], where 0..3 indexes TT, EE, BB, TE.
+        """
+        assert self.Params.DoLensing
+        lmax_unlens = self.Params.max_l
+        if clpp.shape[0] < lmax_unlens + 1:
+            raise CAMBValueError('clpp must go to at least Params.max_l (zero based)')
+        res = np.empty((lmax_unlens + 1, 4), dtype=np.float64)
+        lmax_lensed = c_int(0)
+        lensClsWithSpectrum = lib_import('lensing', '', 'lensclswithspectrum')
+        lensClsWithSpectrum.argtypes = [POINTER(CAMBdata), numpy_1d, numpy_2d, int_arg]
+        clpp = np.array(clpp, dtype=np.float64)
+        lensClsWithSpectrum(byref(self), clpp, res, byref(lmax_lensed))
+        res = res[:min(lmax_lensed.value, lmax or lmax_lensed.value) + 1, :]
+        self._scale_cls(res, CMB_unit, raw_cl)
+        return res
+
+    def get_partially_lensed_cls(self, Alens: float, lmax=None, CMB_unit=None, raw_cl=False):
+        r"""
+           Get lensed CMB power spectra using curved-sky correlation function method, using
+           true lensing spectrum scaled by Alens.
+           Note that if Params.Alens is also set, the result is scaled by the product of both
+
+           :param Alens: scaling of the lensing relative to true, with Alens=1 being the standard result
+           :param lmax: lmax to output to
+           :param CMB_unit: scale results from dimensionless. Use 'muK' for :math:`\mu K^2` units for CMB :math:`C_\ell`
+           :param raw_cl: return :math:`C_\ell` rather than :math:`\ell(\ell+1)C_\ell/2\pi`
+           :return: numpy array CL[0:lmax+1,0:4], where 0..3 indexes TT, EE, BB, TE.
+           """
+        clpp = self.get_lens_potential_cls()[:, 0] * Alens
+        return self.get_lensed_cls_with_spectrum(clpp, lmax, CMB_unit, raw_cl)
+
     def angular_diameter_distance(self, z):
         """
         Get (non-comoving) angular diameter distance to redshift z.

camb/tests/camb_test.py

@@ -88,9 +88,9 @@ def testAssigments(self):
         self.assertTrue(len(pars.SourceWindows) == 0)
         params = camb.get_valid_numerical_params()
         self.assertEqual(params, {'ombh2', 'deltazrei', 'omnuh2', 'tau', 'omk', 'zrei', 'thetastar', 'nrunrun',
-                                  'meffsterile', 'nnu', 'ntrun', 'HMCode_A_baryon', 'HMCode_eta_baryon', 'HMCode_logT_AGN',
-                                  'cosmomc_theta', 'YHe', 'wa', 'cs2', 'H0', 'mnu', 'Alens', 'TCMB', 'ns',
-                                  'nrun', 'As', 'nt', 'r', 'w', 'omch2'})
+                                  'meffsterile', 'nnu', 'ntrun', 'HMCode_A_baryon', 'HMCode_eta_baryon',
+                                  'HMCode_logT_AGN', 'cosmomc_theta', 'YHe', 'wa', 'cs2', 'H0', 'mnu', 'Alens',
+                                  'TCMB', 'ns', 'nrun', 'As', 'nt', 'r', 'w', 'omch2'})
         params2 = camb.get_valid_numerical_params(dark_energy_model='AxionEffectiveFluid')
         self.assertEqual(params2.difference(params), {'fde_zc', 'w_n', 'zc', 'theta_i'})
 
@@ -352,11 +352,18 @@ def testPowers(self):
         pars.set_for_lmax(lmax)
         cls = data.get_cmb_power_spectra(pars)
         data.get_total_cls(2000)
-        data.get_unlensed_scalar_cls(2500)
+        cls_unlensed = data.get_unlensed_scalar_cls(2500)
         data.get_tensor_cls(2000)
         cls_lensed = data.get_lensed_scalar_cls(3000)
         cphi = data.get_lens_potential_cls(2000)
 
+        cls_lensed2 = data.get_lensed_cls_with_spectrum(cls['lens_potential'][:, 0], lmax=3000)
+        np.testing.assert_allclose(cls_lensed2[2:, :], cls_lensed[2:, :], rtol=1e-4)
+        cls_lensed2 = data.get_partially_lensed_cls(1, lmax=3000)
+        np.testing.assert_allclose(cls_lensed2[2:, :], cls_lensed[2:, :], rtol=1e-4)
+        cls_lensed2 = data.get_partially_lensed_cls(0, lmax=2500)
+        np.testing.assert_allclose(cls_lensed2[2:, :], cls_unlensed[2:, :], rtol=1e-4)
+
         # check lensed CL against python; will only agree well for high lmax as python has no extrapolation template
         cls_lensed2 = correlations.lensed_cls(cls['unlensed_scalar'], cls['lens_potential'][:, 0], delta_cls=False)
         np.testing.assert_allclose(cls_lensed2[2:2000, 2], cls_lensed[2:2000, 2], rtol=1e-3)

fortran/lensing.f90

@@ -65,7 +65,7 @@ module lensing
 
     public lens_Cls, lensing_includes_tensors, lensing_method, lensing_method_flat_corr,&
         lensing_method_curv_corr,lensing_method_harmonic, BessI, ALens_Fiducial, &
-        lensing_sanity_check_amplitude, GetFlatSkyCGrads
+        lensing_sanity_check_amplitude, GetFlatSkyCGrads, lensClsWithSpectrum
     contains
 
 
@@ -83,29 +83,58 @@ subroutine lens_Cls(State)
     end if
     end subroutine lens_Cls
 
-
     subroutine CorrFuncFullSky(State)
     class(CAMBdata) :: State
-    integer :: lmax_extrap
+    integer :: lmax_extrap,l
+    real(dl) CPP(0:State%CP%max_l)
 
     lmax_extrap = State%CP%Max_l - lensed_convolution_margin + 750
     lmax_extrap = min(lmax_extrap_highl,lmax_extrap)
-    call CorrFuncFullSkyImpl(State,State%CP%min_l,max(lmax_extrap,State%CP%max_l))
+    do l= State%CP%min_l,State%CP%max_l
+        ! Cl_scalar(l,1,C_Phi) is l^4 C_phi_phi
+        CPP(l) = State%CLdata%Cl_scalar(l,C_Phi)*(l+1)**2/real(l,dl)**2/const_twopi
+    end do
+    call CorrFuncFullSkyImpl(State, State%ClData, State%ClData, CPP, &
+        State%CP%min_l,max(lmax_extrap,State%CP%max_l))
 
     end subroutine CorrFuncFullSky
 
+    subroutine lensClsWithSpectrum(State, CPP, lensedCls, lmax_lensed)
+    !Get lensed CL using CPP as the lensing specturm
+    !CPP is [L(L+1)]^2C_phi_phi/2/pi
+    type(CAMBdata) :: State
+    real(dl), intent(in) :: CPP(0:State%CP%max_l)
+    real(dl) :: lensedCls(4, 0:State%CP%Max_l)
+    integer :: lmax_lensed
+    Type(TCLData) :: CLout
+    integer :: lmax_extrap, l
+
+    lmax_extrap = State%CP%Max_l - lensed_convolution_margin + 750
+    lmax_extrap = min(lmax_extrap_highl,lmax_extrap)
+    call CorrFuncFullSkyImpl(State, State%ClData, CLout, CPP, &
+        State%CP%min_l,max(lmax_extrap,State%CP%max_l))
+    lmax_lensed = CLout%lmax_lensed
+
+    do l=State%CP%min_l, lmax_lensed
+        lensedCls(:,l) = CLout%Cl_lensed(l,:)
+    end do
+
+    end subroutine lensClsWithSpectrum
+
     subroutine AmplitudeError
 
     call GlobalError('You need to normalize realistically to use lensing. ' &
         //'See http://cosmocoffee.info/viewtopic.php?t=94')
 
     end subroutine AmplitudeError
 
-    subroutine CorrFuncFullSkyImpl(State,lmin,lmax)
+    subroutine CorrFuncFullSkyImpl(State,CL,CLout,CPP,lmin,lmax)
     !Accurate curved sky correlation function method
     !Uses non-perturbative isotropic term with 2nd order expansion in C_{gl,2}
     !Neglects C_{gl}(theta) terms (very good approx)
     class(CAMBdata), target :: State
+    Type(TCLData) :: CL, CLout
+    real(dl) :: CPP(0:State%CP%max_l) ! [L(L+1)]^2 C_L_phi_phi/2pi
     integer, intent(in) :: lmin,lmax
     integer l, i
     integer :: npoints
@@ -138,14 +167,12 @@ subroutine CorrFuncFullSkyImpl(State,lmin,lmax)
     real(dl) range_fac
     logical, parameter :: approx = .false.
     real(dl) theta_cut(lmax), LensAccuracyBoost
-    Type(TCLData), pointer :: CL
     Type(TTimer) :: Timer
 
     !$ integer  OMP_GET_THREAD_NUM, OMP_GET_MAX_THREADS
     !$ external OMP_GET_THREAD_NUM, OMP_GET_MAX_THREADS
 
     if (lensing_includes_tensors) call MpiStop('Haven''t implemented tensor lensing')
-    CL =>  State%ClData
     associate(lSamp => State%CLData%CTransScal%ls, CP=>State%CP)
 
         LensAccuracyBoost = CP%Accuracy%AccuracyBoost*CP%Accuracy%LensingBoost
@@ -154,11 +181,11 @@ subroutine CorrFuncFullSkyImpl(State,lmin,lmax)
             max_lensed_ix = max_lensed_ix -1
         end do
         !150 is the default margin added in python by set_for_lmax
-        CL%lmax_lensed = max(lSamp%l(max_lensed_ix), CP%Max_l - 150)
-
-        if (allocated(CL%Cl_lensed)) deallocate(CL%Cl_lensed)
-        allocate(CL%Cl_lensed(lmin:CL%lmax_lensed,1:4), source = 0._dl)
+        CLout%lmax_lensed = max(lSamp%l(max_lensed_ix), CP%Max_l - 150)
 
+        if (allocated(CLout%Cl_lensed)) deallocate(CLout%Cl_lensed)
+        allocate(CLout%Cl_lensed(lmin:CLout%lmax_lensed,1:4), source = 0._dl)
+
         npoints = CP%Max_l  * 2 *LensAccuracyBoost
         short_integral_range = .not. CP%Accuracy%AccurateBB
         dtheta = const_pi / npoints
@@ -199,15 +226,14 @@ subroutine CorrFuncFullSkyImpl(State,lmin,lmax)
             roots(l) = sqrt(real(l,dl))
         end do
 
-
         thread_ix = 1
         !$ thread_ix = OMP_GET_MAX_THREADS()
-        allocate(lens_contrib(4,CL%lmax_lensed,thread_ix))
+        allocate(lens_contrib(4,CLout%lmax_lensed,thread_ix))
         allocate(ddcontribs(lmax,4),corrcontribs(lmax,4))
 
         do l=lmin,CP%Max_l
             ! (2*l+1)l(l+1)/4pi C_phi_phi: Cl_scalar(l,1,C_Phi) is l^4 C_phi_phi
-            Cphil3(l) = CL%Cl_scalar(l,C_Phi)*(2*l+1)*(l+1)/real(l,dl)**3/const_fourpi
+            Cphil3(l) = CPP(l)*(l+0.5_dl)/real((l+1)*l, dl)
             fac = (2*l+1)/const_fourpi * const_twopi/(l*(l+1))
             CTT(l) =  CL%Cl_scalar(l,C_Temp)*fac
             CEE(l) =  CL%Cl_scalar(l,C_E)*fac
@@ -246,7 +272,7 @@ subroutine CorrFuncFullSkyImpl(State,lmin,lmax)
 
         !$OMP PARALLEL DO DEFAULT(PRIVATE),  &
         !$OMP SHARED(lfacs,lfacs2,lrootfacs,Cphil3,CTT,CTE,CEE,lens_contrib,theta_cut), &
-        !$OMP SHARED(lmin, lmax,dtheta,CL,roots, npoints,interp_fac), &
+        !$OMP SHARED(lmin, lmax,dtheta,CL,CLout,roots, npoints,interp_fac), &
         !$OMP SHARED(jmax,ls,xl,short_integral_range,apodize_point_width)
         do i=1,npoints-1
 
@@ -465,7 +491,7 @@ subroutine CorrFuncFullSkyImpl(State,lmin,lmax)
 
             !$  thread_ix = OMP_GET_THREAD_NUM()+1
 
-            do l=lmin, CL%lmax_lensed
+            do l=lmin, CLout%lmax_lensed
                 !theta factors were put in earlier (already in corr)
 
                 lens_contrib(C_Temp, l, thread_ix)= lens_contrib(C_Temp,l, thread_ix) + &
@@ -488,15 +514,15 @@ subroutine CorrFuncFullSkyImpl(State,lmin,lmax)
         end do
         !$OMP END PARALLEL DO
 
-        do l=lmin, CL%lmax_lensed
+        do l=lmin, CLout%lmax_lensed
             !sign from d(cos theta) = -sin theta dtheta
             fac = l*(l+1)/OutputDenominator*dtheta*const_twopi
-            CL%Cl_lensed(l,CT_Temp) = sum(lens_contrib(CT_Temp,l,:))*fac &
+            CLout%Cl_lensed(l,CT_Temp) = sum(lens_contrib(CT_Temp,l,:))*fac &
                 + CL%Cl_scalar(l,C_Temp)
-            CL%Cl_lensed(l,CT_E) = sum(lens_contrib(CT_E,l,:))*fac &
+            CLout%Cl_lensed(l,CT_E) = sum(lens_contrib(CT_E,l,:))*fac &
                 + CL%Cl_scalar(l,C_E)
-            CL%Cl_lensed(l,CT_B) = sum(lens_contrib(CT_B,l,:))*fac
-            CL%Cl_lensed(l,CT_Cross) = sum(lens_contrib(CT_Cross,l,:))*fac &
+            CLout%Cl_lensed(l,CT_B) = sum(lens_contrib(CT_B,l,:))*fac
+            CLout%Cl_lensed(l,CT_Cross) = sum(lens_contrib(CT_Cross,l,:))*fac &
                 + CL%Cl_scalar(l,C_Cross)
 
         end do