support for calculating time transfer functions and then power spectr…

…a with different initial power/non-linear model parameters

camb/__init__.py

@@ -7,7 +7,7 @@
 __author__ = "Antony Lewis"
 __contact__ = "antony at cosmologist dot info"
 __url__ = "https://camb.readthedocs.io"
-__version__ = "1.0.12.1"
+__version__ = "1.0.12.2"
 
 from . import baseconfig
 

camb/camb.py

@@ -42,17 +42,20 @@ def get_results(params):
     return res
 
 
-def get_transfer_functions(params):
+def get_transfer_functions(params, only_time_sources=False):
     """
     Calculate transfer functions for specified parameters and return :class:`~.results.CAMBdata` instance for
     getting results and subsequently calculating power spectra.
 
     :param params: :class:`.model.CAMBparams` instance
+    :param only_time_sources: does not calculate the CMB l,k transfer functions and does not apply any non-linear
+                              correction scaling. Results with only_time_sources=True can therefore be used with
+                              different initial power spectra to get consistent non-linear lensed spectra.
     :return: :class:`~.results.CAMBdata` instance
     """
 
     res = CAMBdata()
-    res.calc_transfers(params)
+    res.calc_transfers(params, only_transfers=True, only_time_sources=only_time_sources)
     return res
 
 

camb/model.py

