Updates for v0.1 to CrossCorrelation likelihood (#116)

* made folder, moved likelihood

* made folder, moved likelihood and added!

* changed inits

* fix for use of gaussianlikelihood

* restored ccl in init, took .yaml file into test env

* use sacc file in galaxykappalikelihood

* fixed tracer selection

* updated yaml files

* implement spectra (not tracer) selection

* test on two bin selection

* codestyle

* added galaxykappa parms to test

* updated test data file to use

* updates to kappag model

* added test wrt Darwish cmass x act-dr4

* remove comments

* remove debugs

* fix to crosscorrelation test params

* added input checks in crosscorrelation

notebooks/dev/cross_correlation_simulatedata.ipynb

@@ -478,7 +478,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.7"
+   "version": "3.10.8"
   }
  },
  "nbformat": 4,

setup.cfg

@@ -27,7 +27,7 @@ install_requires =
     syslibrary
 
 [options.package_data]
-soliket = *.yaml,*.bibtex,clusters/data/*,clusters/data/selFn_equD56/*,lensing/data/*.txt,lensing/data/*.fits,mflike/*.yaml,tests/*.yaml,data/xcorr_simulated/*.txt,data/CosmoPower/CP_paper/CMB/*.pkl,tests/data/test_bandpass/*,cross_correlation/*.yaml
+soliket = *.yaml,*.bibtex,clusters/data/*,clusters/data/selFn_equD56/*,lensing/data/*.txt,lensing/data/*.fits,mflike/*.yaml,tests/*.yaml,data/xcorr_simulated/*.txt,data/CosmoPower/CP_paper/CMB/*.pkl,tests/data/test_bandpass/*,cross_correlation/*.yaml,tests/data/*.fits
 testpaths = "soliket"
 text_file_format = rst
 

soliket/__init__.py

@@ -4,7 +4,7 @@
 from .ps import PSLikelihood, BinnedPSLikelihood  # noqa: F401
 from .mflike import MFLike  # noqa: F401
 from .mflike import TheoryForge_MFLike
-from .cross_correlation import GalaxyKappaLikelihood, ShearKappaLikelihood  # noqa: F401, E501
+from .cross_correlation import CrossCorrelationLikelihood, GalaxyKappaLikelihood, ShearKappaLikelihood  # noqa: F401, E501
 from .xcorr import XcorrLikelihood  # noqa: F401
 from .foreground import Foreground
 from .bandpass import BandPass
@@ -15,4 +15,4 @@
     from .clusters import ClusterLikelihood  # noqa: F401
 except ImportError:
     print('Skipping cluster likelihood (is pyCCL installed?)')
-    pass
+    pass

soliket/cross_correlation/CrossCorrelationLikelihood.yaml

@@ -1,60 +1,3 @@
-auto_file: soliket/data/xcorr_simulated/clgg_noiseless.txt
-cross_file: soliket/data/xcorr_simulated/clkg_noiseless.txt
-dndz_file: soliket/data/xcorr_simulated/dndz.txt
-
 datapath: soliket/tests/data/des_s-act_kappa.toy-sim.sacc.fits
-use_tracers: all
-
-z_nuisance_mode: null
-m_nuisance_mode: null
-ia_mode: null
-
+use_spectra: all
 ncovsims: null
-
-params:
-  b1: 
-    prior:
-      min: 0.
-      max: 10.
-    latex: b_1
-  s1: 
-    value: 0.4
-    latex: s_1
-  A_IA:
-    value: 0.54
-    latex: A_{IA}
-  gs_des_bin1_A_IA:
-    value: 0.50
-  gs_des_bin2_A_IA:
-    value: 0.52
-  gs_des_bin3_A_IA:
-    value: 0.54
-  gs_des_bin4_A_IA:
-    value: 0.56
-  eta_IA:
-    value: 1.7
-    latex: \eta_{IA}
-  gs_des_bin1_deltaz:
-    value: -0.02
-  gs_des_bin2_deltaz:
-    value: -0.02
-  gs_des_bin3_deltaz:
-    value: -0.02
-  gs_des_bin4_deltaz:
-    value: -0.02
-  gs_des_bin1_m:
-    value: -0.04
-  gs_des_bin2_m:
-    value: -0.04
-  gs_des_bin3_m:
-    value: -0.04
-  gs_des_bin4_m:
-    value: -0.04
-  gs_des_bin1_sigmam:
-    value: 0.01
-  gs_des_bin2_sigmam:
-    value: 0.01
-  gs_des_bin3_sigmam:
-    value: 0.01
-  gs_des_bin4_sigmam:
-    value: 0.01

soliket/cross_correlation/GalaxyKappaLikelihood.yaml

@@ -0,0 +1,9 @@
+params:
+  b1: 
+    prior:
+      min: 0.
+      max: 10.
+    latex: b_1
+  s1: 
+    value: 0.4
+    latex: s_1

soliket/cross_correlation/ShearKappaLikelihood.yaml

@@ -0,0 +1,43 @@
+z_nuisance_mode: null
+m_nuisance_mode: null
+ia_mode: null
+
+params:
+  A_IA:
+    value: 0.54
+    latex: A_{IA}
+  gs_des_bin1_A_IA:
+    value: 0.50
+  gs_des_bin2_A_IA:
+    value: 0.52
+  gs_des_bin3_A_IA:
+    value: 0.54
+  gs_des_bin4_A_IA:
+    value: 0.56
+  eta_IA:
+    value: 1.7
+    latex: \eta_{IA}
+  gs_des_bin1_deltaz:
+    value: -0.02
+  gs_des_bin2_deltaz:
+    value: -0.02
+  gs_des_bin3_deltaz:
+    value: -0.02
+  gs_des_bin4_deltaz:
+    value: -0.02
+  gs_des_bin1_m:
+    value: -0.04
+  gs_des_bin2_m:
+    value: -0.04
+  gs_des_bin3_m:
+    value: -0.04
+  gs_des_bin4_m:
+    value: -0.04
+  gs_des_bin1_sigmam:
+    value: 0.01
+  gs_des_bin2_sigmam:
+    value: 0.01
+  gs_des_bin3_sigmam:
+    value: 0.01
+  gs_des_bin4_sigmam:
+    value: 0.01

soliket/cross_correlation/__init__.py

@@ -1 +1 @@
-from .cross_correlation import CrossCorrelationLikelihood, GalaxyKappaLikelihood, ShearKappaLikelihood   # noqa: F401, E501
+from .cross_correlation import CrossCorrelationLikelihood, GalaxyKappaLikelihood, ShearKappaLikelihood   # noqa: F401, E501

soliket/cross_correlation/cross_correlation.py

@@ -19,18 +19,36 @@ class CrossCorrelationLikelihood(GaussianLikelihood):
     """
     def initialize(self):
 
-        if self.datapath is None:
-            self.dndz = np.loadtxt(self.dndz_file)
-
-            x, y, dy = self._get_data()
-            cov = np.diag(dy**2)
-            self.data = GaussianData("CrossCorrelation", x, y, cov, self.ncovsims)
-        else:
-            self._get_sacc_data()
+        self._get_sacc_data()
+        self._check_tracers()
 
     def get_requirements(self):
         return {"CCL": {"kmax": 10, "nonlinear": True}}
 
+    def _check_tracers(self):
+
+        # check correct tracers
+        for tracer_comb in self.sacc_data.get_tracer_combinations():
+
+            if (self.sacc_data.tracers[tracer_comb[0]].quantity ==
+                    self.sacc_data.tracers[tracer_comb[1]].quantity):
+                raise LoggedError(self.log,
+                                  'You have tried to use {0} to calculate an \
+                                   autocorrelation, but it is a cross-correlation \
+                                   likelihood. Please check your tracer selection in the \
+                                   ini file.'.format(self.__class__.__name__))
+
+            for tracer in tracer_comb:
+                if self.sacc_data.tracers[tracer].quantity not in self._allowable_tracers:
+                    raise LoggedError(self.log,
+                                      'You have tried to use a {0} tracer in \
+                                       {1}, which only allows {2}. Please check your \
+                                       tracer selection in the ini file.\
+                                       '.format(self.sacc_data.tracers[tracer].quantity,
+                                                self.__class__.__name__,
+                                                self._allowable_tracers))
+
+
     def _get_nz(self, z, tracer, tracer_name, **params_values):
 
         if self.z_nuisance_mode == 'deltaz':
@@ -46,11 +64,12 @@ def _get_sacc_data(self, **params_values):
 
         self.sacc_data = sacc.Sacc.load_fits(self.datapath)
 
-        if self.use_tracers == 'all':
+        if self.use_spectra == 'all':
             pass
         else:
-            raise LoggedError('Tracer selection not implemented yet!')
-            # self.sacc_data.keep_selection(tracers=self.use_tracers.split(','))
+            for tracer_comb in self.sacc_data.get_tracer_combinations():
+                if tracer_comb not in self.use_spectra:
+                    self.sacc_data.remove_selection(tracers=tracer_comb)
 
         self.x = self._construct_ell_bins()
         self.y = self.sacc_data.mean
@@ -88,6 +107,16 @@ def _get_data(self, **params_values):
 
         return x, y, dy
 
+    def get_binning(self, tracer_comb):
+
+        bpw_idx = self.sacc_data.indices(tracers=tracer_comb)
+        bpw = self.sacc_data.get_bandpower_windows(bpw_idx)
+        ells_theory = bpw.values
+        ells_theory = np.asarray(ells_theory, dtype=int)
+        w_bins = bpw.weight.T
+
+        return ells_theory, w_bins
+
     def logp(self, **params_values):
         theory = self._get_theory(**params_values)
         return self.data.loglike(theory)
@@ -97,31 +126,51 @@ class GalaxyKappaLikelihood(CrossCorrelationLikelihood):
     r"""
     Likelihood for cross-correlations of galaxy and CMB lensing data.
     """
+    _allowable_tracers = ['cmb_convergence', 'galaxy_density']
+
     def _get_theory(self, **params_values):
+
         cosmo = self.provider.get_CCL()["cosmo"]
 
+        tracer_comb = self.sacc_data.get_tracer_combinations()
+
+        for tracer in np.unique(tracer_comb):
+
+            if self.sacc_data.tracers[tracer].quantity == "cmb_convergence":
+                cmbk_tracer = tracer
+            elif self.sacc_data.tracers[tracer].quantity == "galaxy_density":
+                gal_tracer = tracer
+
+        z_gal_tracer = self.sacc_data.tracers[gal_tracer].z
+        nz_gal_tracer = self.sacc_data.tracers[gal_tracer].nz
+
+        # this should use the bias theory!
         tracer_g = ccl.NumberCountsTracer(cosmo,
                                           has_rsd=False,
-                                          dndz=self.dndz.T,
-                                          bias=(self.dndz[:, 0],
+                                          dndz=(z_gal_tracer, nz_gal_tracer),
+                                          bias=(z_gal_tracer,
                                                 params_values["b1"] *
-                                                    np.ones(len(self.dndz[:, 0]))),
-                                          mag_bias=(self.dndz[:, 0],
+                                                    np.ones(len(z_gal_tracer))),
+                                          mag_bias=(z_gal_tracer,
                                                     params_values["s1"] *
-                                                        np.ones(len(self.dndz[:, 0])))
+                                                        np.ones(len(z_gal_tracer)))
                                           )
         tracer_k = ccl.CMBLensingTracer(cosmo, z_source=1060)
 
-        cl_gg = ccl.cells.angular_cl(cosmo, tracer_g, tracer_g, self.ell_auto)  # + 1e-7
-        cl_kg = ccl.cells.angular_cl(cosmo, tracer_k, tracer_g, self.ell_cross)
+        ells_theory_gk, w_bins_gk = self.get_binning((gal_tracer, cmbk_tracer))
+
+        cl_gk_unbinned = ccl.cells.angular_cl(cosmo, tracer_k, tracer_g, ells_theory_gk)
 
-        return np.concatenate([cl_gg, cl_kg])
+        cl_gk_binned = np.dot(w_bins_gk, cl_gk_unbinned)
+
+        return cl_gk_binned
 
 
 class ShearKappaLikelihood(CrossCorrelationLikelihood):
     r"""
     Likelihood for cross-correlations of galaxy weak lensing shear and CMB lensing data.
     """
+    _allowable_tracers = ["cmb_convergence", "galaxy_shear"]
 
     def _get_theory(self, **params_values):
 
@@ -131,14 +180,6 @@ def _get_theory(self, **params_values):
 
         for tracer_comb in self.sacc_data.get_tracer_combinations():
 
-            if (self.sacc_data.tracers[tracer_comb[0]].quantity
-                    == self.sacc_data.tracers[tracer_comb[1]].quantity):
-                self.log.warning('You requested auto-correlation in '\
-                                 'ShearKappaLikelihood but it is only implemented for '\
-                                 'cross-correlations and resulting bias calculations '\
-                                 'will be incorrect. Please check your tracer '\
-                                 'combinations in the sacc file.')
-
             if self.sacc_data.tracers[tracer_comb[0]].quantity == "cmb_convergence":
                 tracer1 = ccl.CMBLensingTracer(cosmo, z_source=1060)