Work on CCL module for cross-correlations (#16)

* Added CCL CosmologyCalculator functionality

* Added CCL CosmologyCalculator functionality

* Added notebook for testing

* Added ccl likelihood

* small changes

* Dynamically importing ccl

* Changed import back to top level

* Renamed Tester likelihood and added data files

* Change creation of CCL cosmo object

* Change likelihood to use GaussianLikelihood structure

* Removed old comments

* Separate likelihood file and tests

* Moved test file and other small changes

* Removed python path

* Updated example notebook

* Updated assert statement

* Changed max_tries

* Increased tolerance of assert

* fixed bug

soliket/CrossCorrelationLikelihood.yaml

@@ -0,0 +1,13 @@
+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
+
+params:
+  b1: 
+    prior:
+      min: 0.
+      max: 10.
+    latex: b_1
+  s1: 
+    value: 0.4
+    latex: s_1

soliket/__init__.py

@@ -3,4 +3,10 @@
 from .ps import PSLikelihood, BinnedPSLikelihood
 from .clusters import ClusterLikelihood
 from .mflike import MFLike
-from .ccl import CCL
+try:
+    import pyccl as ccl
+    from .ccl import CCL
+    from .cross_correlation import CrossCorrelationLikelihood
+except ImportError:
+    print('Skipping CCL module as pyCCL is not installed')
+    pass

soliket/ccl.py

@@ -26,6 +26,10 @@
 Note that this approach preclude sharing results other than the cosmo object itself between different likelihoods.
 
 Also note lots of things still cannot be done consistently in CCL, so this is far from general.
+
+April 2021:
+-----------
+Second version by PL. Using CCL's newly implemented cosmology calculator. 
 """
 
 # For Cobaya docs see
@@ -35,7 +39,7 @@
 import numpy as np
 from typing import Sequence, Union
 from cobaya.theory import Theory
-
+import pyccl as ccl
 
 class CCL(Theory):
     # Options for Pk.
@@ -45,7 +49,9 @@ class CCL(Theory):
     z: Union[Sequence, np.ndarray] = []  # redshift sampling
     extra_args: dict = {}  # extra (non-parameter) arguments passed to ccl.Cosmology()
 
-    _default_z_sampling = np.linspace(0, 5, 100)
+    _logz = np.linspace(-3, np.log10(1100), 150)
+    _default_z_sampling = 10**_logz
+    _default_z_sampling[0] = 0
 
     def initialize(self):
         self._var_pairs = set()
@@ -54,15 +60,12 @@ def initialize(self):
     def get_requirements(self):
         # These are currently required to construct a CCL cosmology object.
         # Ultimately CCL should depend only on observable not parameters
-        # 'As' could be substituted by sigma8.
-        return {'omch2', 'ombh2', 'ns', 'As'}
+        return {'omch2', 'ombh2'}
 
     def must_provide(self, **requirements):
         # requirements is dictionary of things requested by likelihoods
         # Note this may be called more than once
 
-        # CCL currently has no way to infer the required inputs from the required outputs
-        # So a lot of this is fixed
         if 'CCL' not in requirements:
             return {}
         options = requirements.get('CCL') or {}
@@ -94,6 +97,9 @@ def must_provide(self, **requirements):
         needs['Hubble'] = {'z': self.z}
         needs['comoving_radial_distance'] = {'z': self.z}
 
+        needs['fsigma8'] = {'z': self.z}
+        needs['sigma8_z'] = {'z': self.z}
+
         assert len(self._var_pairs) < 2, "CCL doesn't support other Pk yet"
         return needs
 
@@ -109,69 +115,64 @@ def calculate(self, state, want_derived=True, **params_values_dict):
         distance = self.provider.get_comoving_radial_distance(self.z)
         hubble_z = self.provider.get_Hubble(self.z)
         H0 = hubble_z[0]
+        h = H0/100
         E_of_z = hubble_z / H0
+
+        Omega_c = self.provider.get_param('omch2') / h ** 2
+        Omega_b = self.provider.get_param('ombh2') / h ** 2
         # Array z is sorted in ascending order. CCL requires an ascending scale factor
         # as input
         # Flip the arrays to make them a function of the increasing scale factor.
         # If redshift sampling is changed, check that it is monotonically increasing
         distance = np.flip(distance)
         E_of_z = np.flip(E_of_z)
 
-        # Create a CCL cosmology object
-        import pyccl as ccl
-        h = H0 / 100.
-        Omega_c = self.provider.get_param('omch2') / h ** 2
-        Omega_b = self.provider.get_param('ombh2') / h ** 2
-
-        # Currently, CCL requires the (ill-defined) linear "matter" perturbation
-        # growth factor and rate. Because it's ill-defined, we can't get it from
-        # Boltzmann code in general; ultimately CCL should use more physical
-        # inputs for anything of use in general models.
-        # For now just compute from CCL itself to keep it happy:
-        cosmo = ccl.Cosmology(Omega_c=Omega_c, Omega_b=Omega_b, h=h,
-                              n_s=self.provider.get_param('ns'),
-                              A_s=self.provider.get_param('As'))
         # Array z is sorted in ascending order. CCL requires an ascending scale
         # factor as input
         a = 1. / (1 + self.z[::-1])
-        growth = ccl.background.growth_factor(cosmo, a)
-        fgrowth = ccl.background.growth_rate(cosmo, a)
-        # In order to use CCL with input arrays, the cosmology object needs
-        # to be reset. This should be improved...
-        cosmo = ccl.Cosmology(Omega_c=Omega_c, Omega_b=Omega_b, h=h,
-                              n_s=self.provider.get_param('ns'),
-                              A_s=self.provider.get_param('As'), **self.extra_args)
-        cosmo._set_background_from_arrays(a_array=a,
-                                          chi_array=distance,
-                                          hoh0_array=E_of_z,
-                                          growth_array=growth,
-                                          fgrowth_array=fgrowth)
+        #growth = ccl.background.growth_factor(cosmo, a)
+        #fgrowth = ccl.background.growth_rate(cosmo, a)
+
 
         if self.kmax:
             for pair in self._var_pairs:
                 # Get the matter power spectrum:
                 k, z, Pk_lin = self.provider.get_Pk_grid(var_pair=pair, nonlinear=False)
-
-                # np.flip(arr, axis=0) flips the rows of arr, thus making Pk with z
-                # in descending order.
                 Pk_lin = np.flip(Pk_lin, axis=0)
-                cosmo._set_linear_power_from_arrays(a, k, Pk_lin)
 
                 if self.nonlinear:
-                    k, z, Pk_nonlin = self.provider.get_Pk_grid(var_pair=pair, nonlinear=True)
+                    _, z, Pk_nonlin = self.provider.get_Pk_grid(var_pair=pair, nonlinear=True)
                     Pk_nonlin = np.flip(Pk_nonlin, axis=0)
-                    cosmo._set_nonlin_power_from_arrays(a, k, Pk_nonlin)
+
+                    # Create a CCL cosmology object. Because we are giving it background 
+                    # quantities, it should not depend on the cosmology parameters given
+                    cosmo = ccl.CosmologyCalculator(
+                        Omega_c=Omega_c, Omega_b=Omega_b, h=h, sigma8=0.8, n_s=0.96,
+                        background={'a': a,
+                                'chi': distance,
+                                'h_over_h0': E_of_z},
+                        pk_linear={'a': a,
+                        'k': k,
+                        'delta_matter:delta_matter': Pk_lin}, 
+                        pk_nonlin={'a': a,
+                        'k': k,
+                        'delta_matter:delta_matter': Pk_nonlin}
+                    )
+
+                else:
+                    cosmo = ccl.CosmologyCalculator(
+                        Omega_c=Omega_c, Omega_b=Omega_b, h=h, sigma8=0.8, n_s=0.96,
+                        background={'a': a,
+                                'chi': distance,
+                                'h_over_h0': E_of_z},
+                        pk_linear={'a': a,
+                        'k': k,
+                        'delta_matter:delta_matter': Pk_lin}
+                    )                    
 
         state['CCL'] = {'cosmo': cosmo}
         for required_result, method in self._required_results.items():
             state['CCL'][required_result] = method(cosmo)
 
     def get_CCL(self):
-        """
-        Get dictionary of CCL computed quantities.
-        results['cosmo'] contains the initialized CCL Cosmology object.
-        Other entries are computed by methods passed in as the requirements
-
-        :return: dict of results
-        """
-        return self._current_state['CCL']
+        return self._current_state['CCL']

soliket/cross_correlation.py

@@ -0,0 +1,59 @@
+"""
+Simple likelihood for CMB-galaxy cross-correlations using the cobaya
+CCL module. 
+
+First version by Pablo Lemos
+"""
+
+import numpy as np
+from .gaussian import GaussianData, GaussianLikelihood
+import pyccl as ccl
+
+class CrossCorrelationLikelihood(GaussianLikelihood):
+    def initialize(self):
+        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)
+
+    def get_requirements(self):
+        return {'CCL': {"kmax": 10,
+                        "nonlinear": True}}
+
+    def _get_data(self, **params_values):
+        data_auto = np.loadtxt(self.auto_file)
+        data_cross = np.loadtxt(self.cross_file)
+
+        # Get data
+        self.ell_auto = data_auto[0]
+        cl_auto = data_auto[1]
+        cl_auto_err = data_auto[2]
+
+        self.ell_cross = data_cross[0]
+        cl_cross = data_cross[1]
+        cl_cross_err = data_cross[2]  
+
+        x = np.concatenate([self.ell_auto, self.ell_cross])
+        y = np.concatenate([cl_auto, cl_cross])
+        dy = np.concatenate([cl_auto_err, cl_cross_err])
+
+        return x, y, dy
+
+    def _get_theory(self, **params_values):
+        cosmo = self.provider.get_CCL()['cosmo']
+
+        tracer_g = ccl.NumberCountsTracer(cosmo, has_rsd=False, dndz = self.dndz.T, 
+                                            bias =(self.dndz[:,0], params_values['b1']*np.ones(len(self.dndz[:,0]))), 
+                                            mag_bias = (self.dndz[:,0], params_values['s1']*np.ones(len(self.dndz[:,0])))
+                                            )
+        tracer_k = ccl.CMBLensingTracer(cosmo, z_source = 1060)
+
+        cl_gg = ccl.cls.angular_cl(cosmo, tracer_g, tracer_g, self.ell_auto) #+ 1e-7
+        cl_kg = ccl.cls.angular_cl(cosmo, tracer_k, tracer_g, self.ell_cross)    
+
+        return np.concatenate([cl_gg, cl_kg])
+
+    def logp(self, **params_values):
+        theory = self._get_theory(**params_values)
+        return self.data.loglike(theory)

soliket/data/xcorr_simulated/clgg_noiseless.txt

@@ -0,0 +1,3 @@
+3.550000000000000000e+01 7.650000000000000000e+01 1.265000000000000000e+02 1.765000000000000000e+02 2.265000000000000000e+02 2.765000000000000000e+02 3.265000000000000000e+02 3.765000000000000000e+02 4.265000000000000000e+02 4.765000000000000000e+02 5.265000000000000000e+02 5.760000000000000000e+02
+2.171386555992278359e-07 2.074818774537804928e-07 1.760569576742365607e-07 1.556208570083821835e-07 1.434590789832146224e-07 1.347070083607898911e-07 1.277731731363569596e-07 1.224832434049751676e-07 1.186677876068538790e-07 1.159080212760721562e-07 1.137746016890491462e-07 1.120196152233489823e-07
+2.171386555992278359e-08 2.074818774537804994e-08 1.760569576742365674e-08 1.556208570083821835e-08 1.434590789832146224e-08 1.347070083607899044e-08 1.277731731363569629e-08 1.224832434049751775e-08 1.186677876068538790e-08 1.159080212760721661e-08 1.137746016890491529e-08 1.120196152233489823e-08

soliket/data/xcorr_simulated/clkg_noiseless.txt

@@ -0,0 +1,3 @@
+3.550000000000000000e+01 7.650000000000000000e+01 1.265000000000000000e+02 1.765000000000000000e+02 2.265000000000000000e+02 2.765000000000000000e+02 3.265000000000000000e+02 3.765000000000000000e+02 4.265000000000000000e+02 4.765000000000000000e+02 5.265000000000000000e+02 5.760000000000000000e+02
+1.424982171873059867e-07 1.234732647730792034e-07 8.611508401982061149e-08 6.248438515784781822e-08 4.837007692483073317e-08 3.852305256406458136e-08 3.090458671774709724e-08 2.505965588344152354e-08 2.078012214338600926e-08 1.768136302352318810e-08 1.532330752048001726e-08 1.340720288202413845e-08
+1.424982171873059966e-08 1.234732647730792068e-08 8.611508401982062142e-09 6.248438515784782318e-09 4.837007692483073648e-09 3.852305256406458467e-09 3.090458671774709973e-09 2.505965588344152520e-09 2.078012214338600843e-09 1.768136302352318975e-09 1.532330752048001851e-09 1.340720288202413845e-09