Symbolic Pofk implementation as boltzmann solver

CHANGELOG.md

@@ -24,3 +24,4 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 1.0.6 : Restructured default parameters for EBS and Survey Specifications
 
 1.0.7 : Added a new installation script to download external data from the github
+1.1.0 : Implemented symbolic as a new boltzmann solver class. Fixes of feedback prints. More docstrings.

boltzmann_yaml_files/camb/default.yaml

@@ -11,7 +11,8 @@ COSMO_SETTINGS:
 # Massive neutrino baseline model
   'kmax' : 50
   'share_delta_neff' : True
-  'halofit_version' : 'takahashi'
+  #'halofit_version' : 'takahashi'
+  'halofit_version' : 'mead2020'
   'dark_energy_model' : 'ppf'
   'num_nu_massive' : 1
   'TCMB': 2.7255

cosmicfishpie/LSSsurvey/photo_obs.py

@@ -85,8 +85,9 @@ def faster_integral_efficiency(i, ngal_func, comoving_func, zarr):
     callable
         callable function that receives a numpy.ndarray of requested redshifts and returns the lensing efficiency for the i-th bin as a numpy.ndarray
     """
-    zprime = zarr[:, None]
-    wintgd = ngal_func(zprime, i) * (1.0 - comoving_func(zarr) / comoving_func(zprime))
+    wintgd = ngal_func(zarr, i)[:, None] * (
+        1.0 - comoving_func(zarr)[None, :] / comoving_func(zarr)[:, None]
+    )
     witri = np.tril(wintgd)
     wint = integrate.trapezoid(witri, zarr, axis=0)
     intp = interp1d(zarr, wint, kind="cubic")
@@ -154,7 +155,7 @@ def __init__(
         self.feed_lvl = cfg.settings["feedback"]
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=0,
+            min_level=2,
             text="-> Started Cls calculation",
             instance=self,
         )
@@ -168,7 +169,7 @@ def __init__(
         tcosmo2 = time()
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=1,
+            min_level=2,
             text="---> Cosmological functions obtained in ",
             instance=self,
             time_ini=tcosmo1,
@@ -189,7 +190,7 @@ def __init__(
         tnuis2 = time()
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=2,
+            min_level=3,
             text="---> Nuisance functions obtained in ",
             instance=self,
             time_ini=tnuis1,
@@ -228,19 +229,8 @@ def __init__(
         self.ell = np.logspace(
             np.log10(cfg.specs["ellmin"]), np.log10(cfg.specs["ellmax"] + 10), num=self.ellsamp
         )
-        # self.ell      = np.linspace(cfg.specs['ellmin'],cfg.specs['ellmax']+10,num=self.ellsamp)
-
-        if cfg.input_type == "camb":
-            self.z_min = cfg.specs["z_bins"][0]
-            self.z_max = cfg.specs["z_bins"][-1]  # +0.5
-        if cfg.input_type == "class":
-            self.z_min = cfg.specs["z_bins"][0]
-            self.z_max = cfg.specs["z_bins"][-1]  # +0.5
-        elif cfg.input_type == "external":
-            self.z_min = np.max([cfg.specs["z_bins"][0], self.cosmo.results.zgrid[0]])
-            # self.z_max = np.min([cfg.specs['z_bins'][-1]+0.5, self.cosmo.results.zgrid[-1]])
-            self.z_max = np.min([cfg.specs["z_bins"][-1], self.cosmo.results.zgrid[-1]])
-        # +1 to go beyond the bin limit
+        self.z_min = cfg.specs["z_bins"][0]
+        self.z_max = cfg.specs["z_bins"][-1]
         self.z = np.linspace(self.z_min, self.z_max, self.zsamp)
         self.dz = np.mean(np.diff(self.z))
 
@@ -263,7 +253,7 @@ def compute_all(self):
         tini = time()
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=0,
+            min_level=2,
             text="-> Computing power spectra and kernels ",
             instance=self,
         )
@@ -274,7 +264,7 @@ def compute_all(self):
         tplim2 = time()
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=2,
+            min_level=3,
             text="---> Computed P_limber in ",
             time_ini=tplim1,
             time_fin=tplim2,
@@ -286,7 +276,7 @@ def compute_all(self):
         tkern2 = time()
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=2,
+            min_level=3,
             text="---> Computed Kernels in: ",
             time_ini=tkern1,
             time_fin=tkern2,
@@ -298,7 +288,7 @@ def compute_all(self):
         tcls2 = time()
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=2,
+            min_level=3,
             text="---> Computed Cls in: ",
             time_ini=tcls1,
             time_fin=tcls2,
@@ -308,7 +298,7 @@ def compute_all(self):
 
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=1,
+            min_level=2,
             text="--> Total Cls computation performed in : ",
             time_ini=tini,
             time_fin=tend,
@@ -317,19 +307,20 @@ def compute_all(self):
 
     def print_numerical_specs(self):
         """prints the numerical specifications of the internal computations"""
-        print("***")
-        print("Numerical specifications: ")
-        print("WL ell max = ", str(cfg.specs["lmax_WL"]))
-        print("GCph ell max = ", str(cfg.specs["lmax_GCph"]))
-        print("ell min = ", str(cfg.specs["ellmin"]))
-        print("ell max = ", str(cfg.specs["ellmax"]))
-        print("ell sampling: ", str(self.ellsamp))
-        print("z sampling: ", str(self.zsamp))
-        print("z_min : ", str(self.z_min))
-        print("z_max : ", str(self.z_max))
-        print("z_max : ", str(self.z_max))
-        print("delta_z : ", str(self.dz))
-        print("***")
+        if self.feed_lvl >= 2:
+            print("***")
+            print("Numerical specifications: ")
+            print("WL ell max = ", str(cfg.specs["lmax_WL"]))
+            print("GCph ell max = ", str(cfg.specs["lmax_GCph"]))
+            print("ell min = ", str(cfg.specs["ellmin"]))
+            print("ell max = ", str(cfg.specs["ellmax"]))
+            print("ell sampling: ", str(self.ellsamp))
+            print("z sampling: ", str(self.zsamp))
+            print("z_min : ", str(self.z_min))
+            print("z_max : ", str(self.z_max))
+            print("z_max : ", str(self.z_max))
+            print("delta_z : ", str(self.dz))
+            print("***")
 
     def P_limber(self):
         """
