From fa0b8a08c744d404642959e09517f6370ec0369e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 3 Jul 2026 22:01:44 +0200 Subject: [PATCH] feat(cosmo): consolidate sp_validation cosmology theory-curve machinery into cs_util.cosmo Move the cosmology theory-curve machinery from sp_validation.cosmology into cs_util.cosmo so cs_util is the single home for the collaboration's cosmology primitives (per a conversation with M. Kilbinger). Additive: no existing cs_util.cosmo symbol changes, so shear_psf_leakage (get_cosmo_default, xipm_theo) is unaffected. Adds: - PLANCK18 fiducial-cosmology dict - get_cosmo: flexible Planck18-default cosmology builder (individual, CCL, CAMB, or CosmoCov parameter formats) - get_theo_c_ell, c_ell_to_xi, get_theo_xi: angular power spectra and shear correlation functions, CCL and CAMB backends - _ccl_to_camb / _camb_to_ccl / _cosmocov_to_ccl parameter converters - tests/test_cosmology.py, ported from sp_validation Caps scipy < 1.18 (mirroring sp_validation): scipy 1.18's FITPACK-to-C port changed a RectBivariateSpline return shape that breaks camb's BBN Y_He predictor in set_cosmology (and pyccl's Boltzmann spline path likewise). With the cap the full cosmo suite passes for real on py3.10 AND py3.12 (numpy 2.5): 31/31. This also removes the numpy>=2.5 xfail on test_xipm_theo -- a stopgap for the same root cause (#70/#71) -- which now passes. Bumps version 0.2.1 -> 0.2.2. No new dependency (camb, pyccl already required). Part of #75; the sp_validation delete + repoint follows in a companion PR. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_013NCc4BHtZ4Lx6n9zuUwEeP --- cs_util/cosmo.py | 523 +++++++++++++++++++++++++++++++- cs_util/tests/test_cosmo.py | 11 - cs_util/tests/test_cosmology.py | 519 +++++++++++++++++++++++++++++++ pyproject.toml | 17 +- 4 files changed, 1053 insertions(+), 17 deletions(-) create mode 100644 cs_util/tests/test_cosmology.py diff --git a/cs_util/cosmo.py b/cs_util/cosmo.py index ebe874e..421a5a8 100644 --- a/cs_util/cosmo.py +++ b/cs_util/cosmo.py @@ -2,17 +2,28 @@ :Name: cosmo.py -:Description: This file contains methods for cosmological - quantities. - -:Author: Martin Kilbinger +:Description: This file contains methods for cosmological quantities: + lensing critical surface mass density (``sigma_crit`` family) + and shear two-point theory predictions. The theory-curve + machinery — flexible cosmology construction (``get_cosmo``), + angular power spectra and correlation functions + (``get_theo_c_ell``, ``c_ell_to_xi``, ``get_theo_xi``), and + CCL/CAMB/CosmoCov parameter conversion — was consolidated here + from ``sp_validation.cosmology`` so that ``cs_util`` is the + single home for the collaboration's cosmology primitives. + +:Authors: Martin Kilbinger + Axel Guinot + Cail Daley """ +import camb import numpy as np from astropy import constants from astropy import units +from astropy.cosmology import Planck18 import pyccl as ccl @@ -307,3 +318,507 @@ def xipm_theo( ) return xipm["GG+"], xipm["GG-"] + + +# ============================================================================ +# Theory-curve machinery consolidated from sp_validation.cosmology +# (flexible cosmology construction, C_ell / xi_pm predictions, and +# CCL/CAMB/CosmoCov parameter conversion). cs_util is the single home. +# ============================================================================ +# ============================================================================= +# 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 + # Flat-universe dark-energy density (1 - Omega_m = 0.69034, matching the + # paper-era planck18.json and CosmoCov's flat-LCDM convention; astropy's + # Ode0 = 0.68885 differs by the radiation/neutrino contributions). + "Omega_v": 1.0 - Planck18.Om0, # 0.69034 + "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/cs_util/tests/test_cosmo.py b/cs_util/tests/test_cosmo.py index bf058c6..7380ba2 100644 --- a/cs_util/tests/test_cosmo.py +++ b/cs_util/tests/test_cosmo.py @@ -279,17 +279,6 @@ def test_get_cosmo_default(self): npt.assert_equal(pickle.dumps(cos_def), pickle.dumps(self._cos_def)) - @pytest.mark.xfail( - np.lib.NumpyVersion(np.__version__) >= "2.5.0", - reason=( - "pyccl 3.3.4 calls float() on a non-0-d array in boltzmann.py, " - "which numpy>=2.5 turns into a hard TypeError; not a cs_util bug. " - "Stopgap, not the real fix. Tracked in CosmoStat/cs_util#71 " - "(revisit the pyccl/numpy pin and where cosmology code should " - "live)." - ), - strict=False, - ) def test_xipm_theo(self): """Test ``cs_util.xipm_theo``: xi+/- theory against stored values.""" xip, xim = cosmo.xipm_theo( diff --git a/cs_util/tests/test_cosmology.py b/cs_util/tests/test_cosmology.py new file mode 100644 index 0000000..019f055 --- /dev/null +++ b/cs_util/tests/test_cosmology.py @@ -0,0 +1,519 @@ +"""UNIT TESTS FOR COSMOLOGY SUBPACKAGE. + +This module contains unit tests for the module package +cs_util.cosmo (moved from sp_validation.cosmology). + +:Author: Claude Code Assistant + +""" + +import numpy as np +import numpy.testing as npt +import pyccl as ccl +import pytest + +from cs_util import cosmo as 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/pyproject.toml b/pyproject.toml index 48b1764..dabdae3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "cs_util" -version = "0.2.1" +version = "0.2.2" description = "Utility library for CosmoStat" authors = [ { name = "Martin Kilbinger", email = "martin.kilbinger@cea.fr" }, @@ -15,7 +15,14 @@ dependencies = [ "numpy>=1.26.4", "matplotlib>=3.8.4", "swig>=4.2.1", - "scipy>=1.13.0", + # Cap scipy < 1.18: scipy 1.18 ported FITPACK from Fortran to C, changing + # RectBivariateSpline(scalar, scalar, grid=False) from a 0-d array to a + # shape-(1,) array. camb's BBN Y_He predictor then fails set_cosmology with + # "only 0-dimensional arrays can be converted to Python scalars", breaking + # every get_cosmo path (and pyccl's Boltzmann spline path likewise). camb + # 1.6.6 is unfixed; lift the cap when camb ships the scipy-1.18 fix. Mirrors + # the same cap in sp_validation. + "scipy>=1.13.0,<1.18", "vos>=3.6.1", "keyring>=25.2.0", "pyccl>=3.0.2", @@ -53,3 +60,9 @@ docs = [ "nbsphinx-link>=1.3", ] +[tool.pytest.ini_options] +markers = [ + "fast: quick unit test (the default for the cosmo theory-curve suite)", + "slow: slower integration / precision test", +] +