From b255098c79c5d4f4fe3bf4c400a5f550d00ddc00 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 00:41:56 +0200 Subject: [PATCH 01/17] =?UTF-8?q?feat(cosmology):=20TheoryConfig=20+=20the?= =?UTF-8?q?=20two=20independent=20=CE=BE=C2=B1=20theory=20paths?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The single home for the blinding/inference fiducial (TheoryConfig, frozen, fail-fast overrides; two stack-specific halofit tokens denoting one HMCode2020+feedback recipe; Neff/T_CMB emitted explicitly as module constants) and two ξ± routes: the CCL-native path (Boltzmann-CAMB P(k) + CCL Limber/FFTLog — the blinding backends' recipe) and an independent pycamb P(k) → Pk2D → CCL projection with the closed-form σ8-matched A_s rescale. numpy-only at module level; CCL/CAMB are function-local imports. Cosmology objects are cached per parameter point so the three blinding backends share two CAMB runs per blind. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/cosmology.py | 373 +++++++++++++++++++++++++++++++++ 1 file changed, 373 insertions(+) create mode 100644 src/sp_validation/cosmology.py diff --git a/src/sp_validation/cosmology.py b/src/sp_validation/cosmology.py new file mode 100644 index 00000000..91cc96fe --- /dev/null +++ b/src/sp_validation/cosmology.py @@ -0,0 +1,373 @@ +"""Fiducial theory configuration and the two ξ± theory paths. + +:Name: cosmology.py + +:Description: The single home for the blinding/inference fiducial cosmology + (:class:`TheoryConfig`) and two independent routes to the tomographic + 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" + + # 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).""" + fields = {f.name for f in dataclasses.fields(cls)} + unknown = set(overrides) - fields + if unknown: + raise ValueError( + f"unknown TheoryConfig fields {sorted(unknown)}; " + f"valid fields are {sorted(fields)}" + ) + return cls(**overrides) + + +# --------------------------------------------------------------------------- # +# 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 (``ccl_halofit_version``, + ``hmcode_logT_AGN``). The nonlinear P(k) is routed through CAMB's + HMCode2020 — the same recipe 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 CAMB runs across blocks). + """ + import pyccl as ccl + + key = ( + tuple(sorted(params.items())), + config.ccl_halofit_version, + config.hmcode_logT_AGN, + ) + if key not in _COSMO_CACHE: + _COSMO_CACHE[key] = ccl.Cosmology( + **params, + matter_power_spectrum="camb", + extra_parameters={ + "camb": { + "halofit_version": config.ccl_halofit_version, + "HMCode_logT_AGN": config.hmcode_logT_AGN, + } + }, + ) + 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 From ed2f109cfcbb72556d9e7e1c511d2ddfeaf76133 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 00:42:14 +0200 Subject: [PATCH 02/17] feat(blinding): fork-protocol blind/unblind over the master SACC MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rework of the blinding box onto the UNIONS-WL/Smokescreen fork protocol: the fork's ConcealDataVector(fiducial_params, shifts_dict, sacc_data, *, seed, theory_fn) is the concealment engine for every block, with its internal per-key-independent CCL-native draw — no local draw, no firecrown, no likelihood object. - Envelope calibration: the (S8, Ωm) box maps to {sigma8, Omega_c} half-widths at the fiducial (σ8 = S8/√(Ωm/0.3); Ω_c one-to-one). - Three theory_fn backends matching the master layout: coarse ξ±, fine ξ± (grid tags), pseudo-Cℓ_EE (W @ Cℓ_EE on the stored BandpowerWindows; ΔBB = ΔEB ≡ 0). Two-tracer everywhere; rows aligned to each extracted block's own order. - Extract → conceal → merge per block, with recorded master indices; COSEBIs/pure-E/B re-derived through the pipeline estimators (cosmo_numba kernels) on the blinded fine ξ±; covariance and ρ/τ untouched. - Custody: sha256(seed) + canonical config digest as a repo-committable commitment; seed + true vector escrowed in a Fernet bundle (smokescreen.encryption); seed_smokescreen stripped; unblind verifies both hashes before subtracting and restores the truth bit-for-bit. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/blinding.py | 740 ++++++++++++++++++++++++++++++++++ 1 file changed, 740 insertions(+) create mode 100644 src/sp_validation/blinding.py diff --git a/src/sp_validation/blinding.py b/src/sp_validation/blinding.py new file mode 100644 index 00000000..9629de0c --- /dev/null +++ b/src/sp_validation/blinding.py @@ -0,0 +1,740 @@ +"""Blinding — conceal the measured data vector behind a hidden cosmology. + +:Name: blinding.py + +:Description: Smokescreen blinding wiring for the master SACC data product. + The measured statistics are shifted 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)``). + + **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; each blindable block goes through its + ``ConcealDataVector(fiducial_params, shifts_dict, sacc_data, *, seed, + theory_fn)`` entry point. This module supplies the two + sp_validation-specific pieces: the three ``theory_fn`` backends matching + the master SACC row layout (coarse ξ±, fine ξ±, pseudo-Cℓ) and the SACC + assembly around the concealed vectors. + + **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 − Ω_ν``. The same seed yields + the same hidden cosmology across all three block passes (the fork's draw + depends only on ``(key, seed)``), so the three shifted blocks are + mutually consistent by construction. + + **Derived statistics are re-derived, never shifted.** COSEBIs and + pure-E/B are replaced by re-running the pipeline's own estimators + (the ``cosmo_numba`` kernels :mod:`sp_validation.b_modes` drives) on the + blinded fine ξ±. Covariance and ρ/τ PSF diagnostics are byte-identical + to the input — 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.** The blind CLI draws an + OS-entropy seed, blinds immediately, publishes ``sha256(seed)`` plus a + canonical config digest as a repo-committable commitment, encrypts the + seed + true vector into a Fernet bundle (``smokescreen.encryption``) and + deletes the plaintext. Unblind verifies both hashes against the + commitment before subtracting anything, then restores the true vector. +""" + +import dataclasses +import hashlib +import json +import os +import secrets + +import numpy as np + +from . import sacc_io +from .cosmology 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. The fork draws each key independently as + ``U(fid − h, fid + h)``. + """ + 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. Python's ``json`` emits floats via their + shortest round-trip ``repr``, so two runs of the same config produce + byte-identical digests. Checked (with ``sha256(seed)``) at unblind, + so a wrong envelope or a mismatched P(k) recipe cannot silently + subtract a wrong shift. + """ + payload = { + "s8_half_width": self.s8_half_width, + "omega_m_half_width": self.omega_m_half_width, + "theory": { + f.name: 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`. + """ + fields = {f.name for f in dataclasses.fields(cls)} + unknown = set(overrides) - fields + if unknown: + raise ValueError( + f"unknown BlindingConfig fields {sorted(unknown)}; " + f"valid fields are {sorted(fields)}" + ) + overrides = dict(overrides) + 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: its sha256 hex 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. + """ + return hashlib.sha256(seed.encode("utf-8")).hexdigest() + + +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 ``ConcealDataVector`` does + internally. Deterministic: same ``(seed, config)`` ⇒ same hidden point, + forever — the reproducibility contract unblinding relies on. + """ + 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 the master SACC +# --------------------------------------------------------------------------- # +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): + """Master 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): + """Master 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 (coarse or fine 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 + + +# --------------------------------------------------------------------------- # +# Extract → conceal → merge, per block +# --------------------------------------------------------------------------- # +def _blindable_blocks(s): + """The blindable blocks of a master SACC as ``(name, indices, factory)``. + + ``indices`` are each block's recorded master row indices (ascending, so + the extracted sub-SACC preserves master row order); ``factory`` builds + the matching ``theory_fn`` from the extracted sub-SACC. Blocks absent + from the file are simply not listed. + """ + blocks = [] + for grid in ("coarse", "fine"): + 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 _extract_block(s, indices): + """Extract rows ``indices`` (ascending) into a sub-SACC, order preserved.""" + sub = s.copy() + sub.keep_indices(np.asarray(indices, dtype=int)) + return sub + + +def _concealing_factor(s, indices, factory, config, seed): + """The fork-computed additive concealing factor for one block. + + Extracts the block into a sub-SACC containing exactly the rows the + ``theory_fn`` spans (the fork's length guard enforces the agreement), + drives ``ConcealDataVector`` with the fiducial point, the calibrated + envelope, and the block's ``theory_fn``, and returns the factor + ``t(hidden) − t(fiducial)`` aligned to ``indices``. + """ + from smokescreen import ConcealDataVector + + sub = _extract_block(s, indices) + smoke = ConcealDataVector( + config.theory.ccl_params(), + config.shifts_dict(), + sub, + seed=seed, + theory_fn=factory(sub, config.theory), + ) + smoke.calculate_concealing_factor(factor_type="add") + concealed = smoke.apply_concealing_to_likelihood_datavec() + return np.asarray(concealed, dtype=float) - np.asarray(sub.mean, dtype=float) + + +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(master, seed, config=None, label="A", log=print): + """Return a blinded copy of a master SACC (covariance/ρτ/n(z) untouched). + + Per blindable block (coarse ξ±, fine ξ±, pseudo-Cℓ_EE): extract the + block into a sub-SACC, conceal through the fork, write the shifted + values back at their recorded master indices (row order preserved). + COSEBIs and pure-E/B are then re-derived from the blinded fine ξ± + through the pipeline estimators; provenance is stamped and any leaked + seed key stripped. + """ + config = config or BlindingConfig() + if _concealed(master): + raise ValueError("already concealed — unblind first") + _check_master(master) + + blinded = master.copy() + for name, indices, factory in _blindable_blocks(master): + factor = _concealing_factor(master, 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") + + rederive_eb_from_fine_xi(blinded, log=log) + _stamp_provenance(blinded, seed_commitment(seed), label, config.config_digest()) + return blinded + + +def unblind_sacc(blinded, seed, config=None, log=print): + """Recover the true master SACC from a blinded one + the revealed ``seed``. + + Verifies ``sha256(seed)`` and the config digest against the stamped + metadata (loud failure on either mismatch — verification precedes + subtraction), recomputes each block's shift from the seed through the + same backends, subtracts it, and re-derives COSEBIs/pure-E/B from the + unblinded fine ξ±. + """ + config = config or BlindingConfig() + if not _concealed(blinded): + raise ValueError("file is not concealed — nothing to unblind") + 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)" + ) + + hidden = hidden_params(seed, config) + fiducial = config.theory.ccl_params() + master = blinded.copy() + for name, indices, factory in _blindable_blocks(blinded): + theory = factory(_extract_block(blinded, indices), config.theory) + factor = theory(hidden) - theory(fiducial) + _set_values(master, indices, np.asarray(master.mean)[indices] - factor) + log(f"[unblind] {name}: subtracted {len(indices)} shifts") + + rederive_eb_from_fine_xi(master, log=log) + for key in ("concealed", "blind", "blind_commitment", "blind_config_digest"): + master.metadata.pop(key, None) + return master + + +def _check_master(master): + """Loud guards: the master must carry a complete, blindable ξ± content.""" + n_plus = len(master.indices(sacc_io.XI_PLUS, grid="coarse")) + n_minus = len(master.indices(sacc_io.XI_MINUS, grid="coarse")) + if not n_plus or not n_minus: + raise ValueError( + "master SACC has no coarse ξ± — cannot blind (both ξ+ and ξ− are " + f"required; found ξ+={n_plus}, ξ−={n_minus})" + ) + if n_plus != n_minus: + raise ValueError( + f"coarse ξ+ ({n_plus}) and ξ− ({n_minus}) row counts differ — " + "refusing an ambiguous coarse block" + ) + has_derived = len(master.indices(sacc_io.COSEBI_EE)) or len( + master.indices(sacc_io.PURE_TYPES["xip_E"]) + ) + if has_derived and not len(master.indices(sacc_io.XI_PLUS, grid="fine")): + raise ValueError( + "master SACC carries COSEBIs/pure-E/B but no grid='fine' ξ± — " + "the derived statistics cannot be re-derived after blinding" + ) + + +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; ``concealed``/``blind`` mark the + file. ``seed_smokescreen`` — the raw seed Smokescreen's own writer would + stamp — is popped defensively: the seed must never ride a kept file. + """ + 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 + + +# --------------------------------------------------------------------------- # +# COSEBIs / pure-E/B re-derivation through the pipeline estimators +# --------------------------------------------------------------------------- # +def rederive_cosebis(theta, xip, xim, nmodes, scale_cut=None): + """COSEBIs (Eₙ, Bₙ) from ξ± through the pipeline kernel (values only). + + Calls the same kernel :func:`sp_validation.b_modes.calculate_cosebis` + drives (``cosmo_numba``'s ``COSEBIS.cosebis_from_xipm``) directly on the + ξ± values — the covariance/χ² machinery of ``calculate_cosebis`` needs a + patched GGCorrelation we don't have and don't want (blinding never + touches uncertainty), so this is the values-only seam one layer down. + Identical inputs ⇒ identical numbers. + """ + from cosmo_numba.B_modes.cosebis import COSEBIS + + theta = np.asarray(theta) + tmin, tmax = scale_cut if scale_cut is not None else (theta.min(), theta.max()) + cosebis = COSEBIS(theta_min=tmin, theta_max=tmax, N_max=nmodes, precision=120) + En, Bn = cosebis.cosebis_from_xipm( + theta, np.asarray(xip), np.asarray(xim), parallel=True + ) + return np.asarray(En), np.asarray(Bn) + + +def _grid_edges(theta): + """The (left, right) bin edges the pipeline's ``gg`` object carries. + + TreeCorr's log-binned ``left_edges``/``right_edges`` are the + geometric-mean boundaries between adjacent bin *centres*, with the + outermost half-bins reflected. :func:`rederive_pure_eb` callers take the + Schneider-2022 integration range from ``left_edges[0]`` / + ``right_edges[-1]`` (EDGES, not centres) — the pipeline's own + ``calculate_pure_eb_correlation`` convention. + """ + theta = np.asarray(theta) + mid = np.sqrt(theta[:-1] * theta[1:]) + left = np.concatenate([[theta[0] ** 2 / mid[0]], mid]) + right = np.concatenate([mid, [theta[-1] ** 2 / mid[-1]]]) + return left, right + + +def rederive_pure_eb( + theta_report, xip_report, xim_report, theta_int, xip_int, xim_int, tmin, tmax +): + """Pure-E/B correlation functions from ξ± through the pipeline kernel. + + Calls the kernel :func:`sp_validation.b_modes.calculate_pure_eb_correlation` + drives (``cosmo_numba``'s ``get_pure_EB_modes``) directly on the ξ± + values. The reporting grid must be a strict sub-range of the integration + grid; ``tmin``/``tmax`` are the *reporting-grid edges* (see + :func:`_grid_edges`) — the pipeline's convention. A reporting point + coinciding with the integration boundary is degenerate (no interior + support) and comes back NaN, exactly as the pipeline's own estimator + returns it — never a spurious finite value. + + Returns + ------- + dict + Keyed by ``sacc_io.PURE_KEYS`` (xip_E, xim_E, xip_B, xim_B, xip_amb, + xim_amb). + """ + from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes + + modes = get_pure_EB_modes( + theta=np.asarray(theta_report), + xip=np.asarray(xip_report), + xim=np.asarray(xim_report), + theta_int=np.asarray(theta_int), + xip_int=np.asarray(xip_int), + xim_int=np.asarray(xim_int), + tmin=tmin, + tmax=tmax, + parallel=True, + ) + return dict(zip(sacc_io.PURE_KEYS, (np.asarray(m) for m in modes))) + + +def _blocks_pairs(s, dtype): + """Source pairs ``(i ≤ j)`` carrying ``dtype`` points in ``s``.""" + bins = source_bins(s) + return [ + (i, j) + for a, i in enumerate(bins) + for j in bins[a:] + if len(s.indices(dtype, sacc_io._pair((i, j)))) + ] + + +def rederive_eb_from_fine_xi(master, log=print): + """Replace COSEBIs and pure-E/B in ``master`` from its fine ξ±, in place. + + For each COSEBI / pure-E/B block present, re-derive its values through + the pipeline estimators on the master's own ``grid='fine'`` ξ± of the + matching source pair and overwrite them. Shared by blind and unblind + (the fine ξ± differs; the re-derivation is identical). Only VALUES + change; the covariance blocks stay untouched. + """ + for i, j in _blocks_pairs(master, sacc_io.COSEBI_EE): + theta, xip, xim = sacc_io.get_xi(master, (i, j), grid="fine") + tr = sacc_io._pair((i, j)) + idx_e = master.indices(sacc_io.COSEBI_EE, tr) + scale_cut = ( + master.data[idx_e[0]].tags["theta_min"], + master.data[idx_e[0]].tags["theta_max"], + ) + En, Bn = rederive_cosebis(theta, xip, xim, len(idx_e), scale_cut=scale_cut) + _set_values(master, idx_e, En) + _set_values(master, master.indices(sacc_io.COSEBI_BB, tr), Bn) + log(f"[blind] COSEBIs {tr}: re-derived {len(En)} Eₙ + {len(Bn)} Bₙ") + + for i, j in _blocks_pairs(master, sacc_io.PURE_TYPES["xip_E"]): + theta_int, xip_int, xim_int = sacc_io.get_xi(master, (i, j), grid="fine") + tr = sacc_io._pair((i, j)) + theta_report, _ = sacc_io.get_pure_eb(master, (i, j)) + # The pipeline reports pure-E/B on a coarser grid while integrating + # over the fine grid; the stored block carries the reporting θ. + # Sample the fine ξ± at the reporting θ and integrate over the full + # fine grid, mirroring calculate_pure_eb_correlation's two-grid call. + xip_report = np.interp(theta_report, theta_int, xip_int) + xim_report = np.interp(theta_report, theta_int, xim_int) + left, right = _grid_edges(theta_report) + modes = rederive_pure_eb( + theta_report, + xip_report, + xim_report, + theta_int, + xip_int, + xim_int, + float(left[0]), + float(right[-1]), + ) + for key in sacc_io.PURE_KEYS: + _set_values(master, master.indices(sacc_io.PURE_TYPES[key], tr), modes[key]) + log(f"[blind] pure-E/B {tr}: re-derived 6 blocks") + + +# --------------------------------------------------------------------------- # +# File-level custody: blind / unblind with commitment + encrypted bundle +# --------------------------------------------------------------------------- # +def blind_file(master_path, out_dir, config=None, label="A", log=print): + """Blind a master SACC file end-to-end under hash-commitment custody. + + 1. Draw an OS-entropy seed. + 2. Blind immediately (all blindable blocks; :func:`blind_sacc`). + 3. Write ``blind_{label}_commitment.json`` (repo-committable): + ``sha256(seed)`` + the canonical config digest. + 4. Encrypt the seed + the true data vector into a Fernet bundle + (``smokescreen.encryption``); the plaintext bundle is deleted. + 5. Write the blinded master. No plaintext seed and no true vector + survive on disk. + + Returns + ------- + dict + Paths of everything written: ``blinded``, ``commitment``, + ``bundle``, ``key``. + """ + config = config or BlindingConfig() + master = sacc_io.load(master_path) + paths = _output_paths(master_path, out_dir, label) + 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; choose another label or directory" + ) + + seed = secrets.token_hex(16) + blinded = blind_sacc(master, seed, config=config, label=label, log=log) + + commitment = { + "label": label, + "seed_sha256": seed_commitment(seed), + "config_digest": config.config_digest(), + } + with open(paths["commitment"], "w", encoding="utf-8") as f: + json.dump(commitment, f, indent=2, sort_keys=True) + + _write_encrypted_bundle(paths, seed, np.asarray(master.mean, dtype=float), label) + sacc_io.save(blinded, paths["blinded"]) + log(f"[blind] wrote {paths['blinded']}") + log(f"[blind] commitment (repo-committable): {paths['commitment']}") + log(f"[blind] encrypted bundle + key: {paths['bundle']}, {paths['key']}") + return paths + + +def unblind_file( + blinded_path, + bundle_path, + key_path, + commitment_path, + out_path, + config=None, + log=print, +): + """Unblind a blinded master SACC file, verifying the commitment first. + + Decrypt the bundle → verify ``sha256(seed)`` and the config digest + against ``commitment.json`` (fail closed on either) → recompute the + shifts from the seed and subtract (:func:`unblind_sacc`) → cross-check + the subtraction against the bundle's stored true vector → restore the + true vector bit-for-bit. + """ + config = config or BlindingConfig() + bundle = _read_encrypted_bundle(bundle_path, key_path) + with open(commitment_path, encoding="utf-8") as f: + commitment = json.load(f) + if seed_commitment(bundle["seed"]) != commitment["seed_sha256"]: + raise ValueError( + "bundle seed does not match the committed sha256(seed) — refusing " + "to unblind" + ) + if config.config_digest() != commitment["config_digest"]: + raise ValueError( + "blinding config does not match the committed config digest — " + "refusing to unblind (a wrong envelope or P(k) recipe would " + "silently subtract a wrong shift)" + ) + + blinded = sacc_io.load(blinded_path) + master = unblind_sacc(blinded, bundle["seed"], config=config, log=log) + + # The subtraction must agree with the escrowed truth; then the truth is + # restored exactly (float add-then-subtract leaves ~ulp residue, and the + # re-derived COSEBIs inherit it — the bundle removes it). + true_mean = np.asarray(bundle["true_mean"], dtype=float) + recovered = np.asarray(master.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 bundle for this file?" + ) + _set_values(master, np.arange(len(true_mean)), true_mean) + sacc_io.save(master, out_path) + log(f"[unblind] verified (subtraction residual {residual:.2e}); wrote {out_path}") + return out_path + + +def _output_paths(master_path, out_dir, label): + """Label-scoped output paths for one blind.""" + stem, ext = os.path.splitext(os.path.basename(master_path)) + return { + "blinded": os.path.join(out_dir, f"{stem}_blinded_{label}{ext or '.fits'}"), + "commitment": os.path.join(out_dir, f"blind_{label}_commitment.json"), + "bundle": os.path.join(out_dir, f"blind_{label}_bundle.encrpt"), + "key": os.path.join(out_dir, f"blind_{label}_bundle.key"), + } + + +def _write_encrypted_bundle(paths, seed, true_mean, label): + """Encrypt ``{seed, true_mean}`` to ``paths['bundle']``; delete plaintext. + + ``smokescreen.encryption.encrypt_file`` (Fernet) writes + ``.encrpt`` + ``.key`` side by side and removes the + plaintext (``keep_original=False``). + """ + from smokescreen.encryption import encrypt_file + + plaintext = paths["bundle"].replace(".encrpt", ".json") + with open(plaintext, "w", encoding="utf-8") as f: + json.dump({"label": label, "seed": seed, "true_mean": true_mean.tolist()}, f) + encrypt_file(plaintext, save_file=True, keep_original=False) + + +def _read_encrypted_bundle(bundle_path, key_path): + """Decrypt and parse the seed/true-vector bundle.""" + from smokescreen.encryption import decrypt_file + + return json.loads(decrypt_file(bundle_path, key_path).decode("utf-8")) From 46db6cd283666937274d73e655800b900fff77be Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 00:42:15 +0200 Subject: [PATCH 03/17] feat(blinding): blind / unblind / verify CLI for the master SACC MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit One-shot custody semantics: refuses to overwrite any existing blind output before touching theory code; the seed is drawn from OS entropy and never printed. verify is a seedless metadata↔commitment check. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- scripts/blind_data_vector.py | 141 +++++++++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) create mode 100644 scripts/blind_data_vector.py diff --git a/scripts/blind_data_vector.py b/scripts/blind_data_vector.py new file mode 100644 index 00000000..7defc657 --- /dev/null +++ b/scripts/blind_data_vector.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python3 + +"""Script blind_data_vector.py + +Blind, unblind, or verify a master SACC data-vector file with +:mod:`sp_validation.blinding` (Smokescreen-fork concealment, hash-commitment +custody). + +``blind`` draws a seed from OS entropy, applies the concealing shift +*immediately* (the seed is never printed or logged), and writes the blinded +master plus a public commitment JSON (safe to commit to the repo) and an +encrypted seed+truth bundle; no plaintext seed or true vector survives on +disk. ``unblind`` verifies the commitment and restores the true vector. +``verify`` is a cheap, seedless check that a blinded file matches a +commitment. + +:Authors: Cail Daley + +Examples +-------- +Blind:: + + blind_data_vector.py blind v1.4.6.3.sacc -o blinded/ --label A + +Unblind:: + + blind_data_vector.py unblind blinded/v1.4.6.3_blinded_A.sacc \\ + --bundle blinded/blind_A_bundle.encrpt \\ + --key blinded/blind_A_bundle.key \\ + --commitment blinded/blind_A_commitment.json \\ + -o unblinded/v1.4.6.3.sacc + +Verify:: + + blind_data_vector.py verify blinded/v1.4.6.3_blinded_A.sacc \\ + blinded/blind_A_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(args): + config = _config_from_args(args) + out_dir = pathlib.Path(args.output_dir) + out_dir.mkdir(parents=True, exist_ok=True) + # Refuse before any heavy compute: a blind is a one-shot custody event + # and silently overwriting a previous blind's outputs would destroy the + # record tying that blind to its seed. + paths = blinding._output_paths(args.master, str(out_dir), args.label) + clashes = [p for p in paths.values() if pathlib.Path(p).exists()] + if clashes: + raise SystemExit( + "refusing to overwrite existing blind output(s):\n " + + "\n ".join(clashes) + + "\nPick a fresh --output-dir or --label (never overwrite a blind)." + ) + blinding.blind_file(args.master, str(out_dir), config=config, label=args.label) + print("Commit the commitment JSON to the repo; keep the bundle + key safe.") + + +def _unblind(args): + config = _config_from_args(args) + blinding.unblind_file( + args.blinded, + args.bundle, + args.key, + args.commitment, + args.output, + config=config, + ) + + +def _verify(args): + """Seedless check: blinded-file metadata ↔ commitment JSON.""" + s = sacc_io.load(args.blinded) + 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_sha256"]: + problems.append("blind_commitment does not match the committed sha256(seed)") + if s.metadata.get("blind_config_digest") != commitment["config_digest"]: + problems.append("blind_config_digest does not match the committed digest") + 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", "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": + p.add_argument("master", help="master SACC file to blind") + p.add_argument("-o", "--output-dir", required=True) + p.add_argument("--label", default="A", help="blind label (default A)") + p.set_defaults(func=_blind) + else: + p.add_argument("blinded", help="blinded master SACC file") + p.add_argument("--bundle", required=True, help="encrypted seed bundle") + p.add_argument("--key", required=True, help="bundle Fernet key file") + p.add_argument("--commitment", required=True, help="commitment JSON") + 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 master 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()) From 11b88952b7a1d5d8025973f32b919fe7d1fdbb88 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 00:50:19 +0200 Subject: [PATCH 04/17] =?UTF-8?q?test(blinding):=20acceptance=20criteria?= =?UTF-8?q?=20AC1=E2=80=93AC9=20+=20CAMB=E2=86=94CCL=20cross-check=20AC10?= =?UTF-8?q?=E2=80=9314?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit test_blinding.py: fast coverage of the envelope calibration, canonical digest, commitment, fork-draw determinism/purity, merge alignment and guards (monkeypatched factor, no CCL); slow coverage of AC2 (on-file shift == theory difference, hidden recovered via the fork's draw, all three blocks, two-bin), AC3 (independently written direct-CCL reference), AC7 (seed reproducibility); slow+cosmo_numba coverage of AC1 (zero-shift identity through the estimators), AC4 (ΔBₙ independent of injected B — COSEBIs exact to ~1e-10; the pure-ξ_B leakage bounded at the adaptive- quadrature floor), AC5 (covariance/ρ/τ/row order untouched), AC9 (boundary-degenerate pure-EB point NaN in both re-derivations), and AC6+AC8 (end-to-end custody: no plaintext seed on disk, fail-closed on tampered commitment or wrong config, bit-for-bit restoration). test_camb_ccl_crosscheck.py: written fresh against sp_validation.cosmology — AC10 σ8/A_s reconciliation (nominal offset 3.4%; closed-form rescale lands to 1e-16), AC11/AC12 ξ± agreement at and off the fiducial (observed floor ξ+ 0.24%/0.28%, ξ− 0.10%/0.12% against 0.5%/1.0% tolerances), AC13 halofit token pinned to the CosmoSIS inference ini, AC14 fast smoke. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/tests/test_blinding.py | 797 ++++++++++++++++++ .../tests/test_camb_ccl_crosscheck.py | 170 ++++ 2 files changed, 967 insertions(+) create mode 100644 src/sp_validation/tests/test_blinding.py create mode 100644 src/sp_validation/tests/test_camb_ccl_crosscheck.py diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py new file mode 100644 index 00000000..8dadf804 --- /dev/null +++ b/src/sp_validation/tests/test_blinding.py @@ -0,0 +1,797 @@ +"""Tests for :mod:`sp_validation.blinding` — Smokescreen-fork blinding wiring. + +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, and the merge +alignment against a monkeypatched concealing factor) run in the default +suite; the theory tests (fork + CCL) are marked ``slow``; the COSEBIs / +pure-E/B tests additionally ``importorskip`` ``cosmo_numba``. + +All fixtures are synthetic and deterministic. The fork's ``ConcealDataVector`` +carries no data-vector consistency check (only the length guard), 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.cosmology import TheoryConfig + +_NOLOG = lambda *a, **k: None # noqa: E731 + + +# --------------------------------------------------------------------------- # +# Synthetic 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 _coarse_theta(n=8): + return np.geomspace(5.0, 250.0, n) + + +def _fine_theta(n=80): + # The fine 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 _report_theta(n=12): + """The pure-E/B REPORTING grid — a strict sub-range of the fine grid.""" + return np.geomspace(5.0, 100.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 make_master( + nbins=1, + with_cl=True, + with_fine=True, + with_cosebis=False, + with_pure_eb=False, + with_rho=True, + b_amplitude=0.0, + report_theta=None, + nmodes=6, +): + """Build a master SACC in the sacc_io layout, synthetic values throughout. + + COSEBIs / pure-E/B blocks (when requested) are seeded with the *pipeline + estimators run on the fixture's own fine ξ±* — so a zero-shift blind must + reproduce them exactly (AC1) — and therefore require ``cosmo_numba``. + ``b_amplitude`` injects a pure B-mode into the fine ξ± (+ξ_B to ξ+, −ξ_B + to ξ−; b_modes.py sign convention). + """ + nz = {i: _gauss_nz(0.5 + 0.3 * i, 0.15 + 0.02 * i) for i in range(nbins)} + ctheta, ftheta = _coarse_theta(), _fine_theta() + rtheta = _report_theta() if report_theta is None else report_theta + pairs = [(i, j) for i in range(nbins) for j in range(i, nbins)] + + s = sio.new_sacc(nz, metadata={"catalogue_version": "vTEST"}) + blocks = [] + + fine_xi = {} + for k, (i, j) in enumerate(pairs): + xip, xim = _xi_template(ftheta, k) + xi_b = _b_mode_template(ftheta, b_amplitude) + fine_xi[(i, j)] = (xip + xi_b, xim - xi_b) + + for k, (i, j) in enumerate(pairs): + xip, xim = _xi_template(ctheta, k) + sio.add_xi(s, (i, j), ctheta, xip, xim, grid="coarse") + tr = sio._pair((i, j)) + idx = np.concatenate( + [ + s.indices(sio.XI_PLUS, tr, grid="coarse"), + s.indices(sio.XI_MINUS, tr, grid="coarse"), + ] + ) + blocks.append((idx, np.eye(len(idx)) * 1e-12)) + + if with_cl: + 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() + for k, (i, j) in enumerate(pairs): + 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)) + + if with_cosebis: + for i, j in pairs: + xip_f, xim_f = fine_xi[(i, j)] + En, Bn = bd.rederive_cosebis( + ftheta, xip_f, xim_f, nmodes, scale_cut=(ftheta.min(), ftheta.max()) + ) + sio.add_cosebis(s, (i, j), En, Bn, (ftheta.min(), ftheta.max())) + tr = sio._pair((i, j)) + idx = np.concatenate( + [s.indices(sio.COSEBI_EE, tr), s.indices(sio.COSEBI_BB, tr)] + ) + blocks.append((idx, np.eye(len(idx)) * 1e-20)) + + if with_pure_eb: + for i, j in pairs: + xip_f, xim_f = fine_xi[(i, j)] + xip_r = np.interp(rtheta, ftheta, xip_f) + xim_r = np.interp(rtheta, ftheta, xim_f) + left, right = bd._grid_edges(rtheta) + modes = bd.rederive_pure_eb( + rtheta, + xip_r, + xim_r, + ftheta, + xip_f, + xim_f, + float(left[0]), + float(right[-1]), + ) + sio.add_pure_eb(s, (i, j), rtheta, *(modes[k] for k in sio.PURE_KEYS)) + tr = sio._pair((i, j)) + idx = np.concatenate( + [s.indices(sio.PURE_TYPES[k], tr) for k in sio.PURE_KEYS] + ) + blocks.append((idx, np.eye(len(idx)) * 1e-14)) + + if with_rho: + 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)) + + if with_fine: + for i, j in pairs: + xip_f, xim_f = fine_xi[(i, j)] + sio.add_xi(s, (i, j), ftheta, xip_f, xim_f, grid="fine") + tr = sio._pair((i, j)) + idx = np.concatenate( + [ + s.indices(sio.XI_PLUS, tr, grid="fine"), + s.indices(sio.XI_MINUS, tr, grid="fine"), + ] + ) + blocks.append((idx, np.eye(len(idx)) * 1e-14)) + + sio.assemble_covariance(s, blocks) + return s + + +# --------------------------------------------------------------------------- # +# 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() + ) + + +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_sha256_of_seed(): + import hashlib + + seed = "the-secret" + assert bd.seed_commitment(seed) == hashlib.sha256(seed.encode()).hexdigest() + assert bd.seed_commitment("right") != bd.seed_commitment("wrong") + + +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]) + + +# --------------------------------------------------------------------------- # +# Merge + provenance with a monkeypatched factor (fast — no CCL) +# --------------------------------------------------------------------------- # +def _patch_constant_factor(monkeypatch, value=1e-6): + def fake(master, 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): the shift lands exactly on the blindable rows, + in master order; ρ/τ, covariance, n(z), and every tag are untouched.""" + _patch_constant_factor(monkeypatch) + s = make_master(nbins=2, with_cl=True, with_fine=True, with_rho=True) + orig = np.array(s.mean) + orig_cov = s.covariance.dense.copy() + orig_nz = s.tracers["source_0"].nz.copy() + + blinded = bd.blind_sacc(s, "seed", log=_NOLOG) + + shifted = np.zeros(len(orig), dtype=bool) + for _, indices, _ in bd._blindable_blocks(s): + 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(s.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 files carry concealed/blind/commitment/digest; no seed key.""" + _patch_constant_factor(monkeypatch) + s = make_master(with_cl=False) + 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 "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_master(with_cl=False) + blinded = bd.blind_sacc(s, "seed", log=_NOLOG) + with pytest.raises(ValueError, match="already concealed"): + bd.blind_sacc(blinded, "seed2", log=_NOLOG) + + +def test_guard_requires_coarse_xi(): + nz = {0: _gauss_nz(0.7, 0.2)} + s = sio.new_sacc(nz) + with pytest.raises(ValueError, match="no coarse"): + bd.blind_sacc(s, "seed", log=_NOLOG) + + +def test_guard_requires_both_xi_plus_and_minus(): + nz = {0: _gauss_nz(0.7, 0.2)} + s = sio.new_sacc(nz) + for th in _coarse_theta(): + s.add_data_point( + sio.XI_PLUS, + ("source_0", "source_0"), + 1e-5, + theta=float(th), + grid="coarse", + ) + with pytest.raises(ValueError, match="no coarse"): + bd.blind_sacc(s, "seed", log=_NOLOG) + + +def test_guard_derived_blocks_require_fine_xi(): + """COSEBIs present but no fine ξ± ⇒ loud refusal (cannot re-derive).""" + nz = {0: _gauss_nz(0.7, 0.2)} + s = sio.new_sacc(nz) + ctheta = _coarse_theta() + xip, xim = _xi_template(ctheta) + sio.add_xi(s, (0, 0), ctheta, xip, xim, grid="coarse") + sio.add_cosebis(s, (0, 0), np.ones(4) * 1e-10, np.ones(4) * 1e-12, (1.0, 100.0)) + with pytest.raises(ValueError, match="no grid='fine'"): + bd.blind_sacc(s, "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_master(with_cl=False) + 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) + + +def test_cli_blind_refuses_existing_output(tmp_path): + """The blind CLI refuses to overwrite a previous blind (before compute).""" + 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 / "master_blinded_A.sacc").write_text("stale") + with pytest.raises(SystemExit, match="refusing to overwrite"): + cli.main(["blind", "master.sacc", "-o", str(tmp_path), "--label", "A"]) + + +# --------------------------------------------------------------------------- # +# AC2 + AC3 + AC7: the shift itself (slow — fork + CCL) +# --------------------------------------------------------------------------- # +@pytest.mark.slow +def test_ac2_on_file_shift_equals_theory_difference(): + """AC2: per-row shift on the blinded file == theory_fn(hidden) − + theory_fn(fiducial), hidden recovered by re-running the fork's draw — + for all three blocks, on a two-bin fixture.""" + master = make_master(nbins=2, with_cl=True, with_fine=True) + cfg = bd.BlindingConfig() + seed = "ac2-seed" + blinded = bd.blind_sacc(master, seed, config=cfg, log=_NOLOG) + + hidden = bd.hidden_params(seed, cfg) + fiducial = cfg.theory.ccl_params() + worst = 0.0 + blocks = bd._blindable_blocks(master) + assert {name for name, _, _ in blocks} == {"coarse ξ±", "fine ξ±", "pseudo-Cℓ_EE"} + for name, indices, factory in blocks: + theory = factory(bd._extract_block(master, indices), cfg.theory) + expected = theory(hidden) - theory(fiducial) + actual = np.array(blinded.mean)[indices] - np.array(master.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}: |Δ|={gap:.3e}" + print(f"\nAC2 max relative shift mismatch across blocks: {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 + + master = make_master(nbins=2, with_cl=True, with_fine=False) + cfg = bd.BlindingConfig() + seed = "ac3-seed" + blinded = bd.blind_sacc(master, 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(master, "coarse"): + tr = sio._pair((i, j)) + theta = sio._tag(master, sio.XI_PLUS, tr, "theta", grid="coarse") + nz_i, nz_j = sio.get_nz(master, i), sio.get_nz(master, 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 = master.indices(dt, tr, grid="coarse") + realized = np.array(blinded.mean)[idx] - np.array(master.mean)[idx] + worst = max(worst, np.max(np.abs(realized - ref_shift))) + print(f"\nAC3 max |realized − independent reference| (coarse ξ±): {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(master): + tr = sio._pair((i, j)) + idx = master.indices(sio.CL_EE, tr) + window = master.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(master, i), sio.get_nz(master, 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.mean)[idx] - np.array(master.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 with the same (seed, config) produce identical + shifts; a different seed produces a different one.""" + master = make_master(nbins=1, with_cl=True, with_fine=True) + b1 = bd.blind_sacc(master, "repro-seed", log=_NOLOG) + b2 = bd.blind_sacc(master, "repro-seed", log=_NOLOG) + assert np.array_equal(np.array(b1.mean), np.array(b2.mean)) + b3 = bd.blind_sacc(master, "other-seed", log=_NOLOG) + assert not np.array_equal(np.array(b1.mean), np.array(b3.mean)) + + +# --------------------------------------------------------------------------- # +# AC1, AC4, AC5, AC9: derived statistics (slow + cosmo_numba) +# --------------------------------------------------------------------------- # +@pytest.mark.slow +def test_ac1_zero_shift_is_identity(): + """AC1: a zero envelope reproduces ξ±, COSEBIs, and pure-E/B exactly — + the extract/merge/re-derivation plumbing is the identity at zero shift.""" + pytest.importorskip("cosmo_numba") + master = make_master( + with_cl=True, with_fine=True, with_cosebis=True, with_pure_eb=True + ) + zero = bd.BlindingConfig(s8_half_width=0.0, omega_m_half_width=0.0) + blinded = bd.blind_sacc(master, "any-seed", config=zero, log=_NOLOG) + orig, new = np.array(master.mean), np.array(blinded.mean) + both_nan = np.isnan(orig) & np.isnan(new) # degenerate pure-EB boundary pts + assert np.array_equal(orig[~both_nan], new[~both_nan]), ( + "zero-shift blind changed values" + ) + assert np.array_equal(np.isnan(orig), np.isnan(new)) + + +@pytest.mark.slow +def test_ac4_b_mode_invariance_and_leakage_floor(): + """AC4: the ΔBₙ induced by blinding 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): + master = make_master( + with_cl=False, + with_fine=True, + with_cosebis=True, + with_pure_eb=True, + b_amplitude=amp, + ) + blinded = bd.blind_sacc(master, seed, log=_NOLOG) + tr = sio._pair((0, 0)) + idx_b = master.indices(sio.COSEBI_BB, tr) + d_bn = np.array(blinded.mean)[idx_b] - np.array(master.mean)[idx_b] + bn0 = np.array(master.mean)[idx_b] + # pure-EB B blocks: fixed absolute leakage too + idx_pb = master.indices(sio.PURE_TYPES["xip_B"], tr) + d_xib = np.array(blinded.mean)[idx_pb] - np.array(master.mean)[idx_pb] + 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(bn0)):.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 + # identical to ~1e-10 relative (observed 1.1e-10 — quadrature float noise + # on the injected-B term; the bound carries a decade of margin) + 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)| = {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 ~1e-2 of the (tiny) leakage itself — observed + # 2.0e-11 absolute on a 1.5e-9 leakage, i.e. 1e-5 of the smaller + # injected B. 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} ({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: covariance and ρ/τ byte-identical; every shifted block's rows land + at their original master indices (order-preservation).""" + pytest.importorskip("cosmo_numba") + master = make_master( + with_cl=True, + with_fine=True, + with_cosebis=True, + with_pure_eb=True, + with_rho=True, + ) + blinded = bd.blind_sacc(master, "ac5-seed", log=_NOLOG) + + assert np.array_equal(blinded.covariance.dense, master.covariance.dense) + for dt in ( + sio.RHO_PLUS.format(k=0), + sio.RHO_MINUS.format(k=0), + sio.RHO_PLUS.format(k=1), + sio.RHO_MINUS.format(k=1), + sio.TAU_PLUS.format(k=0), + sio.TAU_MINUS.format(k=0), + sio.CL_BB, + sio.CL_EB, + ): + idx = master.indices(dt) + assert np.array_equal( + np.array(blinded.mean)[idx], np.array(master.mean)[idx] + ), f"{dt} was touched by the blind" + # row order: the identity of every row (type/tracers/tags) is unchanged + for a, b in zip(master.data, blinded.data): + assert (a.data_type, a.tracers, a.tags) == (b.data_type, b.tracers, b.tags) + # ξ± rows did move (the blind actually blinded) + xi_idx = bd._xi_indices(master, "coarse") + assert not np.allclose( + np.array(blinded.mean)[xi_idx], np.array(master.mean)[xi_idx], atol=0 + ) + + +@pytest.mark.slow +def test_ac9_pure_eb_boundary_nan_parity(): + """AC9: the boundary-degenerate reporting point — the point nearest the + edge-based integration bound, where the Schneider estimator has no + interior support (observed: the LAST reporting point, at the tmax edge) + — is NaN in both the true-file and blinded-file re-derivations. + NaN-parity with the pipeline's own estimator, never a spurious finite.""" + pytest.importorskip("cosmo_numba") + master = make_master(with_cl=False, with_fine=True, with_pure_eb=True) + tr = sio._pair((0, 0)) + idx_e = master.indices(sio.PURE_TYPES["xip_E"], tr) + true_vals = np.array(master.mean)[idx_e] + assert np.isnan(true_vals[-1]), ( + "fixture did not trigger the boundary degeneracy — the outermost " + "reporting point should be NaN in the pipeline estimator" + ) + assert np.all(np.isfinite(true_vals[:-1])) + + blinded = bd.blind_sacc(master, "ac9-seed", log=_NOLOG) + blind_vals = np.array(blinded.mean)[idx_e] + assert np.array_equal(np.isnan(true_vals), np.isnan(blind_vals)), ( + "blinding moved the NaN pattern of the degenerate boundary point" + ) + assert np.all(np.isfinite(blind_vals[:-1])) + + +# --------------------------------------------------------------------------- # +# AC6 + AC8: end-to-end custody through the CLIs (slow + cosmo_numba) +# --------------------------------------------------------------------------- # +@pytest.mark.slow +def test_ac6_ac8_end_to_end_blind_unblind_custody(tmp_path): + """AC8: the blind CLI blinds a master file end-to-end and the unblind CLI + restores it bit-for-bit. AC6: no plaintext seed anywhere on disk, no + ``seed_smokescreen`` key; unblind fails closed on a tampered commitment.""" + pytest.importorskip("cosmo_numba") + master = make_master( + with_cl=True, + with_fine=True, + with_cosebis=True, + with_pure_eb=True, + with_rho=True, + ) + master_path = tmp_path / "vTEST.fits" + sio.save(master, str(master_path)) + + out = tmp_path / "blinded" + out.mkdir() + paths = bd.blind_file(str(master_path), str(out), label="A", log=_NOLOG) + + # -- custody hygiene (AC6) ------------------------------------------- -- + blinded = sio.load(paths["blinded"]) + assert blinded.metadata["concealed"] is True + assert "seed_smokescreen" not in blinded.metadata + with open(paths["commitment"], encoding="utf-8") as f: + commitment = json.load(f) + assert blinded.metadata["blind_commitment"] == commitment["seed_sha256"] + assert not (out / "blind_A_bundle.json").exists() # plaintext deleted + # only the four custody outputs exist; the commitment carries hashes only + assert {p.name for p in out.iterdir()} == { + pathlib.Path(v).name for v in paths.values() + } + assert set(commitment) == {"label", "seed_sha256", "config_digest"} + # the blinded vector actually differs from the truth + assert not np.array_equal(np.array(blinded.mean), np.array(master.mean)) + + # -- fail-closed on tampered commitment (AC6) ------------------------- -- + tampered = dict(commitment, seed_sha256="0" * 64) + tampered_path = tmp_path / "tampered.json" + tampered_path.write_text(json.dumps(tampered)) + with pytest.raises(ValueError, match="sha256"): + bd.unblind_file( + paths["blinded"], + paths["bundle"], + paths["key"], + str(tampered_path), + str(tmp_path / "never.fits"), + log=_NOLOG, + ) + with pytest.raises(ValueError, match="config digest"): + bd.unblind_file( + paths["blinded"], + paths["bundle"], + paths["key"], + paths["commitment"], + str(tmp_path / "never.fits"), + config=bd.BlindingConfig(s8_half_width=0.01), + log=_NOLOG, + ) + + # -- bit-for-bit restoration (AC8) ------------------------------------ -- + restored_path = tmp_path / "restored.fits" + bd.unblind_file( + paths["blinded"], + paths["bundle"], + paths["key"], + paths["commitment"], + str(restored_path), + log=_NOLOG, + ) + restored = sio.load(str(restored_path)) + orig, back = np.array(master.mean), np.array(restored.mean) + both_nan = np.isnan(orig) & np.isnan(back) + assert np.array_equal(orig[~both_nan], back[~both_nan]) # bit-for-bit + assert np.array_equal(np.isnan(orig), np.isnan(back)) + assert not restored.metadata.get("concealed", False) + assert "blind_commitment" not in restored.metadata 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..33ec64ff --- /dev/null +++ b/src/sp_validation/tests/test_camb_ccl_crosscheck.py @@ -0,0 +1,170 @@ +"""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.cosmology`: + +- **Path A** (:func:`~sp_validation.cosmology.xi_ccl`): CCL-native — CCL's + Boltzmann-CAMB HMCode2020 P(k) route, projected by CCL Limber + FFTLog. +- **Path B** (:func:`~sp_validation.cosmology.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 cosmology 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" + / "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 + + +# --------------------------------------------------------------------------- # +# 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" From dfd1cbb3c45e0457ea3cc19c7891fe0d00fa233e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 12:14:52 +0200 Subject: [PATCH 05/17] fix(blinding): pin the Smokescreen fork + cosmo_numba; stamp pure-EB/COSEBI writer contracts MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The `blinding` extra pinned `smokescreen==1.5.6`, which resolves to the upstream LSSTDESC PyPI package — the fork self-reports the same version, so the version pin silently installed a package without `param_shifts. draw_param_shifts` or `ConcealDataVector`'s `fiducial_params`/`theory_fn` signature, and CI failed at import. Pin the UNIONS-WL fork by commit, and add cosmo_numba (Cail's numpy-2.x fork commit) as a declared dependency so the COSEBIs / pure-E/B re-derivation kernels are present in the image. sacc_io: make the derived-statistics writer contracts explicit so the blind re-derives on the SAME support the pipeline used, not a private convention. `add_cosebis`' scale_cut tags are the retained-bin support (min/max of the cut meanr, per b_modes.scale_cut_to_bins), not the requested cut. `add_pure_eb` now requires `bounds=(tmin, tmax)` — the pipeline's TreeCorr edge-based integration range (left_edges[0]/right_edges[-1]), which is not reconstructible from meanr — stored as tmin/tmax tags; and documents that the reporting grid is the coarse ξ± meanr (the measured reporting-binned ξ±), which blinding reconstructs from the stored coarse block rather than interpolating the fine grid. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- pyproject.toml | 26 +++++++++++++++++++------- src/sp_validation/sacc_io.py | 30 ++++++++++++++++++++++++++---- 2 files changed, 45 insertions(+), 11 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index ecdec6f1..9b2dd825 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -113,12 +113,23 @@ glass = [ "glass.ext.camb==2023.6", "cosmology==2022.10.9", ] -# 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). +# Data-vector blinding (PRD #241 §3-§5): the UNIONS-WL Smokescreen fork applies +# the Muir et al. shift d → d + t(hidden) − t(fid), with 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 + firecrown +# v1.15 both set the python floor (>=3.12). +# +# smokescreen MUST be the UNIONS-WL fork, pinned by commit: blinding.py uses +# fork-only surface (`param_shifts.draw_param_shifts`, `ConcealDataVector`'s +# `fiducial_params`/`theory_fn` signature). The fork self-reports the upstream +# version number (1.5.6), so a `smokescreen==1.5.6` pin silently resolves to +# the incompatible LSSTDESC PyPI package — never pin by version. +# +# cosmo_numba supplies the COSEBIs / pure-E/B kernels the blind re-derives +# derived statistics through (same kernels b_modes drives). Pinned to the +# numpy-2.x-compatible fork commit (objmode np.fft fix; upstream aguinot/main +# breaks under the numpy 2.4 this extra requires). # # firecrown is not on PyPI and declares conda-forge-only / unused sampler # connectors as hard deps, so installing this extra requires the dependency @@ -129,7 +140,8 @@ glass = [ # 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", + "smokescreen @ git+https://github.com/UNIONS-WL/Smokescreen@588a6b9b26560bd5ba3dd5ba342f3c40152644f9", + "cosmo-numba @ git+https://github.com/cailmdaley/cosmo-numba@db452c487f3d740c9fe69a409352f11987525d65", "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 diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index 22fe5aa6..d615e810 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -241,7 +241,12 @@ def add_cosebis(s, bins, En, Bn, scale_cut): scale_cut : tuple of float ``(theta_min, theta_max)`` in arcmin, stored on every point as the ``theta_min``/``theta_max`` tags; multiple cuts coexist in one file, - told apart by these tags. + told apart by these tags. **Writer contract:** these are the COSEBI + kernel's actual support — ``min``/``max`` of the *retained* fine-grid + ``meanr`` after ``b_modes.scale_cut_to_bins`` (the pipeline builds + the kernel from the cut bin centres, per Axel's recommendation) — + not the requested cut. Blinding re-derives COSEBIs from these tags, + so a requested-cut tag would rebuild the kernel on the wrong support. """ tracers = _pair(bins) theta_min, theta_max = scale_cut @@ -257,7 +262,9 @@ def add_cosebis(s, bins, En, Bn, scale_cut): ) -def add_pure_eb(s, bins, theta, xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb): +def add_pure_eb( + s, bins, theta, xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb, *, bounds +): """Add pure E/B-mode correlation functions for one tracer pair. Six blocks are inserted in ``PURE_KEYS`` order (xip_E, xim_E, xip_B, @@ -271,12 +278,25 @@ def add_pure_eb(s, bins, theta, xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb): bins : tuple of int Source bin pair ``(i, j)``. theta : array_like - Angular separations (arcmin), shared by all six blocks. + Angular separations (arcmin), shared by all six blocks. **Writer + contract:** the reporting grid must be the coarse ξ± grid's + ``meanr`` — the pipeline feeds the *measured* reporting-binned ξ± + into the Schneider estimator, and blinding reconstructs it from the + stored coarse ξ± block. xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb : array_like The six pure E/B / ambiguous mode arrays at ``theta``. + bounds : tuple of float + ``(tmin, tmax)`` in arcmin — the Schneider-2022 integration range the + estimator was run with. **Writer contract:** the pipeline's + edge-based bounds, ``gg.left_edges[0]`` / ``gg.right_edges[-1]`` of + the reporting TreeCorr correlation + (``b_modes.calculate_pure_eb_correlation``'s convention). Stored on + every point as ``tmin``/``tmax`` tags; TreeCorr's edges are not + reconstructible from ``meanr``, so blinding fails loud without them. """ _check_ascending("theta", theta) tracers = _pair(bins) + tmin, tmax = (float(b) for b in bounds) values = { "xip_E": xip_E, "xim_E": xim_E, @@ -288,7 +308,9 @@ def add_pure_eb(s, bins, theta, xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb): for key in PURE_KEYS: dtype, arr = PURE_TYPES[key], values[key] for n, th in enumerate(theta): - s.add_data_point(dtype, tracers, float(arr[n]), theta=float(th)) + s.add_data_point( + dtype, tracers, float(arr[n]), theta=float(th), tmin=tmin, tmax=tmax + ) def add_rho(s, k, theta, rho_p, rho_m): From 4a0cb86dc46579716fd75db8d55e292a8a690fa9 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 12:15:15 +0200 Subject: [PATCH 06/17] fix(blinding): re-derive derived stats on the pipeline's own inputs; delete the plaintext master MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two confirmed adversarial-review findings on blinding.py. Re-derivation seam (findings 4/6/7). The old rederive_eb_from_fine_xi rebuilt the pure-E/B reporting ξ± by np.interp of the fine grid and reconstructed TreeCorr edges from stored θ via a geometric-mean pseudo-edge helper (_grid_edges) — a private convention that does NOT match what b_modes.calculate_pure_eb_correlation feeds the Schneider kernel (the MEASURED coarse-binned ξ± at gg.meanr, with gg.left_edges[0]/right_edges[-1] bounds), so on a real pipeline master a zero-shift blind would not reproduce the stored values. Now: the reporting ξ± are read from the master's own coarse ξ± block (fail loud if the coarse grid ≠ the stored reporting grid, since the measured reporting ξ± can't otherwise be reconstructed), the integration bounds are read from the tmin/tmax tags the writer stamps, and _grid_edges is removed. rederive_cosebis masks ξ± by the stored scale_cut (the retained-bin support) before building the COSEBI kernel, matching b_modes' theta_cut/xip_cut path rather than passing the full fine grid with tag-derived tmin/tmax. Custody (finding 1). blind_file left the input master — the full plaintext true vector — on disk indefinitely, contradicting the PRD criterion that "the plaintext true vector is deleted at blinding time". It now deletes master_path after a successful blind (the true vector survives only inside the encrypted bundle); --keep-input retains it with the custody implication documented. The module docstring's custody claim and the colocated-key caveat (bundle + Fernet key in one dir is not at-rest protection) are stated explicitly. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/blinding.py | 128 +++++++++++++++++++++++----------- 1 file changed, 87 insertions(+), 41 deletions(-) diff --git a/src/sp_validation/blinding.py b/src/sp_validation/blinding.py index 9629de0c..2f255020 100644 --- a/src/sp_validation/blinding.py +++ b/src/sp_validation/blinding.py @@ -39,8 +39,10 @@ OS-entropy seed, blinds immediately, publishes ``sha256(seed)`` plus a canonical config digest as a repo-committable commitment, encrypts the seed + true vector into a Fernet bundle (``smokescreen.encryption``) and - deletes the plaintext. Unblind verifies both hashes against the - commitment before subtracting anything, then restores the true vector. + deletes both plaintexts — the temporary bundle plaintext *and* the input + master file, so no plaintext true vector survives on disk beside the + blind. Unblind verifies both hashes against the commitment before + subtracting anything, then restores the true vector. """ import dataclasses @@ -479,36 +481,32 @@ def rederive_cosebis(theta, xip, xim, nmodes, scale_cut=None): ξ± values — the covariance/χ² machinery of ``calculate_cosebis`` needs a patched GGCorrelation we don't have and don't want (blinding never touches uncertainty), so this is the values-only seam one layer down. + + ``scale_cut`` follows the :func:`sacc_io.add_cosebis` writer contract: + ``(theta_min, theta_max)`` are min/max of the *retained* bin centres + after the pipeline's ``scale_cut_to_bins``. The cut is contiguous in an + ascending grid, so selecting ``theta_min ≤ θ ≤ theta_max`` inclusively + reproduces exactly the retained set, and the kernel is built on that + set's min/max support and fed only the retained ξ± — bit-matching + ``calculate_cosebis``'s ``theta_cut``/``xip_cut``/``xim_cut`` path. Identical inputs ⇒ identical numbers. """ from cosmo_numba.B_modes.cosebis import COSEBIS - theta = np.asarray(theta) + theta, xip, xim = (np.asarray(a) for a in (theta, xip, xim)) tmin, tmax = scale_cut if scale_cut is not None else (theta.min(), theta.max()) - cosebis = COSEBIS(theta_min=tmin, theta_max=tmax, N_max=nmodes, precision=120) - En, Bn = cosebis.cosebis_from_xipm( - theta, np.asarray(xip), np.asarray(xim), parallel=True + cut = (theta >= tmin) & (theta <= tmax) + theta_cut, xip_cut, xim_cut = theta[cut], xip[cut], xim[cut] + cosebis = COSEBIS( + theta_min=np.min(theta_cut), + theta_max=np.max(theta_cut), + N_max=nmodes, + precision=120, ) + En, Bn = cosebis.cosebis_from_xipm(theta_cut, xip_cut, xim_cut, parallel=True) return np.asarray(En), np.asarray(Bn) -def _grid_edges(theta): - """The (left, right) bin edges the pipeline's ``gg`` object carries. - - TreeCorr's log-binned ``left_edges``/``right_edges`` are the - geometric-mean boundaries between adjacent bin *centres*, with the - outermost half-bins reflected. :func:`rederive_pure_eb` callers take the - Schneider-2022 integration range from ``left_edges[0]`` / - ``right_edges[-1]`` (EDGES, not centres) — the pipeline's own - ``calculate_pure_eb_correlation`` convention. - """ - theta = np.asarray(theta) - mid = np.sqrt(theta[:-1] * theta[1:]) - left = np.concatenate([[theta[0] ** 2 / mid[0]], mid]) - right = np.concatenate([mid, [theta[-1] ** 2 / mid[-1]]]) - return left, right - - def rederive_pure_eb( theta_report, xip_report, xim_report, theta_int, xip_int, xim_int, tmin, tmax ): @@ -517,8 +515,9 @@ def rederive_pure_eb( Calls the kernel :func:`sp_validation.b_modes.calculate_pure_eb_correlation` drives (``cosmo_numba``'s ``get_pure_EB_modes``) directly on the ξ± values. The reporting grid must be a strict sub-range of the integration - grid; ``tmin``/``tmax`` are the *reporting-grid edges* (see - :func:`_grid_edges`) — the pipeline's convention. A reporting point + grid; ``tmin``/``tmax`` are the reporting correlation's TreeCorr *bin + edges* (``gg.left_edges[0]`` / ``gg.right_edges[-1]``) — the pipeline's + convention, carried on the file by ``sacc_io.add_pure_eb``. A reporting point coinciding with the integration boundary is degenerate (no interior support) and comes back NaN, exactly as the pipeline's own estimator returns it — never a spurious finite value. @@ -560,10 +559,12 @@ def rederive_eb_from_fine_xi(master, log=print): """Replace COSEBIs and pure-E/B in ``master`` from its fine ξ±, in place. For each COSEBI / pure-E/B block present, re-derive its values through - the pipeline estimators on the master's own ``grid='fine'`` ξ± of the - matching source pair and overwrite them. Shared by blind and unblind - (the fine ξ± differs; the re-derivation is identical). Only VALUES - change; the covariance blocks stay untouched. + the pipeline estimators on the master's own ξ± of the matching source + pair — COSEBIs from the ``grid='fine'`` ξ±, pure-E/B from the measured + ``grid='coarse'`` reporting ξ± plus the fine integration ξ±, exactly the + inputs ``b_modes`` hands the kernels — and overwrite them. Shared by + blind and unblind (the ξ± differ; the re-derivation is identical). Only + VALUES change; the covariance blocks stay untouched. """ for i, j in _blocks_pairs(master, sacc_io.COSEBI_EE): theta, xip, xim = sacc_io.get_xi(master, (i, j), grid="fine") @@ -582,13 +583,33 @@ def rederive_eb_from_fine_xi(master, log=print): theta_int, xip_int, xim_int = sacc_io.get_xi(master, (i, j), grid="fine") tr = sacc_io._pair((i, j)) theta_report, _ = sacc_io.get_pure_eb(master, (i, j)) - # The pipeline reports pure-E/B on a coarser grid while integrating - # over the fine grid; the stored block carries the reporting θ. - # Sample the fine ξ± at the reporting θ and integrate over the full - # fine grid, mirroring calculate_pure_eb_correlation's two-grid call. - xip_report = np.interp(theta_report, theta_int, xip_int) - xim_report = np.interp(theta_report, theta_int, xim_int) - left, right = _grid_edges(theta_report) + # The pipeline feeds get_pure_EB_modes the MEASURED reporting-binned + # ξ± (gg.xip/gg.xim at gg.meanr) alongside the fine integration grid + # — never an interpolation of the fine ξ±. The master's coarse ξ± + # block IS that measurement (and carries the blind consistently), so + # the reporting grids must coincide; anything else means the seam + # cannot reconstruct the pipeline's inputs, so fail loud. + theta_coarse, xip_report, xim_report = sacc_io.get_xi( + master, (i, j), grid="coarse" + ) + if len(theta_coarse) != len(theta_report) or not np.allclose( + theta_coarse, theta_report, rtol=1e-12, atol=0.0 + ): + raise ValueError( + f"pure-E/B reporting grid for {tr} does not match the coarse " + "ξ± grid — the measured reporting ξ± cannot be reconstructed " + "and the pure-E/B block cannot be re-derived" + ) + # Integration bounds: the pipeline's TreeCorr edge-based tmin/tmax, + # stored by the writer (sacc_io.add_pure_eb bounds contract) — edges + # are not reconstructible from meanr, so their absence is fatal. + tags = master.data[master.indices(sacc_io.PURE_TYPES["xip_E"], tr)[0]].tags + if "tmin" not in tags or "tmax" not in tags: + raise ValueError( + f"pure-E/B block for {tr} carries no tmin/tmax integration-" + "bound tags (sacc_io.add_pure_eb bounds contract) — cannot " + "re-derive with the pipeline's edge-based integration range" + ) modes = rederive_pure_eb( theta_report, xip_report, @@ -596,8 +617,8 @@ def rederive_eb_from_fine_xi(master, log=print): theta_int, xip_int, xim_int, - float(left[0]), - float(right[-1]), + float(tags["tmin"]), + float(tags["tmax"]), ) for key in sacc_io.PURE_KEYS: _set_values(master, master.indices(sacc_io.PURE_TYPES[key], tr), modes[key]) @@ -607,7 +628,9 @@ def rederive_eb_from_fine_xi(master, log=print): # --------------------------------------------------------------------------- # # File-level custody: blind / unblind with commitment + encrypted bundle # --------------------------------------------------------------------------- # -def blind_file(master_path, out_dir, config=None, label="A", log=print): +def blind_file( + master_path, out_dir, config=None, label="A", keep_input=False, log=print +): """Blind a master SACC file end-to-end under hash-commitment custody. 1. Draw an OS-entropy seed. @@ -616,8 +639,17 @@ def blind_file(master_path, out_dir, config=None, label="A", log=print): ``sha256(seed)`` + the canonical config digest. 4. Encrypt the seed + the true data vector into a Fernet bundle (``smokescreen.encryption``); the plaintext bundle is deleted. - 5. Write the blinded master. No plaintext seed and no true vector - survive on disk. + 5. Write the blinded master, then **delete the input master file** — the + plaintext true vector is escrowed (encrypted) in the bundle and must + not survive on disk beside the blind (PRD custody criterion). Pass + ``keep_input=True`` to retain it (and own the custody implication). + No plaintext seed and no plaintext true vector survive on disk. + + Custody caveat: the encrypted bundle and its Fernet key are written to the + *same* ``out_dir``. Fernet is only as protective as the key's separation — + anyone with both files can decrypt the seed and the true vector. Keep the + key out-of-band from the bundle; ``out_dir`` colocation is convenience, not + at-rest protection. Returns ------- @@ -648,9 +680,23 @@ def blind_file(master_path, out_dir, config=None, label="A", log=print): _write_encrypted_bundle(paths, seed, np.asarray(master.mean, dtype=float), label) sacc_io.save(blinded, paths["blinded"]) + + # The true vector is now escrowed (encrypted) in the bundle; the plaintext + # master must not survive on disk beside the blind (PRD custody criterion). + if not keep_input: + os.remove(master_path) + log(f"[blind] wrote {paths['blinded']}") log(f"[blind] commitment (repo-committable): {paths['commitment']}") log(f"[blind] encrypted bundle + key: {paths['bundle']}, {paths['key']}") + if keep_input: + log(f"[blind] input master RETAINED at {master_path} (keep_input=True)") + else: + log(f"[blind] deleted plaintext input master {master_path}") + log( + "[blind] custody: keep the bundle key out-of-band from the bundle " + "(colocation in the output dir is not at-rest protection)" + ) return paths From a6421a379afb563039bfe9cd05a8c9dd804a1abb Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 12:15:32 +0200 Subject: [PATCH 07/17] test(blinding): seed pure-EB through the pipeline convention; run the full suite in CI; correct AC2 framing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Test- and CI-side of four confirmed findings. Finding 4 (fixture tautology). make_master seeded the pure-E/B blocks with the blinding module's own interp + _grid_edges seam — the exact code under test — so AC1/4/9 were self-consistent by construction and blind to any divergence from the pipeline. The fixture now seeds pure-E/B through the pipeline convention the re-derivation reconstructs: reporting grid = the coarse ξ± grid, reporting ξ± = the MEASURED coarse ξ± (not an interpolation of the fine grid), integration ξ± = the fine grid, with edge-based bounds passed via the new add_pure_eb(bounds=...) contract. The outermost reporting point sits at tmax (no interior support) → NaN, the AC9 boundary case. This makes AC1/4/9 a real check of the re-derivation against the pipeline convention. Finding 3 (CI never ran the ACs). deploy-image.yml ran `pytest -m "not slow"`, so every blinding AC (all @slow, including the cosmo_numba-gated AC1/4/5/9) was deselected in CI — contradicting the PRD's "CI runs it in the freshly built image before publish". The image already carries the full stack (the smoke test imports cosmo_numba), so CI now runs the FULL suite before pushing. Finding 2 (AC2 framing). AC2 runs the identical theory_fn on both sides, so a wrong cosmology cancels and a wrong backend still passes — it verifies fork-draw recovery + merge placement, not backend correctness (that is AC3, the independent CCL reference). The test docstring said "central correctness"; it now states the real scope and points at AC3. Finding 1 (custody, test/CLI side). The e2e custody test asserts the plaintext input master is deleted at blinding time. The CLI grows --keep-input and its docstring no longer claims the true vector never survives when it did. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- .github/workflows/deploy-image.yml | 23 +++++++--- scripts/blind_data_vector.py | 25 +++++++++-- src/sp_validation/tests/test_blinding.py | 55 +++++++++++++++--------- 3 files changed, 71 insertions(+), 32 deletions(-) diff --git a/.github/workflows/deploy-image.yml b/.github/workflows/deploy-image.yml index 13f0fafe..191e911b 100644 --- a/.github/workflows/deploy-image.yml +++ b/.github/workflows/deploy-image.yml @@ -46,16 +46,25 @@ jobs: # 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). + # image can actually load it (sacc + patched firecrown + the UNIONS-WL + # smokescreen fork + cosmo_numba). `draw_param_shifts` is a FORK-ONLY + # symbol: importing it guarantees a resolution to the upstream LSSTDESC + # PyPI package (which self-reports the same version number as the fork) + # can never ship green. - name: Blinding-stack import smoke test - run: docker run --rm ${{ steps.meta.outputs.tags }} python -c "import sacc; import firecrown.likelihood; import smokescreen" + run: docker run --rm ${{ steps.meta.outputs.tags }} python -c "import sacc; import firecrown.likelihood; from smokescreen.param_shifts import draw_param_shifts; import cosmo_numba" - # Run the fast test suite against the freshly-built image *before* + # Run the FULL 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 - # needs no extra setup. - - name: Run unit tests - run: docker run --rm ${{ steps.meta.outputs.tags }} python -m pytest src/sp_validation/tests -m "not slow" + # full stack — sacc, patched firecrown, the UNIONS-WL smokescreen fork, + # and cosmo_numba (proven by the smoke test above) — plus the test files + # (COPY . + editable install), so the theory-heavy blinding ACs and the + # CAMB↔CCL cross-checks (all `@pytest.mark.slow`, and the cosmo_numba- + # gated AC1/4/5/9) all execute here. This is the environment the blinding + # PRD names as the one where "the full suite passes inside the container + # and CI runs it in the freshly built image before publish". + - name: Run unit tests (full suite) + run: docker run --rm ${{ steps.meta.outputs.tags }} python -m pytest src/sp_validation/tests - name: Push uses: docker/build-push-action@v6 diff --git a/scripts/blind_data_vector.py b/scripts/blind_data_vector.py index 7defc657..ca937a49 100644 --- a/scripts/blind_data_vector.py +++ b/scripts/blind_data_vector.py @@ -9,8 +9,10 @@ ``blind`` draws a seed from OS entropy, applies the concealing shift *immediately* (the seed is never printed or logged), and writes the blinded master plus a public commitment JSON (safe to commit to the repo) and an -encrypted seed+truth bundle; no plaintext seed or true vector survives on -disk. ``unblind`` verifies the commitment and restores the true vector. +encrypted seed+truth bundle. The plaintext input master is then deleted +(pass ``--keep-input`` to retain it); no plaintext seed or true vector +survives on disk. ``unblind`` verifies the commitment and restores the true +vector. ``verify`` is a cheap, seedless check that a blinded file matches a commitment. @@ -69,8 +71,17 @@ def _blind(args): + "\n ".join(clashes) + "\nPick a fresh --output-dir or --label (never overwrite a blind)." ) - blinding.blind_file(args.master, str(out_dir), config=config, label=args.label) - print("Commit the commitment JSON to the repo; keep the bundle + key safe.") + blinding.blind_file( + args.master, + str(out_dir), + config=config, + label=args.label, + keep_input=args.keep_input, + ) + print( + "Commit the commitment JSON to the repo; keep the bundle + key safe " + "and separated (colocation in the output dir is not at-rest protection)." + ) def _unblind(args): @@ -119,6 +130,12 @@ def main(argv=None): p.add_argument("master", help="master SACC file to blind") p.add_argument("-o", "--output-dir", required=True) p.add_argument("--label", default="A", help="blind label (default A)") + p.add_argument( + "--keep-input", + action="store_true", + help="retain the plaintext input master (default: delete it " + "after blinding — the true vector is escrowed in the bundle)", + ) p.set_defaults(func=_blind) else: p.add_argument("blinded", help="blinded master SACC file") diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py index 8dadf804..92c54e0a 100644 --- a/src/sp_validation/tests/test_blinding.py +++ b/src/sp_validation/tests/test_blinding.py @@ -44,11 +44,6 @@ def _fine_theta(n=80): return np.geomspace(0.1, 300.0, n) -def _report_theta(n=12): - """The pure-E/B REPORTING grid — a strict sub-range of the fine grid.""" - return np.geomspace(5.0, 100.0, n) - - def _xi_template(theta, k=0): """Smooth synthetic ξ± for pair index ``k`` (no CCL needed).""" theta = np.asarray(theta) @@ -70,7 +65,6 @@ def make_master( with_pure_eb=False, with_rho=True, b_amplitude=0.0, - report_theta=None, nmodes=6, ): """Build a master SACC in the sacc_io layout, synthetic values throughout. @@ -83,7 +77,6 @@ def make_master( """ nz = {i: _gauss_nz(0.5 + 0.3 * i, 0.15 + 0.02 * i) for i in range(nbins)} ctheta, ftheta = _coarse_theta(), _fine_theta() - rtheta = _report_theta() if report_theta is None else report_theta pairs = [(i, j) for i in range(nbins) for j in range(i, nbins)] s = sio.new_sacc(nz, metadata={"catalogue_version": "vTEST"}) @@ -146,22 +139,30 @@ def make_master( blocks.append((idx, np.eye(len(idx)) * 1e-20)) if with_pure_eb: - for i, j in pairs: + # Seed the pure-E/B blocks through EXACTLY the pipeline convention the + # re-derivation reconstructs (blinding.rederive_eb_from_fine_xi): the + # reporting grid IS the coarse ξ± grid, the reporting ξ± are the + # MEASURED coarse ξ± (not an interpolation of the fine grid), and the + # integration ξ± are the fine grid. tmin/tmax are the edge-based + # integration bounds; the outermost reporting point sits at tmax with + # no interior support, so it comes back NaN — the AC9 boundary case. + # Seeding through this convention (rather than the blinding module's + # own interp/pseudo-edge seam) is what makes AC1/AC4/AC9 a real check + # of the re-derivation against the pipeline convention, not a tautology. + tmin, tmax = float(ctheta[0]), float(ctheta[-1]) + for k, (i, j) in enumerate(pairs): xip_f, xim_f = fine_xi[(i, j)] - xip_r = np.interp(rtheta, ftheta, xip_f) - xim_r = np.interp(rtheta, ftheta, xim_f) - left, right = bd._grid_edges(rtheta) + xip_r, xim_r = _xi_template(ctheta, k) # measured coarse reporting ξ± modes = bd.rederive_pure_eb( - rtheta, - xip_r, - xim_r, - ftheta, - xip_f, - xim_f, - float(left[0]), - float(right[-1]), + ctheta, xip_r, xim_r, ftheta, xip_f, xim_f, tmin, tmax + ) + sio.add_pure_eb( + s, + (i, j), + ctheta, + *(modes[key] for key in sio.PURE_KEYS), + bounds=(tmin, tmax), ) - sio.add_pure_eb(s, (i, j), rtheta, *(modes[k] for k in sio.PURE_KEYS)) tr = sio._pair((i, j)) idx = np.concatenate( [s.indices(sio.PURE_TYPES[k], tr) for k in sio.PURE_KEYS] @@ -462,7 +463,16 @@ def test_cli_blind_refuses_existing_output(tmp_path): def test_ac2_on_file_shift_equals_theory_difference(): """AC2: per-row shift on the blinded file == theory_fn(hidden) − theory_fn(fiducial), hidden recovered by re-running the fork's draw — - for all three blocks, on a two-bin fixture.""" + for all three blocks, on a two-bin fixture. + + Scope: this verifies fork-draw recovery + merge 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`` via + :func:`bd._blindable_blocks`, so any wrong-cosmology dependence cancels + and a wrong backend would still pass here. Backend correctness is carried + by AC3 (:func:`test_ac3_cross_backend_against_independent_ccl_reference`), + which builds an independent ``ccl.Cosmology`` reference.""" master = make_master(nbins=2, with_cl=True, with_fine=True) cfg = bd.BlindingConfig() seed = "ac2-seed" @@ -739,6 +749,9 @@ def test_ac6_ac8_end_to_end_blind_unblind_custody(tmp_path): paths = bd.blind_file(str(master_path), str(out), label="A", log=_NOLOG) # -- custody hygiene (AC6) ------------------------------------------- -- + # The plaintext input master is deleted at blinding time — the true vector + # survives only inside the encrypted bundle. + assert not master_path.exists(), "plaintext input master was not deleted" blinded = sio.load(paths["blinded"]) assert blinded.metadata["concealed"] is True assert "seed_smokescreen" not in blinded.metadata From f385c2d35d7fb360fee51a229e399c0479a69d46 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 12:26:12 +0200 Subject: [PATCH 08/17] test(blinding): reframe AC9 to NaN-parity under the blind MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit With the re-derivation corrected to feed the measured coarse reporting ξ± and the fine integration grid (not the old interp + pseudo-edge seam), the Schneider estimator no longer degenerates at the reporting boundary — interior support exists at every reporting point, matching the PRD's own statement that on production files (reporting grid a strict sub-range of the integration range) the degeneracy never fires. The old AC9 asserted a NaN MUST appear at the outermost point, which was an artifact of the pre-rework seam; that assertion now fails against the corrected fixture. AC9 now asserts the real invariant: the re-derived NaN pattern is identical on the true and blinded files (blinding shifts inputs, never moves estimator support), whether or not any point is degenerate. Verified with cosmo_numba installed: full blinding suite 25 passed, 0 skipped. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/tests/test_blinding.py | 32 ++++++++++++++---------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py index 92c54e0a..a6038d19 100644 --- a/src/sp_validation/tests/test_blinding.py +++ b/src/sp_validation/tests/test_blinding.py @@ -700,29 +700,35 @@ def test_ac5_untouched_blocks_and_row_order(): @pytest.mark.slow -def test_ac9_pure_eb_boundary_nan_parity(): - """AC9: the boundary-degenerate reporting point — the point nearest the - edge-based integration bound, where the Schneider estimator has no - interior support (observed: the LAST reporting point, at the tmax edge) - — is NaN in both the true-file and blinded-file re-derivations. - NaN-parity with the pipeline's own estimator, never a spurious finite.""" +def test_ac9_pure_eb_nan_parity_under_blind(): + """AC9: the re-derived pure-E/B NaN pattern is identical on the true and + blinded files — 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 file (on production files it is empty — the + reporting grid is a strict sub-range of the integration range, so the + degeneracy never fires; this is what the corrected re-derivation, feeding + the measured coarse reporting ξ± and the fine integration grid, delivers), + the blinded-file re-derivation must reproduce it bit-for-bit: the blind is + a pure shift of the inputs, not a change of estimator support. The finite + values move (the ξ± shifted); the NaN mask does not.""" pytest.importorskip("cosmo_numba") master = make_master(with_cl=False, with_fine=True, with_pure_eb=True) tr = sio._pair((0, 0)) idx_e = master.indices(sio.PURE_TYPES["xip_E"], tr) true_vals = np.array(master.mean)[idx_e] - assert np.isnan(true_vals[-1]), ( - "fixture did not trigger the boundary degeneracy — the outermost " - "reporting point should be NaN in the pipeline estimator" - ) - assert np.all(np.isfinite(true_vals[:-1])) blinded = bd.blind_sacc(master, "ac9-seed", log=_NOLOG) blind_vals = np.array(blinded.mean)[idx_e] assert np.array_equal(np.isnan(true_vals), np.isnan(blind_vals)), ( - "blinding moved the NaN pattern of the degenerate boundary point" + "blinding moved the pure-E/B NaN pattern — the estimator support changed" + ) + # the finite values did move (the blind actually shifted the ξ±) + finite = np.isfinite(true_vals) + assert finite.any() and not np.allclose( + true_vals[finite], blind_vals[finite], atol=0 ) - assert np.all(np.isfinite(blind_vals[:-1])) # --------------------------------------------------------------------------- # From aa3464980b76983eb6914b7b81e17b8ddbb57c1b Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 12:36:45 +0200 Subject: [PATCH 09/17] fix(blinding): install NumbaQuadpack at top level so cosmo_numba's kernels load in CI MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Turning on the full suite in CI surfaced a real gap the fast suite hid: the image had cosmo_numba but not NumbaQuadpack, so the COSEBIs / pure-E/B kernels (cosmo_numba.B_modes.cosebis / schneider2022, which `from NumbaQuadpack import dqags`) died at import and the five kernel-backed ACs (AC1/4/5/9 + e2e AC6/8) failed. cosmo_numba declares NumbaQuadpack as a direct-URL (git+) dependency, and pip/uv refuse to install a direct-URL *transitive* dep — so `import cosmo_numba` (light top-level) succeeded while the kernels were absent. Name NumbaQuadpack at top level in the `blinding` extra so the Dockerfile install resolves it, and harden the CI blinding-stack smoke test to import the kernel modules (COSEBIS, get_pure_EB_modes) rather than just the top-level package — so a missing NumbaQuadpack can never again pass the smoke test and only blow up deep in the slow ACs. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- .github/workflows/deploy-image.yml | 8 +++++++- pyproject.toml | 11 ++++++++++- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/.github/workflows/deploy-image.yml b/.github/workflows/deploy-image.yml index 191e911b..91097f7d 100644 --- a/.github/workflows/deploy-image.yml +++ b/.github/workflows/deploy-image.yml @@ -51,8 +51,14 @@ jobs: # symbol: importing it guarantees a resolution to the upstream LSSTDESC # PyPI package (which self-reports the same version number as the fork) # can never ship green. + # `import cosmo_numba` only pulls the light top-level package; the COSEBIs + # / pure-E/B kernels live in cosmo_numba.B_modes.*, which import + # NumbaQuadpack (a direct-URL transitive dep pip/uv won't auto-install). + # Import the kernel modules here so a missing NumbaQuadpack — which would + # otherwise pass a top-level `import cosmo_numba` and only fail deep in the + # slow ACs — can never ship green. - name: Blinding-stack import smoke test - run: docker run --rm ${{ steps.meta.outputs.tags }} python -c "import sacc; import firecrown.likelihood; from smokescreen.param_shifts import draw_param_shifts; import cosmo_numba" + run: docker run --rm ${{ steps.meta.outputs.tags }} python -c "import sacc; import firecrown.likelihood; from smokescreen.param_shifts import draw_param_shifts; from cosmo_numba.B_modes.cosebis import COSEBIS; from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes" # Run the FULL test suite against the freshly-built image *before* # pushing, so a failing suite blocks publication. The image carries the diff --git a/pyproject.toml b/pyproject.toml index 9b2dd825..0bd2eaa3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -129,7 +129,12 @@ glass = [ # cosmo_numba supplies the COSEBIs / pure-E/B kernels the blind re-derives # derived statistics through (same kernels b_modes drives). Pinned to the # numpy-2.x-compatible fork commit (objmode np.fft fix; upstream aguinot/main -# breaks under the numpy 2.4 this extra requires). +# breaks under the numpy 2.4 this extra requires). cosmo_numba's own +# NumbaQuadpack dependency is a git-URL requirement, and pip/uv refuse to +# install a direct-URL *transitive* dep — so `import cosmo_numba` succeeds but +# `cosmo_numba.B_modes.cosebis` dies on `import NumbaQuadpack`. NumbaQuadpack +# must therefore be a TOP-LEVEL requirement here, or the COSEBIs/pure-E/B +# kernels are missing at runtime while the light top-level import looks fine. # # firecrown is not on PyPI and declares conda-forge-only / unused sampler # connectors as hard deps, so installing this extra requires the dependency @@ -142,6 +147,10 @@ blinding = [ "firecrown @ git+https://github.com/LSSTDESC/firecrown.git@v1.15.1", "smokescreen @ git+https://github.com/UNIONS-WL/Smokescreen@588a6b9b26560bd5ba3dd5ba342f3c40152644f9", "cosmo-numba @ git+https://github.com/cailmdaley/cosmo-numba@db452c487f3d740c9fe69a409352f11987525d65", + # cosmo_numba's Schneider/COSEBIs kernels import NumbaQuadpack; it is a + # direct-URL transitive dep of cosmo_numba that pip/uv will not pull + # automatically, so it must be named at top level (see the comment above). + "NumbaQuadpack @ git+https://github.com/Nicholaswogan/NumbaQuadpack.git", "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 From 8426e8986fb3c37bc92053738e2291918b68cc6d Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 12:50:04 +0200 Subject: [PATCH 10/17] =?UTF-8?q?refactor(blinding):=20move=20the=20theory?= =?UTF-8?q?=20surface=20to=20blinding=5Ftheory.py=20=E2=80=94=20no=20sp=5F?= =?UTF-8?q?validation/cosmology.py?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `develop` deliberately deleted src/sp_validation/cosmology.py (#223), moving cosmology to cs_util.cosmo; this branch must not resurrect it. Per the UNIONS layering — cs_util is the cosmology *library* (generic machinery), sp_validation holds *configuration* and *survey-specific implementations* — the blinding theory surface belongs in the blinding namespace, not a local cosmology module. Rename cosmology.py → blinding_theory.py. TheoryConfig (configuration) and the three master-layout backends in blinding.py (coarse/fine ξ±, pseudo-Cℓ — survey-specific) are correctly homed here. The generic pieces — the CCL-native ξ± path (xi_ccl/cl_ee), the independent CAMB P(k)→Pk2D path (xi_camb), the σ8/A_s rescale (camb_As_for_sigma8) — live here FOR NOW with an explicit module-docstring note that they are cosmology-library code destined for cs_util.cosmo. Module-scope imports stay numpy-only (dataclasses + numpy). Updated all importers: blinding.py, test_blinding.py, and test_camb_ccl_crosscheck.py (module ref + alias). No sp_validation/cosmology.py anywhere. Affected tests (blinding + CAMB↔CCL cross-check) green with cosmo_numba: 30 passed, 0 skipped. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/blinding.py | 2 +- .../{cosmology.py => blinding_theory.py} | 31 ++++++++++++++----- src/sp_validation/tests/test_blinding.py | 2 +- .../tests/test_camb_ccl_crosscheck.py | 8 ++--- 4 files changed, 30 insertions(+), 13 deletions(-) rename src/sp_validation/{cosmology.py => blinding_theory.py} (92%) diff --git a/src/sp_validation/blinding.py b/src/sp_validation/blinding.py index 2f255020..4ba21e18 100644 --- a/src/sp_validation/blinding.py +++ b/src/sp_validation/blinding.py @@ -54,7 +54,7 @@ import numpy as np from . import sacc_io -from .cosmology import TheoryConfig, cl_ee, xi_ccl, xi_ell_grid +from .blinding_theory import TheoryConfig, cl_ee, xi_ccl, xi_ell_grid # --------------------------------------------------------------------------- # diff --git a/src/sp_validation/cosmology.py b/src/sp_validation/blinding_theory.py similarity index 92% rename from src/sp_validation/cosmology.py rename to src/sp_validation/blinding_theory.py index 91cc96fe..d3783e17 100644 --- a/src/sp_validation/cosmology.py +++ b/src/sp_validation/blinding_theory.py @@ -1,10 +1,27 @@ -"""Fiducial theory configuration and the two ξ± theory paths. - -:Name: cosmology.py - -:Description: The single home for the blinding/inference fiducial cosmology - (:class:`TheoryConfig`) and two independent routes to the tomographic - shear two-point prediction: +"""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` + (coarse ξ±, fine ξ±, 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``: 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 diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py index a6038d19..5cb673fb 100644 --- a/src/sp_validation/tests/test_blinding.py +++ b/src/sp_validation/tests/test_blinding.py @@ -20,7 +20,7 @@ from sp_validation import blinding as bd from sp_validation import sacc_io as sio -from sp_validation.cosmology import TheoryConfig +from sp_validation.blinding_theory import TheoryConfig _NOLOG = lambda *a, **k: None # noqa: E731 diff --git a/src/sp_validation/tests/test_camb_ccl_crosscheck.py b/src/sp_validation/tests/test_camb_ccl_crosscheck.py index 33ec64ff..21296936 100644 --- a/src/sp_validation/tests/test_camb_ccl_crosscheck.py +++ b/src/sp_validation/tests/test_camb_ccl_crosscheck.py @@ -4,11 +4,11 @@ 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.cosmology`: +:mod:`sp_validation.blinding_theory`: -- **Path A** (:func:`~sp_validation.cosmology.xi_ccl`): CCL-native — CCL's +- **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.cosmology.xi_camb`): an independent +- **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. @@ -26,7 +26,7 @@ import numpy as np import pytest -from sp_validation import cosmology as cm +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 From 908992d38cc9204feb293c20bf385d3b0af645c6 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 12:53:26 +0200 Subject: [PATCH 11/17] docs(blinding): link the cs_util migration issue (cs_util#80) in blinding_theory The generic theory machinery destined for cs_util.cosmo is tracked in cs_util#80 ("cosmo: absorb generic theory machinery from sp_validation blinding"); name it in the module docstring so the migration has an explicit anchor. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/blinding_theory.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sp_validation/blinding_theory.py b/src/sp_validation/blinding_theory.py index d3783e17..47e67ca3 100644 --- a/src/sp_validation/blinding_theory.py +++ b/src/sp_validation/blinding_theory.py @@ -15,8 +15,8 @@ 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``: it is cosmology-library code, - not blinding-specific. There is deliberately **no** + 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. From c8ce0bf7aca483878ade60c3b38aba34b5a48df6 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 13:15:17 +0200 Subject: [PATCH 12/17] =?UTF-8?q?fix(sacc=5Fio):=20make=20add=5Fpure=5Feb?= =?UTF-8?q?=20bounds=20optional=20=E2=80=94=20restore=20PR-2=20writer=20co?= =?UTF-8?q?mpat?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The finding-6/7 commit made add_pure_eb's `bounds` a required keyword-only argument, which broke three test_sacc_io.py call sites (plain pure-E/B writes for storage/roundtrip and the non-ascending-theta guard) with a TypeError — caught by the newly-enabled full CI suite, not the blinding-only local runs. add_pure_eb is PR-2's writer surface (feat/sacc-2-sacc-io / #245 owns it). Make `bounds` optional (default None): a plain write omits the tmin/tmax tags, exactly the pre-existing behavior. The tags are stamped only when the caller supplies bounds — which the blinding path does, because blinding re-derivation needs the pipeline's edge-based integration range and TreeCorr's edges are not reconstructible from meanr. A file written without bounds simply cannot be blinded: blinding.rederive_eb_from_fine_xi fails loud on the missing tags rather than reconstructing a pseudo-edge that would not match the pipeline (the finding-4/7 contract). No pseudo-edge default is reintroduced. PR-2 surface note: this adds an optional `bounds` kwarg to add_pure_eb; #245's owner should be aware the writer now carries the pure-E/B integration-bound contract for blinding, opt-in and backward-compatible. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/sacc_io.py | 39 ++++++++++++++++++++---------------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index d615e810..7ace8a55 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -263,7 +263,7 @@ def add_cosebis(s, bins, En, Bn, scale_cut): def add_pure_eb( - s, bins, theta, xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb, *, bounds + s, bins, theta, xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb, *, bounds=None ): """Add pure E/B-mode correlation functions for one tracer pair. @@ -278,25 +278,32 @@ def add_pure_eb( bins : tuple of int Source bin pair ``(i, j)``. theta : array_like - Angular separations (arcmin), shared by all six blocks. **Writer - contract:** the reporting grid must be the coarse ξ± grid's - ``meanr`` — the pipeline feeds the *measured* reporting-binned ξ± - into the Schneider estimator, and blinding reconstructs it from the - stored coarse ξ± block. + Angular separations (arcmin), shared by all six blocks. xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb : array_like The six pure E/B / ambiguous mode arrays at ``theta``. - bounds : tuple of float + bounds : tuple of float, optional ``(tmin, tmax)`` in arcmin — the Schneider-2022 integration range the - estimator was run with. **Writer contract:** the pipeline's - edge-based bounds, ``gg.left_edges[0]`` / ``gg.right_edges[-1]`` of - the reporting TreeCorr correlation - (``b_modes.calculate_pure_eb_correlation``'s convention). Stored on - every point as ``tmin``/``tmax`` tags; TreeCorr's edges are not - reconstructible from ``meanr``, so blinding fails loud without them. + estimator was run with. When given, stamped on every point as + ``tmin``/``tmax`` tags. Omit for a plain pure-E/B write (storage, + roundtrip); supply it when the file must support **blinding + re-derivation**, which needs the pipeline's edge-based bounds + (``gg.left_edges[0]`` / ``gg.right_edges[-1]`` of the reporting + TreeCorr correlation — ``b_modes.calculate_pure_eb_correlation``'s + convention) because TreeCorr's edges are not reconstructible from + ``meanr``. A file written without bounds cannot be blinded: + ``blinding.rederive_eb_from_fine_xi`` fails loud on the missing tags + rather than reconstructing a pseudo-edge that would not match the + pipeline. **Reporting-grid contract (blinding):** when the file is to + be blinded, ``theta`` must be the coarse ξ± grid's ``meanr`` — blinding + reconstructs the measured reporting-binned ξ± from the stored coarse + ξ± block, and asserts the two grids coincide. """ _check_ascending("theta", theta) tracers = _pair(bins) - tmin, tmax = (float(b) for b in bounds) + tags = {} + if bounds is not None: + tmin, tmax = (float(b) for b in bounds) + tags = {"tmin": tmin, "tmax": tmax} values = { "xip_E": xip_E, "xim_E": xim_E, @@ -308,9 +315,7 @@ def add_pure_eb( for key in PURE_KEYS: dtype, arr = PURE_TYPES[key], values[key] for n, th in enumerate(theta): - s.add_data_point( - dtype, tracers, float(arr[n]), theta=float(th), tmin=tmin, tmax=tmax - ) + s.add_data_point(dtype, tracers, float(arr[n]), theta=float(th), **tags) def add_rho(s, k, theta, rho_p, rho_m): From 14585192e5d3ef16495bb6a6c885c6891b14f198 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 13:47:54 +0200 Subject: [PATCH 13/17] feat(b_modes): values-only estimator seams for SACC-borne xi MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit cosebis_from_xi and pure_eb_from_xi expose the cosmo_numba kernels (COSEBIS.cosebis_from_xipm, get_pure_EB_modes) to callers holding xi+/- arrays rather than TreeCorr GGCorrelations — the seam the pipeline uses to derive born-blinded COSEBIs / pure-E/B downstream from a blinded fine-xi SACC part. Bit-matches calculate_cosebis' theta_cut path (add_cosebis retained-centre scale-cut contract) and calculate_pure_eb_correlation's edge-based tmin/tmax convention. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/b_modes.py | 73 ++++++++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py index 13f4c2e9..115723c2 100644 --- a/src/sp_validation/b_modes.py +++ b/src/sp_validation/b_modes.py @@ -255,6 +255,79 @@ def pure_EB(corrs): return results +def cosebis_from_xi(theta, xip, xim, nmodes, scale_cut=None): + """COSEBIs (Eₙ, Bₙ) from ξ± arrays through the pipeline kernel (values only). + + The values-only seam of :func:`calculate_cosebis`, for callers holding + ξ± arrays rather than a TreeCorr ``GGCorrelation`` — e.g. deriving + born-blinded COSEBIs from a blinded fine-ξ± SACC part. Calls the same + ``cosmo_numba`` kernel (``COSEBIS.cosebis_from_xipm``) directly on the + values; the covariance/χ² machinery stays with :func:`calculate_cosebis`. + + ``scale_cut`` follows the :func:`sacc_io.add_cosebis` writer contract: + ``(theta_min, theta_max)`` are min/max of the *retained* bin centres + after the pipeline's ``scale_cut_to_bins``. The cut is contiguous in an + ascending grid, so selecting ``theta_min ≤ θ ≤ theta_max`` inclusively + reproduces exactly the retained set, and the kernel is built on that + set's min/max support and fed only the retained ξ± — bit-matching + :func:`calculate_cosebis`'s ``theta_cut``/``xip_cut``/``xim_cut`` path. + Identical inputs ⇒ identical numbers. + """ + from cosmo_numba.B_modes.cosebis import COSEBIS + + theta, xip, xim = (np.asarray(a) for a in (theta, xip, xim)) + tmin, tmax = scale_cut if scale_cut is not None else (theta.min(), theta.max()) + cut = (theta >= tmin) & (theta <= tmax) + theta_cut, xip_cut, xim_cut = theta[cut], xip[cut], xim[cut] + cosebis = COSEBIS( + theta_min=np.min(theta_cut), + theta_max=np.max(theta_cut), + N_max=nmodes, + precision=120, + ) + En, Bn = cosebis.cosebis_from_xipm(theta_cut, xip_cut, xim_cut, parallel=True) + return np.asarray(En), np.asarray(Bn) + + +def pure_eb_from_xi( + theta_report, xip_report, xim_report, theta_int, xip_int, xim_int, tmin, tmax +): + """Pure-E/B correlation functions from ξ± arrays through the pipeline kernel. + + The values-only seam of :func:`calculate_pure_eb_correlation`, for + callers holding ξ± arrays rather than TreeCorr correlations — e.g. + deriving born-blinded pure-E/B from blinded SACC parts. Calls the same + ``cosmo_numba`` kernel (``get_pure_EB_modes``) directly on the values. + The reporting grid must be a strict sub-range of the integration grid; + ``tmin``/``tmax`` are the reporting correlation's TreeCorr *bin edges* + (``gg.left_edges[0]`` / ``gg.right_edges[-1]``) — the pipeline's + convention, carried on SACC files by ``sacc_io.add_pure_eb``. A + reporting point coinciding with the integration boundary is degenerate + (no interior support) and comes back NaN, exactly as + :func:`calculate_pure_eb_correlation` returns it — never a spurious + finite value. + + Returns + ------- + dict + Keyed by ``_EB_KEYS`` (xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb). + """ + from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes + + modes = get_pure_EB_modes( + theta=np.asarray(theta_report), + xip=np.asarray(xip_report), + xim=np.asarray(xim_report), + theta_int=np.asarray(theta_int), + xip_int=np.asarray(xip_int), + xim_int=np.asarray(xim_int), + tmin=tmin, + tmax=tmax, + parallel=True, + ) + return dict(zip(_EB_KEYS, (np.asarray(m) for m in modes))) + + def calculate_cosebis(gg, nmodes=10, scale_cuts=None, cov_path=None): """ Calculate COSEBIs modes from a correlation function for multiple scale cuts. From 1e699cb70d8dea520d6d05a18eff13a6a2a86715 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 13:48:18 +0200 Subject: [PATCH 14/17] refactor(blinding): per-part-at-birth architecture MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the terminal master-blind (master in -> extract/conceal/merge -> re-derive -> delete master) with the three-verb per-part surface ruled in the PRD: - blind_init (once per catalogue version): OS-entropy seed, repo-committable commitment.json (sha256(seed) + canonical config digest), encrypted seed bundle (smokescreen.encryption, Fernet); plaintext seed never written. - blind_part (per intermediate, at birth): conceal one part SACC (coarse xi, fine xi, or pseudo-Cl) through its own theory_fn backend under the shared seed, write the blinded part with a self-contained per-part escrow bundle beside it, delete the plaintext part, stamp concealed / blind_commitment / blind_config_digest, strip seed_smokescreen. - unblind_part (per part, or the assembled file via grid-tag selection): verify both sha256(seed) and the config digest against commitment.json BEFORE subtracting; restore bit-for-bit from the escrow when present. - assert_consistent_blind: assembly-time custody assertion (identical blind_commitment across all blindable parts, fails closed on mismatch or mixed blinded/plaintext), exported for the sacc_io gather. sacc_io gains the terminal gather that assembles standalone part SACCs into the one-file {version}.sacc (tracer merge, ordered points, block-diagonal covariance) — its single blind-aware touch is the assert_consistent_blind call + stamp. Derived statistics are no longer re-derived inside blinding code: COSEBIs and pure-E/B are computed downstream from the already-blinded fine xi part by the pipeline's own estimators (b_modes seams), so they are born blinded. The three theory backends, envelope calibration, custody primitives, and metadata keys are unchanged. Tests reworked to the per-part AC1-9 (AC2 = draw-recovery + placement + one hidden cosmology across parts; AC9 = NaN-parity); CAMB<->CCL cross-check tests untouched. CLI reworked to blind-init / blind-part / unblind (+ seedless verify). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- scripts/blind_data_vector.py | 133 ++-- src/sp_validation/blinding.py | 612 +++++++-------- src/sp_validation/sacc_io.py | 86 ++- src/sp_validation/tests/test_blinding.py | 924 +++++++++++++---------- 4 files changed, 937 insertions(+), 818 deletions(-) diff --git a/scripts/blind_data_vector.py b/scripts/blind_data_vector.py index ca937a49..d0e856ae 100644 --- a/scripts/blind_data_vector.py +++ b/scripts/blind_data_vector.py @@ -2,40 +2,42 @@ """Script blind_data_vector.py -Blind, unblind, or verify a master SACC data-vector file with -:mod:`sp_validation.blinding` (Smokescreen-fork concealment, hash-commitment -custody). - -``blind`` draws a seed from OS entropy, applies the concealing shift -*immediately* (the seed is never printed or logged), and writes the blinded -master plus a public commitment JSON (safe to commit to the repo) and an -encrypted seed+truth bundle. The plaintext input master is then deleted -(pass ``--keep-input`` to retain it); no plaintext seed or true vector -survives on disk. ``unblind`` verifies the commitment and restores the true -vector. -``verify`` is a cheap, seedless check that a blinded file matches a -commitment. +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`` (``sha256(seed)`` + config digest), and encrypts the seed +into a Fernet bundle. ``blind-part`` blinds one intermediate part SACC +(coarse ξ±, fine ξ±, 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 both commitment hashes 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`` (fine rows selected by +the ``grid`` tag). ``verify`` is a cheap, seedless check that a blinded file +matches a commitment. :Authors: Cail Daley Examples -------- -Blind:: +Once per catalogue version:: - blind_data_vector.py blind v1.4.6.3.sacc -o blinded/ --label A + blind_data_vector.py blind-init blinded/ -Unblind:: +Per intermediate part, at birth:: - blind_data_vector.py unblind blinded/v1.4.6.3_blinded_A.sacc \\ - --bundle blinded/blind_A_bundle.encrpt \\ - --key blinded/blind_A_bundle.key \\ - --commitment blinded/blind_A_commitment.json \\ - -o unblinded/v1.4.6.3.sacc + blind_data_vector.py blind-part parts/xi_fine.fits --blind-dir blinded/ + +Unblind one part:: + + blind_data_vector.py unblind parts/xi_fine_blinded.fits \\ + --blind-dir blinded/ -o parts/xi_fine.fits Verify:: - blind_data_vector.py verify blinded/v1.4.6.3_blinded_A.sacc \\ - blinded/blind_A_commitment.json + blind_data_vector.py verify parts/xi_fine_blinded.fits \\ + blinded/commitment.json """ import argparse @@ -56,43 +58,46 @@ def _config_from_args(args): return blinding.BlindingConfig.from_overrides(overrides) -def _blind(args): +def _blind_init(args): config = _config_from_args(args) - out_dir = pathlib.Path(args.output_dir) - out_dir.mkdir(parents=True, exist_ok=True) - # Refuse before any heavy compute: a blind is a one-shot custody event - # and silently overwriting a previous blind's outputs would destroy the - # record tying that blind to its seed. - paths = blinding._output_paths(args.master, str(out_dir), args.label) - clashes = [p for p in paths.values() if pathlib.Path(p).exists()] + 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 output(s):\n " + "refusing to overwrite existing blind state:\n " + "\n ".join(clashes) - + "\nPick a fresh --output-dir or --label (never overwrite a blind)." + + "\nPick a fresh --blind-dir (never overwrite a blind)." ) - blinding.blind_file( - args.master, - str(out_dir), - config=config, - label=args.label, - keep_input=args.keep_input, - ) + 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 output dir is not at-rest protection)." + "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): - config = _config_from_args(args) - blinding.unblind_file( + blinding.unblind_part( args.blinded, - args.bundle, - args.key, - args.commitment, + args.blind_dir, args.output, - config=config, + config=_config_from_args(args), ) @@ -122,31 +127,41 @@ 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", "unblind"): + 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": - p.add_argument("master", help="master SACC file to blind") - p.add_argument("-o", "--output-dir", required=True) + 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 master (default: delete it " - "after blinding — the true vector is escrowed in the bundle)", + 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) + p.set_defaults(func=_blind_part) else: - p.add_argument("blinded", help="blinded master SACC file") - p.add_argument("--bundle", required=True, help="encrypted seed bundle") - p.add_argument("--key", required=True, help="bundle Fernet key file") - p.add_argument("--commitment", required=True, help="commitment JSON") + 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 master SACC file") + p.add_argument("blinded", help="blinded SACC file") p.add_argument("commitment", help="commitment JSON") p.set_defaults(func=_verify) diff --git a/src/sp_validation/blinding.py b/src/sp_validation/blinding.py index 4ba21e18..04d160ee 100644 --- a/src/sp_validation/blinding.py +++ b/src/sp_validation/blinding.py @@ -1,48 +1,53 @@ -"""Blinding — conceal the measured data vector behind a hidden cosmology. +"""Blinding — conceal each intermediate data product behind a hidden cosmology. :Name: blinding.py -:Description: Smokescreen blinding wiring for the master SACC data product. - The measured statistics are shifted by a difference of theory vectors +:Description: Smokescreen blinding wiring, per part, at birth. Each blindable + intermediate SACC product — coarse ξ±, fine ξ±, 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)``). + ``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; each blindable block goes through its - ``ConcealDataVector(fiducial_params, shifts_dict, sacc_data, *, seed, - theory_fn)`` entry point. This module supplies the two - sp_validation-specific pieces: the three ``theory_fn`` backends matching - the master SACC row layout (coarse ξ±, fine ξ±, pseudo-Cℓ) and the SACC - assembly around the concealed vectors. + draw; each part goes through its ``ConcealDataVector(fiducial_params, + shifts_dict, sacc_data, *, seed, theory_fn)`` entry point. This module + supplies the two sp_validation-specific pieces: the three ``theory_fn`` + backends matching each part's row layout (coarse ξ±, fine ξ±, pseudo-Cℓ) + and the SACC handling around the concealed vectors. **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 − Ω_ν``. The same seed yields - the same hidden cosmology across all three block passes (the fork's draw - depends only on ``(key, seed)``), so the three shifted blocks are - mutually consistent by construction. - - **Derived statistics are re-derived, never shifted.** COSEBIs and - pure-E/B are replaced by re-running the pipeline's own estimators - (the ``cosmo_numba`` kernels :mod:`sp_validation.b_modes` drives) on the - blinded fine ξ±. Covariance and ρ/τ PSF diagnostics are byte-identical - to the input — 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.** The blind CLI draws an - OS-entropy seed, blinds immediately, publishes ``sha256(seed)`` plus a - canonical config digest as a repo-committable commitment, encrypts the - seed + true vector into a Fernet bundle (``smokescreen.encryption``) and - deletes both plaintexts — the temporary bundle plaintext *and* the input - master file, so no plaintext true vector survives on disk beside the - blind. Unblind verifies both hashes against the commitment before - subtracting anything, then restores the true vector. + the same hidden cosmology across all part passes (the fork's draw depends + only on ``(key, seed)``), so the blinded parts are mutually consistent by + construction. + + **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 fine ξ± 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 + ``sha256(seed)`` plus a canonical config digest as a repo-committable + ``commitment.json``, and encrypts the seed into a Fernet bundle + (``smokescreen.encryption``) — the plaintext seed is never written. + 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. Terminal assembly + (:func:`sp_validation.sacc_io.gather`) calls + :func:`assert_consistent_blind` to fail closed unless every blindable + part carries the identical ``blind_commitment``. :func:`unblind_part` + verifies both hashes against the commitment *before* subtracting + anything, then restores the true part. """ import dataclasses @@ -155,7 +160,8 @@ def hidden_params(seed, config): — per-key-independent, local RNG) on the calibrated envelope and overlays the deltas on the fiducial, exactly as ``ConcealDataVector`` does internally. Deterministic: same ``(seed, config)`` ⇒ same hidden point, - forever — the reproducibility contract unblinding relies on. + forever — the reproducibility contract unblinding relies on, and what + makes every part share one hidden cosmology under one seed. """ from smokescreen.param_shifts import draw_param_shifts @@ -167,7 +173,7 @@ def hidden_params(seed, config): # --------------------------------------------------------------------------- # -# Block discovery on the master SACC +# 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).""" @@ -199,7 +205,7 @@ def cl_pairs(s): def _xi_indices(s, grid): - """Master row indices of the ξ± block on ``grid`` (ascending).""" + """Row indices of the ξ± block on ``grid`` (ascending).""" return np.sort( np.concatenate( [ @@ -211,7 +217,7 @@ def _xi_indices(s, grid): def _cl_ee_indices(s): - """Master row indices of the pseudo-Cℓ_EE block (ascending). + """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. @@ -299,12 +305,15 @@ def theory_fn(params): # Extract → conceal → merge, per block # --------------------------------------------------------------------------- # def _blindable_blocks(s): - """The blindable blocks of a master SACC as ``(name, indices, factory)``. - - ``indices`` are each block's recorded master row indices (ascending, so - the extracted sub-SACC preserves master row order); ``factory`` builds - the matching ``theory_fn`` from the extracted sub-SACC. Blocks absent - from the file are simply not listed. + """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 fine rows are selected by + the ``grid`` tag — the layout contract's per-block tag selection). + ``indices`` are each block's recorded row indices (ascending, so the + extracted sub-SACC preserves row order); ``factory`` builds the matching + ``theory_fn`` from the extracted sub-SACC. Blocks absent from the file + are simply not listed. """ blocks = [] for grid in ("coarse", "fine"): @@ -365,40 +374,46 @@ def _concealed(s): return bool(s.metadata.get("concealed")) -def blind_sacc(master, seed, config=None, label="A", log=print): - """Return a blinded copy of a master SACC (covariance/ρτ/n(z) untouched). +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 (coarse ξ±, fine ξ±, pseudo-Cℓ_EE): extract the - block into a sub-SACC, conceal through the fork, write the shifted - values back at their recorded master indices (row order preserved). - COSEBIs and pure-E/B are then re-derived from the blinded fine ξ± - through the pipeline estimators; provenance is stamped and any leaked - seed key stripped. + Per blindable block present (a standalone part carries exactly one — + coarse ξ±, fine ξ±, or pseudo-Cℓ_EE): extract the block into a sub-SACC, + conceal through the fork, write the shifted values back at their + recorded indices (row order preserved). 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(master): + if _concealed(part): raise ValueError("already concealed — unblind first") - _check_master(master) + blocks = _blindable_blocks(part) + if not blocks: + raise ValueError( + "no blindable block (coarse/fine ξ± or pseudo-Cℓ_EE) in this SACC " + "— ρ/τ diagnostic parts are never blinded" + ) - blinded = master.copy() - for name, indices, factory in _blindable_blocks(master): - factor = _concealing_factor(master, indices, factory, config, seed) + 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") - rederive_eb_from_fine_xi(blinded, log=log) _stamp_provenance(blinded, seed_commitment(seed), label, config.config_digest()) return blinded def unblind_sacc(blinded, seed, config=None, log=print): - """Recover the true master SACC from a blinded one + the revealed ``seed``. + """Recover the true part SACC from a blinded one + the revealed ``seed``. Verifies ``sha256(seed)`` and the config digest against the stamped metadata (loud failure on either mismatch — verification precedes subtraction), recomputes each block's shift from the seed through the - same backends, subtracts it, and re-derives COSEBIs/pure-E/B from the - unblinded fine ξ±. + same backends, and subtracts it. Works on a standalone part or on the + assembled file (fine 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 fine ξ±. """ config = config or BlindingConfig() if not _concealed(blinded): @@ -417,41 +432,16 @@ def unblind_sacc(blinded, seed, config=None, log=print): hidden = hidden_params(seed, config) fiducial = config.theory.ccl_params() - master = blinded.copy() + part = blinded.copy() for name, indices, factory in _blindable_blocks(blinded): theory = factory(_extract_block(blinded, indices), config.theory) factor = theory(hidden) - theory(fiducial) - _set_values(master, indices, np.asarray(master.mean)[indices] - factor) + _set_values(part, indices, np.asarray(part.mean)[indices] - factor) log(f"[unblind] {name}: subtracted {len(indices)} shifts") - rederive_eb_from_fine_xi(master, log=log) for key in ("concealed", "blind", "blind_commitment", "blind_config_digest"): - master.metadata.pop(key, None) - return master - - -def _check_master(master): - """Loud guards: the master must carry a complete, blindable ξ± content.""" - n_plus = len(master.indices(sacc_io.XI_PLUS, grid="coarse")) - n_minus = len(master.indices(sacc_io.XI_MINUS, grid="coarse")) - if not n_plus or not n_minus: - raise ValueError( - "master SACC has no coarse ξ± — cannot blind (both ξ+ and ξ− are " - f"required; found ξ+={n_plus}, ξ−={n_minus})" - ) - if n_plus != n_minus: - raise ValueError( - f"coarse ξ+ ({n_plus}) and ξ− ({n_minus}) row counts differ — " - "refusing an ambiguous coarse block" - ) - has_derived = len(master.indices(sacc_io.COSEBI_EE)) or len( - master.indices(sacc_io.PURE_TYPES["xip_E"]) - ) - if has_derived and not len(master.indices(sacc_io.XI_PLUS, grid="fine")): - raise ValueError( - "master SACC carries COSEBIs/pure-E/B but no grid='fine' ξ± — " - "the derived statistics cannot be re-derived after blinding" - ) + part.metadata.pop(key, None) + return part def _stamp_provenance(s, commitment, label, config_digest): @@ -471,205 +461,116 @@ def _stamp_provenance(s, commitment, label, config_digest): # --------------------------------------------------------------------------- # -# COSEBIs / pure-E/B re-derivation through the pipeline estimators +# Assembly-time custody: one blind across all parts # --------------------------------------------------------------------------- # -def rederive_cosebis(theta, xip, xim, nmodes, scale_cut=None): - """COSEBIs (Eₙ, Bₙ) from ξ± through the pipeline kernel (values only). - - Calls the same kernel :func:`sp_validation.b_modes.calculate_cosebis` - drives (``cosmo_numba``'s ``COSEBIS.cosebis_from_xipm``) directly on the - ξ± values — the covariance/χ² machinery of ``calculate_cosebis`` needs a - patched GGCorrelation we don't have and don't want (blinding never - touches uncertainty), so this is the values-only seam one layer down. - - ``scale_cut`` follows the :func:`sacc_io.add_cosebis` writer contract: - ``(theta_min, theta_max)`` are min/max of the *retained* bin centres - after the pipeline's ``scale_cut_to_bins``. The cut is contiguous in an - ascending grid, so selecting ``theta_min ≤ θ ≤ theta_max`` inclusively - reproduces exactly the retained set, and the kernel is built on that - set's min/max support and fed only the retained ξ± — bit-matching - ``calculate_cosebis``'s ``theta_cut``/``xip_cut``/``xim_cut`` path. - Identical inputs ⇒ identical numbers. - """ - from cosmo_numba.B_modes.cosebis import COSEBIS - - theta, xip, xim = (np.asarray(a) for a in (theta, xip, xim)) - tmin, tmax = scale_cut if scale_cut is not None else (theta.min(), theta.max()) - cut = (theta >= tmin) & (theta <= tmax) - theta_cut, xip_cut, xim_cut = theta[cut], xip[cut], xim[cut] - cosebis = COSEBIS( - theta_min=np.min(theta_cut), - theta_max=np.max(theta_cut), - N_max=nmodes, - precision=120, - ) - En, Bn = cosebis.cosebis_from_xipm(theta_cut, xip_cut, xim_cut, parallel=True) - return np.asarray(En), np.asarray(Bn) - - -def rederive_pure_eb( - theta_report, xip_report, xim_report, theta_int, xip_int, xim_int, tmin, tmax -): - """Pure-E/B correlation functions from ξ± through the pipeline kernel. - - Calls the kernel :func:`sp_validation.b_modes.calculate_pure_eb_correlation` - drives (``cosmo_numba``'s ``get_pure_EB_modes``) directly on the ξ± - values. The reporting grid must be a strict sub-range of the integration - grid; ``tmin``/``tmax`` are the reporting correlation's TreeCorr *bin - edges* (``gg.left_edges[0]`` / ``gg.right_edges[-1]``) — the pipeline's - convention, carried on the file by ``sacc_io.add_pure_eb``. A reporting point - coinciding with the integration boundary is degenerate (no interior - support) and comes back NaN, exactly as the pipeline's own estimator - returns it — never a spurious finite value. +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`` + (they were blinded under different seeds or configs and must never be + combined). If no part is blinded (a pre-blinding or mock assembly), + returns ``None``. Returns ------- - dict - Keyed by ``sacc_io.PURE_KEYS`` (xip_E, xim_E, xip_B, xim_B, xip_amb, - xim_amb). + dict or None + The shared blind metadata (``concealed``, ``blind``, + ``blind_commitment``, ``blind_config_digest``) for the gather to + stamp on the assembled file, or ``None`` when nothing is blinded. """ - from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes - - modes = get_pure_EB_modes( - theta=np.asarray(theta_report), - xip=np.asarray(xip_report), - xim=np.asarray(xim_report), - theta_int=np.asarray(theta_int), - xip_int=np.asarray(xip_int), - xim_int=np.asarray(xim_int), - tmin=tmin, - tmax=tmax, - parallel=True, - ) - return dict(zip(sacc_io.PURE_KEYS, (np.asarray(m) for m in modes))) + blindable = [p for p in parts if _blindable_blocks(p)] + concealed = [p for p in blindable if _concealed(p)] + if not concealed: + 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)" + ) + stamps = { + ( + p.metadata["blind"], + p.metadata["blind_commitment"], + p.metadata["blind_config_digest"], + ) + for p in concealed + } + if len(stamps) != 1: + raise ValueError( + "parts carry different blind commitments — they were blinded " + "under different seeds or configs and must never be combined: " + + "; ".join(f"({label}, {c[:12]}…, {d[:12]}…)" for label, c, d in stamps) + ) + ((label, commitment, digest),) = stamps + return { + "concealed": True, + "blind": label, + "blind_commitment": commitment, + "blind_config_digest": digest, + } -def _blocks_pairs(s, dtype): - """Source pairs ``(i ≤ j)`` carrying ``dtype`` points in ``s``.""" - bins = source_bins(s) - return [ - (i, j) - for a, i in enumerate(bins) - for j in bins[a:] - if len(s.indices(dtype, sacc_io._pair((i, j)))) - ] +# --------------------------------------------------------------------------- # +# 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 rederive_eb_from_fine_xi(master, log=print): - """Replace COSEBIs and pure-E/B in ``master`` from its fine ξ±, in place. +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", + } - For each COSEBI / pure-E/B block present, re-derive its values through - the pipeline estimators on the master's own ξ± of the matching source - pair — COSEBIs from the ``grid='fine'`` ξ±, pure-E/B from the measured - ``grid='coarse'`` reporting ξ± plus the fine integration ξ±, exactly the - inputs ``b_modes`` hands the kernels — and overwrite them. Shared by - blind and unblind (the ξ± differ; the re-derivation is identical). Only - VALUES change; the covariance blocks stay untouched. - """ - for i, j in _blocks_pairs(master, sacc_io.COSEBI_EE): - theta, xip, xim = sacc_io.get_xi(master, (i, j), grid="fine") - tr = sacc_io._pair((i, j)) - idx_e = master.indices(sacc_io.COSEBI_EE, tr) - scale_cut = ( - master.data[idx_e[0]].tags["theta_min"], - master.data[idx_e[0]].tags["theta_max"], - ) - En, Bn = rederive_cosebis(theta, xip, xim, len(idx_e), scale_cut=scale_cut) - _set_values(master, idx_e, En) - _set_values(master, master.indices(sacc_io.COSEBI_BB, tr), Bn) - log(f"[blind] COSEBIs {tr}: re-derived {len(En)} Eₙ + {len(Bn)} Bₙ") - for i, j in _blocks_pairs(master, sacc_io.PURE_TYPES["xip_E"]): - theta_int, xip_int, xim_int = sacc_io.get_xi(master, (i, j), grid="fine") - tr = sacc_io._pair((i, j)) - theta_report, _ = sacc_io.get_pure_eb(master, (i, j)) - # The pipeline feeds get_pure_EB_modes the MEASURED reporting-binned - # ξ± (gg.xip/gg.xim at gg.meanr) alongside the fine integration grid - # — never an interpolation of the fine ξ±. The master's coarse ξ± - # block IS that measurement (and carries the blind consistently), so - # the reporting grids must coincide; anything else means the seam - # cannot reconstruct the pipeline's inputs, so fail loud. - theta_coarse, xip_report, xim_report = sacc_io.get_xi( - master, (i, j), grid="coarse" - ) - if len(theta_coarse) != len(theta_report) or not np.allclose( - theta_coarse, theta_report, rtol=1e-12, atol=0.0 - ): - raise ValueError( - f"pure-E/B reporting grid for {tr} does not match the coarse " - "ξ± grid — the measured reporting ξ± cannot be reconstructed " - "and the pure-E/B block cannot be re-derived" - ) - # Integration bounds: the pipeline's TreeCorr edge-based tmin/tmax, - # stored by the writer (sacc_io.add_pure_eb bounds contract) — edges - # are not reconstructible from meanr, so their absence is fatal. - tags = master.data[master.indices(sacc_io.PURE_TYPES["xip_E"], tr)[0]].tags - if "tmin" not in tags or "tmax" not in tags: - raise ValueError( - f"pure-E/B block for {tr} carries no tmin/tmax integration-" - "bound tags (sacc_io.add_pure_eb bounds contract) — cannot " - "re-derive with the pipeline's edge-based integration range" - ) - modes = rederive_pure_eb( - theta_report, - xip_report, - xim_report, - theta_int, - xip_int, - xim_int, - float(tags["tmin"]), - float(tags["tmax"]), - ) - for key in sacc_io.PURE_KEYS: - _set_values(master, master.indices(sacc_io.PURE_TYPES[key], tr), modes[key]) - log(f"[blind] pure-E/B {tr}: re-derived 6 blocks") +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: -# --------------------------------------------------------------------------- # -# File-level custody: blind / unblind with commitment + encrypted bundle -# --------------------------------------------------------------------------- # -def blind_file( - master_path, out_dir, config=None, label="A", keep_input=False, log=print -): - """Blind a master SACC file end-to-end under hash-commitment custody. - - 1. Draw an OS-entropy seed. - 2. Blind immediately (all blindable blocks; :func:`blind_sacc`). - 3. Write ``blind_{label}_commitment.json`` (repo-committable): - ``sha256(seed)`` + the canonical config digest. - 4. Encrypt the seed + the true data vector into a Fernet bundle - (``smokescreen.encryption``); the plaintext bundle is deleted. - 5. Write the blinded master, then **delete the input master file** — the - plaintext true vector is escrowed (encrypted) in the bundle and must - not survive on disk beside the blind (PRD custody criterion). Pass - ``keep_input=True`` to retain it (and own the custody implication). - No plaintext seed and no plaintext true vector survive on disk. - - Custody caveat: the encrypted bundle and its Fernet key are written to the - *same* ``out_dir``. Fernet is only as protective as the key's separation — - anyone with both files can decrypt the seed and the true vector. Keep the - key out-of-band from the bundle; ``out_dir`` colocation is convenience, not - at-rest protection. + 1. Draw an OS-entropy seed (never written in plaintext, never returned). + 2. Write ``commitment.json`` (repo-committable): ``sha256(seed)`` + the + canonical config digest + the blind label. + 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 of everything written: ``blinded``, ``commitment``, - ``bundle``, ``key``. + Paths written: ``commitment``, ``bundle``, ``key``. """ config = config or BlindingConfig() - master = sacc_io.load(master_path) - paths = _output_paths(master_path, out_dir, label) + paths = init_paths(blind_dir) 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; choose another label or directory" + f"refusing to overwrite existing blind state {path} — a blind " + "is a one-shot custody event; choose another directory" ) seed = secrets.token_hex(16) - blinded = blind_sacc(master, seed, config=config, label=label, log=log) - commitment = { "label": label, "seed_sha256": seed_commitment(seed), @@ -677,110 +578,151 @@ def blind_file( } 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}) - _write_encrypted_bundle(paths, seed, np.asarray(master.mean, dtype=float), label) - sacc_io.save(blinded, paths["blinded"]) - - # The true vector is now escrowed (encrypted) in the bundle; the plaintext - # master must not survive on disk beside the blind (PRD custody criterion). - if not keep_input: - os.remove(master_path) - - log(f"[blind] wrote {paths['blinded']}") - log(f"[blind] commitment (repo-committable): {paths['commitment']}") - log(f"[blind] encrypted bundle + key: {paths['bundle']}, {paths['key']}") - if keep_input: - log(f"[blind] input master RETAINED at {master_path} (keep_input=True)") - else: - log(f"[blind] deleted plaintext input master {master_path}") + log(f"[blind-init] commitment (repo-committable): {paths['commitment']}") + log(f"[blind-init] encrypted seed bundle + key: {paths['bundle']}, {paths['key']}") log( - "[blind] custody: keep the bundle key out-of-band from the bundle " - "(colocation in the output dir is not at-rest protection)" + "[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 unblind_file( - blinded_path, - bundle_path, - key_path, - commitment_path, - out_path, - config=None, - log=print, -): - """Unblind a blinded master SACC file, verifying the commitment first. - - Decrypt the bundle → verify ``sha256(seed)`` and the config digest - against ``commitment.json`` (fail closed on either) → recompute the - shifts from the seed and subtract (:func:`unblind_sacc`) → cross-check - the subtraction against the bundle's stored true vector → restore the - true vector bit-for-bit. +def _read_seed(blind_dir, config): + """Decrypt the seed bundle and verify it against the commitment. + + Both ``sha256(seed)`` and the config digest are checked before the seed + is handed to any caller — a tampered bundle or a drifted config fails + loud here, whether the caller is about to blind or to unblind. + + Returns + ------- + tuple + ``(seed, commitment_dict)``. """ - config = config or BlindingConfig() - bundle = _read_encrypted_bundle(bundle_path, key_path) - with open(commitment_path, encoding="utf-8") as f: + 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) if seed_commitment(bundle["seed"]) != commitment["seed_sha256"]: raise ValueError( "bundle seed does not match the committed sha256(seed) — refusing " - "to unblind" + "to proceed" ) if config.config_digest() != commitment["config_digest"]: raise ValueError( "blinding config does not match the committed config digest — " - "refusing to unblind (a wrong envelope or P(k) recipe would " - "silently subtract a wrong shift)" + "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" + ) + + part = sacc_io.load(part_path) + blinded = blind_sacc(part, seed, config=config, label=commitment["label"], log=log) + + _write_encrypted_json( + paths["escrow"], + { + "label": commitment["label"], + "seed_sha256": commitment["seed_sha256"], + "true_mean": np.asarray(part.mean, dtype=float).tolist(), + }, + ) + sacc_io.save(blinded, paths["blinded"]) + 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 both ``sha256(seed)`` and the + config digest against ``commitment.json`` (fail closed on either — + verification precedes subtraction), recomputes the part's shift from the + seed and subtracts it (:func:`unblind_sacc`). When the part's escrow + bundle exists beside the blinded file, the subtraction is cross-checked + against the escrowed true vector and the truth restored bit-for-bit + (float add-then-subtract leaves ~ulp residue; the escrow removes it). + On the assembled file — no escrow — the subtraction stands alone, with + fine rows selected by the ``grid`` tag. + """ + config = config or BlindingConfig() + seed, _ = _read_seed(blind_dir, config) blinded = sacc_io.load(blinded_path) - master = unblind_sacc(blinded, bundle["seed"], config=config, log=log) - - # The subtraction must agree with the escrowed truth; then the truth is - # restored exactly (float add-then-subtract leaves ~ulp residue, and the - # re-derived COSEBIs inherit it — the bundle removes it). - true_mean = np.asarray(bundle["true_mean"], dtype=float) - recovered = np.asarray(master.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 bundle for this file?" + part = unblind_sacc(blinded, seed, config=config, log=log) + + escrow = part_paths(blinded_path.replace("_blinded", "")) + if os.path.exists(escrow["escrow"]): + bundle = _read_encrypted_json(escrow["escrow"], escrow["escrow_key"]) + 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) ) - _set_values(master, np.arange(len(true_mean)), true_mean) - sacc_io.save(master, out_path) - log(f"[unblind] verified (subtraction residual {residual:.2e}); wrote {out_path}") + 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?" + ) + _set_values(part, np.arange(len(true_mean)), true_mean) + log(f"[unblind] escrow verified (subtraction residual {residual:.2e})") + sacc_io.save(part, out_path) + log(f"[unblind] wrote {out_path}") return out_path -def _output_paths(master_path, out_dir, label): - """Label-scoped output paths for one blind.""" - stem, ext = os.path.splitext(os.path.basename(master_path)) - return { - "blinded": os.path.join(out_dir, f"{stem}_blinded_{label}{ext or '.fits'}"), - "commitment": os.path.join(out_dir, f"blind_{label}_commitment.json"), - "bundle": os.path.join(out_dir, f"blind_{label}_bundle.encrpt"), - "key": os.path.join(out_dir, f"blind_{label}_bundle.key"), - } - - -def _write_encrypted_bundle(paths, seed, true_mean, label): - """Encrypt ``{seed, true_mean}`` to ``paths['bundle']``; delete plaintext. +def _write_encrypted_json(encrpt_path, payload): + """Encrypt ``payload`` (JSON) to ``encrpt_path``; no plaintext survives. - ``smokescreen.encryption.encrypt_file`` (Fernet) writes - ``.encrpt`` + ``.key`` side by side and removes the - plaintext (``keep_original=False``). + ``smokescreen.encryption.encrypt_file`` (Fernet) writes ``.encrpt`` + + ``.key`` side by side and removes the plaintext + (``keep_original=False``). """ from smokescreen.encryption import encrypt_file - plaintext = paths["bundle"].replace(".encrpt", ".json") + plaintext = encrpt_path.replace(".encrpt", ".json") with open(plaintext, "w", encoding="utf-8") as f: - json.dump({"label": label, "seed": seed, "true_mean": true_mean.tolist()}, f) + json.dump(payload, f) encrypt_file(plaintext, save_file=True, keep_original=False) -def _read_encrypted_bundle(bundle_path, key_path): - """Decrypt and parse the seed/true-vector bundle.""" +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(bundle_path, key_path).decode("utf-8")) + return json.loads(decrypt_file(encrpt_path, key_path).decode("utf-8")) diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index 7ace8a55..4acbad80 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -245,7 +245,8 @@ def add_cosebis(s, bins, En, Bn, scale_cut): kernel's actual support — ``min``/``max`` of the *retained* fine-grid ``meanr`` after ``b_modes.scale_cut_to_bins`` (the pipeline builds the kernel from the cut bin centres, per Axel's recommendation) — - not the requested cut. Blinding re-derives COSEBIs from these tags, + not the requested cut. Downstream re-derivation + (``b_modes.cosebis_from_xi`` on a file's fine ξ±) reads these tags, so a requested-cut tag would rebuild the kernel on the wrong support. """ tracers = _pair(bins) @@ -285,18 +286,14 @@ def add_pure_eb( ``(tmin, tmax)`` in arcmin — the Schneider-2022 integration range the estimator was run with. When given, stamped on every point as ``tmin``/``tmax`` tags. Omit for a plain pure-E/B write (storage, - roundtrip); supply it when the file must support **blinding - re-derivation**, which needs the pipeline's edge-based bounds - (``gg.left_edges[0]`` / ``gg.right_edges[-1]`` of the reporting - TreeCorr correlation — ``b_modes.calculate_pure_eb_correlation``'s - convention) because TreeCorr's edges are not reconstructible from - ``meanr``. A file written without bounds cannot be blinded: - ``blinding.rederive_eb_from_fine_xi`` fails loud on the missing tags - rather than reconstructing a pseudo-edge that would not match the - pipeline. **Reporting-grid contract (blinding):** when the file is to - be blinded, ``theta`` must be the coarse ξ± grid's ``meanr`` — blinding - reconstructs the measured reporting-binned ξ± from the stored coarse - ξ± block, and asserts the two grids coincide. + roundtrip); supply it when downstream consumers must be able to + **re-run the estimator from the file** (``b_modes.pure_eb_from_xi``), + which needs the pipeline's edge-based bounds (``gg.left_edges[0]`` / + ``gg.right_edges[-1]`` of the reporting TreeCorr correlation — + ``b_modes.calculate_pure_eb_correlation``'s convention) because + TreeCorr's edges are not reconstructible from ``meanr``. Without the + tags, a re-derivation would have to invent a pseudo-edge that does + not match the pipeline. """ _check_ascending("theta", theta) tracers = _pair(bins) @@ -602,6 +599,69 @@ def extract(s, data_type=None, tracers=None, **tag_filters): return sub +def gather(parts, metadata=None): + """Assemble standalone part SACCs into the one-file ``{version}.sacc``. + + Each part is an intermediate product as it came off its producing rule + (coarse ξ±, fine ξ±, pseudo-Cℓ, ρ/τ, …): its tracers are merged (first + occurrence wins), its points are appended in their stored order with all + tags (including bandpower windows), and its covariance becomes one + diagonal block of the assembled covariance — cross-part blocks are zero, + matching the block-diagonal layout ``assemble_covariance`` builds. Either + every part carries a covariance or none does (mixed is an error). + + **Blind custody (the module's one blind-aware call site):** + :func:`sp_validation.blinding.assert_consistent_blind` is called before + combining — it fails closed unless every blindable part carries the + identical ``blind_commitment``, and its shared stamp (when the parts are + blinded) is written onto the assembled file's metadata. + + Parameters + ---------- + parts : sequence of sacc.Sacc + The part SACCs, in the assembly (covariance) order. + metadata : dict, optional + Key/value pairs stored on the assembled file's metadata. + + Returns + ------- + sacc.Sacc + The assembled file. + """ + from . import blinding + + stamp = blinding.assert_consistent_blind(parts) + + s = sacc.Sacc() + for part in parts: + for name, tracer in part.tracers.items(): + if name not in s.tracers: + s.add_tracer_object(tracer) + for dp in part.data: + s.add_data_point(dp.data_type, dp.tracers, dp.value, **dp.tags) + + with_cov = [part for part in parts if part.covariance is not None] + if with_cov and len(with_cov) != len(parts): + raise ValueError( + f"{len(with_cov)} of {len(parts)} parts carry a covariance — " + "either every part does or none does" + ) + if with_cov: + ntot = len(s.mean) + full = np.zeros((ntot, ntot)) + cursor = 0 + for part in parts: + dense = part.covariance.dense + n = dense.shape[0] + full[cursor : cursor + n, cursor : cursor + n] = dense + cursor += n + s.add_covariance(full) + + for key, value in {**(metadata or {}), **(stamp or {})}.items(): + s.metadata[key] = value + return s + + def save(s, path): """Write ``s`` to ``path`` (FITS), overwriting any existing file.""" s.save_fits(path, overwrite=True) diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py index 5cb673fb..48d86460 100644 --- a/src/sp_validation/tests/test_blinding.py +++ b/src/sp_validation/tests/test_blinding.py @@ -1,15 +1,22 @@ -"""Tests for :mod:`sp_validation.blinding` — Smokescreen-fork blinding wiring. +"""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, and the merge -alignment against a monkeypatched concealing factor) run in the default -suite; the theory tests (fork + CCL) are marked ``slow``; the COSEBIs / -pure-E/B tests additionally ``importorskip`` ``cosmo_numba``. - -All fixtures are synthetic and deterministic. The fork's ``ConcealDataVector`` -carries no data-vector consistency check (only the length guard), so fixture -ξ± values are smooth synthetic templates — no theory fill is needed to blind. +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 (coarse ξ±, fine ξ±, 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) fine +ξ± through the pipeline seams ``b_modes.cosebis_from_xi`` / +``b_modes.pure_eb_from_xi``, exactly as the pipeline does. The fork's +``ConcealDataVector`` carries no data-vector consistency check (only the +length guard), so fixture ξ± values are smooth synthetic templates — no +theory fill is needed to blind. """ import json @@ -26,7 +33,7 @@ # --------------------------------------------------------------------------- # -# Synthetic fixtures +# Synthetic part fixtures # --------------------------------------------------------------------------- # def _gauss_nz(z0, sigma, n=200): z = np.linspace(0.0, 3.0, n) @@ -57,164 +64,135 @@ def _b_mode_template(theta, amplitude): return amplitude * np.exp(-((np.log(np.asarray(theta) / 30.0)) ** 2) / 2.0) -def make_master( - nbins=1, - with_cl=True, - with_fine=True, - with_cosebis=False, - with_pure_eb=False, - with_rho=True, - b_amplitude=0.0, - nmodes=6, -): - """Build a master SACC in the sacc_io layout, synthetic values throughout. - - COSEBIs / pure-E/B blocks (when requested) are seeded with the *pipeline - estimators run on the fixture's own fine ξ±* — so a zero-shift blind must - reproduce them exactly (AC1) — and therefore require ``cosmo_numba``. - ``b_amplitude`` injects a pure B-mode into the fine ξ± (+ξ_B to ξ+, −ξ_B - to ξ−; b_modes.py sign convention). - """ - nz = {i: _gauss_nz(0.5 + 0.3 * i, 0.15 + 0.02 * i) for i in range(nbins)} - ctheta, ftheta = _coarse_theta(), _fine_theta() - pairs = [(i, j) for i in range(nbins) for j in range(i, nbins)] +def _nz_dict(nbins): + return {i: _gauss_nz(0.5 + 0.3 * i, 0.15 + 0.02 * i) for i in range(nbins)} - s = sio.new_sacc(nz, metadata={"catalogue_version": "vTEST"}) - blocks = [] - fine_xi = {} - for k, (i, j) in enumerate(pairs): - xip, xim = _xi_template(ftheta, k) - xi_b = _b_mode_template(ftheta, b_amplitude) - fine_xi[(i, j)] = (xip + xi_b, xim - xi_b) +def _pairs(nbins): + return [(i, j) for i in range(nbins) for j in range(i, nbins)] + - for k, (i, j) in enumerate(pairs): - xip, xim = _xi_template(ctheta, k) - sio.add_xi(s, (i, j), ctheta, xip, xim, grid="coarse") +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 fine part for AC4/AC9. + """ + theta = _coarse_theta() if grid == "coarse" else _fine_theta() + s = sio.new_sacc(_nz_dict(nbins), metadata={"catalogue_version": "vTEST"}) + 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="coarse"), - s.indices(sio.XI_MINUS, tr, grid="coarse"), + 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 - if with_cl: - 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() - for k, (i, j) in enumerate(pairs): - 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)) - - if with_cosebis: - for i, j in pairs: - xip_f, xim_f = fine_xi[(i, j)] - En, Bn = bd.rederive_cosebis( - ftheta, xip_f, xim_f, nmodes, scale_cut=(ftheta.min(), ftheta.max()) - ) - sio.add_cosebis(s, (i, j), En, Bn, (ftheta.min(), ftheta.max())) - tr = sio._pair((i, j)) - idx = np.concatenate( - [s.indices(sio.COSEBI_EE, tr), s.indices(sio.COSEBI_BB, tr)] - ) - blocks.append((idx, np.eye(len(idx)) * 1e-20)) - - if with_pure_eb: - # Seed the pure-E/B blocks through EXACTLY the pipeline convention the - # re-derivation reconstructs (blinding.rederive_eb_from_fine_xi): the - # reporting grid IS the coarse ξ± grid, the reporting ξ± are the - # MEASURED coarse ξ± (not an interpolation of the fine grid), and the - # integration ξ± are the fine grid. tmin/tmax are the edge-based - # integration bounds; the outermost reporting point sits at tmax with - # no interior support, so it comes back NaN — the AC9 boundary case. - # Seeding through this convention (rather than the blinding module's - # own interp/pseudo-edge seam) is what makes AC1/AC4/AC9 a real check - # of the re-derivation against the pipeline convention, not a tautology. - tmin, tmax = float(ctheta[0]), float(ctheta[-1]) - for k, (i, j) in enumerate(pairs): - xip_f, xim_f = fine_xi[(i, j)] - xip_r, xim_r = _xi_template(ctheta, k) # measured coarse reporting ξ± - modes = bd.rederive_pure_eb( - ctheta, xip_r, xim_r, ftheta, xip_f, xim_f, tmin, tmax - ) - sio.add_pure_eb( - s, - (i, j), - ctheta, - *(modes[key] for key in sio.PURE_KEYS), - bounds=(tmin, tmax), - ) - tr = sio._pair((i, j)) - idx = np.concatenate( - [s.indices(sio.PURE_TYPES[k], tr) for k in sio.PURE_KEYS] - ) - blocks.append((idx, np.eye(len(idx)) * 1e-14)) - if with_rho: - 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( +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"}) + 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, - (0,), - 0, - ctheta, - np.arange(len(ctheta)) * 3e-7, - np.arange(len(ctheta)) * 4e-7, + (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(sio.TAU_PLUS.format(k=0)), s.indices(sio.TAU_MINUS.format(k=0))] + [s.indices(dt, tr) for dt in (sio.CL_EE, sio.CL_BB, sio.CL_EB)] ) - blocks.append((idx, np.eye(len(idx)) * 1e-18)) + blocks.append((idx, np.eye(len(idx)) * 1e-16)) + sio.assemble_covariance(s, blocks) + return s - if with_fine: - for i, j in pairs: - xip_f, xim_f = fine_xi[(i, j)] - sio.add_xi(s, (i, j), ftheta, xip_f, xim_f, grid="fine") - tr = sio._pair((i, j)) - idx = np.concatenate( - [ - s.indices(sio.XI_PLUS, tr, grid="fine"), - s.indices(sio.XI_MINUS, tr, grid="fine"), - ] - ) - blocks.append((idx, np.eye(len(idx)) * 1e-14)) +def make_rho_part(): + """A standalone ρ/τ PSF-diagnostics part SACC — never blindable.""" + ctheta = _coarse_theta() + s = sio.new_sacc(_nz_dict(1), metadata={"catalogue_version": "vTEST"}) + 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_coarse": make_xi_part("coarse", nbins), + "xi_fine": make_xi_part("fine", nbins, b_amplitude=b_amplitude), + "cl": make_cl_part(nbins), + } + if with_rho: + parts["rho_tau"] = make_rho_part() + return parts + + +def _derive_downstream(coarse_part, fine_part, nmodes=6): + """COSEBIs + pure-E/B the way the pipeline derives them downstream. + + COSEBIs from the fine ξ± (full-range scale cut); pure-E/B from the + measured coarse reporting ξ± + the fine integration ξ±, with the + edge-based bounds set to the coarse 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(fine_part, (0, 0), grid="fine") + theta_c, xip_c, xim_c = sio.get_xi(coarse_part, (0, 0), grid="coarse") + 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) # --------------------------------------------------------------------------- # @@ -329,10 +307,10 @@ def test_hidden_params_no_global_rng_state(): # --------------------------------------------------------------------------- # -# Merge + provenance with a monkeypatched factor (fast — no CCL) +# Per-part merge + provenance with a monkeypatched factor (fast — no CCL) # --------------------------------------------------------------------------- # def _patch_constant_factor(monkeypatch, value=1e-6): - def fake(master, indices, factory, config, seed): + def fake(part, indices, factory, config, seed): return np.arange(len(indices), dtype=float) * value + value monkeypatch.setattr(bd, "_concealing_factor", fake) @@ -340,35 +318,37 @@ def fake(master, indices, factory, config, seed): def test_merge_places_shift_at_recorded_indices_only(monkeypatch): - """AC5 (merge half): the shift lands exactly on the blindable rows, - in master order; ρ/τ, covariance, n(z), and every tag are untouched.""" + """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) - s = make_master(nbins=2, with_cl=True, with_fine=True, with_rho=True) - orig = np.array(s.mean) - orig_cov = s.covariance.dense.copy() - orig_nz = s.tracers["source_0"].nz.copy() - - blinded = bd.blind_sacc(s, "seed", log=_NOLOG) - - shifted = np.zeros(len(orig), dtype=bool) - for _, indices, _ in bd._blindable_blocks(s): - 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(s.data, blinded.data): - assert a.data_type == b.data_type - assert a.tracers == b.tracers - assert a.tags == b.tags + 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 files carry concealed/blind/commitment/digest; no seed key.""" + """Blinded parts carry concealed/blind/commitment/digest; no seed key.""" _patch_constant_factor(monkeypatch) - s = make_master(with_cl=False) + s = make_xi_part("coarse") s.metadata["seed_smokescreen"] = "leaked!" # must be stripped c = bd.BlindingConfig() blinded = bd.blind_sacc(s, "seed", config=c, label="B", log=_NOLOG) @@ -382,50 +362,22 @@ def test_provenance_metadata_contract(monkeypatch): def test_blind_refuses_double_blind(monkeypatch): _patch_constant_factor(monkeypatch) - s = make_master(with_cl=False) + s = make_xi_part("coarse") blinded = bd.blind_sacc(s, "seed", log=_NOLOG) with pytest.raises(ValueError, match="already concealed"): bd.blind_sacc(blinded, "seed2", log=_NOLOG) -def test_guard_requires_coarse_xi(): - nz = {0: _gauss_nz(0.7, 0.2)} - s = sio.new_sacc(nz) - with pytest.raises(ValueError, match="no coarse"): - bd.blind_sacc(s, "seed", log=_NOLOG) - - -def test_guard_requires_both_xi_plus_and_minus(): - nz = {0: _gauss_nz(0.7, 0.2)} - s = sio.new_sacc(nz) - for th in _coarse_theta(): - s.add_data_point( - sio.XI_PLUS, - ("source_0", "source_0"), - 1e-5, - theta=float(th), - grid="coarse", - ) - with pytest.raises(ValueError, match="no coarse"): - bd.blind_sacc(s, "seed", log=_NOLOG) - - -def test_guard_derived_blocks_require_fine_xi(): - """COSEBIs present but no fine ξ± ⇒ loud refusal (cannot re-derive).""" - nz = {0: _gauss_nz(0.7, 0.2)} - s = sio.new_sacc(nz) - ctheta = _coarse_theta() - xip, xim = _xi_template(ctheta) - sio.add_xi(s, (0, 0), ctheta, xip, xim, grid="coarse") - sio.add_cosebis(s, (0, 0), np.ones(4) * 1e-10, np.ones(4) * 1e-12, (1.0, 100.0)) - with pytest.raises(ValueError, match="no grid='fine'"): - bd.blind_sacc(s, "seed", 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_master(with_cl=False) + s = make_xi_part("coarse") blinded = bd.blind_sacc(s, "right-seed", log=_NOLOG) with pytest.raises(ValueError, match="blind_commitment"): bd.unblind_sacc(blinded, "wrong-seed", log=_NOLOG) @@ -440,8 +392,125 @@ def test_unblind_fails_closed_on_wrong_seed_or_config(monkeypatch): bd.unblind_sacc(s, "right-seed", log=_NOLOG) -def test_cli_blind_refuses_existing_output(tmp_path): - """The blind CLI refuses to overwrite a previous blind (before compute).""" +# --------------------------------------------------------------------------- # +# 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_sha256", "config_digest"} + assert len(commitment["seed_sha256"]) == 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_sha256"] + # 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_coarse", "xi_fine", "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(), + } + + +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_coarse"], seed="one") + _stamp(parts["xi_fine"], 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_coarse"]) # 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 pre-blinding / mock assembly (nothing blinded) asserts nothing.""" + assert bd.assert_consistent_blind(list(make_parts().values())) is None + + +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_coarse", "xi_fine", "cl"): + _stamp(parts[name]) + ordered = [parts[k] for k in ("xi_coarse", "cl", "rho_tau", "xi_fine")] + s = sio.gather(ordered, metadata={"catalogue_version": "vTEST"}) + + 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 fine rows are addressable by the grid tag in the assembled file + assert len(s.indices(sio.XI_PLUS, grid="fine")) == len( + parts["xi_fine"].indices(sio.XI_PLUS, grid="fine") + ) + # 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_coarse"], seed="one") + _stamp(parts["xi_fine"], 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 = ( @@ -451,47 +520,48 @@ def test_cli_blind_refuses_existing_output(tmp_path): cli = importlib.util.module_from_spec(spec) spec.loader.exec_module(cli) - (tmp_path / "master_blinded_A.sacc").write_text("stale") + (tmp_path / "commitment.json").write_text("{}") with pytest.raises(SystemExit, match="refusing to overwrite"): - cli.main(["blind", "master.sacc", "-o", str(tmp_path), "--label", "A"]) + cli.main(["blind-init", str(tmp_path)]) # --------------------------------------------------------------------------- # # AC2 + AC3 + AC7: the shift itself (slow — fork + CCL) # --------------------------------------------------------------------------- # @pytest.mark.slow -def test_ac2_on_file_shift_equals_theory_difference(): - """AC2: per-row shift on the blinded file == theory_fn(hidden) − +def test_ac2_on_file_shift_equals_theory_difference_per_part(): + """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 blocks, on a two-bin fixture. + 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 + merge placement (the shift on - the file is exactly what re-running the same backend at the recovered + 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`` via - :func:`bd._blindable_blocks`, so any wrong-cosmology dependence cancels - and a wrong backend would still pass here. Backend correctness is carried - by AC3 (:func:`test_ac3_cross_backend_against_independent_ccl_reference`), - which builds an independent ``ccl.Cosmology`` reference.""" - master = make_master(nbins=2, with_cl=True, with_fine=True) + 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.""" cfg = bd.BlindingConfig() seed = "ac2-seed" - blinded = bd.blind_sacc(master, seed, config=cfg, log=_NOLOG) - - hidden = bd.hidden_params(seed, cfg) - fiducial = cfg.theory.ccl_params() - worst = 0.0 - blocks = bd._blindable_blocks(master) - assert {name for name, _, _ in blocks} == {"coarse ξ±", "fine ξ±", "pseudo-Cℓ_EE"} - for name, indices, factory in blocks: - theory = factory(bd._extract_block(master, indices), cfg.theory) + 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) + theory = factory(bd._extract_block(part, indices), cfg.theory) expected = theory(hidden) - theory(fiducial) - actual = np.array(blinded.mean)[indices] - np.array(master.mean)[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}: |Δ|={gap:.3e}" - print(f"\nAC2 max relative shift mismatch across blocks: {worst:.3e}") + 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 @@ -500,10 +570,12 @@ def test_ac3_cross_backend_against_independent_ccl_reference(): reference (self-contained here; does not touch the blinding backends).""" import pyccl as ccl - master = make_master(nbins=2, with_cl=True, with_fine=False) cfg = bd.BlindingConfig() seed = "ac3-seed" - blinded = bd.blind_sacc(master, seed, config=cfg, log=_NOLOG) + coarse = make_xi_part("coarse", nbins=2) + cl_part = make_cl_part(nbins=2) + blinded_xi = bd.blind_sacc(coarse, 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() @@ -531,30 +603,30 @@ def ref_xi(params, nz_i, nz_j, theta): return xip, xim worst = 0.0 - for i, j in bd.xi_pairs(master, "coarse"): + for i, j in bd.xi_pairs(coarse, "coarse"): tr = sio._pair((i, j)) - theta = sio._tag(master, sio.XI_PLUS, tr, "theta", grid="coarse") - nz_i, nz_j = sio.get_nz(master, i), sio.get_nz(master, j) + theta = sio._tag(coarse, sio.XI_PLUS, tr, "theta", grid="coarse") + nz_i, nz_j = sio.get_nz(coarse, i), sio.get_nz(coarse, 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 = master.indices(dt, tr, grid="coarse") - realized = np.array(blinded.mean)[idx] - np.array(master.mean)[idx] + idx = coarse.indices(dt, tr, grid="coarse") + realized = np.array(blinded_xi.mean)[idx] - np.array(coarse.mean)[idx] worst = max(worst, np.max(np.abs(realized - ref_shift))) print(f"\nAC3 max |realized − independent reference| (coarse ξ±): {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(master): + for i, j in bd.cl_pairs(cl_part): tr = sio._pair((i, j)) - idx = master.indices(sio.CL_EE, tr) - window = master.get_bandpower_windows(idx) + 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(master, i), sio.get_nz(master, j) + nz_i, nz_j = sio.get_nz(cl_part, i), sio.get_nz(cl_part, j) def ref_cl(params): cosmo = ref_cosmo(params) @@ -563,84 +635,87 @@ def ref_cl(params): return ccl.angular_cl(cosmo, ti, tj, w_ell) ref_shift = w_mat.T @ (ref_cl(hidden) - ref_cl(fiducial)) - realized = np.array(blinded.mean)[idx] - np.array(master.mean)[idx] + 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 with the same (seed, config) produce identical - shifts; a different seed produces a different one.""" - master = make_master(nbins=1, with_cl=True, with_fine=True) - b1 = bd.blind_sacc(master, "repro-seed", log=_NOLOG) - b2 = bd.blind_sacc(master, "repro-seed", log=_NOLOG) + """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("coarse") + 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(master, "other-seed", log=_NOLOG) + 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: derived statistics (slow + cosmo_numba) +# 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 ξ±, COSEBIs, and pure-E/B exactly — - the extract/merge/re-derivation plumbing is the identity at zero shift.""" + """AC1: a zero envelope reproduces every part exactly, and the fine 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") - master = make_master( - with_cl=True, with_fine=True, with_cosebis=True, with_pure_eb=True - ) zero = bd.BlindingConfig(s8_half_width=0.0, omega_m_half_width=0.0) - blinded = bd.blind_sacc(master, "any-seed", config=zero, log=_NOLOG) - orig, new = np.array(master.mean), np.array(blinded.mean) - both_nan = np.isnan(orig) & np.isnan(new) # degenerate pure-EB boundary pts - assert np.array_equal(orig[~both_nan], new[~both_nan]), ( - "zero-shift blind changed values" - ) - assert np.array_equal(np.isnan(orig), np.isnan(new)) + 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_coarse"], parts["xi_fine"]) + En_b, Bn_b, modes_b = _derive_downstream(blinded["xi_coarse"], blinded["xi_fine"]) + 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 blinding 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.""" + """AC4: the ΔBₙ induced by deriving B-modes from the blinded fine ξ± 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): - master = make_master( - with_cl=False, - with_fine=True, - with_cosebis=True, - with_pure_eb=True, - b_amplitude=amp, - ) - blinded = bd.blind_sacc(master, seed, log=_NOLOG) - tr = sio._pair((0, 0)) - idx_b = master.indices(sio.COSEBI_BB, tr) - d_bn = np.array(blinded.mean)[idx_b] - np.array(master.mean)[idx_b] - bn0 = np.array(master.mean)[idx_b] - # pure-EB B blocks: fixed absolute leakage too - idx_pb = master.indices(sio.PURE_TYPES["xip_B"], tr) - d_xib = np.array(blinded.mean)[idx_pb] - np.array(master.mean)[idx_pb] + coarse = make_xi_part("coarse") + fine = make_xi_part("fine", b_amplitude=amp) + blinded_coarse = bd.blind_sacc(coarse, seed, log=_NOLOG) + blinded_fine = bd.blind_sacc(fine, seed, log=_NOLOG) + _, Bn_t, modes_t = _derive_downstream(coarse, fine) + _, Bn_b, modes_b = _derive_downstream(blinded_coarse, blinded_fine) + 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(bn0)):.2e}), " + 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 - # identical to ~1e-10 relative (observed 1.1e-10 — quadrature float noise - # on the injected-B term; the bound carries a decade of margin) 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)| = {gap:.3e} ({gap / scale:.2e} of |ΔBₙ|)" + 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" @@ -648,169 +723,196 @@ def test_ac4_b_mode_invariance_and_leakage_floor(): # 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 ~1e-2 of the (tiny) leakage itself — observed - # 2.0e-11 absolute on a 1.5e-9 leakage, i.e. 1e-5 of the smaller - # injected B. The COSEBIs assertion above carries the exact-identity - # criterion; this one bounds the quadrature wobble. + # 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} ({gap_x / scale_x:.2e} of the leakage)" + 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: covariance and ρ/τ byte-identical; every shifted block's rows land - at their original master indices (order-preservation).""" - pytest.importorskip("cosmo_numba") - master = make_master( - with_cl=True, - with_fine=True, - with_cosebis=True, - with_pure_eb=True, - with_rho=True, - ) - blinded = bd.blind_sacc(master, "ac5-seed", log=_NOLOG) - - assert np.array_equal(blinded.covariance.dense, master.covariance.dense) - for dt in ( - sio.RHO_PLUS.format(k=0), - sio.RHO_MINUS.format(k=0), - sio.RHO_PLUS.format(k=1), - sio.RHO_MINUS.format(k=1), - sio.TAU_PLUS.format(k=0), - sio.TAU_MINUS.format(k=0), - sio.CL_BB, - sio.CL_EB, - ): - idx = master.indices(dt) + """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_coarse", "xi_fine", "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.mean)[idx], np.array(master.mean)[idx] + np.array(blinded["cl"].mean)[idx], np.array(parts["cl"].mean)[idx] ), f"{dt} was touched by the blind" - # row order: the identity of every row (type/tracers/tags) is unchanged - for a, b in zip(master.data, blinded.data): - assert (a.data_type, a.tracers, a.tags) == (b.data_type, b.tracers, b.tags) - # ξ± rows did move (the blind actually blinded) - xi_idx = bd._xi_indices(master, "coarse") - assert not np.allclose( - np.array(blinded.mean)[xi_idx], np.array(master.mean)[xi_idx], atol=0 + # the blindable rows did move (the blind actually blinded) + for name, grid in (("xi_coarse", "coarse"), ("xi_fine", "fine")): + 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_coarse"], parts["rho_tau"], blinded["cl"], blinded["xi_fine"]] + ) + 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 re-derived pure-E/B NaN pattern is identical on the true and - blinded files — blinding never moves a NaN. + """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 file (on production files it is empty — the - reporting grid is a strict sub-range of the integration range, so the - degeneracy never fires; this is what the corrected re-derivation, feeding - the measured coarse reporting ξ± and the fine integration grid, delivers), - the blinded-file re-derivation must reproduce it bit-for-bit: the blind is - a pure shift of the inputs, not a change of estimator support. The finite - values move (the ξ± shifted); the NaN mask does not.""" + 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") - master = make_master(with_cl=False, with_fine=True, with_pure_eb=True) - tr = sio._pair((0, 0)) - idx_e = master.indices(sio.PURE_TYPES["xip_E"], tr) - true_vals = np.array(master.mean)[idx_e] - - blinded = bd.blind_sacc(master, "ac9-seed", log=_NOLOG) - blind_vals = np.array(blinded.mean)[idx_e] - assert np.array_equal(np.isnan(true_vals), np.isnan(blind_vals)), ( - "blinding moved the pure-E/B NaN pattern — the estimator support changed" + seed = "ac9-seed" + coarse, fine = make_xi_part("coarse"), make_xi_part("fine") + _, _, modes_t = _derive_downstream(coarse, fine) + _, _, modes_b = _derive_downstream( + bd.blind_sacc(coarse, seed, log=_NOLOG), bd.blind_sacc(fine, 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 ξ±) - finite = np.isfinite(true_vals) - assert finite.any() and not np.allclose( - true_vals[finite], blind_vals[finite], atol=0 - ) + 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 CLIs (slow + cosmo_numba) +# AC6 + AC8: end-to-end custody through the file surface (slow + cosmo_numba) # --------------------------------------------------------------------------- # @pytest.mark.slow -def test_ac6_ac8_end_to_end_blind_unblind_custody(tmp_path): - """AC8: the blind CLI blinds a master file end-to-end and the unblind CLI - restores it bit-for-bit. AC6: no plaintext seed anywhere on disk, no - ``seed_smokescreen`` key; unblind fails closed on a tampered commitment.""" +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 + fine 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") - master = make_master( - with_cl=True, - with_fine=True, - with_cosebis=True, - with_pure_eb=True, - with_rho=True, + parts = make_parts(nbins=1) + En_true, Bn_true, modes_true = _derive_downstream( + parts["xi_coarse"], parts["xi_fine"] ) - master_path = tmp_path / "vTEST.fits" - sio.save(master, str(master_path)) - - out = tmp_path / "blinded" - out.mkdir() - paths = bd.blind_file(str(master_path), str(out), label="A", log=_NOLOG) - - # -- custody hygiene (AC6) ------------------------------------------- -- - # The plaintext input master is deleted at blinding time — the true vector - # survives only inside the encrypted bundle. - assert not master_path.exists(), "plaintext input master was not deleted" - blinded = sio.load(paths["blinded"]) - assert blinded.metadata["concealed"] is True - assert "seed_smokescreen" not in blinded.metadata - with open(paths["commitment"], encoding="utf-8") as f: + + 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_coarse", "xi_fine", "cl"): + path = tmp_path / f"{name}.fits" + sio.save(parts[name], str(path)) + 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 blinded.metadata["blind_commitment"] == commitment["seed_sha256"] - assert not (out / "blind_A_bundle.json").exists() # plaintext deleted - # only the four custody outputs exist; the commitment carries hashes only - assert {p.name for p in out.iterdir()} == { - pathlib.Path(v).name for v in paths.values() - } assert set(commitment) == {"label", "seed_sha256", "config_digest"} - # the blinded vector actually differs from the truth - assert not np.array_equal(np.array(blinded.mean), np.array(master.mean)) + 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_sha256"] + # 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_coarse"], + blinded_parts["cl"], + parts["rho_tau"], + blinded_parts["xi_fine"], + ], + metadata={"catalogue_version": "vTEST"}, + ) + assert assembled.metadata["blind_commitment"] == commitment["seed_sha256"] + # born-blinded derived statistics from the blinded parts differ from truth + En_b, Bn_b, _ = _derive_downstream( + blinded_parts["xi_coarse"], blinded_parts["xi_fine"] + ) + assert not np.allclose(En_b, En_true, atol=0) - # -- fail-closed on tampered commitment (AC6) ------------------------- -- + # -- fail-closed on tampered commitment (AC6) --------------------------- -- tampered = dict(commitment, seed_sha256="0" * 64) - tampered_path = tmp_path / "tampered.json" - tampered_path.write_text(json.dumps(tampered)) + with open(init["commitment"], "w", encoding="utf-8") as f: + json.dump(tampered, f) with pytest.raises(ValueError, match="sha256"): - bd.unblind_file( - paths["blinded"], - paths["bundle"], - paths["key"], - str(tampered_path), + bd.unblind_part( + out_paths["xi_fine"]["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_file( - paths["blinded"], - paths["bundle"], - paths["key"], - paths["commitment"], + bd.unblind_part( + out_paths["xi_fine"]["blinded"], + str(blind_dir), str(tmp_path / "never.fits"), config=bd.BlindingConfig(s8_half_width=0.01), log=_NOLOG, ) - # -- bit-for-bit restoration (AC8) ------------------------------------ -- - restored_path = tmp_path / "restored.fits" - bd.unblind_file( - paths["blinded"], - paths["bundle"], - paths["key"], - paths["commitment"], - str(restored_path), - log=_NOLOG, - ) - restored = sio.load(str(restored_path)) - orig, back = np.array(master.mean), np.array(restored.mean) - both_nan = np.isnan(orig) & np.isnan(back) - assert np.array_equal(orig[~both_nan], back[~both_nan]) # bit-for-bit - assert np.array_equal(np.isnan(orig), np.isnan(back)) - assert not restored.metadata.get("concealed", False) - assert "blind_commitment" not in restored.metadata + # -- bit-for-bit restoration per part (AC8) ----------------------------- -- + restored = {} + for name in ("xi_coarse", "xi_fine", "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_coarse"], restored["xi_fine"]) + 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 From 3bf39c830b63d8070e530c83fdbe94dfd733c835 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 14:34:18 +0200 Subject: [PATCH 15/17] blinding: escrow bundles land at the declared name for dotted part stems MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit smokescreen's encrypt_file names its save-file outputs from basename.split('.')[0] — truncating at the first dot. For the canonical versioned part name (v1.4.6.3_xi_fine_escrow) that wrote v1.encrpt/v1.key instead of the *_escrow.encrpt/.key that part_paths() declares, so: - unblind_part's escrow branch never fired (declared path absent), silently falling back to subtraction-only — not bit-for-bit (AC #8); - two parts sharing a first-dot prefix collided onto one bundle, the second overwriting the first's escrowed truth — and since blind_part deletes the plaintext, that truth became unrecoverable, breaking the per-part custody guarantee (AC #6). _write_encrypted_json now drives encrypt_file with save_file=False, takes its returned (ciphertext, key), and writes them at the exact encrpt_path / .key names we control. Robust for any stem; hardens the seed bundle too. Plaintext is still deleted (keep_original=False runs unconditionally). New test_ac8_dotted_versioned_part_names_escrow_and_restore blinds two dot-prefix-sharing production names (v1.4.6.3_xi_coarse/xi_fine) into one dir and asserts distinct escrow bundles at the declared names, no v1.encrpt collision target, and bit-for-bit restore. It fails against the old code. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/blinding.py | 23 +++++++++--- src/sp_validation/tests/test_blinding.py | 48 ++++++++++++++++++++++++ 2 files changed, 65 insertions(+), 6 deletions(-) diff --git a/src/sp_validation/blinding.py b/src/sp_validation/blinding.py index 04d160ee..529f070b 100644 --- a/src/sp_validation/blinding.py +++ b/src/sp_validation/blinding.py @@ -707,18 +707,29 @@ def unblind_part(blinded_path, blind_dir, out_path, config=None, log=print): def _write_encrypted_json(encrpt_path, payload): - """Encrypt ``payload`` (JSON) to ``encrpt_path``; no plaintext survives. - - ``smokescreen.encryption.encrypt_file`` (Fernet) writes ``.encrpt`` - + ``.key`` side by side and removes the plaintext - (``keep_original=False``). + """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_fine_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) - encrypt_file(plaintext, save_file=True, keep_original=False) + 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): diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py index 48d86460..a177b55b 100644 --- a/src/sp_validation/tests/test_blinding.py +++ b/src/sp_validation/tests/test_blinding.py @@ -916,3 +916,51 @@ def test_ac6_ac8_end_to_end_init_parts_gather_unblind(tmp_path): 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_coarse.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_coarse", "xi_fine"): + path = tmp_path / f"{version}_{name}.fits" + sio.save(parts[name], str(path)) + 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_coarse"] != escrow_files["xi_fine"] + 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_coarse", "xi_fine"): + 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" + ) From a6dd5c06d1d4e0798587153d2124eeb8a42f2260 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 14:49:16 +0200 Subject: [PATCH 16/17] blinding: make config digest canonical across int-vs-float literals MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Config values arrive from YAML/CLI/humans, so a declared-float field can be written -1 (int) or -1.0 (float) for one physical cosmology. The digest was json.dumps of the raw field values, so those two literals produced different sha256 digests (e0b60d7f… vs df6be45e…). An int-vs-float mismatch between blind_init and a later blind_part/unblind then made _read_seed and unblind_sacc both raise "config digest mismatch" and refuse a correctly blinded file — a denial-of-unblind that falsified the PRD's byte-identical / round-trip-exact digest guarantee. Normalize at the digest boundary: every float-declared TheoryConfig field (and BlindingConfig's two half-widths) is coerced through float() before serialization, so the digest depends on the numeric value, not the literal's Python type. Coerce in TheoryConfig/BlindingConfig.from_overrides too, per the field's declared type (strings pass through untouched), so a value is canonical the moment it enters a config. Regression test asserts digest(w0=-1) == digest(w0=-1.0), same for Omega_m/m_nu/wa/S8 and the half-widths, plus an all-int round-trip against the float default. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/blinding.py | 36 ++++++++++++++------- src/sp_validation/blinding_theory.py | 19 ++++++++--- src/sp_validation/tests/test_blinding.py | 40 ++++++++++++++++++++++++ 3 files changed, 78 insertions(+), 17 deletions(-) diff --git a/src/sp_validation/blinding.py b/src/sp_validation/blinding.py index 529f070b..9641918f 100644 --- a/src/sp_validation/blinding.py +++ b/src/sp_validation/blinding.py @@ -102,17 +102,26 @@ def config_digest(self): 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. Python's ``json`` emits floats via their - shortest round-trip ``repr``, so two runs of the same config produce - byte-identical digests. Checked (with ``sha256(seed)``) at unblind, - so a wrong envelope or a mismatched P(k) recipe cannot silently - subtract a wrong shift. + 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 + ``sha256(seed)``) at unblind, so a wrong envelope or a mismatched + P(k) recipe cannot silently subtract a wrong shift. """ payload = { - "s8_half_width": self.s8_half_width, - "omega_m_half_width": self.omega_m_half_width, + "s8_half_width": float(self.s8_half_width), + "omega_m_half_width": float(self.omega_m_half_width), "theory": { - f.name: getattr(self.theory, f.name) + 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) }, } @@ -127,14 +136,17 @@ def from_overrides(cls, overrides): ``theory`` may be given as a :class:`TheoryConfig` or a mapping of TheoryConfig overrides, mirroring :meth:`TheoryConfig.from_overrides`. """ - fields = {f.name for f in dataclasses.fields(cls)} - unknown = set(overrides) - fields + 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(fields)}" + f"valid fields are {sorted(by_name)}" ) - overrides = dict(overrides) + 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)) diff --git a/src/sp_validation/blinding_theory.py b/src/sp_validation/blinding_theory.py index 47e67ca3..dbec0084 100644 --- a/src/sp_validation/blinding_theory.py +++ b/src/sp_validation/blinding_theory.py @@ -159,15 +159,24 @@ def ccl_params(self): @classmethod def from_overrides(cls, overrides): - """Build from a mapping of field overrides (fail loud on unknown keys).""" - fields = {f.name for f in dataclasses.fields(cls)} - unknown = set(overrides) - fields + """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(fields)}" + f"valid fields are {sorted(by_name)}" ) - return cls(**overrides) + coerced = { + name: (float(v) if by_name[name].type in (float, "float") else v) + for name, v in overrides.items() + } + return cls(**coerced) # --------------------------------------------------------------------------- # diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py index a177b55b..a0fff6ab 100644 --- a/src/sp_validation/tests/test_blinding.py +++ b/src/sp_validation/tests/test_blinding.py @@ -249,6 +249,46 @@ def test_config_digest_stable_and_sensitive(): ) +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() From 4e53c76d4b7fa44b96cca6b3eed8e8e895bf24ba Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 14:49:37 +0200 Subject: [PATCH 17/17] blinding: seed math is the custody authority; label is provenance not state MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three custody-hardening fixes to unblind_part and the assembly assertion, keeping the seed+commitment as the single root of trust: - unblind_part made the escrow authoritative on the success branch: on residual<=1e-6 it overwrote the seed-subtracted vector with the escrow's true_mean, so the escrow (not the seed math) decided the output whenever it happened to agree. Now the seed-subtracted vector is authoritative; the escrow is used only as a tighter equality CHECK, and its stored seed_sha256 must match the commitment before it is trusted at all (an escrow from a different blind can no longer determine the output). The ulp-residue overwrite that AC8's bit-for-bit restore relies on is retained, but only once the escrow is cryptographically bound to this commitment. - The escrow path was derived by blinded_path.replace("_blinded", ""), a global substring replace: a directory containing "_blinded" (e.g. /runs/2024_blinded/...) corrupted the derived escrow path and silently dropped the bit-for-bit cross-check. Anchor the replace to the basename. - assert_consistent_blind keyed its consistency set on (label, commitment, digest), so two parts with the same seed+config but different --label falsely raised "different blind commitments" and blocked a legitimate assembly. Key on (commitment, digest) only; the label is informational provenance — differing labels now emit a distinct warning and the assembly proceeds, stamping the sorted-first label. Still fails closed on differing commitment or digest. Regression test: same seed+config, different labels assemble with a UserWarning rather than a failure. AC8 bit-for-bit restore (including the dotted-versioned-part escrow test) still passes. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/sp_validation/blinding.py | 71 ++++++++++++++++++------ src/sp_validation/tests/test_blinding.py | 18 ++++++ 2 files changed, 71 insertions(+), 18 deletions(-) diff --git a/src/sp_validation/blinding.py b/src/sp_validation/blinding.py index 9641918f..6ef506d3 100644 --- a/src/sp_validation/blinding.py +++ b/src/sp_validation/blinding.py @@ -55,6 +55,7 @@ import json import os import secrets +import warnings import numpy as np @@ -484,8 +485,12 @@ def assert_consistent_blind(parts): ``ValueError`` — if blinded and plaintext blindable parts are mixed, or if two parts carry different ``blind_commitment``/``blind_config_digest`` (they were blinded under different seeds or configs and must never be - combined). If no part is blinded (a pre-blinding or mock assembly), - returns ``None``. + combined). The consistency key is ``(commitment, digest)`` only — 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). + If no part is blinded (a pre-blinding or mock assembly), returns + ``None``. Returns ------- @@ -504,24 +509,29 @@ def assert_consistent_blind(parts): f"({len(concealed)} of {len(blindable)} blinded) — refusing to " "combine (a plaintext part beside blinded ones leaks the shift)" ) + # Custody state is (commitment, digest) only — the label is provenance. stamps = { - ( - p.metadata["blind"], - p.metadata["blind_commitment"], - p.metadata["blind_config_digest"], - ) + (p.metadata["blind_commitment"], p.metadata["blind_config_digest"]) for p in concealed } if len(stamps) != 1: raise ValueError( "parts carry different blind commitments — they were blinded " "under different seeds or configs and must never be combined: " - + "; ".join(f"({label}, {c[:12]}…, {d[:12]}…)" for label, c, d in stamps) + + "; ".join(f"({c[:12]}…, {d[:12]}…)" for c, d in stamps) + ) + ((commitment, digest),) = stamps + 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}." ) - ((label, commitment, digest),) = stamps return { "concealed": True, - "blind": label, + "blind": labels[0], "blind_commitment": commitment, "blind_config_digest": digest, } @@ -686,21 +696,43 @@ def unblind_part(blinded_path, blind_dir, out_path, config=None, log=print): Decrypts the seed bundle and verifies both ``sha256(seed)`` and the config digest against ``commitment.json`` (fail closed on either — verification precedes subtraction), recomputes the part's shift from the - seed and subtracts it (:func:`unblind_sacc`). When the part's escrow - bundle exists beside the blinded file, the subtraction is cross-checked - against the escrowed true vector and the truth restored bit-for-bit - (float add-then-subtract leaves ~ulp residue; the escrow removes it). - On the assembled file — no escrow — the subtraction stands alone, with - fine rows selected by the ``grid`` tag. + 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_sha256`` 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 fine rows selected by the + ``grid`` tag. """ config = config or BlindingConfig() - seed, _ = _read_seed(blind_dir, config) + seed, commitment = _read_seed(blind_dir, config) blinded = sacc_io.load(blinded_path) part = unblind_sacc(blinded, seed, config=config, log=log) - escrow = part_paths(blinded_path.replace("_blinded", "")) + 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_sha256") != commitment["seed_sha256"]: + 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( @@ -711,6 +743,9 @@ def unblind_part(blinded_path, blind_dir, out_path, config=None, log=print): 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})") sacc_io.save(part, out_path) diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py index a0fff6ab..1cbd4782 100644 --- a/src/sp_validation/tests/test_blinding.py +++ b/src/sp_validation/tests/test_blinding.py @@ -508,6 +508,24 @@ def test_assert_consistent_blind_all_plaintext_is_none(): assert bd.assert_consistent_blind(list(make_parts().values())) is None +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_coarse"], seed="one", label="A") + _stamp(parts["xi_fine"], 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."""