diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 76c6e51..4727120 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -20,34 +20,4 @@ jobs: - name: Upload to Codecov uses: codecov/codecov-action@v5 with: - token: ${{ secrets.CODECOV_TOKEN }} - - name: coverage.py badge - uses: tj-actions/coverage-badge-py@v2.0.3 - - name: Verify Changed files - uses: tj-actions/verify-changed-files@v20 - id: verify-changed-files - with: - files: coverage.svg - - name: Commit files - if: steps.verify-changed-files.outputs.files_changed == 'true' - run: | - git config --local user.email "github-actions[bot]@users.noreply.github.com" - git config --local user.name "github-actions[bot]" - git add coverage.svg - git commit -m "Updated coverage.svg" - - name: Push changes - if: steps.verify-changed-files.outputs.files_changed == 'true' - uses: ad-m/github-push-action@master - with: - github_token: ${{ secrets.github_token }} - branch: ${{ github.ref }} - - name: Re-pull on failure - if: ${{ failure() }} - run: git pull origin ${{ github.ref }} --autostash --rebase -X ours - - name: Re-push on failure - if: ${{ failure() }} - uses: ad-m/github-push-action@master - with: - branch: ${{ github.ref }} - force: true - github_token: ${{ secrets.GITHUB_TOKEN }} \ No newline at end of file + token: ${{ secrets.CODECOV_TOKEN }} \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md index c5f03b8..88c29d1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -34,4 +34,5 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 1.2.1 : Updating configs of other surveys: SKAO, DESI, LSST to work in new config file structure 1.2.2 : Galaxy sample split for sources and lenses. Feedback prints more consistent. 1.2.3 : Added a new likelihood function for photometric and spectroscopic surveys. -1.2.4 : Discontinued support for python3.8. Fixed the style for likelihood scripts. \ No newline at end of file +1.2.4 : Discontinued support for python3.8. Fixed the style for likelihood scripts. +1.2.5 : Likelihood modules tested for photo and spectro diff --git a/README.md b/README.md index c9c74e5..eae33e1 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,10 @@ # cosmicfishpie [![Documentation Status](https://readthedocs.org/projects/cosmicfishpie/badge/?version=latest)](https://cosmicfishpie.readthedocs.io/en/latest/?badge=latest) -![Coverage](./coverage.svg) [![codecov](https://codecov.io/github/santiagocasas/cosmicfishpie/graph/badge.svg?token=BXTVDPXPUO)](https://codecov.io/github/santiagocasas/cosmicfishpie) +[![codecov](https://codecov.io/github/santiagocasas/cosmicfishpie/graph/badge.svg?token=BXTVDPXPUO)](https://codecov.io/github/santiagocasas/cosmicfishpie) +[![37.50 % FAIR](https://img.shields.io/badge/FAIR_assessment-37.50_%25-red)](https://fair-checker.france-bioinformatique.fr/assessment/68da9e3cc49e421b3e2cf501) - -See CITATION.cff for citation metadata. +See [`CITATION.cff`](CITATION.cff) for citation metadata or click on the "Cite this repository" button on the sidebar.
diff --git a/cosmicfishpie/analysis/fishconsumer.py b/cosmicfishpie/analysis/fishconsumer.py index e45951a..e64bb85 100644 --- a/cosmicfishpie/analysis/fishconsumer.py +++ b/cosmicfishpie/analysis/fishconsumer.py @@ -204,7 +204,7 @@ def fishtable_to_pandas( if ff.name in filter_names: newname = ff.name.split(" ")[0] # print(newname) - ltxname = "$\\mathrm{" + str(newname).replace(" ", "\ ") + "}$" + ltxname = "$\\mathrm{" + str(newname).replace(" ", "\\ ") + "}$" barplot_data[ltxname] = rel_err if return_data_bar: return barplot_data @@ -586,7 +586,7 @@ def prepare_settings_plot( .replace(r"\,", "-") .replace("$", "") .replace(r"\\", "") - .replace("\mathrm", "") + .replace("\\mathrm", "") .replace("___", "_") .replace("__", "_") .replace("----", "-") diff --git a/cosmicfishpie/likelihood/__init__.py b/cosmicfishpie/likelihood/__init__.py new file mode 100644 index 0000000..5541017 --- /dev/null +++ b/cosmicfishpie/likelihood/__init__.py @@ -0,0 +1,7 @@ +""" +Likelihood subpackage: tools to compute data/theory spectra and chi2. +""" + +from .base import Likelihood # noqa: F401 +from .photo_like import PhotometricLikelihood # noqa: F401 +from .spectro_like import SpectroLikelihood # noqa: F401 diff --git a/cosmicfishpie/likelihood/base.py b/cosmicfishpie/likelihood/base.py new file mode 100644 index 0000000..a514d81 --- /dev/null +++ b/cosmicfishpie/likelihood/base.py @@ -0,0 +1,220 @@ +"""Base infrastructure for likelihood modules in Cosmicfishpie. + +This module provides the core interface and utilities for implementing +likelihood calculations in Cosmicfishpie, particularly for cosmological +parameter estimation using Fisher matrices. +""" + +from __future__ import annotations + +from abc import ABC, abstractmethod +from collections.abc import Sequence +from typing import Any, Dict, Iterable, Optional + +import numpy as np +from nautilus import Prior, Sampler + +from cosmicfishpie.fishermatrix.cosmicfish import FisherMatrix + + +def is_indexable_iterable(var: Any) -> bool: + """Check if a variable is an indexable iterable. + + Args: + var: The variable to check + + Returns: + bool: True if the variable is an indexable iterable (list, numpy array, or Sequence) + and not a string or bytes object, False otherwise. + """ + return isinstance(var, (list, np.ndarray, Sequence)) and not isinstance(var, (str, bytes)) + + +class Likelihood(ABC): + """Common interface for likelihood evaluations used in Cosmicfishpie. + + This abstract base class defines the core interface for likelihood calculations + in Cosmicfishpie. It provides methods for computing data representations, + theory predictions, and chi-squared values, as well as utilities for parameter + handling and running Nautilus samplers. + """ + + def __init__( + self, + *, + cosmo_data: FisherMatrix, + cosmo_theory: Optional[FisherMatrix] = None, + leg_flag: str = "wedges", + ) -> None: + """Initialize the Likelihood object with Fisher matrices. + + Args: + cosmo_data: FisherMatrix object containing the observed data + cosmo_theory: Optional FisherMatrix object for theory predictions. + If None, uses cosmoFM_data + leg_flag: Flag indicating the type of data representation to use + ("wedges" or other supported types) + """ + self.cosmo_data = cosmo_data + self.cosmo_theory = cosmo_theory or cosmo_data + self.leg_flag = leg_flag + self.data_obs = self.compute_data() + + @abstractmethod + def compute_data(self) -> Any: + """Compute and return the observed data representation. + + This method should be implemented by subclasses to return the data + in the appropriate format (e.g., wedges or multipoles) for likelihood + calculations. + + Returns: + The observed data representation + """ + + @abstractmethod + def compute_theory(self, param_dict: Dict[str, Any]) -> Any: + """Compute and return the theory prediction. + + This method should be implemented by subclasses to return the theory + prediction in the same representation as the observed data. + + Args: + param_dict: Dictionary of cosmological parameters + + Returns: + The theory prediction in the same format as the observed data + """ + + @abstractmethod + def compute_chi2(self, theory_obs: Any) -> float: + """Compute and return the chi-squared value. + + This method should be implemented by subclasses to compute the chi-squared + value between the observed data and theory prediction. + + Args: + theory_obs: Theory prediction in the same format as observed data + + Returns: + The chi-squared value + """ + + def build_param_dict( + self, + *, + param_vec: Optional[Iterable[float]] = None, + param_dict: Optional[Dict[str, Any]] = None, + prior: Optional[Any] = None, + ) -> Dict[str, Any]: + """Create a parameter dictionary from vector inputs when needed. + + This method converts a parameter vector to a dictionary using the prior + information, or returns the provided parameter dictionary directly. + + Args: + param_vec: Optional iterable of parameter values + param_dict: Optional dictionary of parameters + prior: Optional prior information for parameter mapping + + Returns: + Dictionary of parameters + + Raises: + ValueError: If neither param_dict nor (param_vec and prior) are provided + TypeError: If param_vec is not an indexable iterable + AttributeError: If prior object doesn't expose an ordered 'keys' attribute + """ + if param_dict is not None: + return dict(param_dict) + + if param_vec is None or prior is None: + raise ValueError( + "Provide either param_dict or (param_vec and prior) to build the parameter mapping" + ) + + if not is_indexable_iterable(param_vec): + raise TypeError( + "param_vec must be an indexable iterable when no param_dict is supplied" + ) + + prior_keys = getattr(prior, "keys", None) + if prior_keys is None: + raise AttributeError( + "prior object must expose an ordered 'keys' attribute to map the vector to parameters" + ) + + return {key: param_vec[i] for i, key in enumerate(prior_keys)} + + def loglike( + self, + param_vec: Optional[Iterable[float]] = None, + *, + param_dict: Optional[Dict[str, Any]] = None, + prior: Optional[Any] = None, + ) -> float: + """Compute the log-likelihood value. + + This method computes the log-likelihood value (-0.5 * χ²) for the + supplied parameters. It can accept either a parameter vector or a + parameter dictionary. + + Args: + param_vec: Optional iterable of parameter values + param_dict: Optional dictionary of parameters + prior: Optional prior information for parameter mapping + + Returns: + The log-likelihood value (-0.5 * χ²) + """ + params = self.build_param_dict(param_vec=param_vec, param_dict=param_dict, prior=prior) + theory_obs = self.compute_theory(params) + chi2 = self.compute_chi2(theory_obs) + return -0.5 * chi2 + + +class NautilusMixin: + """Mixin class for running Nautilus samplers.""" + + def create_nautilus_prior(self, prior_dict: Dict[str, Any]) -> Prior: + """Create a Nautilus prior object from a dictionary of parameter names and their prior ranges. + + Args: + prior_dict: Dictionary of parameter names and their prior ranges + + Returns: + Nautilus Prior object + """ + + prior = Prior() + for par, (lower, upper) in prior_dict.items(): + prior.add_parameter(par, (lower, upper)) + return prior + + def run_nautilus( + self, + *, + prior: Any, + sampler_kwargs: Optional[Dict[str, Any]] = None, + run_kwargs: Optional[Dict[str, Any]] = None, + ) -> Sampler: + """Convenience wrapper to launch a Nautilus sampler using this likelihood. + + This method provides a convenient interface to run Nautilus samplers + with the current likelihood function. + + Args: + prior: Prior information for the parameters + sampler_kwargs: Optional dictionary of keyword arguments for the Sampler + run_kwargs: Optional dictionary of keyword arguments for the run method + + Returns: + The Nautilus Sampler object + """ + + sampler_kwargs = dict(sampler_kwargs or {}) + run_kwargs = dict(run_kwargs or {}) + + sampler = Sampler(prior, self.loglike, **sampler_kwargs, likelihood_kwargs={"prior": prior}) + sampler.run(**run_kwargs) + return sampler diff --git a/cosmicfishpie/likelihood/photo_like.py b/cosmicfishpie/likelihood/photo_like.py new file mode 100644 index 0000000..e5fa46c --- /dev/null +++ b/cosmicfishpie/likelihood/photo_like.py @@ -0,0 +1,256 @@ +"""Photometric likelihood implementation for Cosmicfishpie.""" + +from __future__ import annotations + +import logging +from copy import copy, deepcopy +from itertools import product +from typing import Any, Dict, Iterable, Optional + +import numpy as np + +from cosmicfishpie.LSSsurvey import photo_cov as pcov +from cosmicfishpie.LSSsurvey import photo_obs as pobs + +from .base import Likelihood, NautilusMixin + +logger = logging.getLogger("cosmicfishpie.likelihood.photo") + + +def _dict_with_updates(template: Dict[str, Any], pool: Dict[str, Any]) -> Dict[str, Any]: + updated = deepcopy(template) + for key in template: + if key in pool: + updated[key] = pool.pop(key) + return updated + + +def _cells_from_cls( + photo_cls: pobs.ComputeCls, + photo_cov: pcov.PhotoCov, + observables: Iterable[str], +) -> Dict[str, np.ndarray]: + photo_cls.compute_all() + + observables = tuple(observables) + ells = np.array(photo_cls.result["ells"], copy=True) + data: Dict[str, np.ndarray] = {"ells": ells} + + if "WL" in observables: + wl_bins = list(photo_cls.binrange_WL) + n_wl = len(wl_bins) + cell_ll = np.empty((len(ells), n_wl, n_wl), dtype=np.float64) + for i, j in product(wl_bins, repeat=2): + base = photo_cls.result[f"WL {i}xWL {j}"] + noise = 0.0 + if hasattr(photo_cov, "ellipt_error") and hasattr(photo_cov, "ngalbin_WL"): + noise = ( + (photo_cov.ellipt_error**2.0) / photo_cov.ngalbin_WL[i - 1] if i == j else 0.0 + ) + cell_ll[:, i - 1, j - 1] = base + noise + data["Cell_LL"] = cell_ll + else: + wl_bins = [] + + if "GCph" in observables: + gc_bins = list(photo_cls.binrange_GCph) + n_gc = len(gc_bins) + cell_gg = np.empty((len(ells), n_gc, n_gc), dtype=np.float64) + for i, j in product(gc_bins, repeat=2): + base = photo_cls.result[f"GCph {i}xGCph {j}"] + noise = 0.0 + if hasattr(photo_cov, "ngalbin_GCph"): + noise = (1.0 / photo_cov.ngalbin_GCph[i - 1]) if i == j else 0.0 + cell_gg[:, i - 1, j - 1] = base + noise + data["Cell_GG"] = cell_gg + else: + gc_bins = [] + + if "WL" in observables and "GCph" in observables: + cell_gl = np.empty((len(ells), len(gc_bins), len(wl_bins)), dtype=np.float64) + for i, j in product(gc_bins, wl_bins): + cell_gl[:, i - 1, j - 1] = photo_cls.result[f"GCph {i}xWL {j}"] + data["Cell_GL"] = cell_gl + + return data + + +def _chi2_per_obs( + cell_fid: np.ndarray, cell_th: np.ndarray, ells: np.ndarray, dells: np.ndarray +) -> float: + dfid = np.linalg.det(cell_fid) + dth = np.linalg.det(cell_th) + + dmix = 0.0 + for idx in range(cell_fid.shape[-1]): + mix = copy(cell_th) + mix[:, idx, :] = cell_fid[:, idx, :] + dmix += np.linalg.det(mix) + + integrand = (2 * ells + 1) * (dmix / dth + np.log(dth / dfid) - cell_fid.shape[-1]) + integrand = np.array(integrand, copy=False) + result = np.sum( + np.concatenate( + [((integrand[1:] + integrand[:-1]) / 2) * dells[:-1], integrand[-1:] * dells[-1:]] + ) + ) + return float(result) + + +class PhotometricLikelihood(Likelihood, NautilusMixin): + """Likelihood built from photometric clusterings (WL / GCph).""" + + def __init__( + self, + *, + cosmo_data, + cosmo_theory=None, + observables: Optional[Iterable[str]] = None, + data_cells: Optional[Dict[str, np.ndarray]] = None, + ) -> None: + self.observables = tuple(observables or cosmo_data.observables) + self._preloaded_cells = ( + None if data_cells is None else {k: np.array(v) for k, v in data_cells.items()} + ) + self.photo_cov_data: Optional[pcov.PhotoCov] = None + self._ells = None + self._ellmax_WL = None + self._ellmax_GC = None + self._ellmax_XC = None + super().__init__(cosmo_data=cosmo_data, cosmo_theory=cosmo_theory, leg_flag="cells") + + def compute_data(self) -> Dict[str, np.ndarray]: + if self._preloaded_cells is not None: + self._ells = np.array(self._preloaded_cells.get("ells"), copy=True) + return self._preloaded_cells + + photo_cls = getattr(self.cosmo_data, "photo_obs_fid", None) + if photo_cls is None: + photo_cls = pobs.ComputeCls( + cosmopars=self.cosmo_data.fiducialcosmopars, + photopars=self.cosmo_data.photopars, + IApars=self.cosmo_data.IApars, + biaspars=self.cosmo_data.photobiaspars, + fiducial_cosmo=self.cosmo_data.fiducialcosmo, + ) + + self.photo_cov_data = getattr(self.cosmo_data, "photo_LSS", None) + if self.photo_cov_data is None: + self.photo_cov_data = pcov.PhotoCov( + cosmopars=self.cosmo_data.fiducialcosmopars, + photopars=self.cosmo_data.photopars, + IApars=self.cosmo_data.IApars, + biaspars=self.cosmo_data.photobiaspars, + fiducial_Cls=photo_cls, + ) + + cells = _cells_from_cls(photo_cls, self.photo_cov_data, self.observables) + self._ells = cells["ells"] + self._ellmax_WL = self.cosmo_data.specs.get("lmax_WL") + self._ellmax_GC = self.cosmo_data.specs.get("lmax_GCph") + if self._ellmax_WL is not None and self._ellmax_GC is not None: + self._ellmax_XC = min(self._ellmax_WL, self._ellmax_GC) + else: + self._ellmax_XC = None + return cells + + def compute_theory(self, param_dict: Dict[str, Any]) -> Dict[str, np.ndarray]: + params = deepcopy(param_dict) + + cosmopars = _dict_with_updates(self.cosmo_theory.fiducialcosmopars, params) + photopars = _dict_with_updates(self.cosmo_theory.photopars, params) + IApars = _dict_with_updates(self.cosmo_theory.IApars, params) + photobias = _dict_with_updates(self.cosmo_theory.photobiaspars, params) + + if params: + logger.debug( + "PhotometricLikelihood received unused parameters: %s", + ", ".join(sorted(params.keys())), + ) + + photo_cls = pobs.ComputeCls( + cosmopars=cosmopars, + photopars=photopars, + IApars=IApars, + biaspars=photobias, + fiducial_cosmo=None, + ) + return _cells_from_cls(photo_cls, self.photo_cov_data, self.observables) + + def compute_chi2(self, theory_obs: Dict[str, np.ndarray]) -> float: + if self._ells is None or self.photo_cov_data is None: + raise RuntimeError("Data ells not initialised") + + ells = self._ells + chi2 = 0.0 + fsky_wl = getattr(self.photo_cov_data, "fsky_WL", 1.0) + fsky_gc = getattr(self.photo_cov_data, "fsky_GCph", 1.0) + + def ell_subset(limit: Optional[int]) -> tuple[int, np.ndarray, np.ndarray]: + working = np.array(ells, copy=True) + if limit is not None: + insert_idx = np.searchsorted(working, limit) + working = np.insert(working, insert_idx, limit) + else: + insert_idx = len(working) + + deltas = np.diff(working) + if deltas.size < insert_idx: + last = deltas[-1] if deltas.size else 1.0 + deltas = np.concatenate([deltas, [last]]) + else: + deltas = deltas[:insert_idx] + + return insert_idx, working[:insert_idx], deltas + + if "Cell_LL" in self.data_obs and "Cell_LL" in theory_obs: + n_wl, ell_wl, d_ell_wl = ell_subset(self._ellmax_WL) + chi2 += fsky_wl * _chi2_per_obs( + self.data_obs["Cell_LL"][:n_wl], theory_obs["Cell_LL"][:n_wl], ell_wl, d_ell_wl + ) + + if "Cell_GG" in self.data_obs and "Cell_GG" in theory_obs: + n_gc, ell_gc, d_ell_gc = ell_subset(self._ellmax_GC) + chi2 += fsky_gc * _chi2_per_obs( + self.data_obs["Cell_GG"][:n_gc], theory_obs["Cell_GG"][:n_gc], ell_gc, d_ell_gc + ) + + if ( + "Cell_GL" in self.data_obs + and "Cell_GL" in theory_obs + and "Cell_GG" in theory_obs + and "Cell_LL" in theory_obs + ): + n_xc, ell_xc, d_ell_xc = ell_subset(self._ellmax_XC) + big_th = np.block( + [ + [ + theory_obs["Cell_LL"][:n_xc], + np.transpose(theory_obs["Cell_GL"], (0, 2, 1))[:n_xc], + ], + [ + theory_obs["Cell_GL"][:n_xc], + theory_obs["Cell_GG"][:n_xc], + ], + ] + ) + big_fid = np.block( + [ + [ + self.data_obs["Cell_LL"][:n_xc], + np.transpose(self.data_obs["Cell_GL"], (0, 2, 1))[:n_xc], + ], + [ + self.data_obs["Cell_GL"][:n_xc], + self.data_obs["Cell_GG"][:n_xc], + ], + ] + ) + chi2 += np.sqrt( + self.photo_cov_data.fsky_WL * self.photo_cov_data.fsky_GCph + ) * _chi2_per_obs(big_fid, big_th, ell_xc, d_ell_xc) + chi2 += fsky_wl * _chi2_per_obs( + self.data_obs["Cell_LL"][:n_xc], theory_obs["Cell_LL"][:n_xc], ell_xc, d_ell_xc + ) + + return float(chi2) diff --git a/cosmicfishpie/likelihood/spectro_like.py b/cosmicfishpie/likelihood/spectro_like.py index b5f070e..e1e9e74 100644 --- a/cosmicfishpie/likelihood/spectro_like.py +++ b/cosmicfishpie/likelihood/spectro_like.py @@ -1,6 +1,8 @@ +from __future__ import annotations + import logging -from collections.abc import Sequence from copy import deepcopy +from typing import Any, Dict, Iterable, Optional import numpy as np from scipy.integrate import simpson @@ -13,15 +15,13 @@ from cosmicfishpie.utilities import legendre_tools as lgt from cosmicfishpie.utilities.utils import printing as upr +from .base import Likelihood, NautilusMixin + logger = logging.getLogger("cosmicfishpie.cosmology.nuisance") logger.setLevel(logging.INFO) upr.debug = False -def is_indexable_iterable(var): - return isinstance(var, (list, np.ndarray, Sequence)) and not isinstance(var, (str, bytes)) - - def observable_Pgg(theory_spectro, cosmoFM: cosmicfish.FisherMatrix, nuisance_shot=None): """ cosmoFM: cosmicfish.FisherMatrix object @@ -56,6 +56,16 @@ def legendre_Pgg(obs_Pgg, cosmoFM: cosmicfish.FisherMatrix): return P_ell +def _dict_with_updates(template: Dict[str, Any], pool: Dict[str, Any]) -> Dict[str, Any]: + """Return a deepcopy of ``template`` with any matching keys updated from ``pool``.""" + + updated = deepcopy(template) + for key in template: + if key in pool: + updated[key] = pool.pop(key) + return updated + + def compute_covariance_legendre(P_ell, cosmoFM: cosmicfish.FisherMatrix): """ Compute covariance matrix for power spectrum multipoles @@ -179,29 +189,31 @@ def compute_wedge_chi2(P_obs_data, P_obs_theory, cosmoFM_data: cosmicfish.Fisher return chi2 -def compute_theory_spectro(param_dict, cosmoFM_theory: cosmicfish.FisherMatrix, leg_flag="wedges"): +def compute_theory_spectro( + param_dict: Dict[str, Any], + cosmoFM_theory: cosmicfish.FisherMatrix, + leg_flag: str = "wedges", +) -> np.ndarray: + params = deepcopy(param_dict) z_bins = cosmoFM_theory.pk_cov.global_z_bin_mids nuisance_shot = np.zeros(len(z_bins)) - pshotpars = deepcopy(cosmoFM_theory.PShotpars) - for ii, pp in enumerate(pshotpars.keys()): - nuisance_shot[ii] = param_dict.pop(pp, cosmoFM_theory.PShotpars[pp]) - - spectrobiaspars = deepcopy(cosmoFM_theory.Spectrobiaspars) - for ii, pp in enumerate(spectrobiaspars.keys()): - spectrobiaspars[pp] = param_dict.pop(pp, cosmoFM_theory.Spectrobiaspars[pp]) + for index, key in enumerate(cosmoFM_theory.PShotpars.keys()): + nuisance_shot[index] = params.pop(key, cosmoFM_theory.PShotpars[key]) - spectrononlinearpars = deepcopy(cosmoFM_theory.Spectrononlinpars) - for ii, pp in enumerate(spectrononlinearpars.keys()): - spectrononlinearpars[pp] = param_dict.pop(pp, cosmoFM_theory.Spectrononlinpars[pp]) + spectrobiaspars = _dict_with_updates(cosmoFM_theory.Spectrobiaspars, params) + spectrononlinearpars = _dict_with_updates(cosmoFM_theory.Spectrononlinpars, params) + IMbiaspars = _dict_with_updates(cosmoFM_theory.IMbiasparams, params) + cosmological = _dict_with_updates(cosmoFM_theory.fiducialcosmopars, params) - IMbiaspars = deepcopy(cosmoFM_theory.IMbiasparams) - for pp in IMbiaspars.keys(): - IMbiaspars[pp] = param_dict.pop(pp, cosmoFM_theory.IMbiasparams[pp]) + if params: + logger.debug( + "SpectroLikelihood received unused parameters: %s", ", ".join(sorted(params.keys())) + ) spectro_vary = spobs.ComputeGalSpectro( - param_dict, + cosmological, spectrobiaspars=spectrobiaspars, spectrononlinearpars=spectrononlinearpars, PShotpars=cosmoFM_theory.PShotpars, @@ -211,14 +223,80 @@ def compute_theory_spectro(param_dict, cosmoFM_theory: cosmicfish.FisherMatrix, configuration=cosmoFM_theory, ) spectro_cov_vary = spcov.SpectroCov( - fiducialpars=param_dict, configuration=cosmoFM_theory, fiducial_specobs=spectro_vary + fiducialpars=cosmological, configuration=cosmoFM_theory, fiducial_specobs=spectro_vary ) obsPgg_vary = observable_Pgg(spectro_cov_vary, cosmoFM_theory, nuisance_shot=nuisance_shot) if leg_flag == "wedges": return obsPgg_vary - elif leg_flag == "legendre": - P_ell_vary = legendre_Pgg(obsPgg_vary, cosmoFM_theory) - return P_ell_vary + if leg_flag == "legendre": + return legendre_Pgg(obsPgg_vary, cosmoFM_theory) + raise ValueError(f"Unknown leg_flag '{leg_flag}'. Use 'wedges' or 'legendre'.") + + +class SpectroLikelihood(Likelihood, NautilusMixin): + """Likelihood for spectroscopic clustering using CosmicFish Fisher matrices.""" + + def __init__( + self, + *, + cosmoFM_data: cosmicfish.FisherMatrix, + cosmoFM_theory: Optional[cosmicfish.FisherMatrix] = None, + leg_flag: str = "wedges", + data_obs: Optional[np.ndarray] = None, + nuisance_shot: Optional[Iterable[float]] = None, + ) -> None: + self._preloaded_data = None if data_obs is None else np.array(data_obs) + self._nuisance_shot = ( + None if nuisance_shot is None else np.array(nuisance_shot, dtype=float) + ) + self._inv_cov_legendre = None + self._data_wedges = None + super().__init__(cosmo_data=cosmoFM_data, cosmo_theory=cosmoFM_theory, leg_flag=leg_flag) + + @property + def data_wedges(self) -> Optional[np.ndarray]: + """Return wedge data even when the likelihood operates in Legendre space.""" + + return self._data_wedges + + def compute_data(self) -> np.ndarray: + if self._preloaded_data is not None: + data = np.array(self._preloaded_data, copy=False) + if self.leg_flag == "legendre": + _, self._inv_cov_legendre = compute_covariance_legendre(data, self.cosmoFM_data) + else: + self._data_wedges = data + return data + + if not hasattr(self.cosmoFM_data, "pk_cov") or self.cosmoFM_data.pk_cov is None: + raise AttributeError( + "cosmoFM_data.pk_cov is not available. Ensure the FisherMatrix was initialised for spectroscopic probes." + ) + + obsPgg = observable_Pgg( + self.cosmoFM_data.pk_cov, + self.cosmoFM_data, + nuisance_shot=self._nuisance_shot, + ) + self._data_wedges = obsPgg + if self.leg_flag == "legendre": + p_ell = legendre_Pgg(obsPgg, self.cosmoFM_data) + _, self._inv_cov_legendre = compute_covariance_legendre(p_ell, self.cosmoFM_data) + return p_ell + return obsPgg + + def compute_theory(self, param_dict: Dict[str, Any]) -> np.ndarray: + return compute_theory_spectro(param_dict, self.cosmoFM_theory, self.leg_flag) + + def compute_chi2(self, theory_obs: np.ndarray) -> float: + if self.leg_flag == "wedges": + return compute_wedge_chi2(self.data_obs, theory_obs, self.cosmoFM_data) + + if self._inv_cov_legendre is None: + _, self._inv_cov_legendre = compute_covariance_legendre( + self.data_obs, self.cosmoFM_data + ) + return compute_chi2_legendre(self.data_obs, theory_obs, self._inv_cov_legendre) def loglike( @@ -230,26 +308,25 @@ def loglike( cosmoFM_data: cosmicfish.FisherMatrix = None, data_obsPgg: np.ndarray = None, ): - if theory_obsPgg is None and param_vec is not None: + if cosmoFM_data is None: + return -np.inf + + likelihood = SpectroLikelihood( + cosmoFM_data=cosmoFM_data, + cosmoFM_theory=cosmoFM_theory, + leg_flag=leg_flag, + data_obs=data_obsPgg, + ) + + if theory_obsPgg is not None: + return -0.5 * likelihood.compute_chi2(theory_obsPgg) + + try: if isinstance(param_vec, dict): - param_dict = deepcopy(param_vec) - elif is_indexable_iterable(param_vec) and prior is not None: - # print(f'Loading prior with keys: {prior.keys}') - param_dict = {key: param_vec[i] for i, key in enumerate(prior.keys)} - theory_obsPgg = compute_theory_spectro(param_dict, cosmoFM_theory, leg_flag) - elif theory_obsPgg is not None: - pass - else: - upr.debug_print("No theory_obsPgg provided and no param_vec provided") - return 1e23 - if leg_flag == "wedges": - chi2 = compute_wedge_chi2( - P_obs_data=data_obsPgg, P_obs_theory=theory_obsPgg, cosmoFM_data=cosmoFM_data - ) - elif leg_flag == "legendre": - P_ell_data = legendre_Pgg(data_obsPgg, cosmoFM_data) - covariance_leg, inv_covariance_leg = compute_covariance_legendre( - P_ell=P_ell_data, cosmoFM=cosmoFM_data - ) - chi2 = compute_chi2_legendre(P_ell_data, theory_obsPgg, inv_covariance_leg) - return -0.5 * chi2 + params = dict(param_vec) + else: + params = likelihood.build_param_dict(param_vec=param_vec, prior=prior) + except (TypeError, ValueError, AttributeError): + return -np.inf + + return likelihood.loglike(param_dict=params) diff --git a/cosmicfishpie/version.py b/cosmicfishpie/version.py index b1cfdab..703a8e9 100644 --- a/cosmicfishpie/version.py +++ b/cosmicfishpie/version.py @@ -2,7 +2,7 @@ _MINOR = "2" # On main and in a nightly release the patch should be one ahead of the last # released build. -_PATCH = "4" +_PATCH = "5" # This is mainly for nightly builds which have the suffix ".dev$DATE". See # https://semver.org/#is-v123-a-semantic-version for the semantics. _SUFFIX = "" diff --git a/coverage-badge.svg b/coverage-badge.svg deleted file mode 100644 index f11b9fb..0000000 --- a/coverage-badge.svg +++ /dev/null @@ -1 +0,0 @@ -coverage: 23.96%coverage23.96% \ No newline at end of file diff --git a/coverage.svg b/coverage.svg deleted file mode 100644 index f78bcfc..0000000 --- a/coverage.svg +++ /dev/null @@ -1,21 +0,0 @@ - - - - - - - - - - - - - - - - coverage - coverage - 35% - 35% - - diff --git a/notebooks/likelihod_photometric.ipynb b/notebooks/likelihod_photometric.ipynb index b68db0a..2a5c9d8 100644 --- a/notebooks/likelihod_photometric.ipynb +++ b/notebooks/likelihod_photometric.ipynb @@ -71,7 +71,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -135,7 +135,7 @@ "\n", " -> Computing cosmology at the fiducial point\n", "\n", - " ---> Cosmological functions obtained in: 0.16 s\n" + " ---> Cosmological functions obtained in: 0.15 s\n" ] } ], @@ -222,16 +222,16 @@ " 'betaIA': 2.17,\n", " 'etaIA': -0.41,\n", " 'bias_model': 'binned',\n", - " 'b1': 1.0997727037892875,\n", - " 'b2': 1.220245876862528,\n", - " 'b3': 1.2723993083933989,\n", - " 'b4': 1.316624471897739,\n", - " 'b5': 1.35812370570578,\n", - " 'b6': 1.3998214171814918,\n", - " 'b7': 1.4446452851824907,\n", - " 'b8': 1.4964959071110084,\n", - " 'b9': 1.5652475842498528,\n", - " 'b10': 1.7429859437184225,\n", + " 'b1': np.float64(1.0997727037892875),\n", + " 'b2': np.float64(1.220245876862528),\n", + " 'b3': np.float64(1.2723993083933989),\n", + " 'b4': np.float64(1.316624471897739),\n", + " 'b5': np.float64(1.35812370570578),\n", + " 'b6': np.float64(1.3998214171814918),\n", + " 'b7': np.float64(1.4446452851824907),\n", + " 'b8': np.float64(1.4964959071110084),\n", + " 'b9': np.float64(1.5652475842498528),\n", + " 'b10': np.float64(1.7429859437184225),\n", " 'fout': 0.1,\n", " 'co': 1,\n", " 'cb': 1,\n", @@ -261,12 +261,12 @@ "text": [ "range(1, 11)\n", "range(1, 11)\n", - "[35454308.58012683 35454308.58012683 35454308.58012683 35454308.58012683\n", - " 35454308.58012683 35454308.58012683 35454308.58012683 35454308.58012683\n", - " 35454308.58012683 35454308.58012683]\n", - "[35454308.58012683 35454308.58012683 35454308.58012683 35454308.58012683\n", - " 35454308.58012683 35454308.58012683 35454308.58012683 35454308.58012683\n", - " 35454308.58012683 35454308.58012683]\n" + "[35454308.58012684 35454308.58012684 35454308.58012684 35454308.58012684\n", + " 35454308.58012684 35454308.58012684 35454308.58012684 35454308.58012684\n", + " 35454308.58012684 35454308.58012684]\n", + "[35454308.58012684 35454308.58012684 35454308.58012684 35454308.58012684\n", + " 35454308.58012684 35454308.58012684 35454308.58012684 35454308.58012684\n", + " 35454308.58012684 35454308.58012684]\n" ] } ], @@ -584,16 +584,16 @@ " 'mnu': 0.06,\n", " 'Neff': 3.044,\n", " 'bias_model': 'binned',\n", - " 'b1': 1.0997727037892875,\n", - " 'b2': 1.220245876862528,\n", - " 'b3': 1.2723993083933989,\n", - " 'b4': 1.316624471897739,\n", - " 'b5': 1.35812370570578,\n", - " 'b6': 1.3998214171814918,\n", - " 'b7': 1.4446452851824907,\n", - " 'b8': 1.4964959071110084,\n", - " 'b9': 1.5652475842498528,\n", - " 'b10': 1.7429859437184225,\n", + " 'b1': np.float64(1.0997727037892875),\n", + " 'b2': np.float64(1.220245876862528),\n", + " 'b3': np.float64(1.2723993083933989),\n", + " 'b4': np.float64(1.316624471897739),\n", + " 'b5': np.float64(1.35812370570578),\n", + " 'b6': np.float64(1.3998214171814918),\n", + " 'b7': np.float64(1.4446452851824907),\n", + " 'b8': np.float64(1.4964959071110084),\n", + " 'b9': np.float64(1.5652475842498528),\n", + " 'b10': np.float64(1.7429859437184225),\n", " 'fout': 0.1,\n", " 'co': 1,\n", " 'cb': 1,\n", @@ -625,42 +625,22 @@ "name": "stdout", "output_type": "stream", "text": [ - "Sample likelihood 4.038266295364123e-11\n" + "Sample likelihood 4.330972755260946e-11\n", + "Sample likelihood -1568.703678015382\n", + "Sample likelihood -1522.4374216102938\n" ] } ], "source": [ - "samp1dic = {'Omegam': 0.3145714273,\n", - " 'Omegab': 0.0491989,\n", - " 'h': 0.6737,\n", - " 'ns': 0.96605,\n", - " 'sigma8': 0.81,\n", - " 'w0': -1.0,\n", - " 'wa': 0.0,\n", - " 'mnu': 0.06,\n", - " 'Neff': 3.044,\n", - " 'bias_model': 'binned',\n", - " 'b1': 1.0997727037892875,\n", - " 'b2': 1.220245876862528,\n", - " 'b3': 1.2723993083933989,\n", - " 'b4': 1.316624471897739,\n", - " 'b5': 1.35812370570578,\n", - " 'b6': 1.3998214171814918,\n", - " 'b7': 1.4446452851824907,\n", - " 'b8': 1.4964959071110084,\n", - " 'b9': 1.5652475842498528,\n", - " 'b10': 1.7429859437184225,\n", - " 'fout': 0.1,\n", - " 'co': 1,\n", - " 'cb': 1,\n", - " 'sigma_o': 0.05,\n", - " 'sigma_b': 0.05,\n", - " 'zo': 0.1,\n", - " 'zb': 0.0,\n", - " 'IA_model': 'eNLA',\n", - " 'AIA': 1.72,\n", - " 'betaIA': 2.17,\n", - " 'etaIA': -0.41*1.1}\n", + "samp1dic = cosmoFM_fid.allparams.copy()\n", + "print(\"Sample likelihood\", loglike(samp1dic))\n", + "\n", + "samp1dic['b1'] = 1.0\n", + "\n", + "print(\"Sample likelihood\", loglike(samp1dic))\n", + "\n", + "samp1dic['h'] = 0.7\n", + "\n", "print(\"Sample likelihood\", loglike(samp1dic))" ] }, @@ -672,7 +652,7 @@ { "data": { "text/plain": [ - "4.038266295364123e-11" + "np.float64(4.330972755260946e-11)" ] }, "execution_count": 25, @@ -706,16 +686,16 @@ " 'betaIA': 2.17,\n", " 'etaIA': -0.41,\n", " 'bias_model': 'binned',\n", - " 'b1': 1.0997727037892875,\n", - " 'b2': 1.220245876862528,\n", - " 'b3': 1.2723993083933989,\n", - " 'b4': 1.316624471897739,\n", - " 'b5': 1.35812370570578,\n", - " 'b6': 1.3998214171814918,\n", - " 'b7': 1.4446452851824907,\n", - " 'b8': 1.4964959071110084,\n", - " 'b9': 1.5652475842498528,\n", - " 'b10': 1.7429859437184225,\n", + " 'b1': np.float64(1.0997727037892875),\n", + " 'b2': np.float64(1.220245876862528),\n", + " 'b3': np.float64(1.2723993083933989),\n", + " 'b4': np.float64(1.316624471897739),\n", + " 'b5': np.float64(1.35812370570578),\n", + " 'b6': np.float64(1.3998214171814918),\n", + " 'b7': np.float64(1.4446452851824907),\n", + " 'b8': np.float64(1.4964959071110084),\n", + " 'b9': np.float64(1.5652475842498528),\n", + " 'b10': np.float64(1.7429859437184225),\n", " 'fout': 0.1,\n", " 'co': 1,\n", " 'cb': 1,\n", @@ -734,112 +714,6 @@ "photo_cov_fid.allparsfid" ] }, - { - "cell_type": "code", - "execution_count": 27, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "In class: FisherMatrix ----> Computing photo Fisher matrix\n", - "\n", - "Computing fiducial\n", - "\n", - "Fiducial generated in 0.89 s\n", - "\n", - "Noise added to fiducial\n", - "\n", - "Noisy Cls generated in 0.00 s\n", - "\n", - "Computed covariance matrix\n", - "\n", - "Covmat of Cls generated in 0.20 s\n", - "\n", - "Total calculation in 1.08 s\n", - ">> computing derivs >>\n", - "\n", - " +++ Computing derivative on Omegam\n", - "\n", - "In class: derivatives Derivative on Omegam done! in : 1.99 s\n", - "\n", - " +++ Computing derivative on Omegab\n", - "\n", - "In class: derivatives Derivative on Omegab done! in : 1.95 s\n", - "\n", - " +++ Computing derivative on h\n", - "\n", - "In class: derivatives Derivative on h done! in : 1.97 s\n", - "\n", - " +++ Computing derivative on ns\n", - "\n", - "In class: derivatives Derivative on ns done! in : 2.02 s\n", - "\n", - " +++ Computing derivative on sigma8\n", - "\n", - "In class: derivatives Derivative on sigma8 done! in : 1.96 s\n", - "\n", - " +++ Computing derivative on b1\n", - "\n", - "In class: derivatives Derivative on b1 done! in : 2.46 s\n", - "\n", - " +++ Computing derivative on b2\n", - "\n", - "In class: derivatives Derivative on b2 done! in : 2.05 s\n", - "\n", - " +++ Computing derivative on b3\n", - "\n", - "In class: derivatives Derivative on b3 done! in : 2.01 s\n", - "\n", - " +++ Computing derivative on b4\n", - "\n", - "In class: derivatives Derivative on b4 done! in : 1.86 s\n", - "\n", - " +++ Computing derivative on b5\n", - "\n", - "In class: derivatives Derivative on b5 done! in : 1.80 s\n", - "\n", - " +++ Computing derivative on b6\n", - "\n", - "In class: derivatives Derivative on b6 done! in : 1.79 s\n", - "\n", - " +++ Computing derivative on b7\n", - "\n", - "In class: derivatives Derivative on b7 done! in : 1.83 s\n", - "\n", - " +++ Computing derivative on b8\n", - "\n", - "In class: derivatives Derivative on b8 done! in : 1.81 s\n", - "\n", - " +++ Computing derivative on b9\n", - "\n", - "In class: derivatives Derivative on b9 done! in : 1.73 s\n", - "\n", - " +++ Computing derivative on b10\n", - "\n", - "In class: derivatives Derivative on b10 done! in : 1.81 s\n", - "\n", - "In class: PhotoCov -->> Derivatives computed in 29.05 s\n", - "\n", - "In class: FisherMatrix Computing Fisher matrix\n", - "\n", - " Finished calculation of Fisher Matrix for ['GCph'] in: 0.02 s\n", - "\n", - " Fisher matrix calculation finished in 30.50 s\n", - "\n", - " Fisher matrix exported: ./results/CosmicFish_v1.2.4_cosmicjellyfish_Euclid-ISTF-Pess-GGphoto_symb_withnuis_GCph_fishermatrix.txt\n", - "\n", - "\n", - "In class: FisherMatrix CosmicFish settings and Survey specifications exported: ./results/CosmicFish_v1.2.4_cosmicjellyfish_Euclid-ISTF-Pess-GGphoto_symb_withnuis_GCph_fishermatrix_specifications.dat\n" - ] - } - ], - "source": [ - "fishmat_photo = cosmoFM_fid.compute()" - ] - }, { "cell_type": "code", "execution_count": 28, @@ -980,7 +854,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Starting sampler at 2025-09-25 22:18:55\n" + "Starting sampler at 2025-09-29 12:16:05\n" ] } ], @@ -991,27 +865,9 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "ename": "OSError", - "evalue": "Unable to synchronously open file (truncated file: eof = 24309760, sblock->base_addr = 0, stored_eof = 24340472)", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mOSError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[35], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m sampler \u001b[38;5;241m=\u001b[39m \u001b[43mSampler\u001b[49m\u001b[43m(\u001b[49m\u001b[43mprior_chosen\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 2\u001b[0m \u001b[43m \u001b[49m\u001b[43mloglike\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 3\u001b[0m \u001b[43m \u001b[49m\u001b[43mn_live\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msampler_settings\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mn_live\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 4\u001b[0m \u001b[43m \u001b[49m\u001b[43mn_networks\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msampler_settings\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mn_networks\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 5\u001b[0m \u001b[43m \u001b[49m\u001b[43mn_batch\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msampler_settings\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mn_batch\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 6\u001b[0m \u001b[43m \u001b[49m\u001b[43mpool\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msampler_settings\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mpool\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 7\u001b[0m \u001b[43m \u001b[49m\u001b[43mpass_dict\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m 8\u001b[0m \u001b[43m \u001b[49m\u001b[43mfilepath\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43moutroot\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43m.hdf5\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 9\u001b[0m \u001b[43m \u001b[49m\u001b[43mresume\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m 10\u001b[0m \u001b[43m \u001b[49m\u001b[43mlikelihood_kwargs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m{\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mprior\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m:\u001b[49m\u001b[43m \u001b[49m\u001b[43mprior_chosen\u001b[49m\u001b[43m}\u001b[49m\n\u001b[1;32m 11\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 12\u001b[0m sampler\u001b[38;5;241m.\u001b[39mrun(verbose\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m, discard_exploration\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[1;32m 13\u001b[0m log_z_all \u001b[38;5;241m=\u001b[39m sampler\u001b[38;5;241m.\u001b[39mevidence()\n", - "File \u001b[0;32m~/anaconda3/envs/cosmicfishpie/lib/python3.10/site-packages/nautilus/sampler.py:328\u001b[0m, in \u001b[0;36mSampler.__init__\u001b[0;34m(self, prior, likelihood, n_dim, n_live, n_update, enlarge_per_dim, n_points_min, split_threshold, n_networks, neural_network_kwargs, prior_args, prior_kwargs, likelihood_args, likelihood_kwargs, n_batch, n_like_new_bound, vectorized, pass_dict, pool, seed, blobs_dtype, filepath, resume)\u001b[0m\n\u001b[1;32m 326\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfilepath \u001b[38;5;241m=\u001b[39m filepath\n\u001b[1;32m 327\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m resume \u001b[38;5;129;01mand\u001b[39;00m filepath \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;129;01mand\u001b[39;00m Path(filepath)\u001b[38;5;241m.\u001b[39mexists():\n\u001b[0;32m--> 328\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[43mh5py\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mFile\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfilepath\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mr\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mas\u001b[39;00m fstream:\n\u001b[1;32m 330\u001b[0m group \u001b[38;5;241m=\u001b[39m fstream[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msampler\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[1;32m 332\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mrng\u001b[38;5;241m.\u001b[39mbit_generator\u001b[38;5;241m.\u001b[39mstate \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mdict\u001b[39m(\n\u001b[1;32m 333\u001b[0m bit_generator\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mPCG64\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m 334\u001b[0m state\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mdict\u001b[39m(\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 337\u001b[0m has_uint32\u001b[38;5;241m=\u001b[39mgroup\u001b[38;5;241m.\u001b[39mattrs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mrng_has_uint32\u001b[39m\u001b[38;5;124m'\u001b[39m],\n\u001b[1;32m 338\u001b[0m uinteger\u001b[38;5;241m=\u001b[39mgroup\u001b[38;5;241m.\u001b[39mattrs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mrng_uinteger\u001b[39m\u001b[38;5;124m'\u001b[39m])\n", - "File \u001b[0;32m~/anaconda3/envs/cosmicfishpie/lib/python3.10/site-packages/h5py/_hl/files.py:561\u001b[0m, in \u001b[0;36mFile.__init__\u001b[0;34m(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, fs_strategy, fs_persist, fs_threshold, fs_page_size, page_buf_size, min_meta_keep, min_raw_keep, locking, alignment_threshold, alignment_interval, meta_block_size, **kwds)\u001b[0m\n\u001b[1;32m 552\u001b[0m fapl \u001b[38;5;241m=\u001b[39m make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0,\n\u001b[1;32m 553\u001b[0m locking, page_buf_size, min_meta_keep, min_raw_keep,\n\u001b[1;32m 554\u001b[0m alignment_threshold\u001b[38;5;241m=\u001b[39malignment_threshold,\n\u001b[1;32m 555\u001b[0m alignment_interval\u001b[38;5;241m=\u001b[39malignment_interval,\n\u001b[1;32m 556\u001b[0m meta_block_size\u001b[38;5;241m=\u001b[39mmeta_block_size,\n\u001b[1;32m 557\u001b[0m \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwds)\n\u001b[1;32m 558\u001b[0m fcpl \u001b[38;5;241m=\u001b[39m make_fcpl(track_order\u001b[38;5;241m=\u001b[39mtrack_order, fs_strategy\u001b[38;5;241m=\u001b[39mfs_strategy,\n\u001b[1;32m 559\u001b[0m fs_persist\u001b[38;5;241m=\u001b[39mfs_persist, fs_threshold\u001b[38;5;241m=\u001b[39mfs_threshold,\n\u001b[1;32m 560\u001b[0m fs_page_size\u001b[38;5;241m=\u001b[39mfs_page_size)\n\u001b[0;32m--> 561\u001b[0m fid \u001b[38;5;241m=\u001b[39m \u001b[43mmake_fid\u001b[49m\u001b[43m(\u001b[49m\u001b[43mname\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43muserblock_size\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfapl\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfcpl\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mswmr\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mswmr\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 563\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(libver, \u001b[38;5;28mtuple\u001b[39m):\n\u001b[1;32m 564\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_libver \u001b[38;5;241m=\u001b[39m libver\n", - "File \u001b[0;32m~/anaconda3/envs/cosmicfishpie/lib/python3.10/site-packages/h5py/_hl/files.py:235\u001b[0m, in \u001b[0;36mmake_fid\u001b[0;34m(name, mode, userblock_size, fapl, fcpl, swmr)\u001b[0m\n\u001b[1;32m 233\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m swmr \u001b[38;5;129;01mand\u001b[39;00m swmr_support:\n\u001b[1;32m 234\u001b[0m flags \u001b[38;5;241m|\u001b[39m\u001b[38;5;241m=\u001b[39m h5f\u001b[38;5;241m.\u001b[39mACC_SWMR_READ\n\u001b[0;32m--> 235\u001b[0m fid \u001b[38;5;241m=\u001b[39m \u001b[43mh5f\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen\u001b[49m\u001b[43m(\u001b[49m\u001b[43mname\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mflags\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfapl\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfapl\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 236\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m mode \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mr+\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 237\u001b[0m fid \u001b[38;5;241m=\u001b[39m h5f\u001b[38;5;241m.\u001b[39mopen(name, h5f\u001b[38;5;241m.\u001b[39mACC_RDWR, fapl\u001b[38;5;241m=\u001b[39mfapl)\n", - "File \u001b[0;32mh5py/_objects.pyx:54\u001b[0m, in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n", - "File \u001b[0;32mh5py/_objects.pyx:55\u001b[0m, in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n", - "File \u001b[0;32mh5py/h5f.pyx:102\u001b[0m, in \u001b[0;36mh5py.h5f.open\u001b[0;34m()\u001b[0m\n", - "\u001b[0;31mOSError\u001b[0m: Unable to synchronously open file (truncated file: eof = 24309760, sblock->base_addr = 0, stored_eof = 24340472)" - ] - } - ], + "outputs": [], "source": [ "sampler = Sampler(prior_chosen, \n", " loglike, \n", @@ -1034,17 +890,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "ename": "", - "evalue": "", - "output_type": "error", - "traceback": [ - "\u001b[1;31mnotebook controller is DISPOSED. \n", - "\u001b[1;31mView Jupyter log for further details." - ] - } - ], + "outputs": [], "source": [ "tfin = time.time()\n", "elapsed = tfin - tini\n", @@ -1059,17 +905,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "ename": "", - "evalue": "", - "output_type": "error", - "traceback": [ - "\u001b[1;31mnotebook controller is DISPOSED. \n", - "\u001b[1;31mView Jupyter log for further details." - ] - } - ], + "outputs": [], "source": [ "sample_wghlkl = (np.vstack((points_all.T, np.exp(log_w_all), log_l_all)).T)" ] @@ -1078,17 +914,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "ename": "", - "evalue": "", - "output_type": "error", - "traceback": [ - "\u001b[1;31mnotebook controller is DISPOSED. \n", - "\u001b[1;31mView Jupyter log for further details." - ] - } - ], + "outputs": [], "source": [ "outfile_chain = options[\"outroot\"]+\".txt\"\n", "print(f\"Saving chain to text file {outfile_chain}\")" @@ -1098,17 +924,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "ename": "", - "evalue": "", - "output_type": "error", - "traceback": [ - "\u001b[1;31mnotebook controller is DISPOSED. \n", - "\u001b[1;31mView Jupyter log for further details." - ] - } - ], + "outputs": [], "source": [ "headerlist = ['loglike', 'weights'] + list(prior_chosen.keys)\n", "header = \" \".join(headerlist)\n", @@ -1119,17 +935,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "ename": "", - "evalue": "", - "output_type": "error", - "traceback": [ - "\u001b[1;31mnotebook controller is DISPOSED. \n", - "\u001b[1;31mView Jupyter log for further details." - ] - } - ], + "outputs": [], "source": [ "np.savetxt(outfile_chain, sample_wghlkl, header=header)" ] diff --git a/notebooks/likelihod_photometric_module.ipynb b/notebooks/likelihod_photometric_module.ipynb new file mode 100644 index 0000000..90e4a29 --- /dev/null +++ b/notebooks/likelihod_photometric_module.ipynb @@ -0,0 +1,1147 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import seaborn as sns\n", + "import numpy as np\n", + "import logging\n", + "from itertools import product\n", + "from copy import deepcopy, copy\n", + "from collections.abc import Sequence" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "from cosmicfishpie.fishermatrix import cosmicfish\n", + "from cosmicfishpie.LSSsurvey import photo_obs as pobs\n", + "from cosmicfishpie.LSSsurvey import photo_cov as pcov\n", + "from cosmicfishpie.likelihood import PhotometricLikelihood\n", + "from cosmicfishpie.utilities.utils import printing as upr\n", + "from nautilus import Prior\n", + "from nautilus import Sampler\n", + "import re\n", + "import time" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "snscolors = sns.color_palette(\"colorblind\")\n", + "def is_indexable_iterable(var):\n", + " return isinstance(var, (list, np.ndarray, Sequence)) and not isinstance(var, (str, bytes))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "logger = logging.getLogger(\"cosmicfishpie.cosmology.nuisance\")\n", + "logger.setLevel(logging.INFO)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "upr.debug = False" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "upr.debug_print(\"test\")" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "outroot = \"cosmicjellyfish_Euclid-ISTF-Pess-GCphoto_symb_simple\"" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "sampler_settings = {\n", + " \"n_live\": 1000,\n", + " \"n_networks\": 8,\n", + " \"n_batch\": 256,\n", + " \"pool\": 8,\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "fiducial = {\n", + " \"Omegam\": 0.3145714273,\n", + " \"Omegab\": 0.0491989,\n", + " \"h\": 0.6737,\n", + " \"ns\": 0.96605,\n", + " \"sigma8\": 0.81,\n", + " \"w0\": -1.0,\n", + " \"wa\": 0.0,\n", + " \"mnu\": 0.06,\n", + " \"Neff\": 3.044,\n", + "}\n", + "observables = ['GCph']" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "****************************************************************\n", + " _____ _ _____ __ \n", + " / ___/__ ___ __ _ (_)___/ __(_)__ / / \n", + " / /__/ _ \\(_- Survey loaded: Euclid-Photometric-ISTF-Pessimistic\n", + "\n", + " -> Computing cosmology at the fiducial point\n", + "\n", + " ---> Cosmological functions obtained in: 0.11 s\n" + ] + } + ], + "source": [ + "options = {\n", + " \"accuracy\": 1,\n", + " \"feedback\": 1,\n", + " \"code\": \"symbolic\",\n", + " \"outroot\": outroot,\n", + " \"survey_name\": \"Euclid\",\n", + " \"survey_name_photo\": \"Euclid-Photometric-ISTF-Pessimistic\",\n", + " \"survey_name_spectro\": False,\n", + " \"specs_dir\": \"../cosmicfishpie/configs/default_survey_specifications/\",\n", + " \"cosmo_model\": \"LCDM\",\n", + "}\n", + "cosmoFM_fid = cosmicfish.FisherMatrix(\n", + " fiducialpars=fiducial,\n", + " options=options,\n", + " observables=observables,\n", + " cosmoModel=options[\"cosmo_model\"],\n", + " surveyName=options[\"survey_name\"],\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'Omegam': 0.01,\n", + " 'Omegab': 0.01,\n", + " 'h': 0.01,\n", + " 'ns': 0.01,\n", + " 'sigma8': 0.01,\n", + " 'b1': 0.06,\n", + " 'b2': 0.06,\n", + " 'b3': 0.06,\n", + " 'b4': 0.06,\n", + " 'b5': 0.06,\n", + " 'b6': 0.06,\n", + " 'b7': 0.06,\n", + " 'b8': 0.06,\n", + " 'b9': 0.06,\n", + " 'b10': 0.06}" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cosmoFM_fid.freeparams" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "# #photo_fid = pobs.ComputeCls(cosmopars=cosmoFM_fid.fiducialcosmopars,\n", + "# photopars=cosmoFM_fid.photopars,\n", + "# IApars=cosmoFM_fid.IApars,\n", + "# biaspars=cosmoFM_fid.photobiaspars)\n", + "\n", + "# #photo_fid.compute_all()\n", + "\n", + "# #photo_cov_fid = pcov.PhotoCov(cosmopars=cosmoFM_fid.fiducialcosmopars,\n", + "# #photopars=cosmoFM_fid.photopars,\n", + "# #IApars=cosmoFM_fid.IApars,\n", + "# #biaspars=cosmoFM_fid.photobiaspars,\n", + "# #fiducial_Cls=photo_fid)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'Omegam': 0.3145714273,\n", + " 'Omegab': 0.0491989,\n", + " 'h': 0.6737,\n", + " 'ns': 0.96605,\n", + " 'sigma8': 0.81,\n", + " 'w0': -1.0,\n", + " 'wa': 0.0,\n", + " 'mnu': 0.06,\n", + " 'Neff': 3.044,\n", + " 'bias_model': 'binned',\n", + " 'b1': np.float64(1.0997727037892875),\n", + " 'b2': np.float64(1.220245876862528),\n", + " 'b3': np.float64(1.2723993083933989),\n", + " 'b4': np.float64(1.316624471897739),\n", + " 'b5': np.float64(1.35812370570578),\n", + " 'b6': np.float64(1.3998214171814918),\n", + " 'b7': np.float64(1.4446452851824907),\n", + " 'b8': np.float64(1.4964959071110084),\n", + " 'b9': np.float64(1.5652475842498528),\n", + " 'b10': np.float64(1.7429859437184225),\n", + " 'fout': 0.1,\n", + " 'co': 1,\n", + " 'cb': 1,\n", + " 'sigma_o': 0.05,\n", + " 'sigma_b': 0.05,\n", + " 'zo': 0.1,\n", + " 'zb': 0.0,\n", + " 'IA_model': 'eNLA',\n", + " 'AIA': 1.72,\n", + " 'betaIA': 2.17,\n", + " 'etaIA': -0.41}" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cosmoFM_fid.allparams" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['GCph']" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cosmoFM_fid.observables" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "In class: FisherMatrix ----> Computing photo Fisher matrix\n", + "\n", + "Computing fiducial\n", + "\n", + "Fiducial generated in 0.55 s\n", + "\n", + "Noise added to fiducial\n", + "\n", + "Noisy Cls generated in 0.00 s\n", + "\n", + "Computed covariance matrix\n", + "\n", + "Covmat of Cls generated in 0.13 s\n", + "\n", + "Total calculation in 0.69 s\n", + ">> computing derivs >>\n", + "\n", + " +++ Computing derivative on Omegam\n", + "\n", + "In class: derivatives Derivative on Omegam done! in : 1.31 s\n", + "\n", + " +++ Computing derivative on Omegab\n", + "\n", + "In class: derivatives Derivative on Omegab done! in : 1.30 s\n", + "\n", + " +++ Computing derivative on h\n", + "\n", + "In class: derivatives Derivative on h done! in : 1.24 s\n", + "\n", + " +++ Computing derivative on ns\n", + "\n", + "In class: derivatives Derivative on ns done! in : 1.21 s\n", + "\n", + " +++ Computing derivative on sigma8\n", + "\n", + "In class: derivatives Derivative on sigma8 done! in : 1.22 s\n", + "\n", + " +++ Computing derivative on b1\n", + "\n", + "In class: derivatives Derivative on b1 done! in : 1.04 s\n", + "\n", + " +++ Computing derivative on b2\n", + "\n", + "In class: derivatives Derivative on b2 done! in : 1.02 s\n", + "\n", + " +++ Computing derivative on b3\n", + "\n", + "In class: derivatives Derivative on b3 done! in : 1.08 s\n", + "\n", + " +++ Computing derivative on b4\n", + "\n", + "In class: derivatives Derivative on b4 done! in : 1.04 s\n", + "\n", + " +++ Computing derivative on b5\n", + "\n", + "In class: derivatives Derivative on b5 done! in : 1.04 s\n", + "\n", + " +++ Computing derivative on b6\n", + "\n", + "In class: derivatives Derivative on b6 done! in : 1.03 s\n", + "\n", + " +++ Computing derivative on b7\n", + "\n", + "In class: derivatives Derivative on b7 done! in : 1.04 s\n", + "\n", + " +++ Computing derivative on b8\n", + "\n", + "In class: derivatives Derivative on b8 done! in : 1.03 s\n", + "\n", + " +++ Computing derivative on b9\n", + "\n", + "In class: derivatives Derivative on b9 done! in : 1.03 s\n", + "\n", + " +++ Computing derivative on b10\n", + "\n", + "In class: derivatives Derivative on b10 done! in : 1.02 s\n", + "\n", + "In class: PhotoCov -->> Derivatives computed in 16.63 s\n", + "\n", + "In class: FisherMatrix Computing Fisher matrix\n", + "\n", + " Finished calculation of Fisher Matrix for ['GCph'] in: 0.01 s\n", + "\n", + " Fisher matrix calculation finished in 17.57 s\n", + "\n", + " Fisher matrix exported: ./results/CosmicFish_v1.2.4_cosmicjellyfish_Euclid-ISTF-Pess-GCphoto_symb_simple_GCph_fishermatrix.txt\n", + "\n", + "\n", + "In class: FisherMatrix CosmicFish settings and Survey specifications exported: ./results/CosmicFish_v1.2.4_cosmicjellyfish_Euclid-ISTF-Pess-GCphoto_symb_simple_GCph_fishermatrix_specifications.dat\n" + ] + } + ], + "source": [ + "#fishmat_photo = cosmoFM_fid.compute()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "photo_like = PhotometricLikelihood(cosmo_data=cosmoFM_fid, cosmo_theory=cosmoFM_fid)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "#print(Cells_fid[\"Cell_LL\"].shape)\n", + "#print(Cells_fid[\"Cell_GG\"].shape)\n", + "#print(Cells_fid[\"Cell_GL\"].shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "750\n" + ] + } + ], + "source": [ + "ellmax_GC = cosmoFM_fid.specs[\"lmax_GCph\"]\n", + "print(ellmax_GC)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "samp1dic = cosmoFM_fid.allparams.copy()" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Sample likelihood 4.330972755260946e-11\n", + "Sample likelihood -1568.703678015382\n", + "Sample likelihood -1522.4374216102938\n" + ] + } + ], + "source": [ + "print(\"Sample likelihood\", photo_like.loglike(param_dict=samp1dic))\n", + "samp1dic = cosmoFM_fid.allparams\n", + "\n", + "samp1dic['b1'] = 1.0\n", + "\n", + "print(\"Sample likelihood\", photo_like.loglike(param_dict=samp1dic))\n", + "\n", + "samp1dic['h'] = 0.7\n", + "\n", + "print(\"Sample likelihood\", photo_like.loglike(param_dict=samp1dic))" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "prior_nonuis = Prior()\n", + "prior_withnuis = Prior()" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "prior_dict ={\n", + " 'Omegam': [0.24, 0.4],\n", + " 'Omegab': [0.04, 0.06],\n", + " 'h': [0.61, 0.75],\n", + " 'ns': [0.92, 1.00],\n", + " 'sigma8': [0.79, 0.83],\n", + " 'AIA': [1.0, 3.0],\n", + " 'etaIA' :[-6.0, 6.0],\n", + " 'b1': [1.0, 3.0],\n", + " 'b2': [1.0, 3.0],\n", + " 'b3': [1.0, 3.0],\n", + " 'b4': [1.0, 3.0],\n", + " 'b5': [1.0, 3.0],\n", + " 'b6': [1.0, 3.0],\n", + " 'b7': [1.0, 3.0],\n", + " 'b8': [1.0, 3.0],\n", + " 'b9': [1.0, 3.0],\n", + " 'b10': [1.0, 3.0]\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'Omegam': 0.01,\n", + " 'Omegab': 0.01,\n", + " 'h': 0.01,\n", + " 'ns': 0.01,\n", + " 'sigma8': 0.01,\n", + " 'b1': 0.06,\n", + " 'b2': 0.06,\n", + " 'b3': 0.06,\n", + " 'b4': 0.06,\n", + " 'b5': 0.06,\n", + " 'b6': 0.06,\n", + " 'b7': 0.06,\n", + " 'b8': 0.06,\n", + " 'b9': 0.06,\n", + " 'b10': 0.06}" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cosmoFM_fid.freeparams" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "for par in prior_dict.keys():\n", + " if par in cosmoFM_fid.freeparams.keys():\n", + " dist_prior = (prior_dict[par][0], prior_dict[par][1])\n", + " if re.match(r'b\\d+', par):\n", + " prior_withnuis.add_parameter(par, dist_prior)\n", + " elif re.search(r'IA', par):\n", + " prior_withnuis.add_parameter(par, dist_prior)\n", + " else:\n", + " prior_nonuis.add_parameter(par, dist_prior)\n", + " prior_withnuis.add_parameter(par, dist_prior)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['Omegam', 'Omegab', 'h', 'ns', 'sigma8']\n", + "['Omegam', 'Omegab', 'h', 'ns', 'sigma8', 'b1', 'b2', 'b3', 'b4', 'b5', 'b6', 'b7', 'b8', 'b9', 'b10']\n" + ] + } + ], + "source": [ + "print(prior_nonuis.keys)\n", + "print(prior_withnuis.keys)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "use_nuisance = True" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Loading prior with keys: ['Omegam', 'Omegab', 'h', 'ns', 'sigma8', 'b1', 'b2', 'b3', 'b4', 'b5', 'b6', 'b7', 'b8', 'b9', 'b10']\n" + ] + } + ], + "source": [ + "if use_nuisance: \n", + " prior_chosen = prior_withnuis\n", + "else:\n", + " prior_chosen = prior_nonuis\n", + "print(\"Loading prior with keys: \", prior_chosen.keys)" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Starting sampler at 2025-09-29 13:00:47\n" + ] + } + ], + "source": [ + "tini = time.time()\n", + "print(\"Starting sampler at\", time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(tini)))" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Starting the nautilus sampler...\n", + "Please report issues at github.com/johannesulf/nautilus.\n", + "Status | Bounds | Ellipses | Networks | Calls | f_live | N_eff | log Z \n", + "Computing | 10 | 1 | 8 | 12800 | 1.0000 | 1 | -7351.37 \r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Process ForkPoolWorker-7:\n", + "Process ForkPoolWorker-8:\n", + "Process ForkPoolWorker-2:\n", + "Process ForkPoolWorker-5:\n", + "Process ForkPoolWorker-1:\n", + "Process ForkPoolWorker-3:\n", + "Process ForkPoolWorker-4:\n", + "Process ForkPoolWorker-6:\n", + "Traceback (most recent call last):\n", + "Traceback (most recent call last):\n", + "Traceback (most recent call last):\n", + "Traceback (most recent call last):\n", + "Traceback (most recent call last):\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 314, in _bootstrap\n", + " self.run()\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 314, in _bootstrap\n", + " self.run()\n", + "Traceback (most recent call last):\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 314, in _bootstrap\n", + " self.run()\n", + "Traceback (most recent call last):\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 108, in run\n", + " self._target(*self._args, **self._kwargs)\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 314, in _bootstrap\n", + " self.run()\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 108, in run\n", + " self._target(*self._args, **self._kwargs)\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 314, in _bootstrap\n", + " self.run()\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 108, in run\n", + " self._target(*self._args, **self._kwargs)\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 314, in _bootstrap\n", + " self.run()\n", + "Traceback (most recent call last):\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 125, in worker\n", + " result = (True, func(*args, **kwds))\n", + " ^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 125, in worker\n", + " result = (True, func(*args, **kwds))\n", + " ^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 108, in run\n", + " self._target(*self._args, **self._kwargs)\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 314, in _bootstrap\n", + " self.run()\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 108, in run\n", + " self._target(*self._args, **self._kwargs)\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 108, in run\n", + " self._target(*self._args, **self._kwargs)\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 48, in mapstar\n", + " return list(map(*args))\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 125, in worker\n", + " result = (True, func(*args, **kwds))\n", + " ^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 48, in mapstar\n", + " return list(map(*args))\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 108, in run\n", + " self._target(*self._args, **self._kwargs)\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 314, in _bootstrap\n", + " self.run()\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 125, in worker\n", + " result = (True, func(*args, **kwds))\n", + " ^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 125, in worker\n", + " result = (True, func(*args, **kwds))\n", + " ^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 48, in mapstar\n", + " return list(map(*args))\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/pool.py\", line 33, in likelihood_worker\n", + " return LIKELIHOOD(*args)\n", + " ^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/pool.py\", line 33, in likelihood_worker\n", + " return LIKELIHOOD(*args)\n", + " ^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 125, in worker\n", + " result = (True, func(*args, **kwds))\n", + " ^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 125, in worker\n", + " result = (True, func(*args, **kwds))\n", + " ^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 48, in mapstar\n", + " return list(map(*args))\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/process.py\", line 108, in run\n", + " self._target(*self._args, **self._kwargs)\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 48, in mapstar\n", + " return list(map(*args))\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/base.py\", line 83, in loglike\n", + " The theory prediction in the same format as the observed data\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/base.py\", line 83, in loglike\n", + " The theory prediction in the same format as the observed data\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 48, in mapstar\n", + " return list(map(*args))\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/pool.py\", line 33, in likelihood_worker\n", + " return LIKELIHOOD(*args)\n", + " ^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/pool.py\", line 33, in likelihood_worker\n", + " return LIKELIHOOD(*args)\n", + " ^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 125, in worker\n", + " result = (True, func(*args, **kwds))\n", + " ^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 167, in compute_theory\n", + " photo_cls = pobs.ComputeCls(\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 174, in compute_theory\n", + " return _cells_from_cls(photo_cls, self.photo_cov_data, self.observables)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/pool.py\", line 33, in likelihood_worker\n", + " return LIKELIHOOD(*args)\n", + " ^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/base.py\", line 83, in loglike\n", + " The theory prediction in the same format as the observed data\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/pool.py\", line 33, in likelihood_worker\n", + " return LIKELIHOOD(*args)\n", + " ^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/base.py\", line 83, in loglike\n", + " The theory prediction in the same format as the observed data\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 209, in __init__\n", + " self.window = photo_window.GalaxyPhotoDist(self.photopars)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 48, in mapstar\n", + " return list(map(*args))\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 34, in _cells_from_cls\n", + " photo_cls.compute_all()\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/base.py\", line 83, in loglike\n", + " The theory prediction in the same format as the observed data\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 174, in compute_theory\n", + " return _cells_from_cls(photo_cls, self.photo_cov_data, self.observables)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/pool.py\", line 33, in likelihood_worker\n", + " return LIKELIHOOD(*args)\n", + " ^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/base.py\", line 83, in loglike\n", + " The theory prediction in the same format as the observed data\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_window.py\", line 57, in __init__\n", + " self.normalization = {\"GCph\": self.norm(\"GCph\"), \"WL\": self.norm(\"WL\")}\n", + " ^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 276, in compute_all\n", + " self.P_limber()\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 167, in compute_theory\n", + " photo_cls = pobs.ComputeCls(\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/miniforge3/lib/python3.12/multiprocessing/pool.py\", line 48, in mapstar\n", + " return list(map(*args))\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 167, in compute_theory\n", + " photo_cls = pobs.ComputeCls(\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 34, in _cells_from_cls\n", + " photo_cls.compute_all()\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/base.py\", line 83, in loglike\n", + " The theory prediction in the same format as the observed data\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_window.py\", line 213, in norm\n", + " trapezoid([self.ngal_photoz(z, i, obs) for z in zint], dx=dz)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 174, in compute_theory\n", + " return _cells_from_cls(photo_cls, self.photo_cov_data, self.observables)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/pool.py\", line 33, in likelihood_worker\n", + " return LIKELIHOOD(*args)\n", + " ^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 358, in P_limber\n", + " * self.cosmo.matpow(\n", + " ^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 165, in __init__\n", + " self.cosmo = cosmology.cosmo_functions(cosmopars, cfg.input_type)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 277, in compute_all\n", + " self.sqrtP_limber()\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 165, in __init__\n", + " self.cosmo = cosmology.cosmo_functions(cosmopars, cfg.input_type)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_window.py\", line 157, in ngal_photoz\n", + " (np.sqrt(1 / 2) * (z - self.photo[\"zo\"] - self.photo[\"co\"] * z_bins[i]))\n", + " ^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1562, in __init__\n", + " symbresults = boltzmann_code(cosmopars, code=\"symbolic\")\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1637, in matpow\n", + " def matpow(self, z, k, nonlinear=False, tracer=\"matter\"):\n", + " \n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 174, in compute_theory\n", + " return _cells_from_cls(photo_cls, self.photo_cov_data, self.observables)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/base.py\", line 83, in loglike\n", + " The theory prediction in the same format as the observed data\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 34, in _cells_from_cls\n", + " photo_cls.compute_all()\n", + "KeyboardInterrupt\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 277, in compute_all\n", + " self.sqrtP_limber()\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1562, in __init__\n", + " symbresults = boltzmann_code(cosmopars, code=\"symbolic\")\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 34, in _cells_from_cls\n", + " photo_cls.compute_all()\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 158, in __init__\n", + " self.symbolic_results()\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/likelihood/photo_like.py\", line 167, in compute_theory\n", + " photo_cls = pobs.ComputeCls(\n", + " ^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 397, in sqrtP_limber\n", + " self.cosmo.SigmaMG(self.z[zin], kappa[ell, zin])\n", + "KeyboardInterrupt\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 301, in compute_all\n", + " self.result = self.computecls_vectorized()\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1934, in SigmaMG\n", + " def SigmaMG(self, z, k):\n", + " \n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 399, in sqrtP_limber\n", + " self.cosmo.matpow(\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/utilities/utils.py\", line 85, in wrapper\n", + " return func(*args, **kwargs)\n", + " ^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 165, in __init__\n", + " self.cosmo = cosmology.cosmo_functions(cosmopars, cfg.input_type)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 158, in __init__\n", + " self.symbolic_results()\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/LSSsurvey/photo_obs.py\", line 730, in computecls_vectorized\n", + " hub = self.cosmo.Hubble(self.z)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "KeyboardInterrupt\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1092, in symbolic_results\n", + " D_kz = np.array([self.symbcosmo.growthFactor(self.zgrid) for kk in self.kgrid_1Mpc]).T\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1562, in __init__\n", + " symbresults = boltzmann_code(cosmopars, code=\"symbolic\")\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1675, in matpow\n", + " return self.Pmm(z, k, nonlinear=nonlinear)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/utilities/utils.py\", line 85, in wrapper\n", + " return func(*args, **kwargs)\n", + " ^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1592, in Hubble\n", + " hubble = prefactor * self.results.h_of_z(z)\n", + " ^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 2077, in growthFactor\n", + " D = self._zFunction('lnzp1_growthfactor', self._growthFactorExact, z, derivative = derivative,\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 158, in __init__\n", + " self.symbolic_results()\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1135, in symbolic_results\n", + " Pk_nl = vectorized_run_halofit(self.zgrid)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1691, in Pmm\n", + " power = self.results.Pk_nl(z, k, grid=False)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 1123, in _zFunction\n", + " ret = interpolator(x, nu = derivative)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2544, in __call__\n", + " return self._call_as_normal(*args, **kwargs)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/utilities/utils.py\", line 85, in wrapper\n", + " return func(*args, **kwargs)\n", + " ^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2544, in __call__\n", + " return self._call_as_normal(*args, **kwargs)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/scipy/interpolate/_fitpack2.py\", line 974, in __call__\n", + " def __call__(self, x, y, dx=0, dy=0, grid=True):\n", + " \n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2537, in _call_as_normal\n", + " return self._vectorize_call(func=func, args=vargs)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1135, in symbolic_results\n", + " Pk_nl = vectorized_run_halofit(self.zgrid)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/scipy/interpolate/_fitpack2.py\", line 369, in __call__\n", + " def __call__(self, x, nu=0, ext=None):\n", + " \n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2537, in _call_as_normal\n", + " return self._vectorize_call(func=func, args=vargs)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "KeyboardInterrupt\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2544, in __call__\n", + " return self._call_as_normal(*args, **kwargs)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2625, in _vectorize_call\n", + " outputs = ufunc(*args, out=...)\n", + " ^^^^^^^^^^^^^^^^^^^^^\n", + "KeyboardInterrupt\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2618, in _vectorize_call\n", + " res = self._vectorize_call_with_signature(func, args)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1081, in \n", + " lambda zz: self.symbcosmo.Hz(zz) / cosmo_functions.c\n", + " ^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2537, in _call_as_normal\n", + " return self._vectorize_call(func=func, args=vargs)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2656, in _vectorize_call_with_signature\n", + " results = func(*(arg[index] for arg in args))\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2618, in _vectorize_call\n", + " res = self._vectorize_call_with_signature(func, args)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 927, in Hz\n", + " H = self.Ez(z) * self.H0\n", + " ^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1119, in vectorized_halofit\n", + " return ((1 / self.h_now) ** 3) * self.symbfit.run_halofit(\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 899, in Ez\n", + " if np.any(t_array < 0.0):\n", + " ^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py\", line 2656, in _vectorize_call_with_signature\n", + " results = func(*(arg[index] for arg in args))\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/symbolic_pofk/syrenhalofit.py\", line 201, in run_halofit\n", + " plin = linear.plin_emulated(\n", + " ^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/_core/fromnumeric.py\", line 2580, in any\n", + " return _wrapreduction_any_all(a, np.logical_or, 'any', axis, out,\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/cosmicfishpie/cosmicfishpie/cosmology/cosmology.py\", line 1119, in vectorized_halofit\n", + " return ((1 / self.h_now) ** 3) * self.symbfit.run_halofit(\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/symbolic_pofk/linear.py\", line 383, in plin_emulated\n", + " D = cosmo.growthFactor(1/a - 1)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/symbolic_pofk/syrenhalofit.py\", line 201, in run_halofit\n", + " plin = linear.plin_emulated(\n", + " ^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/numpy/_core/fromnumeric.py\", line 102, in _wrapreduction_any_all\n", + " return ufunc.reduce(obj, axis, bool, out, **passkwargs)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/symbolic_pofk/linear.py\", line 383, in plin_emulated\n", + " D = cosmo.growthFactor(1/a - 1)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "KeyboardInterrupt\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 2077, in growthFactor\n", + " D = self._zFunction('lnzp1_growthfactor', self._growthFactorExact, z, derivative = derivative,\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 2077, in growthFactor\n", + " D = self._zFunction('lnzp1_growthfactor', self._growthFactorExact, z, derivative = derivative,\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 1077, in _zFunction\n", + " interpolator = self._zInterpolator(table_name, func, inverse = inverse, future = future)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 1077, in _zFunction\n", + " interpolator = self._zInterpolator(table_name, func, inverse = inverse, future = future)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 1028, in _zInterpolator\n", + " def _zInterpolator(self, table_name, func, inverse = False, future = True):\n", + " \n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 1044, in _zInterpolator\n", + " x_table = func(z_table)\n", + " ^^^^^^^^^^^^^\n", + "KeyboardInterrupt\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 2035, in _growthFactorExact\n", + " D = self.growthFactorUnnormalized(z) / self.growthFactorUnnormalized(0.0)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 2003, in growthFactorUnnormalized\n", + " D[mask1] = 5.0 / 2.0 * self.Om0 * Ez_D(z1) * self._integral(integrand, z1, np.inf)\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 996, in _integral\n", + " integ[i], _ = scipy.integrate.quad(integrand, z_min_use[i], z_max_use[i])\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/scipy/integrate/_quadpack_py.py\", line 460, in quad\n", + " retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/scipy/integrate/_quadpack_py.py\", line 609, in _quad\n", + " return _quadpack._qagie(func, bound, infbounds, args, full_output,\n", + " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + " File \"/home/casas/Cosmo/.cfpenv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py\", line 1920, in integrand\n", + " def integrand(z):\n", + " \n", + "KeyboardInterrupt\n" + ] + }, + { + "ename": "KeyboardInterrupt", + "evalue": "", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[29], line 12\u001b[0m\n\u001b[1;32m 1\u001b[0m sampler \u001b[38;5;241m=\u001b[39m Sampler(prior_chosen, \n\u001b[1;32m 2\u001b[0m photo_like\u001b[38;5;241m.\u001b[39mloglike, \n\u001b[1;32m 3\u001b[0m n_live\u001b[38;5;241m=\u001b[39msampler_settings[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mn_live\u001b[39m\u001b[38;5;124m\"\u001b[39m], \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 10\u001b[0m likelihood_kwargs\u001b[38;5;241m=\u001b[39m{\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mprior\u001b[39m\u001b[38;5;124m'\u001b[39m: prior_chosen}\n\u001b[1;32m 11\u001b[0m )\n\u001b[0;32m---> 12\u001b[0m \u001b[43msampler\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43mverbose\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdiscard_exploration\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m 13\u001b[0m log_z_all \u001b[38;5;241m=\u001b[39m sampler\u001b[38;5;241m.\u001b[39mevidence()\n\u001b[1;32m 14\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mEvidence:\u001b[39m\u001b[38;5;124m'\u001b[39m, log_z_all)\n", + "File \u001b[0;32m~/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/sampler.py:441\u001b[0m, in \u001b[0;36mSampler.run\u001b[0;34m(self, f_live, n_shell, n_eff, n_like_max, discard_exploration, timeout, verbose)\u001b[0m\n\u001b[1;32m 438\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfilepath \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m 439\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mwrite(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfilepath, overwrite\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[0;32m--> 441\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn_update_iter \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43madd_samples\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mverbose\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mverbose\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 442\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn_like_iter \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn_batch\n\u001b[1;32m 443\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfilepath \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m 444\u001b[0m \u001b[38;5;66;03m# Write the complete file if this is the first batch.\u001b[39;00m\n", + "File \u001b[0;32m~/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/sampler.py:1119\u001b[0m, in \u001b[0;36mSampler.add_samples\u001b[0;34m(self, shell, verbose)\u001b[0m\n\u001b[1;32m 1116\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mprint_status(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mComputing\u001b[39m\u001b[38;5;124m'\u001b[39m, end\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;130;01m\\r\u001b[39;00m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 1118\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mshell_n_sample[shell] \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m n_bound\n\u001b[0;32m-> 1119\u001b[0m log_l, blobs \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mevaluate_likelihood\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpoints\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 1120\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpoints[shell] \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mappend(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpoints[shell], points, axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m)\n\u001b[1;32m 1121\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlog_l[shell] \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mappend(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlog_l[shell], log_l, axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m)\n", + "File \u001b[0;32m~/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/sampler.py:869\u001b[0m, in \u001b[0;36mSampler.evaluate_likelihood\u001b[0;34m(self, points)\u001b[0m\n\u001b[1;32m 867\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlist\u001b[39m(\u001b[38;5;28mzip\u001b[39m(\u001b[38;5;241m*\u001b[39mresult))\n\u001b[1;32m 868\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpool_l \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m--> 869\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlist\u001b[39m(\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpool_l\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmap\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlikelihood\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43margs\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[1;32m 870\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 871\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlist\u001b[39m(\u001b[38;5;28mmap\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlikelihood, args))\n", + "File \u001b[0;32m~/Cosmo/.cfpenv/lib/python3.12/site-packages/nautilus/pool.py:85\u001b[0m, in \u001b[0;36mNautilusPool.map\u001b[0;34m(self, func, iterable)\u001b[0m\n\u001b[1;32m 83\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mlist\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpool\u001b[38;5;241m.\u001b[39mgather(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpool\u001b[38;5;241m.\u001b[39mmap(func, iterable)))\n\u001b[1;32m 84\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m---> 85\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mlist\u001b[39m(\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpool\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmap\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43miterable\u001b[49m\u001b[43m)\u001b[49m)\n", + "File \u001b[0;32m~/miniforge3/lib/python3.12/multiprocessing/pool.py:367\u001b[0m, in \u001b[0;36mPool.map\u001b[0;34m(self, func, iterable, chunksize)\u001b[0m\n\u001b[1;32m 362\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21mmap\u001b[39m(\u001b[38;5;28mself\u001b[39m, func, iterable, chunksize\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m 363\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m'''\u001b[39;00m\n\u001b[1;32m 364\u001b[0m \u001b[38;5;124;03m Apply `func` to each element in `iterable`, collecting the results\u001b[39;00m\n\u001b[1;32m 365\u001b[0m \u001b[38;5;124;03m in a list that is returned.\u001b[39;00m\n\u001b[1;32m 366\u001b[0m \u001b[38;5;124;03m '''\u001b[39;00m\n\u001b[0;32m--> 367\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_map_async\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43miterable\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmapstar\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mchunksize\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/miniforge3/lib/python3.12/multiprocessing/pool.py:768\u001b[0m, in \u001b[0;36mApplyResult.get\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m 767\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21mget\u001b[39m(\u001b[38;5;28mself\u001b[39m, timeout\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[0;32m--> 768\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwait\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtimeout\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 769\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mready():\n\u001b[1;32m 770\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTimeoutError\u001b[39;00m\n", + "File \u001b[0;32m~/miniforge3/lib/python3.12/multiprocessing/pool.py:765\u001b[0m, in \u001b[0;36mApplyResult.wait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m 764\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21mwait\u001b[39m(\u001b[38;5;28mself\u001b[39m, timeout\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[0;32m--> 765\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_event\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwait\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtimeout\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/miniforge3/lib/python3.12/threading.py:655\u001b[0m, in \u001b[0;36mEvent.wait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m 653\u001b[0m signaled \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_flag\n\u001b[1;32m 654\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m signaled:\n\u001b[0;32m--> 655\u001b[0m signaled \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_cond\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwait\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtimeout\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 656\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m signaled\n", + "File \u001b[0;32m~/miniforge3/lib/python3.12/threading.py:355\u001b[0m, in \u001b[0;36mCondition.wait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m 353\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m: \u001b[38;5;66;03m# restore state no matter what (e.g., KeyboardInterrupt)\u001b[39;00m\n\u001b[1;32m 354\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m timeout \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m--> 355\u001b[0m \u001b[43mwaiter\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43macquire\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 356\u001b[0m gotit \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[1;32m 357\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n", + "\u001b[0;31mKeyboardInterrupt\u001b[0m: " + ] + } + ], + "source": [ + "sampler = Sampler(prior_chosen, \n", + " photo_like.loglike, \n", + " n_live=sampler_settings[\"n_live\"], \n", + " n_networks=sampler_settings[\"n_networks\"], \n", + " n_batch=sampler_settings[\"n_batch\"], \n", + " pool=sampler_settings[\"pool\"], \n", + " pass_dict=False,\n", + " filepath=options[\"outroot\"]+\".hdf5\", \n", + " resume=True,\n", + " likelihood_kwargs={'prior': prior_chosen}\n", + " )\n", + "sampler.run(verbose=True, discard_exploration=True)\n", + "log_z_all = sampler.evidence()\n", + "print('Evidence:', log_z_all)\n", + "points_all, log_w_all, log_l_all = sampler.posterior()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tfin = time.time()\n", + "elapsed = tfin - tini\n", + "hours = int(elapsed // 3600)\n", + "minutes = int((elapsed % 3600) // 60)\n", + "seconds = int(elapsed % 60)\n", + "print(\"Sampler finished at\", time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(tfin)))\n", + "print(f\"Total time elapsed: {hours:02d}:{minutes:02d}:{seconds:02d}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sample_wghlkl = (np.vstack((points_all.T, np.exp(log_w_all), log_l_all)).T)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "outfile_chain = options[\"outroot\"]+\".txt\"\n", + "print(f\"Saving chain to text file {outfile_chain}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "headerlist = ['loglike', 'weights'] + list(prior_chosen.keys)\n", + "header = \" \".join(headerlist)\n", + "print(\"Saving header: \", header)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.savetxt(outfile_chain, sample_wghlkl, header=header)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".cfpenv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.11" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tests/test_photo_like.py b/tests/test_photo_like.py new file mode 100644 index 0000000..1301c44 --- /dev/null +++ b/tests/test_photo_like.py @@ -0,0 +1,89 @@ +import math + +import numpy as np +import pytest + +from cosmicfishpie.fishermatrix import cosmicfish +from cosmicfishpie.likelihood import PhotometricLikelihood + + +@pytest.fixture(scope="module") +def photometric_fiducial_obs(): + options = { + "accuracy": 1, + "feedback": 1, + "code": "symbolic", + "outroot": "test_photo_like", + "survey_name": "Euclid", + "survey_name_photo": "Euclid-Photometric-ISTF-Pessimistic", + "survey_name_spectro": False, + "specs_dir": "../cosmicfishpie/configs/default_survey_specifications/", + "cosmo_model": "LCDM", + } + + fiducial = { + "Omegam": 0.3145714273, + "Omegab": 0.0491989, + "h": 0.6737, + "ns": 0.96605, + "sigma8": 0.81, + "w0": -1.0, + "wa": 0.0, + "mnu": 0.06, + "Neff": 3.044, + } + + observables = ["GCph"] + + return cosmicfish.FisherMatrix( + fiducialpars=fiducial, + options=options, + observables=observables, + cosmoModel=options["cosmo_model"], + surveyName=options["survey_name"], + ) + + +@pytest.fixture(scope="module") +def photometric_likelihood(photometric_fiducial_obs): + return PhotometricLikelihood( + cosmo_data=photometric_fiducial_obs, + cosmo_theory=photometric_fiducial_obs, + observables=photometric_fiducial_obs.observables, + ) + + +def _sample_params(photometric_likelihood): + fiducial_obs = photometric_likelihood.cosmo_data + samp_pars = fiducial_obs.allparams.copy() + return samp_pars + + +def test_photometric_cells_have_expected_shape(photometric_likelihood): + cells = photometric_likelihood.data_obs + assert "ells" in cells + assert "Cell_GG" in cells + # assert "Cell_LL" in cells + # assert "Cell_GL" in cells + assert cells["Cell_GG"].shape[0] == len(cells["ells"]) + + +def test_photometric_loglike_matches_notebook_value(photometric_likelihood): + sample_params = _sample_params(photometric_likelihood) + loglike_value = photometric_likelihood.loglike(param_dict=sample_params) + expected = 4.3309e-11 + assert math.isclose(loglike_value, expected, rel_tol=1e-2) + + +def test_photometric_cell_entry_matches_theory(photometric_likelihood): + sample_params = _sample_params(photometric_likelihood) + ells = photometric_likelihood.data_obs["ells"] + target_ell = 300 + idx = int(np.argmin(np.abs(ells - target_ell))) + + data_val = photometric_likelihood.data_obs["Cell_GG"][idx, 1, 8] + theory_cells = photometric_likelihood.compute_theory(dict(sample_params)) + theory_val = theory_cells["Cell_GG"][idx, 1, 8] + nb_val = 9.974414863747675e-14 + assert theory_val == data_val + assert theory_val == pytest.approx(nb_val, rel=1e-12, abs=1e-18) diff --git a/tests/test_spectro_like.py b/tests/test_spectro_like.py new file mode 100644 index 0000000..f9d9753 --- /dev/null +++ b/tests/test_spectro_like.py @@ -0,0 +1,79 @@ +import numpy as np +import pytest + +from cosmicfishpie.likelihood.spectro_like import ( + compute_chi2_legendre, + compute_wedge_chi2, + legendre_Pgg, + loglike, +) + + +class DummyFMLegendre: + def __init__(self, mu_grid, ksamp): + # mu_grid array over which to integrate + self.Pk_mugrid = np.array(mu_grid) + self.Pk_musamp = len(mu_grid) + # number of k-samples (unused except for shape checks) + self.Pk_ksamp = ksamp + + +def test_legendre_Pgg_constant_field(): + # Create dummy FisherMatrix with mu grid symmetric [-1, 1] + n_mu = 101 + n_k = 3 + mu = np.linspace(-1.0, 1.0, n_mu) + fm = DummyFMLegendre(mu_grid=mu, ksamp=n_k) + + # constant field C=1 for two redshift bins + n_z = 2 + obs_Pgg = np.ones((n_z, n_mu, n_k)) + # Compute multipoles + P_ell = legendre_Pgg(obs_Pgg, fm) + # Expect P0 = (2*0+1) * ∫1 dµ = 2, P2 and P4 ~ 0 + # P_ell shape: (n_k, n_z, n_ell=3) + assert P_ell.shape == (n_k, n_z, 3) + # Monopole ~2, quadrupole and hexadecapole ~0 + assert np.allclose(P_ell[:, :, 0], 2.0, atol=1e-3) + assert np.allclose(P_ell[:, :, 1:], 0.0, atol=1e-3) + + +def test_compute_chi2_legendre_simple(): + # Single k, single z, single ell + P_data = np.array([[[2.0]]]) + P_th = np.array([[[1.0]]]) + inv_cov = np.array([[[[0.25]]]]) # variance = 4 + chi2 = compute_chi2_legendre(P_data, P_th, inv_cov) + # (2-1)^2 * 0.25 = 0.25 + assert chi2 == pytest.approx(0.25) + + +class DummyFMWedge: + def __init__(self, k_grid, mu_grid, vol_array): + # k and mu grids for integration + self.Pk_kgrid = np.array(k_grid) + self.Pk_mugrid = np.array(mu_grid) + # volume_survey_array indexed by redshift bins + self.pk_cov = type("pc", (), {"volume_survey_array": np.array(vol_array)}) + + +def test_compute_wedge_chi2_constant_simple(): + # Set up simple case with one redshift bin, two mu and two k samples + k_grid = np.array([1.0, 2.0]) + mu_grid = np.array([0.0, 1.0]) + vol = np.array([2.0]) # one redshift bin + fm = DummyFMWedge(k_grid=k_grid, mu_grid=mu_grid, vol_array=vol) + + # Observed = 1, Theory = 0 + data = np.ones((1, len(mu_grid), len(k_grid))) + theory = np.zeros_like(data) + chi2 = compute_wedge_chi2(data, theory, fm) + # Analytical result: sum_z ∫_mu ∫_k delta^2 / sigma = 5/(4*pi^2) + expected = 5.0 / (4.0 * np.pi**2) + assert chi2 == pytest.approx(expected, rel=1e-3) + + +def test_loglike_no_inputs_returns_minus_inf(): + # loglike with no theory and no params gives -inf + ll = loglike() + assert ll == -np.inf