From 65517bd7265a0acbe8dc9f42ac370fb7bd2fc6e9 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 03:27:02 +0200 Subject: [PATCH 1/2] =?UTF-8?q?feat(blinding):=20firecrown=20likelihood=20?= =?UTF-8?q?+=20CAMB=E2=86=94CCL=20cross-check=20(PR=205)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The blinding box's theory engine and its validation (PRD #241 §5). blinding_likelihood.py — a minimal firecrown Gaussian likelihood over the tomographic coarse ξ± subset of a sacc_io analysis file. Smokescreen (PR 6) calls only its compute_theory_vector to form the Muir et al. shift t(hidden)−t(fid). Generic over N source bins; extracts the coarse ξ± block (past COSEBIs/pseudo-Cℓ/ρ/τ) with its aligned covariance before handing to firecrown, as firecrown models every point it reads and Smokescreen requires sacc.mean == likelihood data vector. TheoryConfig is the single config surface for the fiducial cosmology + model choices (HMCode2020 via CAMB, NLA IA off by default, normal-hierarchy neutrinos), defaults mirroring the CosmoSIS v1.4.6.3 _A_cell fiducial. theory_vector_sacc_order realigns firecrown's statistic-ordered vector to SACC insertion order via sacc_indices, independent of statistic order. test_camb_ccl_crosscheck.py — firecrown+CCL vs an independent CAMB-P(k)→CCL path predict the same ξ± at fixed cosmology. Settles σ8-vs-A_s inside the test (CAMB A_s iterated to share CCL's σ8; nominal A_s is ~3% off). Observed: ξ+ 0.21% / ξ− 0.10% at fiducial, 0.25%/0.12% at a blind-sized offset; tolerances 0.5%/1.0%. Validates P(k)+amplitude consistency, not CCL's Limber+Hankel independently (documented). Observed green: likelihood 7 fast + 6 slow; cross-check 1 fast + 4 slow. Co-Authored-By: Claude Fable 5 --- src/sp_validation/blinding_likelihood.py | 364 ++++++++++++++++++ .../tests/test_blinding_likelihood.py | 306 +++++++++++++++ .../tests/test_camb_ccl_crosscheck.py | 301 +++++++++++++++ 3 files changed, 971 insertions(+) create mode 100644 src/sp_validation/blinding_likelihood.py create mode 100644 src/sp_validation/tests/test_blinding_likelihood.py create mode 100644 src/sp_validation/tests/test_camb_ccl_crosscheck.py diff --git a/src/sp_validation/blinding_likelihood.py b/src/sp_validation/blinding_likelihood.py new file mode 100644 index 00000000..b69e807b --- /dev/null +++ b/src/sp_validation/blinding_likelihood.py @@ -0,0 +1,364 @@ +"""Firecrown cosmic-shear likelihood — the blinding box's theory engine. + +:Name: blinding_likelihood.py + +:Description: A minimal `firecrown `_ Gaussian + likelihood over the tomographic real-space shear 2PCF (ξ±) written by + :mod:`sp_validation.sacc_io`. It exists for one job: give Smokescreen + (PR 6) a ``compute_theory_vector`` aligned element-for-element with the + measured ξ± data vector, so the Muir et al. blinding shift + ``d → d + t(hidden) − t(fiducial)`` can be evaluated. The sampling + likelihood stays with CosmoSIS; only the theory-vector predict step of this + module is ever called. + + **Scope.** Only the *coarse analysis* ξ± subset (``grid='coarse'``). A + full analysis SACC file also carries pseudo-Cℓ, COSEBIs, pure-E/B and ρ/τ; + this likelihood extracts the coarse ξ± block (with its aligned covariance + sub-block) via :func:`sp_validation.sacc_io.extract` before handing it to + firecrown, because firecrown's ``read`` models every point it is given and + Smokescreen's ``_verify_sacc_consistency`` requires the SACC ``mean`` to + equal the likelihood's data vector element-for-element. The fine-grid ξ± + shift is computed in PR 6 as a direct CCL theory difference, not through + this likelihood. + + **N tomographic bins.** Generic over any number of ``source_i`` NZ tracers: + it builds one firecrown ``WeakLensing`` source per bin and one ``TwoPoint`` + (ξ+ and ξ−) statistic per unordered bin pair ``i ≤ j`` that the file + actually contains. Today's single-bin data is the ``n = 1`` case. + + **Configuration surface.** :class:`TheoryConfig` is the single, explicit + home for the fiducial cosmology and model choices (nonlinear P(k) recipe, + intrinsic-alignment model, neutrino treatment). Its defaults mirror the + ``cosmo_inference`` CosmoSIS fiducial (v1.4.6.3 ``_A_cell`` pipeline) so the + two theory stacks agree by construction; PRD #241's open question about the + fiducial choice changes these *values*, not this code. Pass a + :class:`TheoryConfig` (or a mapping of overrides) through the firecrown + ``build_parameters`` under the ``theory_config`` key to change it. +""" + +import dataclasses +import pathlib + +import numpy as np +import sacc + +from . import sacc_io + + +# --------------------------------------------------------------------------- # +# Configuration surface — the ONE place fiducial cosmology + model choices live +# --------------------------------------------------------------------------- # +@dataclasses.dataclass(frozen=True) +class TheoryConfig: + """Fiducial cosmology and model configuration for the theory engine. + + Every field is a deliberate, configurable choice. The defaults reproduce + the ``cosmo_inference`` CosmoSIS fiducial (the ``SP_v1.4.6.3_A_cell`` + pipeline + ``values_ia.ini`` central values): so the CCL theory computed + here and the CAMB theory CosmoSIS computes agree to the level the + CAMB↔CCL cross-check test asserts. PRD #241 §5 leaves the *named* fiducial + choice open for the group; adopting a different one is a change to these + values, not to any code. + + Cosmology is parametrised by ``S8`` and ``Omega_m`` — the axes the blind + shifts (PRD §3) — and converted to CCL's native ``sigma8``/``Omega_c`` at + build time. ``sigma8`` (not ``A_s``) is the CCL amplitude parameter; the + ~0.5% σ8-vs-A_s convention offset against a CAMB ``A_s`` input is settled + inside the cross-check test, not here. + """ + + # Cosmological parameters (blind axes S8, Omega_m + the rest). + S8: float = 0.80 # values_ia.ini S_8_input central + Omega_m: float = 0.30 # (omch2+ombh2)/h^2 with the fiducial baryons ≈ 0.29–0.31 + Omega_b: float = 0.0469 # ombh2=0.023 at h=0.7 -> 0.023/0.7^2 + h: float = 0.70 # h0 central + n_s: float = 0.96 # n_s central + m_nu: float = 0.06 # single massive neutrino, eV (mnu in values_ia.ini) + w0: float = -1.0 + wa: float = 0.0 + + # Nonlinear matter power spectrum: CAMB's HMCode2020 with baryonic + # feedback, matching the CosmoSIS `[camb] halofit_version=mead2020_feedback` + # line. CCL routes P(k) through CAMB when `matter_power_spectrum` resolves + # to CAMB (via firecrown's CCLFactory camb_extra_params), so both stacks + # share the identical nonlinear recipe. + halofit_version: str = "mead2020_feedback" + hmcode_logT_AGN: float = 7.5 # values_ia.ini logT_AGN central + + # Neutrino mass split: normal hierarchy (CosmoSIS `neutrino_hierarchy=normal`). + mass_split: str = "normal" + + # Intrinsic alignments: NLA (linear alignment). The theory engine's sole + # purpose is the (Omega_m, S8) blinding shift t(hidden) − t(fiducial); IA is + # a nuisance that enters identically at both cosmologies and nearly cancels, + # so the fiducial engine defaults IA OFF (ia_bias=0) to isolate the lensing + # prediction — this is also what makes the CAMB↔CCL cross-check a clean test + # of the shear calculation. Set `ia_bias` nonzero (CosmoSIS central A=1.0) + # to include NLA; `ia_z_piv`/`ia_alphaz` set the pivot and z-power. + ia_bias: float = 0.0 + ia_z_piv: float = 0.62 + ia_alphaz: float = 0.0 + + # Derived / fixed CCL background parameters. + Neff: float = 3.046 + T_CMB: float = 2.7255 + + def sigma8(self): + """CCL ``sigma8`` implied by ``S8`` and ``Omega_m``. + + ``S8 ≡ σ8 √(Ωm / 0.3)`` — the standard weak-lensing definition — so + ``σ8 = S8 / √(Ωm / 0.3)``. + """ + return self.S8 / np.sqrt(self.Omega_m / 0.3) + + def omega_c(self): + """CCL cold-dark-matter density ``Omega_c = Omega_m − Omega_b − Omega_nu``. + + The neutrino density ``Ω_ν h² = Σm_ν / 93.14 eV`` is subtracted so that + the *total* matter density is exactly ``Omega_m`` (CCL treats massive + neutrinos as a separate species, not part of ``Omega_c``). + """ + omega_nu = self.m_nu / (93.14 * self.h**2) + return self.Omega_m - self.Omega_b - omega_nu + + def ccl_params(self): + """Cosmological sampling parameters as firecrown expects them. + + Returns the ``{name: value}`` overrides to feed a firecrown + ``ParamsMap`` so the modeling tools build the CCL cosmology at this + fiducial. Names match firecrown's CCLFactory parameter set + (``Omega_c``, ``Omega_b``, ``h``, ``n_s``, ``sigma8``, ``m_nu``, …). + """ + return { + "Omega_c": self.omega_c(), + "Omega_b": self.Omega_b, + "h": self.h, + "n_s": self.n_s, + "sigma8": self.sigma8(), + "m_nu": self.m_nu, + "w0": self.w0, + "wa": self.wa, + "Neff": self.Neff, + "T_CMB": self.T_CMB, + "Omega_k": 0.0, + } + + def ia_params(self): + """NLA intrinsic-alignment sampling parameters (firecrown names).""" + return { + "ia_bias": self.ia_bias, + "z_piv": self.ia_z_piv, + "alphaz": self.ia_alphaz, + } + + @classmethod + def from_overrides(cls, overrides): + """Build from a mapping of field overrides (fail loud on unknown keys).""" + fields = {f.name for f in dataclasses.fields(cls)} + unknown = set(overrides) - fields + if unknown: + raise ValueError( + f"unknown TheoryConfig fields {sorted(unknown)}; " + f"valid fields are {sorted(fields)}" + ) + return cls(**overrides) + + +# --------------------------------------------------------------------------- # +# SACC preparation: extract the coarse ξ± analysis subset +# --------------------------------------------------------------------------- # +def source_bins(s): + """Sorted source-bin indices present in ``s`` (from ``source_i`` tracers).""" + bins = [] + for name in s.tracers: + if name.startswith("source_"): + bins.append(int(name.split("_", 1)[1])) + return sorted(bins) + + +def xi_pairs(s, grid="coarse"): + """Unordered source-bin pairs ``(i, j)`` with ``i ≤ j`` carrying coarse ξ+. + + Discovered from the file rather than assumed, so the likelihood models + exactly the pairs the measurement contains (single-bin ⇒ ``[(0, 0)]``). + """ + bins = source_bins(s) + pairs = [] + for a, i in enumerate(bins): + for j in bins[a:]: + tracers = (sacc_io.source_name(i), sacc_io.source_name(j)) + if len(s.indices(sacc_io.XI_PLUS, tracers, grid=grid)): + pairs.append((i, j)) + return pairs + + +def extract_coarse_xi(s): + """Return a SACC holding only the coarse ξ± subset (+ aligned covariance). + + firecrown models every point handed to ``read`` and Smokescreen requires + ``sacc.mean == likelihood.get_data_vector()``, so the likelihood must + receive a file containing *exactly* the ξ± rows it predicts — no COSEBIs, + pseudo-Cℓ, ρ/τ. This keeps ξ+ and ξ− on the coarse grid for all tracers + (in insertion order, so the retained covariance sub-block stays aligned) + and drops everything else. ``keep_indices`` slices the covariance to match, + exactly as :func:`sp_validation.sacc_io.extract` relies on. + """ + keep = np.sort( + np.concatenate( + [ + s.indices(sacc_io.XI_PLUS, grid="coarse"), + s.indices(sacc_io.XI_MINUS, grid="coarse"), + ] + ) + ) + sub = s.copy() + sub.keep_indices(keep) + return sub + + +# --------------------------------------------------------------------------- # +# firecrown build_likelihood entry point +# --------------------------------------------------------------------------- # +def _resolve_sacc(sacc_data): + """Load a SACC from a path, or copy an in-memory Sacc.""" + if isinstance(sacc_data, (str, pathlib.Path)): + return sacc.Sacc.load_fits(str(sacc_data)) + return sacc_data.copy() + + +def _resolve_config(build_parameters): + """Extract a TheoryConfig from firecrown build parameters (default if absent).""" + try: + cfg = build_parameters["theory_config"] + except (KeyError, TypeError): + return TheoryConfig() + if isinstance(cfg, TheoryConfig): + return cfg + return TheoryConfig.from_overrides(dict(cfg)) + + +def build_likelihood(build_parameters): + """Firecrown factory: a Gaussian ξ± likelihood over the coarse analysis grid. + + Parameters + ---------- + build_parameters : firecrown NamedParameters + Must carry ``sacc_data`` (a SACC path or in-memory ``sacc.Sacc`` — the + *analysis* file; the coarse ξ± subset is extracted here). May carry + ``theory_config`` (a :class:`TheoryConfig` or a mapping of overrides); + absent ⇒ the fiducial defaults. + + Returns + ------- + (firecrown.likelihood.ConstGaussian, firecrown.modeling_tools.ModelingTools) + The Gaussian likelihood and the modeling tools whose CCL cosmology is + configured from the :class:`TheoryConfig`. + """ + from firecrown.likelihood import ConstGaussian, TwoPoint + from firecrown.likelihood import weak_lensing as wl + from firecrown.modeling_tools import ( + CAMBExtraParams, + CCLCreationMode, + CCLFactory, + CCLPureModeTransferFunction, + ModelingTools, + PoweSpecAmplitudeParameter, + ) + + try: + raw = build_parameters["sacc_data"] + except TypeError: + raw = build_parameters.data["sacc_data"] + s = _resolve_sacc(raw) + cfg = _resolve_config(build_parameters) + + coarse = extract_coarse_xi(s) + pairs = xi_pairs(coarse, grid="coarse") + if not pairs: + raise ValueError( + "no coarse ξ± data found in the SACC file — the blinding " + "likelihood needs galaxy_shear_xi_plus/minus points tagged grid='coarse'" + ) + + # One WeakLensing source per bin (NLA systematic wired, off by default). + bins = sorted({b for pair in pairs for b in pair}) + sources = { + i: wl.WeakLensing( + sacc_tracer=sacc_io.source_name(i), + systematics=[wl.LinearAlignmentSystematic(sacc_tracer=None)], + ) + for i in bins + } + + # One (ξ+, ξ−) TwoPoint per pair. Statistic order does not constrain + # correctness: each statistic exposes ``sacc_indices`` mapping its theory + # sub-vector back to SACC rows, so callers realign by those indices. + stats = [] + for i, j in pairs: + stats.append(TwoPoint(sacc_io.XI_PLUS, sources[i], sources[j])) + stats.append(TwoPoint(sacc_io.XI_MINUS, sources[i], sources[j])) + + likelihood = ConstGaussian(statistics=stats) + likelihood.read(coarse) + + # Route the nonlinear P(k) through CAMB's HMCode2020 (matching CosmoSIS's + # `halofit_version=mead2020_feedback`). firecrown gates camb_extra_params + # behind PURE_CCL_MODE + the BOLTZMANN_CAMB transfer function — the pure-CCL + # creation path that builds the ccl.Cosmology directly with CAMB as its + # power-spectrum backend. + tools = ModelingTools( + ccl_factory=CCLFactory( + require_nonlinear_pk=True, + amplitude_parameter=PoweSpecAmplitudeParameter.SIGMA8, + mass_split=cfg.mass_split, + creation_mode=CCLCreationMode.PURE_CCL_MODE, + pure_ccl_transfer_function=CCLPureModeTransferFunction.BOLTZMANN_CAMB, + camb_extra_params=CAMBExtraParams( + halofit_version=cfg.halofit_version, + HMCode_logT_AGN=cfg.hmcode_logT_AGN, + ), + ) + ) + return likelihood, tools + + +def compute_theory_vector(likelihood, tools, config=None): + """Predict the ξ± theory vector at ``config`` (default fiducial). + + A thin, deterministic driver over firecrown's update/prepare/predict cycle, + for the cross-check test and callers that want a theory vector without + Smokescreen. Smokescreen itself drives ``compute_theory_vector`` internally; + this mirrors that path so tests exercise the same code. + + Returns the theory vector in firecrown's statistic order; use + :func:`theory_vector_sacc_order` to realign it to SACC insertion order. + """ + from firecrown.updatable import ParamsMap, get_default_params_map + + cfg = config or TheoryConfig() + params = get_default_params_map(tools, likelihood) + params = ParamsMap({**dict(params), **cfg.ccl_params(), **cfg.ia_params()}) + tools.update(params) + tools.prepare() + likelihood.update(params) + tv = np.asarray(likelihood.compute_theory_vector(tools)) + likelihood.reset() + tools.reset() + return tv + + +def theory_vector_sacc_order(likelihood, theory_vector, n_sacc): + """Scatter a firecrown theory vector into SACC insertion order. + + firecrown returns the theory vector in the order its statistics were listed; + each statistic's ``sacc_indices`` say which SACC rows its block maps to. + Scattering by those indices yields a vector aligned to the SACC data vector + (and hence to the covariance) regardless of statistic order. + """ + out = np.full(n_sacc, np.nan) + cursor = 0 + for stat in likelihood.statistics: + idx = np.asarray(stat.statistic.sacc_indices) + out[idx] = theory_vector[cursor : cursor + len(idx)] + cursor += len(idx) + return out diff --git a/src/sp_validation/tests/test_blinding_likelihood.py b/src/sp_validation/tests/test_blinding_likelihood.py new file mode 100644 index 00000000..db010201 --- /dev/null +++ b/src/sp_validation/tests/test_blinding_likelihood.py @@ -0,0 +1,306 @@ +"""Tests for :mod:`sp_validation.blinding_likelihood`. + +The likelihood is the blinding box's theory engine (PRD #241 §5). These tests +prove it (1) builds from a :mod:`sp_validation.sacc_io`-written analysis file, +extracting the coarse ξ± subset past COSEBIs/ρ/τ; (2) returns a theory vector +of the right shape, aligned to SACC insertion order via ``sacc_indices`` +regardless of statistic order; (3) is deterministic; (4) is generic over N +tomographic bins. All synthetic — deterministic Gaussian n(z), no real data. + +The theory-computing tests import CCL/CAMB and are marked ``slow`` (each CCL +cosmology build + ξ± predict runs CAMB's HMCode). The construction/alignment +tests are fast and run in the default suite. +""" + +import numpy as np +import pytest + +from sp_validation import blinding_likelihood as bl +from sp_validation import sacc_io as sio + + +# --------------------------------------------------------------------------- # +# Synthetic fixtures — deterministic Gaussian n(z), sacc_io-written files +# --------------------------------------------------------------------------- # +def _gauss_nz(z0, sigma, n=300): + """Normalised Gaussian n(z) on a fixed grid (deterministic).""" + z = np.linspace(0.0, 3.0, n) + nz = np.exp(-0.5 * ((z - z0) / sigma) ** 2) + return z, nz / np.trapezoid(nz, z) + + +def _theta(n=10): + return np.geomspace(5.0, 250.0, n) + + +def _analysis_sacc(nbins=1, with_extras=True, ntheta=10): + """A sacc_io analysis file: coarse ξ± for every i≤j pair, plus (optionally) + COSEBIs/ρ so the extractor is exercised against a realistic mixed file. + + A FullCovariance is assembled block-diagonally (ξ block, then any extras), + so the file is a valid firecrown ``ConstGaussian`` input. + """ + theta = _theta(ntheta) + nz = {i: _gauss_nz(0.5 + 0.3 * i, 0.15 + 0.03 * i) for i in range(nbins)} + s = sio.new_sacc(nz) + pairs = [(i, j) for i in range(nbins) for j in range(i, nbins)] + for k, (i, j) in enumerate(pairs): + xip = (np.arange(ntheta) + 1) * (k + 1) * 1e-5 + xim = (np.arange(ntheta) + 1) * (k + 1) * 0.5e-5 + sio.add_xi(s, (i, j), theta, xip, xim, grid="coarse") + + blocks = [] + cursor_pairs = pairs + for i, j in cursor_pairs: + tr = (sio.source_name(i), sio.source_name(j)) + idx = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)]) + blocks.append((idx, np.eye(len(idx)) * 1e-12)) + + if with_extras: + sio.add_cosebis( + s, (0, 0), np.arange(1, 5) * 1e-6, np.arange(1, 5) * 1e-7, (1.0, 100.0) + ) + for kk in range(3): + sio.add_rho( + s, kk, theta, np.arange(ntheta) * 1e-7, np.arange(ntheta) * 2e-7 + ) + tr = ("source_0", "source_0") + co = np.concatenate( + [s.indices(sio.COSEBI_EE, tr), s.indices(sio.COSEBI_BB, tr)] + ) + blocks.append((co, np.eye(len(co)) * 1e-14)) + for kk in range(3): + ridx = np.concatenate( + [ + s.indices(sio.RHO_PLUS.format(k=kk)), + s.indices(sio.RHO_MINUS.format(k=kk)), + ] + ) + blocks.append((ridx, np.eye(len(ridx)) * 1e-16)) + + sio.assemble_covariance(s, blocks) + return s + + +# --------------------------------------------------------------------------- # +# 1. TheoryConfig: the single config surface +# --------------------------------------------------------------------------- # +def test_theory_config_defaults_match_cosmosis_fiducial(): + """Defaults reproduce the CosmoSIS v1.4.6.3 `_A_cell` fiducial values.""" + c = bl.TheoryConfig() + assert c.S8 == 0.80 + assert c.n_s == 0.96 + assert c.h == 0.70 + assert c.m_nu == 0.06 + assert c.halofit_version == "mead2020_feedback" + assert c.mass_split == "normal" + # σ8 = S8 / √(Ωm/0.3); at the fiducial Ωm this is ~0.80. + assert c.sigma8() == pytest.approx(c.S8 / np.sqrt(c.Omega_m / 0.3)) + # Ωc = Ωm − Ωb − Ων, all positive and summing back to Ωm. + omega_nu = c.m_nu / (93.14 * c.h**2) + assert c.omega_c() == pytest.approx(c.Omega_m - c.Omega_b - omega_nu) + assert 0 < c.omega_c() < c.Omega_m + + +def test_theory_config_overrides_fail_loud(): + """Unknown override keys raise (fail-fast config).""" + cfg = bl.TheoryConfig.from_overrides({"S8": 0.75, "Omega_m": 0.32}) + assert cfg.S8 == 0.75 and cfg.Omega_m == 0.32 + with pytest.raises(ValueError, match="unknown TheoryConfig fields"): + bl.TheoryConfig.from_overrides({"S8": 0.75, "sig8": 0.8}) + + +def test_theory_config_is_frozen(): + """The config is immutable — one authoritative fiducial per build.""" + c = bl.TheoryConfig() + with pytest.raises(Exception): + c.S8 = 0.9 # frozen dataclass -> FrozenInstanceError + + +# --------------------------------------------------------------------------- # +# 2. SACC preparation: bin/pair discovery + coarse-ξ extraction +# --------------------------------------------------------------------------- # +def test_source_bins_and_pairs_discovered_from_file(): + s = _analysis_sacc(nbins=3, with_extras=False) + assert bl.source_bins(s) == [0, 1, 2] + assert bl.xi_pairs(s, grid="coarse") == [ + (0, 0), + (0, 1), + (0, 2), + (1, 1), + (1, 2), + (2, 2), + ] + + +def test_extract_coarse_xi_drops_everything_else(): + """Extraction keeps only coarse ξ±; COSEBIs/ρ and any fine grid are gone. + + Built without a covariance so the fine grid can be appended (SACC forbids + adding points after ``add_covariance``); extraction logic is + covariance-independent — the covariance-slicing path is covered by + :func:`test_extract_coarse_xi_covariance_sliced`. + """ + theta = _theta(10) + s = sio.new_sacc({i: _gauss_nz(0.5 + 0.3 * i, 0.15) for i in range(2)}) + for k, (i, j) in enumerate([(0, 0), (0, 1), (1, 1)]): + sio.add_xi( + s, + (i, j), + theta, + (np.arange(10) + 1) * (k + 1) * 1e-5, + (np.arange(10) + 1) * (k + 1) * 0.5e-5, + grid="coarse", + ) + sio.add_cosebis( + s, (0, 0), np.arange(1, 5) * 1e-6, np.arange(1, 5) * 1e-7, (1.0, 100.0) + ) + for kk in range(3): + sio.add_rho(s, kk, theta, np.arange(10) * 1e-7, np.arange(10) * 2e-7) + # a fine grid for (0,0) to confirm it is dropped + tf = np.geomspace(0.1, 250.0, 40) + sio.add_xi(s, (0, 0), tf, np.arange(40) * 1e-5, np.arange(40) * 2e-5, grid="fine") + + coarse = bl.extract_coarse_xi(s) + dtypes = {d.data_type for d in coarse.data} + grids = {d.tags.get("grid") for d in coarse.data} + assert dtypes == {sio.XI_PLUS, sio.XI_MINUS} + assert grids == {"coarse"} + # 3 coarse pairs x 2 types x 10 theta + assert len(coarse.mean) == 3 * 2 * 10 + + +def test_extract_coarse_xi_covariance_sliced(): + """Extraction carries an aligned covariance sub-block for the ξ± rows.""" + s = _analysis_sacc(nbins=2, with_extras=True) + coarse = bl.extract_coarse_xi(s) + n = len(coarse.mean) + assert coarse.covariance.dense.shape == (n, n) + # the ξ blocks in the source file were identity*1e-12; the extracted + # sub-block is exactly that (diagonal, 1e-12), COSEBI/ρ rows excluded + assert np.allclose(np.diag(coarse.covariance.dense), 1e-12) + + +def test_extract_preserves_original(): + s = _analysis_sacc(nbins=1, with_extras=True) + n_before = len(s.mean) + _ = bl.extract_coarse_xi(s) + assert len(s.mean) == n_before + + +# --------------------------------------------------------------------------- # +# 3. build_likelihood: construction from a sacc_io file +# --------------------------------------------------------------------------- # +def _named(sacc_data, **kw): + from firecrown.likelihood import NamedParameters + + return NamedParameters({"sacc_data": sacc_data, **kw}) + + +@pytest.mark.slow +def test_build_likelihood_single_bin(): + s = _analysis_sacc(nbins=1, with_extras=True) + lik, tools = bl.build_likelihood(_named(s)) + # one (ξ+, ξ−) TwoPoint per pair -> 2 statistics for the single-bin case + assert len(lik.statistics) == 2 + + +@pytest.mark.slow +def test_build_likelihood_requires_coarse_xi(): + """A file with no coarse ξ± is a loud error, not a silent empty likelihood.""" + s = sio.new_sacc({0: _gauss_nz(0.7, 0.2)}) + sio.add_cosebis( + s, (0, 0), np.arange(1, 4) * 1e-6, np.arange(1, 4) * 1e-7, (1.0, 100.0) + ) + with pytest.raises(ValueError, match="no coarse"): + bl.build_likelihood(_named(s)) + + +# --------------------------------------------------------------------------- # +# 4. Theory vector: shape, finiteness, determinism, SACC alignment +# --------------------------------------------------------------------------- # +@pytest.mark.slow +def test_theory_vector_shape_and_finite(): + s = _analysis_sacc(nbins=1, with_extras=True) + lik, tools = bl.build_likelihood(_named(s)) + tv = bl.compute_theory_vector(lik, tools) + assert tv.shape == (2 * 10,) # ξ+ and ξ− on 10 θ + assert np.all(np.isfinite(tv)) + assert np.any(tv != 0.0) + # cosmic-shear ξ+ is positive and falls with θ (monotone over this range) + xip = tv[:10] + assert np.all(xip > 0) + assert xip[0] > xip[-1] + + +@pytest.mark.slow +def test_theory_vector_deterministic(): + """Same inputs -> bitwise-identical theory vector.""" + s = _analysis_sacc(nbins=1, with_extras=True) + lik, tools = bl.build_likelihood(_named(s)) + tv1 = bl.compute_theory_vector(lik, tools) + tv2 = bl.compute_theory_vector(lik, tools) + assert np.array_equal(tv1, tv2) + + +@pytest.mark.slow +def test_theory_vector_realigns_to_sacc_order_regardless_of_statistic_order(): + """The `sacc_indices` scatter recovers SACC alignment for any statistic order. + + Two likelihoods over the same file — one with the default (pair-major) + statistic order, one deliberately reversed — must, after scattering by + ``sacc_indices``, produce the SAME SACC-aligned theory vector. This is the + property that lets Smokescreen and the covariance stay aligned without the + likelihood pinning a particular statistic order. + """ + from firecrown.likelihood import ConstGaussian, TwoPoint + from firecrown.likelihood import weak_lensing as wl + + s = _analysis_sacc(nbins=2, with_extras=False) + coarse = bl.extract_coarse_xi(s) + n = len(coarse.mean) + + lik_default, tools = bl.build_likelihood(_named(s)) + tv_default = bl.compute_theory_vector(lik_default, tools) + aligned_default = bl.theory_vector_sacc_order(lik_default, tv_default, n) + + # Rebuild with statistics in reversed order. + pairs = bl.xi_pairs(coarse, grid="coarse") + bins = sorted({b for p in pairs for b in p}) + sources = { + i: wl.WeakLensing( + sacc_tracer=sio.source_name(i), + systematics=[wl.LinearAlignmentSystematic(sacc_tracer=None)], + ) + for i in bins + } + stats = [] + for i, j in pairs: + stats.append(TwoPoint(sio.XI_PLUS, sources[i], sources[j])) + stats.append(TwoPoint(sio.XI_MINUS, sources[i], sources[j])) + lik_rev = ConstGaussian(statistics=list(reversed(stats))) + lik_rev.read(coarse) + _, tools_rev = bl.build_likelihood(_named(s)) + tv_rev = bl.compute_theory_vector(lik_rev, tools_rev) + aligned_rev = bl.theory_vector_sacc_order(lik_rev, tv_rev, n) + + assert np.all(np.isfinite(aligned_default)) + assert np.array_equal(aligned_default, aligned_rev) + + +# --------------------------------------------------------------------------- # +# 5. N-bin generality: a 2-bin file yields 3 pairs / 6 statistics +# --------------------------------------------------------------------------- # +@pytest.mark.slow +def test_two_bin_generality(): + s = _analysis_sacc(nbins=2, with_extras=False) + lik, tools = bl.build_likelihood(_named(s)) + # pairs (0,0),(0,1),(1,1) -> 3 pairs x {ξ+, ξ−} = 6 statistics + assert len(lik.statistics) == 6 + tv = bl.compute_theory_vector(lik, tools) + assert tv.shape == (6 * 10,) + assert np.all(np.isfinite(tv)) and np.any(tv != 0.0) + # realignment covers exactly the coarse ξ rows + coarse = bl.extract_coarse_xi(s) + aligned = bl.theory_vector_sacc_order(lik, tv, len(coarse.mean)) + assert np.all(np.isfinite(aligned)) diff --git a/src/sp_validation/tests/test_camb_ccl_crosscheck.py b/src/sp_validation/tests/test_camb_ccl_crosscheck.py new file mode 100644 index 00000000..4a34698c --- /dev/null +++ b/src/sp_validation/tests/test_camb_ccl_crosscheck.py @@ -0,0 +1,301 @@ +"""CAMB↔CCL theory cross-check (PRD #241 §5). + +The blinding box computes theory through firecrown + CCL, while cosmological +inference runs CosmoSIS + CAMB. The Muir et al. shift ``t(hidden) − t(fid)`` +only means what we intend if the two stacks predict the same ξ± at a fixed +cosmology. This module asserts that agreement. + +**What is compared.** For one cosmology and n(z), ξ± on our θ grid computed two +ways: + +- *Path A — the blinding stack:* :func:`sp_validation.blinding_likelihood` + → firecrown → CCL, with the nonlinear P(k) routed through CAMB's HMCode2020 + (PURE_CCL_MODE + ``camb_extra_params``), the exact recipe the likelihood + uses in production. +- *Path B — an independent CAMB path:* a direct pycamb run producing the + HMCode2020 P(k, z), wrapped in a ``ccl.Pk2D`` and projected to ξ± through + CCL's Limber (``angular_cl``) + FFTLog (``correlation``). + +**The σ8-vs-A_s subtlety, settled here.** Our conventions set σ8 for CCL but +A_s for CAMB. A nominal A_s (2.1e-9) gives CAMB a σ8 that is ~3% off our +target — the amplitudes would silently disagree. Path B therefore iterates +CAMB's A_s until CAMB's σ8 equals the CCL target (σ8² ∝ A_s, so one Newton +step nails it), so both stacks share σ8 to <1e-4 before any ξ± is computed. +That reconciliation is the point of doing this in the test rather than +hoping the amplitudes line up. + +**What this validates — and does not.** It validates that the blinding +theory agrees with an independent CAMB-sourced ξ± at the sub-percent level: +the P(k) recipe (HMCode2020) is consistent across stacks and the amplitude +convention is correctly reconciled. It does *not* independently validate +CCL's Limber+Hankel machinery — both paths use ``ccl.angular_cl`` + +``ccl.correlation`` (a standalone FFTLog library is not in the environment), +so a common bug there would cancel. The residual we observe (~0.2–0.6%, +larger for ξ− which weights small scales) is dominated by P(k)-representation +and ℓ-grid differences between two independent-but-both-CCL projections, not +by cosmology; it is the realistic floor of cross-implementation agreement, +and the tolerances sit comfortably above it. +""" + +import numpy as np +import pytest + +from sp_validation import blinding_likelihood as bl +from sp_validation import sacc_io as sio + + +# --------------------------------------------------------------------------- # +# Deterministic fixture: one Gaussian source bin +# --------------------------------------------------------------------------- # +def _gauss_nz(z0=0.7, sigma=0.2, n=400): + z = np.linspace(0.01, 3.0, n) + nz = np.exp(-0.5 * ((z - z0) / sigma) ** 2) + return z, nz / np.trapezoid(nz, z) + + +# --------------------------------------------------------------------------- # +# Path A — firecrown + CCL (the blinding stack) +# --------------------------------------------------------------------------- # +def _firecrown_xi(z, nz, theta_arcmin, config): + """ξ± from the production likelihood at ``config`` — the blinding stack.""" + from firecrown.likelihood import NamedParameters + + s = sio.new_sacc({0: (z, nz)}) + zeros = np.zeros_like(theta_arcmin) + sio.add_xi(s, (0, 0), theta_arcmin, zeros, zeros, grid="coarse") + tr = ("source_0", "source_0") + xi = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)]) + sio.assemble_covariance(s, [(xi, np.eye(len(xi)) * 1e-12)]) + + lik, tools = bl.build_likelihood( + NamedParameters({"sacc_data": s, "theory_config": config}) + ) + tv = bl.compute_theory_vector(lik, tools, config) + n = len(theta_arcmin) + return tv[:n], tv[n:] + + +# --------------------------------------------------------------------------- # +# Path B — independent CAMB P(k) → CCL Pk2D → ξ± +# --------------------------------------------------------------------------- # +def _camb_As_for_sigma8(config, sigma8_target, n_z=48, kmax=20.0): + """CAMB ``A_s`` giving ``sigma8_target``, and the CAMB param builder. + + Iterates A_s (σ8² ∝ A_s ⇒ Newton converges in one step) so the CAMB stack + shares σ8 with CCL. Returns ``(A_s, make_params)`` where ``make_params(As, + nonlinear)`` builds a CAMBparams at the fiducial background. + """ + import camb + + h, ns, mnu = config.h, config.n_s, config.m_nu + omch2 = config.omega_c() * h**2 + ombh2 = config.Omega_b * h**2 + zlist = list(np.linspace(0.0, 3.0, n_z)) + + def make_params(As, nonlinear): + p = camb.CAMBparams() + p.set_cosmology( + H0=h * 100, + ombh2=ombh2, + omch2=omch2, + mnu=mnu, + num_massive_neutrinos=1, + neutrino_hierarchy=config.mass_split, + ) + p.InitPower.set_params(As=As, ns=ns) + p.set_matter_power(redshifts=zlist, kmax=kmax) + if nonlinear: + p.NonLinear = camb.model.NonLinear_both + p.NonLinearModel.set_params( + halofit_version=config.halofit_version, + HMCode_logT_AGN=config.hmcode_logT_AGN, + ) + else: + p.NonLinear = camb.model.NonLinear_none + return p + + As = 2.1e-9 + for _ in range(3): # linear σ8 for the amplitude solve; 1 step suffices + r = camb.get_results(make_params(As, nonlinear=False)) + As *= (sigma8_target / r.get_sigma8_0()) ** 2 + return As, make_params + + +def _independent_camb_xi( + z, nz, theta_arcmin, config, n_ell=300, ell_max=60000, kmax=20.0 +): + """ξ± from an independent CAMB P(k), projected through CCL (σ8-matched).""" + import camb + import pyccl as ccl + + sigma8 = config.sigma8() + As, make_params = _camb_As_for_sigma8(config, sigma8, kmax=kmax) + results = camb.get_results(make_params(As, nonlinear=True)) + kh, zc, pk = results.get_matter_power_spectrum(minkh=1e-4, maxkh=kmax, npoints=400) + zc = np.asarray(zc) + # CAMB uses h-units (k in h/Mpc, P in (Mpc/h)^3); CCL wants 1/Mpc, Mpc^3. + k = kh * config.h + Pk = pk / config.h**3 + a = 1.0 / (1.0 + zc) + order = np.argsort(a) # Pk2D wants ascending scale factor + pk2d = ccl.Pk2D( + a_arr=a[order], + lk_arr=np.log(k), + pk_arr=np.log(Pk[order]), + is_logp=True, + ) + + cosmo = ccl.Cosmology( + Omega_c=config.omega_c(), + Omega_b=config.Omega_b, + h=config.h, + n_s=config.n_s, + sigma8=sigma8, + m_nu=config.m_nu, + mass_split="single", + matter_power_spectrum="camb", + extra_parameters={ + "camb": { + "halofit_version": config.halofit_version, + "HMCode_logT_AGN": config.hmcode_logT_AGN, + } + }, + ) + lens = ccl.WeakLensingTracer(cosmo, dndz=(z, nz)) + ells = np.unique(np.geomspace(2, ell_max, n_ell).astype(int)).astype(float) + cl = ccl.angular_cl(cosmo, lens, lens, ells, p_of_k_a=pk2d) + theta_deg = theta_arcmin / 60.0 + xip = ccl.correlation(cosmo, ell=ells, C_ell=cl, theta=theta_deg, type="GG+") + xim = ccl.correlation(cosmo, ell=ells, C_ell=cl, theta=theta_deg, type="GG-") + return xip, xim, As + + +# --------------------------------------------------------------------------- # +# The σ8/A_s reconciliation — fast, no shear projection +# --------------------------------------------------------------------------- # +@pytest.mark.slow +def test_sigma8_As_reconciliation(): + """CAMB's A_s must be iterated to share σ8 with CCL — the ~3% offset is real. + + Proves (a) a nominal A_s leaves CAMB's σ8 percent-level off our σ8 target, + so the convention offset is not negligible; (b) the Newton iteration lands + CAMB's σ8 on the target to <1e-4. This is what makes the ξ± cross-check a + same-amplitude comparison rather than an amplitude-mismatch artifact. + """ + import camb + + cfg = bl.TheoryConfig() + target = cfg.sigma8() + _, make_params = _camb_As_for_sigma8(cfg, target) + + nominal = camb.get_results(make_params(2.1e-9, nonlinear=False)).get_sigma8_0() + assert abs(nominal - target) / target > 0.01 # >1% off at nominal A_s + + As, _ = _camb_As_for_sigma8(cfg, target) + matched = camb.get_results(make_params(As, nonlinear=False)).get_sigma8_0() + assert abs(matched - target) < 1e-4 + + +# --------------------------------------------------------------------------- # +# The cross-check: firecrown+CCL vs independent CAMB, σ8-matched +# --------------------------------------------------------------------------- # +# Observed agreement floor on this fixture (12 θ ∈ [5, 250] arcmin, HMCode2020): +# ξ+ ≈ 0.3%, ξ− ≈ 0.6% (converged in ℓ). Tolerances sit above the floor with +# headroom; ξ− is looser because it weights smaller scales, where the +# P(k)-representation and FFTLog differences between the two paths are largest. +XIP_RTOL = 0.005 # 0.5% +XIM_RTOL = 0.010 # 1.0% + + +@pytest.mark.slow +def test_camb_ccl_xi_agreement(): + """firecrown+CCL and independent-CAMB ξ± agree within the stated tolerance.""" + cfg = bl.TheoryConfig() + z, nz = _gauss_nz() + theta_arcmin = np.geomspace(5.0, 250.0, 12) + + fc_xip, fc_xim = _firecrown_xi(z, nz, theta_arcmin, cfg) + ind_xip, ind_xim, _ = _independent_camb_xi(z, nz, theta_arcmin, cfg) + + rel_p = np.abs(ind_xip - fc_xip) / np.abs(fc_xip) + rel_m = np.abs(ind_xim - fc_xim) / np.abs(fc_xim) + assert np.all(fc_xip > 0) and np.all(ind_xip > 0) # sensible cosmic shear + assert rel_p.max() < XIP_RTOL, f"ξ+ max rel diff {rel_p.max():.3%} ≥ {XIP_RTOL:.1%}" + assert rel_m.max() < XIM_RTOL, f"ξ− max rel diff {rel_m.max():.3%} ≥ {XIM_RTOL:.1%}" + + +@pytest.mark.slow +def test_camb_ccl_xi_agreement_shifted_cosmology(): + """Agreement holds away from the fiducial (a blind-sized S8/Ωm offset). + + The blinding shift moves the cosmology by up to ±0.075 in S8 and ±0.1 in + Ωm (PRD §3); the two stacks must still agree there, or the *shift itself* + (a difference of two theory vectors) inherits a stack-disagreement bias. + """ + cfg = bl.TheoryConfig.from_overrides({"S8": 0.80 + 0.075, "Omega_m": 0.30 - 0.05}) + z, nz = _gauss_nz() + theta_arcmin = np.geomspace(5.0, 250.0, 12) + + fc_xip, fc_xim = _firecrown_xi(z, nz, theta_arcmin, cfg) + ind_xip, ind_xim, _ = _independent_camb_xi(z, nz, theta_arcmin, cfg) + + rel_p = np.abs(ind_xip - fc_xip) / np.abs(fc_xip) + rel_m = np.abs(ind_xim - fc_xim) / np.abs(fc_xim) + assert rel_p.max() < XIP_RTOL, f"ξ+ max rel diff {rel_p.max():.3%} ≥ {XIP_RTOL:.1%}" + assert rel_m.max() < XIM_RTOL, f"ξ− max rel diff {rel_m.max():.3%} ≥ {XIM_RTOL:.1%}" + + +@pytest.mark.slow +def test_blinding_shift_consistent_across_stacks(): + """The *shift* d(hidden) − d(fiducial) agrees across stacks. + + What blinding actually adds to the data is a difference of two theory + vectors. Even if each vector carries a small common-mode stack offset, the + difference should agree more tightly than the absolute vectors — this is + the quantity that matters. Assert the shift (fiducial → a blind-sized + hidden cosmology) matches between firecrown+CCL and independent-CAMB to + within the ξ± absolute tolerance. + """ + fid = bl.TheoryConfig() + hidden = bl.TheoryConfig.from_overrides({"S8": 0.80 - 0.05, "Omega_m": 0.30 + 0.03}) + z, nz = _gauss_nz() + theta_arcmin = np.geomspace(5.0, 250.0, 12) + + fc_fid_p, fc_fid_m = _firecrown_xi(z, nz, theta_arcmin, fid) + fc_hid_p, fc_hid_m = _firecrown_xi(z, nz, theta_arcmin, hidden) + ind_fid_p, ind_fid_m, _ = _independent_camb_xi(z, nz, theta_arcmin, fid) + ind_hid_p, ind_hid_m, _ = _independent_camb_xi(z, nz, theta_arcmin, hidden) + + shift_fc_p = fc_hid_p - fc_fid_p + shift_ind_p = ind_hid_p - ind_fid_p + # the shift is O(few %) of ξ+, so compare relative to the shift magnitude + rel = np.abs(shift_fc_p - shift_ind_p) / np.abs(shift_fc_p).max() + assert rel.max() < 0.05, ( + f"blinding shift disagrees by {rel.max():.3%} across stacks" + ) + + +# --------------------------------------------------------------------------- # +# Fast smoke variant — the same machinery at coarse resolution, ~seconds +# --------------------------------------------------------------------------- # +def test_camb_ccl_crosscheck_smoke(): + """Fast smoke: both paths run and roughly agree at coarse resolution. + + Not the precision assertion (that is the ``slow`` test) — this guards that + the cross-check machinery imports, runs, and returns finite same-sign ξ± + that agree to a few percent, so a broken wiring is caught in the fast suite. + """ + cfg = bl.TheoryConfig() + z, nz = _gauss_nz(n=150) + theta_arcmin = np.geomspace(10.0, 100.0, 4) + + fc_xip, fc_xim = _firecrown_xi(z, nz, theta_arcmin, cfg) + ind_xip, ind_xim, As = _independent_camb_xi( + z, nz, theta_arcmin, cfg, n_ell=120, ell_max=30000, kmax=10.0 + ) + + assert np.all(np.isfinite(fc_xip)) and np.all(np.isfinite(ind_xip)) + assert np.all(fc_xip > 0) and np.all(ind_xip > 0) + assert 1e-9 < As < 3e-9 # σ8-matched A_s in a sane range + rel_p = np.abs(ind_xip - fc_xip) / np.abs(fc_xip) + assert rel_p.max() < 0.05, f"smoke ξ+ rel diff {rel_p.max():.3%} unexpectedly large" From b0e7470339e8b078184367f55fbae02c870ab9b7 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 07:45:48 +0200 Subject: [PATCH 2/2] docs(crosscheck): floor comments carry the observed numbers (0.21%/0.10%) Adversarial review measured the fixture floor at xi+ 0.213% / xi- 0.102% - xi- tighter than xi+, inverting the old comment's ~0.3%/~0.6% claim and its rationale. Comments now state the observed floor; the looser xi- tolerance stays as prudential margin, said honestly. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_016Ko89uZF84ez6HEDcTWbJe --- src/sp_validation/tests/test_camb_ccl_crosscheck.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/sp_validation/tests/test_camb_ccl_crosscheck.py b/src/sp_validation/tests/test_camb_ccl_crosscheck.py index 4a34698c..415c2153 100644 --- a/src/sp_validation/tests/test_camb_ccl_crosscheck.py +++ b/src/sp_validation/tests/test_camb_ccl_crosscheck.py @@ -30,8 +30,8 @@ convention is correctly reconciled. It does *not* independently validate CCL's Limber+Hankel machinery — both paths use ``ccl.angular_cl`` + ``ccl.correlation`` (a standalone FFTLog library is not in the environment), -so a common bug there would cancel. The residual we observe (~0.2–0.6%, -larger for ξ− which weights small scales) is dominated by P(k)-representation +so a common bug there would cancel. The residual we observe (0.21% for ξ+, +0.10% for ξ− on the test fixture) is dominated by P(k)-representation and ℓ-grid differences between two independent-but-both-CCL projections, not by cosmology; it is the realistic floor of cross-implementation agreement, and the tolerances sit comfortably above it. @@ -200,9 +200,9 @@ def test_sigma8_As_reconciliation(): # The cross-check: firecrown+CCL vs independent CAMB, σ8-matched # --------------------------------------------------------------------------- # # Observed agreement floor on this fixture (12 θ ∈ [5, 250] arcmin, HMCode2020): -# ξ+ ≈ 0.3%, ξ− ≈ 0.6% (converged in ℓ). Tolerances sit above the floor with -# headroom; ξ− is looser because it weights smaller scales, where the -# P(k)-representation and FFTLog differences between the two paths are largest. +# ξ+ ≈ 0.21%, ξ− ≈ 0.10% (converged in ℓ). Tolerances sit above the floor with +# headroom; ξ− keeps the looser tolerance as prudential margin for small-scale +# P(k)/FFTLog sensitivity on other fixtures, even though it measures tighter here. XIP_RTOL = 0.005 # 0.5% XIM_RTOL = 0.010 # 1.0%