@@ -452,8 +443,14 @@ def galaxy_kernel(self, z, i):
         #                                          self.cosmo.Hubble(z)
 
         tgcend = time()
-        if self.feed_lvl >= 3:
-            print("    ...done bin {} for GCph in {:.2f} s".format(i, tgcend - tgcstart))
+        upt.time_print(
+            feedback_level=self.feed_lvl,
+            min_level=3,
+            text="    ...done bin {} for GCph in: ",
+            instance=self,
+            time_ini=tgcstart,
+            time_fin=tgcend,
+        )
 
         return Wgc
 
@@ -498,8 +495,14 @@ def lensing_kernel(self, z, i):
 
         twlend = time()
 
-        if self.feed_lvl >= 3:
-            print("    ...done bin {} for WL in {:.2f} s".format(i, twlend - twlstart))
+        upt.time_print(
+            feedback_level=self.feed_lvl,
+            min_level=3,
+            text="    ...done bin {} for WL in:",
+            instance=self,
+            time_ini=twlstart,
+            time_fin=twlend,
+        )
 
         # Sakr Fix June 2023
         # return Wwl
@@ -540,8 +543,14 @@ def lensing_efficiency(self):
         efficiency = [self.integral_efficiency(i) for i in self.binrange]
         efficiency.insert(0, None)
         teffend = time()
-        if self.feed_lvl >= 3:
-            print("    ...lensing efficiency computation took {:.2f} s".format(teffend - teffstart))
+        upt.time_print(
+            feedback_level=self.feed_lvl,
+            min_level=3,
+            text="    ...lensing efficiency computation done in:",
+            instance=self,
+            time_ini=teffstart,
+            time_fin=teffend,
+        )
         return efficiency
 
     def compute_kernels(self):
