Merge jelly

CHANGELOG.md

@@ -33,3 +33,4 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 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.
+1.2.3 : Added a new likelihood function for photometric and spectroscopic surveys.

cosmicfishpie/LSSsurvey/spectro_cov.py

@@ -92,14 +92,20 @@ def __init__(
             self.global_z_bins = self.z_bins
             self.inter_z_bin_mids = self.z_bin_mids
             self.inter_z_bins = self.z_bins
+            self.inter_z_bin_mids = self.z_bin_mids
+            self.inter_z_bins = self.z_bins
         if "IM" in self.pk_obs.observables:
+            self.IM_z_bins = self.pk_obs.nuisance.IM_zbins
+            self.IM_z_bin_mids = self.pk_obs.nuisance.IM_zbins_mids
             self.IM_z_bins = self.pk_obs.nuisance.IM_zbins
             self.IM_z_bin_mids = self.pk_obs.nuisance.IM_zbins_mids
             self.Tsys_interp = self.pk_obs.nuisance.IM_THI_noise()
             self.global_z_bin_mids = self.IM_z_bin_mids
             self.global_z_bins = self.IM_z_bins
             self.inter_z_bin_mids = self.IM_z_bin_mids
             self.inter_z_bins = self.IM_z_bins
+            self.inter_z_bin_mids = self.IM_z_bin_mids
+            self.inter_z_bins = self.IM_z_bins
         # Choose longest zbins array to loop in Fisher matrix
         if "GCsp" in self.pk_obs.observables and "IM" in self.pk_obs.observables:
             self.global_z_bin_mids = np.union1d(self.z_bin_mids, self.IM_z_bin_mids)
@@ -140,8 +146,8 @@ def volume_bin(self, zi, zj):
         float, numpy.ndarray
             Volume of the comoving spherical shell between zj and zi
         """
-        d1 = self.pk_obs.cosmo.angdist(zi)
-        d2 = self.pk_obs.cosmo.angdist(zj)
+        d1 = self.pk_obs.fiducialcosmo.angdist(zi)
+        d2 = self.pk_obs.fiducialcosmo.angdist(zj)
         sphere_vol = (4 * np.pi / 3) * (pow((1 + zj) * d2, 3) - pow((1 + zi) * d1, 3))
         return sphere_vol
 
@@ -288,6 +294,8 @@ def veff_II(self, zi, k, mu):
 
         Parameters
         ----------
+        zi : float
+             Redshift
         zi : float
              Redshift
         k  : float, numpy.ndarray
@@ -301,8 +309,12 @@ def veff_II(self, zi, k, mu):
         """
         pobs = self.pk_obs.observed_P_ij(zi, k, mu, si="I", sj="I")
         pnoisy = self.noisy_P_ij(zi, k, mu, si="I", sj="I")
+        pobs = self.pk_obs.observed_P_ij(zi, k, mu, si="I", sj="I")
+        pnoisy = self.noisy_P_ij(zi, k, mu, si="I", sj="I")
         prefactor = 1 / (8 * (np.pi**2))
         covterm = prefactor * (pobs / pnoisy) ** 2
+        if zi < self.inter_z_bin_mids[0] or zi > self.inter_z_bin_mids[-1]:
+            covterm = np.zeros_like(covterm)
         if zi < self.inter_z_bin_mids[0] or zi > self.inter_z_bin_mids[-1]:
             covterm = np.zeros_like(covterm)
         return covterm
@@ -312,6 +324,8 @@ def veff_Ig(self, zi, k, mu):
 
         Parameters
         ----------
+        zi : float
+             Redshift
         zi : float
              Redshift
         k  : float, numpy.ndarray
@@ -330,9 +344,16 @@ def veff_Ig(self, zi, k, mu):
         pnoisy_Ig = self.noisy_P_ij(zi, k, mu, si="I", sj="g")
         pnoisy_II = self.noisy_P_ij(zi, k, mu, si="I", sj="I")
         pnoisy_gg = self.noisy_P_ij(zi, k, mu, si="g", sj="g")
+        # the si, sj indices will be overwritten by the self.bias_samples in the observed_P_Ig function
+        pobs_Ig = self.pk_obs.observed_P_ij(zi, k, mu, si="I", sj="g")
+        pnoisy_Ig = self.noisy_P_ij(zi, k, mu, si="I", sj="g")
+        pnoisy_II = self.noisy_P_ij(zi, k, mu, si="I", sj="I")
+        pnoisy_gg = self.noisy_P_ij(zi, k, mu, si="g", sj="g")
         covterm = pobs_Ig**2 / (pnoisy_gg * pnoisy_II + pnoisy_Ig * pnoisy_Ig)
         prefactor = 1 / (4 * (np.pi**2))
         covterm = prefactor * covterm
+        if zi < self.inter_z_bin_mids[0] or zi > self.inter_z_bin_mids[-1]:
+            covterm = np.zeros_like(covterm)
         if zi < self.inter_z_bin_mids[0] or zi > self.inter_z_bin_mids[-1]:
             covterm = np.zeros_like(covterm)
         return covterm
@@ -504,6 +525,7 @@ def exact_derivs(self, par):
             for ii, zzi in enumerate(self.z_array):
                 pgg_obs = self.pobs.observed_Pgg(zzi, self.pk_kmesh, self.pk_mumesh)
                 z_bin_mids = self.z_array
+                z_bin_mids = self.z_array
                 kron = self.kronecker_bins(par, z_bin_mids, zzi)
                 deriv_i = 1 / pgg_obs
                 deriv[ii] = kron * deriv_i
@@ -681,6 +703,7 @@ def dlnpobs_dnuisp(self, zi, k, mu, par):
             deriv = 0
         # get index in bin
         ii = np.where(np.isclose(self.global_z_bin_mids, zi))
+        ii = np.where(np.isclose(self.global_z_bin_mids, zi))
         ii = ii[0][0] + 1
         pi = par.split("_")
         pi = int(pi[-1])

cosmicfishpie/LSSsurvey/spectro_obs.py

@@ -39,10 +39,16 @@ def __init__(
         ----------
         cosmopars            : dict
                               A dictionary containing the cosmological parameters of the sample cosmology
+                              A dictionary containing the cosmological parameters of the sample cosmology
         fiducial_cosmopars   : dict, optional
                               A dictionary containing the cosmological parameters of the fiducial/reference cosmology
         spectrobiaspars      : dict, optional
                               A dictionary containing the specifications for the galaxy biases
+        IMbiaspars          : dict, optional
+                              A dictionary containing the specifications for the intensity mapping biases
+                              A dictionary containing the cosmological parameters of the fiducial/reference cosmology
+        spectrobiaspars      : dict, optional
+                              A dictionary containing the specifications for the galaxy biases
         IMbiaspars          : dict, optional
                               A dictionary containing the specifications for the intensity mapping biases
         spectrononlinearpars : dict, optional
@@ -64,34 +70,51 @@ def __init__(
         ----------
         feed_lvl                      : int
                                        Number indicating the verbosity of the output. Higher numbers mean more output
+                                       Number indicating the verbosity of the output. Higher numbers mean more output
         observables                   : list
                                        A list of the observables that the observed power spectrum is computed for
+                                       A list of the observables that the observed power spectrum is computed for
         s8terms                       : bool
                                        If True will expand the observed power spectrum with :math:`\\sigma_8` to match the IST:F recipe
+                                       If True will expand the observed power spectrum with :math:`\\sigma_8` to match the IST:F recipe
         tracer                        : str
                                        What Power spectrum should be used when calculating the power spectrum. Either "matter" or "clustering"
+                                       What Power spectrum should be used when calculating the power spectrum. Either "matter" or "clustering"
         cosmo                         : cosmicfishpie.cosmology.cosmology.cosmo_functions
                                        An instance of `cosmo_functions` of the sample cosmology
+                                       An instance of `cosmo_functions` of the sample cosmology
         nuisance                      : cosmicfishpie.cosmology.Nuisance.Nuisance
                                        An instance of `nuisance` that contains the relevant modeling of nuisance parameters
+                                       An instance of `nuisance` that contains the relevant modeling of nuisance parameters
         extraPshot                    : dict
                                        A dictionary containing the values of the additional shot noise per bin
         gcsp_z_bin_mids              : numpy.ndarray
                                        Lists the redshift bin centers for galaxy clustering
+        IM_z_bin_mids                : numpy.ndarray
+                                       Lists the redshift bin centers for intensity mapping
+                                       A dictionary containing the values of the additional shot noise per bin
+        gcsp_z_bin_mids              : numpy.ndarray
+                                       Lists the redshift bin centers for galaxy clustering
         IM_z_bin_mids                : numpy.ndarray
                                        Lists the redshift bin centers for intensity mapping
         k_grid                        : numpy.ndarray
                                        Lists all wavenumbers used in the internal calculations
+                                       Lists all wavenumbers used in the internal calculations
         dk_grid                       : numpy.ndarray
                                        Lists the numerical distance between all wavenumbers used in the internal calculations
+                                       Lists the numerical distance between all wavenumbers used in the internal calculations
         linear_switch                 : bool
                                        If False all nonlinear effects will be included in the computation
+                                       If False all nonlinear effects will be included in the computation
         FoG_switch                    : bool
                                        If True and `linear_switch` is True, then the finger of god effect will be modelled
+                                       If True and `linear_switch` is True, then the finger of god effect will be modelled
         APbool                        : bool
                                        If True and `linear_switch` is True, then the Alcock-Paczynski effect will be considered
+                                       If True and `linear_switch` is True, then the Alcock-Paczynski effect will be considered
         fix_cosmo_nl_terms            : bool
                                        If True and the nonlinear modeling parameters are not varied, then they will be fixed to the fiducial cosmology values
+                                       If True and the nonlinear modeling parameters are not varied, then they will be fixed to the fiducial cosmology values
         dz_err                        : float
                                        Value of the spectroscopic redshift error
         Dd                            : float
@@ -185,17 +208,25 @@ def __init__(
             **self.cosmopars,
             **self.spectrobiaspars,
             **self.IMbiaspars,
+            **self.IMbiaspars,
             **self.PShotpars,
             **self.spectrononlinearpars,
         }
         self.fiducial_allpars = {
             **self.fiducial_cosmopars,
             **self.fiducial_spectrobiaspars,
             **self.fiducial_IMbiaspars,
+            **self.fiducial_IMbiaspars,
             **self.fiducial_PShotpars,
             **self.fiducial_spectrononlinearpars,
         }
 
+        if "IM" in self.observables and "GCsp" in self.observables:
+            self.obs_spectrum = ["I", "g"]
+        elif "IM" in self.observables:  # compute  in the case of IM only
+            self.obs_spectrum = ["I", "I"]
+        elif "GCsp" in self.observables:
+            self.obs_spectrum = ["g", "g"]
         if "IM" in self.observables and "GCsp" in self.observables:
             self.obs_spectrum = ["I", "g"]
         elif "IM" in self.observables:  # compute  in the case of IM only
@@ -503,24 +534,29 @@ def BAO_term(self, z):
         return bao
 
     def bterm_fid(self, z, k=None, bias_sample="g"):
-        """
-        Calculates the fiducial bias term at a given redshift z,
-        and an optional wavenumber k.
-        of either galaxies or intensity mapping.
+        """Calculate the fiducial bias term for galaxies or intensity mapping.
 
-        Parameters:
-        -----------
-            z           : float, numpy.ndarray
-                      The redshifts value at which to evaluate the bias term.
-            k           : float, numpy.ndarray, optional
-                      The wavenumber at which to evaluate the bias term.
-            bias_sample : str, optional
-                      Specifies whether to compute the galaxy ('g') or intensity mapping ('I') bias term. (default='g')
+        Parameters
+        ----------
+        z : float or numpy.ndarray
+            Redshift value(s) at which to evaluate the bias term.
+        k : float or numpy.ndarray, optional
+            Wavenumber(s) at which to evaluate the bias term.
+        bias_sample : str, optional
+            Type of bias to compute. Must be either:
+            - 'g' for galaxy clustering bias (default)
+            - 'I' for intensity mapping bias
 
         Returns
-        --------
-        float
-        The value of the bias term at `z` and `k`, if provided.
+        -------
+        float or numpy.ndarray
+            The computed bias term at the specified z and k values.
+
+        Note
+        -----
+        The total bias is computed as the product of a redshift-dependent term and
+        a scale-dependent term: b_total = b(z) * b(k)
+
         """
         if bias_sample == "g":
             if bias_sample != self.sp_bias_sample: