Merge here the split branch

CHANGELOG.md

@@ -32,3 +32,4 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 1.1.4 : Coverage reports with Codecov
 1.2.0 : Big update of configuration, specification yamls and nuisance parameter interface. No backwards compatibility of yamls!
 1.2.1 : Updating configs of other surveys: SKAO, DESI, LSST to work in new config file structure
+1.2.2 : Galaxy sample split for sources and lenses. Feedback prints more consistent.

cosmicfishpie/LSSsurvey/photo_cov.py

@@ -70,15 +70,24 @@ def __init__(self, cosmopars, photopars, IApars, biaspars, fiducial_Cls=None):
         self.allparsfid.update(self.photopars)
         self.fiducial_Cls = fiducial_Cls
         self.observables = []
+        self.binrange = {}
         for key in cfg.obs:
             if key in ["GCph", "WL"]:
                 self.observables.append(key)
-        self.binrange = cfg.specs["binrange"]
+                if key == "GCph":
+                    self.binrange[key] = cfg.specs["binrange_GCph"]
+                elif key == "WL":
+                    self.binrange[key] = cfg.specs["binrange_WL"]
+
+        self.binrange_WL = cfg.specs["binrange_WL"]
+        self.binrange_GCph = cfg.specs["binrange_GCph"]
         self.feed_lvl = cfg.settings["feedback"]
         self.fsky_WL = cfg.specs.get("fsky_WL")
         self.fsky_GCph = cfg.specs.get("fsky_GCph")
-        self.ngalbin = np.array(cfg.specs["ngalbin"])
-        self.numbins = len(cfg.specs["z_bins_ph"]) - 1
+        self.ngalbin_WL = np.array(cfg.specs["ngalbin_WL"])
+        self.ngalbin_GCph = np.array(cfg.specs["ngalbin_GCph"])
+        self.numbins_WL = len(cfg.specs["z_bins_WL"]) - 1
+        self.numbins_GCph = len(cfg.specs["z_bins_GCph"]) - 1
         self.ellipt_error = cfg.specs["ellipt_error"]
 
     def getcls(self, allpars):
@@ -134,16 +143,17 @@ def getclsnoise(self, cls):
         """
         noisy_cls = copy.deepcopy(cls)
 
-        for ind in self.binrange:
-            for obs in self.observables:
-                if obs == "GCph":
+        for obs in self.observables:
+            if obs == "GCph":
+                for ind in self.binrange_GCph:
                     noisy_cls[obs + " " + str(ind) + "x" + obs + " " + str(ind)] += (
-                        1.0 / self.ngalbin[ind - 1]
+                        1.0 / self.ngalbin_GCph[ind - 1]
                     )
-                elif obs == "WL":
+            elif obs == "WL":
+                for ind in self.binrange_WL:
                     noisy_cls[obs + " " + str(ind) + "x" + obs + " " + str(ind)] += (
                         self.ellipt_error**2.0
-                    ) / self.ngalbin[ind - 1]
+                    ) / self.ngalbin_WL[ind - 1]
 
         return noisy_cls
 
@@ -169,21 +179,26 @@ def get_covmat(self, noisy_cls):
         # Create indexes for data frame
         cols = []
         for o in self.observables:
-            for ind in range(self.numbins):
-                cols.append(o + " " + str(ind + 1))
+            if o == "WL":
+                for ind in range(self.numbins_WL):
+                    cols.append(o + " " + str(ind + 1))
+            elif o == "GCph":
+                for ind in range(self.numbins_GCph):
+                    cols.append(o + " " + str(ind + 1))
 
         for ind, ell in enumerate(noisy_cls["ells"]):
             covdf = pd.DataFrame(index=cols, columns=cols)
             covdf = covdf.fillna(0.0)
 
-            for obs1, obs2, bin1, bin2 in product(
-                self.observables, self.observables, self.binrange, self.binrange
-            ):
-                covdf.at[obs1 + " " + str(bin1), obs2 + " " + str(bin2)] = noisy_cls[
-                    obs1 + " " + str(bin1) + "x" + obs2 + " " + str(bin2)
-                ][ind] / np.sqrt(
-                    np.sqrt(getattr(self, "fsky_" + obs1) * getattr(self, "fsky_" + obs1))
-                )
+            for obs1, obs2 in product(self.observables, self.observables):
+                for bin1, bin2 in product(
+                    self.binrange[obs1], self.binrange[obs2]
+                ):  # MMmod: BEWARE!!! THIS IS VERY UGLY!
+                    covdf.at[obs1 + " " + str(bin1), obs2 + " " + str(bin2)] = noisy_cls[
+                        obs1 + " " + str(bin1) + "x" + obs2 + " " + str(bin2)
+                    ][ind] / np.sqrt(
+                        np.sqrt(getattr(self, "fsky_" + obs1) * getattr(self, "fsky_" + obs1))
+                    )
 
             covvec.append(covdf)
 

cosmicfishpie/LSSsurvey/photo_obs.py

@@ -55,7 +55,7 @@ def memo_integral_efficiency(i, ngal_func, comoving_func, z, zint_mat, diffz):
     intg_mat = np.array(
         [
             (
-                ngal_func(zint_mat[zii], i)
+                ngal_func(zint_mat[zii], i, "WL")
                 * (1 - comoving_func(zint_mat[zii, 0]) / comoving_func(zint_mat[zii]))
             )
             for zii in range(len(zint_mat))
@@ -85,9 +85,8 @@ 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
     """