@@ -622,10 +631,12 @@ def computecls(self):
                 P_{\\delta \\delta} \\big[\\frac{\\ell + 1/2}{r(z)} , z \\big]
         """
 
-        if self.feed_lvl > 1:
-            print("")
-        if self.feed_lvl > 1:
-            print("    Computing Cls integral for {}".format(self.observables))
+        upt.time_print(
+            feedback_level=self.feed_lvl,
+            min_level=2,
+            text="Computing Cls integral for {}".format(self.observables),
+            instance=self,
+        )
         # full_ell = np.linspace(cfg.specs['ellmin'],cfg.specs['ellmax'],cfg.specs['ellmax']-cfg.specs['ellmin'])
         # full_ell = np.round(full_ell, 0)
         full_ell = self.ell
@@ -653,21 +664,25 @@ def computecls(self):
             cls[obs1 + " " + str(bin1) + "x" + obs2 + " " + str(bin2)] = finalcls
 
             tbin = time()
-            if self.feed_lvl > 2:
-                print("")
-            if self.feed_lvl > 2:
-                print(
-                    "    ...{} {} x {} {} done in {:.2f} s".format(
-                        obs1, bin1, obs2, bin2, tbin - tcell
-                    )
-                )
+            upt.time_print(
+                feedback_level=self.feed_lvl,
+                min_level=3,
+                text="    ...{} {} x {} {} done in: ".format(obs1, bin1, obs2, bin2),
+                instance=self,
+                time_ini=tcell,
+                time_fin=tbin,
+            )
             tcell = tbin
 
         tend = time()
-        if self.feed_lvl > 1:
-            print("")
-        if self.feed_lvl > 1:
-            print("    Cls integral computation took {:.2f} s".format(tend - tstart))
+        upt.time_print(
+            feedback_level=self.feed_lvl,
+            min_level=3,
+            text="Cls integral computation done in: ",
+            instance=self,
+            time_ini=tstart,
+            time_fin=tend,
+        )
 
         return cls
 

cosmicfishpie/LSSsurvey/spectro_cov.py

@@ -63,6 +63,7 @@ def __init__(self, fiducialpars, fiducial_specobs=None, bias_samples=["g", "g"])
         # initializing the class only with fiducial parameters
         # if fiducial_specobs is None:
 
+        self.feed_lvl = cfg.settings["feedback"]
         try:
             self.fsky_spectro = cfg.specs["fsky_spectro"]
             self.area_survey = self.fsky_spectro * upm.areasky()
@@ -71,7 +72,12 @@ def __init__(self, fiducialpars, fiducial_specobs=None, bias_samples=["g", "g"])
             self.fsky_spectro = self.area_survey / upm.areasky()
         if "IM" in cfg.obs and "GCsp" in cfg.obs:
             bias_samples = ["I", "g"]
-            print("Entering Cov cross XC IM,g term")
+            upt.time_print(
+                feedback_level=self.feed_lvl,
+                min_level=2,
+                text="Entering Cov cross XC IM,g term",
+                instance=self,
+            )
             self.pk_obs = spec_obs.ComputeGalIM(
                 fiducialpars, fiducialpars, bias_samples=bias_samples
             )
@@ -83,13 +89,23 @@ def __init__(self, fiducialpars, fiducial_specobs=None, bias_samples=["g", "g"])
             )
         elif "IM" in cfg.obs and "I" in bias_samples:
             bias_samples = ["I", "I"]
-            print("Entering Cov IM term")
+            upt.time_print(
+                feedback_level=self.feed_lvl,
+                min_level=2,
+                text="Entering Cov IM term",
+                instance=self,
+            )
             self.pk_obs = spec_obs.ComputeGalIM(
                 fiducialpars, fiducialpars, bias_samples=bias_samples
             )
         elif "GCsp" in cfg.obs and "g" in bias_samples:
             bias_samples = ["g", "g"]
-            print("Entering Cov gg term")
+            upt.time_print(
+                feedback_level=self.feed_lvl,
+                min_level=2,
+                text="Entering Cov gg term",
+                instance=self,
+            )
             self.pk_obs = spec_obs.ComputeGalSpectro(
                 fiducialpars, fiducialpars, bias_samples=bias_samples
             )
@@ -406,7 +422,7 @@ def __init__(self, z_array, pk_kmesh, pk_mumesh, fiducial_spectro_obj, bias_samp
         self.pk_kmesh = pk_kmesh
         self.pk_mumesh = pk_mumesh
         self.freeparams = None  # cfg.freeparams
-        self.feed_lvl = 1  # cfg.settings['feedback']
+        self.feed_lvl = cfg.settings["feedback"]
         # self.get_obs = memory.cache(self.getobs)
 
     def initialize_obs(self, allpars):
@@ -547,7 +563,7 @@ def compute_derivs(self, freeparams=dict()):
             tder2 = time()
             upt.time_print(
                 feedback_level=self.feed_lvl,
-                min_level=1,
+                min_level=3,
                 text="-->> Derivatives computed in ",
                 time_ini=tder1,
                 time_fin=tder2,

cosmicfishpie/LSSsurvey/spectro_obs.py

@@ -121,7 +121,7 @@ def __init__(
         self.feed_lvl = cfg.settings["feedback"]
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=1,
+            min_level=2,
             text="Entered ComputeGalSpectro",
             instance=self,
         )
@@ -241,7 +241,7 @@ def __init__(
         tend = time()
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=1,
+            min_level=2,
             text="GalSpec initialization done in: ",
             time_ini=tini,
             time_fin=tend,
@@ -847,7 +847,7 @@ def observed_Pgg(self, z, k, mu, b_i=None):
         tend = time()
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=1,
+            min_level=2,
             text="observed P_gg computation took: ",
             time_ini=tstart,
             time_fin=tend,
@@ -902,7 +902,9 @@ def __init__(
 
         tini = time()
         self.feed_lvl = cfg.settings["feedback"]
-        print("Entered ComputeGalIM")
+        upt.time_print(
+            feedback_level=self.feed_lvl, min_level=2, text="Entered ComputeGalIM", instance=self
+        )
 
         if "IM" not in self.observables:
             raise AttributeError("Observables list not defined properly")
@@ -935,7 +937,7 @@ def __init__(
         tend = time()
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=1,
+            min_level=2,
             text="GalIM initialization done in: ",
             time_ini=tini,
             time_fin=tend,