"""
:Synopsis: Likelihood for cross-correlation of CMB lensing with Large Scale Structure
data. Makes use of the cobaya CCL module for handling tracers and Limber integration.
:Authors: Pablo Lemos, Ian Harrison.
"""
from typing import ClassVar
import numpy as np
try:
from numpy import trapezoid
except ImportError:
from numpy import trapz as trapezoid
from cobaya.log import LoggedError
from cobaya.theory import Provider
from soliket.ccl import CCL
from soliket.gaussian import GaussianLikelihood
[docs]
class CCLTracersLikelihood(GaussianLikelihood):
r"""
Generic likelihood for CCL tracer objects.
"""
datapath: str
use_spectra: str | list[tuple[str, str]]
ncovsims: int | None
provider: Provider
def initialize(self):
super().initialize()
[docs]
def get_requirements(self) -> dict:
return {"CCL": {"kmax": 10, "nonlinear": True}, "zstar": None}
def _get_CCL_results(self) -> tuple[CCL, dict]:
cosmo_dict = self.provider.get_CCL()
return cosmo_dict["ccl"], cosmo_dict["cosmo"]
def _get_nz(
self, z: np.ndarray, tracer, tracer_name: str, **params_values
) -> np.ndarray:
if self.z_nuisance_mode == "deltaz":
bias = params_values[f"{tracer_name}_deltaz"]
nz_biased = tracer.get_dndz(z - bias)
# nz_biased /= np.trapezoid(nz_biased, z)
return nz_biased
def _get_ia_bias(
self,
z_tracer: np.ndarray,
nz_tracer: np.ndarray,
tracer_name: str,
params_values: dict,
) -> tuple[np.ndarray, np.ndarray] | None:
if self.ia_mode is None:
return None
elif self.ia_mode == "nla":
A_IA = params_values["A_IA"]
eta_IA = params_values["eta_IA"]
z0_IA = trapezoid(z_tracer * nz_tracer)
return (z_tracer, A_IA * ((1 + z_tracer) / (1 + z0_IA)) ** eta_IA)
elif self.ia_mode == "nla-perbin":
A_IA = params_values[f"{tracer_name}_A_IA"]
return (z_tracer, A_IA * np.ones_like(z_tracer))
elif self.ia_mode == "nla-noevo":
A_IA = params_values["A_IA"]
return (z_tracer, A_IA * np.ones_like(z_tracer))
return None
class CCLTracersCrossLikelihood(CCLTracersLikelihood):
r"""
Generic likelihood for cross-correlations of CCL tracer objects.
"""
datapath: str
use_spectra: str | list[tuple[str, str]]
ncovsims: int | None
provider: Provider
def initialize(self):
super().initialize()
self._check_is_cross()
def _check_is_cross(self):
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 {} to calculate an \
autocorrelation, but it is a cross-correlation \
likelihood. Please check your tracer selection in the \
ini file.".format(self.__class__.__name__),
)
class CCLTracersAutoLikelihood(CCLTracersLikelihood):
r"""
Generic likelihood for auto-correlations of CCL tracer objects.
"""
datapath: str
use_spectra: str | list[tuple[str, str]]
ncovsims: int | None
provider: Provider
def initialize(self):
super().initialize()
self._check_is_auto()
def _check_is_auto(self):
for tracer_comb in self.sacc_data.get_tracer_combinations():
if not (
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 {} to calculate a \
cross-correlation, but it is an auto-correlation \
likelihood. Please check your tracer selection in the \
ini file.".format(self.__class__.__name__),
)
[docs]
class GalaxyKappaLikelihood(CCLTracersCrossLikelihood):
r"""
Likelihood for cross-correlations of galaxy and CMB lensing data.
"""
name: str = "GalaxyKappa"
_allowable_tracers: ClassVar[list[str]] = ["cmb_convergence", "galaxy_density"]
params: dict
def _get_theory(self, **params_values) -> np.ndarray:
ccl, cosmo = self._get_CCL_results()
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=(z_gal_tracer, nz_gal_tracer),
bias=(z_gal_tracer, params_values["b1"] * np.ones(len(z_gal_tracer))),
mag_bias=(z_gal_tracer, params_values["s1"] * np.ones(len(z_gal_tracer))),
)
tracer_k = ccl.CMBLensingTracer(cosmo, z_source=self.provider.get_param("zstar"))
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)
cl_gk_binned = np.dot(w_bins_gk, cl_gk_unbinned)
return cl_gk_binned
[docs]
class ShearKappaLikelihood(CCLTracersCrossLikelihood):
r"""
Likelihood for cross-correlations of galaxy weak lensing shear and CMB lensing data.
"""
name: str = "ShearKappa"
_allowable_tracers: ClassVar[list[str]] = ["cmb_convergence", "galaxy_shear"]
z_nuisance_mode: str | bool | None
m_nuisance_mode: str | bool | None
ia_mode: str | None
params: dict
def _get_theory(self, **params_values) -> np.ndarray:
ccl, cosmo = self._get_CCL_results()
cl_binned_list: list[np.ndarray] = []
for tracer_comb in self.sacc_data.get_tracer_combinations():
if self.sacc_data.tracers[tracer_comb[0]].quantity == "cmb_convergence":
tracer1 = ccl.CMBLensingTracer(
cosmo, z_source=self.provider.get_param("zstar")
)
elif self.sacc_data.tracers[tracer_comb[0]].quantity == "galaxy_shear":
sheartracer_name = tracer_comb[0]
z_tracer1 = self.sacc_data.tracers[tracer_comb[0]].z
nz_tracer1 = self.sacc_data.tracers[tracer_comb[0]].nz
ia_z = self._get_ia_bias(
z_tracer1, nz_tracer1, sheartracer_name, params_values
)
tracer1 = ccl.WeakLensingTracer(
cosmo, dndz=(z_tracer1, nz_tracer1), ia_bias=ia_z
)
if self.z_nuisance_mode is not None:
nz_tracer1 = self._get_nz(
z_tracer1, tracer1, tracer_comb[0], **params_values
)
tracer1 = ccl.WeakLensingTracer(
cosmo, dndz=(z_tracer1, nz_tracer1), ia_bias=ia_z
)
if self.sacc_data.tracers[tracer_comb[1]].quantity == "cmb_convergence":
tracer2 = ccl.CMBLensingTracer(
cosmo, z_source=self.provider.get_param("zstar")
)
elif self.sacc_data.tracers[tracer_comb[1]].quantity == "galaxy_shear":
sheartracer_name = tracer_comb[1]
z_tracer2 = self.sacc_data.tracers[tracer_comb[1]].z
nz_tracer2 = self.sacc_data.tracers[tracer_comb[1]].nz
ia_z = self._get_ia_bias(
z_tracer2, nz_tracer2, sheartracer_name, params_values
)
tracer2 = ccl.WeakLensingTracer(
cosmo, dndz=(z_tracer2, nz_tracer2), ia_bias=ia_z
)
if self.z_nuisance_mode is not None:
nz_tracer2 = self._get_nz(
z_tracer2, tracer2, tracer_comb[1], **params_values
)
tracer2 = ccl.WeakLensingTracer(
cosmo, dndz=(z_tracer2, nz_tracer2), ia_bias=ia_z
)
bpw_idx = self.sacc_data.indices(tracers=tracer_comb)
bpw = self.sacc_data.get_bandpower_windows(bpw_idx)
ells_theory = np.asarray(bpw.values, dtype=int)
w_bins = bpw.weight.T
cl_unbinned = ccl.cells.angular_cl(cosmo, tracer1, tracer2, ells_theory)
if self.m_nuisance_mode is not None:
# note this allows wrong calculation, as we can do
# shear x shear if the spectra are in the sacc
# but then we would want (1 + m1) * (1 + m2)
m_bias = params_values[f"{sheartracer_name}_m"]
cl_unbinned = (1 + m_bias) * cl_unbinned
cl_binned = np.dot(w_bins, cl_unbinned)
cl_binned_list.append(cl_binned)
cl_binned_total = np.concatenate(cl_binned_list)
return cl_binned_total
def _get_tracer(self, ccl: CCL, cosmo: dict, tracer_name: str, params_values: dict):
tracer_data = self.sacc_data.tracers[tracer_name]
if tracer_data.quantity == "cmb_convergence":
return ccl.CMBLensingTracer(cosmo, z_source=self.provider.get_param("zstar"))
elif tracer_data.quantity == "galaxy_shear":
z_tracer = self.sacc_data.tracers[tracer_name].z
nz_tracer = self.sacc_data.tracers[tracer_name].nz
ia_z = self._get_ia_bias(z_tracer, nz_tracer, tracer_name, params_values)
tracer = ccl.WeakLensingTracer(
cosmo, dndz=(z_tracer, nz_tracer), ia_bias=ia_z
)
if self.z_nuisance_mode is not None:
nz_tracer = self._get_nz(z_tracer, tracer, tracer_name, **params_values)
tracer = ccl.WeakLensingTracer(
cosmo, dndz=(z_tracer, nz_tracer), ia_bias=ia_z
)
return tracer
return None