diff --git a/docs/conf.py b/docs/conf.py index cee05aa9..c38cd049 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -5,42 +5,43 @@ # Import SOLikeT (for autodoc) import sys + sys.path.insert(0, "..") # Create some mock imports from unittest import mock + MOCK_MODULES = ["cosmopower", "tensorflow", "pyccl", "camb"] for module in MOCK_MODULES: sys.modules[module] = mock.Mock() import soliket -__all__ = ['soliket'] +__all__ = ["soliket"] # -- Project information ----------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information -project = 'SOLikeT' -copyright = '2023, The SO Collaboration' -author = 'The SO Collaboration' +project = "SOLikeT" +copyright = "2023, The SO Collaboration" +author = "The SO Collaboration" # -- General configuration --------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration extensions = [ - "sphinx.ext.autodoc", # Generate doc pages from source docstrings - "sphinx.ext.viewcode", # Generate links to source code - "sphinx.ext.mathjax", # Mathematical symbols - "sphinx_rtd_theme", # readthedocs theme + "sphinx.ext.autodoc", # Generate doc pages from source docstrings + "sphinx.ext.viewcode", # Generate links to source code + "sphinx.ext.mathjax", # Mathematical symbols + "sphinx_rtd_theme", # readthedocs theme ] -templates_path = ['_templates'] -exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] - +templates_path = ["_templates"] +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] # -- Options for HTML output ------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output -html_theme = 'sphinx_rtd_theme' -html_static_path = ['_static'] +html_theme = "sphinx_rtd_theme" +html_static_path = ["_static"] diff --git a/examples/example_2.yaml b/examples/example_2.yaml new file mode 100644 index 00000000..f0fbcdfa --- /dev/null +++ b/examples/example_2.yaml @@ -0,0 +1,124 @@ +debug: False + +# Choose a sampler. +# Running the evaluate sampler will compute the likelihood +# only at a single value of parameters +sampler: + evaluate: null + +# Specify the output location and file prefix +output: output/example_2 + +# Specify which Likelihoods we would like to use +# Any options for the Likelihood will use their default values unless specified here +# In this case the options and defaults are specified in soliket/mflike/MFLike.yaml +# Note that MFLike is a Cobaya `installable likelihood`. +# When running this yaml file, or calling `cobaya-install example_1.yaml` the required +# installable components will automatically be downloaded and installed. +# Note that for the soliket MFLike likelihood we are required to calculate: +# - CMB theory power spectra (from CAMB theory below) +# - Multi-frequency bandpass calibrations (from soliket.BandPass theory below) +# - Multi-frequency foregrounds (from soliket.Foreground theory below) +# - The combination of the above components (from soliket.TheoryForge_MFLike theory below) +likelihood: + soliket.MFLike: + data_folder: MFLike/v0.8 + input_file: LAT_simu_sacc_00000.fits + cov_Bbl_file: data_sacc_w_covar_and_Bbl.fits + lmax_theory: 1500 + defaults: + polarizations: ['TT', 'TE', 'ET', 'EE'] + scales: + TT: [30, 1500] + TE: [30, 1500] + ET: [30, 1500] + EE: [30, 1500] + symmetrize: False + lcuts: 1500 + +# Specify the Theory codes which will compute observables to be compared with the data +# in the Likelihood. +# Here we specify the CAMB Einstein-Boltzmann code, with a number of choices made on +# precision and neutrino model. +theory: + soliket.CosmoPower: + network_path: soliket/cosmopower/data/CP_paper + network_settings: + tt: + type: NN + log: True + filename: cmb_TT_NN + has_ell_factor: False + ee: + type: NN + log: True + filename: cmb_EE_NN + has_ell_factor: False + te: + type: PCAplusNN + log: False + filename: cmb_TE_PCAplusNN + has_ell_factor: False + + soliket.BandPass: + stop_at_error: True + + soliket.Foreground: + stop_at_error: True + spectra: + polarizations: ["tt", "te", "ee"] + lmin: 2 + lmax: 1500 + exp_ch: ["LAT_150"] + eff_freqs: [150.] + + soliket.TheoryForge_MFLike: + spectra: + polarizations: ['tt', 'te', 'ee'] + lmin: 2 + lmax: 1500 + stop_at_error: True + +# Specify the parameter values at which to compute the likelihood +params: + H0: + latex: H_0 + derived: 'lambda h: h * 100.' + h: + latex: h_0 + value: 0.677 + logA: + # Dropped parameters are sampled but not passed to Likelihood/Theory codes + # Here it is As (specified below) which is passed to the Likelihood/Theory code. + # drop: true + latex: \log(10^{10} A_\mathrm{s}) + value: 3.05 + As: + value: 'lambda logA: 1e-10*np.exp(logA)' + ombh2: + latex: \Omega_\mathrm{b} h^2 + value: 0.0224 + omch2: + latex: \Omega_c h^2 + value: 0.1202 + ns: + latex: n_s + value: 0.9649 + Alens: + latex: A_lens + value: 1.0 + tau: + latex: \tau_\mathrm{reio} + value: 0.0554 + mnu1: + value: 0.0 + drop: True + mnu2: + value: 'lambda: np.sqrt(7.5e-5)' + drop: True + mnu3: + value: 'lambda: np.sqrt(2.5e-3)' + drop: True + mnu: + value: 'lambda mnu1, mnu2 ,mnu3: mnu1 + mnu2 + mnu3' + latex: \sum m_\nu diff --git a/setup.py b/setup.py index 4296ae2c..27e88b68 100644 --- a/setup.py +++ b/setup.py @@ -26,7 +26,7 @@ https://soliket.readthedocs.io/en/latest/developers.html#checking-code-in-development """ -if 'test' in sys.argv: +if "test" in sys.argv: print(TEST_HELP) sys.exit(1) @@ -44,7 +44,7 @@ https://soliket.readthedocs.io/en/latest/developers.html#documentation """ -if 'build_docs' in sys.argv or 'build_sphinx' in sys.argv: +if "build_docs" in sys.argv or "build_sphinx" in sys.argv: print(DOCS_HELP) sys.exit(1) @@ -59,5 +59,9 @@ version = '{version}' """.lstrip() -setup(use_scm_version={'write_to': os.path.join('soliket', 'version.py'), - 'write_to_template': VERSION_TEMPLATE}) +setup( + use_scm_version={ + "write_to": os.path.join("soliket", "version.py"), + "write_to_template": VERSION_TEMPLATE, + } +) diff --git a/soliket/__init__.py b/soliket/__init__.py index f58c1e77..6e9631fb 100644 --- a/soliket/__init__.py +++ b/soliket/__init__.py @@ -11,12 +11,16 @@ from .ccl import CCL from .clusters import ClusterLikelihood from .cosmopower import CosmoPower, CosmoPowerDerived -from .cross_correlation import (CrossCorrelationLikelihood, - GalaxyKappaLikelihood, ShearKappaLikelihood) +from .cross_correlation import ( + CrossCorrelationLikelihood, + GalaxyKappaLikelihood, + ShearKappaLikelihood, +) from .foreground import Foreground from .gaussian import GaussianLikelihood, MultiGaussianLikelihood from .lensing import LensingLikelihood, LensingLiteLikelihood from .mflike import MFLike, TheoryForge_MFLike + # from .studentst import StudentstLikelihood from .ps import BinnedPSLikelihood, PSLikelihood from .xcorr import XcorrLikelihood diff --git a/soliket/bandpass/bandpass.py b/soliket/bandpass/bandpass.py index 149a53b3..d2c6ddf1 100644 --- a/soliket/bandpass/bandpass.py +++ b/soliket/bandpass/bandpass.py @@ -127,10 +127,11 @@ class BandPass(Theory): external_bandpass: dict def initialize(self): - - self.expected_params_bp = ["bandint_shift_LAT_93", - "bandint_shift_LAT_145", - "bandint_shift_LAT_225"] + self.expected_params_bp = [ + "bandint_shift_LAT_93", + "bandint_shift_LAT_145", + "bandint_shift_LAT_225", + ] self.exp_ch = None # self.eff_freqs = None @@ -148,28 +149,36 @@ def initialize(self): self.bandint_external_bandpass = bool(self.external_bandpass) if self.bandint_external_bandpass: - path = os.path.normpath(os.path.join(self.data_folder, - self.external_bandpass["path"])) + path = os.path.normpath( + os.path.join(self.data_folder, self.external_bandpass["path"]) + ) arrays = os.listdir(path) self._init_external_bandpass_construction(path, arrays) - if (not self.read_from_sacc and not self.use_top_hat_band - and not self.bandint_external_bandpass): + if ( + not self.read_from_sacc + and not self.use_top_hat_band + and not self.bandint_external_bandpass + ): raise LoggedError( - self.log, "fill the dictionaries in the yaml file for" - "either reading the passband from sacc file (mflike default)" - "or an external passband or building a top-hat one!" + self.log, + "fill the dictionaries in the yaml file for" + "either reading the passband from sacc file (mflike default)" + "or an external passband or building a top-hat one!", ) def initialize_with_params(self): # Check that the parameters are the right ones differences = are_different_params_lists( - self.input_params, self.expected_params_bp, - name_A="given", name_B="expected") + self.input_params, + self.expected_params_bp, + name_A="given", + name_B="expected", + ) if differences: raise LoggedError( - self.log, "Configuration error in parameters: %r.", - differences) + self.log, "Configuration error in parameters: %r.", differences + ) def must_provide(self, **requirements): # bandint_freqs is required by Foreground @@ -177,8 +186,9 @@ def must_provide(self, **requirements): # Assign those from Foreground if "bandint_freqs" in requirements: self.bands = requirements["bandint_freqs"]["bands"] - self.exp_ch = [k.replace("_s0", "") for k in self.bands.keys() - if "_s0" in k] + self.exp_ch = [ + k.replace("_s0", "") for k in self.bands.keys() if "_s0" in k + ] def calculate(self, state, want_derived=False, **params_values_dict): r""" @@ -231,7 +241,7 @@ def _bandpass_construction(self, **params): data_are_monofreq = False bandint_freqs = [] for ifr, fr in enumerate(self.exp_ch): - bandpar = 'bandint_shift_' + str(fr) + bandpar = "bandint_shift_" + str(fr) bands = self.bands[f"{fr}_s0"] nu_ghz, bp = np.asarray(bands["nu"]), np.asarray(bands["bandpass"]) if self.use_top_hat_band: @@ -240,19 +250,25 @@ def _bandpass_construction(self, **params): self.bandint_width = np.full_like( self.exp_ch, self.bandint_width, dtype=float ) - if self.bandint_nsteps > 1 and np.any(np.array(self.bandint_width) == 0): + if self.bandint_nsteps > 1 and np.any( + np.array(self.bandint_width) == 0 + ): raise LoggedError( - self.log, "One band has width = 0, \ - set a positive width and run again" + self.log, + "One band has width = 0, \ + set a positive width and run again", ) # Compute central frequency given bandpass fr = nu_ghz @ bp / bp.sum() if self.bandint_nsteps > 1: - bandlow = fr * (1 - self.bandint_width[ifr] * .5) - bandhigh = fr * (1 + self.bandint_width[ifr] * .5) - nub = np.linspace(bandlow + params[bandpar], - bandhigh + params[bandpar], - self.bandint_nsteps, dtype=float) + bandlow = fr * (1 - self.bandint_width[ifr] * 0.5) + bandhigh = fr * (1 + self.bandint_width[ifr] * 0.5) + nub = np.linspace( + bandlow + params[bandpar], + bandhigh + params[bandpar], + self.bandint_nsteps, + dtype=float, + ) tranb = _cmb2bb(nub) tranb_norm = np.trapz(_cmb2bb(nub), nub) bandint_freqs.append([nub, tranb / tranb_norm]) @@ -313,7 +329,9 @@ def _external_bandpass_construction(self, **params): if not hasattr(bp, "__len__"): bandint_freqs.append(nub) bandint_freqs = np.asarray(bandint_freqs) - self.log.info("bandpass is delta function, no band integration performed") + self.log.info( + "bandpass is delta function, no band integration performed" + ) else: trans_norm = np.trapz(bp * _cmb2bb(nub), nub) trans = bp / trans_norm * _cmb2bb(nub) diff --git a/soliket/bias/bias.py b/soliket/bias/bias.py index 94a14258..8e4a2b8a 100644 --- a/soliket/bias/bias.py +++ b/soliket/bias/bias.py @@ -34,7 +34,7 @@ class Bias(Theory): """Parent class for bias models.""" _logz = np.linspace(-3, np.log10(1100), 150) - _default_z_sampling = 10 ** _logz + _default_z_sampling = 10**_logz _default_z_sampling[0] = 0 def initialize(self): @@ -47,32 +47,40 @@ def must_provide(self, **requirements): options = requirements.get("linear_bias") or {} self.kmax = max(self.kmax, options.get("kmax", self.kmax)) - self.z = np.unique(np.concatenate( - (np.atleast_1d(options.get("z", self._default_z_sampling)), - np.atleast_1d(self.z)))) + self.z = np.unique( + np.concatenate( + ( + np.atleast_1d(options.get("z", self._default_z_sampling)), + np.atleast_1d(self.z), + ) + ) + ) # Dictionary of the things needed from CAMB/CLASS needs = {} self.nonlinear = self.nonlinear or options.get("nonlinear", False) self._var_pairs.update( - set((x, y) for x, y in - options.get("vars_pairs", [("delta_tot", "delta_tot")]))) + set( + (x, y) + for x, y in options.get("vars_pairs", [("delta_tot", "delta_tot")]) + ) + ) needs["Pk_grid"] = { "vars_pairs": self._var_pairs or [("delta_tot", "delta_tot")], "nonlinear": (True, False) if self.nonlinear else False, "z": self.z, - "k_max": self.kmax + "k_max": self.kmax, } assert len(self._var_pairs) < 2, "Bias doesn't support other Pk yet" return needs def _get_Pk_mm(self): - self.k, self.z, Pk_mm = \ - self.provider.get_Pk_grid(var_pair=list(self._var_pairs)[0], - nonlinear=self.nonlinear) + self.k, self.z, Pk_mm = self.provider.get_Pk_grid( + var_pair=list(self._var_pairs)[0], nonlinear=self.nonlinear + ) return Pk_mm def get_Pk_gg_grid(self) -> dict: @@ -89,9 +97,8 @@ class Linear_bias(Bias): Has one free parameter, :math:`b_\mathrm{lin}` (``b_lin``). """ - def calculate(self, state: dict, want_derived: bool = True, - **params_values_dict): + def calculate(self, state: dict, want_derived: bool = True, **params_values_dict): Pk_mm = self._get_Pk_mm() - state["Pk_gg_grid"] = params_values_dict["b_lin"] ** 2. * Pk_mm + state["Pk_gg_grid"] = params_values_dict["b_lin"] ** 2.0 * Pk_mm state["Pk_gm_grid"] = params_values_dict["b_lin"] * Pk_mm diff --git a/soliket/cash/cash.py b/soliket/cash/cash.py index 3f071a93..8a0368b1 100644 --- a/soliket/cash/cash.py +++ b/soliket/cash/cash.py @@ -7,12 +7,12 @@ # Likelihood for independent Poisson-distributed data # (here called Cash-C, see https://arxiv.org/abs/1912.05444) + class CashCLikelihood(Likelihood): name: str = "Cash-C" datapath = Optional[str] def initialize(self): - x, N = self._get_data() self.data = CashCData(self.name, N) diff --git a/soliket/cash/cash_data.py b/soliket/cash/cash_data.py index 4075c122..c25fbcfb 100644 --- a/soliket/cash/cash_data.py +++ b/soliket/cash/cash_data.py @@ -8,12 +8,15 @@ def cash_c_logpdf(theory, data, usestirling=True): ln_fac = np.zeros_like(data, dtype=float) if usestirling: # use Stirling's approximation for N > 10 - ln_fac[data > 10] = 0.918939 + (data[data > 10] + 0.5) \ - * np.log(data[data > 10]) - data[data > 10] + ln_fac[data > 10] = ( + 0.918939 + + (data[data > 10] + 0.5) * np.log(data[data > 10]) + - data[data > 10] + ) ln_fac[data <= 10] = np.log(factorial(data[data <= 10])) else: ln_fac[data > 0] = np.log(factorial(data[data > 0])) - ln_fac[data == 0] = 0. + ln_fac[data == 0] = 0.0 loglike = data * np.log(theory) - theory - ln_fac @@ -21,8 +24,7 @@ def cash_c_logpdf(theory, data, usestirling=True): class CashCData: - """Named multi-dimensional Cash-C distributed data - """ + """Named multi-dimensional Cash-C distributed data""" def __init__(self, name, N, usestirling=True): self.name = str(name) diff --git a/soliket/ccl/ccl.py b/soliket/ccl/ccl.py index ec34bbeb..ef6db3e2 100644 --- a/soliket/ccl/ccl.py +++ b/soliket/ccl/ccl.py @@ -86,8 +86,9 @@ class CCL(Theory): """A theory code wrapper for CCL.""" + _logz = np.linspace(-3, np.log10(1100), 150) - _default_z_sampling = 10 ** _logz + _default_z_sampling = 10**_logz _default_z_sampling[0] = 0 kmax: float z: np.ndarray @@ -97,7 +98,9 @@ def initialize(self) -> None: try: import pyccl as ccl except ImportError: - raise LoggedError(self.log, "Could not import ccl. Install pyccl to use ccl.") + raise LoggedError( + self.log, "Could not import ccl. Install pyccl to use ccl." + ) else: self.ccl = ccl @@ -107,46 +110,54 @@ def initialize(self) -> None: def get_requirements(self) -> set: # These are currently required to construct a CCL cosmology object. # Ultimately CCL should depend only on observable not parameters - return {'omch2', 'ombh2'} + return {"omch2", "ombh2"} def must_provide(self, **requirements) -> dict: # requirements is dictionary of things requested by likelihoods # Note this may be called more than once - if 'CCL' not in requirements: + if "CCL" not in requirements: return {} - options = requirements.get('CCL') or {} - if 'methods' in options: - self._required_results.update(options['methods']) - - self.kmax = max(self.kmax, options.get('kmax', self.kmax)) - self.z = np.unique(np.concatenate( - (np.atleast_1d(options.get("z", self._default_z_sampling)), - np.atleast_1d(self.z)))) + options = requirements.get("CCL") or {} + if "methods" in options: + self._required_results.update(options["methods"]) + + self.kmax = max(self.kmax, options.get("kmax", self.kmax)) + self.z = np.unique( + np.concatenate( + ( + np.atleast_1d(options.get("z", self._default_z_sampling)), + np.atleast_1d(self.z), + ) + ) + ) # Dictionary of the things CCL needs from CAMB/CLASS needs = {} if self.kmax: - self.nonlinear = self.nonlinear or options.get('nonlinear', False) + self.nonlinear = self.nonlinear or options.get("nonlinear", False) # CCL currently only supports ('delta_tot', 'delta_tot'), but call allow # general as placeholder self._var_pairs.update( - set((x, y) for x, y in - options.get('vars_pairs', [('delta_tot', 'delta_tot')]))) - - needs['Pk_grid'] = { - 'vars_pairs': self._var_pairs or [('delta_tot', 'delta_tot')], - 'nonlinear': (True, False) if self.nonlinear else False, - 'z': self.z, - 'k_max': self.kmax + set( + (x, y) + for x, y in options.get("vars_pairs", [("delta_tot", "delta_tot")]) + ) + ) + + needs["Pk_grid"] = { + "vars_pairs": self._var_pairs or [("delta_tot", "delta_tot")], + "nonlinear": (True, False) if self.nonlinear else False, + "z": self.z, + "k_max": self.kmax, } - needs['Hubble'] = {'z': self.z} - needs['comoving_radial_distance'] = {'z': self.z} + needs["Hubble"] = {"z": self.z} + needs["comoving_radial_distance"] = {"z": self.z} - needs['fsigma8'] = {'z': self.z} - needs['sigma8_z'] = {'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 @@ -155,8 +166,7 @@ def get_can_support_params(self) -> Sequence[str]: # return any nuisance parameters that CCL can support return [] - def calculate(self, state: dict, want_derived: bool = True, - **params_values_dict): + def calculate(self, state: dict, want_derived: bool = True, **params_values_dict): # calculate the general CCL cosmo object which likelihoods can then use to get # what they need (likelihoods should cache results appropriately) # get our requirements from self.provider @@ -167,8 +177,8 @@ def calculate(self, state: dict, want_derived: bool = True, 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 + 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. @@ -178,7 +188,7 @@ def calculate(self, state: dict, want_derived: bool = True, # Array z is sorted in ascending order. CCL requires an ascending scale # factor as input - a = 1. / (1 + self.z[::-1]) + a = 1.0 / (1 + self.z[::-1]) # growth = ccl.background.growth_factor(cosmo, a) # fgrowth = ccl.background.growth_rate(cosmo, a) @@ -189,8 +199,9 @@ def calculate(self, state: dict, want_derived: bool = True, Pk_lin = np.flip(Pk_lin, axis=0) if self.nonlinear: - _, 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) # Create a CCL cosmology object. Because we are giving it background @@ -201,15 +212,17 @@ def calculate(self, state: dict, want_derived: bool = True, 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}, # noqa E501 - pk_nonlin={'a': a, - 'k': k, - 'delta_matter:delta_matter': Pk_nonlin} # noqa E501 + background={"a": a, "chi": distance, "h_over_h0": E_of_z}, + pk_linear={ + "a": a, + "k": k, + "delta_matter:delta_matter": Pk_lin, + }, # noqa E501 + pk_nonlin={ + "a": a, + "k": k, + "delta_matter:delta_matter": Pk_nonlin, + }, # noqa E501 ) else: @@ -219,17 +232,17 @@ def calculate(self, state: dict, want_derived: bool = True, 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} # noqa E501 + background={"a": a, "chi": distance, "h_over_h0": E_of_z}, + pk_linear={ + "a": a, + "k": k, + "delta_matter:delta_matter": Pk_lin, + }, # noqa E501 ) - state['CCL'] = {'cosmo': cosmo, 'ccl': self.ccl} + state["CCL"] = {"cosmo": cosmo, "ccl": self.ccl} for required_result, method in self._required_results.items(): - state['CCL'][required_result] = method(cosmo) + state["CCL"][required_result] = method(cosmo) def get_CCL(self): - return self._current_state['CCL'] + return self._current_state["CCL"] diff --git a/soliket/clusters/clusters.py b/soliket/clusters/clusters.py index ae7a5390..8c4913b4 100644 --- a/soliket/clusters/clusters.py +++ b/soliket/clusters/clusters.py @@ -40,6 +40,7 @@ class ClusterLikelihood(PoissonLikelihood): """ Poisson Likelihood for un-binned :math:`y`-map galaxy cluster counts. """ + name = "Clusters" columns = ["tsz_signal", "z", "tsz_signal_err"] @@ -48,8 +49,11 @@ class ClusterLikelihood(PoissonLikelihood): def initialize(self): self.data_path = self.data_path or os.path.join( - self.get_class_path(), 'data', 'selFn_equD56') - self.data_name = os.path.join(self.get_class_path(), 'data', 'E-D56Clusters.fits') + self.get_class_path(), "data", "selFn_equD56" + ) + self.data_name = os.path.join( + self.get_class_path(), "data", "E-D56Clusters.fits" + ) self.zarr = np.arange(0, 2, 0.05) self.k = np.logspace(-4, np.log10(5), 200) @@ -58,8 +62,10 @@ def initialize(self): try: import pyccl as ccl except ImportError: - raise LoggedError(self.log, "Could not import ccl. " - "Install pyccl to use ClusterLikelihood.") + raise LoggedError( + self.log, + "Could not import ccl. " "Install pyccl to use ClusterLikelihood.", + ) else: self.ccl = ccl super().initialize() @@ -114,12 +120,12 @@ def _get_catalog(self): def _get_om(self): return (self.provider.get_param("omch2") + self.provider.get_param("ombh2")) / ( - (self.provider.get_param("H0") / 100.0) ** 2 + (self.provider.get_param("H0") / 100.0) ** 2 ) def _get_ob(self): return (self.provider.get_param("ombh2")) / ( - (self.provider.get_param("H0") / 100.0) ** 2 + (self.provider.get_param("H0") / 100.0) ** 2 ) def _get_Ez(self): @@ -149,7 +155,7 @@ def _get_HMF(self): ) # self.provider.get_Hubble(self.zarr) / self.provider.get_param("H0") om = self._get_om() - hmf = mf.HMF(om, Ez, pk=pks * h ** 3, kh=self.k / h, zarr=self.zarr) + hmf = mf.HMF(om, Ez, pk=pks * h**3, kh=self.k / h, zarr=self.zarr) return hmf @@ -199,7 +205,9 @@ def Prob_per_cluster(z, tsz_signal, tsz_signal_err): HMF.M, c_z, c_y * 1e-4, c_yerr * 1e-4, param_vals, Ez_fn, DA_fn ) - dn_dzdm = 10 ** np.squeeze(dn_dzdm_interp((np.log10(HMF.M), c_z))) * h ** 4.0 + dn_dzdm = ( + 10 ** np.squeeze(dn_dzdm_interp((np.log10(HMF.M), c_z))) * h**4.0 + ) ans = trapezoid(dn_dzdm * Pfunc_ind, dx=np.diff(HMF.M, axis=0), axis=0) return ans @@ -211,9 +219,9 @@ def _get_dVdz(self): DA_z = self.provider.get_angular_diameter_distance(self.zarr) dV_dz = ( - DA_z ** 2 - * (1.0 + self.zarr) ** 2 - / (self.provider.get_Hubble(self.zarr) / C_KM_S) + DA_z**2 + * (1.0 + self.zarr) ** 2 + / (self.provider.get_Hubble(self.zarr) / C_KM_S) ) # dV_dz *= (self.provider.get_param("H0") / 100.0) ** 3.0 # was h0 @@ -236,7 +244,7 @@ def _get_n_expected(self, **kwargs): Ntot = 0 dVdz = self._get_dVdz() - dn_dzdm = HMF.dn_dM(HMF.M, 500.0) * h ** 4.0 # getting rid of hs + dn_dzdm = HMF.dn_dM(HMF.M, 500.0) * h**4.0 # getting rid of hs for Yt, frac in zip(self.survey.Ythresh, self.survey.frac_of_survey): Pfunc = self.szutils.PfuncY(Yt, HMF.M, z_arr, param_vals, Ez_fn, DA_fn) @@ -244,11 +252,11 @@ def _get_n_expected(self, **kwargs): dn_dzdm * Pfunc, dx=np.diff(HMF.M[:, None] / h, axis=0), axis=0 ) Ntot += ( - trapezoid(N_z * dVdz, x=z_arr) - * 4.0 - * np.pi - * self.survey.fskytotal - * frac + trapezoid(N_z * dVdz, x=z_arr) + * 4.0 + * np.pi + * self.survey.fskytotal + * frac ) return Ntot @@ -265,15 +273,15 @@ def _test_n_tot(self, **kwargs): Ntot = 0 dVdz = self._get_dVdz() - dn_dzdm = HMF.dn_dM(HMF.M, 500.0) * h ** 4.0 # getting rid of hs + dn_dzdm = HMF.dn_dM(HMF.M, 500.0) * h**4.0 # getting rid of hs # Test Mass function against Nemo. Pfunc = 1.0 N_z = trapezoid(dn_dzdm * Pfunc, dx=np.diff(HMF.M[:, None] / h, axis=0), axis=0) Ntot = ( - trapezoid(N_z * dVdz, x=z_arr) - * 4.0 - * np.pi - * (600.0 / (4 * np.pi * (180 / np.pi) ** 2)) + trapezoid(N_z * dVdz, x=z_arr) + * 4.0 + * np.pi + * (600.0 / (4 * np.pi * (180 / np.pi) ** 2)) ) return Ntot diff --git a/soliket/clusters/massfunc.py b/soliket/clusters/massfunc.py index 0f944391..5553ec97 100644 --- a/soliket/clusters/massfunc.py +++ b/soliket/clusters/massfunc.py @@ -13,15 +13,15 @@ from .tinker import dn_dlogM -np.seterr(divide='ignore', invalid='ignore') +np.seterr(divide="ignore", invalid="ignore") class HMF: """ Build halo mass function """ - def __init__(self, om, Ez, pk=None, kh=None, zarr=None): + def __init__(self, om, Ez, pk=None, kh=None, zarr=None): # Initialize redshift and mass ranges if zarr is None: self.zarr = np.arange(0.05, 1.95, 0.1) @@ -32,7 +32,7 @@ def __init__(self, om, Ez, pk=None, kh=None, zarr=None): # self.M = 10**np.arange(13.5, 15.7, 0.02) M_edges = 10 ** np.arange(13.5, 15.72, 0.02) - self.M = (M_edges[1:] + M_edges[:-1]) / 2. # 10**np.arange(13.5, 15.7, 0.02) + self.M = (M_edges[1:] + M_edges[:-1]) / 2.0 # 10**np.arange(13.5, 15.7, 0.02) assert len(Ez) == len(zarr), "Ez and z arrays do not match" @@ -40,12 +40,13 @@ def __init__(self, om, Ez, pk=None, kh=None, zarr=None): # Initialize rho critical values for usage self.om = om - self.rho_crit0H100 = (3. / (8. * np.pi) * (100 * 1.e5) ** 2.) \ - / G_CGS * MPC2CM / MSUN_CGS + self.rho_crit0H100 = ( + (3.0 / (8.0 * np.pi) * (100 * 1.0e5) ** 2.0) / G_CGS * MPC2CM / MSUN_CGS + ) self.rhoc0om = self.rho_crit0H100 * self.om if pk is None: - print('this will not work') + print("this will not work") else: self.pk = pk self.kh = kh @@ -55,7 +56,7 @@ def rhoc(self): """ Critical density as a function of z """ - ans = self.rho_crit0H100 * self.E_z ** 2. + ans = self.rho_crit0H100 * self.E_z**2.0 return ans def rhom(self): @@ -77,8 +78,9 @@ def dn_dM(self, M, delta): :param delta: Threshold for critical density """ delts = self.critdensThreshold(delta) - dn_dlnm = dn_dlogM(M, self.zarr, self.rhoc0om, delts, self.kh, self.pk, - 'comoving') + dn_dlnm = dn_dlogM( + M, self.zarr, self.rhoc0om, delts, self.kh, self.pk, "comoving" + ) dn_dm = dn_dlnm / M[:, None] return dn_dm @@ -89,6 +91,7 @@ def inter_dndmLogm(self, delta, M=None): if M is None: M = self.M dndM = self.dn_dM(M, delta) - ans = RegularGridInterpolator((np.log10(M), self.zarr), - np.log10(dndM), method='cubic', fill_value=0) + ans = RegularGridInterpolator( + (np.log10(M), self.zarr), np.log10(dndM), method="cubic", fill_value=0 + ) return ans diff --git a/soliket/clusters/survey.py b/soliket/clusters/survey.py index 33617028..1f7ca18c 100644 --- a/soliket/clusters/survey.py +++ b/soliket/clusters/survey.py @@ -13,6 +13,7 @@ import astropy.table as atpy import numpy as np from astropy.io import fits + # from astLib import astWCS from astropy.wcs import WCS from scipy import interpolate diff --git a/soliket/clusters/sz_utils.py b/soliket/clusters/sz_utils.py index 800ba181..5b24471d 100644 --- a/soliket/clusters/sz_utils.py +++ b/soliket/clusters/sz_utils.py @@ -8,6 +8,7 @@ """ import numpy as np + try: from numpy import trapezoid except ImportError: @@ -15,23 +16,35 @@ from scipy import interpolate # from nemo import signals -from soliket.constants import (C_M_S, G_CGS, MPC2CM, MSUN_CGS, T_CMB, - electron_mass_kg, elementary_charge, h_Planck, - k_Boltzmann) +from soliket.constants import ( + C_M_S, + G_CGS, + MPC2CM, + MSUN_CGS, + T_CMB, + electron_mass_kg, + elementary_charge, + h_Planck, + k_Boltzmann, +) # from astropy.cosmology import FlatLambdaCDM # from .clusters import C_KM_S as C_in_kms -rho_crit0H100 = (3. / (8. * np.pi) * (100. * 1.e5) ** 2.) \ - / G_CGS * MPC2CM / MSUN_CGS +rho_crit0H100 = ( + (3.0 / (8.0 * np.pi) * (100.0 * 1.0e5) ** 2.0) / G_CGS * MPC2CM / MSUN_CGS +) def gaussian(xx, mu, sig, noNorm=False): if noNorm: - return np.exp(-1.0 * (xx - mu) ** 2 / (2.0 * sig ** 2.0)) + return np.exp(-1.0 * (xx - mu) ** 2 / (2.0 * sig**2.0)) else: - return 1.0 / (sig * np.sqrt(2 * np.pi)) \ - * np.exp(-1.0 * (xx - mu) ** 2 / (2.0 * sig ** 2.0)) + return ( + 1.0 + / (sig * np.sqrt(2 * np.pi)) + * np.exp(-1.0 * (xx - mu) ** 2 / (2.0 * sig**2.0)) + ) class szutils: @@ -55,9 +68,9 @@ def P_Yo(self, LgY, M, z, param_vals, Ez_fn, Da_fn): B0=param_vals["B0"], H0=param_vals["H0"], Ez_fn=Ez_fn, - Da_fn=Da_fn + Da_fn=Da_fn, ) - Y = 10 ** LgY + Y = 10**LgY # Ytilde = np.repeat(Ytilde[:, :, np.newaxis], LgY.shape[2], axis=2) @@ -66,8 +79,9 @@ def P_Yo(self, LgY, M, z, param_vals, Ez_fn, Da_fn): numer = -1.0 * (np.log(Y / Ytilde)) ** 2 ans = ( - 1.0 / (param_vals["scat"] * np.sqrt(2 * np.pi)) * - np.exp(numer / (2.0 * param_vals["scat"] ** 2)) + 1.0 + / (param_vals["scat"] * np.sqrt(2 * np.pi)) + * np.exp(numer / (2.0 * param_vals["scat"] ** 2)) ) return ans @@ -85,14 +99,15 @@ def P_Yo_vec(self, LgY, M, z, param_vals, Ez_fn, Da_fn): Ez_fn=Ez_fn, Da_fn=Da_fn, ) - Y = 10 ** LgY + Y = 10**LgY Ytilde = np.repeat(Ytilde[:, :, np.newaxis], LgY.shape[2], axis=2) numer = -1.0 * (np.log(Y / Ytilde)) ** 2 ans = ( - 1.0 / (param_vals["scat"] * np.sqrt(2 * np.pi)) * - np.exp(numer / (2.0 * param_vals["scat"] ** 2)) + 1.0 + / (param_vals["scat"] * np.sqrt(2 * np.pi)) + * np.exp(numer / (2.0 * param_vals["scat"] ** 2)) ) return ans @@ -103,11 +118,12 @@ def Y_erf(self, Y, Ynoise): return ans def P_of_gt_SN(self, LgY, MM, zz, Ynoise, param_vals, Ez_fn, Da_fn): - Y = 10 ** LgY + Y = 10**LgY sig_tr = np.outer(np.ones([MM.shape[0], MM.shape[1]]), self.Y_erf(Y, Ynoise)) - sig_thresh = np.reshape(sig_tr, - (MM.shape[0], MM.shape[1], len(self.Y_erf(Y, Ynoise)))) + sig_thresh = np.reshape( + sig_tr, (MM.shape[0], MM.shape[1], len(self.Y_erf(Y, Ynoise))) + ) LgYa = np.outer(np.ones([MM.shape[0], MM.shape[1]]), LgY) LgYa2 = np.reshape(LgYa, (MM.shape[0], MM.shape[1], len(LgY))) @@ -138,7 +154,7 @@ def P_of_Y_per(self, LgY, MM, zz, Y_c, Y_err, param_vals): return ans def Y_prob(self, Y_c, LgY, YNoise): - Y = 10 ** LgY + Y = 10**LgY ans = gaussian(Y, Y_c, YNoise) return ans @@ -196,9 +212,7 @@ def Pfunc_per_zarr(self, MM, z_c, Y_c, Y_err, int_HMF, param_vals): # ---------------------------------------------------------------------------------------- def calcR500Mpc(z, M500, Ez_fn, H0): - """Given z, M500 (in MSun), returns R500 in Mpc, with respect to critical density. - - """ + """Given z, M500 (in MSun), returns R500 in Mpc, with respect to critical density.""" if type(M500) == str: raise Exception( @@ -208,7 +222,7 @@ def calcR500Mpc(z, M500, Ez_fn, H0): Ez = Ez_fn(z) - criticalDensity = rho_crit0H100 * (H0 / 100.) ** 2 * Ez ** 2 + criticalDensity = rho_crit0H100 * (H0 / 100.0) ** 2 * Ez**2 R500Mpc = np.power((3 * M500) / (4 * np.pi * 500 * criticalDensity), 1.0 / 3.0) return R500Mpc @@ -231,9 +245,7 @@ def calcTheta500Arcmin(z, M500, Ez_fn, Da_fn, H0): # ---------------------------------------------------------------------------------------- def calcQ(theta500Arcmin, tck): - """Returns Q, given theta500Arcmin, and a set of spline fit knots for (theta, Q). - - """ + """Returns Q, given theta500Arcmin, and a set of spline fit knots for (theta, Q).""" # Q=np.poly1d(coeffs)(theta500Arcmin) Q = interpolate.splev(theta500Arcmin, tck) @@ -261,7 +273,7 @@ def calcFRel(z, M500, obsFreqGHz=148.0, Ez_fn=None): TKelvin = TkeV * ((1000 * elementary_charge) / k_Boltzmann) # Itoh et al. (1998) eqns. 2.25 - 2.30 - thetae = (k_Boltzmann * TKelvin) / (electron_mass_kg * C_M_S ** 2) + thetae = (k_Boltzmann * TKelvin) / (electron_mass_kg * C_M_S**2) X = (h_Planck * obsFreqGHz * 1e9) / (k_Boltzmann * T_CMB) Xtw = X * (np.cosh(X / 2.0) / np.sinh(X / 2.0)) Stw = X / np.sinh(X / 2.0) @@ -269,92 +281,103 @@ def calcFRel(z, M500, obsFreqGHz=148.0, Ez_fn=None): Y0 = -4 + Xtw Y1 = ( - -10.0 - + (47 / 2.0) * Xtw - - (42 / 5.0) * Xtw ** 2 - + (7 / 10.0) * Xtw ** 3 - + np.power(Stw, 2) * (-(21 / 5.0) + (7 / 5.0) * Xtw) + -10.0 + + (47 / 2.0) * Xtw + - (42 / 5.0) * Xtw**2 + + (7 / 10.0) * Xtw**3 + + np.power(Stw, 2) * (-(21 / 5.0) + (7 / 5.0) * Xtw) ) Y2 = ( - -(15 / 2.0) - + (1023 / 8.0) * Xtw - - (868 / 5.0) * Xtw ** 2 - + (329 / 5.0) * Xtw ** 3 - - (44 / 5.0) * Xtw ** 4 - + (11 / 30.0) * Xtw ** 5 - + np.power(Stw, 2) - * (-(434 / 5.0) + (658 / 5.0) * Xtw - - (242 / 5.0) * Xtw ** 2 - + (143 / 30.0) * Xtw ** 3) - + np.power(Stw, 4) * (-(44 / 5.0) + (187 / 60.0) * Xtw) + -(15 / 2.0) + + (1023 / 8.0) * Xtw + - (868 / 5.0) * Xtw**2 + + (329 / 5.0) * Xtw**3 + - (44 / 5.0) * Xtw**4 + + (11 / 30.0) * Xtw**5 + + np.power(Stw, 2) + * ( + -(434 / 5.0) + + (658 / 5.0) * Xtw + - (242 / 5.0) * Xtw**2 + + (143 / 30.0) * Xtw**3 + ) + + np.power(Stw, 4) * (-(44 / 5.0) + (187 / 60.0) * Xtw) ) Y3 = ( - (15 / 2.0) - + (2505 / 8.0) * Xtw - - (7098 / 5.0) * Xtw ** 2 - + (14253 / 10.0) * Xtw ** 3 - - (18594 / 35.0) * Xtw ** 4 - + (12059 / 140.0) * Xtw ** 5 - - (128 / 21.0) * Xtw ** 6 - + (16 / 105.0) * Xtw ** 7 - + np.power(Stw, 2) - * ( - -(7098 / 10.0) - + (14253 / 5.0) * Xtw - - (102267 / 35.0) * Xtw ** 2 - + (156767 / 140.0) * Xtw ** 3 - - (1216 / 7.0) * Xtw ** 4 - + (64 / 7.0) * Xtw ** 5 - ) - + np.power(Stw, 4) - * (-(18594 / 35.0) + (205003 / 280.0) * Xtw - - (1920 / 7.0) * Xtw ** 2 + (1024 / 35.0) * Xtw ** 3) - + np.power(Stw, 6) * (-(544 / 21.0) + (992 / 105.0) * Xtw) + (15 / 2.0) + + (2505 / 8.0) * Xtw + - (7098 / 5.0) * Xtw**2 + + (14253 / 10.0) * Xtw**3 + - (18594 / 35.0) * Xtw**4 + + (12059 / 140.0) * Xtw**5 + - (128 / 21.0) * Xtw**6 + + (16 / 105.0) * Xtw**7 + + np.power(Stw, 2) + * ( + -(7098 / 10.0) + + (14253 / 5.0) * Xtw + - (102267 / 35.0) * Xtw**2 + + (156767 / 140.0) * Xtw**3 + - (1216 / 7.0) * Xtw**4 + + (64 / 7.0) * Xtw**5 + ) + + np.power(Stw, 4) + * ( + -(18594 / 35.0) + + (205003 / 280.0) * Xtw + - (1920 / 7.0) * Xtw**2 + + (1024 / 35.0) * Xtw**3 + ) + + np.power(Stw, 6) * (-(544 / 21.0) + (992 / 105.0) * Xtw) ) Y4 = ( - -(135 / 32.0) - + (30375 / 128.0) * Xtw - - (62391 / 10.0) * Xtw ** 2 - + (614727 / 40.0) * Xtw ** 3 - - (124389 / 10.0) * Xtw ** 4 - + (355703 / 80.0) * Xtw ** 5 - - (16568 / 21.0) * Xtw ** 6 - + (7516 / 105.0) * Xtw ** 7 - - (22 / 7.0) * Xtw ** 8 - + (11 / 210.0) * Xtw ** 9 - + np.power(Stw, 2) - * ( - -(62391 / 20.0) - + (614727 / 20.0) * Xtw - - (1368279 / 20.0) * Xtw ** 2 - + (4624139 / 80.0) * Xtw ** 3 - - (157396 / 7.0) * Xtw ** 4 - + (30064 / 7.0) * Xtw ** 5 - - (2717 / 7.0) * Xtw ** 6 - + (2761 / 210.0) * Xtw ** 7 - ) - + np.power(Stw, 4) - * ( - -(124389 / 10.0) - + (6046951 / 160.0) * Xtw - - (248520 / 7.0) * Xtw ** 2 - + (481024 / 35.0) * Xtw ** 3 - - (15972 / 7.0) * Xtw ** 4 - + (18689 / 140.0) * Xtw ** 5 - ) - + np.power(Stw, 6) - * (-(70414 / 21.0) + (465992 / 105.0) * Xtw - - (11792 / 7.0) * Xtw ** 2 + (19778 / 105.0) * Xtw ** 3) - + np.power(Stw, 8) * (-(682 / 7.0) + (7601 / 210.0) * Xtw) + -(135 / 32.0) + + (30375 / 128.0) * Xtw + - (62391 / 10.0) * Xtw**2 + + (614727 / 40.0) * Xtw**3 + - (124389 / 10.0) * Xtw**4 + + (355703 / 80.0) * Xtw**5 + - (16568 / 21.0) * Xtw**6 + + (7516 / 105.0) * Xtw**7 + - (22 / 7.0) * Xtw**8 + + (11 / 210.0) * Xtw**9 + + np.power(Stw, 2) + * ( + -(62391 / 20.0) + + (614727 / 20.0) * Xtw + - (1368279 / 20.0) * Xtw**2 + + (4624139 / 80.0) * Xtw**3 + - (157396 / 7.0) * Xtw**4 + + (30064 / 7.0) * Xtw**5 + - (2717 / 7.0) * Xtw**6 + + (2761 / 210.0) * Xtw**7 + ) + + np.power(Stw, 4) + * ( + -(124389 / 10.0) + + (6046951 / 160.0) * Xtw + - (248520 / 7.0) * Xtw**2 + + (481024 / 35.0) * Xtw**3 + - (15972 / 7.0) * Xtw**4 + + (18689 / 140.0) * Xtw**5 + ) + + np.power(Stw, 6) + * ( + -(70414 / 21.0) + + (465992 / 105.0) * Xtw + - (11792 / 7.0) * Xtw**2 + + (19778 / 105.0) * Xtw**3 + ) + + np.power(Stw, 8) * (-(682 / 7.0) + (7601 / 210.0) * Xtw) ) deltaSZE = ( - ((X ** 3) / (np.exp(X) - 1)) - * ((thetae * X * np.exp(X)) / (np.exp(X) - 1)) - * (Y0 + Y1 * thetae + Y2 * thetae ** 2 + Y3 * thetae ** 3 + Y4 * thetae ** 4) + ((X**3) / (np.exp(X) - 1)) + * ((thetae * X * np.exp(X)) / (np.exp(X) - 1)) + * (Y0 + Y1 * thetae + Y2 * thetae**2 + Y3 * thetae**3 + Y4 * thetae**4) ) fRel = 1 + deltaSZE @@ -364,17 +387,17 @@ def calcFRel(z, M500, obsFreqGHz=148.0, Ez_fn=None): # ---------------------------------------------------------------------------------------- def y0FromLogM500( - log10M500, - z, - tckQFit, - tenToA0=4.95e-5, - B0=0.08, - Mpivot=3e14, - sigma_int=0.2, - fRelWeightsDict={148.0: 1.0}, - H0=70., - Ez_fn=None, - Da_fn=None + log10M500, + z, + tckQFit, + tenToA0=4.95e-5, + B0=0.08, + Mpivot=3e14, + sigma_int=0.2, + fRelWeightsDict={148.0: 1.0}, + H0=70.0, + Ez_fn=None, + Da_fn=None, ): """Predict y0~ given logM500 (in MSun) and redshift. Default scaling relation parameters are A10 (as in H13). diff --git a/soliket/clusters/tinker.py b/soliket/clusters/tinker.py index 55d160bb..872a0cbc 100644 --- a/soliket/clusters/tinker.py +++ b/soliket/clusters/tinker.py @@ -11,9 +11,10 @@ # Tinker stuff -tinker_data = np.transpose([[float(x) for x in line.split()] - for line in - """200 0.186 1.47 2.57 1.19 +tinker_data = np.transpose( + [ + [float(x) for x in line.split()] + for line in """200 0.186 1.47 2.57 1.19 300 0.200 1.52 2.25 1.27 400 0.212 1.56 2.05 1.34 600 0.218 1.61 1.87 1.45 @@ -21,7 +22,11 @@ 1200 0.255 2.13 1.51 1.80 1600 0.260 2.30 1.46 1.97 2400 0.260 2.53 1.44 2.24 - 3200 0.260 2.66 1.41 2.44""".split('\n')]) + 3200 0.260 2.66 1.41 2.44""".split( + "\n" + ) + ] +) tinker_splines = None @@ -34,7 +39,7 @@ def tinker_params_spline(delta, z=None): for y in data: # Extend to large Delta p = np.polyfit(D[-2:], y[-2:], 1) - x = np.hstack((D, D[-1] + 3.)) + x = np.hstack((D, D[-1] + 3.0)) y = np.hstack((y, np.polyval(p, x[-1]))) tinker_splines.append(iuSpline(x, y, k=2)) A0, a0, b0, c0 = [ts(np.log(delta)) for ts in tinker_splines] @@ -42,9 +47,9 @@ def tinker_params_spline(delta, z=None): return A0, a0, b0, c0 z = np.asarray(z) - A = A0 * (1 + z) ** -.14 - a = a0 * (1 + z) ** -.06 - alpha = 10. ** (-((.75 / np.log10(delta / 75.)) ** 1.2)) + A = A0 * (1 + z) ** -0.14 + a = a0 * (1 + z) ** -0.06 + alpha = 10.0 ** (-((0.75 / np.log10(delta / 75.0)) ** 1.2)) b = b0 * (1 + z) ** -alpha c = np.zeros(np.shape(z)) + c0 return A, a, b, c @@ -53,13 +58,12 @@ def tinker_params_spline(delta, z=None): def tinker_params_analytic(delta, z=None): alpha = None if np.asarray(delta).ndim == 0: # scalar delta. - A0, a0, b0, c0 = [p[0] for p in - tinker_params(np.array([delta]), z=None)] + A0, a0, b0, c0 = [p[0] for p in tinker_params(np.array([delta]), z=None)] if z is not None: - if delta < 75.: - alpha = 1. + if delta < 75.0: + alpha = 1.0 else: - alpha = 10. ** (-((.75 / np.log10(delta / 75.)) ** 1.2)) + alpha = 10.0 ** (-((0.75 / np.log10(delta / 75.0)) ** 1.2)) else: log_delta = np.log10(delta) @@ -67,18 +71,18 @@ def tinker_params_analytic(delta, z=None): a0 = 1.43 + (log_delta - 2.3) ** 1.5 b0 = 1.0 + (log_delta - 1.6) ** (-1.5) c0 = log_delta - 2.35 - A0[delta > 1600] = .26 + A0[delta > 1600] = 0.26 a0[log_delta < 2.3] = 1.43 b0[log_delta < 1.6] = 1.0 - c0[c0 < 0] = 0. - c0 = 1.2 + c0 ** 1.6 + c0[c0 < 0] = 0.0 + c0 = 1.2 + c0**1.6 if z is None: return A0, a0, b0, c0 - A = A0 * (1 + z) ** -.14 - a = a0 * (1 + z) ** -.06 + A = A0 * (1 + z) ** -0.14 + a = a0 * (1 + z) ** -0.06 if alpha is None: - alpha = 10. ** (-((.75 / np.log10(delta / 75.)) ** 1.2)) - alpha[delta < 75.] = 1. + alpha = 10.0 ** (-((0.75 / np.log10(delta / 75.0)) ** 1.2)) + alpha[delta < 75.0] = 1.0 b = b0 * (1 + z) ** -alpha c = np.zeros(np.shape(z)) + c0 return A, a, b, c @@ -89,16 +93,17 @@ def tinker_params_analytic(delta, z=None): def tinker_f(sigma, params): A, a, b, c = params - return A * ((sigma / b) ** -a + 1) * np.exp(-c / sigma ** 2) + return A * ((sigma / b) ** -a + 1) * np.exp(-c / sigma**2) # Sigma-evaluation, and top-hat functions. + def radius_from_mass(M, rho): """ Convert mass M to radius R assuming density rho. """ - return (3. * M / (4. * np.pi * rho)) ** (1 / 3.) + return (3.0 * M / (4.0 * np.pi * rho)) ** (1 / 3.0) def top_hatf(kR): @@ -110,7 +115,7 @@ def top_hatf(kR): * This is called many times and costs a lot of runtime. * For small values, use Taylor series. """ - out = np.nan_to_num(3 * (np.sin(kR) - kR * np.cos(kR))) / (kR ** 3) + out = np.nan_to_num(3 * (np.sin(kR) - kR * np.cos(kR))) / (kR**3) return out @@ -123,13 +128,18 @@ def sigma_sq_integral(R_grid, power_spt, k_val): smarter way using numpy arrays. """ to_integ = np.array( - [top_hatf(R_grid * k) ** 2 * np.tile( - power_spt[:, i], - (R_grid.shape[0], 1), - ) * k ** 2 for k, i in zip(k_val, np.arange(len(k_val)))] + [ + top_hatf(R_grid * k) ** 2 + * np.tile( + power_spt[:, i], + (R_grid.shape[0], 1), + ) + * k**2 + for k, i in zip(k_val, np.arange(len(k_val))) + ] ) - return simpson(to_integ / (2 * np.pi ** 2), x=k_val, axis=0) + return simpson(to_integ / (2 * np.pi**2), x=k_val, axis=0) def dn_dlogM(M, z, rho, delta, k, P, comoving=False): @@ -152,7 +162,7 @@ def dn_dlogM(M, z, rho, delta, k, P, comoving=False): if not comoving: # if you do this make sure rho still has shape of z. R = R * np.transpose(1 + z) # Fluctuations on those scales (P and k are comoving) - sigma = sigma_sq_integral(R, P, k) ** .5 + sigma = sigma_sq_integral(R, P, k) ** 0.5 # d log(sigma^-1) # gradient is broken. if R.shape[-1] == 1: diff --git a/soliket/constants.py b/soliket/constants.py index 8d5f3bab..45dd45c0 100644 --- a/soliket/constants.py +++ b/soliket/constants.py @@ -1,14 +1,14 @@ from scipy import constants C_M_S = constants.c -C_KM_S = constants.c * 1.e-3 -C_HMPC = constants.c * 1.e-5 +C_KM_S = constants.c * 1.0e-3 +C_HMPC = constants.c * 1.0e-5 h_Planck = constants.h k_Boltzmann = constants.k elementary_charge = constants.e electron_mass_kg = constants.m_e -T_CMB = 2.72548 # in K -MSUN_CGS = 1.98840987e+33 +T_CMB = 2.72548 # in K +MSUN_CGS = 1.98840987e33 G_CGS = constants.G * 1e3 MPC2CM = constants.parsec * 1e8 diff --git a/soliket/cosmopower/cosmopower.py b/soliket/cosmopower/cosmopower.py index a4735d71..d14ebbc2 100644 --- a/soliket/cosmopower/cosmopower.py +++ b/soliket/cosmopower/cosmopower.py @@ -124,18 +124,20 @@ def initialize(self) -> None: netpath = os.path.join(self.network_path, nettype["filename"]) if nettype["type"] == "NN": - network = cp.cosmopower_NN( - restore=True, restore_filename=netpath) + network = cp.cosmopower_NN(restore=True, restore_filename=netpath) elif nettype["type"] == "PCAplusNN": network = cp.cosmopower_PCAplusNN( - restore=True, restore_filename=netpath) + restore=True, restore_filename=netpath + ) elif self.stop_at_error: # pragma: no cover raise ValueError( - f"Unknown network type {nettype['type']} for network {spectype}.") + f"Unknown network type {nettype['type']} for network {spectype}." + ) else: # pragma: no cover self.log.warn( f"Unknown network type {nettype['type']}\ - for network {spectype}: skipped!") + for network {spectype}: skipped!" + ) netdata["type"] = nettype["type"] netdata["log"] = nettype.get("log", True) @@ -153,8 +155,7 @@ def initialize(self) -> None: self.extra_args["lmax"] = None self.log.info(f"Loaded CosmoPower from directory {self.network_path}") - self.log.info( - f"CosmoPower will expect the parameters {self.all_parameters}") + self.log.info(f"CosmoPower will expect the parameters {self.all_parameters}") def calculate(self, state: dict, want_derived: bool = True, **params) -> bool: ## sadly, this syntax not valid until python 3.9 @@ -163,22 +164,22 @@ def calculate(self, state: dict, want_derived: bool = True, **params) -> bool: # } | { # self.translate_param(p): [params[p]] for p in params # } - cmb_params = {**{ - p: [params[p]] for p in params - }, **{ - self.translate_param(p): [params[p]] for p in params - }} + cmb_params = { + **{p: [params[p]] for p in params}, + **{self.translate_param(p): [params[p]] for p in params}, + } ells = None for spectype in self.networks: network = self.networks[spectype] - used_params = {par: (cmb_params[par] if par in cmb_params else [ - params[par]]) for par in network["parameters"]} + used_params = { + par: (cmb_params[par] if par in cmb_params else [params[par]]) + for par in network["parameters"] + } if network["log"]: - data = network["network"].ten_to_predictions_np(used_params)[ - 0, :] + data = network["network"].ten_to_predictions_np(used_params)[0, :] else: data = network["network"].predictions_np(used_params)[0, :] @@ -210,12 +211,13 @@ def get_Cl(self, ell_factor: bool = False, units: str = "FIRASmuK2") -> dict: if ell_factor: prefac *= self.ell_factor(ls, k) - cls[k][ls] = cls_old[k] * prefac * \ - self.cmb_unit_factor(k, units, 2.7255) + cls[k][ls] = cls_old[k] * prefac * self.cmb_unit_factor(k, units, 2.7255) cls[k][:2] = 0.0 if np.any(np.isnan(cls[k])): - self.log.warning("CosmoPower used outside of trained " - "{} ell range. Filled in with NaNs.".format(k)) + self.log.warning( + "CosmoPower used outside of trained " + "{} ell range. Filled in with NaNs.".format(k) + ) return cls @@ -245,15 +247,15 @@ def ell_factor(self, ls: np.ndarray, spectra: str) -> np.ndarray: if spectra in ["tt", "te", "tb", "ee", "et", "eb", "bb", "bt", "be"]: ellfac = ls * (ls + 1.0) / (2.0 * np.pi) elif spectra in ["pt", "pe", "pb", "tp", "ep", "bp"]: - ellfac = (ls * (ls + 1.0)) ** (3. / 2.) / (2.0 * np.pi) + ellfac = (ls * (ls + 1.0)) ** (3.0 / 2.0) / (2.0 * np.pi) elif spectra in ["pp"]: ellfac = (ls * (ls + 1.0)) ** 2.0 / (2.0 * np.pi) return ellfac - def cmb_unit_factor(self, spectra: str, - units: str = "FIRASmuK2", - Tcmb: float = 2.7255) -> float: + def cmb_unit_factor( + self, spectra: str, units: str = "FIRASmuK2", Tcmb: float = 2.7255 + ) -> float: """ Calculate the CMB prefactor for going from dimensionless power spectra to CMB units. @@ -270,12 +272,12 @@ def cmb_unit_factor(self, spectra: str, if x == "t" or x == "e" or x == "b": res *= self._cmb_unit_factor(units, Tcmb) elif x == "p": - res *= 1. / np.sqrt(2.0 * np.pi) + res *= 1.0 / np.sqrt(2.0 * np.pi) if y == "t" or y == "e" or y == "b": res *= self._cmb_unit_factor(units, Tcmb) elif y == "p": - res *= 1. / np.sqrt(2.0 * np.pi) + res *= 1.0 / np.sqrt(2.0 * np.pi) return res @@ -308,26 +310,26 @@ def initialize(self) -> None: netpath = os.path.join(self.network_path, self.network_settings["filename"]) if self.network_settings["type"] == "NN": - self.network = cp.cosmopower_NN( - restore=True, restore_filename=netpath) + self.network = cp.cosmopower_NN(restore=True, restore_filename=netpath) elif self.network_settings["type"] == "PCAplusNN": self.network = cp.cosmopower_PCAplusNN( - restore=True, restore_filename=netpath) + restore=True, restore_filename=netpath + ) else: - raise LoggedError( - f"Unknown network type {self.network_settings['type']}.") + raise LoggedError(f"Unknown network type {self.network_settings['type']}.") self.input_parameters = set(self.network.parameters) self.log_data = self.network_settings.get("log", False) + self.log.info(f"Loaded CosmoPowerDerived from directory {self.network_path}") self.log.info( - f"Loaded CosmoPowerDerived from directory {self.network_path}") - self.log.info( - f"CosmoPowerDerived will expect the parameters {self.input_parameters}") + f"CosmoPowerDerived will expect the parameters {self.input_parameters}" + ) self.log.info( f"CosmoPowerDerived can provide the following parameters: \ - {self.get_can_provide()}.") + {self.get_can_provide()}." + ) def translate_param(self, p): return self.renames.get(p, p) @@ -339,11 +341,10 @@ def calculate(self, state: dict, want_derived: bool = True, **params) -> bool: # } | { # self.translate_param(p): [params[p]] for p in params # } - input_params = {**{ - p: [params[p]] for p in params - }, **{ - self.translate_param(p): [params[p]] for p in params - }} + input_params = { + **{p: [params[p]] for p in params}, + **{self.translate_param(p): [params[p]] for p in params}, + } if self.log_data: data = self.network.ten_to_predictions_np(input_params)[0, :] @@ -378,5 +379,10 @@ def get_requirements(self) -> Iterable[Tuple[str, str]]: return requirements def get_can_provide(self) -> Iterable[str]: - return set([par for par in self.derived_parameters - if (len(par) > 0 and not par == "_")]) + return set( + [ + par + for par in self.derived_parameters + if (len(par) > 0 and not par == "_") + ] + ) diff --git a/soliket/cross_correlation/__init__.py b/soliket/cross_correlation/__init__.py index 51085114..a1e56aff 100644 --- a/soliket/cross_correlation/__init__.py +++ b/soliket/cross_correlation/__init__.py @@ -1,5 +1,8 @@ -from .cross_correlation import (CrossCorrelationLikelihood, - GalaxyKappaLikelihood, ShearKappaLikelihood) +from .cross_correlation import ( + CrossCorrelationLikelihood, + GalaxyKappaLikelihood, + ShearKappaLikelihood, +) __all__ = [ "CrossCorrelationLikelihood", diff --git a/soliket/cross_correlation/cross_correlation.py b/soliket/cross_correlation/cross_correlation.py index b6cbccba..3c163471 100644 --- a/soliket/cross_correlation/cross_correlation.py +++ b/soliket/cross_correlation/cross_correlation.py @@ -6,6 +6,7 @@ """ import numpy as np + try: from numpy import trapezoid except ImportError: @@ -22,7 +23,6 @@ class CrossCorrelationLikelihood(GaussianLikelihood): """ def initialize(self): - self._get_sacc_data() self._check_tracers() @@ -34,32 +34,42 @@ def _get_CCL_results(self): return cosmo_dict["ccl"], cosmo_dict["cosmo"] 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 \ + 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__)) + 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 \ + 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)) + ".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': - bias = params_values['{}_deltaz'.format(tracer_name)] + if self.z_nuisance_mode == "deltaz": + bias = params_values["{}_deltaz".format(tracer_name)] nz_biased = tracer.get_dndz(z - bias) # nz_biased /= np.trapezoid(nz_biased, z) @@ -67,10 +77,9 @@ def _get_nz(self, z, tracer, tracer_name, **params_values): return nz_biased def _get_sacc_data(self, **params_values): - self.sacc_data = sacc.Sacc.load_fits(self.datapath) - if self.use_spectra == 'all': + if self.use_spectra == "all": pass else: for tracer_comb in self.sacc_data.get_tracer_combinations(): @@ -84,7 +93,6 @@ def _get_sacc_data(self, **params_values): self.data = GaussianData(self.name, self.x, self.y, self.cov, self.ncovsims) def _construct_ell_bins(self): - ell_eff = [] for tracer_comb in self.sacc_data.get_tracer_combinations(): @@ -114,7 +122,6 @@ 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 @@ -132,16 +139,14 @@ class GalaxyKappaLikelihood(CrossCorrelationLikelihood): r""" Likelihood for cross-correlations of galaxy and CMB lensing data. """ - _allowable_tracers = ['cmb_convergence', 'galaxy_density'] + _allowable_tracers = ["cmb_convergence", "galaxy_density"] def _get_theory(self, **params_values): - 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": @@ -151,17 +156,16 @@ def _get_theory(self, **params_values): 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')) + 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)) @@ -179,19 +183,17 @@ class ShearKappaLikelihood(CrossCorrelationLikelihood): _allowable_tracers = ["cmb_convergence", "galaxy_shear"] def _get_theory(self, **params_values): - ccl, cosmo = self._get_CCL_results() cl_binned_list = [] 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')) + 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 @@ -199,39 +201,38 @@ def _get_theory(self, **params_values): if self.ia_mode is None: ia_z = None - elif self.ia_mode == 'nla': - A_IA = params_values['A_IA'] - eta_IA = params_values['eta_IA'] + elif self.ia_mode == "nla": + A_IA = params_values["A_IA"] + eta_IA = params_values["eta_IA"] z0_IA = trapezoid(z_tracer1 * nz_tracer1) ia_z = (z_tracer1, A_IA * ((1 + z_tracer1) / (1 + z0_IA)) ** eta_IA) - elif self.ia_mode == 'nla-perbin': - A_IA = params_values['{}_A_IA'.format(sheartracer_name)] + elif self.ia_mode == "nla-perbin": + A_IA = params_values["{}_A_IA".format(sheartracer_name)] ia_z = (z_tracer1, A_IA * np.ones_like(z_tracer1)) - elif self.ia_mode == 'nla-noevo': - A_IA = params_values['A_IA'] + elif self.ia_mode == "nla-noevo": + A_IA = params_values["A_IA"] ia_z = (z_tracer1, A_IA * np.ones_like(z_tracer1)) - tracer1 = ccl.WeakLensingTracer(cosmo, - dndz=(z_tracer1, nz_tracer1), - ia_bias=ia_z) + 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) + 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) + 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')) + 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 @@ -239,32 +240,31 @@ def _get_theory(self, **params_values): if self.ia_mode is None: ia_z = None - elif self.ia_mode == 'nla': - A_IA = params_values['A_IA'] - eta_IA = params_values['eta_IA'] + elif self.ia_mode == "nla": + A_IA = params_values["A_IA"] + eta_IA = params_values["eta_IA"] z0_IA = trapezoid(z_tracer2 * nz_tracer2) ia_z = (z_tracer2, A_IA * ((1 + z_tracer2) / (1 + z0_IA)) ** eta_IA) - elif self.ia_mode == 'nla-perbin': - A_IA = params_values['{}_A_IA'.format(sheartracer_name)] + elif self.ia_mode == "nla-perbin": + A_IA = params_values["{}_A_IA".format(sheartracer_name)] ia_z = (z_tracer2, A_IA * np.ones_like(z_tracer2)) - elif self.ia_mode == 'nla-noevo': - A_IA = params_values['A_IA'] + elif self.ia_mode == "nla-noevo": + A_IA = params_values["A_IA"] ia_z = (z_tracer2, A_IA * np.ones_like(z_tracer2)) - tracer2 = ccl.WeakLensingTracer(cosmo, - dndz=(z_tracer2, nz_tracer2), - ia_bias=ia_z) + 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) + 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) + 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) @@ -278,7 +278,7 @@ def _get_theory(self, **params_values): # 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['{}_m'.format(sheartracer_name)] + m_bias = params_values["{}_m".format(sheartracer_name)] cl_unbinned = (1 + m_bias) * cl_unbinned cl_binned = np.dot(w_bins, cl_unbinned) diff --git a/soliket/foreground/foreground.py b/soliket/foreground/foreground.py index bb9fbac7..21095896 100644 --- a/soliket/foreground/foreground.py +++ b/soliket/foreground/foreground.py @@ -73,9 +73,22 @@ def initialize(self): from fgspectra import frequency as fgf from fgspectra import power as fgp - self.expected_params_fg = ["a_tSZ", "a_kSZ", "a_p", "beta_p", - "a_c", "beta_c", "a_s", "a_gtt", "a_gte", "a_gee", - "a_psee", "a_pste", "xi", "T_d"] + self.expected_params_fg = [ + "a_tSZ", + "a_kSZ", + "a_p", + "beta_p", + "a_c", + "beta_c", + "a_s", + "a_gtt", + "a_gte", + "a_gee", + "a_psee", + "a_pste", + "xi", + "T_d", + ] self.requested_cls = self.spectra["polarizations"] self.lmin = self.spectra["lmin"] @@ -86,18 +99,22 @@ def initialize(self): if hasattr(self.eff_freqs, "__len__"): if not len(self.exp_ch) == len(self.eff_freqs): raise LoggedError( - self.log, "list of effective frequency has to have" - "same length as list of channels!" + self.log, + "list of effective frequency has to have" + "same length as list of channels!", ) # self.bands to be filled with passbands read from sacc file # if mflike is used - self.bands = {f"{expc}_s0": {'nu': [self.eff_freqs[iexpc]], 'bandpass': [1.]} - for iexpc, expc in enumerate(self.exp_ch)} + self.bands = { + f"{expc}_s0": {"nu": [self.eff_freqs[iexpc]], "bandpass": [1.0]} + for iexpc, expc in enumerate(self.exp_ch) + } - template_path = os.path.join(os.path.dirname(os.path.abspath(fgp.__file__)), - 'data') - cibc_file = os.path.join(template_path, 'cl_cib_Choi2020.dat') + template_path = os.path.join( + os.path.dirname(os.path.abspath(fgp.__file__)), "data" + ) + cibc_file = os.path.join(template_path, "cl_cib_Choi2020.dat") # set pivot freq and multipole self.fg_nu_0 = self.foregrounds["normalisation"]["nu_0"] @@ -109,8 +126,9 @@ def initialize(self): self.cibp = fgc.FactorizedCrossSpectrum(fgf.ModifiedBlackBody(), fgp.PowerLaw()) self.radio = fgc.FactorizedCrossSpectrum(fgf.PowerLaw(), fgp.PowerLaw()) self.tsz = fgc.FactorizedCrossSpectrum(fgf.ThermalSZ(), fgp.tSZ_150_bat()) - self.cibc = fgc.FactorizedCrossSpectrum(fgf.CIB(), - fgp.PowerSpectrumFromFile(cibc_file)) + self.cibc = fgc.FactorizedCrossSpectrum( + fgf.CIB(), fgp.PowerSpectrumFromFile(cibc_file) + ) self.dust = fgc.FactorizedCrossSpectrum(fgf.ModifiedBlackBody(), fgp.PowerLaw()) self.tSZ_and_CIB = fgc.SZxCIB_Choi2020() @@ -119,21 +137,26 @@ def initialize(self): def initialize_with_params(self): # Check that the parameters are the right ones differences = are_different_params_lists( - self.input_params, self.expected_params_fg, - name_A="given", name_B="expected") + self.input_params, + self.expected_params_fg, + name_A="given", + name_B="expected", + ) if differences: raise LoggedError( - self.log, "Configuration error in parameters: %r.", - differences) + self.log, "Configuration error in parameters: %r.", differences + ) # Gets the actual power spectrum of foregrounds given the passed parameters - def _get_foreground_model(self, - requested_cls=None, - ell=None, - exp_ch=None, - bandint_freqs=None, - eff_freqs=None, - **fg_params): + def _get_foreground_model( + self, + requested_cls=None, + ell=None, + exp_ch=None, + bandint_freqs=None, + eff_freqs=None, + **fg_params, + ): r""" Gets the foreground power spectra for each component computed by ``fgspectra``. The computation assumes the bandpass transmissions from the ``BandPass`` class @@ -166,26 +189,26 @@ def _get_foreground_model(self, requested_cls = self.requested_cls # if ell = None, it uses ell from yaml, otherwise the ell array provided # useful to make tests at different l_max than the data - if not hasattr(ell, '__len__'): + if not hasattr(ell, "__len__"): ell = self.ell ell_0 = self.fg_ell_0 nu_0 = self.fg_nu_0 # Normalisation of radio sources - ell_clp = ell * (ell + 1.) - ell_0clp = ell_0 * (ell_0 + 1.) + ell_clp = ell * (ell + 1.0) + ell_0clp = ell_0 * (ell_0 + 1.0) # Set component spectra self.fg_component_list = {s: self.components[s] for s in requested_cls} # Set exp_ch list - if not hasattr(exp_ch, '__len__'): + if not hasattr(exp_ch, "__len__"): exp_ch = self.exp_ch # Set array of freqs to use if bandint_freqs is None - if not hasattr(bandint_freqs, '__len__'): - if hasattr(eff_freqs, '__len__'): + if not hasattr(bandint_freqs, "__len__"): + if hasattr(eff_freqs, "__len__"): bandint_freqs = np.asarray(eff_freqs) else: raise LoggedError( @@ -193,88 +216,89 @@ def _get_foreground_model(self, ) model = {} - model["tt", "kSZ"] = fg_params["a_kSZ"] * self.ksz({"nu": bandint_freqs}, - {"ell": ell, - "ell_0": ell_0}) - - model["tt", "cibp"] = fg_params["a_p"] * self.cibp({"nu": bandint_freqs, - "nu_0": nu_0, - "temp": fg_params["T_d"], - "beta": fg_params["beta_p"]}, - {"ell": ell_clp, - "ell_0": ell_0clp, - "alpha": 1}) - - model["tt", "radio"] = fg_params["a_s"] * self.radio({"nu": bandint_freqs, - "nu_0": nu_0, - "beta": -0.5 - 2.}, - {"ell": ell_clp, - "ell_0": ell_0clp, - "alpha": 1}) - - model["tt", "tSZ"] = fg_params["a_tSZ"] * self.tsz({"nu": bandint_freqs, - "nu_0": nu_0}, - {"ell": ell, - "ell_0": ell_0}) - - model["tt", "cibc"] = fg_params["a_c"] * self.cibc({"nu": bandint_freqs, - "nu_0": nu_0, - "temp": fg_params["T_d"], - "beta": fg_params["beta_c"]}, - {'ell': ell, - 'ell_0': ell_0}) - - model["tt", "dust"] = fg_params["a_gtt"] * self.dust({"nu": bandint_freqs, - "nu_0": nu_0, - "temp": 19.6, - "beta": 1.5}, - {"ell": ell, - "ell_0": 500., - "alpha": -0.6}) - - model["tt", "tSZ_and_CIB"] = \ - self.tSZ_and_CIB({'kwseq': ({'nu': bandint_freqs, 'nu_0': nu_0}, - {'nu': bandint_freqs, 'nu_0': nu_0, - 'temp': fg_params['T_d'], - 'beta': fg_params["beta_c"]})}, - {'kwseq': ({'ell': ell, 'ell_0': ell_0, - 'amp': fg_params['a_tSZ']}, - {'ell': ell, 'ell_0': ell_0, - 'amp': fg_params['a_c']}, - {'ell': ell, 'ell_0': ell_0, - 'amp': - fg_params['xi'] \ - * np.sqrt(fg_params['a_tSZ'] * - fg_params['a_c'])})}) - - model["ee", "radio"] = fg_params["a_psee"] * self.radio({"nu": bandint_freqs, - "nu_0": nu_0, - "beta": -0.5 - 2.}, - {"ell": ell_clp, - "ell_0": ell_0clp, - "alpha": 1}) - - model["ee", "dust"] = fg_params["a_gee"] * self.dust({"nu": bandint_freqs, - "nu_0": nu_0, - "temp": 19.6, - "beta": 1.5}, - {"ell": ell, - "ell_0": 500., - "alpha": -0.4}) - - model["te", "radio"] = fg_params["a_pste"] * self.radio({"nu": bandint_freqs, - "nu_0": nu_0, - "beta": -0.5 - 2.}, - {"ell": ell_clp, - "ell_0": ell_0clp, - "alpha": 1}) - - model["te", "dust"] = fg_params["a_gte"] * self.dust({"nu": bandint_freqs, - "nu_0": nu_0, - "temp": 19.6, - "beta": 1.5}, - {"ell": ell, - "ell_0": 500., - "alpha": -0.4}) + model["tt", "kSZ"] = fg_params["a_kSZ"] * self.ksz( + {"nu": bandint_freqs}, {"ell": ell, "ell_0": ell_0} + ) + + model["tt", "cibp"] = fg_params["a_p"] * self.cibp( + { + "nu": bandint_freqs, + "nu_0": nu_0, + "temp": fg_params["T_d"], + "beta": fg_params["beta_p"], + }, + {"ell": ell_clp, "ell_0": ell_0clp, "alpha": 1}, + ) + + model["tt", "radio"] = fg_params["a_s"] * self.radio( + {"nu": bandint_freqs, "nu_0": nu_0, "beta": -0.5 - 2.0}, + {"ell": ell_clp, "ell_0": ell_0clp, "alpha": 1}, + ) + + model["tt", "tSZ"] = fg_params["a_tSZ"] * self.tsz( + {"nu": bandint_freqs, "nu_0": nu_0}, {"ell": ell, "ell_0": ell_0} + ) + + model["tt", "cibc"] = fg_params["a_c"] * self.cibc( + { + "nu": bandint_freqs, + "nu_0": nu_0, + "temp": fg_params["T_d"], + "beta": fg_params["beta_c"], + }, + {"ell": ell, "ell_0": ell_0}, + ) + + model["tt", "dust"] = fg_params["a_gtt"] * self.dust( + {"nu": bandint_freqs, "nu_0": nu_0, "temp": 19.6, "beta": 1.5}, + {"ell": ell, "ell_0": 500.0, "alpha": -0.6}, + ) + + model["tt", "tSZ_and_CIB"] = self.tSZ_and_CIB( + { + "kwseq": ( + {"nu": bandint_freqs, "nu_0": nu_0}, + { + "nu": bandint_freqs, + "nu_0": nu_0, + "temp": fg_params["T_d"], + "beta": fg_params["beta_c"], + }, + ) + }, + { + "kwseq": ( + {"ell": ell, "ell_0": ell_0, "amp": fg_params["a_tSZ"]}, + {"ell": ell, "ell_0": ell_0, "amp": fg_params["a_c"]}, + { + "ell": ell, + "ell_0": ell_0, + "amp": -fg_params["xi"] + * np.sqrt(fg_params["a_tSZ"] * fg_params["a_c"]), + }, + ) + }, + ) + + model["ee", "radio"] = fg_params["a_psee"] * self.radio( + {"nu": bandint_freqs, "nu_0": nu_0, "beta": -0.5 - 2.0}, + {"ell": ell_clp, "ell_0": ell_0clp, "alpha": 1}, + ) + + model["ee", "dust"] = fg_params["a_gee"] * self.dust( + {"nu": bandint_freqs, "nu_0": nu_0, "temp": 19.6, "beta": 1.5}, + {"ell": ell, "ell_0": 500.0, "alpha": -0.4}, + ) + + model["te", "radio"] = fg_params["a_pste"] * self.radio( + {"nu": bandint_freqs, "nu_0": nu_0, "beta": -0.5 - 2.0}, + {"ell": ell_clp, "ell_0": ell_0clp, "alpha": 1}, + ) + + model["te", "dust"] = fg_params["a_gte"] * self.dust( + {"nu": bandint_freqs, "nu_0": nu_0, "temp": 19.6, "beta": 1.5}, + {"ell": ell, "ell_0": 500.0, "alpha": -0.4}, + ) fg_dict = {} for c1, f1 in enumerate(exp_ch): @@ -286,9 +310,9 @@ def _get_foreground_model(self, fg_dict[s, "tSZ", f1, f2] = model[s, "tSZ"][c1, c2] fg_dict[s, "cibc", f1, f2] = model[s, "cibc"][c1, c2] fg_dict[s, "tSZxCIB", f1, f2] = ( - model[s, comp][c1, c2] - - model[s, "tSZ"][c1, c2] - - model[s, "cibc"][c1, c2] + model[s, comp][c1, c2] + - model[s, "tSZ"][c1, c2] + - model[s, "cibc"][c1, c2] ) fg_dict[s, "all", f1, f2] += model[s, comp][c1, c2] else: @@ -333,8 +357,14 @@ def calculate(self, state, want_derived=False, **params_values_dict): # compute bandpasses at each step only if bandint_shift params are not null # and bandint_freqs has been computed at least once if np.all( - np.array([params_values_dict[k] for k in params_values_dict.keys() - if "bandint_shift_" in k]) == 0.0 + np.array( + [ + params_values_dict[k] + for k in params_values_dict.keys() + if "bandint_shift_" in k + ] + ) + == 0.0 ): if not hasattr(self, "bandint_freqs"): self.log.info("Computing bandpass at first step, no shifts") @@ -343,11 +373,13 @@ def calculate(self, state, want_derived=False, **params_values_dict): self.bandint_freqs = self.get_bandpasses(**params_values_dict) fg_params = {k: params_values_dict[k] for k in self.expected_params_fg} - state["fg_dict"] = self._get_foreground_model(requested_cls=self.requested_cls, - ell=self.ell, - exp_ch=self.exp_ch, - bandint_freqs=self.bandint_freqs, - **fg_params) + state["fg_dict"] = self._get_foreground_model( + requested_cls=self.requested_cls, + ell=self.ell, + exp_ch=self.exp_ch, + bandint_freqs=self.bandint_freqs, + **fg_params, + ) def get_fg_dict(self): """ diff --git a/soliket/gaussian/gaussian.py b/soliket/gaussian/gaussian.py index 012f2ce1..fbf363b2 100644 --- a/soliket/gaussian/gaussian.py +++ b/soliket/gaussian/gaussian.py @@ -55,10 +55,10 @@ class MultiGaussianLikelihood(GaussianLikelihood): cross_cov_path: Optional[str] = None def __init__(self, info=empty_dict, **kwargs): - - if 'components' in info: - self.likelihoods = [get_likelihood(*kv) for kv in zip(info['components'], - info['options'])] + if "components" in info: + self.likelihoods = [ + get_likelihood(*kv) for kv in zip(info["components"], info["options"]) + ] default_info = merge_info(*[like.get_defaults() for like in self.likelihoods]) default_info.update(info) @@ -71,7 +71,7 @@ def initialize(self): data_list = [like.data for like in self.likelihoods] self.data = MultiGaussianData(data_list, self.cross_cov) - self.log.info('Initialized.') + self.log.info("Initialized.") def initialize_with_provider(self, provider): # pragma: no cover for like in self.likelihoods: @@ -88,8 +88,7 @@ def get_helper_theories(self): # pragma: no cover def _get_theory(self, **kwargs): return np.concatenate([like._get_theory(**kwargs) for like in self.likelihoods]) - def get_requirements(self): # pragma: no cover - + def get_requirements(self): # pragma: no cover # Reqs with arguments like 'lmax', etc. may have to be carefully treated here to # merge reqs = {} diff --git a/soliket/gaussian/gaussian_data.py b/soliket/gaussian/gaussian_data.py index aa98b99b..0913310f 100644 --- a/soliket/gaussian/gaussian_data.py +++ b/soliket/gaussian/gaussian_data.py @@ -5,8 +5,9 @@ class GaussianData: """ - Named multivariate gaussian data + Named multivariate gaussian data """ + name: str # name identifier for the data x: Sequence # labels for each data point y: np.ndarray # data point values @@ -16,15 +17,22 @@ class GaussianData: _fast_chi_squared = _fast_chi_square() - def __init__(self, name, x: Sequence, y: Sequence[float], cov: np.ndarray, - ncovsims: Optional[int] = None): - + def __init__( + self, + name, + x: Sequence, + y: Sequence[float], + cov: np.ndarray, + ncovsims: Optional[int] = None, + ): self.name = str(name) self.ncovsims = ncovsims if not (len(x) == len(y) and cov.shape == (len(x), len(x))): - raise ValueError(f"Incompatible shapes! x={len(x)}, y={len(y)}, \ - cov={cov.shape}") + raise ValueError( + f"Incompatible shapes! x={len(x)}, y={len(y)}, \ + cov={cov.shape}" + ) self.x = x self.y = np.ascontiguousarray(y) @@ -61,7 +69,6 @@ class MultiGaussianData(GaussianData): """ def __init__(self, data_list, cross_covs=None): - if cross_covs is None: cross_covs = {} @@ -120,8 +127,11 @@ def norm_const(self): @property def labels(self): - return [x for y in [[name] * len(d) for - name, d in zip(self.names, self.data_list)] for x in y] + return [ + x + for y in [[name] * len(d) for name, d in zip(self.names, self.data_list)] + for x in y + ] def _index_range(self, name): if name not in self.names: @@ -163,5 +173,6 @@ def plot_cov(self, **kwargs): for j, lj in zip(range(len(self.data)), self.labels) ] - return hv.HeatMap(data).opts(tools=["hover"], width=800, height=800, - invert_yaxis=True, xrotation=90) + return hv.HeatMap(data).opts( + tools=["hover"], width=800, height=800, invert_yaxis=True, xrotation=90 + ) diff --git a/soliket/halo_model/__init__.py b/soliket/halo_model/__init__.py index c54ba9bb..0a66250e 100644 --- a/soliket/halo_model/__init__.py +++ b/soliket/halo_model/__init__.py @@ -1,6 +1,3 @@ from .halo_model import HaloModel, HaloModel_pyhm -__all__ = [ - "HaloModel", - "HaloModel_pyhm" -] +__all__ = ["HaloModel", "HaloModel_pyhm"] diff --git a/soliket/halo_model/halo_model.py b/soliket/halo_model/halo_model.py index 6cb4719a..eca958ac 100644 --- a/soliket/halo_model/halo_model.py +++ b/soliket/halo_model/halo_model.py @@ -35,6 +35,7 @@ import numpy as np import pyhalomodel as halo from cobaya.theory import Theory + # from cobaya.theories.cosmo.boltzmannbase import PowerSpectrumInterpolator from scipy.interpolate import RectBivariateSpline @@ -55,28 +56,26 @@ def initialize(self): def _get_Pk_mm_lin(self): for pair in self._var_pairs: - self.k, self.z, Pk_mm = \ - self.provider.get_Pk_grid(var_pair=pair, nonlinear=False) + self.k, self.z, Pk_mm = self.provider.get_Pk_grid( + var_pair=pair, nonlinear=False + ) return Pk_mm def get_Pk_mm_grid(self): - return self.current_state["Pk_mm_grid"] def get_Pk_gg_grid(self): - return self.current_state["Pk_gg_grid"] def get_Pk_gm_grid(self): - return self.current_state["Pk_gm_grid"] class HaloModel_pyhm(HaloModel): """Halo Model wrapping the simple pyhalomodel code of Asgari, Mead & Heymans (2023) - We include this simple halo model for the non-linear matter-matter power spectrum with + We include this simple halo model for the non-linear matter-matter power spectrum with NFW profiles. This is calculated via the `pyhalomodel `_ code. """ @@ -86,42 +85,46 @@ def initialize(self): self.Ms = np.logspace(np.log10(self.Mmin), np.log10(self.Mmax), self.nM) def get_requirements(self): - return {"omegam": None} def must_provide(self, **requirements): - options = requirements.get("halo_model") or {} self._var_pairs.update( - set((x, y) for x, y in - options.get("vars_pairs", [("delta_tot", "delta_tot")]))) + set( + (x, y) + for x, y in options.get("vars_pairs", [("delta_tot", "delta_tot")]) + ) + ) self.kmax = max(self.kmax, options.get("kmax", self.kmax)) - self.z = np.unique(np.concatenate( - (np.atleast_1d(options.get("z", self._default_z_sampling)), - np.atleast_1d(self.z)))) + self.z = np.unique( + np.concatenate( + ( + np.atleast_1d(options.get("z", self._default_z_sampling)), + np.atleast_1d(self.z), + ) + ) + ) needs = {} needs["Pk_grid"] = { - "vars_pairs": self._var_pairs, - "nonlinear": (False, False), - "z": self.z, - "k_max": self.kmax - } - - needs["sigma_R"] = {"vars_pairs": self._var_pairs, - "z": self.z, - "k_max": self.kmax, - "R": np.linspace(0.14, 66, 256) # list of radii required - } - + "vars_pairs": self._var_pairs, + "nonlinear": (False, False), + "z": self.z, + "k_max": self.kmax, + } + + needs["sigma_R"] = { + "vars_pairs": self._var_pairs, + "z": self.z, + "k_max": self.kmax, + "R": np.linspace(0.14, 66, 256), # list of radii required + } return needs - def calculate(self, state: dict, want_derived: bool = True, - **params_values_dict): - + def calculate(self, state: dict, want_derived: bool = True, **params_values_dict): Pk_mm_lin = self._get_Pk_mm_lin() # now wish to interpolate sigma_R to these Rs @@ -133,32 +136,45 @@ def calculate(self, state: dict, want_derived: bool = True, # for sure we could avoid the for loop with some thought for iz, zeval in enumerate(self.z): - hmod = halo.model(zeval, self.provider.get_param('omegam'), - name=self.hmf_name, Dv=self.hmf_Dv) + hmod = halo.model( + zeval, + self.provider.get_param("omegam"), + name=self.hmf_name, + Dv=self.hmf_Dv, + ) Rs = hmod.Lagrangian_radius(self.Ms) rvs = hmod.virial_radius(self.Ms) - cs = 7.85 * (self.Ms / 2e12)**-0.081 * (1. + zeval)**-0.71 + cs = 7.85 * (self.Ms / 2e12) ** -0.081 * (1.0 + zeval) ** -0.71 Uk = self._win_NFW(self.k, rvs, cs) - matter_profile = halo.profile.Fourier(self.k, self.Ms, Uk, - amplitude=self.Ms, - normalisation=hmod.rhom, - mass_tracer=True) - - Pk_2h, Pk_1h, Pk_hm = hmod.power_spectrum(self.k, Pk_mm_lin[iz], - self.Ms, sigmaRs(zeval, Rs)[0], - {'m': matter_profile}, - verbose=False) - - output_Pk_hm_mm[iz] = Pk_hm['m-m'] - - state['Pk_mm_grid'] = output_Pk_hm_mm + matter_profile = halo.profile.Fourier( + self.k, + self.Ms, + Uk, + amplitude=self.Ms, + normalisation=hmod.rhom, + mass_tracer=True, + ) + + Pk_2h, Pk_1h, Pk_hm = hmod.power_spectrum( + self.k, + Pk_mm_lin[iz], + self.Ms, + sigmaRs(zeval, Rs)[0], + {"m": matter_profile}, + verbose=False, + ) + + output_Pk_hm_mm[iz] = Pk_hm["m-m"] + + state["Pk_mm_grid"] = output_Pk_hm_mm # state['Pk_gm_grid'] = Pk_hm['g-m'] # state['Pk_gg_grid'] = Pk_hm['g-g'] def _win_NFW(self, k, rv, c): from scipy.special import sici + rs = rv / c kv = np.outer(k, rv) ks = np.outer(k, rs) @@ -167,6 +183,6 @@ def _win_NFW(self, k, rv, c): f1 = np.cos(ks) * (Cisv - Cis) f2 = np.sin(ks) * (Sisv - Sis) f3 = np.sin(kv) / (ks + kv) - f4 = np.log(1. + c) - c / (1. + c) + f4 = np.log(1.0 + c) - c / (1.0 + c) Wk = (f1 + f2 - f3) / f4 return Wk diff --git a/soliket/lensing/lensing.py b/soliket/lensing/lensing.py index e17eeac3..602f76b7 100644 --- a/soliket/lensing/lensing.py +++ b/soliket/lensing/lensing.py @@ -69,13 +69,14 @@ class LensingLikelihood(BinnedPSLikelihood, InstallableLikelihood): def initialize(self): self.log.info("Initialising.") # Set path to data - if ((not getattr(self, "path", None)) and - (not getattr(self, "packages_path", None))): + if (not getattr(self, "path", None)) and ( + not getattr(self, "packages_path", None) + ): raise LoggedError( self.log, "No path given to LensingLikelihood data. " "Set the likelihood property " - "'path' or 'packages_path'" + "'path' or 'packages_path'", ) # If no path specified, use the modules path @@ -108,11 +109,11 @@ def initialize(self): # Set the fiducial spectra self.ls = np.arange(0, self.lmax, dtype=np.longlong) - self.fcltt = Cls["tt"][0: self.lmax] - self.fclpp = Cls["pp"][0: self.lmax] - self.fclee = Cls["ee"][0: self.lmax] - self.fclte = Cls["te"][0: self.lmax] - self.fclbb = Cls["bb"][0: self.lmax] + self.fcltt = Cls["tt"][0:self.lmax] + self.fclpp = Cls["pp"][0:self.lmax] + self.fclee = Cls["ee"][0:self.lmax] + self.fclte = Cls["te"][0:self.lmax] + self.fclbb = Cls["bb"][0:self.lmax] self.thetaclkk = self.fclpp * (self.ls * (self.ls + 1)) ** 2 * 0.25 # load the correction terms generate from the script n1so.py @@ -173,9 +174,8 @@ def get_requirements(self): "ee": self.theory_lmax, "bb": self.theory_lmax, }, - "CCL": {"kmax": 10, - "nonlinear": True}, - "zstar": None + "CCL": {"kmax": 10, "nonlinear": True}, + "zstar": None, } def _get_CCL_results(self): @@ -183,22 +183,24 @@ def _get_CCL_results(self): return cosmo_dict["ccl"], cosmo_dict["cosmo"] def _get_data(self): - bin_centers, bandpowers, cov = \ - self.sacc.get_ell_cl(None, 'ck', 'ck', return_cov=True) + bin_centers, bandpowers, cov = self.sacc.get_ell_cl( + None, "ck", "ck", return_cov=True + ) self.x = bin_centers self.y = bandpowers return bin_centers, self.y def _get_cov(self): - bin_centers, bandpowers, cov = \ - self.sacc.get_ell_cl(None, 'ck', 'ck', return_cov=True) + bin_centers, bandpowers, cov = self.sacc.get_ell_cl( + None, "ck", "ck", return_cov=True + ) self.cov = cov return cov def _get_binning_matrix(self): - - bin_centers, bandpowers, cov, ind = \ - self.sacc.get_ell_cl(None, 'ck', 'ck', return_cov=True, return_ind=True) + bin_centers, bandpowers, cov, ind = self.sacc.get_ell_cl( + None, "ck", "ck", return_cov=True, return_ind=True + ) bpw = self.sacc.get_bandpower_windows(ind) binning_matrix = bpw.weight.T self.binning_matrix = binning_matrix @@ -215,7 +217,7 @@ def _get_theory(self, **params_values): cl = self.provider.get_Cl(ell_factor=False) if self.pp_ccl is False: - Cl_theo = cl["pp"][0: self.lmax] + Cl_theo = cl["pp"][0:self.lmax] ls = self.ls Clkk_theo = (ls * (ls + 1)) ** 2 * Cl_theo * 0.25 else: @@ -224,27 +226,27 @@ def _get_theory(self, **params_values): cmbk = ccl.CMBLensingTracer(cosmo, z_source=zstar) Clkk_theo = ccl.angular_cl(cosmo, cmbk, cmbk, self.ls) - Cl_tt = cl["tt"][0: self.lmax] - Cl_ee = cl["ee"][0: self.lmax] - Cl_te = cl["te"][0: self.lmax] - Cl_bb = cl["bb"][0: self.lmax] + Cl_tt = cl["tt"][0:self.lmax] + Cl_ee = cl["ee"][0:self.lmax] + Cl_te = cl["te"][0:self.lmax] + Cl_bb = cl["bb"][0:self.lmax] Clkk_binned = self.binning_matrix.dot(Clkk_theo) correction = ( - 2 - * (self.thetaclkk / self.n0) - * ( - np.dot(self.N0cltt, Cl_tt - self.fcltt) - + np.dot(self.N0clee, Cl_ee - self.fclee) - + np.dot(self.N0clbb, Cl_bb - self.fclbb) - + np.dot(self.N0clte, Cl_te - self.fclte) - ) - + np.dot(self.N1clpp, Clkk_theo - self.thetaclkk) - + np.dot(self.N1cltt, Cl_tt - self.fcltt) - + np.dot(self.N1clee, Cl_ee - self.fclee) - + np.dot(self.N1clbb, Cl_bb - self.fclbb) - + np.dot(self.N1clte, Cl_te - self.fclte) + 2 + * (self.thetaclkk / self.n0) + * ( + np.dot(self.N0cltt, Cl_tt - self.fcltt) + + np.dot(self.N0clee, Cl_ee - self.fclee) + + np.dot(self.N0clbb, Cl_bb - self.fclbb) + + np.dot(self.N0clte, Cl_te - self.fclte) + ) + + np.dot(self.N1clpp, Clkk_theo - self.thetaclkk) + + np.dot(self.N1cltt, Cl_tt - self.fcltt) + + np.dot(self.N1clee, Cl_ee - self.fclee) + + np.dot(self.N1clbb, Cl_bb - self.fclbb) + + np.dot(self.N1clte, Cl_te - self.fclte) ) # put the correction term into bandpowers @@ -260,14 +262,16 @@ class LensingLiteLikelihood(BinnedPSLikelihood): data. Simply a Gaussian likelihood between a provided binned ``pp`` data vector and covariance matrix, and the appropriate theory vector. """ + kind: str = "pp" lmax: int = 3000 def initialize(self): - data = os.path.join(self.get_class_path(), 'data') - self.datapath = self.datapath or os.path.join(data, 'binnedauto.txt') - self.covpath = self.covpath or os.path.join(data, 'binnedcov.txt') + data = os.path.join(self.get_class_path(), "data") + self.datapath = self.datapath or os.path.join(data, "binnedauto.txt") + self.covpath = self.covpath or os.path.join(data, "binnedcov.txt") - self.binning_matrix_path = self.binning_matrix_path or \ - os.path.join(data, 'binningmatrix.txt') + self.binning_matrix_path = self.binning_matrix_path or os.path.join( + data, "binningmatrix.txt" + ) super().initialize() diff --git a/soliket/mflike/mflike.py b/soliket/mflike/mflike.py index 3993c729..b12689d4 100644 --- a/soliket/mflike/mflike.py +++ b/soliket/mflike/mflike.py @@ -51,17 +51,19 @@ def initialize(self): self.spec_meta = [] # Set path to data - if ((not getattr(self, "path", None)) and - (not getattr(self, "packages_path", None))): - raise LoggedError(self.log, - "No path given to MFLike data. " - "Set the likelihood property " - "'path' or 'packages_path'" - ) + if (not getattr(self, "path", None)) and ( + not getattr(self, "packages_path", None) + ): + raise LoggedError( + self.log, + "No path given to MFLike data. " + "Set the likelihood property " + "'path' or 'packages_path'", + ) # If no path specified, use the modules path - data_file_path = os.path.normpath(getattr(self, "path", None) or - os.path.join(self.packages_path, - "data")) + data_file_path = os.path.normpath( + getattr(self, "path", None) or os.path.join(self.packages_path, "data") + ) self.data_folder = os.path.join(data_file_path, self.data_folder) if not os.path.exists(self.data_folder): @@ -99,11 +101,14 @@ def get_requirements(self): # mflike requires cmbfg_dict from theoryforge # cmbfg_dict requires some params to be computed reqs = dict() - reqs["cmbfg_dict"] = {"ell": self.l_bpws, - "requested_cls": self.requested_cls, - "lcuts": self.lcuts, - "exp_ch": self.experiments, - "bands": self.bands} + reqs["cmbfg_dict"] = { + "ell": np.arange(2, 2001), + # "ell": self.l_bpws, + "requested_cls": self.requested_cls, + "lcuts": self.lcuts, + "exp_ch": self.experiments, + "bands": self.bands, + } return reqs def _get_theory(self, **params_values): @@ -129,7 +134,8 @@ def loglike(self, cmbfg_dict): logp += self.logp_const self.log.debug( "Log-likelihood value computed " - "= {} (Χ² = {})".format(logp, -2 * (logp - self.logp_const))) + "= {} (Χ² = {})".format(logp, -2 * (logp - self.logp_const)) + ) return logp def prepare_data(self, verbose=False): @@ -142,6 +148,7 @@ def prepare_data(self, verbose=False): the shape of the indices array, lmin, lmax. """ import sacc + data = self.data # Read data input_fname = os.path.join(self.data_folder, self.input_file) @@ -152,8 +159,7 @@ def prepare_data(self, verbose=False): s_b = s if self.cov_Bbl_file: if self.cov_Bbl_file != self.input_file: - cov_Bbl_fname = os.path.join(self.data_folder, - self.cov_Bbl_file) + cov_Bbl_fname = os.path.join(self.data_folder, self.cov_Bbl_file) s_b = sacc.Sacc.load_fits(cov_Bbl_fname) cbbl_extra = True @@ -163,19 +169,19 @@ def prepare_data(self, verbose=False): raise KeyError("You must provide a list of default cuts") # Translation betwen TEB and sacc C_ell types - pol_dict = {"T": "0", - "E": "e", - "B": "b"} - ppol_dict = {"TT": "tt", - "EE": "ee", - "TE": "te", - "ET": "te", - "BB": "bb", - "EB": "eb", - "BE": "eb", - "TB": "tb", - "BT": "tb", - "BB": "bb"} + pol_dict = {"T": "0", "E": "e", "B": "b"} + ppol_dict = { + "TT": "tt", + "EE": "ee", + "TE": "te", + "ET": "te", + "BB": "bb", + "EB": "eb", + "BE": "eb", + "TB": "tb", + "BT": "tb", + "BB": "bb", + } def get_cl_meta(spec): """ @@ -190,11 +196,9 @@ def get_cl_meta(spec): # Experiments/frequencies exp_1, exp_2 = spec["experiments"] # Read off polarization channel combinations - pols = spec.get("polarizations", - default_cuts["polarizations"]).copy() + pols = spec.get("polarizations", default_cuts["polarizations"]).copy() # Read off scale cuts - scls = spec.get("scales", - default_cuts["scales"]).copy() + scls = spec.get("scales", default_cuts["scales"]).copy() # For the same two channels, do not include ET and TE, only TE if exp_1 == exp_2: @@ -207,8 +211,7 @@ def get_cl_meta(spec): else: # Symmetrization if ("TE" in pols) and ("ET" in pols): - symm = spec.get("symmetrize", - default_cuts["symmetrize"]) + symm = spec.get("symmetrize", default_cuts["symmetrize"]) else: symm = False @@ -259,16 +262,19 @@ def get_sacc_names(pol, exp_1, exp_2): for pol in pols: tname_1, tname_2, dtype = get_sacc_names(pol, exp_1, exp_2) lmin, lmax = scls[pol] - ind = s.indices(dtype, # Power spectrum type - (tname_1, tname_2), # Channel combinations - ell__gt=lmin, ell__lt=lmax) # Scale cuts + ind = s.indices( + dtype, # Power spectrum type + (tname_1, tname_2), # Channel combinations + ell__gt=lmin, + ell__lt=lmax, + ) # Scale cuts indices += list(ind) # Note that data in the cov_Bbl file may be in different order. if cbbl_extra: - ind_b = s_b.indices(dtype, - (tname_1, tname_2), - ell__gt=lmin, ell__lt=lmax) + ind_b = s_b.indices( + dtype, (tname_1, tname_2), ell__gt=lmin, ell__lt=lmax + ) indices_b += list(ind_b) if symm and pol == "ET": @@ -308,8 +314,7 @@ def get_sacc_names(pol, exp_1, exp_2): # loop over data["spectra"]. ls, cls, ind = s.get_ell_cl(dtype, tname_1, tname_2, return_ind=True) if cbbl_extra: - ind_b = s_b.indices(dtype, - (tname_1, tname_2)) + ind_b = s_b.indices(dtype, (tname_1, tname_2)) ws = s_b.get_bandpower_windows(ind_b) else: ws = s.get_bandpower_windows(ind) @@ -323,10 +328,8 @@ def get_sacc_names(pol, exp_1, exp_2): if (pol in ["TE", "ET"]) and symm: pol2 = pol[::-1] pols.remove(pol2) - tname_1, tname_2, dtype = get_sacc_names(pol2, - exp_1, exp_2) - ind2 = s.indices(dtype, - (tname_1, tname_2)) + tname_1, tname_2, dtype = get_sacc_names(pol2, exp_1, exp_2) + ind2 = s.indices(dtype, (tname_1, tname_2)) cls2 = s.get_ell_cl(dtype, tname_1, tname_2)[1] cls = 0.5 * (cls + cls2) @@ -334,8 +337,7 @@ def get_sacc_names(pol, exp_1, exp_2): mat_compress[index_sofar + i, j1] = 0.5 mat_compress[index_sofar + i, j2] = 0.5 if cbbl_extra: - ind2_b = s_b.indices(dtype, - (tname_1, tname_2)) + ind2_b = s_b.indices(dtype, (tname_1, tname_2)) for i, (j1, j2) in enumerate(zip(ind_b, ind2_b)): mat_compress_b[index_sofar + i, j1] = 0.5 mat_compress_b[index_sofar + i, j2] = 0.5 @@ -347,23 +349,24 @@ def get_sacc_names(pol, exp_1, exp_2): mat_compress_b[index_sofar + i, j1] = 1 # The fields marked with # below aren't really used, but # we store them just in case. - self.spec_meta.append({"ids": (index_sofar + - np.arange(cls.size, - dtype=int)), - "pol": ppol_dict[pol], - "hasYX_xsp": pol in ["ET", "BE", "BT"], # For symm - "t1": exp_1, - "t2": exp_2, - "leff": ls, - "cl_data": cls, - "bpw": ws}) + self.spec_meta.append( + { + "ids": (index_sofar + np.arange(cls.size, dtype=int)), + "pol": ppol_dict[pol], + "hasYX_xsp": pol in ["ET", "BE", "BT"], # For symm + "t1": exp_1, + "t2": exp_2, + "leff": ls, + "cl_data": cls, + "bpw": ws, + } + ) index_sofar += cls.size if not cbbl_extra: mat_compress_b = mat_compress # Put data and covariance in the right order. self.data_vec = np.dot(mat_compress, s.mean) - self.cov = np.dot(mat_compress_b, - s_b.covariance.covmat.dot(mat_compress_b.T)) + self.cov = np.dot(mat_compress_b, s_b.covariance.covmat.dot(mat_compress_b.T)) self.inv_cov = np.linalg.inv(self.cov) self.logp_const = np.log(2 * np.pi) * (-len(self.data_vec) / 2) self.logp_const -= 0.5 * np.linalg.slogdet(self.cov)[1] @@ -399,30 +402,36 @@ def _get_power_spectra(self, cmbfg): ps_vec = np.zeros_like(self.data_vec) DlsObs = dict() # Rescale l_bpws because cmbfg spectra start from first element of l_bpws (l=2) - ell = self.l_bpws - self.l_bpws[0] - - for m in self.spec_meta: - p = m["pol"] - i = m["ids"] - w = m["bpw"].weight.T - - if p in ['tt', 'ee', 'bb']: - DlsObs[p, m['t1'], m['t2']] = cmbfg[p, m['t1'], m['t2']][ell] + if self.l_bpws[-1] <= self.lmax_theory: + ell = self.l_bpws - self.l_bpws[0] + else: + ell = np.arange(self.l_bpws[0], self.lmax_theory) + + for m0 in self.spec_meta: + p0 = m0["pol"] + i0 = m0["ids"] + w0 = m0["bpw"].weight.T + + if p0 in ["tt", "ee", "bb"]: + DlsObs[p0, m0["t1"], m0["t2"]] = cmbfg[p0, m0["t1"], m0["t2"]][ell] else: # ['te','tb','eb'] - if m['hasYX_xsp']: # not symmetrizing - DlsObs[p, m['t2'], m['t1']] = cmbfg[p, m['t2'], m['t1']][ell] + if m0["hasYX_xsp"]: # not symmetrizing + DlsObs[p0, m0["t2"], m0["t1"]] = cmbfg[p0, m0["t2"], m0["t1"]][ell] else: - DlsObs[p, m['t1'], m['t2']] = cmbfg[p, m['t1'], m['t2']][ell] + DlsObs[p0, m0["t1"], m0["t2"]] = cmbfg[p0, m0["t1"], m0["t2"]][ell] # - if self.defaults['symmetrize']: # we average TE and ET (as for data) - DlsObs[p, m['t1'], m['t2']] += cmbfg[p, m['t2'], m['t1']][ell] - DlsObs[p, m['t1'], m['t2']] *= 0.5 - - dls_obs = DlsObs[p, m["t2"], m["t1"]] if m["hasYX_xsp"] \ - else DlsObs[p, m["t1"], m["t2"]] - - clt = w @ dls_obs - ps_vec[i] = clt + if self.defaults["symmetrize"]: # we average TE and ET (as for data) + DlsObs[p0, m0["t1"], m0["t2"]] += cmbfg[p0, m0["t2"], m0["t1"]][ell] + DlsObs[p0, m0["t1"], m0["t2"]] *= 0.5 + + dls_obs = ( + DlsObs[p0, m0["t2"], m0["t1"]] + if m0["hasYX_xsp"] + else DlsObs[p0, m0["t1"], m0["t2"]] + ) + + clt = w0[:, ell] @ dls_obs + ps_vec[i0] = clt return ps_vec diff --git a/soliket/mflike/theoryforge_MFLike.py b/soliket/mflike/theoryforge_MFLike.py index 04c8bd72..960824ad 100644 --- a/soliket/mflike/theoryforge_MFLike.py +++ b/soliket/mflike/theoryforge_MFLike.py @@ -70,7 +70,6 @@ class TheoryForge_MFLike(Theory): systematics_template: dict def initialize(self): - self.lmin = self.spectra["lmin"] self.lmax = self.spectra["lmax"] self.ell = np.arange(self.lmin, self.lmax + 1) @@ -93,23 +92,33 @@ def initialize(self): if hasattr(self.eff_freqs, "__len__"): if not len(self.exp_ch) == len(self.eff_freqs): raise LoggedError( - self.log, "list of effective frequency has to have" - "same length as list of channels!" + self.log, + "list of effective frequency has to have" + "same length as list of channels!", ) # self.bands to be filled with passbands read from sacc file # if mflike is used - self.bands = {f"{expc}_s0": {'nu': [self.eff_freqs[iexpc]], 'bandpass': [1.]} - for iexpc, expc in enumerate(self.exp_ch)} - - self.expected_params_nuis = ["cal_LAT_93", "cal_LAT_145", "cal_LAT_225", - "calT_LAT_93", "calE_LAT_93", - "calT_LAT_145", "calE_LAT_145", - "calT_LAT_225", "calE_LAT_225", - "calG_all", - "alpha_LAT_93", "alpha_LAT_145", - "alpha_LAT_225", - ] + self.bands = { + f"{expc}_s0": {"nu": [self.eff_freqs[iexpc]], "bandpass": [1.0]} + for iexpc, expc in enumerate(self.exp_ch) + } + + self.expected_params_nuis = [ + "cal_LAT_93", + "cal_LAT_145", + "cal_LAT_225", + "calT_LAT_93", + "calE_LAT_93", + "calT_LAT_145", + "calE_LAT_145", + "calT_LAT_225", + "calE_LAT_225", + "calG_all", + "alpha_LAT_93", + "alpha_LAT_145", + "alpha_LAT_225", + ] self.use_systematics_template = bool(self.systematics_template) @@ -121,12 +130,15 @@ def initialize(self): def initialize_with_params(self): # Check that the parameters are the right ones differences = are_different_params_lists( - self.input_params, self.expected_params_nuis, - name_A="given", name_B="expected") + self.input_params, + self.expected_params_nuis, + name_A="given", + name_B="expected", + ) if differences: raise LoggedError( - self.log, "Configuration error in parameters: %r.", - differences) + self.log, "Configuration error in parameters: %r.", differences + ) def must_provide(self, **requirements): # cmbfg_dict is required by mflike @@ -136,8 +148,8 @@ def must_provide(self, **requirements): if "cmbfg_dict" in requirements: req = requirements["cmbfg_dict"] self.ell = req.get("ell", self.ell) - #self.log.info('%d', self.ell[0]) - #self.log.info('%d', self.ell[-1]) + # self.log.info('%d', self.ell[0]) + # self.log.info('%d', self.ell[-1]) self.requested_cls = req.get("requested_cls", self.requested_cls) self.lcuts = req.get("lcuts", self.lcuts) self.exp_ch = req.get("exp_ch", self.exp_ch) @@ -150,9 +162,12 @@ def must_provide(self, **requirements): reqs = dict() # Be sure that CMB is computed at lmax > lmax_data (lcuts from mflike here) reqs["Cl"] = {k: max(c, self.lmax_boltzmann + 1) for k, c in self.lcuts.items()} - reqs["fg_dict"] = {"requested_cls": self.requested_cls, - "ell": self.ell, - "exp_ch": self.exp_ch, "bands": self.bands} + reqs["fg_dict"] = { + "requested_cls": self.requested_cls, + "ell": self.ell, + "exp_ch": self.exp_ch, + "bands": self.bands, + } return reqs def get_cmb_theory(self, **params): @@ -163,13 +178,15 @@ def get_foreground_theory(self, **params): def calculate(self, state, want_derived=False, **params_values_dict): Dls = self.get_cmb_theory(**params_values_dict) - #self.Dls = {s: Dls[s][self.ell] for s, _ in self.lcuts.items()} + # self.Dls = {s: Dls[s][self.ell] for s, _ in self.lcuts.items()} Dls_cut = {s: Dls[s][self.ell] for s, _ in self.lcuts.items()} - params_values_nocosmo = {k: params_values_dict[k] for k in ( - self.expected_params_nuis)} + params_values_nocosmo = { + k: params_values_dict[k] for k in (self.expected_params_nuis) + } fg_dict = self.get_foreground_theory(**params_values_nocosmo) - state["cmbfg_dict"] = self.get_modified_theory(Dls_cut, - fg_dict, **params_values_nocosmo) + state["cmbfg_dict"] = self.get_modified_theory( + Dls_cut, fg_dict, **params_values_nocosmo + ) def get_cmbfg_dict(self): return self.current_state["cmbfg_dict"] @@ -197,8 +214,7 @@ def get_modified_theory(self, Dls, fg_dict, **params): for f1 in self.exp_ch: for f2 in self.exp_ch: for s in self.requested_cls: - cmbfg_dict[s, f1, f2] = (Dls[s] + - fg_dict[s, 'all', f1, f2]) + cmbfg_dict[s, f1, f2] = Dls[s] + fg_dict[s, "all", f1, f2] # Apply alm based calibration factors cmbfg_dict = self._get_calibrated_spectra(cmbfg_dict, **nuis_params) @@ -245,15 +261,21 @@ def _get_calibrated_spectra(self, dls_dict, **nuis_params): cal_pars = {} if "tt" in self.requested_cls or "te" in self.requested_cls: - cal = (nuis_params["calG_all"] * - np.array([nuis_params[f"cal_{exp}"] * nuis_params[f"calT_{exp}"] - for exp in self.exp_ch])) + cal = nuis_params["calG_all"] * np.array( + [ + nuis_params[f"cal_{exp}"] * nuis_params[f"calT_{exp}"] + for exp in self.exp_ch + ] + ) cal_pars["t"] = 1 / cal if "ee" in self.requested_cls or "te" in self.requested_cls: - cal = (nuis_params["calG_all"] * - np.array([nuis_params[f"cal_{exp}"] * nuis_params[f"calE_{exp}"] - for exp in self.exp_ch])) + cal = nuis_params["calG_all"] * np.array( + [ + nuis_params[f"cal_{exp}"] * nuis_params[f"calE_{exp}"] + for exp in self.exp_ch + ] + ) cal_pars["e"] = 1 / cal calib = syl.Calibration_alm(ell=self.ell, spectra=dls_dict) @@ -315,8 +337,9 @@ def _init_template_from_file(self): # decide where to store systematics template. # Currently stored inside syslibrary package - templ_from_file = \ - syl.ReadTemplateFromFile(rootname=self.systematics_template["rootname"]) + templ_from_file = syl.ReadTemplateFromFile( + rootname=self.systematics_template["rootname"] + ) self.dltempl_from_file = templ_from_file(ell=self.ell) def _get_template_from_file(self, dls_dict, **nuis_params): @@ -333,13 +356,16 @@ def _get_template_from_file(self, dls_dict, **nuis_params): # templ_pars=[nuis_params['templ_'+str(fr)] for fr in self.exp_ch] # templ_pars currently hard-coded # but ideally should be passed as input nuisance - templ_pars = {cls: np.zeros((len(self.exp_ch), len(self.exp_ch))) - for cls in self.requested_cls} + templ_pars = { + cls: np.zeros((len(self.exp_ch), len(self.exp_ch))) + for cls in self.requested_cls + } for cls in self.requested_cls: for i1, f1 in enumerate(self.exp_ch): for i2, f2 in enumerate(self.exp_ch): - dls_dict[cls, f1, f2] += (templ_pars[cls][i1][i2] * - self.dltempl_from_file[cls, f1, f2]) + dls_dict[cls, f1, f2] += ( + templ_pars[cls][i1][i2] * self.dltempl_from_file[cls, f1, f2] + ) return dls_dict diff --git a/soliket/poisson/poisson.py b/soliket/poisson/poisson.py index db2b8fdd..80d8da56 100644 --- a/soliket/poisson/poisson.py +++ b/soliket/poisson/poisson.py @@ -23,13 +23,11 @@ def _get_catalog(self): return catalog def _get_rate_fn(self, **kwargs): - """Returns a callable rate function that takes each of 'columns' as kwargs. - """ + """Returns a callable rate function that takes each of 'columns' as kwargs.""" raise NotImplementedError def _get_n_expected(self, **kwargs): - """Computes and returns the integral of the rate function - """ + """Computes and returns the integral of the rate function""" raise NotImplementedError def logp(self, **params_values): diff --git a/soliket/poisson/poisson_data.py b/soliket/poisson/poisson_data.py index c10bc572..c9054af3 100644 --- a/soliket/poisson/poisson_data.py +++ b/soliket/poisson/poisson_data.py @@ -26,12 +26,16 @@ def __init__(self, name, catalog, columns, samples=None): if samples is not None: for c in columns: if c not in samples: - raise ValueError(f"If providing samples, must have samples \ - for all columns: {columns}") + raise ValueError( + f"If providing samples, must have samples \ + for all columns: {columns}" + ) if "prior" not in samples: - raise ValueError('Must provide value of interim prior \ - for all samples, under "prior" key!') + raise ValueError( + 'Must provide value of interim prior \ + for all samples, under "prior" key!' + ) assert all( [samples[k].shape == samples["prior"].shape for k in samples] @@ -59,8 +63,9 @@ def loglike(self, rate_fn, n_expected, broadcastable=False): # Simple case; no uncertainties if self.samples is None: if broadcastable: - rate_densities = rate_fn(**{c: self.catalog[c].values for - c in self.columns}) + rate_densities = rate_fn( + **{c: self.catalog[c].values for c in self.columns} + ) else: rate_densities = np.array( [ @@ -73,8 +78,10 @@ def loglike(self, rate_fn, n_expected, broadcastable=False): else: # Eqn (11) of DFM, Hogg & Morton (https://arxiv.org/pdf/1406.3020.pdf) - summand = rate_fn(**{c: self.samples[c] for - c in self.columns}) / self.samples["prior"] + summand = ( + rate_fn(**{c: self.samples[c] for c in self.columns}) + / self.samples["prior"] + ) l_k = 1 / self.N_k * summand.sum(axis=1) assert l_k.shape == (self._len,) return -n_expected + sum(np.log(l_k)) diff --git a/soliket/ps/ps.py b/soliket/ps/ps.py index cb97538f..3139a96c 100644 --- a/soliket/ps/ps.py +++ b/soliket/ps/ps.py @@ -17,7 +17,7 @@ def _get_Cl(self): def _get_theory(self, **params_values): cl_theory = self._get_Cl() - return cl_theory[self.kind][:self.lmax] + return cl_theory[self.kind][: self.lmax] class BinnedPSLikelihood(PSLikelihood): @@ -25,8 +25,9 @@ class BinnedPSLikelihood(PSLikelihood): def initialize(self): self.binning_matrix = self._get_binning_matrix() - self.bin_centers = \ - self.binning_matrix.dot(np.arange(self.binning_matrix.shape[1])) + self.bin_centers = self.binning_matrix.dot( + np.arange(self.binning_matrix.shape[1]) + ) super().initialize() @classmethod @@ -41,4 +42,4 @@ def _get_data(self): def _get_theory(self, **params_values): cl_theory = self._get_Cl() - return self.binning_matrix.dot(cl_theory[self.kind][:self.lmax]) + return self.binning_matrix.dot(cl_theory[self.kind][: self.lmax]) diff --git a/soliket/tests/__init__.py b/soliket/tests/__init__.py index 62c1a1e1..017a4254 100644 --- a/soliket/tests/__init__.py +++ b/soliket/tests/__init__.py @@ -1,3 +1,3 @@ import soliket -__all__ = ['soliket'] +__all__ = ["soliket"] diff --git a/soliket/tests/conftest.py b/soliket/tests/conftest.py index f255cc44..29e0ec9c 100644 --- a/soliket/tests/conftest.py +++ b/soliket/tests/conftest.py @@ -3,7 +3,7 @@ def pytest_collection_modifyitems(config, items): - if sys.platform.startswith('win'): + if sys.platform.startswith("win"): skip_on_windows = pytest.mark.skip(reason="Skipped on Windows") for item in items: if "require_ccl" in item.keywords: diff --git a/soliket/tests/test_bandpass.py b/soliket/tests/test_bandpass.py index 4616a3b5..a1452359 100644 --- a/soliket/tests/test_bandpass.py +++ b/soliket/tests/test_bandpass.py @@ -6,12 +6,14 @@ bandpass_params = { "bandint_shift_LAT_93": 0.0, "bandint_shift_LAT_145": 0.0, - "bandint_shift_LAT_225": 0.0 + "bandint_shift_LAT_225": 0.0, } -bands = {"LAT_93_s0": {"nu": [93], "bandpass": [1.]}, - "LAT_145_s0": {"nu": [145], "bandpass": [1.]}, - "LAT_225_s0": {"nu": [225], "bandpass": [1.]}} +bands = { + "LAT_93_s0": {"nu": [93], "bandpass": [1.0]}, + "LAT_145_s0": {"nu": [145], "bandpass": [1.0]}, + "LAT_225_s0": {"nu": [225], "bandpass": [1.0]}, +} exp_ch = [k.replace("_s0", "") for k in bands.keys()] @@ -30,9 +32,10 @@ def test_bandpass_model(evaluate_one_info): from soliket.bandpass import BandPass evaluate_one_info["params"] = bandpass_params - evaluate_one_info["theory"] = {"bandpass": { - "external": BandPass, - }, + evaluate_one_info["theory"] = { + "bandpass": { + "external": BandPass, + }, } model = get_model(evaluate_one_info) # noqa F841 @@ -47,18 +50,17 @@ def test_bandpass_read_from_sacc(evaluate_one_info): } model = get_model(evaluate_one_info) - model.add_requirements({"bandint_freqs": {"bands": bands} - }) + model.add_requirements({"bandint_freqs": {"bands": bands}}) - model.logposterior(evaluate_one_info['params']) # force computation of model + model.logposterior(evaluate_one_info["params"]) # force computation of model - lhood = model.likelihood['one'] + lhood = model.likelihood["one"] bandpass = lhood.provider.get_bandint_freqs() bandint_freqs = np.empty_like(exp_ch, dtype=float) for ifr, fr in enumerate(exp_ch): - bandpar = 'bandint_shift_' + fr + bandpar = "bandint_shift_" + fr bandint_freqs[ifr] = ( np.asarray(bands[fr + "_s0"]["nu"]) + evaluate_one_info["params"][bandpar] ) @@ -72,21 +74,19 @@ def test_bandpass_top_hat(evaluate_one_info): # now testing top-hat construction evaluate_one_info["params"] = bandpass_params evaluate_one_info["theory"] = { - "bandpass": {"external": BandPass, - "top_hat_band": { - "nsteps": 3, - "bandwidth": 0.5}, - "external_bandpass": {}, - "read_from_sacc": {}, - }, + "bandpass": { + "external": BandPass, + "top_hat_band": {"nsteps": 3, "bandwidth": 0.5}, + "external_bandpass": {}, + "read_from_sacc": {}, + }, } model = get_model(evaluate_one_info) - model.add_requirements({"bandint_freqs": {"bands": bands} - }) - model.logposterior(evaluate_one_info['params']) # force computation of model + model.add_requirements({"bandint_freqs": {"bands": bands}}) + model.logposterior(evaluate_one_info["params"]) # force computation of model - lhood = model.likelihood['one'] + lhood = model.likelihood["one"] bandpass = lhood.provider.get_bandint_freqs() @@ -94,15 +94,18 @@ def test_bandpass_top_hat(evaluate_one_info): nsteps = evaluate_one_info["theory"]["bandpass"]["top_hat_band"]["nsteps"] bandwidth = evaluate_one_info["theory"]["bandpass"]["top_hat_band"]["bandwidth"] for ifr, fr in enumerate(exp_ch): - bandpar = 'bandint_shift_' + fr + bandpar = "bandint_shift_" + fr bd = bands[f"{fr}_s0"] nu_ghz, bp = np.asarray(bd["nu"]), np.asarray(bd["bandpass"]) fr = nu_ghz @ bp / bp.sum() - bandlow = fr * (1 - bandwidth * .5) - bandhigh = fr * (1 + bandwidth * .5) - nub = np.linspace(bandlow + evaluate_one_info["params"][bandpar], - bandhigh + evaluate_one_info["params"][bandpar], - nsteps, dtype=float) + bandlow = fr * (1 - bandwidth * 0.5) + bandhigh = fr * (1 + bandwidth * 0.5) + nub = np.linspace( + bandlow + evaluate_one_info["params"][bandpar], + bandhigh + evaluate_one_info["params"][bandpar], + nsteps, + dtype=float, + ) tranb = _cmb2bb(nub) tranb_norm = np.trapz(_cmb2bb(nub), nub) bandint_freqs.append([nub, tranb / tranb_norm]) @@ -115,33 +118,34 @@ def test_bandpass_external_file(request, evaluate_one_info): from soliket.bandpass import BandPass - filepath = os.path.join(request.config.rootdir, - "soliket/tests/data/") + filepath = os.path.join(request.config.rootdir, "soliket/tests/data/") # now testing reading from external file evaluate_one_info["params"] = bandpass_params evaluate_one_info["theory"] = { - "bandpass": {"external": BandPass, - "data_folder": f"{filepath}", - "top_hat_band": {}, - "external_bandpass": { - "path": "test_bandpass"}, - "read_from_sacc": {}, - }, + "bandpass": { + "external": BandPass, + "data_folder": f"{filepath}", + "top_hat_band": {}, + "external_bandpass": {"path": "test_bandpass"}, + "read_from_sacc": {}, + }, } model = get_model(evaluate_one_info) - model.add_requirements({"bandint_freqs": {"bands": bands} - }) + model.add_requirements({"bandint_freqs": {"bands": bands}}) - model.logposterior(evaluate_one_info['params']) # force computation of model + model.logposterior(evaluate_one_info["params"]) # force computation of model - lhood = model.likelihood['one'] + lhood = model.likelihood["one"] bandpass = lhood.provider.get_bandint_freqs() - path = os.path.normpath(os.path.join( - evaluate_one_info["theory"]["bandpass"]["data_folder"], - evaluate_one_info["theory"]["bandpass"]["external_bandpass"]["path"])) + path = os.path.normpath( + os.path.join( + evaluate_one_info["theory"]["bandpass"]["data_folder"], + evaluate_one_info["theory"]["bandpass"]["external_bandpass"]["path"], + ) + ) arrays = os.listdir(path) external_bandpass = [] diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 97e6bbab..f891b984 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -1,9 +1,7 @@ import numpy as np from cobaya.model import get_model -bias_params = { - "b_lin": 1.1 -} +bias_params = {"b_lin": 1.1} def test_bias_import(): @@ -15,48 +13,52 @@ def test_linear_bias_import(): def test_linear_bias_model(evaluate_one_info, test_cosmology_params): - from soliket.bias import Linear_bias evaluate_one_info["params"] = test_cosmology_params evaluate_one_info["params"].update(bias_params) evaluate_one_info["theory"] = { - "camb": None, - "linear_bias": {"external": Linear_bias} - } + "camb": None, + "linear_bias": {"external": Linear_bias}, + } model = get_model(evaluate_one_info) # noqa F841 def test_linear_bias_compute_grid(evaluate_one_info, test_cosmology_params): - from soliket.bias import Linear_bias evaluate_one_info["params"] = test_cosmology_params evaluate_one_info["params"].update(bias_params) evaluate_one_info["theory"] = { - "camb": None, - "linear_bias": {"external": Linear_bias} - } + "camb": None, + "linear_bias": {"external": Linear_bias}, + } model = get_model(evaluate_one_info) - model.add_requirements({"Pk_grid": {"z": 0., "k_max": 10., - "nonlinear": False, - "vars_pairs": ('delta_tot', 'delta_tot') - }, - "Pk_gg_grid": None, - "Pk_gm_grid": None - }) - - model.logposterior(evaluate_one_info['params']) # force computation of model - - lhood = model.likelihood['one'] - - k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), - nonlinear=False) + model.add_requirements( + { + "Pk_grid": { + "z": 0.0, + "k_max": 10.0, + "nonlinear": False, + "vars_pairs": ("delta_tot", "delta_tot"), + }, + "Pk_gg_grid": None, + "Pk_gm_grid": None, + } + ) + + model.logposterior(evaluate_one_info["params"]) # force computation of model + + lhood = model.likelihood["one"] + + k, z, Pk_mm_lin = lhood.provider.get_Pk_grid( + var_pair=("delta_tot", "delta_tot"), nonlinear=False + ) Pk_gg = lhood.provider.get_Pk_gg_grid() Pk_gm = lhood.provider.get_Pk_gm_grid() - assert np.allclose(Pk_mm_lin * evaluate_one_info["params"]["b_lin"]**2., Pk_gg) + assert np.allclose(Pk_mm_lin * evaluate_one_info["params"]["b_lin"] ** 2.0, Pk_gg) assert np.allclose(Pk_mm_lin * evaluate_one_info["params"]["b_lin"], Pk_gm) diff --git a/soliket/tests/test_cash.py b/soliket/tests/test_cash.py index 1a202a35..a4000949 100644 --- a/soliket/tests/test_cash.py +++ b/soliket/tests/test_cash.py @@ -5,7 +5,6 @@ class cash_theory_calculator(Theory): - def calculate(self, state, want_derived=False, **params_values_dict): state["cash_theory"] = np.arange(params_values_dict["param_test_cash"]) @@ -31,8 +30,9 @@ def test_cash_read_data(request): from soliket.cash import CashCLikelihood - cash_data_path = os.path.join(request.config.rootdir, - "soliket/tests/data/cash_data.txt") + cash_data_path = os.path.join( + request.config.rootdir, "soliket/tests/data/cash_data.txt" + ) cash_lkl = CashCLikelihood({"datapath": cash_data_path}) cash_data = cash_lkl._get_data() @@ -45,8 +45,9 @@ def test_cash_logp(request): from soliket.cash import CashCLikelihood params = {"cash_test_logp": 20} - cash_data_path = os.path.join(request.config.rootdir, - "soliket/tests/data/cash_data.txt") + cash_data_path = os.path.join( + request.config.rootdir, "soliket/tests/data/cash_data.txt" + ) cash_lkl = CashCLikelihood({"datapath": cash_data_path}) cash_logp = cash_lkl.logp(**params) @@ -54,7 +55,6 @@ def test_cash_logp(request): def test_cash(): - data1d, theory1d, data2d, theory2d = toy_data() cashdata1d = CashCData("toy 1d", data1d) diff --git a/soliket/tests/test_ccl.py b/soliket/tests/test_ccl.py index 152c8431..1ce817d5 100644 --- a/soliket/tests/test_ccl.py +++ b/soliket/tests/test_ccl.py @@ -24,17 +24,8 @@ def get_requirements(self): ccl_like_and_theory = { - "likelihood": { - "checkLike": {"external": CheckLike} - }, - "theory": { - "camb": { - }, - "soliket.CCL": { - "kmax": 10.0, - "nonlinear": True - } - } + "likelihood": {"checkLike": {"external": CheckLike}}, + "theory": {"camb": {}, "soliket.CCL": {"kmax": 10.0, "nonlinear": True}}, } diff --git a/soliket/tests/test_clusters.py b/soliket/tests/test_clusters.py index e8936eb1..65bd0c7d 100644 --- a/soliket/tests/test_clusters.py +++ b/soliket/tests/test_clusters.py @@ -15,7 +15,7 @@ "nonlinear": False, "kmax": 10.0, "dark_energy_model": "ppf", - "bbn_predictor": "PArthENoPE_880.2_standard.dat" + "bbn_predictor": "PArthENoPE_880.2_standard.dat", } }, }, @@ -37,7 +37,7 @@ def test_clusters_loglike(evaluate_one_info, test_cosmology_params): lnl = model_fiducial.loglikes({})[0] - assert np.isclose(lnl, -847.22462272, rtol=1.e-3, atol=1.e-5) + assert np.isclose(lnl, -847.22462272, rtol=1.0e-3, atol=1.0e-5) def test_clusters_n_expected(evaluate_one_info, test_cosmology_params): @@ -50,5 +50,5 @@ def test_clusters_n_expected(evaluate_one_info, test_cosmology_params): like = model_fiducial.likelihood["soliket.ClusterLikelihood"] - assert np.isclose(lnl, -847.22462272, rtol=1.e-3, atol=1.e-5) + assert np.isclose(lnl, -847.22462272, rtol=1.0e-3, atol=1.0e-5) assert like._get_n_expected() > 40 diff --git a/soliket/tests/test_cosmopower.py b/soliket/tests/test_cosmopower.py index fc017527..0c912994 100644 --- a/soliket/tests/test_cosmopower.py +++ b/soliket/tests/test_cosmopower.py @@ -59,29 +59,31 @@ "omch2": "omega_cdm", "ns": "n_s", "logA": "ln10^{10}A_s", - "tau": "tau_reio" - } + "tau": "tau_reio", + }, } }, } -@pytest.mark.skipif(not HAS_COSMOPOWER, reason='test requires cosmopower') +@pytest.mark.skipif(not HAS_COSMOPOWER, reason="test requires cosmopower") def test_cosmopower_import(request): from soliket.cosmopower import CosmoPower # noqa F401 -@pytest.mark.skipif(not HAS_COSMOPOWER, reason='test requires cosmopower') +@pytest.mark.skipif(not HAS_COSMOPOWER, reason="test requires cosmopower") def test_cosmopower_theory(request): - info_dict['theory']['soliket.CosmoPower']['network_path'] = \ - os.path.join(request.config.rootdir, 'soliket/cosmopower/data/CP_paper') - model_fiducial = get_model(info_dict) # noqa F841 + info_dict["theory"]["soliket.CosmoPower"]["network_path"] = os.path.join( + request.config.rootdir, "soliket/cosmopower/data/CP_paper" + ) + model_fiducial = get_model(info_dict) # noqa F841 -@pytest.mark.skipif(not HAS_COSMOPOWER, reason='test requires cosmopower') +@pytest.mark.skipif(not HAS_COSMOPOWER, reason="test requires cosmopower") def test_cosmopower_loglike(request): - info_dict['theory']['soliket.CosmoPower']['network_path'] = \ - os.path.join(request.config.rootdir, 'soliket/cosmopower/data/CP_paper') + info_dict["theory"]["soliket.CosmoPower"]["network_path"] = os.path.join( + request.config.rootdir, "soliket/cosmopower/data/CP_paper" + ) model_cp = get_model(info_dict) logL_cp = float(model_cp.loglikes({})[0]) @@ -89,35 +91,27 @@ def test_cosmopower_loglike(request): assert np.isclose(logL_cp, -295.139) -@pytest.mark.skipif(not HAS_COSMOPOWER, reason='test requires cosmopower') +@pytest.mark.skipif(not HAS_COSMOPOWER, reason="test requires cosmopower") def test_cosmopower_against_camb(request): - - info_dict['theory'] = {'camb': {'stop_at_error': True}} + info_dict["theory"] = {"camb": {"stop_at_error": True}} model_camb = get_model(info_dict) logL_camb = float(model_camb.loglikes({})[0]) - camb_cls = model_camb.theory['camb'].get_Cl() + camb_cls = model_camb.theory["camb"].get_Cl() - info_dict['theory'] = { + info_dict["theory"] = { "soliket.CosmoPower": { "stop_at_error": True, - "extra_args": {'lmax': camb_cls['ell'].max()}, - 'network_path': os.path.join(request.config.rootdir, - 'soliket/cosmopower/data/CP_paper'), + "extra_args": {"lmax": camb_cls["ell"].max()}, + "network_path": os.path.join( + request.config.rootdir, "soliket/cosmopower/data/CP_paper" + ), "network_settings": { - "tt": { - "type": "NN", - "log": True, - "filename": "cmb_TT_NN" - }, - "ee": { - "type": "NN", - "log": True, - "filename": "cmb_EE_NN" - }, + "tt": {"type": "NN", "log": True, "filename": "cmb_TT_NN"}, + "ee": {"type": "NN", "log": True, "filename": "cmb_EE_NN"}, "te": { "type": "PCAplusNN", "log": False, - "filename": "cmb_TE_PCAplusNN" + "filename": "cmb_TE_PCAplusNN", }, }, "renames": { @@ -125,16 +119,16 @@ def test_cosmopower_against_camb(request): "omch2": "omega_cdm", "ns": "n_s", "logA": "ln10^{10}A_s", - "tau": "tau_reio" - } + "tau": "tau_reio", + }, } } model_cp = get_model(info_dict) logL_cp = float(model_cp.loglikes({})[0]) - cp_cls = model_cp.theory['soliket.CosmoPower'].get_Cl() + cp_cls = model_cp.theory["soliket.CosmoPower"].get_Cl() - nanmask = ~np.isnan(cp_cls['tt']) + nanmask = ~np.isnan(cp_cls["tt"]) - assert np.allclose(cp_cls['tt'][nanmask], camb_cls['tt'][nanmask], rtol=1.e-2) - assert np.isclose(logL_camb, logL_cp, rtol=1.e-1) + assert np.allclose(cp_cls["tt"][nanmask], camb_cls["tt"][nanmask], rtol=1.0e-2) + assert np.isclose(logL_camb, logL_cp, rtol=1.0e-1) diff --git a/soliket/tests/test_cross_correlation.py b/soliket/tests/test_cross_correlation.py index 337e9e78..69902784 100644 --- a/soliket/tests/test_cross_correlation.py +++ b/soliket/tests/test_cross_correlation.py @@ -6,8 +6,8 @@ pytestmark = pytest.mark.require_ccl -gammakappa_sacc_file = 'soliket/tests/data/des_s-act_kappa.toy-sim.sacc.fits' -gkappa_sacc_file = 'soliket/tests/data/gc_cmass-actdr4_kappa.sacc.fits' +gammakappa_sacc_file = "soliket/tests/data/des_s-act_kappa.toy-sim.sacc.fits" +gkappa_sacc_file = "soliket/tests/data/gc_cmass-actdr4_kappa.sacc.fits" cross_correlation_params = { "b1": 1.0, @@ -35,9 +35,11 @@ def test_galaxykappa_model(request, evaluate_one_info, test_cosmology_params): evaluate_one_info["theory"] = cross_correlation_theory evaluate_one_info["likelihood"] = { - "GalaxyKappaLikelihood": {"external": GalaxyKappaLikelihood, - "datapath": os.path.join(request.config.rootdir, - gkappa_sacc_file)}} + "GalaxyKappaLikelihood": { + "external": GalaxyKappaLikelihood, + "datapath": os.path.join(request.config.rootdir, gkappa_sacc_file), + } + } model = get_model(evaluate_one_info) # noqa F841 @@ -48,11 +50,12 @@ def test_shearkappa_model(request, evaluate_one_info, test_cosmology_params): evaluate_one_info["params"] = test_cosmology_params evaluate_one_info["theory"] = cross_correlation_theory - evaluate_one_info["likelihood"] = {"ShearKappaLikelihood": - {"external": ShearKappaLikelihood, - "datapath": - os.path.join(request.config.rootdir, - gammakappa_sacc_file)}} + evaluate_one_info["likelihood"] = { + "ShearKappaLikelihood": { + "external": ShearKappaLikelihood, + "datapath": os.path.join(request.config.rootdir, gammakappa_sacc_file), + } + } model = get_model(evaluate_one_info) # noqa F841 @@ -65,10 +68,12 @@ def test_galaxykappa_like(request, evaluate_one_info, test_cosmology_params): evaluate_one_info["theory"] = cross_correlation_theory evaluate_one_info["likelihood"] = { - "GalaxyKappaLikelihood": {"external": GalaxyKappaLikelihood, - "datapath": os.path.join(request.config.rootdir, - gkappa_sacc_file), - "use_spectra": [('gc_cmass', 'ck_actdr4')]}} + "GalaxyKappaLikelihood": { + "external": GalaxyKappaLikelihood, + "datapath": os.path.join(request.config.rootdir, gkappa_sacc_file), + "use_spectra": [("gc_cmass", "ck_actdr4")], + } + } model = get_model(evaluate_one_info) loglikes, derived = model.loglikes() @@ -87,14 +92,16 @@ def test_shearkappa_like(request, evaluate_one_info): test_datapath = os.path.join(rootdir, cs82_file) evaluate_one_info["likelihood"] = { - "ShearKappaLikelihood": {"external": ShearKappaLikelihood, - "datapath": test_datapath} + "ShearKappaLikelihood": { + "external": ShearKappaLikelihood, + "datapath": test_datapath, + } } # Cosmological parameters for the test data, digitized from # Fig. 3 and Eq. 8 of Hall & Taylor (2014). # See https://github.com/simonsobs/SOLikeT/pull/58 for validation plots - evaluate_one_info['params'] = { + evaluate_one_info["params"] = { "omch2": 0.118, # Planck + lensing + WP + highL "ombh2": 0.0222, "H0": 68.0, @@ -102,7 +109,7 @@ def test_shearkappa_like(request, evaluate_one_info): "As": 2.1e-9, "tau": 0.094, "mnu": 0.0, - "nnu": 3.046 + "nnu": 3.046, } model = get_model(evaluate_one_info) @@ -124,36 +131,38 @@ def test_shearkappa_tracerselect(request, evaluate_one_info, test_cosmology_para test_datapath = os.path.join(rootdir, gammakappa_sacc_file) evaluate_one_info["likelihood"] = { - "ShearKappaLikelihood": {"external": ShearKappaLikelihood, - "datapath": test_datapath, - 'use_spectra': 'all'} + "ShearKappaLikelihood": { + "external": ShearKappaLikelihood, + "datapath": test_datapath, + "use_spectra": "all", + } } info_onebin = copy.deepcopy(evaluate_one_info) - info_onebin['likelihood']['ShearKappaLikelihood']['use_spectra'] = \ - [('gs_des_bin1', 'ck_act')] + info_onebin["likelihood"]["ShearKappaLikelihood"]["use_spectra"] = [ + ("gs_des_bin1", "ck_act") + ] info_twobin = copy.deepcopy(evaluate_one_info) - info_twobin['likelihood']['ShearKappaLikelihood']['use_spectra'] = \ - [ - ('gs_des_bin1', 'ck_act'), - ('gs_des_bin3', 'ck_act'), - ] + info_twobin["likelihood"]["ShearKappaLikelihood"]["use_spectra"] = [ + ("gs_des_bin1", "ck_act"), + ("gs_des_bin3", "ck_act"), + ] model = get_model(evaluate_one_info) loglikes, derived = model.loglikes() - lhood = model.likelihood['ShearKappaLikelihood'] + lhood = model.likelihood["ShearKappaLikelihood"] model_onebin = get_model(info_onebin) loglikes_onebin, derived_onebin = model_onebin.loglikes() - lhood_onebin = model_onebin.likelihood['ShearKappaLikelihood'] + lhood_onebin = model_onebin.likelihood["ShearKappaLikelihood"] model_twobin = get_model(info_twobin) loglikes_twobin, derived_twobin = model_twobin.loglikes() - lhood_twobin = model_twobin.likelihood['ShearKappaLikelihood'] + lhood_twobin = model_twobin.likelihood["ShearKappaLikelihood"] n_ell_perbin = len(lhood.data.x) // 4 @@ -161,9 +170,15 @@ def test_shearkappa_tracerselect(request, evaluate_one_info, test_cosmology_para assert np.allclose(lhood.data.y[:n_ell_perbin], lhood_onebin.data.y) assert 2 * n_ell_perbin == len(lhood_twobin.data.x) - assert np.allclose(np.concatenate([lhood.data.y[:n_ell_perbin], - lhood.data.y[2 * n_ell_perbin:3 * n_ell_perbin]]), - lhood_twobin.data.y) + assert np.allclose( + np.concatenate( + [ + lhood.data.y[:n_ell_perbin], + lhood.data.y[2 * n_ell_perbin:3 * n_ell_perbin], + ] + ), + lhood_twobin.data.y, + ) def test_shearkappa_hartlap(request, evaluate_one_info): @@ -177,23 +192,26 @@ def test_shearkappa_hartlap(request, evaluate_one_info): test_datapath = os.path.join(rootdir, cs82_file) evaluate_one_info["likelihood"] = { - "ShearKappaLikelihood": {"external": ShearKappaLikelihood, - "datapath": test_datapath} + "ShearKappaLikelihood": { + "external": ShearKappaLikelihood, + "datapath": test_datapath, + } } # Cosmological parameters for the test data, digitized from # Fig. 3 and Eq. 8 of Hall & Taylor (2014). # See https://github.com/simonsobs/SOLikeT/pull/58 for validation plots - evaluate_one_info['params'] = \ - {"omch2": 0.118, # Planck + lensing + WP + highL - "ombh2": 0.0222, - "H0": 68.0, - "ns": 0.962, - # "As": 2.1e-9, - "As": 2.5e-9, # offset the theory to upweight inv_cov in loglike - "tau": 0.094, - "mnu": 0.0, - "nnu": 3.046} + evaluate_one_info["params"] = { + "omch2": 0.118, # Planck + lensing + WP + highL + "ombh2": 0.0222, + "H0": 68.0, + "ns": 0.962, + # "As": 2.1e-9, + "As": 2.5e-9, # offset the theory to upweight inv_cov in loglike + "tau": 0.094, + "mnu": 0.0, + "nnu": 3.046, + } model = get_model(evaluate_one_info) loglikes, derived = model.loglikes() @@ -203,8 +221,9 @@ def test_shearkappa_hartlap(request, evaluate_one_info): model = get_model(evaluate_one_info) loglikes_hartlap, derived = model.loglikes() - assert np.isclose(np.abs(loglikes - loglikes_hartlap), 0.0010403, - rtol=1.e-5, atol=1.e-5) + assert np.isclose( + np.abs(loglikes - loglikes_hartlap), 0.0010403, rtol=1.0e-5, atol=1.0e-5 + ) def test_shearkappa_deltaz(request, evaluate_one_info, test_cosmology_params): @@ -213,12 +232,13 @@ def test_shearkappa_deltaz(request, evaluate_one_info, test_cosmology_params): evaluate_one_info["params"] = test_cosmology_params evaluate_one_info["theory"] = cross_correlation_theory - evaluate_one_info["likelihood"] = \ - {"ShearKappaLikelihood": - {"external": ShearKappaLikelihood, - "datapath": os.path.join(request.config.rootdir, - gammakappa_sacc_file), - "z_nuisance_mode": "deltaz"}} + evaluate_one_info["likelihood"] = { + "ShearKappaLikelihood": { + "external": ShearKappaLikelihood, + "datapath": os.path.join(request.config.rootdir, gammakappa_sacc_file), + "z_nuisance_mode": "deltaz", + } + } model = get_model(evaluate_one_info) loglikes, derived = model.loglikes() @@ -232,12 +252,13 @@ def test_shearkappa_m(request, evaluate_one_info, test_cosmology_params): evaluate_one_info["params"] = test_cosmology_params evaluate_one_info["theory"] = cross_correlation_theory - evaluate_one_info["likelihood"] = \ - {"ShearKappaLikelihood": - {"external": ShearKappaLikelihood, - "datapath": os.path.join(request.config.rootdir, - gammakappa_sacc_file), - "m_nuisance_mode": True}} + evaluate_one_info["likelihood"] = { + "ShearKappaLikelihood": { + "external": ShearKappaLikelihood, + "datapath": os.path.join(request.config.rootdir, gammakappa_sacc_file), + "m_nuisance_mode": True, + } + } model = get_model(evaluate_one_info) loglikes, derived = model.loglikes() @@ -251,12 +272,13 @@ def test_shearkappa_ia_nla_noevo(request, evaluate_one_info, test_cosmology_para evaluate_one_info["params"] = test_cosmology_params evaluate_one_info["theory"] = cross_correlation_theory - evaluate_one_info["likelihood"] = \ - {"ShearKappaLikelihood": - {"external": ShearKappaLikelihood, - "datapath": os.path.join(request.config.rootdir, - gammakappa_sacc_file), - "ia_mode": 'nla-noevo'}} + evaluate_one_info["likelihood"] = { + "ShearKappaLikelihood": { + "external": ShearKappaLikelihood, + "datapath": os.path.join(request.config.rootdir, gammakappa_sacc_file), + "ia_mode": "nla-noevo", + } + } model = get_model(evaluate_one_info) loglikes, derived = model.loglikes() @@ -270,12 +292,13 @@ def test_shearkappa_ia_nla(request, evaluate_one_info, test_cosmology_params): evaluate_one_info["params"] = test_cosmology_params evaluate_one_info["theory"] = cross_correlation_theory - evaluate_one_info["likelihood"] = \ - {"ShearKappaLikelihood": - {"external": ShearKappaLikelihood, - "datapath": os.path.join(request.config.rootdir, - gammakappa_sacc_file), - "ia_mode": 'nla'}} + evaluate_one_info["likelihood"] = { + "ShearKappaLikelihood": { + "external": ShearKappaLikelihood, + "datapath": os.path.join(request.config.rootdir, gammakappa_sacc_file), + "ia_mode": "nla", + } + } evaluate_one_info["params"]["eta_IA"] = 1.7 @@ -291,12 +314,13 @@ def test_shearkappa_ia_perbin(request, evaluate_one_info, test_cosmology_params) evaluate_one_info["params"] = test_cosmology_params evaluate_one_info["theory"] = cross_correlation_theory - evaluate_one_info["likelihood"] = \ - {"ShearKappaLikelihood": - {"external": ShearKappaLikelihood, - "datapath": os.path.join(request.config.rootdir, - gammakappa_sacc_file), - "ia_mode": 'nla-perbin'}} + evaluate_one_info["likelihood"] = { + "ShearKappaLikelihood": { + "external": ShearKappaLikelihood, + "datapath": os.path.join(request.config.rootdir, gammakappa_sacc_file), + "ia_mode": "nla-perbin", + } + } model = get_model(evaluate_one_info) loglikes, derived = model.loglikes() @@ -310,18 +334,20 @@ def test_shearkappa_hmcode(request, evaluate_one_info, test_cosmology_params): evaluate_one_info["params"] = test_cosmology_params evaluate_one_info["theory"] = cross_correlation_theory - evaluate_one_info["likelihood"] = \ - {"ShearKappaLikelihood": - {"external": ShearKappaLikelihood, - "datapath": os.path.join(request.config.rootdir, - gammakappa_sacc_file)}} + evaluate_one_info["likelihood"] = { + "ShearKappaLikelihood": { + "external": ShearKappaLikelihood, + "datapath": os.path.join(request.config.rootdir, gammakappa_sacc_file), + } + } evaluate_one_info["theory"] = { "camb": { - 'extra_args': { - 'halofit_version': 'mead2020_feedback', 'HMCode_logT_AGN': 7.8 + "extra_args": { + "halofit_version": "mead2020_feedback", + "HMCode_logT_AGN": 7.8, } }, - "ccl": {"external": CCL, "nonlinear": False} + "ccl": {"external": CCL, "nonlinear": False}, } model = get_model(evaluate_one_info) diff --git a/soliket/tests/test_foreground.py b/soliket/tests/test_foreground.py index b3557614..83d96f68 100644 --- a/soliket/tests/test_foreground.py +++ b/soliket/tests/test_foreground.py @@ -95,7 +95,14 @@ def test_foreground_compute(evaluate_one_info): fg_model = lhood.provider.get_fg_dict() fg_model_test = get_fg( - exp_ch, eff_freqs, ell, ell_0, nu_0, requested_cls, components, evaluate_one_info + exp_ch, + eff_freqs, + ell, + ell_0, + nu_0, + requested_cls, + components, + evaluate_one_info, ) for k in fg_model_test.keys(): diff --git a/soliket/tests/test_gaussian.py b/soliket/tests/test_gaussian.py index 1b0517ef..c0864ba5 100644 --- a/soliket/tests/test_gaussian.py +++ b/soliket/tests/test_gaussian.py @@ -24,7 +24,7 @@ def toy_data(): # Generate arbitrary covariance matrix, partition into parts full_cov = make_spd_matrix(n1 + n2 + n3, random_state=1234) cov1 = full_cov[:n1, :n1] - cov2 = full_cov[n1: n1 + n2, n1: n1 + n2] + cov2 = full_cov[n1:n1 + n2, n1:n1 + n2] cov3 = full_cov[n1 + n2:, n1 + n2:] data1 = GaussianData(name1, x1, y1, cov1) @@ -33,9 +33,9 @@ def toy_data(): cross_cov = CrossCov( { - (name1, name2): full_cov[:n1, n1: n1 + n2], + (name1, name2): full_cov[:n1, n1:n1 + n2], (name1, name3): full_cov[:n1, n1 + n2:], - (name2, name3): full_cov[n1: n1 + n2, n1 + n2:], + (name2, name3): full_cov[n1:n1 + n2, n1 + n2:], } ) @@ -50,9 +50,15 @@ def test_gaussian(): name1, name2, name3 = [d.name for d in datalist] data1, data2, data3 = datalist - assert (multi.cross_covs[(name1, name2)] == multi.cross_covs[(name2, name1)].T).all() - assert (multi.cross_covs[(name1, name3)] == multi.cross_covs[(name3, name1)].T).all() - assert (multi.cross_covs[(name2, name3)] == multi.cross_covs[(name3, name2)].T).all() + assert ( + multi.cross_covs[(name1, name2)] == multi.cross_covs[(name2, name1)].T + ).all() + assert ( + multi.cross_covs[(name1, name3)] == multi.cross_covs[(name3, name1)].T + ).all() + assert ( + multi.cross_covs[(name2, name3)] == multi.cross_covs[(name3, name2)].T + ).all() assert (multi.cross_covs[(name1, name1)] == data1.cov).all() assert (multi.cross_covs[(name2, name2)] == data2.cov).all() @@ -60,31 +66,32 @@ def test_gaussian(): def test_gaussian_hartlap(): - np.random.seed(1234) name1 = "A" n1 = 10 x1 = np.arange(n1) - y1th = x1**2. + y1th = x1**2.0 y1 = np.random.random(n1) nsims1 = 50 cov1 = make_spd_matrix(n1, random_state=1234) data1 = GaussianData(name1, x1, y1, cov1) - data1_simcov = GaussianData(name1 + 'simcov', x1, y1, cov1, - ncovsims=nsims1) - data1_manysimcov = GaussianData(name1 + 'simcov', x1, y1, cov1, - ncovsims=(100 * nsims1)) + data1_simcov = GaussianData(name1 + "simcov", x1, y1, cov1, ncovsims=nsims1) + data1_manysimcov = GaussianData( + name1 + "simcov", x1, y1, cov1, ncovsims=(100 * nsims1) + ) hartlap_factor = (nsims1 - n1 - 2) / (nsims1 - 1) hartlap_manyfactor = (100 * nsims1 - n1 - 2) / (100 * nsims1 - 1) - assert np.isclose(data1.loglike(y1th), - data1_simcov.loglike(y1th) / hartlap_factor, - rtol=1.e-3) + assert np.isclose( + data1.loglike(y1th), data1_simcov.loglike(y1th) / hartlap_factor, rtol=1.0e-3 + ) - assert np.isclose(data1.loglike(y1th), - data1_manysimcov.loglike(y1th) / hartlap_manyfactor, - rtol=1.e-5) + assert np.isclose( + data1.loglike(y1th), + data1_manysimcov.loglike(y1th) / hartlap_manyfactor, + rtol=1.0e-5, + ) diff --git a/soliket/tests/test_lensing.py b/soliket/tests/test_lensing.py index 3fb85d0f..d8af9d8c 100644 --- a/soliket/tests/test_lensing.py +++ b/soliket/tests/test_lensing.py @@ -7,17 +7,18 @@ # Cosmological parameters for the test data from SO sims # See https://github.com/simonsobs/SOLikeT/pull/101 for validation plots fiducial_params = { - 'omch2': 0.1203058, - 'ombh2': 0.02219218, - 'H0': 67.02393, - 'ns': 0.9625356, - 'As': 2.15086031154146e-9, - 'mnu': 0.06, - 'tau': 0.06574325, - 'nnu': 3.04} + "omch2": 0.1203058, + "ombh2": 0.02219218, + "H0": 67.02393, + "ns": 0.9625356, + "As": 2.15086031154146e-9, + "mnu": 0.06, + "tau": 0.06574325, + "nnu": 3.04, +} info = {"theory": {"camb": {"extra_args": {"kmax": 0.9}}}} -info['params'] = fiducial_params +info["params"] = fiducial_params def test_lensing_import(request): @@ -26,6 +27,7 @@ def test_lensing_import(request): def test_lensing_like(request): from cobaya.install import install + install( {"likelihood": {"soliket.lensing.LensingLikelihood": None}}, path=packages_path, @@ -51,6 +53,7 @@ def test_lensing_ccl_limber(request): """ from cobaya.install import install + install( {"likelihood": {"soliket.lensing.LensingLikelihood": None}}, path=packages_path, @@ -66,19 +69,20 @@ def test_lensing_ccl_limber(request): info_dict = deepcopy(info) # Neutrino mass put to 0 as far as it is not included in the ccl wrapper - info_dict['params']["mnu"] = 0 - info_dict['params']["omnuh2"] = 0 - info_dict['likelihood'] = {"LensingLikelihood": {"external": LensingLikelihood}} + info_dict["params"]["mnu"] = 0 + info_dict["params"]["omnuh2"] = 0 + info_dict["likelihood"] = {"LensingLikelihood": {"external": LensingLikelihood}} model = get_model(info_dict) model.loglikes({}) - cl_camb = model.likelihood['LensingLikelihood']._get_theory() + cl_camb = model.likelihood["LensingLikelihood"]._get_theory() - info_dict["likelihood"] = {"LensingLikelihood": {"external": LensingLikelihood, - "pp_ccl": True}} + info_dict["likelihood"] = { + "LensingLikelihood": {"external": LensingLikelihood, "pp_ccl": True} + } info_dict["theory"]["soliket.CCL"] = {"kmax": 10, "nonlinear": True} model = get_model(info_dict) model.loglikes({}) - cl_ccl = model.likelihood['LensingLikelihood']._get_theory() + cl_ccl = model.likelihood["LensingLikelihood"]._get_theory() assert np.any(np.not_equal(cl_ccl, cl_camb)) assert np.allclose(cl_ccl, cl_camb, rtol=1e-2, atol=0) diff --git a/soliket/tests/test_mflike.py b/soliket/tests/test_mflike.py index 25d58212..2c2e83d4 100644 --- a/soliket/tests/test_mflike.py +++ b/soliket/tests/test_mflike.py @@ -48,22 +48,15 @@ } -if Version(camb.__version__) >= Version('1.4'): - chi2s = {"tt": 544.9017, - "te": 136.6051, - "ee": 166.1897, - "tt-te-et-ee": 787.9529} +if Version(camb.__version__) >= Version("1.4"): + chi2s = {"tt": 544.9017, "te": 136.6051, "ee": 166.1897, "tt-te-et-ee": 787.9529} else: - chi2s = {"tt": 544.8797, - "te-et": 151.8197, - "ee": 166.2835, - "tt-te-et-ee": 787.9843} + chi2s = {"tt": 544.8797, "te-et": 151.8197, "ee": 166.2835, "tt-te-et-ee": 787.9843} pre = "test_data_sacc_" class Test_mflike: - @classmethod def setup_class(cls): from cobaya.install import install @@ -79,7 +72,6 @@ def setup_class(cls): @pytest.mark.usefixtures("test_cosmology_params") def test_mflike(self, test_cosmology_params): - # As of now, there is not a mechanism # in soliket to ensure there is .loglike that can be called like this # w/out cobaya @@ -89,9 +81,9 @@ def test_mflike(self, test_cosmology_params): pars = camb.set_params(**test_cosmology_params) results = camb.get_results(pars) powers = results.get_cmb_power_spectra(pars, CMB_unit="muK") - cl_dict = {k: powers["total"][:, v] for - k, v in {"tt": 0, "ee": 1, "te": 3}.items()} - + cl_dict = { + k: powers["total"][:, v] for k, v in {"tt": 0, "ee": 1, "te": 3}.items() + } BP = soliket.BandPass() FG = soliket.Foreground() @@ -102,13 +94,11 @@ def test_mflike(self, test_cosmology_params): requested_cls = TF.requested_cls BP.bands = bands - BP.exp_ch = [k.replace("_s0", "") for k in bands.keys() - if "_s0" in k] + BP.exp_ch = [k.replace("_s0", "") for k in bands.keys() if "_s0" in k] bandpass = BP._bandpass_construction(**nuisance_params) for select, chi2 in chi2s.items(): - my_mflike = TestMFLike( { "external": TestMFLike, @@ -130,11 +120,13 @@ def test_mflike(self, test_cosmology_params): ell_cut = my_mflike.l_bpws dls_cut = {s: cl_dict[s][ell_cut] for s, _ in my_mflike.lcuts.items()} - fg_dict = FG._get_foreground_model(requested_cls=requested_cls, - ell=ell_cut, - exp_ch=exp_ch, - bandint_freqs=bandpass, - **nuisance_params) + fg_dict = FG._get_foreground_model( + requested_cls=requested_cls, + ell=ell_cut, + exp_ch=exp_ch, + bandint_freqs=bandpass, + **nuisance_params, + ) dlobs_dict = TF.get_modified_theory(dls_cut, fg_dict, **nuisance_params) loglike = my_mflike.loglike(dlobs_dict) @@ -145,7 +137,6 @@ def test_mflike(self, test_cosmology_params): @pytest.mark.usefixtures("test_cosmology_params") def test_cobaya(self, test_cosmology_params): - info = { "likelihood": { "soliket.mflike.TestMFLike": { @@ -164,16 +155,20 @@ def test_cobaya(self, test_cosmology_params): }, }, }, - "theory": {"camb": {"extra_args": {"lens_potential_accuracy": 1}, - "stop_at_error": True}}, + "theory": { + "camb": { + "extra_args": {"lens_potential_accuracy": 1}, + "stop_at_error": True, + } + }, "params": test_cosmology_params, "modules": packages_path, "debug": True, } - info["theory"]["soliket.TheoryForge_MFLike"] = {'stop_at_error': True} - info["theory"]["soliket.Foreground"] = {'stop_at_error': True} - info["theory"]["soliket.BandPass"] = {'stop_at_error': True} + info["theory"]["soliket.TheoryForge_MFLike"] = {"stop_at_error": True} + info["theory"]["soliket.Foreground"] = {"stop_at_error": True} + info["theory"]["soliket.BandPass"] = {"stop_at_error": True} from cobaya.model import get_model model = get_model(info) diff --git a/soliket/tests/test_multi.py b/soliket/tests/test_multi.py index a00d66cb..181378a1 100644 --- a/soliket/tests/test_multi.py +++ b/soliket/tests/test_multi.py @@ -18,8 +18,10 @@ def test_multi(test_cosmology_params): camb_options = {"extra_args": {"lens_potential_accuracy": 1}} - fg_params = {"a_tSZ": {"prior": {"min": 3.0, "max": 3.6}}, - "a_kSZ": {"prior": {"min": 1.4, "max": 1.8}}} + fg_params = { + "a_tSZ": {"prior": {"min": 3.0, "max": 3.6}}, + "a_kSZ": {"prior": {"min": 1.4, "max": 1.8}}, + } mflike_params = {**test_cosmology_params, **nuisance_params} mflike_params.update(fg_params) @@ -28,24 +30,31 @@ def test_multi(test_cosmology_params): info = { "likelihood": { "soliket.gaussian.MultiGaussianLikelihood": { - "components": ["soliket.mflike.TestMFLike", "soliket.LensingLikelihood"], + "components": [ + "soliket.mflike.TestMFLike", + "soliket.LensingLikelihood", + ], "options": [mflike_options, lensing_options], "stop_at_error": True, } }, - "theory": {"camb": camb_options, - "soliket.TheoryForge_MFLike": {'stop_at_error': True}, - "soliket.Foreground": {"stop_at_error": True}, - "soliket.BandPass": {"stop_at_error": True}}, + "theory": { + "camb": camb_options, + "soliket.TheoryForge_MFLike": {"stop_at_error": True}, + "soliket.Foreground": {"stop_at_error": True}, + "soliket.BandPass": {"stop_at_error": True}, + }, "params": {**mflike_params}, } info1 = { "likelihood": {"soliket.mflike.TestMFLike": mflike_options}, - "theory": {"camb": camb_options, - "soliket.TheoryForge_MFLike": {'stop_at_error': True}, - "soliket.Foreground": {"stop_at_error": True}, - "soliket.BandPass": {"stop_at_error": True}}, + "theory": { + "camb": camb_options, + "soliket.TheoryForge_MFLike": {"stop_at_error": True}, + "soliket.Foreground": {"stop_at_error": True}, + "soliket.BandPass": {"stop_at_error": True}, + }, "params": {**mflike_params}, } diff --git a/soliket/tests/test_poisson.py b/soliket/tests/test_poisson.py index 2cfec387..9a2866a4 100644 --- a/soliket/tests/test_poisson.py +++ b/soliket/tests/test_poisson.py @@ -10,13 +10,12 @@ def rate_density(x, a): - """simple linear rate density - """ + """simple linear rate density""" return a * x def n_expected(a): - return 0.5 * a * (x_max ** 2 - x_min ** 2) # integral(rate_density, x_min, x_max) + return 0.5 * a * (x_max**2 - x_min**2) # integral(rate_density, x_min, x_max) def generate_data(a, with_samples=False, unc=0.3, Nk=64): @@ -27,7 +26,7 @@ def generate_data(a, with_samples=False, unc=0.3, Nk=64): u = np.random.random(n) # From inverting CDF of above normalized density - x = np.sqrt(u * (x_max ** 2 - x_min ** 2) + x_min ** 2) + x = np.sqrt(u * (x_max**2 - x_min**2) + x_min**2) if not with_samples: return x @@ -49,8 +48,9 @@ def test_poisson_experiment(a_true=3, N=100, with_samples=False, Nk=64): data = PoissonData("toy_samples", catalog, ["x"], samples=samples) a_grid = np.arange(0.1, 10, 0.1) - lnl = [data.loglike(partial(rate_density, a=a), - n_expected(a)) for a in a_grid] + lnl = [ + data.loglike(partial(rate_density, a=a), n_expected(a)) for a in a_grid + ] a_maxlike = a_grid[np.argmax(lnl)] a_maxlikes.append(a_maxlike) diff --git a/soliket/tests/test_ps.py b/soliket/tests/test_ps.py index 34d159f8..1429bdd7 100644 --- a/soliket/tests/test_ps.py +++ b/soliket/tests/test_ps.py @@ -23,7 +23,7 @@ def initialize(self): x = np.arange(self.n) if self.cov is None: cov = make_spd_matrix(self.n) * self.off_diag_amp - cov += np.diag(np.ones(self.n) * self.sigma ** 2) + cov += np.diag(np.ones(self.n) * self.sigma**2) else: cov = self.cov @@ -72,7 +72,9 @@ def test_toy(): like2 = get_likelihood(lhood, info2) like3 = get_likelihood(lhood, info3) - assert np.isclose(multilike1.logp(), sum([likex.logp() for - likex in [like1, like2, like3]])) - assert not np.isclose(multilike2.logp(), sum([likex.logp() for - likex in [like1, like2, like3]])) + assert np.isclose( + multilike1.logp(), sum([likex.logp() for likex in [like1, like2, like3]]) + ) + assert not np.isclose( + multilike2.logp(), sum([likex.logp() for likex in [like1, like2, like3]]) + ) diff --git a/soliket/tests/test_runs.py b/soliket/tests/test_runs.py index e3cd5c70..43b982e5 100644 --- a/soliket/tests/test_runs.py +++ b/soliket/tests/test_runs.py @@ -8,41 +8,49 @@ packages_path = resolve_packages_path() -@pytest.mark.parametrize("lhood", - ["mflike", - "lensing", - "lensing_lite", - "multi", - # "galaxykappa", - # "shearkappa" - # "xcorr" - ]) +@pytest.mark.parametrize( + "lhood", + [ + "mflike", + "lensing", + "lensing_lite", + "multi", + # "galaxykappa", + # "shearkappa" + # "xcorr" + ], +) def test_evaluate(lhood): info = yaml_load(pkgutil.get_data("soliket", f"tests/test_{lhood}.yaml")) info["force"] = True - info['sampler'] = {'evaluate': {}} + info["sampler"] = {"evaluate": {}} from cobaya.install import install + install(info, path=packages_path, skip_global=True, no_set_global=True) updated_info, sampler = run(info) -@pytest.mark.parametrize("lhood", - ["mflike", - "lensing", - "lensing_lite", - "multi", - # "galaxykappa", - # "shearkappa" - # "xcorr" - ]) +@pytest.mark.parametrize( + "lhood", + [ + "mflike", + "lensing", + "lensing_lite", + "multi", + # "galaxykappa", + # "shearkappa" + # "xcorr" + ], +) def test_mcmc(lhood): info = yaml_load(pkgutil.get_data("soliket", f"tests/test_{lhood}.yaml")) info["force"] = True - info['sampler'] = {'mcmc': {'max_samples': 10, 'max_tries': 1000}} + info["sampler"] = {"mcmc": {"max_samples": 10, "max_tries": 1000}} from cobaya.install import install + install(info, path=packages_path, skip_global=True, no_set_global=True) updated_info, sampler = run(info) diff --git a/soliket/tests/test_xcorr.py b/soliket/tests/test_xcorr.py index 77df8c12..a13927c7 100644 --- a/soliket/tests/test_xcorr.py +++ b/soliket/tests/test_xcorr.py @@ -77,42 +77,42 @@ def get_demo_xcorr_model(theory): @pytest.mark.skip(reason="Under development") -@pytest.mark.parametrize("theory", ["camb"])# , "classy"]) +@pytest.mark.parametrize("theory", ["camb"]) # , "classy"]) def test_xcorr(theory): - - params = {'b1': 1.0, 's1': 0.4} + params = {"b1": 1.0, "s1": 0.4} model = get_demo_xcorr_model(theory) lnl = model.loglike(params)[0] assert np.isfinite(lnl) - xcorr_lhood = model.likelihood['soliket.XcorrLikelihood'] + xcorr_lhood = model.likelihood["soliket.XcorrLikelihood"] setup_chi_out = xcorr_lhood._setup_chi() - Pk_interpolator = xcorr_lhood.provider.get_Pk_interpolator(("delta_nonu", - "delta_nonu"), - extrap_kmax=1.e8, - nonlinear=False).P + Pk_interpolator = xcorr_lhood.provider.get_Pk_interpolator( + ("delta_nonu", "delta_nonu"), extrap_kmax=1.0e8, nonlinear=False + ).P from soliket.xcorr.limber import do_limber - cl_gg, cl_kappag = do_limber(xcorr_lhood.ell_range, - xcorr_lhood.provider, - xcorr_lhood.dndz, - xcorr_lhood.dndz, - params['s1'], - params['s1'], - Pk_interpolator, - params['b1'], - params['b1'], - xcorr_lhood.alpha_auto, - xcorr_lhood.alpha_cross, - setup_chi_out, - Nchi=xcorr_lhood.Nchi, - dndz1_mag=xcorr_lhood.dndz, - dndz2_mag=xcorr_lhood.dndz) + cl_gg, cl_kappag = do_limber( + xcorr_lhood.ell_range, + xcorr_lhood.provider, + xcorr_lhood.dndz, + xcorr_lhood.dndz, + params["s1"], + params["s1"], + Pk_interpolator, + params["b1"], + params["b1"], + xcorr_lhood.alpha_auto, + xcorr_lhood.alpha_cross, + setup_chi_out, + Nchi=xcorr_lhood.Nchi, + dndz1_mag=xcorr_lhood.dndz, + dndz2_mag=xcorr_lhood.dndz, + ) ell_load = xcorr_lhood.data.x cl_load = xcorr_lhood.data.y @@ -128,36 +128,48 @@ def test_xcorr(theory): # Nell_unwise_g = np.ones_like(cl_gg) \ # / (xcorr_lhood.ngal * (60 * 180 / np.pi)**2) - Nell_obs_unwise_g = np.ones_like(cl_obs_gg) \ - / (xcorr_lhood.ngal * (60 * 180 / np.pi)**2) + Nell_obs_unwise_g = np.ones_like(cl_obs_gg) / ( + xcorr_lhood.ngal * (60 * 180 / np.pi) ** 2 + ) import pyccl as ccl - h2 = (xcorr_lhood.provider.get_param('H0') / 100)**2 - - cosmo = ccl.Cosmology(Omega_c=xcorr_lhood.provider.get_param('omch2') / h2, - Omega_b=xcorr_lhood.provider.get_param('ombh2') / h2, - h=xcorr_lhood.provider.get_param('H0') / 100, - n_s=xcorr_lhood.provider.get_param('ns'), - A_s=xcorr_lhood.provider.get_param('As'), - Omega_k=xcorr_lhood.provider.get_param('omk'), - Neff=xcorr_lhood.provider.get_param('nnu'), - matter_power_spectrum='linear') - - g_bias_zbz = (xcorr_lhood.dndz[:, 0], - params['b1'] * np.ones(len(xcorr_lhood.dndz[:, 0]))) - mag_bias_zbz = (xcorr_lhood.dndz[:, 0], - params['s1'] * np.ones(len(xcorr_lhood.dndz[:, 0]))) - - tracer_g = ccl.NumberCountsTracer(cosmo, - has_rsd=False, - dndz=xcorr_lhood.dndz.T, - bias=g_bias_zbz, - mag_bias=mag_bias_zbz) + + h2 = (xcorr_lhood.provider.get_param("H0") / 100) ** 2 + + cosmo = ccl.Cosmology( + Omega_c=xcorr_lhood.provider.get_param("omch2") / h2, + Omega_b=xcorr_lhood.provider.get_param("ombh2") / h2, + h=xcorr_lhood.provider.get_param("H0") / 100, + n_s=xcorr_lhood.provider.get_param("ns"), + A_s=xcorr_lhood.provider.get_param("As"), + Omega_k=xcorr_lhood.provider.get_param("omk"), + Neff=xcorr_lhood.provider.get_param("nnu"), + matter_power_spectrum="linear", + ) + + g_bias_zbz = ( + xcorr_lhood.dndz[:, 0], + params["b1"] * np.ones(len(xcorr_lhood.dndz[:, 0])), + ) + mag_bias_zbz = ( + xcorr_lhood.dndz[:, 0], + params["s1"] * np.ones(len(xcorr_lhood.dndz[:, 0])), + ) + + tracer_g = ccl.NumberCountsTracer( + cosmo, + has_rsd=False, + dndz=xcorr_lhood.dndz.T, + bias=g_bias_zbz, + mag_bias=mag_bias_zbz, + ) tracer_k = ccl.CMBLensingTracer(cosmo, z_source=1100) cl_gg_ccl = ccl.cells.angular_cl(cosmo, tracer_g, tracer_g, xcorr_lhood.ell_range) - cl_kappag_ccl = ccl.cells.angular_cl(cosmo, tracer_k, tracer_g, xcorr_lhood.ell_range) + cl_kappag_ccl = ccl.cells.angular_cl( + cosmo, tracer_k, tracer_g, xcorr_lhood.ell_range + ) assert np.allclose(cl_gg_ccl, cl_gg) assert np.allclose(cl_kappag_ccl, cl_kappag) diff --git a/soliket/utils.py b/soliket/utils.py index de923fb1..65491d07 100644 --- a/soliket/utils.py +++ b/soliket/utils.py @@ -65,8 +65,12 @@ class OneWithCls(one): lmax = 10000 def get_requirements(self): - return {"Cl": {"pp": self.lmax, - "tt": self.lmax, - "te": self.lmax, - "ee": self.lmax, - "bb": self.lmax, }} + return { + "Cl": { + "pp": self.lmax, + "tt": self.lmax, + "te": self.lmax, + "ee": self.lmax, + "bb": self.lmax, + } + } diff --git a/soliket/xcorr/limber.py b/soliket/xcorr/limber.py index da8f3fef..0e9ff9a9 100644 --- a/soliket/xcorr/limber.py +++ b/soliket/xcorr/limber.py @@ -13,7 +13,7 @@ from numpy import trapz as trapezoid from soliket.constants import C_HMPC -oneover_chmpc = 1. / C_HMPC +oneover_chmpc = 1.0 / C_HMPC def mag_bias_kernel(provider, dndz, s1, zatchi, chi_arr, chiprime_arr, zprime_arr): @@ -23,49 +23,92 @@ def mag_bias_kernel(provider, dndz, s1, zatchi, chi_arr, chiprime_arr, zprime_ar norm = trapezoid(dndz[:, 1], x=dndz[:, 0]) dndzprime = dndzprime / norm # TODO check this norm is right - g_integrand = (chiprime_arr - chi_arr[np.newaxis, :]) / chiprime_arr \ - * (oneover_chmpc * provider.get_param('H0') / 100) \ - * np.sqrt(provider.get_param('omegam') * (1 + zprime_arr) ** 3. - + 1 - provider.get_param('omegam')) \ - * dndzprime + g_integrand = ( + (chiprime_arr - chi_arr[np.newaxis, :]) + / chiprime_arr + * (oneover_chmpc * provider.get_param("H0") / 100) + * np.sqrt( + provider.get_param("omegam") * (1 + zprime_arr) ** 3.0 + + 1 + - provider.get_param("omegam") + ) + * dndzprime + ) g = chi_arr * trapezoid(g_integrand, x=chiprime_arr, axis=0) - W_mu = (5. * s1 - 2.) * 1.5 * provider.get_param('omegam') \ - * (provider.get_param('H0') / 100) ** 2 * oneover_chmpc ** 2 \ - * (1. + zatchi(chi_arr)) * g + W_mu = ( + (5.0 * s1 - 2.0) + * 1.5 + * provider.get_param("omegam") + * (provider.get_param("H0") / 100) ** 2 + * oneover_chmpc**2 + * (1.0 + zatchi(chi_arr)) + * g + ) return W_mu -def do_limber(ell_arr, provider, dndz1, dndz2, s1, s2, pk, b1_HF, b2_HF, - alpha_auto, alpha_cross, - chi_grids, - # use_zeff=True, - Nchi=50, dndz1_mag=None, dndz2_mag=None, normed=False): - zatchi = chi_grids['zatchi'] +def do_limber( + ell_arr, + provider, + dndz1, + dndz2, + s1, + s2, + pk, + b1_HF, + b2_HF, + alpha_auto, + alpha_cross, + chi_grids, + # use_zeff=True, + Nchi=50, + dndz1_mag=None, + dndz2_mag=None, + normed=False, +): + zatchi = chi_grids["zatchi"] # chiatz = chi_grids['chiatz'] - chi_arr = chi_grids['chival'] + chi_arr = chi_grids["chival"] # z_arr = chi_grids['zval'] - chiprime_arr = chi_grids['chivalp'] - zprime_arr = chi_grids['zvalp'] + chiprime_arr = chi_grids["chivalp"] + zprime_arr = chi_grids["zvalp"] - chistar = provider.get_comoving_radial_distance(provider.get_param('zstar')) + chistar = provider.get_comoving_radial_distance(provider.get_param("zstar")) # Galaxy kernels, assumed to be b(z) * dN/dz - W_g1 = np.interp(zatchi(chi_arr), dndz1[:, 0], dndz1[:, 1] \ - * provider.get_Hubble(dndz1[:, 0], units='1/Mpc'), left=0, right=0) + W_g1 = np.interp( + zatchi(chi_arr), + dndz1[:, 0], + dndz1[:, 1] * provider.get_Hubble(dndz1[:, 0], units="1/Mpc"), + left=0, + right=0, + ) if not normed: W_g1 /= trapezoid(W_g1, x=chi_arr) - W_g2 = np.interp(zatchi(chi_arr), dndz2[:, 0], dndz2[:, 1] \ - * provider.get_Hubble(dndz2[:, 0], units='1/Mpc'), left=0, right=0) + W_g2 = np.interp( + zatchi(chi_arr), + dndz2[:, 0], + dndz2[:, 1] * provider.get_Hubble(dndz2[:, 0], units="1/Mpc"), + left=0, + right=0, + ) if not normed: W_g2 /= trapezoid(W_g2, x=chi_arr) - W_kappa = oneover_chmpc ** 2. * 1.5 * provider.get_param('omegam') \ - * (provider.get_param('H0') / 100) ** 2. * (1. + zatchi(chi_arr)) \ - * chi_arr * (chistar - chi_arr) / chistar + W_kappa = ( + oneover_chmpc**2.0 + * 1.5 + * provider.get_param("omegam") + * (provider.get_param("H0") / 100) ** 2.0 + * (1.0 + zatchi(chi_arr)) + * chi_arr + * (chistar - chi_arr) + / chistar + ) # Get effective redshift # if use_zeff: @@ -75,8 +118,9 @@ def do_limber(ell_arr, provider, dndz1, dndz2, s1, s2, pk, b1_HF, b2_HF, # zeff = -1.0 # set up magnification bias kernels - W_mu1 = mag_bias_kernel(provider, dndz1, s1, - zatchi, chi_arr, chiprime_arr, zprime_arr) + W_mu1 = mag_bias_kernel( + provider, dndz1, s1, zatchi, chi_arr, chiprime_arr, zprime_arr + ) c_ell_g1g1 = np.zeros([ell_arr.shape[0], 1, chi_arr.shape[0]]) c_ell_g1kappa = np.zeros([ell_arr.shape[0], 1, chi_arr.shape[0]]) @@ -94,22 +138,22 @@ def do_limber(ell_arr, provider, dndz1, dndz2, s1, s2, pk, b1_HF, b2_HF, p_gg = b1_HF * b1_HF * p_mm_hf # lets just stay at constant linear bias for now p_gm = b1_HF * p_mm_hf - W_g1g1 = W_g1[i_chi] * W_g1[i_chi] / (chi ** 2) * p_gg + W_g1g1 = W_g1[i_chi] * W_g1[i_chi] / (chi**2) * p_gg c_ell_g1g1[:, :, i_chi] = W_g1g1.T - W_g1kappa = W_g1[i_chi] * W_kappa[i_chi] / (chi ** 2) * p_gm + W_g1kappa = W_g1[i_chi] * W_kappa[i_chi] / (chi**2) * p_gm c_ell_g1kappa[:, :, i_chi] = W_g1kappa.T # W_kappakappa = W_kappa[i_chi] * W_kappa[i_chi] / (chi**2) * p_mm # c_ell_kappakappa[:,:,i_chi] = W_kappakappa.T - W_g1mu1 = W_g1[i_chi] * W_mu1[i_chi] / (chi ** 2) * p_gm + W_g1mu1 = W_g1[i_chi] * W_mu1[i_chi] / (chi**2) * p_gm c_ell_g1mu1[:, :, i_chi] = W_g1mu1.T - W_mu1mu1 = W_mu1[i_chi] * W_mu1[i_chi] / (chi ** 2) * p_mm + W_mu1mu1 = W_mu1[i_chi] * W_mu1[i_chi] / (chi**2) * p_mm c_ell_mu1mu1[:, :, i_chi] = W_mu1mu1.T - W_mu1kappa = W_kappa[i_chi] * W_mu1[i_chi] / (chi ** 2) * p_mm + W_mu1kappa = W_kappa[i_chi] * W_mu1[i_chi] / (chi**2) * p_mm c_ell_mu1kappa[:, :, i_chi] = W_mu1kappa.T c_ell_g1g1 = trapezoid(c_ell_g1g1, x=chi_arr, axis=-1) @@ -120,7 +164,7 @@ def do_limber(ell_arr, provider, dndz1, dndz2, s1, s2, pk, b1_HF, b2_HF, c_ell_mu1mu1 = trapezoid(c_ell_mu1mu1, x=chi_arr, axis=-1) c_ell_mu1kappa = trapezoid(c_ell_mu1kappa, x=chi_arr, axis=-1) - clobs_gg = c_ell_g1g1 + 2. * c_ell_g1mu1 + c_ell_mu1mu1 + clobs_gg = c_ell_g1g1 + 2.0 * c_ell_g1mu1 + c_ell_mu1mu1 clobs_kappag = c_ell_g1kappa + c_ell_mu1kappa # clobs_kappakappa = c_ell_kappakappa diff --git a/soliket/xcorr/xcorr.py b/soliket/xcorr/xcorr.py index 6134119c..49c5e6b0 100644 --- a/soliket/xcorr/xcorr.py +++ b/soliket/xcorr/xcorr.py @@ -62,10 +62,9 @@ class XcorrLikelihood(GaussianLikelihood): def initialize(self): name: str = "Xcorr" # noqa F841 - self.log.info('Initialising.') + self.log.info("Initialising.") if self.datapath is None: - dndz_file: Optional[str] # noqa F821 auto_file: Optional[str] # noqa F821 cross_file: Optional[str] # noqa F821 @@ -74,34 +73,33 @@ def initialize(self): self.x, self.y, self.dy = self._get_data() if self.covpath is None: - self.log.info('No xcorr covariance specified. Using diag(dy^2).') + self.log.info("No xcorr covariance specified. Using diag(dy^2).") self.cov = np.diag(self.dy**2) else: self.cov = self._get_cov() else: - self.k_tracer_name: Optional[str] # noqa F821 self.gc_tracer_name: Optional[str] # noqa F821 # tracer_combinations: Optional[str] # TODO: implement with keep_selection self.sacc_data = self._get_sacc_data() - self.x = self.sacc_data['x'] - self.y = self.sacc_data['y'] - self.cov = self.sacc_data['cov'] - self.dndz = self.sacc_data['dndz'] - self.ngal = self.sacc_data['ngal'] + self.x = self.sacc_data["x"] + self.y = self.sacc_data["y"] + self.cov = self.sacc_data["cov"] + self.dndz = self.sacc_data["dndz"] + self.ngal = self.sacc_data["ngal"] # TODO is this resolution limit on zarray a CAMB problem? self.nz: Optional[int] # noqa F821 assert self.nz <= 149, "CAMB limitations requires nz <= 149" self.zarray = np.linspace(self.dndz[:, 0].min(), self.dndz[:, 0].max(), self.nz) - self.zbgdarray = np.concatenate([self.zarray, [1100]]) # TODO: unfix zstar + self.zbgdarray = np.concatenate([self.zarray, [1100]]) # TODO: unfix zstar self.Nchi: Optional[int] # noqa F821 self.Nchi_mag: Optional[int] # noqa F821 - #self.use_zeff: Optional[bool] # noqa F821 + # self.use_zeff: Optional[bool] # noqa F821 self.Pk_interp_kmax: Optional[float] # noqa F821 @@ -114,73 +112,69 @@ def initialize(self): self.data = GaussianData(self.name, self.x, self.y, self.cov) - def get_requirements(self): return { - 'Cl': {'lmax': self.high_ell, - 'pp': self.high_ell}, - "Pk_interpolator": { - "z": self.zarray[:-1], - "k_max": self.Pk_interp_kmax, - #"extrap_kmax": 20.0, - "nonlinear": False, - "hubble_units": False, # cobaya told me to - "k_hunit": False, # cobaya told me to - "vars_pairs": [["delta_nonu", "delta_nonu"]], - }, - "Hubble": {"z": self.zarray}, - "angular_diameter_distance": {"z": self.zbgdarray}, - "comoving_radial_distance": {"z": self.zbgdarray}, - 'H0': None, - 'ombh2': None, - 'omch2': None, - 'omk': None, - 'omegam': None, - 'zstar': None, - 'As': None, - 'ns': None - } + "Cl": {"lmax": self.high_ell, "pp": self.high_ell}, + "Pk_interpolator": { + "z": self.zarray[:-1], + "k_max": self.Pk_interp_kmax, + # "extrap_kmax": 20.0, + "nonlinear": False, + "hubble_units": False, # cobaya told me to + "k_hunit": False, # cobaya told me to + "vars_pairs": [["delta_nonu", "delta_nonu"]], + }, + "Hubble": {"z": self.zarray}, + "angular_diameter_distance": {"z": self.zbgdarray}, + "comoving_radial_distance": {"z": self.zbgdarray}, + "H0": None, + "ombh2": None, + "omch2": None, + "omk": None, + "omegam": None, + "zstar": None, + "As": None, + "ns": None, + } def _bin(self, theory_cl, lmin, lmax): binned_theory_cl = np.zeros_like(lmin) for i in range(len(lmin)): - binned_theory_cl[i] = np.mean(theory_cl[(self.ell_range >= lmin[i]) - & (self.ell_range < lmax[i])]) + binned_theory_cl[i] = np.mean( + theory_cl[(self.ell_range >= lmin[i]) & (self.ell_range < lmax[i])] + ) return binned_theory_cl def _get_sacc_data(self, **params_values): - data_sacc = sacc.Sacc.load_fits(self.datapath) # TODO: would be better to use keep_selection data_sacc.remove_selection(tracers=(self.k_tracer_name, self.k_tracer_name)) - ell_auto, cl_auto = data_sacc.get_ell_cl('cl_00', - self.gc_tracer_name, - self.gc_tracer_name) - ell_cross, cl_cross = data_sacc.get_ell_cl('cl_00', - self.gc_tracer_name, - self.k_tracer_name) # TODO: check order + ell_auto, cl_auto = data_sacc.get_ell_cl( + "cl_00", self.gc_tracer_name, self.gc_tracer_name + ) + ell_cross, cl_cross = data_sacc.get_ell_cl( + "cl_00", self.gc_tracer_name, self.k_tracer_name + ) # TODO: check order cov = data_sacc.covariance.covmat x = np.concatenate([ell_auto, ell_cross]) y = np.concatenate([cl_auto, cl_cross]) - dndz = np.column_stack([data_sacc.tracers[self.gc_tracer_name].z, - data_sacc.tracers[self.gc_tracer_name].nz]) - ngal = data_sacc.tracers[self.gc_tracer_name].metadata['ngal'] + dndz = np.column_stack( + [ + data_sacc.tracers[self.gc_tracer_name].z, + data_sacc.tracers[self.gc_tracer_name].nz, + ] + ) + ngal = data_sacc.tracers[self.gc_tracer_name].metadata["ngal"] - data = {'x': x, - 'y': y, - 'cov': cov, - 'dndz': dndz, - 'ngal': ngal} + data = {"x": x, "y": y, "cov": cov, "dndz": dndz, "ngal": ngal} return data - def _get_data(self, **params_values): - data_auto = np.loadtxt(self.auto_file) data_cross = np.loadtxt(self.cross_file) @@ -200,55 +194,59 @@ def _get_data(self, **params_values): return x, y, dy def _setup_chi(self): - chival = self.provider.get_comoving_radial_distance(self.zarray) zatchi = Spline(chival, self.zarray) chiatz = Spline(self.zarray, chival) - chimin = np.min(chival) + 1.e-5 + chimin = np.min(chival) + 1.0e-5 chimax = np.max(chival) chival = np.linspace(chimin, chimax, self.Nchi) zval = zatchi(chival) - chistar = \ - self.provider.get_comoving_radial_distance(self.provider.get_param('zstar')) - chivalp = \ - np.array(list(map(lambda x: np.linspace(x, chistar, self.Nchi_mag), chival))) + chistar = self.provider.get_comoving_radial_distance( + self.provider.get_param("zstar") + ) + chivalp = np.array( + list(map(lambda x: np.linspace(x, chistar, self.Nchi_mag), chival)) + ) chivalp = chivalp.transpose()[0] zvalp = zatchi(chivalp) - chi_result = {'zatchi': zatchi, - 'chiatz': chiatz, - 'chival': chival, - 'zval': zval, - 'chivalp': chivalp, - 'zvalp': zvalp} + chi_result = { + "zatchi": zatchi, + "chiatz": chiatz, + "chival": chival, + "zval": zval, + "chivalp": chivalp, + "zvalp": zvalp, + } return chi_result def _get_theory(self, **params_values): - setup_chi_out = self._setup_chi() - Pk_interpolator = self.provider.get_Pk_interpolator(("delta_nonu", "delta_nonu"), - extrap_kmax=1.e8, - nonlinear=False).P - - cl_gg, cl_kappag = do_limber(self.ell_range, - self.provider, - self.dndz, - self.dndz, - params_values['s1'], - params_values['s1'], - Pk_interpolator, - params_values['b1'], - params_values['b1'], - self.alpha_auto, - self.alpha_cross, - setup_chi_out, - Nchi=self.Nchi, - #use_zeff=self.use_zeff, - dndz1_mag=self.dndz, - dndz2_mag=self.dndz) + Pk_interpolator = self.provider.get_Pk_interpolator( + ("delta_nonu", "delta_nonu"), extrap_kmax=1.0e8, nonlinear=False + ).P + + cl_gg, cl_kappag = do_limber( + self.ell_range, + self.provider, + self.dndz, + self.dndz, + params_values["s1"], + params_values["s1"], + Pk_interpolator, + params_values["b1"], + params_values["b1"], + self.alpha_auto, + self.alpha_cross, + setup_chi_out, + Nchi=self.Nchi, + # use_zeff=self.use_zeff, + dndz1_mag=self.dndz, + dndz2_mag=self.dndz, + ) # TODO: this is not the correct binning, # but there needs to be a consistent way to specify it @@ -256,6 +254,6 @@ def _get_theory(self, **params_values): ell_gg, clobs_gg = utils.binner(self.ell_range, cl_gg, bin_edges) ell_kappag, clobs_kappag = utils.binner(self.ell_range, cl_kappag, bin_edges) - #ell_kappakappa, clobs_kappakappa = utils.binner(self.ell_range, cl_kappakappa, bin_edges) # noqa E501 + # ell_kappakappa, clobs_kappakappa = utils.binner(self.ell_range, cl_kappakappa, bin_edges) # noqa E501 return np.concatenate([clobs_gg, clobs_kappag])