-    wintgd = ngal_func(zarr, i)[:, None] * (
-        1.0 - comoving_func(zarr)[None, :] / comoving_func(zarr)[:, None]
-    )
+    zprime = zarr[:, None]
+    wintgd = ngal_func(zprime, i, "WL") * (1.0 - comoving_func(zarr) / comoving_func(zprime))
     witri = np.tril(wintgd)
     wint = integrate.trapezoid(witri, zarr, axis=0)
     intp = interp1d(zarr, wint, kind="cubic")
@@ -176,9 +175,17 @@ def __init__(
             time_fin=tcosmo2,
         )
         self.observables = []
+        self.binrange = {}
         for key in cfg.obs:
             if key in ["GCph", "WL"]:
                 self.observables.append(key)
+                if key == "GCph":
+                    self.binrange[key] = cfg.specs["binrange_GCph"]
+                elif key == "WL":
+                    self.binrange[key] = cfg.specs["binrange_WL"]
+
+        self.binrange_WL = cfg.specs["binrange_WL"]
+        self.binrange_GCph = cfg.specs["binrange_GCph"]
 
         tnuis1 = time()
         self.biaspars = biaspars
@@ -203,7 +210,7 @@ def __init__(
         tngal2 = time()
         upt.time_print(
             feedback_level=self.feed_lvl,
-            min_level=2,
+            min_level=3,
             text="---> Galaxy Photometric distributions obtained in ",
             instance=self,
             time_ini=tngal1,
@@ -223,7 +230,7 @@ def __init__(
             raise ValueError("Observables not specified correctly")
 
         self.tracer = cfg.settings["GCph_Tracer"]
-        self.binrange = cfg.specs["binrange"]
+        #        self.binrange = cfg.specs["binrange"]
         self.zsamp = int(round(200 * cfg.settings["accuracy"]))
         if cfg.settings["ell_sampling"] == "accuracy":
             self.ellsamp = int(round(100 * cfg.settings["accuracy"]))
@@ -235,8 +242,9 @@ def __init__(
         self.ell = np.logspace(
             np.log10(cfg.specs["ellmin"]), np.log10(cfg.specs["ellmax"] + 10), num=self.ellsamp
         )
-        self.z_min = cfg.specs["z_bins_ph"][0]
-        self.z_max = cfg.specs["z_bins_ph"][-1]
+
+        self.z_min = np.min([cfg.specs["z_bins_GCph"][0], cfg.specs["z_bins_WL"][0]])
+        self.z_max = np.max([cfg.specs["z_bins_GCph"][-1], cfg.specs["z_bins_WL"][-1]])
         self.z = np.linspace(self.z_min, self.z_max, self.zsamp)
         self.dz = np.mean(np.diff(self.z))
 
@@ -444,7 +452,7 @@ def galaxy_kernel(self, z, i):
         # Wgc = self.window.norm_ngal_photoz(z,i) * np.array([self.nuisance.bias(self.biaspars, i)(z) * \
         # Wgc = self.window.norm_ngal_photoz(z,i) *
         # np.array([self.biaspars['b{:d}'.format(i)] * \
-        Wgc = self.window.norm_ngal_photoz(z, i) * self.cosmo.Hubble(z)
+        Wgc = self.window.norm_ngal_photoz(z, i, "GCph") * self.cosmo.Hubble(z)
         # Wgc = self.window.norm_ngal_photoz(z,i) * self.nuisance.bias(self.biaspars, i)(z) * \
         #                                          self.cosmo.Hubble(z)
 
@@ -493,7 +501,7 @@ def lensing_kernel(self, z, i):
         )
 
         # Adding Intrinsic alignment
-        WIA = self.window.norm_ngal_photoz(z, i) * np.array(
+        WIA = self.window.norm_ngal_photoz(z, i, "WL") * np.array(
             [self.IAvalue(zi) * self.cosmo.Hubble(zi) for zi in z]
         )
         # Sakr Fix June 2023
@@ -546,7 +554,7 @@ def lensing_efficiency(self):
             list of callable functions that give the lensing efficiency for each bin
         """
         teffstart = time()
-        efficiency = [self.integral_efficiency(i) for i in self.binrange]
+        efficiency = [self.integral_efficiency(i) for i in self.binrange_WL]
         efficiency.insert(0, None)
         teffend = time()
         upt.time_print(
@@ -569,7 +577,7 @@ def compute_kernels(self):
         if "GCph" in self.observables:
             self.GC = [
                 interp1d(self.z, self.galaxy_kernel(self.z, ind), kind="cubic")
-                for ind in self.binrange
+                for ind in self.binrange_GCph
             ]
             self.GC.insert(0, None)
         if "WL" in self.observables:
@@ -578,13 +586,13 @@ def compute_kernels(self):
             # self.WL         = [interp1d(self.z,self.lensing_kernel(self.z,ind), kind='cubic') for ind in self.binrange]
             self.WL = [
                 interp1d(self.z, self.lensing_kernel(self.z, ind)[0], kind="cubic")
-                for ind in self.binrange
+                for ind in self.binrange_WL
             ]
             self.WL.insert(0, None)
             # Sakr Fix June 2023
             self.WL_IA = [
                 interp1d(self.z, self.lensing_kernel(self.z, ind)[1], kind="cubic")
-                for ind in self.binrange
+                for ind in self.binrange_WL
             ]
             self.WL_IA.insert(0, None)
         return None
@@ -655,19 +663,21 @@ def computecls(self):
         tcell = time()
 
         # PYTHONIZE THIS HORRIBLE THING
-        for obs1, obs2, bin1, bin2 in product(
-            self.observables, self.observables, self.binrange, self.binrange
-        ):
-            clinterp = self.clsintegral(obs1, obs2, bin1, bin2, hub)
-
-            finalcls = np.zeros((len(full_ell)))
-            for ind, lval in enumerate(full_ell):
-                if (cfg.specs["lmin_" + obs1] <= lval <= cfg.specs["lmax_" + obs1]) and (
-                    cfg.specs["lmin_" + obs2] <= lval <= cfg.specs["lmax_" + obs2]
-                ):
-                    finalcls[ind] = clinterp(lval)
-
-            cls[obs1 + " " + str(bin1) + "x" + obs2 + " " + str(bin2)] = finalcls
+        # for obs1, obs2, bin1, bin2 in product(
+        #        self.observables, self.observables, self.binrange[0], self.binrange[1] #MMmod: BEWARE! THIS IS UGLY!
+        # ):
+        for obs1, obs2 in product(self.observables, self.observables):
+            for bin1, bin2 in product(self.binrange[obs1], self.binrange[obs2]):
+                clinterp = self.clsintegral(obs1, obs2, bin1, bin2, hub)
+
+                finalcls = np.zeros((len(full_ell)))
+                for ind, lval in enumerate(full_ell):
+                    if (cfg.specs["lmin_" + obs1] <= lval <= cfg.specs["lmax_" + obs1]) and (
+                        cfg.specs["lmin_" + obs2] <= lval <= cfg.specs["lmax_" + obs2]
+                    ):
+                        finalcls[ind] = clinterp(lval)
+
+                cls[obs1 + " " + str(bin1) + "x" + obs2 + " " + str(bin2)] = finalcls
 
             tbin = time()
             upt.time_print(
@@ -723,10 +733,13 @@ def computecls_vectorized(self):
 
         # Pre-compute all clsintegrals
         clinterps = {}
-        for obs1, obs2, bin1, bin2 in product(
-            self.observables, self.observables, self.binrange, self.binrange
-        ):
-            clinterps[(obs1, obs2, bin1, bin2)] = self.clsintegral(obs1, obs2, bin1, bin2, hub)
+        for obs1, obs2 in product(self.observables, self.observables):
+            # Get the correct binrange for each observable
+            bins1 = self.binrange[obs1]
+            bins2 = self.binrange[obs2]
+
+            for bin1, bin2 in product(bins1, bins2):
+                clinterps[(obs1, obs2, bin1, bin2)] = self.clsintegral(obs1, obs2, bin1, bin2, hub)
 
         # Vectorized computation of finalcls
         for (obs1, obs2, bin1, bin2), clinterp in clinterps.items():

cosmicfishpie/LSSsurvey/photo_window.py

@@ -44,15 +44,17 @@ def __init__(self, photopars):
         n_i_vec : callable
                   callable function that receives the index of a redshift bin and a numpy.ndarray of redshifts and gives back the binned galaxy redshift distribution without photometric redshift errors
         """
-        self.z_bins = cfg.specs["z_bins_ph"]
-        self.n_bins = len(self.z_bins)
+        self.z_bins_WL = cfg.specs["z_bins_WL"]
+        self.n_bins_WL = len(self.z_bins_WL)
+        self.z_bins_GCph = cfg.specs["z_bins_GCph"]
+        self.n_bins_GCph = len(self.z_bins_GCph)
         self.z0 = cfg.specs["z0"]
         self.z0_p = cfg.specs["z0_p"]
         self.ngamma = cfg.specs["ngamma"]
         self.photo = photopars
-        self.z_min = self.z_bins[0]
-        self.z_max = self.z_bins[-1]
-        self.normalization = self.norm()
+        self.z_min = np.min([cfg.specs["z_bins_GCph"][0], cfg.specs["z_bins_WL"][0]])
+        self.z_max = np.max([cfg.specs["z_bins_GCph"][-1], cfg.specs["z_bins_WL"][-1]])
+        self.normalization = {"GCph": self.norm("GCph"), "WL": self.norm("WL")}
         self.n_i_vec = np.vectorize(self.n_i)
 
     def dNdz(self, z):
@@ -105,7 +107,7 @@ def n_i(self, z, i):
         dNdz_at_z[~mask] = 0.0
         return dNdz_at_z
 
-    def ngal_photoz(self, z, i):
+    def ngal_photoz(self, z, i, obs):
         """Function to compute the binned galaxy redshift distribution convolved with photometric redshift errors n^{ph}_i(z)
 
         Parameters
@@ -128,38 +130,43 @@ def ngal_photoz(self, z, i):
             p_{ph}(z_p|z) = \\frac{1-f_{out}}{\\sqrt{2\\pi}\\sigma_b(1+z)} \\exp\\left\\{-\\frac{1}{2}\\left[\\frac{z-c_bz_p-z_b}{\\sigma_b(1+z)}\\right]^2\\right\\} \\ + \\frac{f_{out}}{\\sqrt{2\\pi}\\sigma_0(1+z)} \\exp\\left\\{-\\frac{1}{2}\\left[\\frac{z-c_0z_p-z_0}{\\sigma_0(1+z)}\\right]^2\\right\\}
         """
 
-        if i == 0 or i >= 11:
-            return None
+        if obs == "GCph":
+            z_bins = self.z_bins_GCph
+        elif obs == "WL":
+            z_bins = self.z_bins_WL
+
+        # if i == 0 or i >= 11:
+        #    return None
 
         term1 = (
             self.photo["cb"]
             * self.photo["fout"]
             * erf(
-                (np.sqrt(1 / 2) * (z - self.photo["zo"] - self.photo["co"] * self.z_bins[i - 1]))
+                (np.sqrt(1 / 2) * (z - self.photo["zo"] - self.photo["co"] * z_bins[i - 1]))
                 / (self.photo["sigma_o"] * (1 + z))
             )
         )
         term2 = (
             -self.photo["cb"]
             * self.photo["fout"]
             * erf(
-                (np.sqrt(1 / 2) * (z - self.photo["zo"] - self.photo["co"] * self.z_bins[i]))
+                (np.sqrt(1 / 2) * (z - self.photo["zo"] - self.photo["co"] * z_bins[i]))
                 / (self.photo["sigma_o"] * (1 + z))
             )
         )
         term3 = (
             self.photo["co"]
             * (1 - self.photo["fout"])
             * erf(
-                (np.sqrt(1 / 2) * (z - self.photo["zb"] - self.photo["cb"] * self.z_bins[i - 1]))
+                (np.sqrt(1 / 2) * (z - self.photo["zb"] - self.photo["cb"] * z_bins[i - 1]))
                 / (self.photo["sigma_b"] * (1 + z))
             )
         )
         term4 = (
             -self.photo["co"]
             * (1 - self.photo["fout"])
             * erf(
-                (np.sqrt(1 / 2) * (z - self.photo["zb"] - self.photo["cb"] * self.z_bins[i]))
+                (np.sqrt(1 / 2) * (z - self.photo["zb"] - self.photo["cb"] * z_bins[i]))
                 / (self.photo["sigma_b"] * (1 + z))
             )
         )
@@ -170,7 +177,7 @@ def ngal_photoz(self, z, i):
             / (2 * self.photo["co"] * self.photo["cb"])
         )
 
-    def norm(self):
+    def norm(self, obs):
         """n^{ph}_i(z)
 
         Parameters
@@ -187,21 +194,26 @@ def norm(self):
 
         """
 
+        if obs == "GCph":
+            z_bins = self.z_bins_GCph
+        elif obs == "WL":
+            z_bins = self.z_bins_WL
+
         # norm = romberg(self.ngal_photoz, self.z_min, self.z_max, args=(i,))
         # Using this as romberg was giving crazy normalizations for the first 2
         # bins
         zint = np.linspace(0.0, self.z_max, 1000)
         dz = self.z_max / 1000
 
         norm = [
-            trapezoid([self.ngal_photoz(z, i) for z in zint], dx=dz)
-            for i in range(1, len(self.z_bins))
+            trapezoid([self.ngal_photoz(z, i, obs) for z in zint], dx=dz)
+            for i in range(1, len(z_bins))
         ]
         norm.insert(0, None)
 
         return norm
 
-    def norm_ngal_photoz(self, z, i):
+    def norm_ngal_photoz(self, z, i, obs):
         """n^{ph}_i(z)
 
         Parameters
@@ -223,4 +235,4 @@ def norm_ngal_photoz(self, z, i):
         # Using this as romberg was giving crazy normalizations for the first 2
         # bins
 
-        return np.array([self.ngal_photoz(zi, i) for zi in z]) / self.normalization[i]
+        return np.array([self.ngal_photoz(zi, i, obs) for zi in z]) / self.normalization[obs][i]