Smokescreen-fork blinding: per-part-at-birth blind, hash-commitment custody, CAMB↔CCL cross-check (PR 6)

.github/workflows/deploy-image.yml

@@ -46,16 +46,31 @@ 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.
+      # `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; 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; from cosmo_numba.B_modes.cosebis import COSEBIS; from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes"
 
-      # 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

pyproject.toml

@@ -113,12 +113,28 @@ 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). 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
@@ -129,7 +145,12 @@ 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",
+    # 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

scripts/blind_data_vector.py

@@ -0,0 +1,173 @@
+#!/usr/bin/env python3
+
+"""Script blind_data_vector.py
+
+Per-part data-vector blinding with :mod:`sp_validation.blinding`
+(Smokescreen-fork concealment, hash-commitment custody).
+
+``blind-init`` runs once per catalogue version: it draws an OS-entropy seed
+(never printed, never written in plaintext), publishes a repo-committable
+``commitment.json`` (``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
+--------
+Once per catalogue version::
+
+    blind_data_vector.py blind-init blinded/
+
+Per intermediate part, at birth::
+
+    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 parts/xi_fine_blinded.fits \\
+        blinded/commitment.json
+"""
+
+import argparse
+import json
+import pathlib
+import sys
+
+from sp_validation import blinding, sacc_io
+
+
+def _config_from_args(args):
+    """A :class:`blinding.BlindingConfig` from optional CLI overrides."""
+    overrides = {}
+    if args.s8_half_width is not None:
+        overrides["s8_half_width"] = args.s8_half_width
+    if args.omega_m_half_width is not None:
+        overrides["omega_m_half_width"] = args.omega_m_half_width
+    return blinding.BlindingConfig.from_overrides(overrides)
+
+
+def _blind_init(args):
+    config = _config_from_args(args)
+    blind_dir = pathlib.Path(args.blind_dir)
+    blind_dir.mkdir(parents=True, exist_ok=True)
+    # Refuse before drawing anything: a blind is a one-shot custody event and
+    # silently overwriting a previous blind's state would destroy the record
+    # tying that blind to its seed.
+    clashes = [
+        p
+        for p in blinding.init_paths(str(blind_dir)).values()
+        if pathlib.Path(p).exists()
+    ]
+    if clashes:
+        raise SystemExit(
+            "refusing to overwrite existing blind state:\n  "
+            + "\n  ".join(clashes)
+            + "\nPick a fresh --blind-dir (never overwrite a blind)."
+        )
+    blinding.blind_init(str(blind_dir), config=config, label=args.label)
+    print(
+        "Commit the commitment JSON to the repo; keep the bundle + key safe "
+        "and separated (colocation in the blind dir is not at-rest protection)."
+    )
+
+
+def _blind_part(args):
+    blinding.blind_part(
+        args.part,
+        args.blind_dir,
+        config=_config_from_args(args),
+        keep_input=args.keep_input,
+    )
+
+
+def _unblind(args):
+    blinding.unblind_part(
+        args.blinded,
+        args.blind_dir,
+        args.output,
+        config=_config_from_args(args),
+    )
+
+
+def _verify(args):
+    """Seedless check: blinded-file metadata ↔ commitment JSON."""
+    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-init", "blind-part", "unblind"):
+        p = sub.add_parser(name)
+        p.add_argument("--s8-half-width", type=float, default=None)
+        p.add_argument("--omega-m-half-width", type=float, default=None)
+        if name == "blind-init":
+            p.add_argument(
+                "blind_dir",
+                help="directory for the blind's fixed state (commitment + "
+                "encrypted seed bundle)",
+            )
+            p.add_argument("--label", default="A", help="blind label (default A)")
+            p.set_defaults(func=_blind_init)
+        elif name == "blind-part":
+            p.add_argument("part", help="intermediate part SACC file to blind")
+            p.add_argument(
+                "--blind-dir", required=True, help="blind-init state directory"
+            )
+            p.add_argument(
+                "--keep-input",
+                action="store_true",
+                help="retain the plaintext input part (default: delete it "
+                "after blinding — the true vector is escrowed beside the "
+                "blinded output)",
+            )
+            p.set_defaults(func=_blind_part)
+        else:
+            p.add_argument("blinded", help="blinded part (or assembled) SACC file")
+            p.add_argument(
+                "--blind-dir", required=True, help="blind-init state directory"
+            )
+            p.add_argument("-o", "--output", required=True, help="output SACC path")
+            p.set_defaults(func=_unblind)
+
+    p = sub.add_parser("verify")
+    p.add_argument("blinded", help="blinded SACC file")
+    p.add_argument("commitment", help="commitment JSON")
+    p.set_defaults(func=_verify)
+
+    args = parser.parse_args(argv)
+    args.func(args)
+
+
+if __name__ == "__main__":
+    sys.exit(main())

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.