diff --git a/papers/bmodes/config/config.yaml b/papers/bmodes/config/config.yaml index 3c8e7a44..47b385bd 100644 --- a/papers/bmodes/config/config.yaml +++ b/papers/bmodes/config/config.yaml @@ -63,7 +63,7 @@ covariance: use_semianalytic: true # Number of Monte Carlo samples for covariance propagation (used by pure E/B mode analysis) n_samples: 2000 - # Cosmology: imported from sp_validation.cosmology.PLANCK18 (astropy Planck18) + # Cosmology: imported from cs_util.cosmo.PLANCK18 (astropy Planck18) # Mask Cl paths: defined in covariance.smk (MASK_CLS_FILES) # Harmonic-space (C_l) analysis parameters diff --git a/papers/bmodes/config/covariance.md b/papers/bmodes/config/covariance.md index 72081e81..54e7ad42 100644 --- a/papers/bmodes/config/covariance.md +++ b/papers/bmodes/config/covariance.md @@ -20,7 +20,7 @@ Gaussian-only approximation: non-Gaussian contributions computationally prohibit |-----------|--------|-------------| | Samples | `config.yaml: covariance.n_samples` | MC samples for propagation (default 2000) | | Use masked | `config.yaml: covariance.default_masked` | Whether to use masked covariance | -| Cosmology | `Snakefile: PLANCK18` | astropy Planck18 via sp_validation.cosmology | +| Cosmology | `Snakefile: PLANCK18` | astropy Planck18 via cs_util.cosmo | | Mask Cls | `covariance.smk: MASK_CLS_FILES` | Per-version mask power spectra | ## Survey Properties diff --git a/papers/bmodes/scripts/precompute_pure_eb_chunk.py b/papers/bmodes/scripts/precompute_pure_eb_chunk.py index 9672c949..363a3192 100644 --- a/papers/bmodes/scripts/precompute_pure_eb_chunk.py +++ b/papers/bmodes/scripts/precompute_pure_eb_chunk.py @@ -79,9 +79,9 @@ def compute_chunk( cosmo_params=None, ): from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes + from cs_util.cosmo import PLANCK18, get_theo_xi from sp_validation.cosmo_val import CosmologyValidation - from sp_validation.cosmology import PLANCK18, get_theo_xi if cosmo_params is None: cosmo_params = dict(PLANCK18) diff --git a/pyproject.toml b/pyproject.toml index 7b8f1ed6..bb1fe413 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,12 @@ dependencies = [ "camb>=1.6", "clmm", "colorama", - "cs_util>=0.2.1", + # Track cs_util's develop branch directly (git dependency) rather than a + # PyPI pin: the two repos are iterating together heavily and cs_util + # releases are infrequent. This PR's cosmology repoint needs get_cosmo / + # get_theo_c_ell / get_theo_xi / PLANCK18 in cs_util.cosmo, which land via + # CosmoStat/cs_util#76 — so this goes green once #76 merges into develop. + "cs_util @ git+https://github.com/CosmoStat/cs_util.git@develop", "emcee", # SHA-pinned snapshot of getdist branch `upper_triangle_whisker`. "getdist @ git+https://github.com/benabed/getdist.git@113cd22a9a0d013b6f72fe734be81f260f3d3be5", diff --git a/scratch/guerrini/namaster_utils.py b/scratch/guerrini/namaster_utils.py index 52dbdb90..43b96083 100644 --- a/scratch/guerrini/namaster_utils.py +++ b/scratch/guerrini/namaster_utils.py @@ -23,10 +23,10 @@ import pymaster as nmt import seaborn as sns from astropy.io import fits +from cs_util.cosmo import get_cosmo, get_theo_c_ell from matplotlib import scale as mscale from tqdm import tqdm -from sp_validation.cosmology import get_cosmo, get_theo_c_ell from sp_validation.rho_tau import SquareRootScale mscale.register_scale(SquareRootScale) diff --git a/src/sp_validation/__init__.py b/src/sp_validation/__init__.py index 3d1ad5f1..d2b502ce 100644 --- a/src/sp_validation/__init__.py +++ b/src/sp_validation/__init__.py @@ -6,7 +6,6 @@ "calibration", "catalog", "catalog_builders", - "cosmology", "format", "galaxy", "glass_mock", diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py index f0da212a..13f4c2e9 100644 --- a/src/sp_validation/b_modes.py +++ b/src/sp_validation/b_modes.py @@ -12,11 +12,10 @@ import seaborn as sns import tqdm import treecorr +from cs_util.cosmo import get_theo_xi from mpl_toolkits.axes_grid1 import make_axes_locatable from scipy import sparse, stats -from .cosmology import get_theo_xi - _EB_KEYS = ("xip_E", "xim_E", "xip_B", "xim_B", "xip_amb", "xim_amb") diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py index 70f41a70..48251bbb 100644 --- a/src/sp_validation/cosmo_val/core.py +++ b/src/sp_validation/cosmo_val/core.py @@ -7,13 +7,13 @@ import colorama import numpy as np import yaml +from cs_util.cosmo import get_cosmo from shear_psf_leakage import run_object, run_scale from ..b_modes import ( _get_pte_from_scale_cut, find_conservative_scale_cut_key, ) -from ..cosmology import get_cosmo from ..statistics import chi2_and_pte from .catalog_characterization import CatalogCharacterizationMixin from .cosebis import CosebisMixin diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index e4733eea..514049a6 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -15,8 +15,8 @@ import numpy as np import pymaster as nmt from astropy.io import fits +from cs_util.cosmo import get_theo_c_ell -from ..cosmology import get_theo_c_ell from ..pseudo_cl import ( apply_random_rotation, get_n_gal_map, diff --git a/src/sp_validation/cosmology.py b/src/sp_validation/cosmology.py deleted file mode 100644 index 2964aa08..00000000 --- a/src/sp_validation/cosmology.py +++ /dev/null @@ -1,510 +0,0 @@ -"""COSMOLOGY. - -:Name: cosmology.py - -:Description: This file contains methods for science - validation of a weak-lensing shape catalogue. - Depends on a cosmological model. - -:Author: Axel Guinot, Martin Kilbinger, Cail Daley - -""" - -import camb -import numpy as np -import pyccl as ccl -from astropy.cosmology import Planck18 - -# ============================================================================= -# Fiducial Cosmology: astropy Planck18 -# ============================================================================= -# Source: Planck 2018 Paper VI, Table 2 (TT,TE,EE+lowE+lensing+BAO) -# Reference: Planck Collaboration 2020, A&A, 641, A6 -# -# Note on sigma_8 / A_s consistency: -# CAMB with A_s=2.1e-9 and m_nu=0.06 eV derives sigma_8 ~ 0.806, not 0.8102. -# This ~0.5% difference arises from Planck's MCMC marginalization details. -# Policy: Use sigma_8=0.8102 for codes taking sigma_8 directly (CosmoCov, CCL); -# use A_s=2.1e-9 for CAMB-based predictions. -# ============================================================================= -PLANCK18 = { - "Omega_m": Planck18.Om0, # 0.30966 - "Omega_b": Planck18.Ob0, # 0.04897 - "h": Planck18.h, # 0.6766 - "n_s": Planck18.meta["n"], # 0.9665 - "sigma_8": Planck18.meta["sigma8"], # 0.8102 - "A_s": 2.1e-9, # ln(10^10 A_s) = 3.047 - "m_nu": 0.06, # eV, sum of neutrino masses - "w0": -1.0, - "wa": 0.0, -} - - -def _ccl_to_camb(cosmo): - """Convert CCL cosmology object to CAMB parameter format. - - Parameters - ---------- - cosmo : ccl.Cosmology - CCL cosmology object - - Returns - ------- - dict - CAMB parameters dictionary with As properly set - """ - - h = cosmo["h"] - camb_params = { - "H0": h * 100, - "ombh2": cosmo["Omega_b"] * h**2, - "omch2": cosmo["Omega_c"] * h**2, - "ns": cosmo["n_s"], - } - - # Handle normalization: prefer As, but convert sigma8 to As if needed - As_val = cosmo.__getitem__("A_s") - sigma8_val = cosmo.__getitem__("sigma8") - - if As_val is not None and not np.isnan(As_val): - # Use As directly - camb_params["As"] = As_val - elif sigma8_val is not None: - # Convert sigma8 to As using iterative CAMB calculation - # see https://cosmocoffee.info/viewtopic.php?t=475 - As_fiducial = 2e-9 # Standard fiducial value - - # Step 1: Calculate current sigma8 with fiducial As - temp_params = camb_params.copy() - temp_params["As"] = As_fiducial - - pars = camb.set_params(**temp_params) - pars.set_matter_power(redshifts=[0.0], kmax=2.0) - results = camb.get_results(pars) - sigma8_current = results.get_sigma8_0() - - # Step 2: Scale As to match target sigma8 - # As scales as sigma8^2 - As_scaled = As_fiducial * (sigma8_val / sigma8_current) ** 2 - camb_params["As"] = As_scaled - - # Step 3: Verify the result - temp_params["As"] = As_scaled - pars = camb.set_params(**temp_params) - pars.set_matter_power(redshifts=[0.0], kmax=2.0) - results = camb.get_results(pars) - sigma8_final = results.get_sigma8_0() - - # Check accuracy (warn if >1% difference) - relative_error = abs(sigma8_final - sigma8_val) / sigma8_val - if relative_error > 0.01: - print( - f"Warning: CAMB sigma8 conversion accuracy: target={sigma8_val:.4f}, " - f"achieved={sigma8_final:.4f}, error={relative_error:.1%}" - ) - else: - # No normalization specified, use CAMB default - pass - - # Add dark energy parameters if they exist - for camb_key, cosmo_key in [("w", "w0"), ("wa", "wa")]: - if hasattr(cosmo._params, cosmo_key): - camb_params[camb_key] = cosmo[cosmo_key] - - return camb_params - - -def _camb_to_ccl(camb_params): - """Convert CAMB parameter format to CCL parameter dictionary. - - Parameters - ---------- - camb_params : dict - CAMB parameters dictionary with required keys: - H0, ombh2, omch2, ns, and either As or sigma8 - - Returns - ------- - dict - CCL parameters dictionary - """ - h = camb_params["H0"] / 100.0 - ccl_params = { - "Omega_c": camb_params["omch2"] / h**2, - "Omega_b": camb_params["ombh2"] / h**2, - "h": h, - "n_s": camb_params["ns"], - **{ - k: camb_params[v] - for k, v in [("w0", "w"), ("wa", "wa")] - if v in camb_params - }, - } - - # CCL accepts either A_s or sigma8 directly - if "As" in camb_params: - ccl_params["A_s"] = camb_params["As"] - elif "sigma8" in camb_params: - ccl_params["sigma8"] = camb_params["sigma8"] - else: - raise ValueError("Must provide either 'As' or 'sigma8' in camb_params") - - return ccl_params - - -def _cosmocov_to_ccl(cosmocov_params): - """Convert CosmoCov parameter format to CCL parameter dictionary. - - Parameters - ---------- - cosmocov_params : dict - CosmoCov parameters dictionary with required keys: - Omega_m, omb, h0, sigma_8, n_spec - - Returns - ------- - dict - CCL parameters dictionary - """ - required_params = ["Omega_m", "omb", "h0", "n_spec", "sigma_8"] - missing_params = [p for p in required_params if p not in cosmocov_params] - if missing_params: - raise KeyError(f"Missing required cosmological parameters: {missing_params}") - - ccl_params = { - "Omega_c": cosmocov_params["Omega_m"] - cosmocov_params["omb"], - "Omega_b": cosmocov_params["omb"], - "h": cosmocov_params["h0"], - "sigma8": cosmocov_params["sigma_8"], - "n_s": cosmocov_params["n_spec"], - **{k: cosmocov_params[k] for k in ["w0", "wa"] if k in cosmocov_params}, - } - - return ccl_params - - -def get_cosmo( - Omega_b=None, - Omega_m=None, - h=None, - sig8=None, - ns=None, - w0=None, - wa=None, - mnu=None, - transfer_function="boltzmann_camb", - matter_power_spectrum="halofit", - cosmocov_params=None, - camb_params=None, - extra_params=None, -): - """Get CCL cosmology object with user-specified parameters. - - Defaults to astropy Planck18 cosmology (Table 2: TT,TE,EE+lowE+lensing+BAO). - Can also use CosmoCov or CAMB parameter formats. - - Parameters - ---------- - Omega_m : float, default=None - Matter density parameter (defaults to Planck18: 0.30966) - Omega_b : float, default=None - Baryon density parameter (defaults to Planck18: 0.04897) - h : float, default=None - Reduced Hubble constant (defaults to Planck18: 0.6766) - sig8 : float, default=None - RMS matter fluctuation amplitude at 8 Mpc/h (defaults to Planck18: 0.8102) - ns : float, default=None - Scalar spectral index (defaults to Planck18: 0.9665) - w0 : float, default=None - Dark energy equation of state parameter (defaults to -1.0) - wa : float, default=None - Dark energy equation of state evolution parameter (defaults to 0.0) - mnu : float, default=None - Total neutrino mass in eV (defaults to 0.06 eV) - transfer_function : str, default='boltzmann_camb' - Transfer function to use - matter_power_spectrum : str, default='halofit' - Matter power spectrum to use - cosmocov_params : dict, optional - Parameters in CosmoCov format (Omega_m, omb, h0, sigma_8, n_spec) - If provided, entries override above parameters. Mutually exclusive with - camb_params. - camb_params : dict, optional - Parameters in CAMB format (H0, ombh2, omch2, ns, sigma8) - If provided, entries override above parameters. Mutually exclusive with - cosmocov_params. - extra_params : dict, optional - Additional parameters to pass to CCL (e.g., for CAMB non-linear settings) - - Returns - ------- - Cosmology - pyccl cosmology object - """ - # Check for parameter format conflicts - if cosmocov_params is not None and camb_params is not None: - raise ValueError( - "Cannot provide both cosmocov_params and camb_params. Choose one format." - ) - - # Convert parameters to CCL format - if cosmocov_params is not None: - print("Using CosmoCov parameters to create CCL cosmology.") - ccl_params = _cosmocov_to_ccl(cosmocov_params) - elif camb_params is not None: - print("Using CAMB parameters to create CCL cosmology.") - ccl_params = _camb_to_ccl(camb_params) - else: - ccl_params = {} - - # Planck 2018 defaults from astropy (see PLANCK18 dict at module level) - planck_defaults = { - "Omega_m": PLANCK18["Omega_m"], - "Omega_b": PLANCK18["Omega_b"], - "h": PLANCK18["h"], - "sig8": PLANCK18["sigma_8"], - "ns": PLANCK18["n_s"], - "w0": PLANCK18["w0"], - "wa": PLANCK18["wa"], - "mnu": PLANCK18["m_nu"], - } - - mnu = ccl_params.get("mnu", mnu or planck_defaults["mnu"]) - h = ccl_params.get("h", h or planck_defaults["h"]) - pars = camb.CAMBparams() - - pars.set_cosmology( - mnu=mnu, - H0=h * 100, - ) - - Omega_nu = pars.omeganu - - combined_params = { - "Omega_c": ccl_params.get( - "Omega_c", - (Omega_m or planck_defaults["Omega_m"]) - - (Omega_b or planck_defaults["Omega_b"]) - - Omega_nu, - ), - "Omega_b": ccl_params.get("Omega_b", Omega_b or planck_defaults["Omega_b"]), - "h": h, - "sigma8": ccl_params.get("sigma8", sig8 or planck_defaults["sig8"]), - "n_s": ccl_params.get("n_s", ns or planck_defaults["ns"]), - "w0": ccl_params.get("w0", w0 or planck_defaults["w0"]), - "wa": ccl_params.get("wa", wa or planck_defaults["wa"]), - "m_nu": mnu, - "extra_parameters": extra_params, - } - - return ccl.Cosmology( - **combined_params, - transfer_function=transfer_function, - matter_power_spectrum=matter_power_spectrum, - ) - - -def get_theo_c_ell( - ell, - z, - nz, - backend="ccl", - cosmo=None, - Omega_b=None, - Omega_m=None, - h=None, - sig8=None, - ns=None, - w0=None, - wa=None, -): - """Calculate theoretical angular power spectrum C_ell for weak lensing. - - Parameters - ---------- - ell : array - Multipole moments (e.g., np.arange(2, 2000)) - z : array - Redshifts for n(z) distribution - nz : array - n(z) redshift distribution - backend : str, default="ccl" - Backend to use: "ccl" or "camb" - cosmo : ccl.Cosmology, optional - CCL cosmology object. If None, will create using individual parameters. - Omega_b : float, optional - Baryon density parameter (defaults to Planck 2018) - Omega_m : float, optional - Matter density parameter (defaults to Planck 2018) - h : float, optional - Reduced Hubble constant (defaults to Planck 2018) - sig8 : float, optional - RMS matter fluctuation amplitude at 8 Mpc/h (defaults to Planck 2018) - ns : float, optional - Scalar spectral index (defaults to Planck 2018) - w0 : float, optional - Dark energy equation of state parameter (defaults to -1.0) - wa : float, optional - Dark energy equation of state evolution parameter (defaults to 0.0) - - Returns - ------- - cl : array - Angular power spectrum - """ - if cosmo is None: - cosmo = get_cosmo( - Omega_b=Omega_b, - Omega_m=Omega_m, - h=h, - sig8=sig8, - ns=ns, - w0=w0, - wa=wa, - ) - - if backend == "ccl": - # Create lensing tracer - lens = ccl.WeakLensingTracer(cosmo, dndz=(z, nz)) - - # Calculate power spectrum - cl = ccl.angular_cl(cosmo, lens, lens, ell) - - elif backend == "camb": - # Convert CCL cosmology to CAMB parameters - import camb - - camb_kwargs = _ccl_to_camb(cosmo) - - # Set up CAMB parameters - pars = camb.set_params( - **camb_kwargs, - WantTransfer=True, - NonLinear=camb.model.NonLinear_both, - ) - - # Adjust for neutrino contribution - if "mnu" in camb_kwargs and camb_kwargs["mnu"] > 0: - omch2_adj = camb_kwargs["omch2"] - pars.omeganu * (pars.H0 / 100) ** 2 - pars.set_cosmology(omch2=omch2_adj) - - # Set up lensing source window - pars.min_l = ell.min() - pars.set_for_lmax(ell.max()) - pars.SourceWindows = [ - camb.sources.SplinedSourceWindow(z=z, W=nz, source_type="lensing") - ] - - # Calculate power spectrum - results = camb.get_results(pars) - theory_cls = results.get_source_cls_dict(lmax=ell.max(), raw_cl=True) - cl_full = theory_cls["W1xW1"] - - # Interpolate to match input ell array - # CAMB returns C_ell for ell = 0, 1, 2, ..., lmax - ell_camb = np.arange(len(cl_full)) - cl = np.interp(ell, ell_camb, cl_full) - - else: - raise ValueError(f"Unknown backend: {backend}. Must be 'ccl' or 'camb'") - - return cl - - -def c_ell_to_xi(cosmo, theta, ell, cl): - """Convert angular power spectrum to correlation functions using CCL. - - Parameters - ---------- - cosmo : ccl.Cosmology - CCL cosmology object (used for correlation function calculation) - theta : array - Angular separations in arcminutes - ell : array - Multipole moments - cl : array - Angular power spectrum - - Returns - ------- - xip, xim : arrays - xi+ and xi- correlation functions - """ - theta_deg = theta / 60.0 # arcmin to degrees - - xip = ccl.correlation( - cosmo, ell=ell, C_ell=cl, theta=theta_deg, type="GG+", method="Bessel" - ) - xim = ccl.correlation( - cosmo, ell=ell, C_ell=cl, theta=theta_deg, type="GG-", method="Bessel" - ) - - return xip, xim - - -def get_theo_xi( - theta, - z, - nz, - Omega_m=None, - h=None, - Omega_b=None, - sig8=None, - ns=None, - ell_min=10, - ell_max=20000, - n_ell=500, - backend="ccl", - cosmo=None, - **cosmo_kwargs, -): - """Calculate theoretical xi+/xi- using individual parameters. - - Parameters - ---------- - theta : array - Angular separations in arcminutes - z : array - Redshift array - nz : array - n(z) redshift distribution - Omega_m : float, default=None - Matter density parameter (defaults to Planck 2018) - h : float, default=None - Reduced Hubble constant (defaults to Planck 2018) - Omega_b : float, default=None - Baryon density parameter (defaults to Planck 2018) - sig8 : float, default=None - RMS matter fluctuation amplitude at 8 Mpc/h (defaults to Planck 2018) - ns : float, default=None - Scalar spectral index (defaults to Planck 2018) - ell_min : int, default=0 - Minimum ell for power spectrum calculation - ell_max : int, default=20000 - Maximum ell for power spectrum calculation - n_ell : int, default=500 - Number of ell bins - backend : str, default="ccl" - Backend to use: "ccl" or "camb" - **cosmo_kwargs - Additional arguments passed to backend - - Returns - ------- - xip, xim : arrays - Theoretical xi+ and xi- correlation functions - """ - # Create ell array for C_ell calculation - ell = np.geomspace(ell_min, ell_max, n_ell) - - # Use provided cosmology or create from parameters - if cosmo is None: - cosmo = get_cosmo( - Omega_m=Omega_m, Omega_b=Omega_b, h=h, sig8=sig8, ns=ns, **cosmo_kwargs - ) - - # Calculate C_ell - cl = get_theo_c_ell(ell, z, nz, backend=backend, cosmo=cosmo) - - # Convert to xi - return c_ell_to_xi(cosmo, theta, ell, cl) diff --git a/src/sp_validation/tests/test_cosmology.py b/src/sp_validation/tests/test_cosmology.py deleted file mode 100644 index 348f1985..00000000 --- a/src/sp_validation/tests/test_cosmology.py +++ /dev/null @@ -1,519 +0,0 @@ -"""UNIT TESTS FOR COSMOLOGY SUBPACKAGE. - -This module contains unit tests for the module package -sp_validation.cosmology. - -:Author: Claude Code Assistant - -""" - -import numpy as np -import numpy.testing as npt -import pyccl as ccl -import pytest - -from sp_validation import cosmology - -# Test performance markers -pytestmark = pytest.mark.fast # Mark all tests as fast by default - -# Standard test cosmology used across all tests -TEST_COSMOLOGY = { - "Omega_m": 0.3, - "Omega_b": 0.05, - "h": 0.7, - "sig8": 0.8, - "ns": 0.96, - "w0": -1.0, - "wa": 0.0, -} - -# Derived CAMB format from TEST_COSMOLOGY -TEST_CAMB = { - "H0": TEST_COSMOLOGY["h"] * 100, # 70.0 - "ombh2": TEST_COSMOLOGY["Omega_b"] * TEST_COSMOLOGY["h"] ** 2, - "omch2": ( - (TEST_COSMOLOGY["Omega_m"] - TEST_COSMOLOGY["Omega_b"]) - * TEST_COSMOLOGY["h"] ** 2 - ), - "ns": TEST_COSMOLOGY["ns"], # 0.96 - "sigma8": TEST_COSMOLOGY["sig8"], # 0.8 - "w": TEST_COSMOLOGY["w0"], # -1.0 - "wa": TEST_COSMOLOGY["wa"], # 0.0 -} - -# Derived CosmoCov format from TEST_COSMOLOGY -TEST_COSMOCOV = { - "Omega_m": TEST_COSMOLOGY["Omega_m"], # 0.3 - "omb": TEST_COSMOLOGY["Omega_b"], # 0.05 - "h0": TEST_COSMOLOGY["h"], # 0.7 - "sigma_8": TEST_COSMOLOGY["sig8"], # 0.8 - "n_spec": TEST_COSMOLOGY["ns"], # 0.96 - "w0": TEST_COSMOLOGY["w0"], # -1.0 - "wa": TEST_COSMOLOGY["wa"], # 0.0 -} - -# Fast test parameters (reduced precision for speed) -FAST_ELL_MAX = 2000 -FAST_N_ELL = 100 -FAST_Z_POINTS = 20 - -# Realistic test parameters (full precision for one integration test) -REALISTIC_ELL_MAX = 10_000 -REALISTIC_N_ELL = 500 - - -class TestCosmologyHelpers: - """Test helper functions for parameter conversion.""" - - @pytest.fixture(autouse=True) - def setup_method(self): - """Set up test fixtures.""" - # Use standard test cosmology definitions - self.camb_params = { - k: v - for k, v in TEST_CAMB.items() - if k in ["H0", "ombh2", "omch2", "ns", "sigma8"] - } - self.camb_params_full = TEST_CAMB.copy() - self.cosmocov_params = { - k: v - for k, v in TEST_COSMOCOV.items() - if k in ["Omega_m", "omb", "h0", "sigma_8", "n_spec"] - } - self.cosmocov_params_full = TEST_COSMOCOV.copy() - - def test_ccl_to_camb_basic(self): - """Test CCL to CAMB parameter conversion.""" - # Create a test cosmology using TEST_COSMOLOGY values - omega_c = TEST_COSMOLOGY["Omega_m"] - TEST_COSMOLOGY["Omega_b"] - cosmo = ccl.Cosmology( - Omega_c=omega_c, - Omega_b=TEST_COSMOLOGY["Omega_b"], - h=TEST_COSMOLOGY["h"], - sigma8=TEST_COSMOLOGY["sig8"], - n_s=TEST_COSMOLOGY["ns"], - ) - - # Convert to CAMB parameters - camb_params = cosmology._ccl_to_camb(cosmo) - - # Check basic parameters against expected TEST_CAMB values - npt.assert_almost_equal(camb_params["H0"], TEST_CAMB["H0"], decimal=2) - npt.assert_almost_equal(camb_params["ombh2"], TEST_CAMB["ombh2"], decimal=3) - npt.assert_almost_equal(camb_params["omch2"], TEST_CAMB["omch2"], decimal=3) - npt.assert_almost_equal(camb_params["ns"], TEST_CAMB["ns"], decimal=4) - # Check that As is present and positive (sigma8 converted to As) - assert "As" in camb_params - assert camb_params["As"] > 0 - - def test_camb_to_ccl_basic(self): - """Test CAMB to CCL parameter conversion.""" - ccl_params = cosmology._camb_to_ccl(self.camb_params) - - # Check conversion against expected TEST_COSMOLOGY values - assert isinstance(ccl_params, dict) - npt.assert_almost_equal( - ccl_params["Omega_b"], TEST_COSMOLOGY["Omega_b"], decimal=3 - ) - npt.assert_almost_equal( - ccl_params["Omega_c"], - TEST_COSMOLOGY["Omega_m"] - TEST_COSMOLOGY["Omega_b"], - decimal=3, - ) - npt.assert_almost_equal(ccl_params["h"], TEST_COSMOLOGY["h"], decimal=4) - npt.assert_almost_equal(ccl_params["n_s"], TEST_COSMOLOGY["ns"], decimal=4) - npt.assert_almost_equal(ccl_params["sigma8"], TEST_COSMOLOGY["sig8"], decimal=4) - - def test_camb_to_ccl_with_optional_params(self): - """Test CAMB to CCL conversion with optional parameters.""" - ccl_params = cosmology._camb_to_ccl(self.camb_params_full) - - # Check optional parameters against TEST_CAMB values - npt.assert_almost_equal(ccl_params["w0"], TEST_CAMB["w"], decimal=4) - npt.assert_almost_equal(ccl_params["wa"], TEST_CAMB["wa"], decimal=4) - - def test_camb_to_ccl_missing_normalization(self): - """Test error when neither As nor sigma8 provided.""" - bad_params = { - "H0": TEST_CAMB["H0"], - "ombh2": TEST_CAMB["ombh2"], - "omch2": TEST_CAMB["omch2"], - "ns": TEST_CAMB["ns"], - } - - with pytest.raises(ValueError): - cosmology._camb_to_ccl(bad_params) - - def test_cosmocov_to_ccl_basic(self): - """Test CosmoCov to CCL parameter conversion.""" - ccl_params = cosmology._cosmocov_to_ccl(self.cosmocov_params) - - # Check conversion against TEST_COSMOCOV values - assert isinstance(ccl_params, dict) - expected_omega_c = TEST_COSMOCOV["Omega_m"] - TEST_COSMOCOV["omb"] - npt.assert_almost_equal(ccl_params["Omega_c"], expected_omega_c, decimal=4) - npt.assert_almost_equal(ccl_params["Omega_b"], TEST_COSMOCOV["omb"], decimal=4) - npt.assert_almost_equal(ccl_params["h"], TEST_COSMOCOV["h0"], decimal=4) - npt.assert_almost_equal( - ccl_params["sigma8"], TEST_COSMOCOV["sigma_8"], decimal=4 - ) - npt.assert_almost_equal(ccl_params["n_s"], TEST_COSMOCOV["n_spec"], decimal=4) - - def test_cosmocov_to_ccl_with_optional_params(self): - """Test CosmoCov to CCL conversion with optional parameters.""" - ccl_params = cosmology._cosmocov_to_ccl(self.cosmocov_params_full) - - # Check optional parameters against TEST_COSMOCOV values - npt.assert_almost_equal(ccl_params["w0"], TEST_COSMOCOV["w0"], decimal=4) - npt.assert_almost_equal(ccl_params["wa"], TEST_COSMOCOV["wa"], decimal=4) - - def test_cosmocov_to_ccl_missing_params(self): - """Test error when required CosmoCov parameters are missing.""" - incomplete_params = { - "Omega_m": TEST_COSMOCOV["Omega_m"], - "omb": TEST_COSMOCOV["omb"], - # Missing h0, sigma_8, n_spec - } - - with pytest.raises(KeyError): - cosmology._cosmocov_to_ccl(incomplete_params) - - -class TestGetCosmo: - """Test get_cosmo function.""" - - def test_default_cosmology(self): - """Test that default parameters are Planck18 (not our test values).""" - cosmo = cosmology.get_cosmo() - - # Check that we get a CCL cosmology object - assert isinstance(cosmo, ccl.Cosmology) - - # Check that defaults are NOT our test values (should be Planck18) - assert cosmo["h"] != TEST_COSMOLOGY["h"] - assert cosmo["Omega_b"] != TEST_COSMOLOGY["Omega_b"] - - def test_individual_parameters(self): - """Test overriding individual parameters.""" - cosmo = cosmology.get_cosmo(h=TEST_COSMOLOGY["h"], sig8=TEST_COSMOLOGY["sig8"]) - - # Check custom parameters - assert cosmo["h"] == pytest.approx(TEST_COSMOLOGY["h"], abs=1e-4) - assert cosmo["sigma8"] == pytest.approx(TEST_COSMOLOGY["sig8"], abs=1e-4) - - # Check other parameters use defaults (not our test values) - assert cosmo["Omega_b"] != TEST_COSMOLOGY["Omega_b"] - - def test_cosmocov_params(self): - """Test using CosmoCov parameter format.""" - cosmo = cosmology.get_cosmo(cosmocov_params=TEST_COSMOCOV) - - # Check parameters against TEST_COSMOCOV constants - assert cosmo["Omega_b"] == pytest.approx(TEST_COSMOCOV["omb"], abs=1e-4) - assert cosmo["h"] == pytest.approx(TEST_COSMOCOV["h0"], abs=1e-4) - assert cosmo["sigma8"] == pytest.approx(TEST_COSMOCOV["sigma_8"], abs=1e-4) - assert cosmo["n_s"] == pytest.approx(TEST_COSMOCOV["n_spec"], abs=1e-4) - - def test_camb_params(self): - """Test using CAMB parameter format.""" - cosmo = cosmology.get_cosmo(camb_params=TEST_CAMB) - - # Check parameters against TEST_CAMB constants - expected_h = TEST_CAMB["H0"] / 100.0 - expected_omega_b = TEST_CAMB["ombh2"] / expected_h**2 - - assert cosmo["h"] == pytest.approx(expected_h, abs=1e-4) - assert cosmo["Omega_b"] == pytest.approx(expected_omega_b, abs=1e-4) - assert cosmo["sigma8"] == pytest.approx(TEST_CAMB["sigma8"], abs=1e-4) - - def test_mutually_exclusive_params(self): - """Test that cosmocov_params and camb_params are mutually exclusive.""" - with pytest.raises(ValueError): - cosmology.get_cosmo(cosmocov_params=TEST_COSMOCOV, camb_params=TEST_CAMB) - - def test_transfer_function_options(self): - """Test different transfer function options.""" - cosmo1 = cosmology.get_cosmo(transfer_function="boltzmann_camb") - cosmo2 = cosmology.get_cosmo(transfer_function="boltzmann_class") - - # Both should create valid cosmologies - assert isinstance(cosmo1, ccl.Cosmology) - assert isinstance(cosmo2, ccl.Cosmology) - - -# Shared test fixtures to avoid recomputation -@pytest.fixture -def fast_redshift_data(): - """Fast redshift distribution for most tests.""" - z = np.linspace(0.1, 2.0, FAST_Z_POINTS) - nz = np.exp(-(((z - 0.8) / 0.2) ** 2)) # Gaussian around z=0.8 - nz /= np.trapezoid(nz, z) # Normalize - return z, nz - - -@pytest.fixture -def fast_ell_array(): - """Fast ell array for most tests.""" - return np.array([10, 100, 1000]) # Reduced from [10, 50, 100, 500, 1000] - - -@pytest.fixture -def realistic_redshift_data(): - """Realistic redshift distribution for integration tests.""" - z = np.linspace(0.1, 2.0, 50) - nz = np.exp(-(((z - 0.8) / 0.2) ** 2)) - nz /= np.trapezoid(nz, z) - return z, nz - - -class TestGetTheoCell: - """Test get_theo_c_ell function.""" - - def test_ccl_backend_basic(self, fast_ell_array, fast_redshift_data): - """Test C_ell calculation with CCL backend - basic functionality.""" - z, nz = fast_redshift_data - cosmo = cosmology.get_cosmo() - - cl = cosmology.get_theo_c_ell(fast_ell_array, z, nz, backend="ccl", cosmo=cosmo) - - # Check output shape and positivity - assert len(cl) == len(fast_ell_array) - assert np.all(cl > 0) - - def test_ccl_backend_with_params(self, fast_ell_array, fast_redshift_data): - """Test C_ell calculation with CCL backend using individual parameters.""" - z, nz = fast_redshift_data - cl = cosmology.get_theo_c_ell( - fast_ell_array, z, nz, backend="ccl", **TEST_COSMOLOGY - ) - - # Check output shape and positivity - assert len(cl) == len(fast_ell_array) - assert np.all(cl > 0) - - def test_ccl_backend_default_params(self, fast_ell_array, fast_redshift_data): - """Test C_ell calculation with CCL backend using default parameters.""" - z, nz = fast_redshift_data - cl = cosmology.get_theo_c_ell(fast_ell_array, z, nz, backend="ccl") - - # Check output shape and positivity - assert len(cl) == len(fast_ell_array) - assert np.all(cl > 0) - - def test_camb_backend(self, fast_ell_array, fast_redshift_data): - """Test C_ell calculation with CAMB backend.""" - z, nz = fast_redshift_data - try: - cl = cosmology.get_theo_c_ell(fast_ell_array, z, nz, backend="camb") - - # Check output shape and positivity - assert len(cl) == len(fast_ell_array) - assert np.all(cl > 0) - except ImportError: - pytest.skip("CAMB not available") - - def test_invalid_backend(self, fast_ell_array, fast_redshift_data): - """Test error with invalid backend.""" - z, nz = fast_redshift_data - with pytest.raises(ValueError): - cosmology.get_theo_c_ell(fast_ell_array, z, nz, backend="invalid") - - def test_small_ell_array(self, fast_redshift_data): - """Test behavior with small ell array.""" - z, nz = fast_redshift_data - small_ell = np.array([10, 100]) - cl = cosmology.get_theo_c_ell(small_ell, z, nz, backend="ccl") - - assert len(cl) == 2 - assert np.all(cl > 0) - - @pytest.mark.slow - def test_backend_consistency_fast(self, fast_ell_array, fast_redshift_data): - """Test that CCL and CAMB backends give consistent results (fast version).""" - z, nz = fast_redshift_data - try: - # Calculate with both backends using fast parameters - cl_ccl = cosmology.get_theo_c_ell(fast_ell_array, z, nz, backend="ccl") - cl_camb = cosmology.get_theo_c_ell(fast_ell_array, z, nz, backend="camb") - - # Check that both have same shape - assert len(cl_ccl) == len(cl_camb) == len(fast_ell_array) - - # Check relative differences between cosmology codes - # CCL and CAMB should agree well since they compute the same physics - # Observed differences: ~4.8% max, ~2.8% average - rel_diff = np.abs(cl_camb - cl_ccl) / cl_ccl - max_rel_diff = rel_diff.max() - - assert max_rel_diff < 0.08, ( - f"CCL and CAMB backends differ by >8%: {max_rel_diff:.1%}" - ) - except ImportError: - pytest.skip("CAMB not available") - - -@pytest.fixture -def fast_theta_array(): - """Fast theta array for most tests.""" - return np.array([5.0, 20.0, 100.0]) # Reduced from logspace(0, 2, 10) - - -class TestGetTheoXi: - """Test get_theo_xi function.""" - - def test_backwards_compatibility(self, fast_theta_array, fast_redshift_data): - """Test backwards compatible parameter interface with fast parameters.""" - z, nz = fast_redshift_data - xip, xim = cosmology.get_theo_xi( - fast_theta_array, - z, - nz, - **TEST_COSMOLOGY, - ell_min=10, - ell_max=FAST_ELL_MAX, - n_ell=FAST_N_ELL, - backend="ccl", - ) - - # Check output shapes - assert len(xip) == len(fast_theta_array) - assert len(xim) == len(fast_theta_array) - - # Check that xi+ is positive for typical scales - assert np.any(xip > 0) - - def test_default_parameters(self, fast_theta_array, fast_redshift_data): - """Test with default cosmological parameters using fast settings.""" - z, nz = fast_redshift_data - xip, xim = cosmology.get_theo_xi( - fast_theta_array, z, nz, ell_min=10, ell_max=FAST_ELL_MAX, n_ell=FAST_N_ELL - ) - - # Check output shapes - assert len(xip) == len(fast_theta_array) - assert len(xim) == len(fast_theta_array) - - @pytest.mark.slow - def test_different_ell_ranges(self, fast_theta_array, fast_redshift_data): - """Test that different ell ranges give consistent results.""" - z, nz = fast_redshift_data - - # Coarse integration (lower ell_max) - xip1, xim1 = cosmology.get_theo_xi( - fast_theta_array, z, nz, ell_min=10, ell_max=1000, n_ell=50 - ) - - # Fine integration (higher ell_max) - xip2, xim2 = cosmology.get_theo_xi( - fast_theta_array, z, nz, ell_min=10, ell_max=FAST_ELL_MAX, n_ell=FAST_N_ELL - ) - - # Scale-dependent tolerance: - # Xi+ converges well at all scales (~4% max differences) - # Xi- has poor convergence at small scales due to oscillatory Hankel - # transforms and sensitivity to high-ell power (up to ~15% at θ<10 arcmin) - # but excellent convergence at large scales (~0.2% at θ>20 arcmin) - - small_scale_mask = fast_theta_array <= 10.0 # arcmin - large_scale_mask = fast_theta_array > 10.0 # arcmin - - # Xi+ has good convergence at all scales - npt.assert_allclose(xip1, xip2, rtol=0.06) # 6% tolerance for xi+ - - # Xi- has scale-dependent convergence - if np.any(small_scale_mask): - # Small scales (θ ≤ 10 arcmin): physics-limited precision - npt.assert_allclose( - xim1[small_scale_mask], xim2[small_scale_mask], rtol=0.18 - ) # 18% tolerance for xi- at small scales - - if np.any(large_scale_mask): - # Large scales (θ > 10 arcmin): excellent convergence - npt.assert_allclose( - xim1[large_scale_mask], xim2[large_scale_mask], rtol=0.05 - ) # 5% tolerance for xi- at large scales - - -class TestCosmologyIntegration: - """Integration tests for cosmology functions.""" - - def test_end_to_end_pipeline_fast(self, fast_redshift_data): - """Test complete pipeline from parameters to xi (fast version).""" - z, nz = fast_redshift_data - # Create cosmology using TEST_CAMB parameters - camb_subset = { - k: v - for k, v in TEST_CAMB.items() - if k in ["H0", "ombh2", "omch2", "ns", "sigma8"] - } - cosmo = cosmology.get_cosmo(camb_params=camb_subset) - - # Calculate C_ell with fast parameters - ell = np.logspace(1, 3, 15) # Reduced from 30 points - cl = cosmology.get_theo_c_ell(ell, z, nz, cosmo=cosmo) - - # Calculate xi - theta = np.array([5.0, 10.0, 20.0]) - xip, xim = cosmology.c_ell_to_xi(cosmo, theta, ell, cl) - - # Sanity checks - assert len(cl) == len(ell) - assert len(xip) == len(theta) - assert len(xim) == len(theta) - assert np.all(cl > 0) - - @pytest.mark.slow - def test_realistic_precision_pipeline(self, realistic_redshift_data): - """Test complete pipeline with realistic precision for production use.""" - z, nz = realistic_redshift_data - - # Create cosmology with realistic parameters - cosmo = cosmology.get_cosmo(**TEST_COSMOLOGY) - - # Calculate C_ell with realistic precision - ell = np.logspace(1, np.log10(REALISTIC_ELL_MAX), REALISTIC_N_ELL) - cl = cosmology.get_theo_c_ell(ell, z, nz, cosmo=cosmo) - - # Test xi calculation with realistic theta range - theta = np.geomspace(1, 300, 20) # 1 to 300 arcmin, 20 points - xip, xim = cosmology.get_theo_xi( - theta, - z, - nz, - ell_min=10, - ell_max=REALISTIC_ELL_MAX, - n_ell=REALISTIC_N_ELL, - **TEST_COSMOLOGY, - ) - - # Comprehensive checks for realistic test - assert len(cl) == len(ell) - assert len(xip) == len(theta) - assert len(xim) == len(theta) - assert np.all(cl > 0) - - # Check that xi+ is positive at small scales (expected for cosmic shear) - small_scale_mask = theta < 10 # arcmin - assert np.any(xip[small_scale_mask] > 0) - - # Check that xi+ amplitude is reasonable (order of magnitude check) - max_xip = np.max(np.abs(xip)) - assert 1e-6 < max_xip < 1e-2, f"xi+ amplitude unrealistic: {max_xip}" - - def test_parameter_format_consistency(self): - """Test that different parameter formats give consistent results.""" - # Individual parameters - cosmo1 = cosmology.get_cosmo(**TEST_COSMOLOGY) - - # CosmoCov format - cosmo2 = cosmology.get_cosmo(cosmocov_params=TEST_COSMOCOV) - - # CAMB format - cosmo3 = cosmology.get_cosmo(camb_params=TEST_CAMB) - - # All should give consistent results - npt.assert_allclose(cosmo1["Omega_b"], cosmo2["Omega_b"], rtol=1e-4) - npt.assert_allclose(cosmo1["Omega_b"], cosmo3["Omega_b"], rtol=1e-4) - npt.assert_allclose(cosmo1["h"], cosmo2["h"], rtol=1e-4) - npt.assert_allclose(cosmo1["h"], cosmo3["h"], rtol=1e-4) diff --git a/workflow/common.py b/workflow/common.py index d8b250fe..3df3413b 100644 --- a/workflow/common.py +++ b/workflow/common.py @@ -22,7 +22,7 @@ BLOCK_PAIRS = [("++", "1"), ("--", "2"), ("+-", "3")] # Fiducial cosmology: Planck 2018 (astropy Planck18, Table 2 + BAO) -# Source of truth: sp_validation.cosmology.PLANCK18 +# Source of truth: cs_util.cosmo.PLANCK18 # Regenerate with: snakemake results/cosmology/planck18.json # Resolved relative to the run directory at configure() time. COSMOLOGY_PARAMS = "results/cosmology/planck18.json" diff --git a/workflow/rules/covariance.smk b/workflow/rules/covariance.smk index d52608e7..49e40201 100644 --- a/workflow/rules/covariance.smk +++ b/workflow/rules/covariance.smk @@ -37,14 +37,14 @@ rule cosmology_params: """Generate cosmology parameters JSON from sp_validation. This decouples snakemake parse-time from sp_validation import. - Source of truth remains sp_validation.cosmology.PLANCK18. + Source of truth remains cs_util.cosmo.PLANCK18. """ output: COSMOLOGY_PARAMS shell: """ python -c " -from sp_validation.cosmology import PLANCK18 +from cs_util.cosmo import PLANCK18 import json from pathlib import Path diff --git a/workflow/scripts/generate_cosmocov_ini.py b/workflow/scripts/generate_cosmocov_ini.py index 631315c2..11e785cf 100644 --- a/workflow/scripts/generate_cosmocov_ini.py +++ b/workflow/scripts/generate_cosmocov_ini.py @@ -2,7 +2,7 @@ CLI refactor of the former ``rule covariance_ini`` heredoc. Cosmology is read from the frozen ``planck18.json`` snapshot (the cosmology_snapshot lc output — -source of truth is sp_validation.cosmology.PLANCK18); survey (area, n_eff, +source of truth is cs_util.cosmo.PLANCK18); survey (area, n_eff, sigma_e) from the catalog config's per-version ``cov_th``; n(z) via the same path convention as workflow/common.build_redshift_path; the footprint mask power spectrum is passed explicitly (empty string for the unmasked variant).