diff --git a/.github/workflows/deploy-image.yml b/.github/workflows/deploy-image.yml index 192e87a6..13f0fafe 100644 --- a/.github/workflows/deploy-image.yml +++ b/.github/workflows/deploy-image.yml @@ -44,6 +44,12 @@ jobs: - name: Import smoke test run: docker run --rm ${{ steps.meta.outputs.tags }} python -c "import sp_validation" + # The fast suite doesn't import the blinding stack, so a broken + # firecrown/smokescreen install would otherwise ship green. Prove the + # image can actually load it (sacc + patched firecrown + smokescreen). + - name: Blinding-stack import smoke test + run: docker run --rm ${{ steps.meta.outputs.tags }} python -c "import sacc; import firecrown.likelihood; import smokescreen" + # Run the fast test suite against the freshly-built image *before* # pushing, so a failing suite blocks publication. The image carries the # full stack and the test files (COPY . + editable install), so this diff --git a/Dockerfile b/Dockerfile index ccc930bf..3caa129f 100644 --- a/Dockerfile +++ b/Dockerfile @@ -33,7 +33,28 @@ COPY pyproject.toml uv.lock /sp_validation/ RUN uv sync --frozen --inexact --no-install-project \ --extra test --extra glass --extra workflow -# Install sp_validation itself (editable) into the same venv; deps are already -# satisfied by the sync above. +# Full source in place before the blinding-extra install below: it is an +# editable install of the project itself and needs uv-overrides.txt and +# scripts/patch_firecrown.py from the tree. COPY . /sp_validation -RUN uv pip install --no-deps -e . + +# The [blinding] extra (SACC/Smokescreen blinding stack: firecrown + smokescreen) +# is not in uv.lock — firecrown is not on PyPI and declares conda-forge-only / +# unused sampler connectors as hard deps, so it needs the override file (see +# uv-overrides.txt for the full story). Installed as a separate pass on top of +# the locked sync. +RUN uv pip install --no-cache-dir --overrides uv-overrides.txt -e '.[blinding]' + +# Same uv gotcha as the cs_util upgrade above (astral-sh/uv #8410): if the base +# image / locked sync already carries a numpy that violates the [blinding] +# extra's `numpy<2.5` cap (firecrown 1.15.1 breaks on numpy 2.5 at import), the +# editable install won't move it. Request the bound explicitly so the image is +# deterministic either way; numpy 2.4.x is ABI-compatible with the compiled +# stack (verified: pyccl/camb/treecorr/healpy/pymaster + fast suite). +RUN uv pip install --no-cache-dir 'numpy>=2.2,<2.5' + +# firecrown is distributed for conda-forge (where NumCosmo always exists) and +# hits NumCosmo at import time in a pip env, on paths unrelated to our use. +# This patches the installed tree (surgical, pinned-version-checked, loud on +# mismatch) and verifies `import firecrown.likelihood; import smokescreen`. +RUN python scripts/patch_firecrown.py diff --git a/cosmo_inference/cosmosis_config/cosmosis_pipeline_A_ia.ini b/cosmo_inference/cosmosis_config/cosmosis_pipeline_A_ia.ini new file mode 100644 index 00000000..eb3ab166 --- /dev/null +++ b/cosmo_inference/cosmosis_config/cosmosis_pipeline_A_ia.ini @@ -0,0 +1,108 @@ +#parameters used elsewhere in this file +[DEFAULT] +COSMOSIS_DIR = /n23data1/n06data/lgoh/scratch/cosmosis-standard-library_lisa + + +[pipeline] +modules = consistency sample_S8 camb load_nz_fits photoz_bias linear_alignment projection add_intrinsic 2pt_shear shear_m_bias 2pt_like +likelihoods = 2pt_like +extra_output = cosmological_parameters/omega_lambda cosmological_parameters/S_8 cosmological_parameters/sigma_8 cosmological_parameters/omega_m +timing = T +debug = T + +[runtime] +sampler = polychord +verbosity = debug + +[polychord] +live_points = 192 +feedback = 3 +resume = T +base_dir = %(SCRATCH)s/polychord + +[test] + +[output] +format = text +lock = F + +[consistency] +file = %(COSMOSIS_DIR)s/utility/consistency/consistency_interface.py +verbose = F + +[sample_S8] +file = %(COSMOSIS_DIR)s/utility/sample_sigma8/sample_S8.py + +[camb] +file = %(COSMOSIS_DIR)s/boltzmann/camb/camb_interface.py +mode=power +lmax=2508 +feedback=0 +do_reionization=F +kmin=1e-5 +kmax=20.0 +nk=200 +zmax=5.0 +zmax_background=5.0 +nz_background=500 +halofit_version=mead2020_feedback +nonlinear=pk +neutrino_hierarchy=normal +kmax_extrapolate = 500.0 + +[load_nz_fits] +file = %(COSMOSIS_DIR)s/number_density/load_nz_fits/load_nz_fits.py +nz_file =%(FITS_FILE)s +data_sets = SOURCE + +[photoz_bias] +file = %(COSMOSIS_DIR)s/number_density/photoz_bias/photoz_bias.py +mode = additive +sample = nz_source +bias_section = nofz_shifts +interpolation = cubic +output_deltaz_section_name = delta_z_out + +[linear_alignment] +file = %(COSMOSIS_DIR)s/intrinsic_alignments/la_model/linear_alignments_interface_znla.py +method = bk_corrected + +[projection] +file = %(COSMOSIS_DIR)s/structure/projection/project_2d.py +ell_min_logspaced = 1.0 +ell_max_logspaced = 25000.0 +n_ell_logspaced = 400 +shear-shear = source-source +shear-intrinsic = source-source +intrinsic-intrinsic = source-source +get_kernel_peaks = F +verbose = F + +[add_intrinsic] +file = %(COSMOSIS_DIR)s/shear/add_intrinsic/add_intrinsic.py +shear-shear=T +position-shear=F +perbin=F + +[2pt_shear] +file = %(COSMOSIS_DIR)s/shear/cl_to_xi_nicaea/nicaea_interface.so +corr_type = 0 ; shear_cl -> shear_xi + +[shear_m_bias] +file = %(COSMOSIS_DIR)s/shear/shear_bias/shear_m_bias.py +m_per_bin = True +; Despite the parameter name, this can operate on xi as well as C_ell. +cl_section = shear_xi_plus shear_xi_minus +verbose = F + +[2pt_like] +file = %(COSMOSIS_DIR)s/likelihood/2pt/2pt_like.py +data_file=%(FITS_FILE)s +gaussian_covariance=F +covmat_name=COVMAT +cut_zeros=F +data_sets=XI_PLUS XI_MINUS +like_name=2pt_like + +angle_range_XI_PLUS_1_1= 10.0 200.0 +angle_range_XI_MINUS_1_1= 20.0 200.0 \ No newline at end of file diff --git a/cosmo_inference/cosmosis_config/cosmosis_pipeline_A_ia_sacc.ini b/cosmo_inference/cosmosis_config/cosmosis_pipeline_A_ia_sacc.ini new file mode 100644 index 00000000..5c31a476 --- /dev/null +++ b/cosmo_inference/cosmosis_config/cosmosis_pipeline_A_ia_sacc.ini @@ -0,0 +1,110 @@ +#parameters used elsewhere in this file +[DEFAULT] +COSMOSIS_DIR = /n23data1/n06data/lgoh/scratch/cosmosis-standard-library_lisa + + +[pipeline] +modules = consistency sample_S8 camb load_nz_sacc photoz_bias linear_alignment projection add_intrinsic 2pt_shear shear_m_bias sacc_like +likelihoods = 2pt_like +extra_output = cosmological_parameters/omega_lambda cosmological_parameters/S_8 cosmological_parameters/sigma_8 cosmological_parameters/omega_m +timing = T +debug = T + +[runtime] +sampler = polychord +verbosity = debug + +[polychord] +live_points = 192 +feedback = 3 +resume = T +base_dir = %(SCRATCH)s/polychord + +[test] + +[output] +format = text +lock = F + +[consistency] +file = %(COSMOSIS_DIR)s/utility/consistency/consistency_interface.py +verbose = F + +[sample_S8] +file = %(COSMOSIS_DIR)s/utility/sample_sigma8/sample_S8.py + +[camb] +file = %(COSMOSIS_DIR)s/boltzmann/camb/camb_interface.py +mode=power +lmax=2508 +feedback=0 +do_reionization=F +kmin=1e-5 +kmax=20.0 +nk=200 +zmax=5.0 +zmax_background=5.0 +nz_background=500 +halofit_version=mead2020_feedback +nonlinear=pk +neutrino_hierarchy=normal +kmax_extrapolate = 500.0 + +[load_nz_sacc] +file = %(COSMOSIS_DIR)s/number_density/load_nz_sacc/load_nz_sacc.py +nz_file = %(SACC_FILE)s +data_sets = source + +[photoz_bias] +file = %(COSMOSIS_DIR)s/number_density/photoz_bias/photoz_bias.py +mode = additive +sample = nz_source +bias_section = nofz_shifts +interpolation = cubic +output_deltaz_section_name = delta_z_out + +[linear_alignment] +file = %(COSMOSIS_DIR)s/intrinsic_alignments/la_model/linear_alignments_interface_znla.py +method = bk_corrected + +[projection] +file = %(COSMOSIS_DIR)s/structure/projection/project_2d.py +ell_min_logspaced = 1.0 +ell_max_logspaced = 25000.0 +n_ell_logspaced = 400 +shear-shear = source-source +shear-intrinsic = source-source +intrinsic-intrinsic = source-source +get_kernel_peaks = F +verbose = F + +[add_intrinsic] +file = %(COSMOSIS_DIR)s/shear/add_intrinsic/add_intrinsic.py +shear-shear=T +position-shear=F +perbin=F + +[2pt_shear] +file = %(COSMOSIS_DIR)s/shear/cl_to_xi_nicaea/nicaea_interface.so +corr_type = 0 ; shear_cl -> shear_xi + +[shear_m_bias] +file = %(COSMOSIS_DIR)s/shear/shear_bias/shear_m_bias.py +m_per_bin = True +; Despite the parameter name, this can operate on xi as well as C_ell. +cl_section = shear_xi_plus shear_xi_minus +verbose = F + +; Native SACC likelihood via the sp_validation shim (arcmin->rad + ordering +; guard over CosmoSIS's SaccClLikelihood). data_sets/angle ranges use the SACC +; grammar (full data-type names + tracer pairs). like_name=2pt_like keeps the +; block keys identical to the 2pt_like path so chain post-processing is unchanged. +[sacc_like] +file = %(SP_VALIDATION_MODULES)s/sacc_like_unions.py +csl_dir = %(COSMOSIS_DIR)s +data_file = %(SACC_FILE)s +data_sets = galaxy_shear_xi_plus galaxy_shear_xi_minus +like_name = 2pt_like + +angle_range_galaxy_shear_xi_plus_source_0_source_0 = 10.0 200.0 +angle_range_galaxy_shear_xi_minus_source_0_source_0 = 20.0 200.0 diff --git a/papers/cosmo_val/config/config.yaml b/papers/cosmo_val/config/config.yaml index 232513c7..16a47647 100644 --- a/papers/cosmo_val/config/config.yaml +++ b/papers/cosmo_val/config/config.yaml @@ -19,6 +19,11 @@ versions: [ # CosmologyValidation suite parameters (cosmo_val.py) # --------------------------------------------------------------------------- cosmo_val: + # Campaign run type: "data" or "mock". One switch: it gates Smokescreen + # blind-at-birth and is stamped as the SACC `type` of every part written + # (custody state at assembly — see blinding.assert_consistent_blind). + type: data + # CosmologyValidation constructor npatch: 100 theta_min: 1.0 @@ -120,11 +125,14 @@ harmonic: binning: powspace nbins: 32 -# Cosmological inference data-product locations (dormant subsystem). +# Cosmological inference data-product locations + tooling. inference: chains_dir: "/n09data/guerrini/output_chains" glass_mock_data_dir: "/n09data/guerrini/glass_mock_v1.4.6/results" glass_mock_chains_dir: "/n09data/guerrini/glass_mock_chains" + # CosmoSIS Standard Library checkout — fills COSMOSIS_DIR in the generated + # pipeline inis (the module `file =` paths and the sacc_like shim's csl_dir). + csl_dir: "/n23data1/n06data/lgoh/scratch/cosmosis-standard-library_lisa" cosebis: theta_min: 1.0 diff --git a/pyproject.toml b/pyproject.toml index 9c067178..3c81460e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,9 +9,10 @@ authors = [ ] license = {text = "MIT"} readme = "README.md" -# 3.12 floor: cosmo-numba (a hard dependency below) requires >=3.12, and the -# production container is Python 3.12 (shapepipe base image). Keeping the floor -# in sync with the container is what makes `uv lock` resolvable. +# 3.12 floor: Smokescreen 1.5.6 (and firecrown v1.15) require it, as does +# cosmo-numba (a hard dependency below); the container base (shapepipe:develop) +# is already python:3.12-slim-bookworm, so this aligns pyproject with the +# actual runtime and keeps `uv lock` resolvable. requires-python = ">=3.12" classifiers = [ "License :: OSI Approved :: MIT License", @@ -162,6 +163,40 @@ workflow = [ # published or resolvable (no public repo found) — left undeclared pending # its source. Same for `unions_wl` (scripts/check_footprint.py). ] +# Data-vector blinding (PRD #241 §3-§5): Smokescreen applies the Muir et al. +# shift d → d + t(hidden) − t(fid), with firecrown + CCL as the theory engine +# (only compute_theory_vector is used; sampling stays with CosmoSIS). The blind +# must be exactly recomputable from the seed at unblinding time, so the whole +# theory stack is pinned exactly, as a set. Smokescreen 1.5.6 + firecrown v1.15 +# both set the python floor (>=3.12). +# +# firecrown is not on PyPI and declares conda-forge-only / unused sampler +# connectors as hard deps, so installing this extra requires the dependency +# override file: `uv pip install --overrides uv-overrides.txt -e '.[blinding]'` +# (see uv-overrides.txt; the Dockerfile does this for the container). +# After installing this extra, run `python scripts/patch_firecrown.py` — it +# makes pip-installed firecrown importable without NumCosmo (conda-forge-only); +# see that script's docstring for the full story. +blinding = [ + "firecrown @ git+https://github.com/LSSTDESC/firecrown.git@v1.15.1", + "smokescreen==1.5.6", + "pyccl==3.3.4", + # firecrown 1.15.1 subclasses npt.NDArray (DataVector); numpy 2.5 turned + # npt.NDArray into a non-subclassable typing alias, breaking firecrown at + # import. firecrown's own env caps numpy<2.4; 2.4.3 is verified against + # the full compiled stack (pyccl/camb/treecorr/healpy/pymaster) + the + # sp_validation fast suite. + "numpy>=2.2,<2.5", +] +# Native SACC inference likelihood (PR 7, sp_validation.sacc_like_unions). The +# shim subclasses CosmoSIS's SaccClLikelihood; this pins the engine. Both engine +# module files are pure-python (no CAMB / compiled CSL), so the equality tests +# (tests/test_sacc_like.py) run with just this extra + a CosmoSIS Standard +# Library checkout, located via the CSL_DIR env var: +# git clone --depth 1 https://github.com/joezuntz/cosmosis-standard-library CSL +# CSL_DIR=/path/to/CSL pytest ... test_sacc_like.py +# Inert in CI (its Dockerfile installs only [test,glass,blinding]). +cosmosis = ["cosmosis>=3.25"] develop = ["sp_validation[test,docs]"] [tool.pytest.ini_options] diff --git a/scripts/blind_data_vector.py b/scripts/blind_data_vector.py new file mode 100644 index 00000000..1050d226 --- /dev/null +++ b/scripts/blind_data_vector.py @@ -0,0 +1,202 @@ +#!/usr/bin/env python3 + +"""Script blind_data_vector.py + +Per-part data-vector blinding with :mod:`sp_validation.blinding` +(Smokescreen-fork concealment, hash-commitment custody). + +``blind-init`` runs once per catalogue version: it draws an OS-entropy seed +(never printed, never written in plaintext), publishes a repo-committable +``commitment.json`` (the seed commitment + the config digest + the installed +Smokescreen fork's draw scheme), and encrypts the seed into a Fernet bundle. +Those three, not the seed alone, are what reproduces a blind (see +``blinding.draw_scheme``). ``blind-part`` blinds one intermediate part SACC (reporting ξ±, integration ξ±, +or pseudo-Cℓ) under that fixed state, escrows the true vector into a per-part +encrypted bundle beside the blinded output, and deletes the plaintext part. +``unblind`` verifies all three commitments and restores a true part +(bit-for-bit when the part's escrow bundle is beside it); it also works on the +assembled ``{version}.sacc`` (integration rows selected by the ``grid`` tag). +``verify`` is a cheap, seedless check that a blinded file matches a commitment +— seedless, but not environment-independent: it also compares both recorded +draw schemes against the installed fork, so it reports a problem on a machine +whose Smokescreen draws differently even when file and commitment agree. + +:Authors: Cail Daley + +Examples +-------- +Once per catalogue version:: + + blind_data_vector.py blind-init blinded/ + +Per intermediate part, at birth:: + + blind_data_vector.py blind-part parts/xi_integration.fits --blind-dir blinded/ + +Unblind one part:: + + blind_data_vector.py unblind parts/xi_integration_blinded.fits \\ + --blind-dir blinded/ -o parts/xi_integration.fits + +Verify:: + + blind_data_vector.py verify parts/xi_integration_blinded.fits \\ + blinded/commitment.json +""" + +import argparse +import json +import pathlib +import sys + +from sp_validation import blinding, sacc_io + + +def _config_from_args(args): + """A :class:`blinding.BlindingConfig` from optional CLI overrides.""" + overrides = {} + if args.s8_half_width is not None: + overrides["s8_half_width"] = args.s8_half_width + if args.omega_m_half_width is not None: + overrides["omega_m_half_width"] = args.omega_m_half_width + return blinding.BlindingConfig.from_overrides(overrides) + + +def _blind_init(args): + config = _config_from_args(args) + blind_dir = pathlib.Path(args.blind_dir) + blind_dir.mkdir(parents=True, exist_ok=True) + # Refuse before drawing anything: a blind is a one-shot custody event and + # silently overwriting a previous blind's state would destroy the record + # tying that blind to its seed. + clashes = [ + p + for p in blinding.init_paths(str(blind_dir)).values() + if pathlib.Path(p).exists() + ] + if clashes: + raise SystemExit( + "refusing to overwrite existing blind state:\n " + + "\n ".join(clashes) + + "\nPick a fresh --blind-dir (never overwrite a blind)." + ) + blinding.blind_init(str(blind_dir), config=config, label=args.label) + print( + "Commit the commitment JSON to the repo; keep the bundle + key safe " + "and separated (colocation in the blind dir is not at-rest protection)." + ) + + +def _blind_part(args): + blinding.blind_part( + args.part, + args.blind_dir, + config=_config_from_args(args), + keep_input=args.keep_input, + ) + + +def _unblind(args): + blinding.unblind_part( + args.blinded, + args.blind_dir, + args.output, + config=_config_from_args(args), + ) + + +def _verify(args): + """Seedless check: blinded-file metadata ↔ commitment JSON ↔ this install. + + No seed is read, so this cannot confirm that the blind is *subtractable* — + only that the file's three custody stamps agree with the commitment, and + that the recorded draw scheme is the one this install implements. That last + check makes the result environment-dependent by design: a machine carrying a + different Smokescreen could not unblind the file, so it reports a problem. + + Loads with ``allow_unblinded=True``: the whole job here is to report on a + file's custody state, including the state where the file is not concealed at + all, which the fail-closed loader would otherwise raise on before any + diagnostic could be assembled. + """ + s = sacc_io.load(args.blinded, allow_unblinded=True) + with open(args.commitment, encoding="utf-8") as f: + commitment = json.load(f) + problems = [] + if not s.metadata.get("concealed"): + problems.append("file is not marked concealed") + if s.metadata.get("blind_commitment") != commitment["seed_commitment"]: + problems.append("blind_commitment does not match the committed seed commitment") + if s.metadata.get("blind_config_digest") != commitment["config_digest"]: + problems.append("blind_config_digest does not match the committed digest") + # The draw scheme is checked two ways: file ↔ commitment, and the file's + # against the installed fork (see blinding.draw_scheme). + file_scheme = s.metadata.get("blind_draw_scheme") + committed = commitment.get("draw_scheme") + if file_scheme != committed: + problems.append( + f"blind_draw_scheme {file_scheme!r} does not match the committed " + f"draw_scheme {committed!r}" + ) + try: + blinding._assert_draw_scheme(file_scheme, "the blinded file") + except ValueError as exc: + problems.append(str(exc)) + if "seed_smokescreen" in s.metadata: + problems.append("PLAINTEXT SEED LEAKED into file metadata (seed_smokescreen)") + if problems: + raise SystemExit("verification FAILED:\n " + "\n ".join(problems)) + print( + f"OK: {args.blinded} matches {args.commitment} " + f"(blind {s.metadata.get('blind')!r})" + ) + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__.split("\n\n")[1]) + sub = parser.add_subparsers(dest="mode", required=True) + + for name in ("blind-init", "blind-part", "unblind"): + p = sub.add_parser(name) + p.add_argument("--s8-half-width", type=float, default=None) + p.add_argument("--omega-m-half-width", type=float, default=None) + if name == "blind-init": + p.add_argument( + "blind_dir", + help="directory for the blind's fixed state (commitment + " + "encrypted seed bundle)", + ) + p.add_argument("--label", default="A", help="blind label (default A)") + p.set_defaults(func=_blind_init) + elif name == "blind-part": + p.add_argument("part", help="intermediate part SACC file to blind") + p.add_argument( + "--blind-dir", required=True, help="blind-init state directory" + ) + p.add_argument( + "--keep-input", + action="store_true", + help="retain the plaintext input part (default: delete it " + "after blinding — the true vector is escrowed beside the " + "blinded output)", + ) + p.set_defaults(func=_blind_part) + else: + p.add_argument("blinded", help="blinded part (or assembled) SACC file") + p.add_argument( + "--blind-dir", required=True, help="blind-init state directory" + ) + p.add_argument("-o", "--output", required=True, help="output SACC path") + p.set_defaults(func=_unblind) + + p = sub.add_parser("verify") + p.add_argument("blinded", help="blinded SACC file") + p.add_argument("commitment", help="commitment JSON") + p.set_defaults(func=_verify) + + args = parser.parse_args(argv) + args.func(args) + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/patch_firecrown.py b/scripts/patch_firecrown.py new file mode 100644 index 00000000..7b528da1 --- /dev/null +++ b/scripts/patch_firecrown.py @@ -0,0 +1,225 @@ +"""Make pip-installed firecrown importable without NumCosmo. + +Run *inside* the target environment, after installing the ``[blinding]`` extra: + + python scripts/patch_firecrown.py + +Why this exists (PRD #241, PR 1): firecrown is the theory engine for +Smokescreen blinding — only ``compute_theory_vector`` on the SACC-read +cosmic-shear path is used. Upstream distributes firecrown via conda-forge, +where NumCosmo (a GObject-introspection C library, absent from PyPI) is always +present; in a pip/uv environment, firecrown 1.15.1 hits NumCosmo at *import +time* through two paths that have nothing to do with cosmic shear: + +1. ``firecrown/generators/__init__.py`` eagerly re-exports the LSST Y1/Y10 + predefined n(z) bin constants, defeating the lazy ``__getattr__`` that + ``_inferred_galaxy_zdist`` already provides — and computing those constants + imports NumCosmo. +2. ``firecrown/likelihood/__init__.py`` eagerly imports the cluster + likelihoods, which import ``crow`` (lsstdesc-crow), which subclasses a + NumCosmo C class at module load (``class CountsIntegralND(Ncm.IntegralND)``). + +This script (a) restores laziness in ``generators``, (b) makes the cluster +imports optional, and (c) installs a *loud* ``numcosmo_py`` shim so that any +genuine NumCosmo use raises immediately instead of being silently faked. +Everything is exact-string surgery against the pinned firecrown v1.15.1: if a +target string is missing (e.g. after a version bump), the script fails loudly +so the pin and the patch get reviewed together. Idempotent — safe to re-run. + +The right long-term fix is upstream (guarded/lazy imports in firecrown); until +then this file is the entire cost of staying pip-installable. +""" + +import importlib.metadata +import importlib.util +import subprocess +import sys +from pathlib import Path + +EXPECTED_FIRECROWN = "1.15.1" + +GENERATORS_OLD = """\ + # Lazy-loaded bins (via __getattr__) + Y1_LENS_BINS, + Y1_SOURCE_BINS, + Y10_LENS_BINS, + Y10_SOURCE_BINS, + LSST_Y1_LENS_HARMONIC_BIN_COLLECTION, + LSST_Y1_SOURCE_HARMONIC_BIN_COLLECTION, + LSST_Y10_LENS_HARMONIC_BIN_COLLECTION, + LSST_Y10_SOURCE_HARMONIC_BIN_COLLECTION, +) +""" + +GENERATORS_NEW = """\ +) + +# NOTE (sp_validation patch, scripts/patch_firecrown.py): the LSST Y1/Y10 +# predefined bin constants are computed lazily in _inferred_galaxy_zdist via a +# module-level __getattr__ that imports NumCosmo. Importing them EAGERLY here +# forced NumCosmo at `import firecrown.generators` (hence at +# `import firecrown.likelihood`), which pip cannot satisfy. Re-expose them +# lazily instead; the SACC-read cosmic-shear path never touches them. +_LAZY_BIN_NAMES = frozenset( + { + "Y1_LENS_BINS", + "Y1_SOURCE_BINS", + "Y10_LENS_BINS", + "Y10_SOURCE_BINS", + "LSST_Y1_LENS_HARMONIC_BIN_COLLECTION", + "LSST_Y1_SOURCE_HARMONIC_BIN_COLLECTION", + "LSST_Y10_LENS_HARMONIC_BIN_COLLECTION", + "LSST_Y10_SOURCE_HARMONIC_BIN_COLLECTION", + } +) + + +def __getattr__(name): + if name in _LAZY_BIN_NAMES: + from . import _inferred_galaxy_zdist as _z + + return getattr(_z, name) + raise AttributeError(f"module {__name__!r} has no attribute {name!r}") + +""" + +LIKELIHOOD_OLD = """\ +# Cluster statistics +from firecrown.likelihood._binned_cluster import BinnedCluster +from firecrown.likelihood._binned_cluster_number_counts import ( + BinnedClusterNumberCounts, +) +from firecrown.likelihood._binned_cluster_number_counts_shear import ( + BinnedClusterShearProfile, +) +""" + +LIKELIHOOD_NEW = """\ +# Cluster statistics. +# NOTE (sp_validation patch, scripts/patch_firecrown.py): the cluster +# likelihoods import `crow` (lsstdesc-crow), which subclasses NumCosmo C +# classes at module load. NumCosmo is conda-forge-only, so in a pip/uv env +# these imports fail. They are NOT on the cosmic-shear (TwoPoint/WeakLensing) +# path, so they become optional: without NumCosmo the cluster classes are +# unavailable but everything else loads. +try: + from firecrown.likelihood._binned_cluster import BinnedCluster + from firecrown.likelihood._binned_cluster_number_counts import ( + BinnedClusterNumberCounts, + ) + from firecrown.likelihood._binned_cluster_number_counts_shear import ( + BinnedClusterShearProfile, + ) +except (ImportError, RuntimeError, TypeError): # pragma: no cover + BinnedCluster = None # type: ignore[assignment,misc] + BinnedClusterNumberCounts = None # type: ignore[assignment,misc] + BinnedClusterShearProfile = None # type: ignore[assignment,misc] +""" + +SHIM = '''\ +"""Minimal loud shim for numcosmo_py (installed by sp_validation). + +NumCosmo is a GObject-introspection C library available only via conda-forge. +With the companion patches to firecrown (scripts/patch_firecrown.py), the +SACC-read cosmic-shear likelihood path never imports it; this shim provides +the import-time names so the patched package loads, and any genuine numerical +use of NumCosmo raises loudly rather than being silently faked. +""" + + +class _Missing: + def __init__(self, path="numcosmo_py"): + self._p = path + + def __getattr__(self, name): + return _Missing(f"{self._p}.{name}") + + def __call__(self, *a, **k): + raise RuntimeError( + f"{self._p} was called, but NumCosmo is not installed (conda-forge " + "only, not on PyPI). It is not needed for the SACC-read " + "cosmic-shear likelihood path." + ) + + def __getitem__(self, item): + return _Missing(f"{self._p}[...]") + + +Ncm = _Missing("numcosmo_py.Ncm") +Nc = _Missing("numcosmo_py.Nc") +GObject = _Missing("numcosmo_py.GObject") + + +def dict_to_var_dict(*a, **k): + raise RuntimeError("numcosmo_py.dict_to_var_dict unavailable (no NumCosmo)") + + +def var_dict_to_dict(*a, **k): + raise RuntimeError("numcosmo_py.var_dict_to_dict unavailable (no NumCosmo)") +''' + + +def patch_file(path: Path, old: str, new: str) -> str: + text = path.read_text() + if new in text: + return "already patched" + if old not in text: + sys.exit( + f"FATAL: expected text not found in {path}.\n" + "firecrown has probably been bumped past the pinned version this " + "patch targets — review scripts/patch_firecrown.py together with " + "the [blinding] pin in pyproject.toml." + ) + path.write_text(text.replace(old, new, 1)) + return "patched" + + +def main() -> None: + spec = importlib.util.find_spec("firecrown") + if spec is None or spec.origin is None: + sys.exit("FATAL: firecrown is not installed in this environment.") + pkg = Path(spec.origin).parent + + # Metadata, not `import firecrown` — pre-patch, importing is what's broken. + version = importlib.metadata.version("firecrown") + if version != EXPECTED_FIRECROWN: + sys.exit( + f"FATAL: firecrown {version} != expected {EXPECTED_FIRECROWN}; " + "review this patch against the new version before bumping " + "EXPECTED_FIRECROWN." + ) + + print( + "generators/__init__.py:", + patch_file(pkg / "generators" / "__init__.py", GENERATORS_OLD, GENERATORS_NEW), + ) + print( + "likelihood/__init__.py:", + patch_file(pkg / "likelihood" / "__init__.py", LIKELIHOOD_OLD, LIKELIHOOD_NEW), + ) + + # Loud numcosmo_py shim — only when no real NumCosmo is present. + if importlib.util.find_spec("numcosmo_py") is None: + shim_dir = pkg.parent / "numcosmo_py" + shim_dir.mkdir(exist_ok=True) + (shim_dir / "__init__.py").write_text(SHIM) + print("numcosmo_py shim: installed") + else: + print("numcosmo_py shim: skipped (numcosmo_py importable)") + + check = subprocess.run( + [ + sys.executable, + "-c", + "import firecrown.likelihood; import smokescreen", + ], + capture_output=True, + text=True, + ) + if check.returncode != 0: + sys.exit(f"FATAL: post-patch import check failed:\n{check.stderr}") + print("post-patch import check: firecrown.likelihood + smokescreen OK") + + +if __name__ == "__main__": + main() diff --git a/src/sp_validation/blinding.py b/src/sp_validation/blinding.py new file mode 100644 index 00000000..247b20e9 --- /dev/null +++ b/src/sp_validation/blinding.py @@ -0,0 +1,985 @@ +"""Blinding — conceal each intermediate data product behind a hidden cosmology. + +:Name: blinding.py + +:Description: Smokescreen blinding wiring, per part, at birth. Each blindable + intermediate SACC product — reporting ξ±, integration ξ±, pseudo-Cℓ — is shifted, + the moment the pipeline computes it, by a difference of theory vectors + between the fiducial cosmology and a *hidden* cosmology drawn inside a + fixed amplitude envelope, so no one can read S8 off the data until the + collaboration agrees to unblind (Muir et al. 2019: + ``d → d + t(hidden) − t(fiducial)``). Only blinded parts persist on disk. + + **The fork is the concealment engine.** The hidden cosmology is drawn by + the ``UNIONS-WL/Smokescreen`` fork's own CCL-native, per-key-independent + draw. Blinding is a vector operation, so this module calls the fork's + vector core directly — ``smokescreen.concealing_factor(fiducial_params, + shifts_dict, *, seed, theory_fn)``, which returns + ``t(hidden) − t(fiducial)`` and never sees a SACC. This module supplies + the two sp_validation-specific pieces: the three ``theory_fn`` backends + matching each part's row layout (reporting ξ±, integration ξ±, pseudo-Cℓ) + and the SACC handling around the returned factors. + + **Envelope calibration.** The blinding intent is an amplitude smear of a + chosen (S8, Ωm) box. The fork draws in CCL-native primitives, so + :meth:`BlindingConfig.shifts_dict` maps the intended box to + ``{sigma8, Omega_c}`` half-widths evaluated at the fiducial — + ``σ8 = S8/√(Ωm/0.3)``, ``Ω_c = Ωm − Ω_b − Ω_ν``. Under the installed draw + scheme the deltas depend only on ``(key, seed, shift_distr)``, so one seed + yields one hidden cosmology across all part passes and the blinded parts + are mutually consistent by construction. That "only" is a property of the + *draw scheme*, not of the seed: :func:`draw_scheme` records which scheme a + blind was drawn under and every custody gate refuses an install that + disagrees (see **Custody** below). + + **Derived statistics are born blinded.** COSEBIs and pure-E/B are never + touched by this module: the pipeline's own estimators + (:mod:`sp_validation.b_modes`) run in their normal downstream place on + the already-blinded integration ξ± part, so their outputs are born blinded. + Covariance and ρ/τ PSF diagnostics are never blinded — blinding hides + the vector, not the uncertainty, and the shift is pure E-mode so the + B-mode null tests stay honest under the blind. + + **Custody: hash commitment, no keyholder.** :func:`blind_init` runs once + per catalogue version: it draws an OS-entropy seed, publishes + the seed commitment (:func:`seed_commitment`), a canonical config digest, + and the installed fork's ``DRAW_SCHEME`` as a repo-committable + ``commitment.json``, and encrypts the seed into a Fernet bundle + (``smokescreen.encryption``) — the plaintext seed is never written. The + three together, not the seed alone, are what reproduces a blind (see + :func:`draw_scheme`). + + Each :func:`blind_part` call reads that fixed state, conceals one part, + escrows the part's true vector into its own encrypted bundle beside the + blinded output, and deletes the plaintext part. Every path that assembles + parts into the one-file product goes through + :func:`sp_validation.sacc_io.gather` — the production assembler is passed + *into* it rather than wrapped around it — and gather calls + :func:`assert_consistent_blind` to fail closed unless every blindable part + carries the identical ``blind_commitment``, ``blind_config_digest`` and + ``blind_draw_scheme``. :func:`unblind_part` verifies all three against the + commitment *before* subtracting anything, then restores the true part. + + Custody has no back door, including for its own history: a blind carrying no + scheme record is refused by :func:`_assert_draw_scheme` ahead of every seed + read, with no CLI override — none exists, the scheme having been bound + before the first real blind. +""" + +import dataclasses +import hashlib +import json +import os +import secrets +import warnings + +import numpy as np + +from . import sacc_io +from .blinding_theory import TheoryConfig, cl_ee, xi_ccl, xi_ell_grid + + +# --------------------------------------------------------------------------- # +# Configuration surface — the blinding envelope +# --------------------------------------------------------------------------- # +@dataclasses.dataclass(frozen=True) +class BlindingConfig: + """Blinding envelope and fiducial. + + The hidden cosmology is drawn (by the fork) uniformly and independently + per key inside ``shifts_dict()``'s half-widths about the fiducial. The + half-widths are the deliberate, configurable size of the blind — config, + not code; the group may resize the envelope. ``theory`` carries the + fiducial :class:`TheoryConfig` whose defaults *are* the blinding + fiducial. + """ + + s8_half_width: float = 0.075 + omega_m_half_width: float = 0.1 + theory: TheoryConfig = dataclasses.field(default_factory=TheoryConfig) + + def shifts_dict(self): + """The (S8, Ωm) envelope as CCL-native ``{sigma8, Omega_c}`` half-widths. + + Evaluated at the fiducial: a ΔS8 half-width maps to + ``ΔS8/√(Ωm_fid/0.3)`` in σ8 (at fixed Ωm), and a ΔΩm half-width maps + one-to-one to Ω_c (Ω_b and Ω_ν are fixed). Exact enough for a + blinding smear — the target is a characteristic amplitude, not a + precise (S8, Ωm) posterior. Under the installed draw scheme the fork + draws each key independently as ``U(fid − h, fid + h)``, so adding or + resizing one key never moves another; :func:`draw_scheme` is what + makes that independence a guarantee rather than a hope, by refusing an + install whose draw semantics differ from the blind's. + """ + return { + "sigma8": self.s8_half_width / np.sqrt(self.theory.Omega_m / 0.3), + "Omega_c": self.omega_m_half_width, + } + + def config_digest(self): + """sha256 of a canonical serialization of the full blinding config. + + Binds the envelope half-widths, the complete fiducial + :class:`TheoryConfig` (cosmology + the two ``halofit_version`` + tokens + IA fields) into one digest: JSON with sorted keys over the + full ordered field set. Every ``float``-declared field is coerced + through :func:`float` before serialization, so the digest depends on + the numeric *value*, not on whether a config wrote ``-1`` or ``-1.0`` + — an int-vs-float literal difference (routine when configs come from + YAML/CLI/humans) can no longer split one physical cosmology into two + digests and deny a legitimate unblind. Python's ``json`` then emits + each float via its shortest round-trip ``repr``, so two runs of the + same config produce byte-identical digests. Checked (with + the seed commitment) at unblind, so a wrong envelope or a mismatched + P(k) recipe cannot silently subtract a wrong shift. + """ + payload = { + "s8_half_width": float(self.s8_half_width), + "omega_m_half_width": float(self.omega_m_half_width), + "theory": { + f.name: ( + float(getattr(self.theory, f.name)) + if f.type in (float, "float") + else getattr(self.theory, f.name) + ) + for f in dataclasses.fields(self.theory) + }, + } + return hashlib.sha256( + json.dumps(payload, sort_keys=True).encode("utf-8") + ).hexdigest() + + @classmethod + def from_overrides(cls, overrides): + """Build from a mapping of field overrides (fail loud on unknown keys). + + ``theory`` may be given as a :class:`TheoryConfig` or a mapping of + TheoryConfig overrides, mirroring :meth:`TheoryConfig.from_overrides`. + """ + by_name = {f.name: f for f in dataclasses.fields(cls)} + unknown = set(overrides) - set(by_name) + if unknown: + raise ValueError( + f"unknown BlindingConfig fields {sorted(unknown)}; " + f"valid fields are {sorted(by_name)}" + ) + overrides = { + name: (float(v) if by_name[name].type in (float, "float") else v) + for name, v in overrides.items() + } + theory = overrides.get("theory") + if theory is not None and not isinstance(theory, TheoryConfig): + overrides["theory"] = TheoryConfig.from_overrides(dict(theory)) + return cls(**overrides) + + +# --------------------------------------------------------------------------- # +# Custody primitives +# --------------------------------------------------------------------------- # +def seed_commitment(seed): + """Public commitment for a seed: the fork's domain-separated sha256 digest. + + Safe to publish and commit to the repo — it ties a blinded file to its + blind without revealing the seed, and lets unblind refuse a wrong seed. + + This delegates to :func:`smokescreen.seed_commitment` so the commitment has + exactly one definition across the blinding stack. That function hashes + ``smokescreen.COMMITMENT_DOMAIN + str(seed)``, *not* the bare seed, and the + domain prefix is load-bearing: the fork derives the RNG base seed from the + bare sha256 of the same string, so an undomained commitment would publish + that base seed verbatim in its first 16 hex characters and the blind would + be recoverable from the artifact meant to protect it. + """ + from smokescreen import seed_commitment as _fork_seed_commitment + + return _fork_seed_commitment(seed) + + +def draw_scheme(): + """The installed Smokescreen fork's shift-draw semantics version. + + ``smokescreen.DRAW_SCHEME`` is the fork's own version number for *how* a + seed becomes a set of parameter deltas — scheme 1 is upstream DESC's one + global RNG stream consumed over the sorted keys, scheme 2 (this fork) a + per-key RNG derived from ``(seed, key)``. Two installs that agree on + ``(seed, config)`` but disagree on this number draw *different* hidden + cosmologies from the same inputs. + + That is the one blinding failure with no loud symptom: a blind added under + one scheme and subtracted under another leaves a residual smooth + cosmological shift in the "unblinded" vector, which every hash, digest and + escrow check would still pass. So the scheme is bound into the blind's + custody state — ``commitment.json`` and every blinded file's + ``blind_draw_scheme`` — and re-checked against this function wherever a + shift is drawn or subtracted (:func:`_assert_draw_scheme`). + """ + from smokescreen import DRAW_SCHEME + + return int(DRAW_SCHEME) + + +def _assert_draw_scheme(recorded, what): + """Fail closed unless ``recorded`` is the installed fork's draw scheme. + + ``what`` names the surface the scheme was read from, for the message. + A missing record (``None``) is a failure, not a pass: a blind whose scheme + is unknown cannot be shown to be reproducible by this install (see + :func:`draw_scheme`). + """ + installed = draw_scheme() + if recorded is None: + raise ValueError( + f"{what} carries no draw-scheme record — refusing to proceed. " + f"It predates draw-scheme binding, so there is no way to tell " + f"whether the installed Smokescreen (DRAW_SCHEME={installed}) " + f"reproduces the shift it was blinded with." + ) + if int(recorded) != installed: + raise ValueError( + f"{what} was drawn under Smokescreen DRAW_SCHEME={int(recorded)} " + f"but the installed fork implements DRAW_SCHEME={installed} — " + f"refusing to proceed. The same seed draws a different hidden " + f"cosmology under a different scheme, so this install would " + f"subtract the wrong shift and pass every other check. Install " + f"the Smokescreen the blind was made with." + ) + + +def hidden_params(seed, config): + """The hidden CCL parameter point the fork realizes for ``(seed, config)``. + + Re-runs the fork's own draw (``smokescreen.param_shifts.draw_param_shifts`` + — per-key-independent, local RNG) on the calibrated envelope and overlays + the deltas on the fiducial, exactly as + ``smokescreen.concealing_factor`` does internally. Deterministic under a + fixed draw scheme: same ``(seed, config)`` ⇒ same hidden point on any + install whose :func:`draw_scheme` matches — the reproducibility contract + unblinding relies on, and what makes every part share one hidden cosmology + under one seed. This is an introspection helper (what *was* the hidden + cosmology, once revealed); the blinding path itself never calls it, and so + does not check the scheme here — every gate that acts on a shift does. + """ + from smokescreen.param_shifts import draw_param_shifts + + deltas = draw_param_shifts(config.shifts_dict(), seed) + params = dict(config.theory.ccl_params()) + for key, delta in deltas.items(): + params[key] += delta + return params + + +# --------------------------------------------------------------------------- # +# Block discovery on a SACC (a standalone part, or the assembled file) +# --------------------------------------------------------------------------- # +def source_bins(s): + """Sorted source-bin indices present in ``s`` (from ``source_i`` tracers).""" + return sorted( + int(name.split("_", 1)[1]) for name in s.tracers if name.startswith("source_") + ) + + +def xi_pairs(s, grid): + """Unordered source-bin pairs ``(i ≤ j)`` carrying ξ+ on ``grid``.""" + bins = source_bins(s) + return [ + (i, j) + for a, i in enumerate(bins) + for j in bins[a:] + if len(s.indices(sacc_io.XI_PLUS, sacc_io._pair((i, j)), grid=grid)) + ] + + +def cl_pairs(s): + """Source-bin pairs ``(i ≤ j)`` carrying pseudo-Cℓ_EE.""" + bins = source_bins(s) + return [ + (i, j) + for a, i in enumerate(bins) + for j in bins[a:] + if len(s.indices(sacc_io.CL_EE, sacc_io._pair((i, j)))) + ] + + +def _xi_indices(s, grid): + """Row indices of the ξ± block on ``grid`` (ascending).""" + return np.sort( + np.concatenate( + [ + s.indices(sacc_io.XI_PLUS, grid=grid), + s.indices(sacc_io.XI_MINUS, grid=grid), + ] + ) + ).astype(int) + + +def _cl_ee_indices(s): + """Row indices of the pseudo-Cℓ_EE block (ascending). + + Only EE: a pure E-mode cosmology shift leaves BB and EB identically + zero, so those blocks are never extracted, never concealed. + """ + return np.sort(s.indices(sacc_io.CL_EE)).astype(int) + + +def _pair_nz(s, i, j): + """The two per-bin n(z) for pair ``(i, j)`` as ``((z_i, n_i), (z_j, n_j))``.""" + z_i, n_i = sacc_io.get_nz(s, i) + z_j, n_j = sacc_io.get_nz(s, j) + return (np.asarray(z_i), np.asarray(n_i)), (np.asarray(z_j), np.asarray(n_j)) + + +# --------------------------------------------------------------------------- # +# The three theory backends — callables aligned to a sub-SACC block's rows +# --------------------------------------------------------------------------- # +def xi_theory_fn(block, theory, grid): + """``theory_fn`` for a ξ± sub-SACC block (reporting or integration grid). + + Reads each bin's n(z) directly from the block's own tracers and lays the + output out to match the block's SACC rows element-for-element: for every + pair, ξ± is computed at that pair's stored θ (the ``theta`` tag, + arcmin) and scattered to the rows ``block.indices`` reports — the + block's own row order, never an assumed pairing. IA enters from + ``theory``'s NLA fields, identically at every parameter point. + """ + pairs = xi_pairs(block, grid) + layout = [] + for i, j in pairs: + tr = sacc_io._pair((i, j)) + idx_p = block.indices(sacc_io.XI_PLUS, tr, grid=grid) + idx_m = block.indices(sacc_io.XI_MINUS, tr, grid=grid) + theta = sacc_io._tag(block, sacc_io.XI_PLUS, tr, "theta", grid=grid) + layout.append(((i, j), idx_p, idx_m, np.asarray(theta, dtype=float))) + ell = xi_ell_grid() + + def theory_fn(params): + out = np.full(len(block.mean), np.nan) + for (i, j), idx_p, idx_m, theta in layout: + nz_i, nz_j = _pair_nz(block, i, j) + xip, xim = xi_ccl(params, theory, nz_i, nz_j, theta, ell) + out[idx_p] = xip + out[idx_m] = xim + return out + + return theory_fn + + +def cl_theory_fn(block, theory): + """``theory_fn`` for the pseudo-Cℓ_EE sub-SACC block. + + Per pair: theory Cℓ_EE on the stored ``BandpowerWindows`` support, + binned by the same window matrix the measurement used (``W @ Cℓ_EE``), + scattered to the block's own rows. The concealing factor the fork forms + from this is ``W @ ΔCℓ_EE`` — the shift lands in the measured + bandpowers; ΔBB = ΔEB ≡ 0 by construction (pure E-mode shift), so those + blocks are simply not part of this backend's rows. + """ + layout = [] + for i, j in cl_pairs(block): + tr = sacc_io._pair((i, j)) + idx = block.indices(sacc_io.CL_EE, tr) + window = block.get_bandpower_windows(idx) + layout.append( + ( + (i, j), + np.asarray(idx), + np.asarray(window.values, dtype=float), # (n_ell,) + np.asarray(window.weight, dtype=float), # (n_ell, n_bp) + ) + ) + + def theory_fn(params): + out = np.full(len(block.mean), np.nan) + for (i, j), idx, w_ell, w_mat in layout: + nz_i, nz_j = _pair_nz(block, i, j) + out[idx] = w_mat.T @ cl_ee(params, theory, nz_i, nz_j, w_ell) + return out + + return theory_fn + + +# --------------------------------------------------------------------------- # +# The concealing factor, per block +# --------------------------------------------------------------------------- # +def _blindable_blocks(s): + """The blindable blocks of a SACC as ``(name, indices, factory)``. + + Works identically on a standalone part (which carries exactly one block) + and on the assembled one-file product (whose integration rows are selected by + the ``grid`` tag — the layout contract's per-block tag selection). + ``indices`` are each block's recorded row indices (ascending); ``factory`` + builds the matching ``theory_fn`` from the SACC those indices point into. + Blocks absent from the file are simply not listed. + """ + blocks = [] + for grid in ("reporting", "integration"): + idx = _xi_indices(s, grid) + if len(idx): + blocks.append( + ( + f"{grid} ξ±", + idx, + lambda sub, theory, grid=grid: xi_theory_fn(sub, theory, grid), + ) + ) + idx = _cl_ee_indices(s) + if len(idx): + blocks.append(("pseudo-Cℓ_EE", idx, cl_theory_fn)) + return blocks + + +def _concealing_factor(s, indices, factory, config, seed): + """The fork-computed additive concealing factor for one block of ``s``. + + Blinding is a vector operation, so this goes straight to the fork's vector + core: ``smokescreen.concealing_factor`` draws the hidden deltas from + ``seed``, overlays them on the fiducial, evaluates the block's + ``theory_fn`` at both points and differences them. No data vector and no + SACC reach the fork — nothing is carved out of ``s``, and ``s`` itself is + not modified. + + The factory reads its layout off ``s`` directly, so the returned vector is + full-length: a block's ``theory_fn`` fills only its own rows and leaves + every other row NaN. Slicing to ``indices`` drops those NaNs by + construction; the finite check then proves the converse — that the + factory filled *all* of this block's rows. A row the block claims but the + factory cannot cover (a pair carrying ξ− without ξ+, say) would otherwise + write NaN into the data vector, silently. + + Both :func:`blind_sacc` and :func:`unblind_sacc` reach the fork through + this one function, so the added and the subtracted shift cannot drift + apart. + + Returns + ------- + np.ndarray + ``t(hidden) − t(fiducial)``, aligned to ``indices``. + """ + from smokescreen import concealing_factor + + full = np.asarray( + concealing_factor( + config.theory.ccl_params(), + config.shifts_dict(), + seed=seed, + theory_fn=factory(s, config.theory), + factor_type="add", + ), + dtype=float, + ) + factor = full[indices] + if not np.all(np.isfinite(factor)): + raise ValueError( + f"the theory backend left {int(np.sum(~np.isfinite(factor)))} of " + f"{len(indices)} blindable rows unfilled — refusing to apply the " + "concealing factor (these rows would be shifted by NaN). The " + "block's row layout is not fully covered by its theory_fn." + ) + return factor + + +def _set_values(s, indices, values): + """Overwrite ``s.data[i].value`` for ``indices`` with ``values`` (aligned).""" + for i, v in zip(indices, values): + s.data[int(i)].value = float(v) + + +def _concealed(s): + """Whether ``s`` is already a blinded file (its ``concealed`` mark is set).""" + return bool(s.metadata.get("concealed")) + + +def blind_sacc(part, seed, config=None, label="A", log=print): + """Return a blinded copy of a part SACC (covariance and tags untouched). + + Per blindable block present (a standalone part carries exactly one — + reporting ξ±, integration ξ±, or pseudo-Cℓ_EE): ask the fork for the + block's concealing factor (:func:`_concealing_factor`) and add it at the + block's recorded indices. Row order, tags, n(z) and covariance are + untouched — only ``value`` changes, and only on blindable rows. + Provenance is stamped and any leaked seed key stripped. A file with no + blindable block (e.g. a ρ/τ diagnostic part) is refused loudly — it should + never see a blind call. + """ + config = config or BlindingConfig() + if _concealed(part): + raise ValueError("already concealed — unblind first") + blocks = _blindable_blocks(part) + if not blocks: + raise ValueError( + "no blindable block (reporting/integration ξ± or pseudo-Cℓ_EE) in this SACC " + "— ρ/τ diagnostic parts are never blinded" + ) + + blinded = part.copy() + for name, indices, factory in blocks: + factor = _concealing_factor(part, indices, factory, config, seed) + _set_values(blinded, indices, np.asarray(blinded.mean)[indices] + factor) + log(f"[blind] {name}: shifted {len(indices)} points via Smokescreen fork") + + _stamp_provenance(blinded, seed_commitment(seed), label, config.config_digest()) + return blinded + + +def unblind_sacc(blinded, seed, config=None, log=print): + """Recover the true part SACC from a blinded one + the revealed ``seed``. + + Verifies the stamped draw scheme against the installed fork, then + the seed commitment and the config digest against the stamped metadata (loud + failure on any of the three — verification precedes subtraction), then + recomputes each block's shift through :func:`_concealing_factor`, the same + call :func:`blind_sacc` added it with, and subtracts it. Works on a + standalone part or on the assembled file (integration rows selected by the + ``grid`` tag). Derived statistics, if present (assembled file), are *not* + recomputed here — the pipeline's own estimators re-derive them from the + unblinded integration ξ±. + """ + config = config or BlindingConfig() + if not _concealed(blinded): + raise ValueError("file is not concealed — nothing to unblind") + _assert_draw_scheme(blinded.metadata.get("blind_draw_scheme"), "this blinded file") + if seed_commitment(seed) != blinded.metadata["blind_commitment"]: + raise ValueError( + "seed does not match blind_commitment — refusing to unblind " + "(a wrong seed would silently produce a wrong data vector)" + ) + if config.config_digest() != blinded.metadata["blind_config_digest"]: + raise ValueError( + "blinding config does not match blind_config_digest — refusing to " + "unblind (this config would subtract a different shift than was " + "added)" + ) + + part = blinded.copy() + for name, indices, factory in _blindable_blocks(blinded): + factor = _concealing_factor(blinded, indices, factory, config, seed) + _set_values(part, indices, np.asarray(part.mean)[indices] - factor) + log(f"[unblind] {name}: subtracted {len(indices)} shifts") + + for key in ( + "concealed", + "blind", + "blind_commitment", + "blind_config_digest", + "blind_draw_scheme", + ): + part.metadata.pop(key, None) + return part + + +def _stamp_provenance(s, commitment, label, config_digest): + """Stamp the blind's public provenance; strip any leaked seed. + + ``blind_commitment`` (sha256 of the seed) ties the file to its blind + without revealing it; ``blind_config_digest`` pins the envelope + + fiducial the shift was drawn against; ``blind_draw_scheme`` pins the + Smokescreen draw semantics that turned the seed into that shift (see + :func:`draw_scheme`); ``concealed``/``blind`` mark the file. The three + together are what an unblind must reproduce. ``seed_smokescreen`` — the + raw seed upstream Smokescreen's writer would stamp — is popped + defensively: the seed must never ride a kept file. + + The stamped scheme is always the installed fork's: every caller has already + checked the blind's recorded scheme against it (:func:`_assert_draw_scheme`), + so the two agree by the time this runs. + """ + s.metadata.pop("seed_smokescreen", None) + s.metadata["concealed"] = True + s.metadata["blind"] = label + s.metadata["blind_commitment"] = commitment + s.metadata["blind_config_digest"] = config_digest + s.metadata["blind_draw_scheme"] = draw_scheme() + + +def stamp_concealed_passthrough(s, commitment_path): + """Stamp a part concealed under an existing blind, values untouched. + + A born-blinded derived statistic (COSEBIs / pure-E/B) has its E-mode vector + re-derived from the *already-blinded* integration ξ± before this call, so its + values are blind by construction and only the provenance stamp is missing. + ρ/τ carries no cosmological vector at all — the stamp merely clears it for + assembly under the blind (values unchanged either way). Both cases need the + custody keys the fail-closed load gate and :func:`assert_consistent_blind` + check: ``concealed``, ``blind`` (label), ``blind_commitment`` + (= ``seed_commitment``), ``blind_config_digest``, ``blind_draw_scheme``. This + reads those from the version's ``commitment.json`` (written by + :func:`blind_init`) and stamps them via :func:`_stamp_provenance`, so a + pass-through part shares the exact same custody state as the blinded + ξ±/pseudo-Cℓ parts. The committed draw scheme is checked against the + installed fork first: a pass-through part is blind because it was derived + from parts blinded under that scheme, so stamping it from an install that + draws differently would mint a custody claim this install cannot honour. + + Unlike :func:`blind_sacc`, this shifts nothing and does not require a + blindable block — it is the seam for parts blinded (or made blind-irrelevant) + upstream of the SACC writer. + """ + with open(commitment_path, encoding="utf-8") as f: + commitment = json.load(f) + _assert_draw_scheme( + commitment.get("draw_scheme"), f"the blind at {commitment_path}" + ) + _stamp_provenance( + s, + commitment["seed_commitment"], + commitment["label"], + commitment["config_digest"], + ) + return s + + +# --------------------------------------------------------------------------- # +# Assembly-time custody: one blind across all parts +# --------------------------------------------------------------------------- # +def assert_consistent_blind(parts): + """Assert every blindable part shares one blind; return the shared stamp. + + Custody logic for the terminal :func:`sp_validation.sacc_io.gather`: a + part is *blindable* if it carries a blindable block (ξ± or pseudo-Cℓ_EE); + ρ/τ diagnostic and covariance-only parts are exempt. Fails closed — + ``ValueError`` — if blinded and plaintext blindable parts are mixed, or + if two parts carry different + ``blind_commitment``/``blind_config_digest``/``blind_draw_scheme`` (they + were blinded under different seeds, configs or draw semantics and must + never be combined), or if the shared draw scheme is not the installed + fork's (this install could not unblind what it is about to assemble). + The consistency key is ``(commitment, digest, scheme)`` — the + ``blind`` *label* is informational provenance, not custody state, so + parts blinded under one seed+config but tagged with different ``--label`` + values assemble cleanly (a distinct warning is logged, not a failure). + + **Unconcealed parts must be declared mocks** (PRD §4, "Mocks vs data"): + when no blindable part is concealed, every blindable part's metadata must + carry ``type == "mock"`` — the tag :mod:`sp_validation.sacc_io` stamps + when the data vector is computed. An unconcealed ``type == "data"`` part, + or one missing the tag, fails assembly closed: skipping the blind can + never silently expose real data. A concealed+plaintext mix fails closed + regardless of ``type`` (see above). A fully-mock plaintext assembly + returns ``None``. + + Returns + ------- + dict or None + The shared blind metadata (``concealed``, ``blind``, + ``blind_commitment``, ``blind_config_digest``, ``blind_draw_scheme``) + for the gather to stamp on the assembled file, or ``None`` when + nothing is blinded. + """ + blindable = [p for p in parts if _blindable_blocks(p)] + concealed = [p for p in blindable if _concealed(p)] + if not concealed: + # `.get` is deliberate here: a missing `type` tag must count as + # not-a-mock and fail closed, not KeyError with less context. + exposed = sorted( + {str(p.metadata.get("type", "")) for p in blindable} - {"mock"} + ) + if exposed: + raise ValueError( + f"unconcealed blindable parts with type {exposed} in assembly " + "— only parts declared `type: mock` may assemble without a " + "blind (an unconcealed data part exposes the real vector)" + ) + return None + if len(concealed) != len(blindable): + raise ValueError( + f"blinded and plaintext blindable parts mixed in one assembly " + f"({len(concealed)} of {len(blindable)} blinded) — refusing to " + "combine (a plaintext part beside blinded ones leaks the shift)" + ) + # Custody state is (commitment, digest, scheme) — the label is provenance. + # `.get` on the scheme so a part predating scheme binding reads as None and + # fails at _assert_draw_scheme with its explanation, not with a KeyError. + stamps = { + ( + p.metadata["blind_commitment"], + p.metadata["blind_config_digest"], + p.metadata.get("blind_draw_scheme"), + ) + for p in concealed + } + if len(stamps) != 1: + raise ValueError( + "parts carry different blind commitments — they were blinded " + "under different seeds, configs or draw schemes and must never be " + "combined: " + + "; ".join(f"({c[:12]}…, {d[:12]}…, scheme {v})" for c, d, v in stamps) + ) + ((commitment, digest, scheme),) = stamps + _assert_draw_scheme(scheme, "the blind these parts share") + labels = sorted({p.metadata["blind"] for p in concealed}) + if len(labels) != 1: + warnings.warn( + f"blindable parts share one blind (commitment {commitment[:12]}…, " + f"config {digest[:12]}…) but carry different labels {labels} — " + "assembling anyway; the label is provenance, not custody state. " + f"Stamping the assembled file with label {labels[0]!r}." + ) + return { + "concealed": True, + "blind": labels[0], + "blind_commitment": commitment, + "blind_config_digest": digest, + "blind_draw_scheme": int(scheme), + } + + +# --------------------------------------------------------------------------- # +# File-level custody: blind-init / blind-part / unblind +# --------------------------------------------------------------------------- # +def init_paths(blind_dir): + """The fixed custody state written by :func:`blind_init` in ``blind_dir``.""" + return { + "commitment": os.path.join(blind_dir, "commitment.json"), + "bundle": os.path.join(blind_dir, "blind_seed.encrpt"), + "key": os.path.join(blind_dir, "blind_seed.key"), + } + + +def part_paths(part_path): + """Blinded-output and escrow-bundle paths beside a part file.""" + stem, ext = os.path.splitext(part_path) + return { + "blinded": f"{stem}_blinded{ext or '.fits'}", + "escrow": f"{stem}_escrow.encrpt", + "escrow_key": f"{stem}_escrow.key", + } + + +def blind_init(blind_dir, config=None, label="A", log=print): + """Fix the blind for one catalogue version: seed, commitment, seed bundle. + + Runs once per catalogue version: + + 1. Draw an OS-entropy seed (never written in plaintext, never returned). + 2. Write ``commitment.json`` (repo-committable): the seed commitment + the + canonical config digest + the installed fork's draw scheme + the blind + label. Those first three are the full reproducibility statement — see + :func:`draw_scheme` for why the seed and config alone are not. + 3. Encrypt the seed into a Fernet bundle (``smokescreen.encryption``); + the temporary plaintext is deleted by the encryptor. + + These outputs are the blind's fixed state; every :func:`blind_part` and + :func:`unblind_part` call reads them. + + Custody caveat: the bundle and its Fernet key land in the *same* + ``blind_dir``. Fernet is only as protective as the key's separation — + anyone with both files can decrypt the seed. Keep the key out-of-band; + colocation is convenience, not at-rest protection. + + Returns + ------- + dict + Paths written: ``commitment``, ``bundle``, ``key``. + """ + config = config or BlindingConfig() + paths = init_paths(blind_dir) + for path in paths.values(): + if os.path.exists(path): + raise FileExistsError( + f"refusing to overwrite existing blind state {path} — a blind " + "is a one-shot custody event; choose another directory" + ) + + seed = secrets.token_hex(16) + commitment = { + "label": label, + "seed_commitment": seed_commitment(seed), + "config_digest": config.config_digest(), + "draw_scheme": draw_scheme(), + } + with open(paths["commitment"], "w", encoding="utf-8") as f: + json.dump(commitment, f, indent=2, sort_keys=True) + _write_encrypted_json(paths["bundle"], {"label": label, "seed": seed}) + + log(f"[blind-init] commitment (repo-committable): {paths['commitment']}") + log(f"[blind-init] encrypted seed bundle + key: {paths['bundle']}, {paths['key']}") + log( + "[blind-init] custody: keep the bundle key out-of-band from the bundle " + "(colocation in the blind dir is not at-rest protection)" + ) + return paths + + +def _read_seed(blind_dir, config): + """Decrypt the seed bundle and verify it against the commitment. + + The seed commitment, the config digest, and the committed draw scheme (see + :func:`draw_scheme`) are all checked before the seed is handed to any caller + — a tampered bundle, a drifted config or a Smokescreen that draws + differently from the one that fixed this blind all fail loud here, whether + the caller is about to blind or to unblind. + + Returns + ------- + tuple + ``(seed, commitment_dict)``. + """ + paths = init_paths(blind_dir) + bundle = _read_encrypted_json(paths["bundle"], paths["key"]) + with open(paths["commitment"], encoding="utf-8") as f: + commitment = json.load(f) + _assert_draw_scheme(commitment.get("draw_scheme"), f"the blind in {blind_dir}") + if seed_commitment(bundle["seed"]) != commitment["seed_commitment"]: + raise ValueError( + "bundle seed does not match the committed seed commitment — refusing " + "to proceed" + ) + if config.config_digest() != commitment["config_digest"]: + raise ValueError( + "blinding config does not match the committed config digest — " + "refusing to proceed (a wrong envelope or P(k) recipe would " + "silently produce a wrong shift)" + ) + return bundle["seed"], commitment + + +def blind_part(part_path, blind_dir, config=None, keep_input=False, log=print): + """Blind one intermediate part SACC at birth, under the fixed blind state. + + Reads the encrypted seed and ``commitment.json`` written by + :func:`blind_init` (verifying both hashes), conceals the part through its + matching backend (:func:`blind_sacc`), writes the blinded part beside the + input with the part's true vector escrowed into a per-part Fernet bundle, + and deletes the plaintext part — only the blinded part persists on disk. + Each part's escrow is self-contained: restoring it needs only that + part's bundle plus the seed bundle; corruption of one bundle loses one + part, not all. Pass ``keep_input=True`` to retain the plaintext part + (and own the custody implication). + + Returns + ------- + dict + Paths written: ``blinded``, ``escrow``, ``escrow_key``. + """ + config = config or BlindingConfig() + seed, commitment = _read_seed(blind_dir, config) + paths = part_paths(part_path) + for path in paths.values(): + if os.path.exists(path): + raise FileExistsError( + f"refusing to overwrite existing blind output {path} — a blind " + "is a one-shot custody event" + ) + + # The plaintext part is real, unblinded data — the blinding step is the one + # legitimate reader of the true vector, so it passes the load escape hatch. + part = sacc_io.load(part_path, allow_unblinded=True) + blinded = blind_sacc(part, seed, config=config, label=commitment["label"], log=log) + + _write_encrypted_json( + paths["escrow"], + { + "label": commitment["label"], + "seed_commitment": commitment["seed_commitment"], + "true_mean": np.asarray(part.mean, dtype=float).tolist(), + }, + ) + # The blinded file inherits the part's provenance (data vs mock); it also + # carries concealed=True (stamped by blind_sacc), so it loads without the + # escape hatch. + sacc_io.save(blinded, paths["blinded"], type=part.metadata["type"]) + if not keep_input: + os.remove(part_path) + log(f"[blind-part] deleted plaintext part {part_path}") + else: + log(f"[blind-part] plaintext part RETAINED at {part_path} (keep_input=True)") + log(f"[blind-part] wrote {paths['blinded']} (escrow beside it)") + return paths + + +def unblind_part(blinded_path, blind_dir, out_path, config=None, log=print): + """Unblind one blinded part (or the assembled file), verifying first. + + Decrypts the seed bundle and verifies the seed commitment, the config digest + and the draw scheme against ``commitment.json``, and the same three + against the blinded file's own stamps (fail closed on any — + verification precedes subtraction), recomputes the part's shift from the + seed and subtracts it (:func:`unblind_sacc`). **The seed-subtracted + vector is the authority** — the seed plus its commitment is the custody + root of trust, and the returned vector is what the seed math produced. + + When the part's escrow bundle exists beside the blinded file it serves + two subordinate roles, and only after its stored ``seed_commitment`` is + verified to match the same commitment (so an escrow from a different + blind can never be trusted): (1) a *tighter equality check* — the seed + subtraction and the escrowed truth must agree to ``1e-6`` relative or the + unblind fails closed; (2) *ulp-residue removal* — float add-then-subtract + leaves ~ulp residue against the pre-blind truth, and once the escrow is + bound to this commitment its exact value clears that residue so the + restore is bit-for-bit. The escrow is never the source of correctness: + if it disagrees materially the guard raises, and a mismatched or + unbound escrow never determines the output. On the assembled file — no + escrow — the subtraction stands alone, with integration rows selected by the + ``grid`` tag. + """ + config = config or BlindingConfig() + seed, commitment = _read_seed(blind_dir, config) + blinded = sacc_io.load(blinded_path) + part = unblind_sacc(blinded, seed, config=config, log=log) + + stem, ext = os.path.splitext(blinded_path) + unblinded_stem = os.path.join( + os.path.dirname(stem), os.path.basename(stem).replace("_blinded", "") + ) + escrow = part_paths(unblinded_stem + ext) + if os.path.exists(escrow["escrow"]): + bundle = _read_encrypted_json(escrow["escrow"], escrow["escrow_key"]) + if bundle.get("seed_commitment") != commitment["seed_commitment"]: + raise ValueError( + "escrow bundle beside the blinded file was written under a " + "different seed than the commitment — refusing to trust it " + "(the seed subtraction is authoritative; this escrow is not " + "bound to this blind)" + ) + true_mean = np.asarray(bundle["true_mean"], dtype=float) + recovered = np.asarray(part.mean, dtype=float) + residual = np.nanmax( + np.abs(recovered - true_mean) / (np.abs(true_mean) + 1e-30) + ) + if residual > 1e-6: + raise ValueError( + f"unblinded vector disagrees with the escrowed true vector " + f"(max rel {residual:.2e}) — wrong escrow for this part?" + ) + # Seed-bound escrow: clear the add-then-subtract ulp residue so the + # restore is bit-for-bit. Correctness already came from the seed + # subtraction above; the guard proved the escrow agrees with it. + _set_values(part, np.arange(len(true_mean)), true_mean) + log(f"[unblind] escrow verified (subtraction residual {residual:.2e})") + # unblind_sacc stripped the concealed/blind stamps, so this is the true + # revealed vector; it inherits the blinded file's provenance (data vs mock). + sacc_io.save(part, out_path, type=blinded.metadata["type"]) + log(f"[unblind] wrote {out_path}") + return out_path + + +def _write_encrypted_json(encrpt_path, payload): + """Encrypt ``payload`` (JSON) to ``encrpt_path`` + sibling ``.key``; no + plaintext survives. + + We drive ``smokescreen.encryption.encrypt_file`` (Fernet) for the crypto + but write the ciphertext and key at the exact names we control. Its + ``save_file`` mode names outputs from ``basename.split('.')[0]`` — it + truncates at the first dot — so for a dotted stem (the canonical + catalogue-version case, e.g. ``v1.4.6.3_xi_integration_escrow``) it would land at + ``v1.encrpt``/``v1.key``, diverging from what :func:`part_paths` declares + and colliding across parts that share a first-dot prefix. Instead we take + the returned ``(ciphertext, key)`` and write them ourselves. + """ + from smokescreen.encryption import encrypt_file + + key_path = encrpt_path.replace(".encrpt", ".key") + plaintext = encrpt_path.replace(".encrpt", ".json") + with open(plaintext, "w", encoding="utf-8") as f: + json.dump(payload, f) + ciphertext, key = encrypt_file(plaintext, save_file=False, keep_original=False) + with open(encrpt_path, "wb") as f: + f.write(ciphertext) + with open(key_path, "wb") as f: + f.write(key) + + +def _read_encrypted_json(encrpt_path, key_path): + """Decrypt and parse a Fernet-encrypted JSON bundle.""" + from smokescreen.encryption import decrypt_file + + return json.loads(decrypt_file(encrpt_path, key_path).decode("utf-8")) diff --git a/src/sp_validation/blinding_theory.py b/src/sp_validation/blinding_theory.py new file mode 100644 index 00000000..b44f7e73 --- /dev/null +++ b/src/sp_validation/blinding_theory.py @@ -0,0 +1,424 @@ +"""Blinding theory: fiducial configuration and the two ξ± theory paths. + +:Name: blinding_theory.py + +:Description: The blinding backend's theory surface — the fiducial + configuration (:class:`TheoryConfig`) and two independent routes to the + tomographic shear two-point prediction. + + **Division of responsibility** (per the UNIONS layering: ``cs_util`` is + the cosmology *library* — generic machinery; ``sp_validation`` holds + *configuration* and *survey-specific implementations*). This module lives + in the blinding namespace because :class:`TheoryConfig` is configuration + and the master-layout theory backends in :mod:`sp_validation.blinding` + (reporting ξ±, integration ξ±, pseudo-Cℓ) are survey-specific. The **generic** + theory machinery below — the CCL-native ξ± path (:func:`xi_ccl`, + :func:`cl_ee`), the independent CAMB P(k)→``Pk2D`` path (:func:`xi_camb`), + and the σ8/A_s rescale (:func:`camb_As_for_sigma8`) — lives here **for + now** but is destined for ``cs_util.cosmo`` (tracked in cs_util#80): it is + cosmology-library code, not blinding-specific. There is deliberately **no** + ``sp_validation/cosmology.py`` — ``develop`` removed the local cosmology + module (#223) and moved cosmology to ``cs_util.cosmo``; this module does + not resurrect it. + + Two independent routes to the shear two-point prediction: + + - **CCL-native path** (:func:`xi_ccl`, :func:`cl_ee`): CCL builds the + nonlinear P(k) through its Boltzmann-CAMB HMCode2020 route + (``matter_power_spectrum='camb'`` + ``extra_parameters``) and projects + to Cℓ/ξ± via its own Limber (``angular_cl``) + FFTLog + (``correlation``). This is the recipe the blinding theory backends use. + - **Independent-CAMB path** (:func:`xi_camb`): a direct ``pycamb`` run + produces the HMCode2020 ``P(k, z)`` (σ8-matched via the closed-form + A_s rescale of :func:`camb_As_for_sigma8`), wrapped in a ``ccl.Pk2D`` + and projected through the same CCL Limber + FFTLog machinery. + + Because both paths route their nonlinear P(k) through CAMB's HMCode2020 + and both project through CCL, a common Limber+FFTLog bug cancels between + them: the CAMB↔CCL cross-check test built on these two paths validates + the **P(k) recipe** and the **σ8/A_s amplitude convention**, not the + projection machinery. + + This module imports only ``numpy`` at module level; CCL and CAMB are + imported inside the functions that need them, so importing + :class:`TheoryConfig` never drags in a theory backend. +""" + +import dataclasses + +import numpy as np + +# Fixed constants of the fiducial — load-bearing for the CAMB↔CCL amplitude +# match, so they are emitted explicitly to both stacks rather than left to +# either stack's default. Not user-facing TheoryConfig fields. +NEFF = 3.046 +T_CMB = 2.7255 + + +# --------------------------------------------------------------------------- # +# Configuration surface — the ONE place fiducial cosmology + model choices live +# --------------------------------------------------------------------------- # +@dataclasses.dataclass(frozen=True) +class TheoryConfig: + """Fiducial cosmology and model configuration for the theory paths. + + Every field is a deliberate, configurable choice. The defaults mirror 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. Adopting a different named group fiducial is a + change to these *values*, not to any code. + + Cosmology is parametrised by the blind axes ``S8`` and ``Omega_m`` and + converted to CCL's native ``sigma8``/``Omega_c`` by :meth:`sigma8` / + :meth:`omega_c`. + + **Nonlinear-model tokens (load-bearing).** The HMCode2020+feedback recipe + is named by a token in each stack's API: CCL takes + ``extra_parameters['camb']['halofit_version']``, CAMB takes + ``NonLinearModel.set_params(halofit_version=…)``. :class:`TheoryConfig` + carries **two** tokens (``ccl_halofit_version``, ``camb_halofit_version``) + denoting one recipe, and each stack is fed its own — a single shared + string invites a silent stack disagreement the moment the two APIs name + the recipe differently (``mead2020`` vs ``mead2020_feedback`` differ by + several % at k ≳ 1/Mpc). Today both stacks accept the same string for the + feedback recipe, so the two defaults coincide; the cross-check test pins + the CCL token against the inference config independently. + """ + + # 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 + Omega_b: float = 0.0469 # ombh2=0.023 at h=0.7 -> 0.023/0.7^2 + h: float = 0.70 + n_s: float = 0.96 + m_nu: float = 0.06 # Σm_ν in eV, distributed under `mass_split` + w0: float = -1.0 + wa: float = 0.0 + + # Neutrino mass split: normal hierarchy (CosmoSIS `neutrino_hierarchy=normal`). + mass_split: str = "normal" + + # Boltzmann/transfer-function backend for the CCL path (#280). The default + # `boltzmann_camb` makes CCL call CAMB for the linear P(k), so blinding + # theory and the CosmoSIS+CAMB inference stack share one power-spectrum + # path. The halofit route follows the choice (see `ccl_cosmology`): only + # under `boltzmann_camb` can the nonlinear P(k) run through CAMB's HMCode + # (`matter_power_spectrum="camb"` + the tokens below); any other backend + # falls back to CCL's own halofit — consistent, but a different recipe, + # so a non-default backend is a deliberate cross-check tool, not a + # production setting. + transfer_function: str = "boltzmann_camb" + + # One nonlinear recipe (CAMB HMCode2020 + baryonic feedback), two + # stack-specific tokens — see the class docstring. + ccl_halofit_version: str = "mead2020_feedback" + camb_halofit_version: str = "mead2020_feedback" + hmcode_logT_AGN: float = 7.5 # values_ia.ini logT_AGN central + + # Intrinsic alignments: NLA. The fiducial defaults IA OFF (ia_bias=0) — + # the blinding shift is a difference of two theory vectors at the same IA, + # so IA nearly cancels there, and IA-off keeps 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_bias: float = 0.0 + ia_z_piv: float = 0.62 + ia_alphaz: float = 0.0 + + 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)``. At the fiducial (S8=0.80, Ωm=0.30), + σ8 = 0.80. + """ + 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 − Ω_ν``. + + The neutrino density ``Ω_ν h² = Σm_ν / 93.14 eV`` is subtracted so + 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): + """The fiducial point as a plain CCL-native parameter mapping. + + Exactly the keys ``Omega_c, Omega_b, h, n_s, sigma8, m_nu, + mass_split, w0, wa, Neff, T_CMB`` and no others — no CCL default + rides along. ``Neff``/``T_CMB`` are the fixed module constants. This + mapping is what the Smokescreen fork receives as ``fiducial_params`` + and what every ``theory_fn`` receives back (possibly with + ``sigma8``/``Omega_c`` overlaid by the hidden draw). + """ + 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, + "mass_split": self.mass_split, + "w0": self.w0, + "wa": self.wa, + "Neff": NEFF, + "T_CMB": T_CMB, + } + + @classmethod + def from_overrides(cls, overrides): + """Build from a mapping of field overrides (fail loud on unknown keys). + + Numeric overrides are coerced to ``float`` per the field's declared + type, so a YAML/CLI ``w0: -1`` (int) yields the same value — and the + same :meth:`config_digest` — as the float default ``-1.0``. + """ + by_name = {f.name: f for f in dataclasses.fields(cls)} + unknown = set(overrides) - set(by_name) + if unknown: + raise ValueError( + f"unknown TheoryConfig fields {sorted(unknown)}; " + f"valid fields are {sorted(by_name)}" + ) + coerced = { + name: (float(v) if by_name[name].type in (float, "float") else v) + for name, v in overrides.items() + } + return cls(**coerced) + + +# --------------------------------------------------------------------------- # +# CCL-native path: cosmology construction, Cℓ_EE, ξ± +# --------------------------------------------------------------------------- # +# The two cosmologies of a blind (fiducial + hidden) are evaluated by three +# theory backends over multiple blocks; caching the ccl.Cosmology per parameter +# point avoids re-running the CAMB P(k) computation for every block. +_COSMO_CACHE = {} + + +def ccl_cosmology(params, config): + """A ``pyccl.Cosmology`` at ``params`` with ``config``'s nonlinear recipe. + + ``params`` is a plain CCL-native mapping (:meth:`TheoryConfig.ccl_params`, + possibly with keys overlaid by the hidden draw); ``config`` supplies only + the non-sampled recipe tokens (``transfer_function``, + ``ccl_halofit_version``, ``hmcode_logT_AGN``). The Boltzmann backend is + ``config.transfer_function`` (#280); the halofit route stays consistent + with that choice: under ``boltzmann_camb`` the nonlinear P(k) runs + through CAMB's HMCode2020 (``matter_power_spectrum="camb"`` + the CAMB + tokens), while any other backend has no CAMB run to hand tokens to, so it + takes CCL's own halofit (``matter_power_spectrum="halofit"``). Either + way, the same recipe sits on both sides of any theory difference. + Cosmology objects are cached per parameter point (CCL memoises its P(k) + on the object, so the cache saves repeated Boltzmann runs across blocks). + """ + import pyccl as ccl + + key = ( + tuple(sorted(params.items())), + config.transfer_function, + config.ccl_halofit_version, + config.hmcode_logT_AGN, + ) + if key not in _COSMO_CACHE: + nonlinear = ( + { + "matter_power_spectrum": "camb", + "extra_parameters": { + "camb": { + "halofit_version": config.ccl_halofit_version, + "HMCode_logT_AGN": config.hmcode_logT_AGN, + } + }, + } + if config.transfer_function == "boltzmann_camb" + else {"matter_power_spectrum": "halofit"} + ) + _COSMO_CACHE[key] = ccl.Cosmology( + **params, + transfer_function=config.transfer_function, + **nonlinear, + ) + return _COSMO_CACHE[key] + + +def xi_ell_grid(): + """The ℓ grid the ξ± Hankel projection integrates over. + + Integers 2…49, then 200 log-spaced multipoles up to 6·10⁴ — dense enough + at low ℓ (where ξ± at large θ lives) and wide enough for the small-θ + tail. ``ccl.correlation`` interpolates C(ℓ) internally, so this fixes the + resolution of every ξ± this module produces. + """ + return np.unique( + np.concatenate([np.arange(2, 50), np.geomspace(50, 6e4, 200)]).astype(float) + ) + + +def _tracer(cosmo, z, nz, config): + """A ``WeakLensingTracer`` for one bin's n(z), NLA from ``config``. + + With the fiducial ``ia_bias = 0`` the tracer is built bare — no IA term. + A nonzero ``ia_bias`` enters as the NLA amplitude + ``A(z) = ia_bias · ((1+z)/(1+z_piv))^alphaz``. + """ + import pyccl as ccl + + z = np.asarray(z) + if config.ia_bias == 0.0: + return ccl.WeakLensingTracer(cosmo, dndz=(z, np.asarray(nz))) + a_ia = config.ia_bias * ((1 + z) / (1 + config.ia_z_piv)) ** config.ia_alphaz + return ccl.WeakLensingTracer( + cosmo, dndz=(z, np.asarray(nz)), ia_bias=(z, a_ia), use_A_ia=True + ) + + +def cl_ee(params, config, nz_i, nz_j, ell): + """Cross Cℓ_EE at ``ell`` for the bin pair with n(z) ``nz_i``, ``nz_j``. + + Two-tracer: one :class:`~pyccl.WeakLensingTracer` per bin from that bin's + own ``(z, nz)``, then ``angular_cl(cosmo, tracer_i, tracer_j, ell)`` — + the cross-spectrum for i ≠ j, the auto-spectrum when the two n(z) are the + same bin. The shear ``angular_cl`` is the E-mode spectrum; B and EB are + zero in theory, which is why only Cℓ_EE ever receives a blinding shift. + """ + import pyccl as ccl + + cosmo = ccl_cosmology(params, config) + tracer_i = _tracer(cosmo, *nz_i, config) + tracer_j = _tracer(cosmo, *nz_j, config) + return ccl.angular_cl(cosmo, tracer_i, tracer_j, np.asarray(ell, dtype=float)) + + +def xi_ccl(params, config, nz_i, nz_j, theta_arcmin, ell=None): + """CCL-native ξ± at ``theta_arcmin`` for one bin pair (Path A). + + Cross Cℓ_EE on :func:`xi_ell_grid` (or ``ell``), then ``ccl.correlation`` + (FFTLog Hankel transform) at θ in degrees, ``type="GG+"`` / ``"GG-"``. + + Returns + ------- + (np.ndarray, np.ndarray) + ``(xip, xim)`` aligned to ``theta_arcmin``. + """ + import pyccl as ccl + + ell = xi_ell_grid() if ell is None else np.asarray(ell, dtype=float) + cosmo = ccl_cosmology(params, config) + cl = cl_ee(params, config, nz_i, nz_j, ell) + theta_deg = np.asarray(theta_arcmin) / 60.0 + xip = ccl.correlation(cosmo, ell=ell, C_ell=cl, theta=theta_deg, type="GG+") + xim = ccl.correlation(cosmo, ell=ell, C_ell=cl, theta=theta_deg, type="GG-") + return xip, xim + + +# --------------------------------------------------------------------------- # +# Independent-CAMB path: A_s reconciliation + P(k) → Pk2D → CCL projection +# --------------------------------------------------------------------------- # +def make_camb_params(config, As, *, nonlinear, zmax=3.0, n_z=48, kmax=20.0): + """A ``CAMBparams`` at ``config``'s background with amplitude ``As``. + + Every :class:`TheoryConfig` field CCL sees is fed to CAMB from the same + source — ``w0``/``wa`` via ``set_dark_energy``, ``Neff``/``T_CMB`` as the + module constants, ``m_nu``/``mass_split`` through ``set_cosmology`` — so + the independent path differs from the CCL path only in who computes P(k), + never in an unmatched background parameter. + """ + import camb + + p = camb.CAMBparams() + p.set_cosmology( + H0=config.h * 100, + ombh2=config.Omega_b * config.h**2, + omch2=config.omega_c() * config.h**2, + mnu=config.m_nu, + num_massive_neutrinos=1, + neutrino_hierarchy=config.mass_split, + nnu=NEFF, + TCMB=T_CMB, + ) + p.set_dark_energy(w=config.w0, wa=config.wa, dark_energy_model="ppf") + p.InitPower.set_params(As=As, ns=config.n_s) + p.set_matter_power(redshifts=list(np.linspace(0.0, zmax, n_z)), kmax=kmax) + if nonlinear: + p.NonLinear = camb.model.NonLinear_both + p.NonLinearModel.set_params( + halofit_version=config.camb_halofit_version, + HMCode_logT_AGN=config.hmcode_logT_AGN, + ) + else: + p.NonLinear = camb.model.NonLinear_none + return p + + +def camb_linear_sigma8(config, As, **kwargs): + """CAMB's linear σ8(z=0) at amplitude ``As``.""" + import camb + + results = camb.get_results(make_camb_params(config, As, nonlinear=False, **kwargs)) + return float(results.get_sigma8_0()) + + +def camb_As_for_sigma8(config, sigma8_target, As_seed=2.1e-9, **kwargs): + """The CAMB ``A_s`` whose linear σ8 equals ``sigma8_target``. + + Closed-form: linear σ8² ∝ A_s exactly, so one CAMB linear-σ8 evaluation + at ``As_seed`` and one rescale ``As_seed · (σ8_target/σ8_seed)²`` land on + the target — no iteration. This settles the convention subtlety that our + fiducial fixes σ8 for CCL but A_s for CAMB: a nominal ``A_s = 2.1e-9`` + leaves CAMB's σ8 ≈3% off target, enough to blow a ξ± comparison to + ~9–10%. + """ + sigma8_seed = camb_linear_sigma8(config, As_seed, **kwargs) + return As_seed * (sigma8_target / sigma8_seed) ** 2 + + +def xi_camb(config, nz, theta_arcmin, *, n_ell=300, ell_max=60000, kmax=20.0, n_k=400): + """Independent-CAMB ξ± for one bin (Path B): CAMB P(k) → Pk2D → CCL. + + A direct pycamb run produces the HMCode2020 nonlinear ``P(k, z)`` at a + σ8-matched ``A_s`` (:func:`camb_As_for_sigma8`), extracted through + ``get_matter_power_interpolator(hubble_units=False, k_hunit=False)`` so + it comes out in CCL's native units (k in 1/Mpc, P in Mpc³) — **no** + ``·h`` / ``/h³`` conversion is applied (applying one would double-count + an h³ amplitude error). Both ``Pk2D`` axes are arranged ascending + (log-k ascending; scale factor ascending, i.e. CAMB's z-ascending grid + reversed). Projection is CCL's own Limber + FFTLog with a bare tracer + (IA off — this path exists for the cross-check). + + Returns + ------- + (np.ndarray, np.ndarray, float) + ``(xip, xim, As)`` — the σ8-matched amplitude is returned for + assertion by the cross-check test. + """ + import camb + import pyccl as ccl + + sigma8 = config.sigma8() + As = camb_As_for_sigma8(config, sigma8, kmax=kmax) + results = camb.get_results(make_camb_params(config, As, nonlinear=True, kmax=kmax)) + interp = results.get_matter_power_interpolator( + nonlinear=True, hubble_units=False, k_hunit=False + ) + k = np.geomspace(1e-4, kmax * config.h, n_k) # 1/Mpc + z = np.linspace(0.0, 3.0, 48) + pk = interp.P(z, k) # (n_z, n_k), Mpc^3 + a = 1.0 / (1.0 + z) + 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(config.ccl_params(), config) + z_nz, nz_vals = nz + lens = ccl.WeakLensingTracer(cosmo, dndz=(np.asarray(z_nz), np.asarray(nz_vals))) + 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 = np.asarray(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 diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py index 0c9273d8..05ded922 100644 --- a/src/sp_validation/cosmo_val/core.py +++ b/src/sp_validation/cosmo_val/core.py @@ -98,6 +98,11 @@ class CosmologyValidation( noise debiasing, making those realizations reproducible run-to-run. cosmo_params : dict, optional Cosmological parameters to pass to get_cosmo(). If None, uses Planck 2018. + run_type : {'data', 'mock'}, default 'data' + The campaign's run type, stamped as the SACC ``type`` of every part this + object writes. Custody state, not decoration: a mock campaign must be + built with ``run_type='mock'`` for its parts to assemble at all (see + ``blinding.assert_consistent_blind``). Attributes ---------- @@ -224,6 +229,7 @@ def __init__( path_onecovariance=None, cosmo_params=None, blind=None, + run_type="data", ): self.rho_tau_method = rho_tau_method self.cov_estimate_method = cov_estimate_method @@ -253,6 +259,7 @@ def __init__( self.nside_mask = nside_mask self.path_onecovariance = path_onecovariance self.blind = blind + self.run_type = run_type assert self.cell_method in ["map", "catalog"], ( "cell_method must be 'map' or 'catalog'" diff --git a/src/sp_validation/cosmo_val/cosebis.py b/src/sp_validation/cosmo_val/cosebis.py index 406f0a11..4fbec2ba 100644 --- a/src/sp_validation/cosmo_val/cosebis.py +++ b/src/sp_validation/cosmo_val/cosebis.py @@ -165,7 +165,13 @@ def _fiducial_cosebis_result(results, fiducial_scale_cut): return results[key], tuple(key) def cosebis_to_sacc_part( - self, version, out_path, results, fiducial_scale_cut=None, en_override=None + self, + version, + out_path, + results, + fiducial_scale_cut=None, + en_override=None, + commitment_path=None, ): """Write the COSEBIs SACC part at the fiducial scale cut. @@ -178,7 +184,9 @@ def cosebis_to_sacc_part( ``en_override`` is the consume-the-part plumbing: the E-mode ``En`` written to the part in place of ``result["En"]`` — re-derived from the integration ξ± SACC part at the fiducial scale cut (Bn and the covariance stay from the - raw estimator ``result``). With ``None`` the behaviour is unchanged. + raw estimator ``result``). ``commitment_path`` stamps the part concealed + under that version's blind (see :func:`sacc_io.save`). With both ``None`` + (mock runs) the behaviour is unchanged. """ result, scale_cut = self._fiducial_cosebis_result(results, fiducial_scale_cut) if en_override is not None: @@ -189,7 +197,7 @@ def cosebis_to_sacc_part( result, scale_cut, ) - sacc_io.save(s, out_path, type="data") + sacc_io.save(s, out_path, type=self.run_type, commitment=commitment_path) def plot_cosebis( self, diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index 2a55934a..cfaa349e 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -708,7 +708,7 @@ def pseudo_cl_to_sacc_part(self, version, out_path, ell_eff, cl_all, wsp): cl_all, wsp, ) - sacc_io.save(s, out_path, type="data") + sacc_io.save(s, out_path, type=self.run_type) def plot_pseudo_cl(self): """ diff --git a/src/sp_validation/cosmo_val/psf_systematics.py b/src/sp_validation/cosmo_val/psf_systematics.py index 6f72ac1f..c7984c26 100644 --- a/src/sp_validation/cosmo_val/psf_systematics.py +++ b/src/sp_validation/cosmo_val/psf_systematics.py @@ -25,7 +25,8 @@ class PSFSystematicsMixin: - def calculate_rho_tau_stats(self): + def calculate_rho_tau_stats(self, commitment_path=None): + """Measure ρ/τ statistics per version and write each version's SACC part.""" out_dir = f"{self.cc['paths']['output']}/rho_tau_stats" if not os.path.exists(out_dir): os.mkdir(out_dir) @@ -44,7 +45,12 @@ def calculate_rho_tau_stats(self): npatch=self.npatch, ) self.rho_tau_to_sacc_part( - ver, out_dir, base, rho_stat_handler, tau_stat_handler + ver, + out_dir, + base, + rho_stat_handler, + tau_stat_handler, + commitment_path=commitment_path, ) self.print_done("Rho stats finished") @@ -52,10 +58,21 @@ def calculate_rho_tau_stats(self): self._tau_stat_handler = tau_stat_handler def rho_tau_to_sacc_part( - self, version, out_dir, base, rho_stat_handler, tau_stat_handler + self, + version, + out_dir, + base, + rho_stat_handler, + tau_stat_handler, + commitment_path=None, ): """Write the ρ/τ SACC part for one version. + ρ/τ is a PSF diagnostic carrying no cosmological vector, so it is never + blinded; ``commitment_path`` stamps the part concealed under the + version's blind, values untouched (see :func:`sacc_io.save`), so a data + run's assembly admits it. A mock run passes ``None``. + ρ_0…ρ_5 autos and τ_0/τ_2/τ_5 leakage from the handler tables. The ``CovTauTh`` theory covariance ``cov_tau_{base}_th.npy`` — a ``(3·nbin, 3·nbin)`` plus-folded k-major block over ``{τ0, τ2, τ5}`` — is @@ -80,7 +97,7 @@ def rho_tau_to_sacc_part( tau_cov_th=tau_cov_th, ) out_path = os.path.join(out_dir, f"rho_tau_{base}.sacc") - sacc_io.save(s, out_path, type="data") + sacc_io.save(s, out_path, type=self.run_type, commitment=commitment_path) @property def rho_stat_handler(self): diff --git a/src/sp_validation/cosmo_val/pure_eb.py b/src/sp_validation/cosmo_val/pure_eb.py index 29e1f888..e7220c47 100644 --- a/src/sp_validation/cosmo_val/pure_eb.py +++ b/src/sp_validation/cosmo_val/pure_eb.py @@ -134,7 +134,9 @@ def calculate_pure_eb( return results - def pure_eb_to_sacc_part(self, version, out_path, results, eb_override=None): + def pure_eb_to_sacc_part( + self, version, out_path, results, eb_override=None, commitment_path=None + ): """Write the pure-E/B SACC part (six ``PURE_KEYS`` blocks + covariance). ``results`` is the dict ``calculate_pure_eb`` returned: the six pure-mode @@ -145,8 +147,10 @@ def pure_eb_to_sacc_part(self, version, out_path, results, eb_override=None): ``eb_override`` is the consume-the-part plumbing: the six pure-mode arrays (a mapping keyed by ``sacc_io.PURE_KEYS``) written in place of ``results``' — re-derived from the reporting + integration ξ± SACC parts (the covariance - stays blind-invariant from the raw estimator ``results``). With ``None`` the - behaviour is unchanged. + stays blind-invariant from the raw estimator ``results``). + ``commitment_path`` stamps the part concealed under that version's blind + (see :func:`sacc_io.save`). With both ``None`` (mock runs) the behaviour + is unchanged. """ theta = results["gg"].meanr source = eb_override if eb_override is not None else results @@ -158,7 +162,7 @@ def pure_eb_to_sacc_part(self, version, out_path, results, eb_override=None): eb, covariance=results["cov"], ) - sacc_io.save(s, out_path, type="data") + sacc_io.save(s, out_path, type=self.run_type, commitment=commitment_path) def plot_pure_eb( self, diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index 99a0efe2..287f1e69 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -942,7 +942,7 @@ def update_statistic(s, sub): s.data[idx[0]].value = point.value -def save(s, path, *, type): +def save(s, path, *, type, commitment=None): """Write ``s`` to ``path`` (FITS), overwriting any existing file. Parameters @@ -955,6 +955,12 @@ def save(s, path, *, type): the pipeline computing the data vector — knows whether its input catalogue is a mock; there is deliberately no default. ``load`` refuses ``type='data'`` files that are not blinded. + commitment : str, optional + Path to the version's ``commitment.json``. When given, the file is + stamped concealed under that blind + (:func:`sp_validation.blinding.stamp_concealed_passthrough`, values + untouched) before writing — the seam every born-blinded or + blind-irrelevant part uses to clear the fail-closed load gate. """ if type not in ("data", "mock"): raise ValueError(f"type must be 'data' or 'mock'; got {type!r}") @@ -964,6 +970,10 @@ def save(s, path, *, type): f"refusing to re-stamp as {type!r}" ) s.metadata["type"] = type + if commitment is not None: + from . import blinding + + blinding.stamp_concealed_passthrough(s, commitment) s.save_fits(path, overwrite=True) @@ -1539,3 +1549,61 @@ def covariance_blocks(cov_list, selectors, *, gaussian=True): return [ (selectors, cov_from_one_covariance(np.asarray(cov_list), gaussian=gaussian)) ] + + +# --------------------------------------------------------------------------- # +# Terminal assembly (PR-6 blinding) — gather() and its blind-custody call site. +# --------------------------------------------------------------------------- # +def gather(parts, metadata=None, assemble=None): + """Assemble standalone part SACCs into the one-file ``{version}.sacc``. + + Each part is an intermediate product as it came off its producing rule + (reporting ξ±, integration ξ±, pseudo-Cℓ, ρ/τ, …). Gather is **the** terminal + seam: every path that combines parts into the one-file product goes + through here, because this is where the one thing an assembler cannot know + about is enforced — **blind custody.** + + **Blind custody.** + :func:`sp_validation.blinding.assert_consistent_blind` runs before the + assembly — it fails closed unless every blindable part carries the + identical ``blind_commitment``/``blind_config_digest``/``blind_draw_scheme`` + (or, when nothing is blinded, every blindable part is declared + ``type='mock'``). Its returned shared stamp is written onto the assembled + file so the one-file product carries the blind it was built from; the + blinded parts already carry those keys, so the assembly preserves them and + this stamp is a consistent (idempotent) re-affirmation. + + **The assembly itself is the caller's.** Two exist and both are real: + :func:`merge` (the default) does the first-wins tracer union, in-order + point concatenation with all tags, and a covariance built from whatever the + parts carry — the right thing when parts are already covariance-bearing. + :func:`sp_validation.cosmo_val.sacc_writers.assemble_analysis_sacc` rebuilds + from a given n(z) + metadata and *requires* one covariance block per part — + the right thing for the production terminal, where the ξ± and pseudo-Cℓ + parts are born cov-less and have their blocks injected first. Passing the + assembler in, rather than duplicating the custody wrapper around each one, + is what keeps the guard un-bypassable. + + Parameters + ---------- + parts : sequence of sacc.Sacc + The part SACCs, in the assembly (covariance) order. + metadata : dict, optional + Extra key/value pairs to store on the assembled file's metadata. + assemble : callable, optional + ``assemble(parts) -> sacc.Sacc``. Defaults to :func:`merge`. Bind any + further arguments (n(z), metadata) into the callable. + + Returns + ------- + sacc.Sacc + The assembled file. + """ + from . import blinding + + parts = list(parts) + stamp = blinding.assert_consistent_blind(parts) + s = (assemble or merge)(parts) + for key, value in {**(metadata or {}), **(stamp or {})}.items(): + s.metadata[key] = value + return s diff --git a/src/sp_validation/sacc_like_unions.py b/src/sp_validation/sacc_like_unions.py new file mode 100644 index 00000000..a25837cb --- /dev/null +++ b/src/sp_validation/sacc_like_unions.py @@ -0,0 +1,225 @@ +"""SACC likelihood shim for CosmoSIS — the sp_validation-owned native path. + +This module is a thin subclass of CosmoSIS's ``SaccClLikelihood`` (from the +CosmoSIS Standard Library, ``likelihood/sacc/sacc_like.py``) that fixes two +upstream defects so the native SACC likelihood matches the PR-3 converter → +``2pt_like`` path bit for bit on our real-space ξ± analysis file. It exists as a +CosmoSIS *module file* (``setup``/``execute``/``cleanup`` at module scope), loaded +via ``file = .../sacc_like_unions.py`` in an ini, and is NOT imported by +``sp_validation/__init__`` (cosmosis is an optional dependency). + +Why the shim exists +------------------- +Two things in upstream ``sacc_like`` break a real-space ξ likelihood; both are +documented and empirically verified (probe: Δχ²=3184 raw vs Δχ²=0 shimmed on the +single-bin analysis file), and both are fixed here by overriding ``build_data``: + +1. **No arcmin→radian conversion (the killer).** + ``sacc_likelihoods/twopoint.py`` (L74-90 at CSL commit 4fd2f1c) builds a + ``SpectrumInterp`` over ``block[section, "theta"]`` — theory θ in *radians* + (CosmoSIS convention) — and evaluates it at each data point's raw ``theta`` + tag, which our SACC files (and firecrown) store in *arcmin*. ``2pt_like.py`` + (L179-183) converts its real-space data to radians for exactly this reason; + ``sacc_like`` never does, so the spline is evaluated ~3437× outside its grid, + ``SpectrumInterp`` returns 0 there, and χ² silently collapses to dᵀC⁻¹d. The + only prior use was the ℓ-space (unit-free) Cℓ path, which is why it was never + caught. We evaluate the theory against a radian-θ *copy* of the loaded SACC, + keeping ``self.sacc_data`` in arcmin so upstream's save_theory / + save_realization paths (which copy it and overwrite only values) never write + radian tags into a file downstream consumers read as arcmin. + +2. **Theory↔data ordering is assumed, never enforced.** + The data vector is ``sacc.get_mean()`` (insertion order); the theory vector is + built by looping ``get_data_types() × get_tracer_combinations() × points``. + A comment in ``twopoint.py`` (L35) claims ``to_canonical_order`` was called on + load — it is not. The two orders agree only when the file is grouped + type-major (all of one data type, then the next). Our single-pair + ``[ξ+; ξ−]`` files satisfy this; a tomographic *pair-major* file (ξ+/ξ− per + pair) would silently misalign theory against data. We reconstruct the + theory-loop index order and require it equals ``arange`` — a hard ValueError + otherwise (see ``test_ordering_guard_raises_pair_major_tomographic``). + +When upstream fixes the units, this shim dies: the tripwire test +``test_upstream_unit_gap_tripwire`` in ``tests/test_sacc_like.py`` fails the day +raw ``sacc_like`` stops producing a wildly different χ², signalling the shim can +be retired. + +Requires a CSL checkout: ``setup`` reads ``csl_dir`` from the module options and +imports the upstream ``sacc_like`` from ``/likelihood/sacc``. The tests +locate it via the ``CSL_DIR`` environment variable (see the test module docstring +for the checkout recipe). +""" + +import os +import sys + +import numpy as np + +# cosmosis is an OPTIONAL dependency: this module is a CosmoSIS module file, but +# importing it must succeed without cosmosis installed (the CI image has the +# science stack but no cosmosis, and test_imports.py bare-imports every module). +# So every cosmosis touch is deferred into setup() — nothing at top level imports +# it. numpy is fine at top level (always present). + +# arcmin → radian: the conversion 2pt_like applies to real-space data and that +# sacc_like omits. Applied only to `theta` tags of `real`-category data types. +ARCMIN_TO_RAD = np.pi / (180.0 * 60.0) + + +def _import_upstream_sacc_like(csl_dir): + """Import the upstream ``sacc_like`` module from a CSL checkout. + + ``/likelihood/sacc`` is prepended to ``sys.path`` so both + ``sacc_like`` and its sibling ``sacc_likelihoods`` package (imported by + ``sacc_like`` for the theory-extraction functions) resolve. Idempotent: the + path is only inserted once. + """ + sacc_dir = os.path.join(csl_dir, "likelihood", "sacc") + if not os.path.isdir(sacc_dir): + raise ValueError( + f"csl_dir={csl_dir!r} has no likelihood/sacc directory; point csl_dir " + "at a CosmoSIS Standard Library checkout (see module docstring)" + ) + if sacc_dir not in sys.path: + sys.path.insert(0, sacc_dir) + import sacc_like # noqa: E402 — resolved from the sys.path insertion above + + return sacc_like + + +def _make_subclass(sacc_like): + """Build ``SaccLikeUnions`` as a subclass of the upstream ``SaccClLikelihood``. + + A factory (not an import-time ``class ... :`` statement) so importing this + module never requires the upstream class — that dependency is deferred to + ``setup``, which has ``csl_dir`` in hand. Only ``build_data`` is overridden; + scale cuts, covariance handling (Sellentin/Hartlap), theory extraction and + ``save_theory`` all ride upstream unmodified. + """ + + class SaccLikeUnions(sacc_like.SaccClLikelihood): + """CSL ``SaccClLikelihood`` with the arcmin→rad + ordering-guard fixes.""" + + def build_data(self): + # Run the upstream build FIRST: it loads the SACC, applies data_sets + # selection and the arcmin-grammar scale cuts (matching 2pt_like's + # angle_range convention and the ini ergonomics), populates + # self.sacc_data / self.sections_for_names, and returns the + # (unit-independent) data vector we pass straight through. + x, data_vector = super().build_data() + + # self.sacc_data STAYS in arcmin — save_theory / save_realization copy + # it and only overwrite point values, so its θ tags must remain the + # units the file was written in (any consumer, incl. re-ingesting the + # saved SACC as a data_file, assumes arcmin). The radian conversion the + # theory spline needs lives on a separate copy, swapped in only for the + # extraction (see extract_theory_points). + self._sacc_data_rad = self._radian_theta_copy(self.sacc_data) + self._assert_theory_order_matches_data() + + return x, data_vector + + def _radian_theta_copy(self, sacc_data): + """A copy of ``sacc_data`` with ``real``-category θ tags in radians. + + The theory spline is built on ``block[section, "theta"]`` in radians, + so the ``theta`` tag each point is evaluated against must be radians + too. Scoped to data types whose category + (``sections_for_names[dt][0]``) is ``real``: cosebis ``n`` tags and + spectrum ``ell`` tags are unit-free and left untouched, and only the + ``theta`` tag is scaled (``theta_nom`` etc. are metadata the likelihood + never evaluates). Operates on a ``.copy()`` so the original stays + arcmin for the save paths. + """ + real_types = { + dt + for dt, (category, _section) in self.sections_for_names.items() + if category == "real" + } + converted = sacc_data.copy() + for point in converted.data: + if point.data_type in real_types and "theta" in point.tags: + point.tags["theta"] = point.tags["theta"] * ARCMIN_TO_RAD + return converted + + def extract_theory_points(self, block): + """Extract theory against the radian-θ copy, then restore the original. + + Upstream ``extract_theory_points`` reads ``self.sacc_data`` (the θ tag + per point) to evaluate the theory spline; that read needs radians. + Swap in ``self._sacc_data_rad`` for the duration of the upstream call + and restore in ``finally`` so everything else — including the + save_theory / save_realization copies that run afterward in + ``do_likelihood`` — sees the untouched arcmin ``self.sacc_data``. + """ + original = self.sacc_data + self.sacc_data = self._sacc_data_rad + try: + return super().extract_theory_points(block) + finally: + self.sacc_data = original + + def _assert_theory_order_matches_data(self): + """Require the theory-loop order to equal the data-vector order. + + The data vector is ``sacc.get_mean()`` (insertion order); upstream + builds theory by looping data types, then tracer combinations, then + points, and concatenating — assuming (never enforcing) that this + reproduces insertion order. It does only for type-major files. We + reconstruct that loop's index order and require it be ``arange``; + otherwise theory and data would be silently misaligned (the same bug + class the PR-2/PR-3 reviews caught for the converter). + """ + order = [ + int(i) + for dt in self.sacc_data.get_data_types() + for tracers in self.sacc_data.get_tracer_combinations(dt) + for i in self.sacc_data.indices(dt, tracers) + ] + expected = np.arange(len(self.sacc_data.mean)) + if not np.array_equal(order, expected): + raise ValueError( + "SACC data/theory ordering mismatch: the theory loop " + "(get_data_types × get_tracer_combinations × points) does not " + "reproduce the get_mean() insertion order, so sacc_like would " + "compare theory against data point-by-point in the WRONG order " + "and return a silently wrong χ². This happens when the file is " + "grouped pair-major (ξ+/ξ− interleaved per tracer pair) rather " + "than type-major (all ξ+, then all ξ−). Upstream assumes " + "to_canonical_order() was applied on load but never calls it; " + "write the SACC type-major, or call to_canonical_order() before " + "saving." + ) + + return SaccLikeUnions + + +def setup(options): + """CosmoSIS ``setup`` — build and instantiate the shimmed likelihood. + + Mirrors ``GaussianLikelihood.build_module``'s setup: wrap the raw options in + ``SectionOptions`` and instantiate the likelihood (whose ``__init__`` calls + ``build_data``). The one addition is reading ``csl_dir`` from the module + options to locate and import the upstream class before subclassing it. The + cosmosis import is deferred to here (call time) so importing this module never + requires cosmosis. + """ + from cosmosis.datablock import SectionOptions, option_section + + csl_dir = options.get_string(option_section, "csl_dir") + sacc_like = _import_upstream_sacc_like(csl_dir) + likelihood_class = _make_subclass(sacc_like) + return likelihood_class(SectionOptions(options)) + + +def execute(block, config): + """CosmoSIS ``execute`` — run the likelihood (mirrors ``build_module``).""" + likelihood_calculator = config + likelihood_calculator.do_likelihood(block) + return 0 + + +def cleanup(config): + """CosmoSIS ``cleanup`` — mirror of ``build_module``'s cleanup.""" + likelihood_calculator = config + likelihood_calculator.cleanup() diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py new file mode 100644 index 00000000..f4aa7d58 --- /dev/null +++ b/src/sp_validation/tests/test_blinding.py @@ -0,0 +1,1312 @@ +"""Tests for :mod:`sp_validation.blinding` — per-part Smokescreen blinding. + +Acceptance criteria AC1–AC9 of the blinding PRD, plus fast unit coverage of +the config/custody surface. Fast tests (no CCL import — the envelope +calibration, digest, commitment, fork-draw determinism, blind-init custody, +the merge alignment against a monkeypatched concealing factor, and the +assembly hash assertion) run in the default suite; the theory tests (fork + +CCL) are marked ``slow``; the derived-statistics tests additionally +``importorskip`` ``cosmo_numba``. + +All fixtures are synthetic and deterministic. Each blindable intermediate is +its own standalone part SACC (reporting ξ±, integration ξ±, pseudo-Cℓ), as in the +per-part-at-birth architecture; derived statistics (COSEBIs, pure-E/B) are +never stored in parts — they are computed downstream from the (blinded) integration +ξ± through the pipeline seams ``b_modes.cosebis_from_xi`` / +``b_modes.pure_eb_from_xi``, exactly as the pipeline does. The blinding path +hands the fork no data vector at all (``smokescreen.concealing_factor`` is a +pure theory difference), so fixture ξ± values are smooth synthetic templates — +no theory fill is needed to blind. +""" + +import json +import pathlib + +import numpy as np +import pytest + +from sp_validation import blinding as bd +from sp_validation import sacc_io as sio +from sp_validation.blinding_theory import TheoryConfig + +_NOLOG = lambda *a, **k: None # noqa: E731 + + +# --------------------------------------------------------------------------- # +# Synthetic part fixtures +# --------------------------------------------------------------------------- # +def _gauss_nz(z0, sigma, n=200): + 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 _reporting_theta(n=8): + return np.geomspace(5.0, 250.0, n) + + +def _integration_theta(n=80): + # The integration grid is the pure-E/B INTEGRATION grid, so it spans wider than + # the reporting range on both ends (production: ~0.08–300 arcmin). + return np.geomspace(0.1, 300.0, n) + + +def _xi_template(theta, k=0): + """Smooth synthetic ξ± for pair index ``k`` (no CCL needed).""" + theta = np.asarray(theta) + xip = 1e-4 * (1 + 0.1 * k) * (theta / 10.0) ** -0.6 + xim = 0.5e-4 * (1 + 0.1 * k) * (theta / 10.0) ** -0.9 + return xip, xim + + +def _b_mode_template(theta, amplitude): + """A smooth ξ_B(θ) template. B contributes +ξ_B to ξ+, −ξ_B to ξ−.""" + return amplitude * np.exp(-((np.log(np.asarray(theta) / 30.0)) ** 2) / 2.0) + + +def _nz_dict(nbins): + return {i: _gauss_nz(0.5 + 0.3 * i, 0.15 + 0.02 * i) for i in range(nbins)} + + +def _pairs(nbins): + return [(i, j) for i in range(nbins) for j in range(i, nbins)] + + +def make_xi_part(grid, nbins=1, b_amplitude=0.0): + """A standalone ξ± part SACC (one grid), synthetic values, eye covariance. + + ``b_amplitude`` injects a pure B-mode (+ξ_B to ξ+, −ξ_B to ξ−; + b_modes.py sign convention) — used on the integration part for AC4/AC9. + """ + theta = _reporting_theta() if grid == "reporting" else _integration_theta() + s = sio.new_sacc( + _nz_dict(nbins), metadata={"catalogue_version": "vTEST", "type": "mock"} + ) + blocks = [] + for k, (i, j) in enumerate(_pairs(nbins)): + xip, xim = _xi_template(theta, k) + xi_b = _b_mode_template(theta, b_amplitude) + sio.add_xi(s, (i, j), theta, xip + xi_b, xim - xi_b, grid=grid) + tr = sio._pair((i, j)) + idx = np.concatenate( + [ + s.indices(sio.XI_PLUS, tr, grid=grid), + s.indices(sio.XI_MINUS, tr, grid=grid), + ] + ) + blocks.append((idx, np.eye(len(idx)) * 1e-12)) + sio.assemble_covariance(s, blocks) + return s + + +def make_cl_part(nbins=1): + """A standalone pseudo-Cℓ part SACC (EE/BB/EB + bandpower windows).""" + s = sio.new_sacc( + _nz_dict(nbins), metadata={"catalogue_version": "vTEST", "type": "mock"} + ) + ell_eff = np.array([30.0, 80.0, 150.0, 280.0, 450.0]) + w_ell = np.arange(2, 501).astype(float) + w_mat = np.zeros((len(w_ell), len(ell_eff))) + for b, le in enumerate(ell_eff): + w_mat[:, b] = np.exp(-0.5 * ((w_ell - le) / 40.0) ** 2) + w_mat[:, b] /= w_mat[:, b].sum() + blocks = [] + for k, (i, j) in enumerate(_pairs(nbins)): + cl_ee = 1e-8 * (1 + 0.1 * k) * (ell_eff / 100.0) ** -1.2 + sio.add_pseudo_cl( + s, + (i, j), + ell_eff, + cl_ee, + np.zeros(5), + np.zeros(5), + window_ells=w_ell, + window_weights=w_mat, + ) + tr = sio._pair((i, j)) + idx = np.concatenate( + [s.indices(dt, tr) for dt in (sio.CL_EE, sio.CL_BB, sio.CL_EB)] + ) + blocks.append((idx, np.eye(len(idx)) * 1e-16)) + sio.assemble_covariance(s, blocks) + return s + + +def make_rho_part(): + """A standalone ρ/τ PSF-diagnostics part SACC — never blindable.""" + ctheta = _reporting_theta() + s = sio.new_sacc( + _nz_dict(1), metadata={"catalogue_version": "vTEST", "type": "mock"} + ) + blocks = [] + for k in range(2): + sio.add_rho( + s, k, ctheta, np.arange(len(ctheta)) * 1e-7, np.arange(len(ctheta)) * 2e-7 + ) + idx = np.concatenate( + [s.indices(sio.RHO_PLUS.format(k=k)), s.indices(sio.RHO_MINUS.format(k=k))] + ) + blocks.append((idx, np.eye(len(idx)) * 1e-18)) + sio.add_tau( + s, (0,), 0, ctheta, np.arange(len(ctheta)) * 3e-7, np.arange(len(ctheta)) * 4e-7 + ) + idx = np.concatenate( + [s.indices(sio.TAU_PLUS.format(k=0)), s.indices(sio.TAU_MINUS.format(k=0))] + ) + blocks.append((idx, np.eye(len(idx)) * 1e-18)) + sio.assemble_covariance(s, blocks) + return s + + +def make_parts(nbins=1, b_amplitude=0.0, with_rho=True): + """All intermediate parts of one catalogue version, keyed by name.""" + parts = { + "xi_reporting": make_xi_part("reporting", nbins), + "xi_integration": make_xi_part("integration", nbins, b_amplitude=b_amplitude), + "cl": make_cl_part(nbins), + } + if with_rho: + parts["rho_tau"] = make_rho_part() + return parts + + +def _derive_downstream(reporting_part, integration_part, nmodes=6): + """COSEBIs + pure-E/B the way the pipeline derives them downstream. + + COSEBIs from the integration ξ± (full-range scale cut); pure-E/B from the + measured reporting reporting ξ± + the integration integration ξ±, with the + edge-based bounds set to the reporting grid's span — the outermost + reporting point sits at tmax with no interior support and comes back + NaN (the AC9 boundary case). + """ + from sp_validation import b_modes + + theta_f, xip_f, xim_f = sio.get_xi(integration_part, (0, 0), grid="integration") + theta_c, xip_c, xim_c = sio.get_xi(reporting_part, (0, 0), grid="reporting") + En, Bn = b_modes.cosebis_from_xi( + theta_f, xip_f, xim_f, nmodes, scale_cut=(theta_f.min(), theta_f.max()) + ) + modes = b_modes.pure_eb_from_xi( + theta_c, + xip_c, + xim_c, + theta_f, + xip_f, + xim_f, + float(theta_c[0]), + float(theta_c[-1]), + ) + return En, Bn, modes + + +# --------------------------------------------------------------------------- # +# BlindingConfig: envelope calibration + digest (fast) +# --------------------------------------------------------------------------- # +def test_blinding_config_defaults(): + c = bd.BlindingConfig() + assert c.s8_half_width == 0.075 + assert c.omega_m_half_width == 0.1 + assert c.theory.S8 == 0.80 # fiducial TheoryConfig defaults + + +def test_blinding_config_overrides_fail_loud(): + c = bd.BlindingConfig.from_overrides({"s8_half_width": 0.05}) + assert c.s8_half_width == 0.05 + with pytest.raises(ValueError, match="unknown BlindingConfig fields"): + bd.BlindingConfig.from_overrides({"s8_half_width": 0.05, "bogus": 1}) + with pytest.raises(ValueError, match="unknown TheoryConfig fields"): + bd.BlindingConfig.from_overrides({"theory": {"nope": 1}}) + + +def test_blinding_config_is_frozen(): + with pytest.raises(Exception): + bd.BlindingConfig().s8_half_width = 0.2 + + +def test_envelope_calibration_maps_s8_box_to_ccl_halfwidths(): + """(S8, Ωm) half-widths → {sigma8, Omega_c} at the fiducial (exact forms).""" + c = bd.BlindingConfig() + shifts = c.shifts_dict() + assert set(shifts) == {"sigma8", "Omega_c"} + assert shifts["sigma8"] == pytest.approx( + c.s8_half_width / np.sqrt(c.theory.Omega_m / 0.3) + ) + assert shifts["Omega_c"] == c.omega_m_half_width + # every shift key must exist in the fiducial point (fork contract) + assert set(shifts) <= set(c.theory.ccl_params()) + + +def test_config_digest_stable_and_sensitive(): + """Canonical digest: byte-stable across runs, moves with every bound field.""" + c = bd.BlindingConfig() + assert c.config_digest() == bd.BlindingConfig().config_digest() + assert len(c.config_digest()) == 64 + assert bd.BlindingConfig(s8_half_width=0.05).config_digest() != c.config_digest() + assert ( + bd.BlindingConfig.from_overrides({"theory": {"S8": 0.79}}).config_digest() + != c.config_digest() + ) + # the P(k) recipe tokens are bound: a different halofit token = new digest + assert ( + bd.BlindingConfig.from_overrides( + {"theory": {"ccl_halofit_version": "takahashi"}} + ).config_digest() + != c.config_digest() + ) + # the Boltzmann backend (#280) is bound too — a different transfer + # function is a different P(k) path, so a different blind + assert ( + bd.BlindingConfig.from_overrides( + {"theory": {"transfer_function": "eisenstein_hu"}} + ).config_digest() + != c.config_digest() + ) + + +def test_theory_config_transfer_function_default_and_override(): + """#280: the Boltzmann backend is one config knob, CAMB by default (the + inference pipeline's Boltzmann code), overridable like any other field.""" + assert TheoryConfig().transfer_function == "boltzmann_camb" + cfg = TheoryConfig.from_overrides({"transfer_function": "eisenstein_hu"}) + assert cfg.transfer_function == "eisenstein_hu" + + +def test_config_digest_int_float_canonical(): + """One physical cosmology has one digest, regardless of int-vs-float literals. + + Configs come from YAML/CLI/humans, so a field can arrive as ``-1`` (int) or + ``-1.0`` (float). The digest must depend on the numeric *value*, not the + literal's Python type — otherwise an int-vs-float mismatch between the blind + and a later unblind would raise "config digest mismatch" and deny a + legitimate unblind. Every declared-float field must be canonical this way. + """ + for field, int_val, float_val in [ + ("w0", -1, -1.0), + ("wa", 0, 0.0), + ("Omega_m", 1, 1.0), + ("m_nu", 0, 0.0), + ("S8", 1, 1.0), + ]: + assert ( + bd.BlindingConfig.from_overrides( + {"theory": {field: int_val}} + ).config_digest() + == bd.BlindingConfig.from_overrides( + {"theory": {field: float_val}} + ).config_digest() + ), f"int-vs-float digest split on theory.{field}" + # the envelope half-widths (BlindingConfig's own float fields) too + assert ( + bd.BlindingConfig.from_overrides({"s8_half_width": 1}).config_digest() + == bd.BlindingConfig.from_overrides({"s8_half_width": 1.0}).config_digest() + ) + # and a full round-trip: an all-int override matches the float default digest + assert ( + bd.BlindingConfig.from_overrides( + {"theory": {"w0": -1, "Omega_m": 0, "wa": 0}} + ).config_digest() + == bd.BlindingConfig.from_overrides( + {"theory": {"w0": -1.0, "Omega_m": 0.0, "wa": 0.0}} + ).config_digest() + ) + + +def test_theory_config_ccl_params_exact_keyset(): + """ccl_params() carries exactly the contracted keys — nothing rides along.""" + params = TheoryConfig().ccl_params() + assert set(params) == { + "Omega_c", + "Omega_b", + "h", + "n_s", + "sigma8", + "m_nu", + "mass_split", + "w0", + "wa", + "Neff", + "T_CMB", + } + assert params["sigma8"] == pytest.approx(0.80) # S8=0.80 at Ωm=0.30 + assert params["Neff"] == 3.046 and params["T_CMB"] == 2.7255 + + +# --------------------------------------------------------------------------- # +# Commitment + the fork's draw (fast; smokescreen import is light) +# --------------------------------------------------------------------------- # +def test_commitment_is_the_forks_domain_separated_digest(): + """One definition of the commitment, and it is the fork's.""" + import smokescreen + + seed = "the-secret" + assert bd.seed_commitment(seed) == smokescreen.seed_commitment(seed) + assert bd.seed_commitment("right") != bd.seed_commitment("wrong") + + +def test_commitment_does_not_embed_the_rng_seed(): + """The published commitment must not carry the effective RNG seed. + + The fork derives the base seed for its per-key RNG from the *undomained* + sha256 of the seed string, taking the digest's first 8 bytes. A commitment + hashed over the bare seed would therefore publish that base seed verbatim in + its first 16 hex characters, and anyone holding the (public) commitment and + the (public) fiducial config could redraw the hidden cosmology and subtract + the blind. The domain prefix is what breaks that identity — this is the + regression guard for it. + """ + import secrets + + from smokescreen.param_shifts import _normalize_seed + + # The third seed is drawn exactly as blind_init draws a production one. + for seed in ("my_secret_seed", "the-secret", secrets.token_hex(16)): + commitment = bd.seed_commitment(seed) + assert int(commitment[:16], 16) != _normalize_seed(seed) + + +def test_hidden_params_deterministic_and_in_envelope(): + """Same (seed, config) ⇒ same hidden point; draws respect the envelope.""" + import secrets + + c = bd.BlindingConfig() + fid = c.theory.ccl_params() + h1, h2 = bd.hidden_params("a-seed", c), bd.hidden_params("a-seed", c) + assert h1 == h2 + assert bd.hidden_params("другой", c) != h1 + shifts = c.shifts_dict() + for _ in range(50): + h = bd.hidden_params(secrets.token_hex(8), c) + for key, half in shifts.items(): + assert abs(h[key] - fid[key]) <= half + # only the enveloped keys move + assert all(h[k] == fid[k] for k in fid if k not in shifts) + + +def test_hidden_params_no_global_rng_state(): + """The fork's draw uses a local RNG — global numpy state is untouched.""" + np.random.seed(0) + before = np.random.get_state()[1].copy() + bd.hidden_params("whatever", bd.BlindingConfig()) + assert np.array_equal(before, np.random.get_state()[1]) + + +# --------------------------------------------------------------------------- # +# Per-part merge + provenance with a monkeypatched factor (fast — no CCL) +# --------------------------------------------------------------------------- # +def _patch_constant_factor(monkeypatch, value=1e-6): + def fake(part, indices, factory, config, seed): + return np.arange(len(indices), dtype=float) * value + value + + monkeypatch.setattr(bd, "_concealing_factor", fake) + return fake + + +def test_merge_places_shift_at_recorded_indices_only(monkeypatch): + """AC5 (merge half), per part: the shift lands exactly on the blindable + rows, in stored order; covariance, n(z), and every tag are untouched.""" + _patch_constant_factor(monkeypatch) + for name, part in make_parts(nbins=2, with_rho=False).items(): + orig = np.array(part.mean) + orig_cov = part.covariance.dense.copy() + orig_nz = part.tracers["source_0"].nz.copy() + + blinded = bd.blind_sacc(part, "seed", log=_NOLOG) + + blocks = bd._blindable_blocks(part) + assert len(blocks) == 1, f"{name}: a part carries exactly one block" + shifted = np.zeros(len(orig), dtype=bool) + for _, indices, _ in blocks: + expected = orig[indices] + (np.arange(len(indices)) * 1e-6 + 1e-6) + assert np.array_equal(np.array(blinded.mean)[indices], expected) + shifted[indices] = True + assert np.array_equal(np.array(blinded.mean)[~shifted], orig[~shifted]) + assert np.array_equal(blinded.covariance.dense, orig_cov) + assert np.array_equal(blinded.tracers["source_0"].nz, orig_nz) + # row-order preservation: type/tracers/tags sequence is bitwise unchanged + for a, b in zip(part.data, blinded.data): + assert a.data_type == b.data_type + assert a.tracers == b.tracers + assert a.tags == b.tags + + +def test_provenance_metadata_contract(monkeypatch): + """Blinded parts carry concealed/blind/commitment/digest; no seed key.""" + _patch_constant_factor(monkeypatch) + s = make_xi_part("reporting") + s.metadata["seed_smokescreen"] = "leaked!" # must be stripped + c = bd.BlindingConfig() + blinded = bd.blind_sacc(s, "seed", config=c, label="B", log=_NOLOG) + assert blinded.metadata["concealed"] is True + assert blinded.metadata["blind"] == "B" + assert blinded.metadata["blind_commitment"] == bd.seed_commitment("seed") + assert blinded.metadata["blind_config_digest"] == c.config_digest() + assert blinded.metadata["blind_draw_scheme"] == bd.draw_scheme() + assert "seed_smokescreen" not in blinded.metadata + assert blinded.metadata["catalogue_version"] == "vTEST" + + +def test_blind_refuses_double_blind(monkeypatch): + _patch_constant_factor(monkeypatch) + s = make_xi_part("reporting") + blinded = bd.blind_sacc(s, "seed", log=_NOLOG) + with pytest.raises(ValueError, match="already concealed"): + bd.blind_sacc(blinded, "seed2", log=_NOLOG) + + +def test_blind_refuses_non_blindable_part(): + """A ρ/τ diagnostic part must never see a blind call — loud refusal.""" + with pytest.raises(ValueError, match="no blindable block"): + bd.blind_sacc(make_rho_part(), "seed", log=_NOLOG) + + +def test_unblind_fails_closed_on_wrong_seed_or_config(monkeypatch): + """AC6 (in-memory half): wrong seed and wrong config both refuse loudly.""" + _patch_constant_factor(monkeypatch) + s = make_xi_part("reporting") + blinded = bd.blind_sacc(s, "right-seed", log=_NOLOG) + with pytest.raises(ValueError, match="blind_commitment"): + bd.unblind_sacc(blinded, "wrong-seed", log=_NOLOG) + with pytest.raises(ValueError, match="blind_config_digest"): + bd.unblind_sacc( + blinded, + "right-seed", + config=bd.BlindingConfig(s8_half_width=0.01), + log=_NOLOG, + ) + with pytest.raises(ValueError, match="not concealed"): + bd.unblind_sacc(s, "right-seed", log=_NOLOG) + + +# --------------------------------------------------------------------------- # +# The vector core: full-length theory in, block slice out (fast — fake backend) +# --------------------------------------------------------------------------- # +def _fake_factory(indices, hole=None): + """A backend that fills ``indices`` with ``sigma8 * arange`` and nothing else. + + Mimics the real backends' contract — a block's ``theory_fn`` reads its + layout off the SACC it is given and returns a full-length vector, NaN on + every row outside its own block — without importing CCL. ``hole`` leaves + one row of the block itself unfilled. + """ + + def factory(s, theory): + def theory_fn(params): + out = np.full(len(s.mean), np.nan) + out[indices] = params["sigma8"] * np.arange(len(indices)) + if hole is not None: + out[hole] = np.nan + return out + + return theory_fn + + return factory + + +def test_concealing_factor_slices_its_own_block_from_a_full_length_vector(): + """The factor is the theory difference at the block's rows, and the rows + the backend does not fill are never read. + + This is the shape of the whole path after the sub-SACC carving came out: + the backend is driven off the assembled SACC directly, returns a + full-length vector that is NaN everywhere but its own block, and + ``_concealing_factor`` returns exactly that block's rows. Checked against + :func:`hidden_params`, which reaches the same hidden point by an + independent route. + """ + ordered = ["xi_reporting", "cl", "rho_tau", "xi_integration"] + s = sio.gather([make_parts(nbins=1)[k] for k in ordered]) + blocks = bd._blindable_blocks(s) + assert len(blocks) == 3, "assembled file carries all three blindable blocks" + + cfg = bd.BlindingConfig() + delta = bd.hidden_params("seed", cfg)["sigma8"] - cfg.theory.ccl_params()["sigma8"] + assert delta != 0.0 + for _, indices, _ in blocks: + factor = bd._concealing_factor(s, indices, _fake_factory(indices), cfg, "seed") + assert factor.shape == (len(indices),) + assert np.allclose(factor, delta * np.arange(len(indices))) + + +def test_concealing_factor_refuses_a_row_its_backend_cannot_fill(): + """A NaN on a row the block *claims* is a layout the backend cannot cover. + + Slicing to the block drops the NaNs outside it by construction; this is + the converse guard, and it has to be explicit — without it a row the + backend silently skipped would be shifted by NaN, destroying that point + with no error anywhere. + """ + part = make_xi_part("reporting") + ((_, indices, _),) = bd._blindable_blocks(part) + with pytest.raises(ValueError, match="unfilled"): + bd._concealing_factor( + part, + indices, + _fake_factory(indices, hole=indices[0]), + bd.BlindingConfig(), + "seed", + ) + + +# --------------------------------------------------------------------------- # +# Draw-scheme binding: the blind is (seed, config, draw semantics) +# --------------------------------------------------------------------------- # +def test_draw_scheme_is_the_installed_fork_constant(): + """The recorded scheme is read from Smokescreen, not hardcoded here.""" + import smokescreen + + assert bd.draw_scheme() == int(smokescreen.DRAW_SCHEME) + assert isinstance(bd.draw_scheme(), int) + + +def test_assert_draw_scheme_message_names_both_versions(monkeypatch): + """A scheme mismatch says which scheme made the blind and which is installed.""" + monkeypatch.setattr(bd, "draw_scheme", lambda: 2) + bd._assert_draw_scheme(2, "the blind") # matching scheme is silent + with pytest.raises(ValueError, match=r"DRAW_SCHEME=1.*DRAW_SCHEME=2"): + bd._assert_draw_scheme(1, "the blind") + with pytest.raises(ValueError, match="no draw-scheme record"): + bd._assert_draw_scheme(None, "the blind") + + +def test_unblind_refuses_a_blind_drawn_under_another_scheme(monkeypatch): + """The finding this closes: a blind added under one draw scheme and + subtracted under another silently produces a wrong data vector, because + the seed hash, the config digest and the escrow check all still pass. The + scheme must be part of the custody state, checked before any subtraction. + """ + _patch_constant_factor(monkeypatch) + blinded = bd.blind_sacc(make_xi_part("reporting"), "seed", log=_NOLOG) + assert blinded.metadata["blind_draw_scheme"] == bd.draw_scheme() + + # everything else about this file is valid — only the draw semantics moved + monkeypatch.setattr( + bd, "draw_scheme", lambda: blinded.metadata["blind_draw_scheme"] + 1 + ) + with pytest.raises(ValueError, match="DRAW_SCHEME"): + bd.unblind_sacc(blinded, "seed", log=_NOLOG) + + +def test_unblind_refuses_a_file_predating_scheme_binding(monkeypatch): + """A blinded file with no scheme record fails closed, not open.""" + _patch_constant_factor(monkeypatch) + blinded = bd.blind_sacc(make_xi_part("reporting"), "seed", log=_NOLOG) + del blinded.metadata["blind_draw_scheme"] + with pytest.raises(ValueError, match="no draw-scheme record"): + bd.unblind_sacc(blinded, "seed", log=_NOLOG) + + +def test_unblind_strips_the_scheme_stamp_with_the_rest(monkeypatch): + _patch_constant_factor(monkeypatch) + blinded = bd.blind_sacc(make_xi_part("reporting"), "seed", log=_NOLOG) + part = bd.unblind_sacc(blinded, "seed", log=_NOLOG) + assert "blind_draw_scheme" not in part.metadata + + +def test_read_seed_fails_closed_on_scheme_drift(tmp_path, monkeypatch): + """blind_part reads the seed through _read_seed, so a scheme change between + blinding part 1 and part 2 is caught before the second part is shifted.""" + bd.blind_init(str(tmp_path), log=_NOLOG) + monkeypatch.setattr(bd, "draw_scheme", lambda: 99) + with pytest.raises(ValueError, match="DRAW_SCHEME"): + bd._read_seed(str(tmp_path), bd.BlindingConfig()) + + +def test_stamp_passthrough_carries_the_committed_scheme(tmp_path, monkeypatch): + """A pass-through part inherits the blind's scheme, and cannot be stamped + from an install that draws differently.""" + paths = bd.blind_init(str(tmp_path), log=_NOLOG) + s = bd.stamp_concealed_passthrough(make_rho_part(), paths["commitment"]) + assert s.metadata["blind_draw_scheme"] == bd.draw_scheme() + + monkeypatch.setattr(bd, "draw_scheme", lambda: 99) + with pytest.raises(ValueError, match="DRAW_SCHEME"): + bd.stamp_concealed_passthrough(make_rho_part(), paths["commitment"]) + + +def test_assert_consistent_blind_refuses_divergent_or_foreign_schemes(monkeypatch): + """Parts blinded under different schemes never assemble; nor does a set + that agrees with itself but not with the installed fork.""" + parts = make_parts(nbins=1, with_rho=False) + for p in parts.values(): + _stamp(p) + parts["cl"].metadata["blind_draw_scheme"] = bd.draw_scheme() + 1 + with pytest.raises(ValueError, match="different blind commitments"): + bd.assert_consistent_blind(list(parts.values())) + + parts = make_parts(nbins=1, with_rho=False) + for p in parts.values(): + _stamp(p) + p.metadata["blind_draw_scheme"] = bd.draw_scheme() + 1 + with pytest.raises(ValueError, match="DRAW_SCHEME"): + bd.assert_consistent_blind(list(parts.values())) + + +# --------------------------------------------------------------------------- # +# blind-init custody + assembly hash assertion (fast — encryption only) +# --------------------------------------------------------------------------- # +def test_blind_init_writes_commitment_and_encrypted_bundle_only(tmp_path): + """AC6 (init): commitment.json + encrypted bundle; never a plaintext seed.""" + paths = bd.blind_init(str(tmp_path), log=_NOLOG) + with open(paths["commitment"], encoding="utf-8") as f: + commitment = json.load(f) + assert set(commitment) == { + "label", + "seed_commitment", + "config_digest", + "draw_scheme", + } + assert len(commitment["seed_commitment"]) == 64 + assert commitment["config_digest"] == bd.BlindingConfig().config_digest() + # exactly the three custody outputs, no plaintext bundle + assert {p.name for p in tmp_path.iterdir()} == { + "commitment.json", + "blind_seed.encrpt", + "blind_seed.key", + } + # the decrypted seed matches the public commitment + bundle = bd._read_encrypted_json(paths["bundle"], paths["key"]) + assert bd.seed_commitment(bundle["seed"]) == commitment["seed_commitment"] + # one-shot custody: a second init in the same dir refuses + with pytest.raises(FileExistsError, match="refusing to overwrite"): + bd.blind_init(str(tmp_path), log=_NOLOG) + + +def test_read_seed_fails_closed_on_drifted_config(tmp_path): + bd.blind_init(str(tmp_path), log=_NOLOG) + with pytest.raises(ValueError, match="config digest"): + bd._read_seed(str(tmp_path), bd.BlindingConfig(s8_half_width=0.01)) + + +def _stamp(s, seed="s", label="A", config=None): + bd._stamp_provenance( + s, + bd.seed_commitment(seed), + label, + (config or bd.BlindingConfig()).config_digest(), + ) + return s + + +def test_assert_consistent_blind_shared_stamp(): + """One commitment across all blindable parts ⇒ the shared stamp returns; + ρ/τ parts are exempt from the assertion.""" + parts = make_parts(nbins=1) + for name in ("xi_reporting", "xi_integration", "cl"): + _stamp(parts[name]) + stamp = bd.assert_consistent_blind(list(parts.values())) + assert stamp == { + "concealed": True, + "blind": "A", + "blind_commitment": bd.seed_commitment("s"), + "blind_config_digest": bd.BlindingConfig().config_digest(), + "blind_draw_scheme": bd.draw_scheme(), + } + + +def test_assert_consistent_blind_fails_closed(): + """AC6 (assembly): mismatched commitments and mixed states both refuse.""" + parts = make_parts(nbins=1, with_rho=False) + _stamp(parts["xi_reporting"], seed="one") + _stamp(parts["xi_integration"], seed="one") + _stamp(parts["cl"], seed="two") # different seed ⇒ different commitment + with pytest.raises(ValueError, match="different blind commitments"): + bd.assert_consistent_blind(list(parts.values())) + + parts = make_parts(nbins=1, with_rho=False) + _stamp(parts["xi_reporting"]) # blinded beside plaintext blindable parts + with pytest.raises(ValueError, match="mixed"): + bd.assert_consistent_blind(list(parts.values())) + + +def test_assert_consistent_blind_all_plaintext_is_none(): + """A declared-mock plaintext assembly (nothing blinded) asserts nothing.""" + assert bd.assert_consistent_blind(list(make_parts().values())) is None + + +def test_assert_consistent_blind_unconcealed_data_fails_closed(): + """PRD §4 "Mocks vs data": an unconcealed blindable part may assemble only + if its metadata declares ``type == "mock"`` — an unconcealed ``data`` part, + or one missing the tag, fails closed (skipping the blind can never + silently expose real data). Mixed concealed/plaintext still fails on the + mixed guard regardless of type.""" + parts = make_parts(nbins=1, with_rho=False) + parts["xi_integration"].metadata["type"] = "data" + with pytest.raises(ValueError, match="type"): + bd.assert_consistent_blind(list(parts.values())) + + parts = make_parts(nbins=1, with_rho=False) + del parts["cl"].metadata["type"] # missing tag counts as not-a-mock + with pytest.raises(ValueError, match=""): + bd.assert_consistent_blind(list(parts.values())) + + # concealed data parts assemble integration (that is the whole point of the blind) + parts = make_parts(nbins=1, with_rho=False) + for p in parts.values(): + p.metadata["type"] = "data" + _stamp(p) + assert bd.assert_consistent_blind(list(parts.values()))["concealed"] is True + + # mixed concealed/plaintext fails on the mixed guard even for mocks + parts = make_parts(nbins=1, with_rho=False) + _stamp(parts["xi_reporting"]) + with pytest.raises(ValueError, match="mixed"): + bd.assert_consistent_blind(list(parts.values())) + + +def test_gather_fails_closed_on_unconcealed_data_part(): + """The type guard reaches the terminal gather surface too.""" + parts = make_parts(nbins=1, with_rho=False) + parts["xi_reporting"].metadata["type"] = "data" + with pytest.raises(ValueError, match="type"): + sio.gather(list(parts.values())) + + +def test_assert_consistent_blind_differing_labels_warn_not_fail(): + """Same seed+config, different --label ⇒ assemble cleanly with a warning. + + The label is provenance, not custody state: parts blinded under one blind + but tagged with different labels must not be misread as different blinds. + The assembly succeeds (keyed on commitment+digest); a distinct warning + surfaces the label divergence rather than a false "different commitments". + """ + parts = make_parts(nbins=1, with_rho=False) + _stamp(parts["xi_reporting"], seed="one", label="A") + _stamp(parts["xi_integration"], seed="one", label="A") + _stamp(parts["cl"], seed="one", label="B") # same blind, different label + with pytest.warns(UserWarning, match="different labels"): + stamp = bd.assert_consistent_blind(list(parts.values())) + assert stamp["blind_commitment"] == bd.seed_commitment("one") + assert stamp["blind"] == "A" # deterministic: sorted-first label + + +def test_gather_assembles_parts_and_stamps_blind(): + """The sacc_io gather combines parts (points, tags, covariance blocks), + calls the assembly assertion, and stamps the shared blind.""" + parts = make_parts(nbins=1) + for name in ("xi_reporting", "xi_integration", "cl"): + _stamp(parts[name]) + ordered = [parts[k] for k in ("xi_reporting", "cl", "rho_tau", "xi_integration")] + s = sio.gather(ordered, metadata={"catalogue_version": "vTEST", "type": "mock"}) + + assert len(s.mean) == sum(len(p.mean) for p in ordered) + assert np.array_equal(np.array(s.mean), np.concatenate([p.mean for p in ordered])) + # covariance: block-diagonal of the parts, in order + cursor = 0 + for p in ordered: + n = len(p.mean) + assert np.array_equal( + s.covariance.dense[cursor : cursor + n, cursor : cursor + n], + p.covariance.dense, + ) + cursor += n + # the integration rows are addressable by the grid tag in the assembled file + assert len(s.indices(sio.XI_PLUS, grid="integration")) == len( + parts["xi_integration"].indices(sio.XI_PLUS, grid="integration") + ) + # bandpower windows survive assembly + assert s.get_bandpower_windows(s.indices(sio.CL_EE)) is not None + # blind stamp on the assembled file + assert s.metadata["concealed"] is True + assert s.metadata["blind_commitment"] == bd.seed_commitment("s") + assert s.metadata["catalogue_version"] == "vTEST" + + +def test_gather_fails_closed_on_mismatched_blinds(): + parts = make_parts(nbins=1, with_rho=False) + _stamp(parts["xi_reporting"], seed="one") + _stamp(parts["xi_integration"], seed="two") + _stamp(parts["cl"], seed="one") + with pytest.raises(ValueError, match="different blind commitments"): + sio.gather(list(parts.values())) + + +def test_cli_blind_init_refuses_existing_state(tmp_path): + """The blind-init CLI refuses to overwrite a previous blind's state.""" + import importlib.util + + script = ( + pathlib.Path(__file__).resolve().parents[3] / "scripts" / "blind_data_vector.py" + ) + spec = importlib.util.spec_from_file_location("_blind_cli", script) + cli = importlib.util.module_from_spec(spec) + spec.loader.exec_module(cli) + + (tmp_path / "commitment.json").write_text("{}") + with pytest.raises(SystemExit, match="refusing to overwrite"): + cli.main(["blind-init", str(tmp_path)]) + + +# --------------------------------------------------------------------------- # +# AC2 + AC3 + AC7: the shift itself (slow — fork + CCL) +# --------------------------------------------------------------------------- # +@pytest.mark.slow +@pytest.mark.parametrize( + "transfer_function", + ["boltzmann_camb", "eisenstein_hu"], + ids=["default-camb", "non-default-eh"], +) +def test_ac2_on_file_shift_equals_theory_difference_per_part(transfer_function): + """AC2: per-row shift on each blinded part == theory_fn(hidden) − + theory_fn(fiducial), hidden recovered by re-running the fork's draw — + for all three parts, on a two-bin fixture; and the recovered hidden + cosmology is identical across the three parts (one seed → one hidden). + + Scope: this verifies fork-draw recovery + placement (the shift on the + file is exactly what re-running the same backend at the recovered + hidden/fiducial points predicts, at the recorded rows). It is NOT a + backend-correctness test: both sides run the identical ``theory_fn``, so + any wrong-cosmology dependence cancels and a wrong backend would still + pass here. Backend correctness is carried by AC3. + + Parametrized over the Boltzmann backend (#280): the default CAMB route + and one non-default (Eisenstein–Hu — cheap, no CAMB run) both thread the + same ``transfer_function`` knob through all three theory backends.""" + cfg = bd.BlindingConfig.from_overrides( + {"theory": {"transfer_function": transfer_function}} + ) + seed = "ac2-seed" + parts = make_parts(nbins=2, with_rho=False) + + hiddens, worst = [], 0.0 + for name, part in parts.items(): + blinded = bd.blind_sacc(part, seed, config=cfg, log=_NOLOG) + hidden = bd.hidden_params(seed, cfg) # recovered per part + hiddens.append(hidden) + fiducial = cfg.theory.ccl_params() + ((block_name, indices, factory),) = bd._blindable_blocks(part) + # Independent of the blinding path: the factory is driven directly off + # the part, at the hidden point recovered by hidden_params, and the two + # theory vectors are differenced here rather than by the fork. The + # factory fills only its own block, so slice to it. + theory = factory(part, cfg.theory) + expected = (theory(hidden) - theory(fiducial))[indices] + actual = np.array(blinded.mean)[indices] - np.array(part.mean)[indices] + gap = np.max(np.abs(actual - expected)) + scale = np.max(np.abs(expected)) + worst = max(worst, gap / scale) + assert gap <= 1e-10 * max(scale, 1e-30), f"{name}/{block_name}: |Δ|={gap:.3e}" + # one seed → one hidden cosmology across all parts + assert hiddens[0] == hiddens[1] == hiddens[2] + print(f"\nAC2 max relative shift mismatch across parts: {worst:.3e}") + + +@pytest.mark.slow +def test_ac3_cross_backend_against_independent_ccl_reference(): + """AC3: the realized shift matches an independently written direct-CCL + reference (self-contained here; does not touch the blinding backends).""" + import pyccl as ccl + + cfg = bd.BlindingConfig() + seed = "ac3-seed" + reporting = make_xi_part("reporting", nbins=2) + cl_part = make_cl_part(nbins=2) + blinded_xi = bd.blind_sacc(reporting, seed, config=cfg, log=_NOLOG) + blinded_cl = bd.blind_sacc(cl_part, seed, config=cfg, log=_NOLOG) + hidden = bd.hidden_params(seed, cfg) + fiducial = cfg.theory.ccl_params() + + # ----- independent reference (from scratch; same fixture n(z), θ, ℓ) ---- + def ref_cosmo(params): + return ccl.Cosmology( + **params, + matter_power_spectrum="camb", + extra_parameters={ + "camb": {"halofit_version": "mead2020_feedback", "HMCode_logT_AGN": 7.5} + }, + ) + + ell = np.unique( + np.concatenate([np.arange(2, 50), np.geomspace(50, 6e4, 200)]).astype(float) + ) + + def ref_xi(params, nz_i, nz_j, theta): + cosmo = ref_cosmo(params) + ti = ccl.WeakLensingTracer(cosmo, dndz=nz_i) + tj = ccl.WeakLensingTracer(cosmo, dndz=nz_j) + cl = ccl.angular_cl(cosmo, ti, tj, ell) + xip = ccl.correlation(cosmo, ell=ell, C_ell=cl, theta=theta / 60.0, type="GG+") + xim = ccl.correlation(cosmo, ell=ell, C_ell=cl, theta=theta / 60.0, type="GG-") + return xip, xim + + worst = 0.0 + for i, j in bd.xi_pairs(reporting, "reporting"): + tr = sio._pair((i, j)) + theta = sio._tag(reporting, sio.XI_PLUS, tr, "theta", grid="reporting") + nz_i, nz_j = sio.get_nz(reporting, i), sio.get_nz(reporting, j) + xip_h, xim_h = ref_xi(hidden, nz_i, nz_j, theta) + xip_f, xim_f = ref_xi(fiducial, nz_i, nz_j, theta) + for dt, ref_shift in ( + (sio.XI_PLUS, xip_h - xip_f), + (sio.XI_MINUS, xim_h - xim_f), + ): + idx = reporting.indices(dt, tr, grid="reporting") + realized = np.array(blinded_xi.mean)[idx] - np.array(reporting.mean)[idx] + worst = max(worst, np.max(np.abs(realized - ref_shift))) + print(f"\nAC3 max |realized − independent reference| (reporting ξ±): {worst:.3e}") + assert worst < 1e-8 # observed ~1e-10; factor magnitudes ~1e-6 + + # pseudo-Cℓ: W @ ΔCℓ_EE against the same independent reference + for i, j in bd.cl_pairs(cl_part): + tr = sio._pair((i, j)) + idx = cl_part.indices(sio.CL_EE, tr) + window = cl_part.get_bandpower_windows(idx) + w_ell = np.asarray(window.values, dtype=float) + w_mat = np.asarray(window.weight, dtype=float) + nz_i, nz_j = sio.get_nz(cl_part, i), sio.get_nz(cl_part, j) + + def ref_cl(params): + cosmo = ref_cosmo(params) + ti = ccl.WeakLensingTracer(cosmo, dndz=nz_i) + tj = ccl.WeakLensingTracer(cosmo, dndz=nz_j) + return ccl.angular_cl(cosmo, ti, tj, w_ell) + + ref_shift = w_mat.T @ (ref_cl(hidden) - ref_cl(fiducial)) + realized = np.array(blinded_cl.mean)[idx] - np.array(cl_part.mean)[idx] + assert np.max(np.abs(realized - ref_shift)) < 1e-12 # bandpowers ~1e-9 + + +@pytest.mark.slow +def test_ac7_reproducibility_same_seed_same_shift(): + """AC7: two blind runs of the same part with the same (seed, config) + produce identical shifts; a different seed produces a different one.""" + part = make_xi_part("reporting") + b1 = bd.blind_sacc(part, "repro-seed", log=_NOLOG) + b2 = bd.blind_sacc(part, "repro-seed", log=_NOLOG) + assert np.array_equal(np.array(b1.mean), np.array(b2.mean)) + b3 = bd.blind_sacc(part, "other-seed", log=_NOLOG) + assert not np.array_equal(np.array(b1.mean), np.array(b3.mean)) + + +# --------------------------------------------------------------------------- # +# AC1, AC4, AC5, AC9: born-blinded derived statistics (slow + cosmo_numba) +# --------------------------------------------------------------------------- # +@pytest.mark.slow +def test_ac1_zero_shift_is_identity(): + """AC1: a zero envelope reproduces every part exactly, and the integration part + run downstream through the b_modes seams yields COSEBIs and pure-E/B + identical to the unblinded run — the per-part plumbing and the + born-blinded derivation path are the identity at zero shift.""" + pytest.importorskip("cosmo_numba") + zero = bd.BlindingConfig(s8_half_width=0.0, omega_m_half_width=0.0) + parts = make_parts(nbins=1, with_rho=False) + blinded = { + name: bd.blind_sacc(p, "any-seed", config=zero, log=_NOLOG) + for name, p in parts.items() + } + for name, part in parts.items(): + assert np.array_equal(np.array(part.mean), np.array(blinded[name].mean)), ( + f"zero-shift blind changed {name}" + ) + # derived statistics downstream: identical inputs ⇒ identical numbers + En_t, Bn_t, modes_t = _derive_downstream( + parts["xi_reporting"], parts["xi_integration"] + ) + En_b, Bn_b, modes_b = _derive_downstream( + blinded["xi_reporting"], blinded["xi_integration"] + ) + assert np.array_equal(En_t, En_b) and np.array_equal(Bn_t, Bn_b) + for key in modes_t: + t, b = modes_t[key], modes_b[key] + both_nan = np.isnan(t) & np.isnan(b) + assert np.array_equal(t[~both_nan], b[~both_nan]), key + assert np.array_equal(np.isnan(t), np.isnan(b)), key + + +@pytest.mark.slow +def test_ac4_b_mode_invariance_and_leakage_floor(): + """AC4: the ΔBₙ induced by deriving B-modes from the blinded integration ξ± is + independent of the injected B amplitude (a fixed absolute E→B leakage + offset, not fractional). The magnitude is measured and reported, never + asserted against a constant.""" + pytest.importorskip("cosmo_numba") + seed = "ac4-seed" + deltas, reports = [], [] + for amp in (2e-6, 2e-5): + reporting = make_xi_part("reporting") + integration = make_xi_part("integration", b_amplitude=amp) + blinded_reporting = bd.blind_sacc(reporting, seed, log=_NOLOG) + blinded_integration = bd.blind_sacc(integration, seed, log=_NOLOG) + _, Bn_t, modes_t = _derive_downstream(reporting, integration) + _, Bn_b, modes_b = _derive_downstream(blinded_reporting, blinded_integration) + d_bn = Bn_b - Bn_t + d_xib = modes_b["xip_B"] - modes_t["xip_B"] + finite = np.isfinite(d_xib) + deltas.append((d_bn, d_xib[finite])) + reports.append( + f"B={amp:.0e}: max|ΔBₙ|={np.max(np.abs(d_bn)):.3e} " + f"(ΔBₙ/Bₙ={np.max(np.abs(d_bn)) / np.max(np.abs(Bn_t)):.2e}), " + f"max|Δξ+_B|={np.max(np.abs(d_xib[finite])):.3e} " + f"({np.max(np.abs(d_xib[finite])) / amp:.2%} of injected B)" + ) + print("\nAC4 " + "\nAC4 ".join(reports)) + + (d_bn_1, d_xib_1), (d_bn_2, d_xib_2) = deltas + scale = max(np.max(np.abs(d_bn_1)), 1e-30) + gap = np.max(np.abs(d_bn_1 - d_bn_2)) + print( + f"AC4 ΔBₙ amplitude-independence: max|ΔBₙ(2e-6) − ΔBₙ(2e-5)| = " + f"{gap:.3e} ({gap / scale:.2e} of |ΔBₙ|)" + ) + assert gap <= 1e-9 * scale + 1e-24, ( + "ΔBₙ depends on the injected B amplitude — the shift is not pure E" + ) + # The pure-ξ_B leakage is also amplitude-independent, but only to the + # adaptive-quadrature floor: cosmo_numba's Schneider integrals subdivide + # adaptively, so the estimator is not bit-linear in its inputs and the + # two runs differ at a small fraction of the (tiny) leakage itself. The + # COSEBIs assertion above carries the exact-identity criterion; this one + # bounds the quadrature wobble. + scale_x = max(np.max(np.abs(d_xib_1)), 1e-30) + gap_x = np.max(np.abs(d_xib_1 - d_xib_2)) + print( + f"AC4 Δξ+_B amplitude-independence: {gap_x:.3e} " + f"({gap_x / scale_x:.2e} of the leakage)" + ) + assert gap_x <= 0.05 * scale_x + + +@pytest.mark.slow +def test_ac5_untouched_blocks_and_row_order(): + """AC5: each part's covariance byte-identical; the Cℓ BB/EB rows and the + ρ/τ part never blinded; every part's shifted rows land at their original + within-part indices (order-preservation).""" + parts = make_parts(nbins=1) + seed = "ac5-seed" + blinded = { + name: bd.blind_sacc(parts[name], seed, log=_NOLOG) + for name in ("xi_reporting", "xi_integration", "cl") + } + for name, b in blinded.items(): + part = parts[name] + assert np.array_equal(b.covariance.dense, part.covariance.dense), name + # row order: the identity of every row (type/tracers/tags) unchanged + for a, c in zip(part.data, b.data): + assert (a.data_type, a.tracers, a.tags) == (c.data_type, c.tracers, c.tags) + # BB/EB rows of the Cℓ part untouched (pure E-mode shift) + for dt in (sio.CL_BB, sio.CL_EB): + idx = parts["cl"].indices(dt) + assert np.array_equal( + np.array(blinded["cl"].mean)[idx], np.array(parts["cl"].mean)[idx] + ), f"{dt} was touched by the blind" + # the blindable rows did move (the blind actually blinded) + for name, grid in ( + ("xi_reporting", "reporting"), + ("xi_integration", "integration"), + ): + idx = bd._xi_indices(parts[name], grid) + assert not np.allclose( + np.array(blinded[name].mean)[idx], np.array(parts[name].mean)[idx], atol=0 + ) + # ρ/τ: refused by blind_sacc (test_blind_refuses_non_blindable_part) and + # exempt in assembly — pass through gather untouched + s = sio.gather( + [ + blinded["xi_reporting"], + parts["rho_tau"], + blinded["cl"], + blinded["xi_integration"], + ] + ) + idx = s.indices(sio.RHO_PLUS.format(k=0)) + assert np.array_equal( + np.array(s.mean)[idx], + np.array(parts["rho_tau"].mean)[ + parts["rho_tau"].indices(sio.RHO_PLUS.format(k=0)) + ], + ) + + +@pytest.mark.slow +def test_ac9_pure_eb_nan_parity_under_blind(): + """AC9: the pure-E/B NaN pattern born from the blinded parts is identical + to the true parts' — blinding never moves a NaN. + + The Schneider estimator returns NaN wherever a reporting point lacks + interior support against the edge-based integration bounds. Whatever + that pattern is on the true parts (the fixture puts the outermost + reporting point at the boundary, so it is non-empty here; on production + files it is empty), the blinded derivation must reproduce it bit-for-bit: + the blind is a pure shift of the estimator's inputs, not a change of + estimator support. The finite values move (the ξ± shifted); the NaN mask + does not.""" + pytest.importorskip("cosmo_numba") + seed = "ac9-seed" + reporting, integration = make_xi_part("reporting"), make_xi_part("integration") + _, _, modes_t = _derive_downstream(reporting, integration) + _, _, modes_b = _derive_downstream( + bd.blind_sacc(reporting, seed, log=_NOLOG), + bd.blind_sacc(integration, seed, log=_NOLOG), + ) + for key in modes_t: + t, b = modes_t[key], modes_b[key] + assert np.array_equal(np.isnan(t), np.isnan(b)), ( + f"blinding moved the pure-E/B NaN pattern for {key}" + ) + # the finite values did move (the blind actually shifted the ξ±) + t, b = modes_t["xip_E"], modes_b["xip_E"] + finite = np.isfinite(t) + assert finite.any() and not np.allclose(t[finite], b[finite], atol=0) + + +# --------------------------------------------------------------------------- # +# AC6 + AC8: end-to-end custody through the file surface (slow + cosmo_numba) +# --------------------------------------------------------------------------- # +@pytest.mark.slow +def test_ac6_ac8_end_to_end_init_parts_gather_unblind(tmp_path): + """AC8: blind-init → blind-part on each intermediate part → terminal + gather (hash assertion passes, stamp lands) → unblind restores each part + bit-for-bit, and the derived statistics re-derived from the unblinded + integration part reproduce the truth. AC6: no plaintext part or seed survives + on disk, no ``seed_smokescreen`` key; unblind fails closed on a tampered + commitment.""" + pytest.importorskip("cosmo_numba") + parts = make_parts(nbins=1) + En_true, Bn_true, modes_true = _derive_downstream( + parts["xi_reporting"], parts["xi_integration"] + ) + + blind_dir = tmp_path / "blind" + blind_dir.mkdir() + init = bd.blind_init(str(blind_dir), log=_NOLOG) + + part_files, out_paths = {}, {} + for name in ("xi_reporting", "xi_integration", "cl"): + path = tmp_path / f"{name}.fits" + sio.save(parts[name], str(path), type="mock") + out_paths[name] = bd.blind_part(str(path), str(blind_dir), log=_NOLOG) + part_files[name] = path + + # -- custody hygiene (AC6) --------------------------------------------- -- + for name, path in part_files.items(): + assert not path.exists(), f"plaintext part {name} was not deleted" + blinded = sio.load(out_paths[name]["blinded"]) + assert blinded.metadata["concealed"] is True + assert "seed_smokescreen" not in blinded.metadata + assert pathlib.Path(out_paths[name]["escrow"]).exists() + assert not np.array_equal(np.array(blinded.mean), np.array(parts[name].mean)) + with open(init["commitment"], encoding="utf-8") as f: + commitment = json.load(f) + assert set(commitment) == { + "label", + "seed_commitment", + "config_digest", + "draw_scheme", + } + blinded_parts = {n: sio.load(p["blinded"]) for n, p in out_paths.items()} + for b in blinded_parts.values(): + assert b.metadata["blind_commitment"] == commitment["seed_commitment"] + # no plaintext json anywhere beside the blind outputs + assert not list(tmp_path.rglob("*escrow.json")) + assert not (blind_dir / "blind_seed.json").exists() + + # -- terminal assembly: hash assertion + stamp (AC8) -------------------- -- + assembled = sio.gather( + [ + blinded_parts["xi_reporting"], + blinded_parts["cl"], + parts["rho_tau"], + blinded_parts["xi_integration"], + ], + metadata={"catalogue_version": "vTEST", "type": "mock"}, + ) + assert assembled.metadata["blind_commitment"] == commitment["seed_commitment"] + # born-blinded derived statistics from the blinded parts differ from truth + En_b, Bn_b, _ = _derive_downstream( + blinded_parts["xi_reporting"], blinded_parts["xi_integration"] + ) + assert not np.allclose(En_b, En_true, atol=0) + + # -- fail-closed on tampered commitment (AC6) --------------------------- -- + tampered = dict(commitment, seed_commitment="0" * 64) + with open(init["commitment"], "w", encoding="utf-8") as f: + json.dump(tampered, f) + with pytest.raises(ValueError, match="committed seed commitment"): + bd.unblind_part( + out_paths["xi_integration"]["blinded"], + str(blind_dir), + str(tmp_path / "never.fits"), + log=_NOLOG, + ) + with open(init["commitment"], "w", encoding="utf-8") as f: + json.dump(commitment, f) + with pytest.raises(ValueError, match="config digest"): + bd.unblind_part( + out_paths["xi_integration"]["blinded"], + str(blind_dir), + str(tmp_path / "never.fits"), + config=bd.BlindingConfig(s8_half_width=0.01), + log=_NOLOG, + ) + + # -- bit-for-bit restoration per part (AC8) ----------------------------- -- + restored = {} + for name in ("xi_reporting", "xi_integration", "cl"): + out = tmp_path / f"{name}_restored.fits" + bd.unblind_part( + out_paths[name]["blinded"], str(blind_dir), str(out), log=_NOLOG + ) + restored[name] = sio.load(str(out)) + assert np.array_equal( + np.array(restored[name].mean), np.array(parts[name].mean) + ), f"{name} not restored bit-for-bit" + assert not restored[name].metadata.get("concealed", False) + assert "blind_commitment" not in restored[name].metadata + + # unblinding then re-deriving reproduces the true derived statistics + En_r, Bn_r, modes_r = _derive_downstream( + restored["xi_reporting"], restored["xi_integration"] + ) + assert np.array_equal(En_r, En_true) and np.array_equal(Bn_r, Bn_true) + for key in modes_true: + t, r = modes_true[key], modes_r[key] + both_nan = np.isnan(t) & np.isnan(r) + assert np.array_equal(t[~both_nan], r[~both_nan]), key + assert np.array_equal(np.isnan(t), np.isnan(r)), key + + +def test_ac8_dotted_versioned_part_names_escrow_and_restore(tmp_path): + """AC8 under the canonical catalogue-version naming (dotted stems). + + Production part files carry the versioned name ``v1.4.6.3_xi_reporting.fits`` + etc. ``smokescreen.encryption.encrypt_file`` names its outputs from + ``basename.split('.')[0]``, so both these parts would misfile onto + ``v1.encrpt``/``v1.key`` and the second would silently overwrite the + first's escrowed truth. Guard: the escrow lands at the exact + :func:`part_paths` name, two dot-prefix-sharing parts do not collide, and + each restores bit-for-bit.""" + parts = make_parts(nbins=1) + blind_dir = tmp_path / "blind" + blind_dir.mkdir() + bd.blind_init(str(blind_dir), log=_NOLOG) + + version = "v1.4.6.3" + out_paths, part_files = {}, {} + for name in ("xi_reporting", "xi_integration"): + path = tmp_path / f"{version}_{name}.fits" + sio.save(parts[name], str(path), type="mock") + out_paths[name] = bd.blind_part(str(path), str(blind_dir), log=_NOLOG) + part_files[name] = path + + # escrow bundles landed at the declared names (no split('.') truncation), + # and the two dot-prefix-sharing parts did not collide onto one bundle. + escrow_files = {n: p["escrow"] for n, p in out_paths.items()} + assert escrow_files["xi_reporting"] != escrow_files["xi_integration"] + for name, path in part_files.items(): + assert not path.exists(), f"plaintext part {name} was not deleted" + assert pathlib.Path(out_paths[name]["escrow"]).exists(), name + assert pathlib.Path(out_paths[name]["escrow_key"]).exists(), name + # the truncated-name collision target must not exist + assert not (tmp_path / "v1.encrpt").exists() + assert not (tmp_path / "v1.key").exists() + assert not list(tmp_path.rglob("*escrow.json")) + + # each part restores bit-for-bit via its own escrow (not subtraction-only) + for name in ("xi_reporting", "xi_integration"): + out = tmp_path / f"{version}_{name}_restored.fits" + bd.unblind_part( + out_paths[name]["blinded"], str(blind_dir), str(out), log=_NOLOG + ) + restored = sio.load(str(out)) + assert np.array_equal(np.array(restored.mean), np.array(parts[name].mean)), ( + f"{name} not restored bit-for-bit" + ) diff --git a/src/sp_validation/tests/test_blinding_wiring.py b/src/sp_validation/tests/test_blinding_wiring.py new file mode 100644 index 00000000..5bf5b0d7 --- /dev/null +++ b/src/sp_validation/tests/test_blinding_wiring.py @@ -0,0 +1,380 @@ +"""Tests for the Snakemake blind-at-birth wiring (issues #247/#252, PR #253). + +Two seams are covered here, both independent of a live cluster: + +1. **The path helpers in ``workflow/common.py``** must stay in lockstep with + ``sp_validation.blinding`` — common mirrors ``init_paths`` / ``part_paths`` by + hand (to keep the DAG build from importing the heavy blinding module), so a + drift between them would silently mis-wire ``blind_part``. These tests are the + guard. +2. **The data-run fail-closed assembly**: ``assemble_sacc`` must refuse an + unblinded ``type='data'`` part and succeed once every part is concealed under + one commitment — the terminal custody gate of #252. + +A candide-only test additionally asserts the blinding subgraph resolves in the +cosmo_val DAG dry-run. +""" + +import importlib.util +import json +import os +import subprocess +import sys +import types +from pathlib import Path + +import numpy as np +import pytest + +from sp_validation import blinding +from sp_validation import sacc_io as sio +from sp_validation.cosmo_val import sacc_writers as sw + + +def _repo_root(): + return next( + p for p in Path(__file__).resolve().parents if (p / "pyproject.toml").exists() + ) + + +def _load_module(rel_path, name): + """Import a workflow module/script by file path (off the package path).""" + path = _repo_root() / rel_path + spec = importlib.util.spec_from_file_location(name, path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +common = _load_module("workflow/common.py", "wf_common") +asm = _load_module("workflow/scripts/assemble_sacc.py", "assemble_sacc") + + +# --------------------------------------------------------------------------- # +# 1. common.py path helpers mirror sp_validation.blinding (drift guard) +# --------------------------------------------------------------------------- # +_STEMS = [ + "SP_v1.4.6.3_xi_reporting_minsep=1.0_maxsep=250.0_nbins=20_npatch=100", + "SP_v1.4.6.3_leak_corr_xi_integration", + "pseudo_cl_SP_v1.4.6.3_blind=A_powspace_nbins=32", + "pseudo_cl_SP_v1.4.6.3_leak_corr_blind=A_powspace_nbins=32", +] + + +@pytest.mark.parametrize("stem", _STEMS) +def test_blinded_path_mirrors_blinding_part_paths(stem): + part = f"/out/{stem}.sacc" + assert common.blinded_path(part) == blinding.part_paths(part)["blinded"] + + +def test_blind_state_paths_mirror_blinding_init_paths(): + version = "SP_v1.4.6.3_leak_corr" + common_paths = common.blind_state_paths(version) + ref = blinding.init_paths(common.blind_state_dir(version)) + assert common_paths == ref + + +@pytest.mark.parametrize( + "stem,expected", + [ + (_STEMS[0], "SP_v1.4.6.3"), + (_STEMS[1], "SP_v1.4.6.3_leak_corr"), + (_STEMS[2], "SP_v1.4.6.3"), + (_STEMS[3], "SP_v1.4.6.3_leak_corr"), + ], +) +def test_version_of_extracts_catalogue_version(stem, expected): + assert common.version_of(stem) == expected + + +def test_version_of_raises_without_version(): + with pytest.raises(ValueError, match="no catalogue version"): + common.version_of("cosebis_no_version_here") + + +def test_blindable_part_switches_on_run_type(monkeypatch): + part = "/out/SP_v1.4.6.3_xi_integration.sacc" + monkeypatch.setattr(common, "RUN_TYPE", "data") + assert common.blindable_part(part) == common.blinded_path(part) + monkeypatch.setattr(common, "RUN_TYPE", "mock") + assert common.blindable_part(part) == part + + +# --------------------------------------------------------------------------- # +# 2. Data-run fail-closed assembly (#252 terminal custody gate) +# --------------------------------------------------------------------------- # +META = {"catalogue_version": "vSYNTH", "npatch": 1} +# Two arbitrary-but-consistent hex stamps standing in for a real blind's +# seed commitment / config digest; the assembly only checks they agree across parts. +_COMMIT = "a" * 64 +_DIGEST = "b" * 64 + + +def _nz(): + return np.linspace(0.01, 2.0, 40), np.random.default_rng(0).uniform(0.1, 1.0, 40) + + +def _spd(n, seed): + a = np.random.default_rng(seed).normal(size=(n, n)) + return a @ a.T + n * np.eye(n) + + +def _data_parts(tmp_path, *, conceal, one_plaintext=False, run_type="data"): + """Write the five per-statistic parts, stamped ``type=run_type``. + + ``conceal`` stamps every part with the shared blind (concealed=True). With + ``one_plaintext`` the ξ± reporting part is left unconcealed — a blinded / + plaintext mix the assembly must refuse. ``run_type='mock'`` writes the parts + as a mock campaign's producers do, which is the only way an unconcealed + blindable part is allowed through the assembly. + """ + nz = {0: _nz()} + theta = np.geomspace(1.0, 100.0, 6) + + xi = sw.xi_to_sacc( + nz, META, theta, np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="reporting" + ) + xi.add_covariance(_spd(len(xi.mean), 1)) + co = sw.cosebis_to_sacc( + nz, + META, + { + "En": np.arange(1, 6) * 1e-6, + "Bn": np.arange(1, 6) * 1e-7, + "cov": _spd(10, 3), + }, + (1.0, 100.0), + ) + eb_arrays = {k: np.arange(6) * (i + 1) * 1e-6 for i, k in enumerate(sio.PURE_KEYS)} + eb = sw.pure_eb_to_sacc(nz, META, theta, eb_arrays, covariance=_spd(36, 4)) + rho = {"theta": theta} + tau = {"theta": theta} + rng = np.random.default_rng(5) + for k in sw.RHO_K: + for s in ("p", "m"): + rho[f"rho_{k}_{s}"] = rng.normal(size=6) * 1e-6 + rho[f"varrho_{k}_{s}"] = rng.uniform(1e-14, 1e-13, 6) + for k in sw.TAU_K: + for s in ("p", "m"): + tau[f"tau_{k}_{s}"] = rng.normal(size=6) * 1e-6 + tau[f"vartau_{k}_{s}"] = rng.uniform(1e-14, 1e-13, 6) + rt = sw.rho_tau_to_sacc(nz, META, rho, tau) + + parts = {"xi_reporting": xi, "cosebis": co, "pure_eb": eb, "rho_tau": rt} + paths = {} + for name, part in parts.items(): + if conceal and not (one_plaintext and name == "xi_reporting"): + blinding._stamp_provenance(part, _COMMIT, "A", _DIGEST) + p = tmp_path / f"{name}.sacc" + sio.save(part, str(p), type=run_type) + paths[name] = str(p) + return paths + + +def test_data_assemble_fails_closed_on_unblinded_part(tmp_path): + """A data run refuses to assemble an unconcealed real part (fail closed).""" + paths = _data_parts(tmp_path, conceal=False) + with pytest.raises(ValueError, match="refusing to load an unblinded"): + asm.assemble_sacc( + "vSYNTH", paths, str(tmp_path / "vSYNTH.sacc"), placeholder_var=1.0 + ) + + +def test_data_assemble_passes_on_blinded_parts(tmp_path): + """With every part concealed under one blind, the data-run assembly succeeds + and stamps the shared commitment on the terminal file.""" + paths = _data_parts(tmp_path, conceal=True) + out = tmp_path / "vSYNTH.sacc" + s = asm.assemble_sacc("vSYNTH", paths, str(out), placeholder_var=1.0) + assert s.metadata["concealed"] is True + assert s.metadata["blind_commitment"] == _COMMIT + assert s.metadata["blind_config_digest"] == _DIGEST + # Round-trips through the fail-closed load gate without an escape hatch. + assert sio.load(str(out)).metadata["concealed"] is True + + +def test_data_assemble_refuses_blinded_plaintext_mix(tmp_path): + """A concealed ξ± beside a plaintext one is a custody violation — refuse.""" + paths = _data_parts(tmp_path, conceal=True, one_plaintext=True) + # The plaintext ξ± reporting part fails the load gate first (data + not + # concealed), so the mix can never even reach assembly. + with pytest.raises(ValueError, match="refusing to load an unblinded"): + asm.assemble_sacc( + "vSYNTH", paths, str(tmp_path / "vSYNTH.sacc"), placeholder_var=1.0 + ) + + +def test_data_assemble_runs_behind_the_custody_guard(tmp_path, monkeypatch): + """The custody guard is reached *through* the production assembly. + + Every part here is concealed, so every part clears the fail-closed load + gate — the earlier tests all stop there. Only + ``blinding.assert_consistent_blind`` can catch what is wrong with this + assembly: the install's Smokescreen draws shifts under a different scheme + than the one that made the blind, so it could never unblind what it is + about to write. That ``assemble_sacc`` raises is the check that it runs + through :func:`sacc_io.gather` like every other assembly path, rather than + reimplementing the custody wrapper around its own assembler. + """ + paths = _data_parts(tmp_path, conceal=True) + monkeypatch.setattr(blinding, "draw_scheme", lambda: 99) + with pytest.raises(ValueError, match="DRAW_SCHEME"): + asm.assemble_sacc( + "vSYNTH", paths, str(tmp_path / "vSYNTH.sacc"), placeholder_var=1.0 + ) + + +def test_data_assemble_stamps_the_draw_scheme_on_the_terminal_file(tmp_path): + """The assembled file carries the blind's draw scheme, like its parts.""" + paths = _data_parts(tmp_path, conceal=True) + out = tmp_path / "vSYNTH.sacc" + asm.assemble_sacc("vSYNTH", paths, str(out), placeholder_var=1.0) + assert sio.load(str(out)).metadata["blind_draw_scheme"] == blinding.draw_scheme() + + +def test_mock_assemble_succeeds_without_a_blind(tmp_path): + """A mock campaign assembles its plaintext parts — the gate's mock branch is + reachable from real producers, not only from test fixtures. + + Every part writer stamps the campaign's run type (CosmologyValidation's + ``run_type``, ``run_2pcf``'s ``run_type=``, ``run_2pcf_highres``'s + ``--run-type``), so a ``mock`` campaign's parts declare themselves mocks and + ``assert_consistent_blind`` lets them through unconcealed. The same parts + stamped ``type='data'`` fail closed — that is + ``test_data_assemble_fails_closed_on_unblinded_part``. + """ + paths = _data_parts(tmp_path, conceal=False, run_type="mock") + out = tmp_path / "vSYNTH.sacc" + s = asm.assemble_sacc("vSYNTH", paths, str(out), placeholder_var=1.0) + assert "concealed" not in s.metadata + # No escape hatch: a mock part is not gated by the fail-closed loader. + assert sio.load(str(out)).metadata["type"] == "mock" + + +def test_rho_tau_part_is_stamped_concealed_from_the_commitment(tmp_path): + """ρ/τ carries no cosmological vector, but a data run's assembly still opens + it through the fail-closed load gate — so its writer must stamp it. + + This drives the real writer (``PSFSystematicsMixin.rho_tau_to_sacc_part``) + with the commitment the ``rho_tau_stats`` rule binds on a data run, and + checks the emitted part loads without the escape hatch. Without the stamp a + data run's ``assemble_sacc`` dies on the ρ/τ part before custody is ever + checked. + """ + from sp_validation.cosmo_val.psf_systematics import PSFSystematicsMixin + + blind_dir = tmp_path / "blind" + blind_dir.mkdir() + commitment = blinding.blind_init(str(blind_dir), log=lambda *_: None)["commitment"] + theta = np.geomspace(1.0, 100.0, 6) + rng = np.random.default_rng(11) + rho = {"theta": theta} + tau = {"theta": theta} + for k in sw.RHO_K: + for sign in ("p", "m"): + rho[f"rho_{k}_{sign}"] = rng.normal(size=6) * 1e-6 + rho[f"varrho_{k}_{sign}"] = rng.uniform(1e-14, 1e-13, 6) + for k in sw.TAU_K: + for sign in ("p", "m"): + tau[f"tau_{k}_{sign}"] = rng.normal(size=6) * 1e-6 + tau[f"vartau_{k}_{sign}"] = rng.uniform(1e-14, 1e-13, 6) + + class _Writer(PSFSystematicsMixin): + """The writer's collaborators, stubbed — the method under test is real.""" + + run_type = "data" + + def sacc_nz(self, version): + return {0: _nz()} + + def sacc_metadata(self, version): + return dict(META) + + def print_magenta(self, *args, **kwargs): + pass + + out_dir = tmp_path / "rho_tau_stats" + out_dir.mkdir() + _Writer().rho_tau_to_sacc_part( + "vSYNTH", + str(out_dir), + "vSYNTH", + types.SimpleNamespace(rho_stats=rho), + types.SimpleNamespace(tau_stats=tau), + commitment_path=commitment, + ) + + written = sio.load(str(out_dir / "rho_tau_vSYNTH.sacc")) + assert written.metadata["concealed"] is True + assert written.metadata["blind_draw_scheme"] == blinding.draw_scheme() + with open(commitment, encoding="utf-8") as f: + committed = json.load(f) + assert written.metadata["blind_commitment"] == committed["seed_commitment"] + assert written.metadata["blind_config_digest"] == committed["config_digest"] + + +def test_assert_consistent_blind_rejects_divergent_commitments(tmp_path): + """Two ξ± parts blinded under different commitments must never combine.""" + nz = {0: _nz()} + theta = np.geomspace(1.0, 100.0, 6) + a = sw.xi_to_sacc( + nz, META, theta, np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="reporting" + ) + b = sw.xi_to_sacc( + nz, META, theta, np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="integration" + ) + blinding._stamp_provenance(a, _COMMIT, "A", _DIGEST) + blinding._stamp_provenance(b, "c" * 64, "A", _DIGEST) + with pytest.raises(ValueError, match="different blind commitments"): + blinding.assert_consistent_blind([a, b]) + + +# --------------------------------------------------------------------------- # +# 3. The blinding subgraph resolves in the cosmo_val DAG (candide-only) +# --------------------------------------------------------------------------- # +requires_candide_data = pytest.mark.skipif( + not Path("/n17data/cdaley/unions").exists(), + reason="candide-local workflow config/data (/n17data) absent — off-cluster", +) + + +@requires_candide_data +def test_blinding_subgraph_in_cosmo_val_dry_run(): + """A data-run cosmo_val assemble pulls blind_init + blind_part, and binds the + ξ± / pseudo-Cℓ parts to their *_blinded siblings.""" + env = os.environ | {"PYTHONNOUSERSITE": "1", "PYTHONUNBUFFERED": "1"} + env.pop("SNAKEMAKE_PROFILE", None) + result = subprocess.run( + [ + sys.executable, + "-m", + "snakemake", + "assemble_sacc_all", + "--dry-run", + "--cores", + "1", + "--configfile", + "config/config.yaml", + ], + cwd=_repo_root() / "papers/cosmo_val", + env=env, + text=True, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + timeout=180, + check=False, + ) + assert result.returncode == 0, result.stdout + out = result.stdout + assert "rule blind_init:" in out, out + assert "rule blind_part:" in out, out + # assemble consumes the blinded ξ± reporting and pseudo-Cℓ. The integration + # ξ± is not gathered into the terminal file (per the #247 ruling), but the + # COSEBIs / pure-E/B consumers now bind the blinded integration part (via + # blindable_part) to re-derive their concealed E-modes, so blind_part enters + # the subgraph for it too. + assert ( + "_xi_reporting_minsep=1.0_maxsep=250.0_nbins=20_npatch=100_blinded.sacc" in out + ) + assert "_xi_integration_blinded.sacc" in out + assert "_blind=A_powspace_nbins=32_blinded.sacc" in out diff --git a/src/sp_validation/tests/test_bmodes_workflow_dry_run.py b/src/sp_validation/tests/test_bmodes_workflow_dry_run.py index cc6a293c..17b18b29 100644 --- a/src/sp_validation/tests/test_bmodes_workflow_dry_run.py +++ b/src/sp_validation/tests/test_bmodes_workflow_dry_run.py @@ -96,3 +96,32 @@ def test_cosmo_val_workflow_assemble_dry_runs(): assert f"pseudo_cl_cov_{version}_blind=A_powspace_nbins=32.fits" in out, out for part in ("_xi_reporting_", "_cosebis.sacc", "_pure_eb.sacc", "rho_tau_"): assert part in out, f"missing {part} part in assemble DAG:\n{out}" + + +@requires_candide_data +def test_cosmo_val_inference_prep_dry_runs(): + """The revived (PR-7) inference_prep DAG resolves end to end from the SACC. + + inference_fiducial must pull inference_prep, which consumes the assembled + {version}.sacc and emits the converter 2pt-FITS plus BOTH generated pipeline + inis (2pt_like and the native sacc_like). The old cosmosis_fitting.py real- + data assembly is retired from this path; the glass-mock rules keep it. + """ + version = "SP_v1.4.6.3_leak_corr" + result = _dry_run(_repo_root() / "papers/cosmo_val", ["inference_fiducial"]) + assert result.returncode == 0, result.stdout + out = result.stdout + assert "rule inference_prep:" in out, out + assert "rule inference_fiducial:" in out, out + # inference_prep consumes the assembled analysis SACC (not per-sign xi FITS). + assert f"{version}.sacc" in out, out + # ...and both ini TEMPLATES, bound as inputs so a template edit regenerates the + # configs (Finding 2 — a template as params gave no DAG edge). + assert "cosmosis_pipeline_A_ia.ini" in out, out + assert "cosmosis_pipeline_A_ia_sacc.ini" in out, out + # It emits the converter FITS + both engine inis. + assert f"cosmosis_{version}.fits" in out, out + assert f"cosmosis_pipeline_{version}_A_ia.ini" in out, out + assert f"cosmosis_pipeline_{version}_A_ia_sacc.ini" in out, out + # The retired real-data assembly script must not appear in this DAG's prep. + assert "cosmosis_fitting.py --cosmosis-root" not in out, out 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..3f286474 --- /dev/null +++ b/src/sp_validation/tests/test_camb_ccl_crosscheck.py @@ -0,0 +1,175 @@ +"""CAMB↔CCL theory cross-check (blinding PRD AC10–14). + +The blinding shift is a difference of CCL theory vectors; downstream +inference runs CAMB (CosmoSIS). The shift only means what it is intended to +mean if CCL and CAMB predict the same ξ± at a fixed cosmology on our θ grid. +This module asserts that agreement between the two independent ξ± paths in +:mod:`sp_validation.blinding_theory`: + +- **Path A** (:func:`~sp_validation.blinding_theory.xi_ccl`): CCL-native — CCL's + Boltzmann-CAMB HMCode2020 P(k) route, projected by CCL Limber + FFTLog. +- **Path B** (:func:`~sp_validation.blinding_theory.xi_camb`): an independent + pycamb run produces the HMCode2020 ``P(k, z)`` (σ8-matched ``A_s``), + wrapped in a ``ccl.Pk2D`` and projected through the same CCL machinery. + +Because both paths route their nonlinear P(k) through CAMB's HMCode2020 and +both project through CCL, a common Limber+FFTLog bug cancels: this test +validates the **P(k) recipe** and the **σ8/A_s amplitude convention**, not +the projection. The one convention subtlety it settles: the fiducial fixes +σ8 for CCL but A_s for CAMB; a nominal ``A_s = 2.1e-9`` leaves CAMB's σ8 +≈3% off target — enough to blow a ξ± comparison to ~9–10%. +""" + +import pathlib +import re + +import numpy as np +import pytest + +from sp_validation import blinding_theory as cm + +# Tolerances (AC11/AC12). Observed floor on this fixture: see the printed +# numbers in the slow tests — the tolerances sit above the floor with +# headroom; version bumps move the floor and that is not a regression. +XIP_RTOL = 0.005 # 0.5 % +XIM_RTOL = 0.010 # 1.0 % +# ξ− crosses zero on this grid: the relative assertion applies only where +# |ξ−| exceeds an absolute floor set from the fixture's peak |ξ−|. +XIM_FLOOR_FRAC = 0.05 + + +# --------------------------------------------------------------------------- # +# Deterministic fixture: one Gaussian source bin, 12-point θ grid +# --------------------------------------------------------------------------- # +def _gauss_nz(n=400): + z = np.linspace(0.01, 3.0, n) + nz = np.exp(-0.5 * ((z - 0.7) / 0.2) ** 2) + return z, nz / np.trapezoid(nz, z) + + +THETA_ARCMIN = np.geomspace(5.0, 250.0, 12) + + +def _both_paths(config, **camb_kwargs): + z, nz = _gauss_nz() + xip_a, xim_a = cm.xi_ccl( + config.ccl_params(), config, (z, nz), (z, nz), THETA_ARCMIN + ) + xip_b, xim_b, As = cm.xi_camb(config, (z, nz), THETA_ARCMIN, **camb_kwargs) + return (xip_a, xim_a), (xip_b, xim_b), As + + +def _assert_xi_agreement(a, b, label): + (xip_a, xim_a), (xip_b, xim_b) = a, b + assert np.all(xip_a > 0) and np.all(xip_b > 0) # sensible cosmic shear + rel_p = np.abs(xip_b - xip_a) / np.abs(xip_a) + assert rel_p.max() < XIP_RTOL, ( + f"{label}: ξ+ max rel diff {rel_p.max():.3%} ≥ {XIP_RTOL:.1%}" + ) + floor = XIM_FLOOR_FRAC * np.max(np.abs(xim_a)) + above = np.abs(xim_a) > floor + rel_m = np.abs(xim_b - xim_a)[above] / np.abs(xim_a)[above] + assert rel_m.max() < XIM_RTOL, ( + f"{label}: ξ− max rel diff {rel_m.max():.3%} ≥ {XIM_RTOL:.1%} (on |ξ−| > floor)" + ) + # near the zero crossing: absolute agreement at the floor scale + abs_m = np.abs(xim_b - xim_a)[~above] + if len(abs_m): + assert abs_m.max() < XIM_RTOL * floor, ( + f"{label}: ξ− absolute diff {abs_m.max():.3e} near zero crossing" + ) + print( + f"\n{label}: ξ+ max rel {rel_p.max():.3%}; " + f"ξ− max rel {rel_m.max():.3%} (above floor, " + f"{above.sum()}/{len(above)} points)" + ) + + +# --------------------------------------------------------------------------- # +# AC10: σ8/A_s reconciliation +# --------------------------------------------------------------------------- # +@pytest.mark.slow +def test_ac10_sigma8_As_reconciliation(): + """(a) nominal A_s leaves CAMB's σ8 >2% off target — the convention + offset is real; (b) the closed-form rescale lands on target to <1e-4.""" + cfg = cm.TheoryConfig() + target = cfg.sigma8() + + nominal = cm.camb_linear_sigma8(cfg, 2.1e-9) + offset = abs(nominal / target - 1) + print(f"\nAC10 nominal-A_s σ8 offset: {offset:.4f}") + assert offset > 0.02 + + As = cm.camb_As_for_sigma8(cfg, target) + matched = cm.camb_linear_sigma8(cfg, As) + print(f"AC10 σ8-matched residual: {abs(matched - target):.2e} (A_s={As:.4e})") + assert abs(matched - target) < 1e-4 + + +# --------------------------------------------------------------------------- # +# AC11 + AC12: ξ± agreement at and off the fiducial +# --------------------------------------------------------------------------- # +@pytest.mark.slow +def test_ac11_xi_agreement_at_fiducial(): + cfg = cm.TheoryConfig() + a, b, _ = _both_paths(cfg) + _assert_xi_agreement(a, b, "AC11 fiducial") + + +@pytest.mark.slow +def test_ac12_xi_agreement_off_fiducial(): + """A representative in-envelope offset — the *shift* (a difference of two + theory vectors) must not inherit a stack-disagreement bias.""" + cfg = cm.TheoryConfig.from_overrides({"S8": 0.80 + 0.075, "Omega_m": 0.30 - 0.05}) + a, b, _ = _both_paths(cfg) + _assert_xi_agreement(a, b, "AC12 off-fiducial") + + +# --------------------------------------------------------------------------- # +# AC13: halofit token pinned to the inference config (fast) +# --------------------------------------------------------------------------- # +def test_ac13_halofit_token_matches_inference_config(): + """The blinding fiducial's CCL halofit token equals the CosmoSIS + inference config's ``halofit_version`` — asserted against the config + file itself. All three blinding backends share one recipe by + construction and would agree with each other while jointly diverging + from the inference stack, so this cannot be caught by the cross-backend + test and is asserted independently here.""" + ini = ( + pathlib.Path(__file__).resolve().parents[3] + / "cosmo_inference" + / "cosmosis_config" + / "templates" + / "cosmosis_pipeline_A_ia_cell.ini" + ) + match = re.search(r"^halofit_version\s*=\s*(\S+)", ini.read_text(), re.MULTILINE) + assert match, f"no halofit_version in {ini}" + inference_token = match.group(1) + cfg = cm.TheoryConfig() + assert cfg.ccl_halofit_version == inference_token + # the two stack tokens denote ONE recipe; a divergence is a config bug + assert cfg.camb_halofit_version == cfg.ccl_halofit_version + # #280: the shipped Boltzmann backend is CAMB-through-CCL, matching the + # CosmoSIS+CAMB inference stack — one power-spectrum path. The cross-check + # tests above (AC10–12, 14) all run at this default configuration. + assert cfg.transfer_function == "boltzmann_camb" + + +# --------------------------------------------------------------------------- # +# AC14: fast smoke — broken wiring caught in the fast suite +# --------------------------------------------------------------------------- # +def test_ac14_crosscheck_smoke(): + """Both paths run at coarse resolution: finite, positive, + few-percent-agreeing ξ+, and a σ8-matched A_s in a sane range.""" + cfg = cm.TheoryConfig() + z, nz = _gauss_nz(n=150) + theta = np.geomspace(10.0, 100.0, 4) + xip_a, _ = cm.xi_ccl(cfg.ccl_params(), cfg, (z, nz), (z, nz), theta) + xip_b, _, As = cm.xi_camb( + cfg, (z, nz), theta, n_ell=120, ell_max=30000, kmax=10.0, n_k=200 + ) + assert np.all(np.isfinite(xip_a)) and np.all(np.isfinite(xip_b)) + assert np.all(xip_a > 0) and np.all(xip_b > 0) + assert 1e-9 < As < 3e-9 + rel = np.abs(xip_b - xip_a) / np.abs(xip_a) + assert rel.max() < 0.05, f"smoke ξ+ rel diff {rel.max():.3%} unexpectedly large" diff --git a/src/sp_validation/tests/test_generate_inference_config.py b/src/sp_validation/tests/test_generate_inference_config.py new file mode 100644 index 00000000..4fbbb98d --- /dev/null +++ b/src/sp_validation/tests/test_generate_inference_config.py @@ -0,0 +1,138 @@ +"""Tests for the CosmoSIS inference-config generator. + +The generator (:mod:`workflow.scripts.generate_inference_config`) fills a +pipeline ini template's ``[DEFAULT]`` section with concrete paths so CosmoSIS's +ConfigParser resolves the template's ``%(KEY)s`` placeholders. These tests need +no cosmosis: they check that the substituted DEFAULT keys land, that the module +file paths the templates reference exist on disk, and that no ``%(...)s`` +placeholder is left unresolved after filling. +""" + +import configparser +import importlib.util +from pathlib import Path + +import pytest + +_REPO = Path(__file__).resolve().parents[3] +_SCRIPT = _REPO / "workflow" / "scripts" / "generate_inference_config.py" +_CONFIG_DIR = _REPO / "cosmo_inference" / "cosmosis_config" +_SACC_TEMPLATE = _CONFIG_DIR / "cosmosis_pipeline_A_ia_sacc.ini" +_FITS_TEMPLATE = _CONFIG_DIR / "cosmosis_pipeline_A_ia.ini" + + +def _load_generator(): + spec = importlib.util.spec_from_file_location("gen_inference_cfg", _SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +gen = _load_generator() + + +def _read_interpolated(ini_path): + """Parse the generated ini with interpolation ON (the way CosmoSIS reads it). + + CosmoSIS uses ``%(KEY)s`` BasicInterpolation with case-*preserved* keys, so + set ``optionxform = str`` (stdlib configparser lowercases keys by default, + which would break the uppercase ``%(FITS_FILE)s`` / ``%(SACC_FILE)s`` + lookups). The default BasicInterpolation then raises if a referenced key is + missing — exactly the failure we want to catch. + """ + parser = configparser.ConfigParser() + parser.optionxform = str + parser.read(ini_path) + return parser + + +def test_sacc_template_defaults_filled(tmp_path): + """The generated sacc ini carries the substituted DEFAULT keys and resolves.""" + out = tmp_path / "gen_sacc.ini" + gen.generate_inference_config( + _SACC_TEMPLATE, + out, + gen._substitutions( + scratch="/scratch/run", + cosmosis_dir="/csl", + sacc_file="/data/v1.sacc", + ), + ) + text = out.read_text() + assert "SCRATCH = /scratch/run" in text + assert "SACC_FILE = /data/v1.sacc" in text + assert "COSMOSIS_DIR = /csl" in text + assert "SP_VALIDATION_MODULES = " in text + # FITS_FILE is None for the sacc path — must be dropped, not written as "None". + assert "FITS_FILE" not in text + + parser = _read_interpolated(out) + # The sacc_like data_file interpolates SACC_FILE; load_nz_sacc too. + assert parser["sacc_like"]["data_file"] == "/data/v1.sacc" + assert parser["load_nz_sacc"]["nz_file"] == "/data/v1.sacc" + assert parser["sacc_like"]["csl_dir"] == "/csl" + + +def test_fits_template_defaults_filled(tmp_path): + """The generated 2pt_like ini carries the substituted DEFAULT keys and resolves.""" + out = tmp_path / "gen_fits.ini" + gen.generate_inference_config( + _FITS_TEMPLATE, + out, + gen._substitutions( + scratch="/scratch/run", + cosmosis_dir="/csl", + fits_file="/data/v1.fits", + ), + ) + text = out.read_text() + assert "SCRATCH = /scratch/run" in text + assert "FITS_FILE = /data/v1.fits" in text + assert "SACC_FILE" not in text + + parser = _read_interpolated(out) + assert parser["2pt_like"]["data_file"] == "/data/v1.fits" + assert parser["load_nz_fits"]["nz_file"] == "/data/v1.fits" + + +def test_sacc_template_module_file_exists(): + """The sacc_like module file the generated ini points at exists on disk. + + ``SP_VALIDATION_MODULES`` resolves to the installed package dir; + ``sacc_like_unions.py`` must live there (it is the shim CosmoSIS loads). + """ + modules = Path(gen._sp_validation_modules()) + assert (modules / "sacc_like_unions.py").is_file() + + +def test_all_placeholders_resolve(tmp_path): + """No ``%(...)s`` placeholder survives interpolation in either template. + + A missing DEFAULT key would make ConfigParser raise on access; iterate every + option in every section to force resolution of all placeholders. + """ + for template, subs in ( + ( + _SACC_TEMPLATE, + gen._substitutions(scratch="/s", cosmosis_dir="/csl", sacc_file="/d.sacc"), + ), + ( + _FITS_TEMPLATE, + gen._substitutions(scratch="/s", cosmosis_dir="/csl", fits_file="/d.fits"), + ), + ): + out = tmp_path / (template.stem + ".gen.ini") + gen.generate_inference_config(template, out, subs) + parser = _read_interpolated(out) + for section in parser.sections(): + for key in parser[section]: + value = parser[section][key] # raises if a placeholder is unresolved + assert "%(" not in value, f"[{section}] {key} = {value}" + + +def test_no_default_section_raises(tmp_path): + """A template with no [DEFAULT] section is a loud error.""" + bad = tmp_path / "bad.ini" + bad.write_text("[pipeline]\nmodules = a b c\n") + with pytest.raises(ValueError, match="DEFAULT"): + gen.generate_inference_config(bad, tmp_path / "out.ini", {"SCRATCH": "/s"}) diff --git a/src/sp_validation/tests/test_pseudo_cl.py b/src/sp_validation/tests/test_pseudo_cl.py index 8c095a71..dfb4d3e6 100644 --- a/src/sp_validation/tests/test_pseudo_cl.py +++ b/src/sp_validation/tests/test_pseudo_cl.py @@ -169,6 +169,7 @@ def cv(tmp_path): power=0.5, n_ell_bins=N_ELL_BINS, pol_factor=True, + run_type="mock", ) cv._test_version = version return cv diff --git a/src/sp_validation/tests/test_sacc_like.py b/src/sp_validation/tests/test_sacc_like.py new file mode 100644 index 00000000..f8f87857 --- /dev/null +++ b/src/sp_validation/tests/test_sacc_like.py @@ -0,0 +1,566 @@ +"""Equality tests: the sp_validation SACC-likelihood shim vs CosmoSIS ``2pt_like``. + +PR 7 adopts CosmoSIS's native ``SaccClLikelihood`` for the ξ± inference path, +through the shim :mod:`sp_validation.sacc_like_unions` (which fixes the upstream +arcmin→rad gap and adds an ordering guard). The contract is *in-process module +equality*: run the shimmed ``sacc_like`` on the analysis SACC and CosmoSIS's +``2pt_like`` on the PR-3 converter's 2pt-FITS against an identical synthetic +theory DataBlock, and require the same χ², log-likelihood, theory vector and +post-cut point count. The prototype (``sacc-like-probe/probe_equality.py``) +observed *exact* equality (Δχ²=0, Δtheory=0), so the equality tests assert +``array_equal`` / rtol=1e-12. + +Environment +----------- +Both engines are pure-python CosmoSIS module files (no CAMB, no compiled CSL +modules), so they run in the shared venv with ``cosmosis`` installed. They need a +checkout of the CosmoSIS Standard Library, located via the ``CSL_DIR`` env var; +the module is skipped when cosmosis is absent (CI image) or ``CSL_DIR`` is unset +or missing. Recipe:: + + git clone --depth 1 https://github.com/joezuntz/cosmosis-standard-library CSL + CSL_DIR=/path/to/CSL pytest ... test_sacc_like.py + +The pyproject ``cosmosis`` extra pins ``cosmosis>=3.25`` for the engine itself +(inert in CI, whose Dockerfile installs only ``[test,glass,blinding]``). +""" + +import importlib.util +import os +from pathlib import Path + +import numpy as np +import pytest + +cosmosis = pytest.importorskip("cosmosis") + +# The upstream CSL checkout carrying likelihood/sacc + likelihood/2pt. Skip the +# whole module (not error) when it is not configured, mirroring the cosmo_numba +# env-gated precedent. +_CSL_DIR = os.environ.get("CSL_DIR") +if not _CSL_DIR or not Path(_CSL_DIR, "likelihood", "sacc").is_dir(): + pytest.skip( + "CSL_DIR unset or has no likelihood/sacc — set CSL_DIR to a CosmoSIS " + "Standard Library checkout to run the sacc_like equality tests", + allow_module_level=True, + ) + +from cosmosis.datablock import DataBlock, option_section # noqa: E402 + +from sp_validation import sacc_io # noqa: E402 + +# Reuse the angular-bin count from the converter tests so the two suites' single- +# bin shapes stay in lockstep. (The χ²-dynamics builders here need a covariance +# commensurate with the data, which the converter's byte-compare _sacc is not.) +from sp_validation.tests.test_sacc_io_twopoint import N_ANG # noqa: E402 + +CSL = Path(_CSL_DIR) +ARCMIN_TO_RAD = np.pi / (180.0 * 60.0) + +# The like_name both engines are configured with, so both write identical block +# keys (_CHI2, _LIKE, _theory) — the parity the design mandates. +LIKE_NAME = "2pt_like" + + +# --------------------------------------------------------------------------- +# Shared engine harness +# --------------------------------------------------------------------------- +def _load_module_file(path, name): + """Import a CosmoSIS module file (``setup``/``execute``/``cleanup``) by path. + + ``likelihood/sacc`` and ``likelihood/2pt`` are added to ``sys.path`` so the + upstream modules resolve their siblings (``sacc_likelihoods``, ``spec_tools``, + ``twopoint_cosmosis``, …). Loading ``2pt_like.py`` runs its + ``build_module()`` at import (module-level ``setup, execute, cleanup``); + ``sacc_like_unions.py`` defines those functions directly. + """ + import sys + + for sub in ("likelihood/sacc", "likelihood/2pt"): + p = str(CSL / sub) + if p not in sys.path: + sys.path.insert(0, p) + spec = importlib.util.spec_from_file_location(name, path) + mod = importlib.util.module_from_spec(spec) + spec.loader.exec_module(mod) + return mod + + +def _theory_block(): + """A DataBlock carrying the synthetic ξ± theory both engines interpolate. + + Mirrors the probe: the theory θ grid is in *radians* (CosmoSIS convention), + and the ξ+/ξ− predictions are smooth power laws sampled on it. Both engines + build a spline over this grid and evaluate it at each data point's angular + tag — so as long as both see the same block, the interpolated theory (and + hence χ²) must agree. + """ + theta_arcmin = np.geomspace(0.3, 400.0, 300) + + def t_xip(th): + return 2e-4 * (th / 10.0) ** -0.8 + + def t_xim(th): + return 1e-4 * (th / 10.0) ** -0.5 + + b = DataBlock() + for section, f in (("shear_xi_plus", t_xip), ("shear_xi_minus", t_xim)): + b[section, "theta"] = theta_arcmin * ARCMIN_TO_RAD + b[section, "bin_1_1"] = f(theta_arcmin) + b[section, "is_auto"] = True + b[section, "nbin_a"] = 1 + b[section, "nbin_b"] = 1 + b[section, "sample_a"] = "nz_source" + b[section, "sample_b"] = "nz_source" + b[section, "sep_name"] = "theta" + b[section, "save_name"] = "" + return b + + +def _run(mod, options): + """Run a CosmoSIS likelihood module file standalone; return (like, chi2, theory, n). + + Builds the options DataBlock (module options live under ``option_section``), + calls ``setup`` → ``execute`` on a fresh theory block, and reads the standard + Gaussian-likelihood outputs back out under the shared ``LIKE_NAME`` keys. + """ + opt = DataBlock() + for key, value in options.items(): + opt[option_section, key] = value + config = mod.setup(opt) + block = _theory_block() + mod.execute(block, config) + like = block["likelihoods", f"{LIKE_NAME}_LIKE"] + chi2 = block["data_vector", f"{LIKE_NAME}_CHI2"] + theory = block["data_vector", f"{LIKE_NAME}_theory"] + return like, chi2, np.asarray(theory), len(theory) + + +# --------------------------------------------------------------------------- +# The two engine module files + their options +# --------------------------------------------------------------------------- +_SHIM_PATH = Path(__file__).resolve().parents[1] / "sacc_like_unions.py" + + +@pytest.fixture(scope="module") +def m_shim(): + return _load_module_file(_SHIM_PATH, "sacc_like_unions_test") + + +@pytest.fixture(scope="module") +def m_2pt(): + return _load_module_file(CSL / "likelihood/2pt/2pt_like.py", "twopt_like_test") + + +@pytest.fixture(scope="module") +def m_raw_sacc(): + return _load_module_file(CSL / "likelihood/sacc/sacc_like.py", "raw_sacc_like_test") + + +def _shim_opts(sacc_path, **extra): + return { + "csl_dir": str(CSL), + "data_file": sacc_path, + "data_sets": "galaxy_shear_xi_plus galaxy_shear_xi_minus", + "like_name": LIKE_NAME, + **extra, + } + + +def _twopt_opts(fits_path, **extra): + return { + "data_file": fits_path, + "data_sets": "XI_PLUS XI_MINUS", + "covmat_name": "COVMAT", + "like_name": LIKE_NAME, + "gaussian_covariance": False, + "cut_zeros": False, + **extra, + } + + +def _realistic_sacc(seed=0, *, xip=None): + """A single-bin ξ± SACC with data + covariance sized like the real product. + + The converter-test ``_sacc`` builder uses a covariance ~14 orders of + magnitude larger than the ξ values (fine for byte-comparing the converter, + where covariance *content* is irrelevant), which makes every χ² collapse to + numerical zero — no teeth for the unit-gap tripwire. This builder instead + lays down realistic ξ± power laws and a covariance ~ ``(0.1·|ξ|)²`` (the + probe's recipe), so χ² is O(1)-scale and the raw-vs-shim gap is visible. + + ``xip`` overrides the ξ+ values (perturbation teeth). + """ + rng = np.random.default_rng(seed) + n = N_ANG + theta = np.geomspace(1.0, 250.0, n) # arcmin + z = np.linspace(0.01, 3.0, 200) + nz = z**2 * np.exp(-((z / 0.5) ** 1.5)) + + xip_vals = 2e-4 * (theta / 10.0) ** -0.8 * (1 + 0.05 * rng.standard_normal(n)) + xim_vals = 1e-4 * (theta / 10.0) ** -0.5 * (1 + 0.05 * rng.standard_normal(n)) + if xip is not None: + xip_vals = xip + + s = sacc_io.new_sacc({0: (z, nz)}, {"catalogue_version": "test"}) + sacc_io.add_xi(s, (0, 0), theta, xip_vals, xim_vals, grid="reporting") + + sig = 0.1 * np.abs(np.concatenate([xip_vals, xim_vals])) + a = rng.standard_normal((2 * n, 3 * n)) + cov = (a @ a.T / (3 * n)) * np.outer(sig, sig) * 0.3 + np.diag(sig**2) + s.add_covariance(cov) + return s, theta + + +def _write_pair(tmp_path, seed=0, name="probe", *, xip=None): + """Write a realistic analysis SACC and its PR-3 converter 2pt-FITS. + + Returns ``(sacc_path, fits_path)``. The SACC carries the arcmin ξ± tags the + shim converts; the FITS is the byte-compatible product ``2pt_like`` reads. + """ + s, _theta = _realistic_sacc(seed, xip=xip) + sacc_path = str(tmp_path / f"{name}.sacc") + sacc_io.save(s, sacc_path, type="mock") + fits_path = str(tmp_path / f"{name}_2pt.fits") + sacc_io.sacc_to_twopoint_fits(sacc_io.load(sacc_path), fits_path, n_bins=1) + return sacc_path, fits_path + + +# --------------------------------------------------------------------------- +# 1. Core equality: shim on SACC ≡ 2pt_like on converter FITS +# --------------------------------------------------------------------------- +def test_shimmed_equals_2pt_like_exact(tmp_path, m_shim, m_2pt): + """The shimmed sacc_like and 2pt_like agree exactly on the same data+theory. + + Same synthetic theory block, the SACC through the shim vs the converter FITS + through 2pt_like: χ², log-likelihood and the theory vector must match to + numerical precision (the probe observed exact equality). This is the PR's + central contract — the native path reproduces the validated converter path. + """ + sacc_path, fits_path = _write_pair(tmp_path, seed=0) + + like_s, chi2_s, theory_s, n_s = _run(m_shim, _shim_opts(sacc_path)) + like_t, chi2_t, theory_t, n_t = _run(m_2pt, _twopt_opts(fits_path)) + + assert n_s == n_t == 2 * N_ANG + np.testing.assert_allclose(chi2_s, chi2_t, rtol=1e-12) + np.testing.assert_allclose(like_s, like_t, rtol=1e-12) + np.testing.assert_allclose(theory_s, theory_t, rtol=1e-12) + + +# --------------------------------------------------------------------------- +# 2. Tripwire: raw upstream sacc_like is broken on arcmin tags +# --------------------------------------------------------------------------- +def test_upstream_unit_gap_tripwire(tmp_path, m_raw_sacc, m_2pt): + """RAW upstream sacc_like (no shim) gives a wildly wrong χ² on arcmin tags. + + Documents and guards the arcmin→rad gap the shim fixes: with raw ``theta`` + tags the theory spline is evaluated outside its (radian) grid and returns 0, + collapsing χ² to dᵀC⁻¹d. We assert the relative χ² difference against + 2pt_like exceeds 10 (the probe saw Δχ²≈3184). + + IF THIS TEST EVER FAILS: upstream ``sacc_like`` has fixed its units — the + shim's arcmin→rad conversion is now redundant and the shim can be retired. + """ + sacc_path, fits_path = _write_pair(tmp_path, seed=0) + + _like_r, chi2_raw, _theory_r, _n_r = _run(m_raw_sacc, _shim_opts(sacc_path)) + _like_t, chi2_t, _theory_t, _n_t = _run(m_2pt, _twopt_opts(fits_path)) + + rel = abs(chi2_raw - chi2_t) / abs(chi2_t) + assert rel > 10, ( + f"raw sacc_like χ²={chi2_raw:.6g} is within 10× of 2pt_like χ²={chi2_t:.6g} " + f"(rel diff {rel:.3g}) — the upstream unit gap appears fixed; retire the shim" + ) + + +# --------------------------------------------------------------------------- +# 3. Scale cuts equivalent through both engines +# --------------------------------------------------------------------------- +def test_scale_cuts_equivalent(tmp_path, m_shim, m_2pt): + """The same arcmin scale cuts give the same post-cut N and χ² on both engines. + + Cuts are expressed in the each engine's grammar but the SAME numeric arcmin + values (shim cuts run before the arcmin→rad conversion, so they take arcmin + just like 2pt_like's angle_range). ξ+ ∈ [10, 200], ξ− ∈ [20, 200] arcmin. + """ + sacc_path, fits_path = _write_pair(tmp_path, seed=1) + + shim_cuts = _shim_opts( + sacc_path, + **{ + "angle_range_galaxy_shear_xi_plus_source_0_source_0": np.array( + [10.0, 200.0] + ), + "angle_range_galaxy_shear_xi_minus_source_0_source_0": np.array( + [20.0, 200.0] + ), + }, + ) + twopt_cuts = _twopt_opts( + fits_path, + **{ + "angle_range_XI_PLUS_1_1": np.array([10.0, 200.0]), + "angle_range_XI_MINUS_1_1": np.array([20.0, 200.0]), + }, + ) + + _like_s, chi2_s, _theory_s, n_s = _run(m_shim, shim_cuts) + _like_t, chi2_t, _theory_t, n_t = _run(m_2pt, twopt_cuts) + + assert n_s == n_t + assert n_s < 2 * N_ANG # the cuts actually removed points + np.testing.assert_allclose(chi2_s, chi2_t, rtol=1e-12) + + +# --------------------------------------------------------------------------- +# 4. Perturbation moves both engines identically +# --------------------------------------------------------------------------- +def test_perturbation_moves_both_identically(tmp_path, m_shim, m_2pt): + """Perturbing one data value shifts both engines' χ² by the identical amount. + + Teeth: build the base pair and a pair whose first ξ+ value is bumped, and + require Δχ²(shim) == Δχ²(2pt_like). If either engine ignored the perturbed + point (e.g. a misaligned data vector), the deltas would diverge. + """ + # Base pair, then a pair whose first ξ+ value is bumped. Rebuild the base ξ+ + # from the same seed so only the one perturbed entry differs. + base_s, theta = _realistic_sacc(seed=2) + base_xip = np.array( + [p.value for p in base_s.data if p.data_type == sacc_io.XI_PLUS] + ) + sacc_b, fits_b = _write_pair(tmp_path, seed=2, name="base") + + pert_xip = base_xip.copy() + pert_xip[0] += 5e-5 + sacc_p, fits_p = _write_pair(tmp_path, seed=2, name="pert", xip=pert_xip) + + _l, chi2_shim_b, _t, _n = _run(m_shim, _shim_opts(sacc_b)) + _l, chi2_shim_p, _t, _n = _run(m_shim, _shim_opts(sacc_p)) + _l, chi2_2pt_b, _t, _n = _run(m_2pt, _twopt_opts(fits_b)) + _l, chi2_2pt_p, _t, _n = _run(m_2pt, _twopt_opts(fits_p)) + + d_shim = chi2_shim_p - chi2_shim_b + d_2pt = chi2_2pt_p - chi2_2pt_b + assert abs(d_shim) > 0 # the perturbation actually moved χ² + np.testing.assert_allclose(d_shim, d_2pt, rtol=1e-10) + + +# --------------------------------------------------------------------------- +# 5. Ordering guard raises on a pair-major tomographic file +# --------------------------------------------------------------------------- +def test_ordering_guard_raises_pair_major_tomographic(tmp_path, m_shim): + """A pair-major 2-bin SACC trips the shim's ordering guard at setup. + + Inserting ξ± per pair — (0,0), then (0,1), then (1,1) — lays the data vector + out pair-major ([ξ+;ξ−] per pair), while ``get_data_types()`` groups the + theory loop type-major (all ξ+ pairs, then all ξ− pairs). The two orders + disagree (verified empirically), so the guard must raise a ValueError + mentioning the ordering hazard rather than silently mis-comparing. + """ + theta = np.geomspace(1.0, 250.0, N_ANG) # arcmin + z = np.linspace(0.01, 3.0, 200) + nz = z**2 * np.exp(-((z / 0.5) ** 1.5)) + xip = np.ones(N_ANG) * 1e-4 + xim = np.ones(N_ANG) * 1e-4 + s = sacc_io.new_sacc({0: (z, nz), 1: (z, nz)}) + for pair in [(0, 0), (0, 1), (1, 1)]: + sacc_io.add_xi(s, pair, theta, xip, xim, grid="reporting") + s.add_covariance(np.eye(len(s.mean))) + sacc_path = str(tmp_path / "pair_major.sacc") + sacc_io.save(s, sacc_path, type="mock") + + with pytest.raises(ValueError, match="ordering"): + m_shim.setup(_as_option_block(_shim_opts(sacc_path))) + + +def _as_option_block(options): + """Build a raw options DataBlock (keys under ``option_section``) for setup().""" + opt = DataBlock() + for key, value in options.items(): + opt[option_section, key] = value + return opt + + +# --------------------------------------------------------------------------- +# 6. theta conversion scoped to real-category types (COSEBIs untouched) +# --------------------------------------------------------------------------- +def test_theta_conversion_scoped_to_real_types(tmp_path, m_shim): + """The shim converts ξ ``theta`` tags but leaves cosebi ``n`` tags untouched. + + Build a SACC carrying both ξ± (real) and COSEBIs (a non-real ``cosebis`` + category with integer ``n`` tags). After the shim's ``build_data``, the ξ + ``theta`` tags must be scaled arcmin→rad (so they equal the original arcmin + values times the conversion factor), and the cosebi ``n`` tags must be + numerically unchanged. + """ + # Build ξ± + COSEBIs, then attach the covariance last (sacc forbids adding + # points after add_covariance). + theta = np.geomspace(1.0, 250.0, N_ANG) # arcmin + z = np.linspace(0.01, 3.0, 200) + nz = z**2 * np.exp(-((z / 0.5) ** 1.5)) + n_modes = 5 + En = np.arange(1.0, n_modes + 1) + Bn = np.arange(1.0, n_modes + 1) * 0.1 + + s = sacc_io.new_sacc({0: (z, nz)}) + sacc_io.add_xi( + s, (0, 0), theta, np.ones(N_ANG) * 1e-4, np.ones(N_ANG) * 1e-4, grid="reporting" + ) + sacc_io.add_cosebis(s, (0, 0), En, (1.0, 250.0), Bn=Bn) + s.add_covariance(np.eye(len(s.mean))) + + sacc_path = str(tmp_path / "with_cosebis.sacc") + sacc_io.save(s, sacc_path, type="mock") + + # data_sets keeps the cosebis in (so we can check its tags survive); cosebi's + # section/category resolve from sacc_like's default_sections, so build_data + # needs no extra ini config. Only setup() runs (build_data); the theory loop + # (which would want a cosebi theory block) runs at execute, not here. + config = m_shim.setup( + _as_option_block( + _shim_opts( + sacc_path, + data_sets=( + "galaxy_shear_xi_plus galaxy_shear_xi_minus " + "galaxy_shear_cosebi_ee galaxy_shear_cosebi_bb" + ), + ) + ) + ) + + # The conversion lives on the radian copy (_sacc_data_rad); the original + # self.sacc_data stays arcmin (so save_theory writes arcmin tags — Finding 1). + rad_xi_thetas = [ + p.tags["theta"] + for p in config._sacc_data_rad.data + if p.data_type == sacc_io.XI_PLUS + ] + rad_cosebi_ns = [ + p.tags["n"] + for p in config._sacc_data_rad.data + if p.data_type == sacc_io.COSEBI_EE + ] + orig_xi_thetas = [ + p.tags["theta"] for p in config.sacc_data.data if p.data_type == sacc_io.XI_PLUS + ] + # rad copy: ξ theta scaled to radians (original arcmin × ARCMIN_TO_RAD). + np.testing.assert_allclose( + np.sort(rad_xi_thetas), np.sort(theta * ARCMIN_TO_RAD), rtol=1e-12 + ) + # rad copy: cosebi n tags untouched (still integer modes 1..n_modes). + np.testing.assert_array_equal(np.sort(rad_cosebi_ns), np.arange(1, n_modes + 1)) + # original sacc_data: ξ theta still in arcmin (unmutated). + np.testing.assert_allclose(np.sort(orig_xi_thetas), np.sort(theta), rtol=1e-12) + + +# --------------------------------------------------------------------------- +# 6b. save_theory writes arcmin tags, and re-execute is stable (Finding 1) +# --------------------------------------------------------------------------- +def test_save_theory_writes_arcmin_and_reexecute_stable(tmp_path, m_shim): + """save_theory must write a SACC whose θ tags are still arcmin, not radians. + + Finding 1: because upstream save_theory copies self.sacc_data and overwrites + only point values, self.sacc_data must stay arcmin — otherwise the saved file + carries radian θ tags and any arcmin-assuming consumer (sacc_io.get_xi, or + re-ingesting it as a data_file, which would double-convert to ~8.5e-8) is + silently off by 3437×. Assert the saved θ tags match the input file's tags + exactly (arcmin) and that the saved values equal the theory vector. Also run + execute() twice and require identical χ² — a guard against any accidental + double-conversion creeping back in. + """ + s, theta = _realistic_sacc(seed=5) + sacc_path = str(tmp_path / "in.sacc") + sacc_io.save(s, sacc_path, type="mock") + input_theta = np.array( + [p.tags["theta"] for p in sacc_io.load(sacc_path).data if "theta" in p.tags] + ) + + save_path = str(tmp_path / "saved_theory.sacc") + opt = _as_option_block(_shim_opts(sacc_path, save_theory=save_path)) + config = m_shim.setup(opt) + + block1 = _theory_block() + m_shim.execute(block1, config) + chi2_1 = block1["data_vector", f"{LIKE_NAME}_CHI2"] + theory = np.asarray(block1["data_vector", f"{LIKE_NAME}_theory"]) + + saved = sacc_io.load(save_path) + saved_theta = np.array([p.tags["theta"] for p in saved.data if "theta" in p.tags]) + saved_values = np.array(saved.mean) + + # θ tags in the saved file are arcmin — identical to the input file's tags. + np.testing.assert_array_equal(saved_theta, input_theta) + # and are NOT the radian conversion (guards against the leak explicitly). + assert not np.allclose(saved_theta, input_theta * ARCMIN_TO_RAD) + # saved values are the theory vector (save_theory overwrites values in order). + np.testing.assert_allclose(saved_values, theory, rtol=1e-12) + + # A second execute() yields the identical χ² — no cumulative mutation. + block2 = _theory_block() + m_shim.execute(block2, config) + chi2_2 = block2["data_vector", f"{LIKE_NAME}_CHI2"] + np.testing.assert_allclose(chi2_2, chi2_1, rtol=1e-12) + + +# --------------------------------------------------------------------------- +# 7. Real-data equality (candide-gated) +# --------------------------------------------------------------------------- +_REALDATA = ( + Path("/automnt/n17data/cdaley/unions/code/sp_validation/cosmo_inference/data") + / "SP_v1.4.6_leak_corr_A_minsep=1.0_maxsep=250.0_nbins=20_npatch=1" + / "cosmosis_SP_v1.4.6_leak_corr_A_minsep=1.0_maxsep=250.0_nbins=20_npatch=1.fits" +) + + +@pytest.mark.skipif( + not _REALDATA.exists(), reason=f"real 2pt-FITS not on disk: {_REALDATA}" +) +def test_realdata_shim_equals_2pt_like(tmp_path, m_shim, m_2pt): + """On a real product, the shim on its SACC equals 2pt_like on the FITS. + + Builds a ξ-only analysis SACC from the real 2pt-FITS's own ξ± values and + covariance sub-block (via the converter test's ``_sacc_from_2pt_fits``, then + strip to ξ±), writes it, converts it to a plain-ξ FITS, and runs both engines + against the synthetic theory block. Scoping to ξ± isolates the shear + likelihood equality on the true (20-point-per-sign) data-vector shape — the + IA-only inference scope this PR targets — and keeps ``2pt_like``'s + ``twopoint.from_fits`` from tripping over the real file's separate + COVMAT_CELL / τ blocks. + """ + from astropy.io import fits + + from sp_validation.tests.test_sacc_io_realdata import _sacc_from_2pt_fits + + with fits.open(_REALDATA) as hdul: + full_s, _rho_hdu, _tau_hdu = _sacc_from_2pt_fits(hdul) + + # Rebuild a ξ-only SACC: same n(z), ξ± values and ξ± covariance sub-block. + source = sacc_io.source_name(0) + z, nz = sacc_io.get_nz(full_s, 0) + theta, xip, xim = sacc_io.get_xi(full_s, (0, 0), grid="reporting") + xi_idx = np.concatenate( + [ + full_s.indices(sacc_io.XI_PLUS, (source, source)), + full_s.indices(sacc_io.XI_MINUS, (source, source)), + ] + ) + xi_cov = full_s.covariance.dense[np.ix_(xi_idx, xi_idx)] + + s = sacc_io.new_sacc({0: (z, nz)}) + sacc_io.add_xi(s, (0, 0), theta, xip, xim, grid="reporting") + s.add_covariance(xi_cov) + + sacc_path = str(tmp_path / "real_xi.sacc") + sacc_io.save(s, sacc_path, type="mock") + fits_path = str(tmp_path / "real_xi_2pt.fits") + sacc_io.sacc_to_twopoint_fits(sacc_io.load(sacc_path), fits_path, n_bins=1) + + like_s, chi2_s, theory_s, n_s = _run(m_shim, _shim_opts(sacc_path)) + like_t, chi2_t, theory_t, n_t = _run(m_2pt, _twopt_opts(fits_path)) + + assert n_s == n_t + assert n_s == 2 * len(theta) + np.testing.assert_allclose(chi2_s, chi2_t, rtol=1e-10) + np.testing.assert_allclose(like_s, like_t, rtol=1e-10) + np.testing.assert_allclose(theory_s, theory_t, rtol=1e-10) diff --git a/uv-overrides.txt b/uv-overrides.txt new file mode 100644 index 00000000..813ad055 --- /dev/null +++ b/uv-overrides.txt @@ -0,0 +1,21 @@ +# uv dependency overrides — pass via `--overrides uv-overrides.txt` (or +# UV_OVERRIDE=uv-overrides.txt) to every `uv pip install` against this project. +# +# Why this file exists: firecrown declares its sampler *connectors* as hard +# dependencies, but we use firecrown only as the theory engine for Smokescreen +# blinding (`compute_theory_vector`); sampling stays with CosmoSIS in +# cosmo_inference. Of the three connector deps: +# +# - numcosmo-py exists only on conda-forge, so pip/uv resolution of firecrown +# is *impossible* without an override; +# - cosmosis ships sdist-only (full Fortran/C build with gsl/cfitsio) — a +# heavy, fragile compile in every CI image build, for a connector we never +# import; +# - cobaya is wheel-clean but equally unused. +# +# Each line below replaces the package's requirement (wherever it appears in +# the graph) with one gated on an always-false marker, dropping it from +# resolution. firecrown's likelihood/CCL core imports none of them. +numcosmo-py; python_version < "3" +cosmosis; python_version < "3" +cobaya; python_version < "3" diff --git a/workflow/Snakefile b/workflow/Snakefile index 6f59827b..0a1f3258 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -37,6 +37,10 @@ wildcard_constraints: # Compute rules (infrastructure — raw outputs, no evidence.json) include: "rules/twopoint.smk" +# Smokescreen blind-at-birth custody (blind_init / blind_part). Generic over the +# blindable parts twopoint.smk produces; dormant unless a data run requests a +# blinded part. Included before covariance/cosmo_val so their consumers resolve. +include: "rules/blinding.smk" include: "rules/covariance.smk" include: "rules/inference.smk" include: "rules/masks.smk" diff --git a/workflow/common.py b/workflow/common.py index 1a168041..78224d02 100644 --- a/workflow/common.py +++ b/workflow/common.py @@ -5,6 +5,8 @@ import re from pathlib import Path +from snakemake.io import temp + # Absolute path to the generic workflow's scripts, anchored on this module's own # location (common.py lives in workflow/, is `from common import *`'d into every # Snakefile, and so resolves to the generic workflow dir of the running checkout @@ -51,6 +53,10 @@ # glass-mock A/B/C variant, not Smokescreen blinding — see BLINDS above. "blind": r"[ABC]", "nbins": r"\d+", + # Constrained so a producer's ξ± output pattern cannot greedily absorb the + # "_blinded" suffix into npatch (which would make it ambiguous with the + # blind_part rule's {stem}_blinded output). npatch is always an integer. + "npatch": r"\d+", "min_sep": r"[0-9.]+", "max_sep": r"[0-9.]+", "gaussian": r"(g|ng)", @@ -65,16 +71,22 @@ DEFAULT_MASK_SUFFIX = "" CATALOG_CONFIG = None PLANCK18 = None +# Run type gates Smokescreen blind-at-birth (see the blind custody section +# below): "data" blinds the three blindable parts and binds ξ-derived consumers +# to the blinded siblings; "mock" bypasses blinding entirely. Set from +# config["cosmo_val"]["type"] in configure(); "data" is the production default. +RUN_TYPE = "data" def configure(workflow_config): """Install config-derived values after Snakemake has loaded configfiles.""" - global CATALOG_CONFIG, DEFAULT_MASK_SUFFIX, FIDUCIAL, PLANCK18 + global CATALOG_CONFIG, DEFAULT_MASK_SUFFIX, FIDUCIAL, PLANCK18, RUN_TYPE CATALOG_CONFIG = workflow_config FIDUCIAL = workflow_config["fiducial"] DEFAULT_MASK_SUFFIX = ( "_masked" if workflow_config["covariance"].get("default_masked", False) else "" ) + RUN_TYPE = workflow_config.get("cosmo_val", {}).get("type", "data") with open(COSMOLOGY_PARAMS) as f: PLANCK18 = json.load(f) @@ -196,6 +208,114 @@ def get_shear_catalog(wildcards): return str(Path(subdir) / shear_path) +# --------------------------------------------------------------------------- +# Smokescreen blind-at-birth custody (issues #247/#252, PR #253) +# --------------------------------------------------------------------------- +# Distinct from the glass-mock A/B/C `blind` wildcard above: this is Smokescreen +# concealment (the concealed=True SACC stamp). Blinding is per part, at birth. +# The three blindable parts (reporting ξ±, integration ξ±, pseudo-Cℓ) are each +# concealed the moment they are computed, so only blinded parts persist on disk. +# A `data` run binds every ξ-derived consumer (the terminal assemble, and the +# born-blinded COSEBIs / pure-E/B) to the *_blinded parts, which pulls the +# blind_part → blind_init subgraph into the DAG. A `mock` run bypasses blinding +# entirely and binds to the plaintext parts, so the subgraph never appears. +# RUN_TYPE is the single switch that flips which files exist. +# +# The path helpers below MIRROR sp_validation.blinding.init_paths / part_paths +# by hand rather than importing blinding (which pulls in numpy + smokescreen) at +# DAG-build time. test_blinding_wiring asserts the two stay in lockstep. + + +def run_type(): + """The campaign's run type, ``"data"`` or ``"mock"``. + + Every part writer stamps this as the SACC ``type`` metadata, which is what + ``blinding.assert_consistent_blind`` reads at assembly. A function + rather than the ``RUN_TYPE`` global because ``from common import *`` binds + names before ``configure()`` runs, so only a call reads the configured + value. + """ + return RUN_TYPE + + +def is_data_run(): + """True when blinding is active (production data runs); False for mocks.""" + return run_type() == "data" + + +def blind_state_dir(version): + """Per-version blind-init custody directory (commitment + encrypted seed).""" + return str(COSMO_VAL / "blind" / version) + + +def blind_state_paths(version): + """The fixed custody-state files blind_init writes for a version. + + Mirrors sp_validation.blinding.init_paths(blind_state_dir(version)). + """ + d = blind_state_dir(version) + return { + "commitment": os.path.join(d, "commitment.json"), + "bundle": os.path.join(d, "blind_seed.encrpt"), + "key": os.path.join(d, "blind_seed.key"), + } + + +def commitment_input(version): + """Input mapping binding a version's commitment.json, on data runs only. + + A part writer stamps its output concealed from that file (sacc_io.save's + `commitment=`), which is what lets a born-blinded or blind-irrelevant part + clear the fail-closed load gate the terminal assembly opens every part + through. A mock run binds nothing and the part stays plaintext. + """ + if not is_data_run(): + return {} + return {"commitment": blind_state_paths(version)["commitment"]} + + +def blinded_path(part_path): + """The *_blinded sibling blind_part writes beside a plaintext part. + + Mirrors sp_validation.blinding.part_paths(part_path)["blinded"]. + """ + stem, ext = os.path.splitext(str(part_path)) + return f"{stem}_blinded{ext or '.fits'}" + + +def version_of(stem): + """Extract the catalogue version embedded in a blindable part's stem. + + Every blindable stem carries the version (as {version}_xi_… or + pseudo_cl_{version}_…); blind_part needs it to locate the version's blind + state. Matches the shared `version` wildcard pattern. + """ + m = re.search(WILDCARD_CONSTRAINTS["version"], stem) + if m is None: + raise ValueError(f"no catalogue version found in part stem {stem!r}") + return m.group(0) + + +def blindable_part(part_path): + """On-disk path a run persists for one blindable part. + + Data run -> the blinded sibling (binding it pulls blind_part + blind_init + into the DAG); mock run -> the plaintext part (blinding bypassed). + """ + return blinded_path(part_path) if is_data_run() else str(part_path) + + +def maybe_temp(part_path): + """Wrap a producer's blindable plaintext part temp() on data runs. + + On a data run the plaintext part's only consumer is blind_part, which + escrows the true vector before Snakemake removes the temp file — so no + plaintext blindable part persists. On a mock run the part is the real + product downstream binds to, so it is left persistent. + """ + return temp(str(part_path)) if is_data_run() else str(part_path) + + # --------------------------------------------------------------------------- # CosmologyValidation diagnostic suite (cosmo_val.py) # --------------------------------------------------------------------------- @@ -254,6 +374,9 @@ def cv_init_params(config, version_list=None): nrandom_cell=cv["nrandom_cell"], cell_method=cv["cell_method"], nside_mask=cv["nside_mask"], + # Stamped as the SACC `type` of every part the cv writes; RUN_TYPE reads + # from this same key (see configure()). + run_type=cv.get("type", "data"), ) if cv.get("path_onecovariance"): params["path_onecovariance"] = cv["path_onecovariance"] diff --git a/workflow/rules/blinding.smk b/workflow/rules/blinding.smk new file mode 100644 index 00000000..a5f0dedb --- /dev/null +++ b/workflow/rules/blinding.smk @@ -0,0 +1,74 @@ +# Smokescreen blind-at-birth custody rules (issues #247/#252, PR #253). +# +# Two rules realise the three-verb custody surface of sp_validation.blinding on +# the DAG. They only enter the graph on a `data` run, and only when a consumer +# binds to a *_blinded part through common.blindable_part — a `mock` run never +# requests a blinded file, so blind_part and blind_init stay dormant. +# +# blind_init (once per catalogue version) draws the seed, publishes +# commitment.json + the encrypted seed bundle. +# blind_part (once per blindable part, at birth) conceals the part, escrows +# the true vector beside the blinded output, and lets Snakemake +# remove the plaintext (temp()) once it is the sole consumer. +# +# The terminal assemble_sacc rule (cosmo_val.smk) asserts the shared commitment +# across parts — the assembly-time custody check of #252. + +# The three blindable stems: reporting ξ± (rule xi), integration ξ± (xi_highres), +# and the analysis pseudo-Cℓ (pseudo_cl). None contains "_blinded", so the +# generic blind_part rule can never blind its own output twice. The version +# pattern is the shared one from common.WILDCARD_CONSTRAINTS. +_V = WILDCARD_CONSTRAINTS["version"] +BLINDABLE_STEM = ( + rf"(?:{_V}_xi_reporting_minsep=[0-9.]+_maxsep=[0-9.]+_nbins=\d+_npatch=\d+" + rf"|{_V}_xi_integration" + rf"|pseudo_cl_{_V}_blind=[ABC]_[a-z]+_nbins=\d+)" +) + + +rule blind_init: + """Fix the blind for one catalogue version (blind-init). + + Draws an OS-entropy seed, writes the repo-committable commitment.json + (seed commitment + config digest + the installed fork's draw scheme) and the + Fernet-encrypted seed bundle. Runs once per version and refuses to overwrite + existing state — a blind is a one-shot custody event. + """ + output: + commitment=str(COSMO_VAL / "blind" / "{version}" / "commitment.json"), + bundle=str(COSMO_VAL / "blind" / "{version}" / "blind_seed.encrpt"), + key=str(COSMO_VAL / "blind" / "{version}" / "blind_seed.key"), + params: + blind_dir=lambda w: blind_state_dir(w.version), + resources: + runtime=5, + script: + "../scripts/blind_init.py" + + +rule blind_part: + """Blind one intermediate part SACC at birth (blind-part). + + Conceals the plaintext part through its matching theory backend, escrows the + true vector into a per-part encrypted bundle beside the blinded output, and + leaves the plaintext for Snakemake to remove (it is a temp() output of the + producing rule, and this is its only consumer on a data run). Generic over + the three blindable stems. + """ + input: + part=str(COSMO_VAL / "{stem}.sacc"), + commitment=lambda w: blind_state_paths(version_of(w.stem))["commitment"], + bundle=lambda w: blind_state_paths(version_of(w.stem))["bundle"], + key=lambda w: blind_state_paths(version_of(w.stem))["key"], + output: + blinded=str(COSMO_VAL / "{stem}_blinded.sacc"), + escrow=str(COSMO_VAL / "{stem}_escrow.encrpt"), + escrow_key=str(COSMO_VAL / "{stem}_escrow.key"), + wildcard_constraints: + stem=BLINDABLE_STEM, + params: + blind_dir=lambda w: blind_state_dir(version_of(w.stem)), + resources: + runtime=10, + script: + "../scripts/blind_part.py" diff --git a/workflow/rules/cosmo_val.smk b/workflow/rules/cosmo_val.smk index 0cac5126..049b30b8 100644 --- a/workflow/rules/cosmo_val.smk +++ b/workflow/rules/cosmo_val.smk @@ -368,12 +368,32 @@ rule cv_pseudo_cl: # Pure E/B modes and COSEBIs (per version), then the B-mode summary # --------------------------------------------------------------------------- +# On a data run the COSEBIs / pure-E/B parts re-derive their E-mode vector from +# the *blinded* integration ξ± (COSEBIs) or blinded reporting + integration ξ± +# (pure-E/B). blindable_part returns the plaintext part on a mock run and the +# blinded part on a data run, so the ξ± inputs bind unconditionally for every +# version; see common.commitment_input for the commitment. +def cv_cosebis_inputs(w): + return { + "xi": cv_xi_txt(w.version), + "xi_integration": blindable_part(cv_xi_integration_sacc(w.version)), + **commitment_input(w.version), + } + + +def cv_pure_eb_inputs(w): + return { + "xi": cv_xi_txt(w.version), + "xi_reporting": blindable_part(cv_xi_reporting_sacc(w.version)), + "xi_integration": blindable_part(cv_xi_integration_sacc(w.version)), + **commitment_input(w.version), + } + + rule cv_pure_eb: """Pure E/B-mode decomposition for one version (config-space).""" input: - xi=lambda w: cv_xi_txt(w.version), - xi_reporting=lambda w: cv_xi_reporting_sacc(w.version), - xi_integration=lambda w: cv_xi_integration_sacc(w.version), + unpack(cv_pure_eb_inputs), output: npz=cv_pure_eb_npz("{version}"), sacc=cv_pure_eb_sacc("{version}"), @@ -396,8 +416,7 @@ rule cv_pure_eb: rule cv_cosebis: """COSEBIs E/B decomposition for one version (config-space, fine binning).""" input: - xi=lambda w: cv_xi_txt(w.version), - xi_integration=lambda w: cv_xi_integration_sacc(w.version), + unpack(cv_cosebis_inputs), output: npz=cv_cosebis_npz("{version}"), sacc=cv_cosebis_sacc("{version}"), @@ -488,14 +507,23 @@ def cv_assemble_inputs(version): fiducial harmonic tag). pseudo_cl (+ its cov) is included only when the config toggles the harmonic-space BB into the analysis. """ + # blindable_part binds the raw-signal parts (reporting ξ±, analysis pseudo-Cℓ) + # to their blinded siblings on a data run and to the plaintext on a mock run. + # COSEBIs and pure-E/B are born blinded (their writers derive them from the + # blinded integration ξ± and stamp concealed=True), so they bind by their own + # name in both cases; ρ/τ is a diagnostic carrying no cosmological vector but + # is stamped concealed pass-through so the fail-closed load gate admits it on + # a data run. The integration-grid ξ± itself is blinded at birth (see rule + # xi_highres) but is NOT gathered into the terminal file — it persists as its + # own per-part intermediate (see #247 ruling), consumed by COSEBIs/pure-E/B. parts = dict( - xi_reporting=cv_xi_reporting_sacc(version), + xi_reporting=blindable_part(cv_xi_reporting_sacc(version)), cosebis=cv_cosebis_sacc(version), pure_eb=cv_pure_eb_sacc(version), rho_tau=cv_rho_tau_sacc(version), ) if CV.get("include_pseudo_cl", False): - parts["pseudo_cl"] = cv_pseudo_cl_analysis_sacc(version) + parts["pseudo_cl"] = blindable_part(cv_pseudo_cl_analysis_sacc(version)) parts["pseudo_cl_cov"] = cv_pseudo_cl_cov(version) return parts diff --git a/workflow/rules/inference.smk b/workflow/rules/inference.smk index ff0c982d..40aba513 100644 --- a/workflow/rules/inference.smk +++ b/workflow/rules/inference.smk @@ -1,15 +1,19 @@ -# Imports from Snakefile: FIDUCIAL, COSMO_INFERENCE, COSMO_VAL, covariance_path, build_redshift_path, fiducial_binning_suffix -# NOTE: dormant subsystem. The file-name plumbing (config-driven paths + the -# producer-tagged pseudo-Cl names) is fixed and the DAG is valid, but it has not -# been run end-to-end. Reviving it still needs the FITS-CONTENT plumbing -# reconciled: cosmosis_fitting.py reads ELL/EE/BB + COVAR_FULL, while the -# producers write PSEUDO_CELL/ELL + COVAR_BB_BB. +# Imports from common (via `from common import *`): FIDUCIAL, COSMO_INFERENCE, +# COSMO_VAL, WORKFLOW_SCRIPTS, covariance_path, build_redshift_path, +# fiducial_binning_suffix. cv_analysis_sacc arrives from cosmo_val.smk (resolved +# lazily at DAG time, since that file is included after this one). +# +# Two paths live here: +# * Real-data inference_prep — LIVE (PR 7): consumes the assembled {version}.sacc +# and emits the converter 2pt-FITS + both engine inis (2pt_like, sacc_like). +# * glass-mock rules — still cosmosis_fitting.py-based (their SACC migration is +# out of scope); the pseudo-Cl file-name plumbing they depend on stays below. # Output root for CosmoSIS data products + configs. COSMO_INFERENCE (common.py) # already resolves to THIS repo's cosmo_inference dir, so the products land # beside the code that builds them rather than in a contributor's home. COSMO_INFERENCE_PROD = COSMO_INFERENCE -# Working directory for the cosmosis_fitting.py invocation — the same repo dir. +# Working directory for the (glass-mock) cosmosis_fitting.py invocation. COSMO_INFERENCE_RUNDIR = str(COSMO_INFERENCE) # External chain/mock locations are deployment-specific, so they live in config. @@ -56,88 +60,115 @@ def pseudo_cl_assets(version): return str(cl_path), str(cov_path) # --------------------------------------------------------------------------- -# DORMANT — pre-SACC cosmosis assembly. Migration to native SACC deferred to -# PR 7 (native-SACC inference consumption); do NOT deep-migrate here. +# Real-data inference prep — LIVE (native SACC, PR 7). Consumes the assembled +# analysis {version}.sacc (cosmo_val.smk's assemble_sacc rule) and emits the two +# file-prep products the A_ia (IA-only, ξ±) fiducial pipeline needs: +# (a) the converter 2pt-FITS (sacc_to_twopoint_fits) + a generated 2pt_like ini +# — the validating/legacy path (retiring cosmosis_fitting.py's assembly), +# (b) a generated sacc_like ini pointing at the SACC directly — the native path +# validated bit-for-bit against (a) (test_sacc_like.py). +# The converter is A_ia-scoped: no rho/tau sidecars, so it emits a pure-ξ FITS +# (it ignores the SACC's extra data types). This is file-prep only — the actual +# CosmoSIS sampling still runs via pipeline.sh against these products. # -# The SACC migration (PR 4) removed the data products several of these inputs -# name, so this rule's DAG no longer resolves and is NOT reachable from the -# cosmo_val suite (cosmo_val_all never requests it). Stale inputs: -# - xi_plus / xi_minus FITS: the `xi` rule now emits the reporting ξ± SACC part -# ({version}_xi_reporting_...sacc), not per-sign FITS. -# - pseudo_cl / pseudo_cl_cov via pseudo_cl_assets(): the `pseudo_cl` rule now -# writes .sacc (pseudo_cl_assets still requests .fits). -# PR 7 rewires this to consume the assembled {version}.sacc (built by -# cosmo_val.smk's assemble_sacc rule) directly, retiring cosmosis_fitting.py's -# per-product FITS assembly. Until then the inference target is knowingly red. +# The glass-mock rules below stay cosmosis_fitting.py-based; their SACC migration +# is out of scope for PR 7. # --------------------------------------------------------------------------- +# Generated per-version configs land in the (env-overridable) output root. +INFERENCE_CONFIG_OUT = COSMO_INFERENCE_PROD / "cosmosis_config" +# The ini TEMPLATES are source files: anchor them on the running checkout (repo +# root = the workflow dir's parent, via WORKFLOW_SCRIPTS), NOT on the output root +# — so a template edit in this checkout drives the DAG even when COSMO_INFERENCE +# points elsewhere. In a normal (non-worktree) run the two roots coincide. +INFERENCE_TEMPLATE_DIR = ( + Path(os.path.dirname(WORKFLOW_SCRIPTS)).parent / "cosmo_inference" / "cosmosis_config" +) + + +def _csl_dir(): + """The CSL checkout that fills COSMOSIS_DIR / sacc_like csl_dir in the inis. + + Read lazily (at DAG time, inside inference_prep's params) rather than at + module parse time: inference.smk is included by every paper workflow, but + only papers that run inference (cosmo_val) carry inference.csl_dir. A missing + key still fails loudly — just when the real-data inference is actually built, + not when an unrelated (bmodes) workflow merely parses this file. + """ + return INFERENCE["csl_dir"] + + rule inference_prep: input: - # Processed covariance matrix - use centralized covariance_path() - cov_matrix=lambda w: covariance_path(w.version, w.blind, min_sep=w.min_sep, max_sep=w.max_sep, nbins=w.nbins), - # Xi FITS files — PRE-SACC (no longer produced; see dormant note above) - xi_plus=str(COSMO_VAL / "xi_plus_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), - xi_minus=str(COSMO_VAL / "xi_minus_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), - # n(z) file (using new location with base version mapping) - nz_file=lambda w: build_redshift_path(w.version, w.blind), - # rho/tau stats - rho_stats=str(COSMO_VAL / "rho_tau_stats/rho_stats_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), - tau_stats=str(COSMO_VAL / "rho_tau_stats/tau_stats_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), - # tau covariance (tracked as dependency) - tau_cov=str(COSMO_VAL / "rho_tau_stats/cov_tau_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}_th.npy"), - # pseudo_cl / pseudo_cl_cov — PRE-SACC (.fits path; producer now writes .sacc) - pseudo_cl=lambda w: pseudo_cl_assets(w.version)[0], - pseudo_cl_cov=lambda w: pseudo_cl_assets(w.version)[1], + # The terminal assembled analysis SACC (cosmo_val.smk assemble_sacc). Bound + # lazily through its helper so the filename tracks that rule, not a literal. + sacc=lambda w: cv_analysis_sacc(w.version), + # The two pipeline ini templates are static repo files, but binding them as + # inputs (not params) puts them in the DAG, so editing a template + # regenerates the configs rather than leaving stale output on disk. + template_2pt=str(INFERENCE_TEMPLATE_DIR / "cosmosis_pipeline_A_ia.ini"), + template_sacc=str(INFERENCE_TEMPLATE_DIR / "cosmosis_pipeline_A_ia_sacc.ini"), output: - fits_file=str( - COSMO_INFERENCE_PROD - / "data/{version}_{blind}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}/cosmosis_{version}_{blind}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits" + fits_file=str(COSMO_INFERENCE_PROD / "data/{version}/cosmosis_{version}.fits"), + config_file_2pt=str( + INFERENCE_CONFIG_OUT / "cosmosis_pipeline_{version}_A_ia.ini" + ), + config_file_sacc=str( + INFERENCE_CONFIG_OUT / "cosmosis_pipeline_{version}_A_ia_sacc.ini" ), - config_file=str( - COSMO_INFERENCE_PROD - / "cosmosis_config/cosmosis_pipeline_{version}_{blind}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.ini" - ) params: - cosmosis_root="{version}_{blind}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}", - data_dir=f"{CHAINS_DIR}/{{version}}_{{blind}}_minsep={{min_sep}}_maxsep={{max_sep}}_nbins={{nbins}}_npatch={{npatch}}", - output_root=str(COSMO_INFERENCE_PROD), + # SCRATCH = the per-version chain output root the generated inis point at. + scratch=lambda w: f"{CHAINS_DIR}/{w.version}", + cosmosis_dir=lambda w: _csl_dir(), threads: 1 resources: mem_mb=8000, runtime=10, - shell: - """ - cd {COSMO_INFERENCE_RUNDIR} + run: + import os + import sys - # Run inference preparation step with cosmosis_fitting.py - python scripts/cosmosis_fitting.py \ - --cosmosis-root {params.cosmosis_root} \ - --nz-file {input.nz_file} \ - --data-dir {params.data_dir} \ - --output-root {params.output_root} \ - --xi {input.xi_plus} {input.xi_minus} \ - --cov-xi {input.cov_matrix} \ - --use-rho-tau \ - --rho-stats {input.rho_stats} \ - --tau-stats {input.tau_stats} \ - --cov-tau {input.tau_cov} \ - --cl-file {input.pseudo_cl} \ - --cov-cl {input.pseudo_cl_cov} - """ + from sp_validation import sacc_io + from sp_validation.sacc_io import sacc_to_twopoint_fits + + os.makedirs(os.path.dirname(output.fits_file), exist_ok=True) + + # (a) converter 2pt-FITS — pure ξ (A_ia scope; no rho/tau sidecars). + sacc_to_twopoint_fits(sacc_io.load(input.sacc), output.fits_file, n_bins=1) + + # (b) + (c) the two generated pipeline inis, from the template inputs. + # WORKFLOW_SCRIPTS (common.py) is the absolute generic-workflow scripts dir. + sys.path.insert(0, WORKFLOW_SCRIPTS) + from generate_inference_config import ( + _substitutions, + generate_inference_config, + ) + + generate_inference_config( + input.template_2pt, + output.config_file_2pt, + _substitutions( + scratch=params.scratch, + cosmosis_dir=params.cosmosis_dir, + fits_file=output.fits_file, + ), + ) + generate_inference_config( + input.template_sacc, + output.config_file_sacc, + _substitutions( + scratch=params.scratch, + cosmosis_dir=params.cosmosis_dir, + sacc_file=input.sacc, + ), + ) rule inference_fiducial: input: - # Use the same output patterns as inference_prep with FIDUCIAL params - rules.inference_prep.output.fits_file.format( - version=FIDUCIAL["version"], blind=FIDUCIAL["blind"], - min_sep=FIDUCIAL["min_sep"], max_sep=FIDUCIAL["max_sep"], - nbins=FIDUCIAL["nbins"], npatch=FIDUCIAL["npatch"] - ), - rules.inference_prep.output.config_file.format( - version=FIDUCIAL["version"], blind=FIDUCIAL["blind"], - min_sep=FIDUCIAL["min_sep"], max_sep=FIDUCIAL["max_sep"], - nbins=FIDUCIAL["nbins"], npatch=FIDUCIAL["npatch"] - ) + # The fiducial version's prep products (both engine inis + the FITS). + rules.inference_prep.output.fits_file.format(version=FIDUCIAL["version"]), + rules.inference_prep.output.config_file_2pt.format(version=FIDUCIAL["version"]), + rules.inference_prep.output.config_file_sacc.format(version=FIDUCIAL["version"]), rule inference_glass_mocks: diff --git a/workflow/rules/twopoint.smk b/workflow/rules/twopoint.smk index 69c516e9..31b1dd5a 100644 --- a/workflow/rules/twopoint.smk +++ b/workflow/rules/twopoint.smk @@ -14,7 +14,9 @@ rule xi: # rule to share one wildcard set, and it keeps the reporting .sacc name # self-describing so requesting it binds the xi job unambiguously. txt=str(COSMO_VAL / "{version}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.txt"), - xi_reporting=str(COSMO_VAL / "{version}_xi_reporting_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.sacc"), + # Blindable part: temp() on a data run so only its blinded sibling + # persists (blind_part escrows the true vector first). See common.maybe_temp. + xi_reporting=maybe_temp(str(COSMO_VAL / "{version}_xi_reporting_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.sacc")), threads: 24 params: ver="{version}", @@ -22,6 +24,9 @@ rule xi: max_sep="{max_sep}", nbins="{nbins}", npatch="{npatch}", + # Stamped as the part's SACC `type` — custody state at assembly + # (see blinding.assert_consistent_blind). + type=run_type(), resources: mem_mb=30000, disk_mb=20000, @@ -69,7 +74,9 @@ rule xi_highres: # [0.08, 300] at 1000 bins) so the single part serves both consumers: # pure-E/B needs it to strictly contain its reporting grid down to 0.08; # COSEBIs scale-cuts on the same part. Decoupled from covariance.smk. - xi_integration=str(COSMO_VAL / "{version}_xi_integration.sacc"), + # Blindable part: temp() on a data run so blind_part produces the _blinded + # sibling the COSEBIs/pure-E/B consumers bind (see rule xi / common.maybe_temp). + xi_integration=maybe_temp(str(COSMO_VAL / "{version}_xi_integration.sacc")), params: version="{version}", cat_config=CAT_CONFIG, @@ -78,6 +85,7 @@ rule xi_highres: nbins=_INTEGRATION["nbins"], out=str(COSMO_VAL), scripts=WORKFLOW_SCRIPTS, + run_type=run_type(), threads: 24 resources: mem_mb=40000, @@ -86,7 +94,8 @@ rule xi_highres: "python {params.scripts}/run_2pcf_highres.py " "--version {params.version} --cat-config {params.cat_config} " "--min-sep {params.min_sep} --max-sep {params.max_sep} " - "--nbins {params.nbins} --npatch 1 --out {params.out}" + "--nbins {params.nbins} --npatch 1 --out {params.out} " + "--run-type {params.run_type}" rule run_cosmo_val: @@ -108,6 +117,10 @@ rule run_cosmo_val: rule rho_tau_stats: + # ρ/τ has no blindable input; it binds only the commitment, to stamp its + # part concealed pass-through (see common.commitment_input). + input: + unpack(lambda w: commitment_input(w.version)), output: rho_stats=str(COSMO_VAL / "rho_tau_stats/rho_stats_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), tau_stats=str(COSMO_VAL / "rho_tau_stats/tau_stats_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), @@ -123,6 +136,7 @@ rule rho_tau_stats: max_sep="{max_sep}", nbins="{nbins}", npatch="{npatch}", + type=run_type(), resources: mem_mb=30000, disk_mb=20000, @@ -145,6 +159,13 @@ rule pseudo_cl: concealed=True stamp is a separate axis on the SACC file. """ output: + # This generic rule produces every pseudo-Cℓ variant — the analysis part + # (blind=A, powspace, nbins=32) folded into {version}.sacc, plus the fine + # (COSEBIS) and glass-mock variants. Only the analysis part is a terminal + # blindable, and a data run blinds it via a requested _blinded sibling + # (blind_part reads this plaintext); the fine/mock variants are B-mode / + # validation intermediates left untouched here. The output is therefore + # not temp()'d — see the PR note on residual unblinded pseudo-Cℓ. pseudo_cl=str(COSMO_VAL / "pseudo_cl_{version}_blind={blind}_{binning}_nbins={nbins}.sacc"), wildcard_constraints: blind="[ABC]", # glass-mock variant, not Smokescreen blinding diff --git a/workflow/scripts/assemble_sacc.py b/workflow/scripts/assemble_sacc.py index cfad6a8f..4c7d41a1 100644 --- a/workflow/scripts/assemble_sacc.py +++ b/workflow/scripts/assemble_sacc.py @@ -7,9 +7,13 @@ Each per-statistic ``*.sacc`` *part* (written born-as-SACC by the mixins and the run_2pcf / generate_pseudo_cl scripts) holds one statistic. The assembler loads them in canonical order — ξ± reporting, pseudo-Cℓ, COSEBIs, pure-E/B, ρ/τ — and -calls :func:`sacc_writers.assemble_analysis_sacc`, which rebuilds one Sacc with a -single ``BlockDiagonalCovariance`` (point-insertion order = block order, -validated by ``sacc_io.assemble_covariance``). +hands them to :func:`sacc_io.gather` along with +:func:`sacc_writers.assemble_analysis_sacc` as the assembly, which rebuilds one +Sacc with a single ``BlockDiagonalCovariance`` (point-insertion order = block +order, validated by ``sacc_io.assemble_covariance``). Going through ``gather`` +rather than calling the assembler directly is what puts this path behind the +blind-custody gate (``blinding.assert_consistent_blind``) — there is one +terminal seam, not one per assembler. Covariance sourcing (the part-by-part decision) ----------------------------------------------- @@ -176,7 +180,17 @@ def assemble_sacc( ) if not parts: raise ValueError(f"no parts found for {version}: {part_paths}") - s = assemble_analysis_sacc(nz, metadata, parts) + # Assembly-time custody assertion (#252) runs inside sacc_io.gather, the one + # terminal seam: every blindable part (ξ± / pseudo-Cℓ EE) must share one + # blind commitment + config digest + draw scheme, or assembly fails closed — + # mixed blinded/plaintext parts and divergent-seed parts both raise. ρ/τ and + # covariance-only parts are exempt. Gather also stamps the shared blind onto + # the assembled file, so the union reads as concealed in its own right. + # The production assembly (n(z) + per-part covariance blocks) is passed in; + # gather's own default merge() is for parts that already carry covariance. + s = sacc_io.gather( + parts, assemble=lambda ordered: assemble_analysis_sacc(nz, metadata, ordered) + ) # Assembly preserves its parts' provenance: every part was written by # sacc_io.save and therefore carries the type=data|mock stamp in its # metadata (copied into the assembled file above). diff --git a/workflow/scripts/blind_init.py b/workflow/scripts/blind_init.py new file mode 100644 index 00000000..f9e4b682 --- /dev/null +++ b/workflow/scripts/blind_init.py @@ -0,0 +1,12 @@ +"""Rule blind_init: fix the blind for one catalogue version. + +Thin wrapper over :func:`sp_validation.blinding.blind_init`. Draws the seed, +writes commitment.json + the encrypted seed bundle into the version's blind +directory. The plaintext seed is never written (the encryptor deletes it). +""" + +from snakemake.script import snakemake + +from sp_validation import blinding + +blinding.blind_init(snakemake.params["blind_dir"]) diff --git a/workflow/scripts/blind_part.py b/workflow/scripts/blind_part.py new file mode 100644 index 00000000..b93313fb --- /dev/null +++ b/workflow/scripts/blind_part.py @@ -0,0 +1,19 @@ +"""Rule blind_part: blind one intermediate part SACC at birth. + +Thin wrapper over :func:`sp_validation.blinding.blind_part`. Conceals the part, +escrows the true vector beside the blinded output, and leaves the plaintext in +place: it is a temp() output of the producing rule, so Snakemake removes it once +this (its only consumer on a data run) finishes. keep_input=True hands that +lifecycle to Snakemake rather than deleting inside the blind step, which keeps +the blinded output and its temp input in one consistent DAG accounting. +""" + +from snakemake.script import snakemake + +from sp_validation import blinding + +blinding.blind_part( + snakemake.input["part"], + snakemake.params["blind_dir"], + keep_input=True, +) diff --git a/workflow/scripts/cv_cosebis.py b/workflow/scripts/cv_cosebis.py index 84991d7b..ffb4fd0d 100644 --- a/workflow/scripts/cv_cosebis.py +++ b/workflow/scripts/cv_cosebis.py @@ -41,6 +41,12 @@ theta, xip, xim = sacc_io.get_xi(integ, (0, 0), grid="integration") en_part, _ = cosebis_from_xi(theta, xip, xim, p["nmodes"], scale_cut=fiducial_scale_cut) +# On a data run the integration ξ± part is blinded (blindable_part), so the +# derived En is blinded too; commitment.json (bound only on a data run) stamps +# the emitted COSEBIs part concealed. A mock run binds no commitment and the +# part stays plaintext/unconcealed. +commitment_path = snakemake.input.get("commitment") + # Born-as-SACC COSEBIs part at the fiducial scale cut (plot_cosebis stored the # multi-cut results on the instance); En comes from the consumed part. cv.cosebis_to_sacc_part( @@ -49,6 +55,7 @@ cv._cosebis_results[version], fiducial_scale_cut=fiducial_scale_cut, en_override=en_part, + commitment_path=commitment_path, ) # Overwrite the raw fiducial-cut En plot_cosebis wrote into the diagnostic npz with diff --git a/workflow/scripts/cv_pure_eb.py b/workflow/scripts/cv_pure_eb.py index d4797922..4d3ec244 100644 --- a/workflow/scripts/cv_pure_eb.py +++ b/workflow/scripts/cv_pure_eb.py @@ -44,8 +44,20 @@ ti, xpi, xmi = sacc_io.get_xi(integ, (0, 0), grid="integration") modes = pure_eb_from_xi(tr, xpr, xmr, ti, xpi, xmi, tmin, tmax) +# On a data run the reporting + integration ξ± parts are blinded (blindable_part), +# so the derived modes are blinded too; commitment.json (bound only on a data run) +# stamps the emitted pure-E/B part concealed. A mock run binds no commitment and +# the part stays plaintext/unconcealed. +commitment_path = snakemake.input.get("commitment") + # Born-as-SACC pure-E/B part; the six pure-mode blocks come from the consumed parts. -cv.pure_eb_to_sacc_part(version, snakemake.output["sacc"], results, eb_override=modes) +cv.pure_eb_to_sacc_part( + version, + snakemake.output["sacc"], + results, + eb_override=modes, + commitment_path=commitment_path, +) # Overwrite the raw pure-mode arrays plot_pure_eb wrote into the diagnostic npz with # the part-derived ones (identical to the SACC part's). theta / cov / PTE fields are diff --git a/workflow/scripts/generate_inference_config.py b/workflow/scripts/generate_inference_config.py new file mode 100644 index 00000000..68a8efe7 --- /dev/null +++ b/workflow/scripts/generate_inference_config.py @@ -0,0 +1,159 @@ +"""Generate a CosmoSIS pipeline ini by filling a template's ``[DEFAULT]`` section. + +Dual-mode, like ``assemble_sacc.py``. Under Snakemake (``script:`` directive) the +injected ``snakemake`` object supplies the template, output path and DEFAULT +substitutions; as a standalone CLI (argparse) the same fill runs from explicit +flags. + +The template carries ``%(KEY)s`` interpolation placeholders (SCRATCH, FITS_FILE +or SACC_FILE, COSMOSIS_DIR, SP_VALIDATION_MODULES) in its module sections; this +script prepends the concrete ``KEY = value`` lines into ``[DEFAULT]`` so +CosmoSIS's ConfigParser resolves them at load. It is deliberately plain text +processing — appending lines after the ``[DEFAULT]`` header, the same idiom as +``pipeline.sh``'s ``sed -i "/^\\[DEFAULT\\]/a\\KEY = value"`` — rather than a +configparser round-trip, which would strip the template's comments and its +``%(...)s`` interpolation. + +``SP_VALIDATION_MODULES`` is resolved from ``sp_validation.__file__``'s parent so +the generated ini points at the installed package's module directory (where +``sacc_like_unions.py`` lives) regardless of checkout location. +""" + +import argparse +from pathlib import Path + + +def _sp_validation_modules(): + """The directory holding the sp_validation CosmoSIS module files. + + Resolved from the installed package so the generated ini finds + ``sacc_like_unions.py`` wherever sp_validation is installed. + """ + import sp_validation + + return str(Path(sp_validation.__file__).resolve().parent) + + +def generate_inference_config(template_path, out_path, substitutions): + """Write ``out_path`` from ``template_path`` with ``substitutions`` in DEFAULT. + + Parameters + ---------- + template_path : str or Path + The pipeline ini template (carries ``%(KEY)s`` placeholders). + out_path : str or Path + Destination ini. + substitutions : dict + ``{KEY: value}`` lines prepended into the template's ``[DEFAULT]`` + section. Every referenced ``%(KEY)s`` in the template must have a value + here (COSMOSIS_DIR already sits in the template's DEFAULT and may be + overridden). ``None`` values are dropped (an absent optional key). + """ + lines = Path(template_path).read_text().splitlines(keepends=True) + + header = "[DEFAULT]" + default_idx = next( + (i for i, line in enumerate(lines) if line.strip() == header), None + ) + if default_idx is None: + raise ValueError(f"template {template_path} has no [DEFAULT] section") + + # The end of the DEFAULT section: the next `[section]` header, or EOF. + section_end = next( + ( + i + for i in range(default_idx + 1, len(lines)) + if lines[i].lstrip().startswith("[") + ), + len(lines), + ) + + # A key the template already declares in DEFAULT is REPLACED in place (e.g. the + # template's placeholder COSMOSIS_DIR); a genuinely-new key is prepended just + # after the header. This avoids a duplicate DEFAULT key, which CosmoSIS's + # ConfigParser (strict) rejects. + wanted = {key: value for key, value in substitutions.items() if value is not None} + remaining = dict(wanted) + for i in range(default_idx + 1, section_end): + stripped = lines[i].lstrip() + if not stripped or stripped.startswith(("#", ";", "[")): + continue + existing_key = stripped.split("=", 1)[0].strip() + if existing_key in remaining: + lines[i] = f"{existing_key} = {remaining.pop(existing_key)}\n" + + prepended = [f"{key} = {value}\n" for key, value in remaining.items()] + out_lines = lines[: default_idx + 1] + prepended + lines[default_idx + 1 :] + + out_path = Path(out_path) + out_path.parent.mkdir(parents=True, exist_ok=True) + out_path.write_text("".join(out_lines)) + print( + f"Wrote {out_path} from {template_path} " + f"({len(wanted)} DEFAULT keys, {len(prepended)} new)" + ) + return str(out_path) + + +def _substitutions(scratch, cosmosis_dir, *, fits_file=None, sacc_file=None): + """Assemble the DEFAULT substitution dict, resolving SP_VALIDATION_MODULES. + + ``fits_file`` (2pt_like path) and ``sacc_file`` (sacc_like path) are mutually + the data-file placeholder for their template; whichever the template + references is filled, the other left absent. + """ + return { + "SCRATCH": scratch, + "FITS_FILE": fits_file, + "SACC_FILE": sacc_file, + "COSMOSIS_DIR": cosmosis_dir, + "SP_VALIDATION_MODULES": _sp_validation_modules(), + } + + +def _from_snakemake(smk): + p = smk.params + generate_inference_config( + template_path=smk.input[0] + if not hasattr(smk.input, "template") + else smk.input.template, + out_path=str(smk.output[0]), + substitutions=_substitutions( + scratch=p["scratch"], + cosmosis_dir=p["cosmosis_dir"], + fits_file=p.get("fits_file", None), + sacc_file=p.get("sacc_file", None), + ), + ) + + +def _from_cli(argv=None): + ap = argparse.ArgumentParser( + description="Generate a CosmoSIS pipeline ini from a template + DEFAULT subs." + ) + ap.add_argument("--template", required=True, help="Pipeline ini template") + ap.add_argument("--out", required=True, help="Output ini path") + ap.add_argument("--scratch", required=True, help="SCRATCH value") + ap.add_argument("--cosmosis-dir", required=True, help="COSMOSIS_DIR value") + ap.add_argument("--fits-file", default=None, help="FITS_FILE (2pt_like path)") + ap.add_argument("--sacc-file", default=None, help="SACC_FILE (sacc_like path)") + a = ap.parse_args(argv) + generate_inference_config( + template_path=a.template, + out_path=a.out, + substitutions=_substitutions( + scratch=a.scratch, + cosmosis_dir=a.cosmosis_dir, + fits_file=a.fits_file, + sacc_file=a.sacc_file, + ), + ) + + +if __name__ == "__main__": + try: + snakemake # noqa: F821 — injected by Snakemake's script: directive + except NameError: + _from_cli() + else: + _from_snakemake(snakemake) # noqa: F821 diff --git a/workflow/scripts/run_2pcf.py b/workflow/scripts/run_2pcf.py index 3e513479..7f39660c 100644 --- a/workflow/scripts/run_2pcf.py +++ b/workflow/scripts/run_2pcf.py @@ -39,6 +39,7 @@ def run_2pcf( cat_config, output_dir, sacc_out=None, + run_type="data", ): """Measure ξ±(θ) for ``ver`` and write its reporting SACC part. @@ -49,7 +50,9 @@ def run_2pcf( ``cat_config['paths']['output']`` so the ``.txt`` byproduct lands where lc expects. ``sacc_out`` is the exact destination for the reporting ξ± SACC part (the Snakemake-declared output); it defaults to ``{ver}_xi_reporting.sacc`` - under the resolved output directory for the CLI path. + under the resolved output directory for the CLI path. ``run_type`` + (``"data"`` or ``"mock"``) is stamped as the part's SACC ``type`` — custody + state at assembly (see ``blinding.assert_consistent_blind``). Returns ------- @@ -85,7 +88,7 @@ def run_2pcf( out_path = sacc_out or os.path.join( output_dir or cv.cc["paths"]["output"], f"{ver}_xi_reporting.sacc" ) - sacc_io.save(s, out_path, type="data") + sacc_io.save(s, out_path, type=run_type) print(f"Wrote reporting ξ± SACC part: {out_path}") return gg @@ -107,6 +110,7 @@ def _from_snakemake(smk): # Write the SACC part exactly where the rule declares it (the .txt # byproduct still lands under the resolved output dir via _output_path). sacc_out=smk.output["xi_reporting"], + run_type=p.get("type", "data"), ) @@ -133,6 +137,12 @@ def _from_cli(argv=None): "--cat-config", required=True, help="Absolute path to cat_config.yaml" ) ap.add_argument("--out", required=True, help="Output directory (lc {output})") + ap.add_argument( + "--run-type", + default="data", + choices=("data", "mock"), + help="Campaign run type stamped as the part's SACC `type`", + ) a = ap.parse_args(argv) run_2pcf( ver=a.ver, @@ -142,6 +152,7 @@ def _from_cli(argv=None): npatch=a.npatch, cat_config=a.cat_config, output_dir=a.out, + run_type=a.run_type, ) diff --git a/workflow/scripts/run_2pcf_highres.py b/workflow/scripts/run_2pcf_highres.py index 6049e105..985c0ef4 100644 --- a/workflow/scripts/run_2pcf_highres.py +++ b/workflow/scripts/run_2pcf_highres.py @@ -91,6 +91,9 @@ NPATCH = None OUTPUT_DIR = None PATCH_FILE = None +# Campaign run type, stamped as the part's SACC `type`. Custody state, not +# decoration (see blinding.assert_consistent_blind). +RUN_TYPE = "data" def parse_args(argv=None): @@ -121,6 +124,12 @@ def parse_args(argv=None): "--max-sep", type=float, default=300.0, help="Max separation [arcmin]" ) ap.add_argument("--out", required=True, help="Output directory (lc {output})") + ap.add_argument( + "--run-type", + default="data", + choices=("data", "mock"), + help="Campaign run type stamped as the part's SACC `type`", + ) return ap.parse_args(argv) @@ -233,7 +242,7 @@ def write_xi_integration_sacc(gg): variances=np.concatenate([gg.varxip, gg.varxim]), ) out_path = os.path.join(OUTPUT_DIR, f"{VERSION}_xi_integration.sacc") - sacc_io.save(s, out_path, type="data") + sacc_io.save(s, out_path, type=RUN_TYPE) log(f" Wrote {out_path}") @@ -297,10 +306,11 @@ def resolve_paths(ver): def main(): global CAT_PATH, VERSION, E1_COL, E2_COL, W_COL, REDSHIFT_PATH - global TMIN, TMAX, NBINS, NPATCH, OUTPUT_DIR, PATCH_FILE + global TMIN, TMAX, NBINS, NPATCH, OUTPUT_DIR, PATCH_FILE, RUN_TYPE args = parse_args() VERSION = args.version + RUN_TYPE = args.run_type NBINS = args.nbins NPATCH = args.npatch TMIN = args.min_sep diff --git a/workflow/scripts/run_rho_tau.py b/workflow/scripts/run_rho_tau.py index 9df0bbfd..7b64c418 100644 --- a/workflow/scripts/run_rho_tau.py +++ b/workflow/scripts/run_rho_tau.py @@ -44,9 +44,16 @@ theta_max=float(params["max_sep"]), nbins=int(params["nbins"]), npatch=int(params["npatch"]), + run_type=params.get("type", "data"), ) -cv.calculate_rho_tau_stats() +# On a data run the rule binds the version's commitment.json, which stamps the +# emitted ρ/τ part concealed pass-through — values untouched (ρ/τ carries no +# cosmological vector) but clearing the fail-closed load gate at assembly. A mock +# run binds no commitment and the part stays plaintext, stamped type='mock'. +commitment_path = snakemake.input.get("commitment") # type: ignore + +cv.calculate_rho_tau_stats(commitment_path=commitment_path) # Confirm CosmologyValidation produced the requested outputs. calculate_rho_tau_stats # writes the rho/tau FITS *and* the born-as-SACC rho_tau part (via