@@ -653,7 +653,7 @@ def set_matter_power(self, redshifts=(0.,), kmax=1.2, k_per_logint=None, nonline
 
         :param redshifts: array of redshifts to calculate
         :param kmax: maximum k to calculate (where k is just k, not k/h)
-        :param k_per_logint: number of k steps per log k. Set to zero to use default optimized spacing.
+        :param k_per_logint: minimum number of k steps per log k. Set to zero to use default optimized spacing.
         :param nonlinear: if None, uses existing setting, otherwise boolean for whether to use non-linear matter power.
         :param accurate_massive_neutrino_transfers: if you want the massive neutrino transfers accurately
         :param silent: if True, don't give warnings about sort order
@@ -828,6 +828,18 @@ def get_custom_source_names(self):
     def clear_custom_scalar_sources(self):
         self.f_SetCustomSourcesFunc(byref(c_int(0)), byref(ctypes.c_void_p(0)), np.zeros(0, dtype=np.int32))
 
+    def diff(self, params):
+        """
+        Print differences between this set of parameters and params
+
+        :param params: another CAMBparams instance
+        """
+        p1 = str(params)
+        p2 = str(self)
+        for line1, line2 in zip(p1.split('\n'), p2.split('\n')):
+            if line1 != line2:
+                print(line1, ' <-> ', line2)
+
 
 def set_default_params(P):
     """

camb/results.py

@@ -176,7 +176,8 @@ class CAMBdata(F2003Class):
                  "Number of calculated redshift outputs for the matter transfer (including those for CMB lensing)"),
                 ("transfer_redshifts", AllocatableArrayDouble, "Calculated output redshifts"),
                 ("PK_redshifts_index", c_int * model.max_transfer_redshifts, "Indices of the requested PK_redshifts"),
-                ("OnlyTransfers", c_bool, "Only calculating transfer functions, not power spectra")]
+                ("OnlyTransfers", c_bool, "Only calculating transfer functions, not power spectra"),
+                ("HasScalarTimeSources", c_bool, "calculate and save time source functions, not power spectra")]
 
     # Note there are many more fields in Fortran. Since F2003Class is memory-managed by Fortran, we don't need
     # need to define them all in python.
@@ -281,19 +282,22 @@ def calc_background(self, params):
         if CAMBdata_CalcBackgroundTheory(byref(self), byref(params)):
             config.check_global_error('calc_background')
 
-    def calc_transfers(self, params, only_transfers=True):
+    def calc_transfers(self, params, only_transfers=True, only_time_sources=False):
         """
         Calculate the transfer functions (for CMB and matter power, as determined
         by params.WantCls, params.WantTransfer).
 
         :param params: :class:`~.model.CAMBparams` instance with parameters to use
         :param only_transfers: only calculate transfer functions, no power spectra
+        :param only_time_sources: only calculate time transfer functions, no (p,l,k) transfer functions or non-
+                                  linear scaling
         :return: non-zero if error, zero if OK
         """
         self._check_params(params)
-        if not only_transfers:
+        if not (only_transfers or only_time_sources):
             self._check_powers(params)
-        if CAMBdata_gettransfers(byref(self), byref(params), byref(c_int(1 if only_transfers else 0))):
+        if CAMBdata_gettransfers(byref(self), byref(params), byref(c_int(1 if only_transfers else 0)),
+                                 byref(c_int(1 if only_time_sources else 0))):
             config.check_global_error('calc_transfer')
 
     def _check_powers(self, params=None):
@@ -318,7 +322,7 @@ def calc_power_spectra(self, params=None):
             CAMBdata_transferstopowers(byref(self))
             config.check_global_error()
 
-    def power_spectra_from_transfer(self, initial_power_params, silent=False):
+    def power_spectra_from_transfer(self, initial_power_params=None, silent=False):
         """
         Assuming :meth:`calc_transfers` or :meth:`calc_power_spectra` have already been used, re-calculate the
         power spectra using a new set of initial power spectrum parameters with otherwise the same cosmology.
@@ -328,15 +332,17 @@ def power_spectra_from_transfer(self, initial_power_params, silent=False):
         same results as doing a full recalculation.
 
         :param initial_power_params: :class:`.initialpower.InitialPowerLaw` or :class:`.initialpower.SplinedInitialPower`
-               instance with new primordial power spectrum parameters
+               instance with new primordial power spectrum parameters, or None to use current power spectrum.
         :param silent: suppress warnings about non-linear corrections not being recalculated
         """
-        if not silent and self.Params.NonLinear in [model.NonLinear_lens, model.NonLinear_both] and \
+        if not silent and not self.HasScalarTimeSources and \
+                self.Params.NonLinear in [model.NonLinear_lens, model.NonLinear_both] and \
                 self.Params.WantScalars and self.Params.WantCls and not getattr(self, '_suppress_power_warn', False):
             logging.warning(
                 'power_spectra_from_transfer with non-linear lensing does not recalculate the non-linear correction')
             self._suppress_power_warn = True
-        self.Params.set_initial_power(initial_power_params)
+        if initial_power_params:
+            self.Params.set_initial_power(initial_power_params)
         self._check_powers()
         CAMBdata_transferstopowers(byref(self))
         config.check_global_error()
@@ -1428,7 +1434,7 @@ def cosmomc_theta(self):
 
 CAMBdata_gettransfers = camblib.__handles_MOD_cambdata_gettransfers
 CAMBdata_gettransfers.argtypes = [POINTER(CAMBdata), POINTER(model.CAMBparams),
-                                  POINTER(c_int)]
+                                  POINTER(c_int), POINTER(c_int)]
 CAMBdata_gettransfers.restype = c_int
 
 CAMBdata_transferstopowers = camblib.__camb_MOD_camb_transferstopowers

camb/tests/camb_test.py

@@ -393,6 +393,39 @@ def testPowers(self):
         self.assertAlmostEqual(cls2[1, 0], 1.30388e-10, places=13)
         self.assertAlmostEqual(cls[1, 0], 0)
 
+    def testTimeTransfers(self):
+        from camb import initialpower
+        pars = camb.set_params(H0=69, YHe=0.22, lmax=2000, lens_potential_accuracy=1, ns=0.96, As=2.5e-9)
+        results1 = camb.get_results(pars)
+        cl1 = results1.get_total_cls()
+
+        pars = camb.set_params(H0=69, YHe=0.22, lmax=2000, lens_potential_accuracy=1)
+        results = camb.get_transfer_functions(pars, only_time_sources=True)
+        inflation_params = initialpower.InitialPowerLaw()
+        inflation_params.set_params(ns=0.96, As=2.5e-9)
+        results.power_spectra_from_transfer(inflation_params)
+        cl2 = results.get_total_cls()
+        np.testing.assert_allclose(cl1, cl2, rtol=1e-4)
+        inflation_params.set_params(ns=0.96, As=1.9e-9)
+        results.power_spectra_from_transfer(inflation_params)
+        inflation_params.set_params(ns=0.96, As=2.5e-9)
+        results.power_spectra_from_transfer(inflation_params)
+        cl2 = results.get_total_cls()
+        np.testing.assert_allclose(cl1, cl2, rtol=1e-4)
+
+        pars = camb.CAMBparams()
+        pars.set_cosmology(H0=78, YHe=0.22)
+        pars.set_for_lmax(2000, lens_potential_accuracy=1)
+        pars.WantTensors = True
+        results = camb.get_transfer_functions(pars, only_time_sources=True)
+        cls = []
+        for r in [0, 0.2, 0.4]:
+            inflation_params = initialpower.InitialPowerLaw()
+            inflation_params.set_params(ns=0.96, r=r, nt=0)
+            results.power_spectra_from_transfer(inflation_params)
+            cls += [results.get_total_cls(CMB_unit='muK')]
+        self.assertTrue(np.allclose((cls[1] - cls[0])[2:300, 2] * 2, (cls[2] - cls[0])[2:300, 2], rtol=1e-3))
+
     def testDarkEnergy(self):
         pars = camb.CAMBparams()
         pars.set_cosmology(H0=71)

fortran/camb.f90

@@ -13,9 +13,11 @@ module CAMB
     contains
 
     subroutine CAMB_TransfersToPowers(CData)
+    !From Delta_p_l_k or time transfers to CMB powers and transfers to P(k)
     use CAMBmain
     use lensing
     type (CAMBdata) :: CData
+    logical :: want_tensors, want_vectors
 
     call SetActiveState(CData)
     CData%OnlyTransfer = .false.
@@ -24,6 +26,17 @@ subroutine CAMB_TransfersToPowers(CData)
     if (allocated(Cdata%CAMB_Pk)) deallocate(Cdata%CAMB_PK)
 
     if (CData%CP%WantCls) then
+        if (allocated(CData%ScalarTimeSources) .and. CData%CP%WantScalars) then
+            want_tensors = CData%CP%WantTensors
+            want_vectors = CData%CP%WantVectors
+            Cdata%OnlyTransfer = .true. !prevent ClTransferToCl
+            Cdata%CP%WantTensors = .false.
+            CData%CP%WantVectors = .false.
+            call TimeSourcesToCl
+            Cdata%CP%WantTensors = want_tensors
+            CData%CP%WantVectors = want_vectors
+            Cdata%OnlyTransfer = .false.
+        end if
         call ClTransferToCl(CData)
         if (State%CP%DoLensing .and. global_error_flag==0) call lens_Cls(Cdata)
         if (global_error_flag/=0) return
@@ -33,42 +46,26 @@ subroutine CAMB_TransfersToPowers(CData)
 
     end subroutine CAMB_TransfersToPowers
 
-
     !Call this routine with a set of parameters to generate the results you want.
-    subroutine CAMB_GetResults(OutData, Params, error, onlytransfer)
+    subroutine CAMB_GetResults(OutData, Params, error, onlytransfer, onlytimesources)
     use CAMBmain
     use lensing
     use Bispectrum
     type(CAMBdata)  :: OutData
     type(CAMBparams) :: Params
     integer, optional :: error !Zero if OK
-    logical, optional :: onlytransfer
+    logical, optional :: onlytransfer, onlytimesources
     type(CAMBparams) P
     logical :: call_again
 
     global_error_flag = 0
     call_again = .false.
     call OutData%Free()
     call SetActiveState(OutData)
-    OutData%OnlyTransfer = DefaultFalse(onlytransfer)
-    if (Params%WantCls .and. Params%WantScalars) then
-        P = Params
-        P%Max_eta_k=max(min(P%max_l,3000)*2.5_dl,P%Max_eta_k)
-        P%WantTensors = .false.
-        P%WantVectors = .false.
-        if ((P%NonLinear==NonLinear_lens .or. P%NonLinear==NonLinear_both) .and. &
-            (P%DoLensing .or. State%num_redshiftwindows > 0)) then
-            P%WantTransfer  = .true.
-        end if
-        call OutData%SetParams(P)
-        if (global_error_flag==0) call cmbmain
-        if (global_error_flag/=0) then
-            if (present(error)) error =global_error_flag
-            return
-        end if
-        call_again = .true.
-    end if
+    OutData%HasScalarTimeSources= DefaultFalse(onlytimesources)
+    OutData%OnlyTransfer = DefaultFalse(onlytransfer) .or. OutData%HasScalarTimeSources
 
+    !Vector and tensors first, so at end time steps in state are for scalars
     if (Params%WantCls .and. Params%WantTensors) then
         P=Params
         P%WantTransfer = .false.
@@ -99,6 +96,25 @@ subroutine CAMB_GetResults(OutData, Params, error, onlytransfer)
         call_again = .true.
     end if
 
+    if (Params%WantCls .and. Params%WantScalars) then
+        P = Params
+        P%Max_eta_k=max(min(P%max_l,3000)*2.5_dl,P%Max_eta_k)
+        P%WantTensors = .false.
+        P%WantVectors = .false.
+        if ((P%NonLinear==NonLinear_lens .or. P%NonLinear==NonLinear_both) .and. &
+            (P%DoLensing .or. State%num_redshiftwindows > 0)) then
+            P%WantTransfer  = .true.
+        end if
+        call OutData%SetParams(P)
+        if (global_error_flag==0) call cmbmain
+        if (global_error_flag/=0) then
+            if (present(error)) error =global_error_flag
+            return
+        end if
+        call_again = .true.
+    end if
+
+
     if (Params%WantTransfer .and. .not. (Params%WantCls .and. Params%WantScalars)) then
         P=Params
         P%WantCls = .false.

fortran/camb_python.f90

@@ -243,17 +243,18 @@ subroutine F2003Class_free(cptr,SelfPtr)
 
     end subroutine F2003Class_free
 
-    function CAMBdata_GetTransfers(State, Params, onlytransfer) result(error)
+    function CAMBdata_GetTransfers(State, Params, onlytransfer, onlytimesources) result(error)
     Type (CAMBdata):: State
     type(CAMBparams) :: Params
-    logical :: onlytransfer
+    logical :: onlytransfer, onlytimesources
     integer :: error
 
     error = 0
-    call CAMB_GetResults(State, Params, error, onlytransfer)
+    call CAMB_GetResults(State, Params, error, onlytransfer, onlytimesources)
 
     end function CAMBdata_GetTransfers
 
+
     function CAMBdata_CalcBackgroundTheory(State, P) result(error)
     use cambmain, only: initvars
     Type (CAMBdata):: State