From d31a6f628c341cc98fa180522abd817980627a18 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 01:49:49 +0200 Subject: [PATCH 1/6] deps: py3.12 floor + SACC/blinding stack (sacc, firecrown, smokescreen) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit PRD #241 row 1. sacc>=2.4 joins core dependencies (sp_validation.sacc_io will be core library code). firecrown v1.15.1 + smokescreen 1.5.6 + an exact pyccl pin form the new [blinding] extra: the blind must be exactly recomputable from the seed at unblinding time, so the theory stack is pinned as a set. firecrown is not on PyPI and hard-depends on numcosmo-py (conda-forge only) plus the cosmosis/cobaya sampler connectors we never import — the new uv-overrides.txt drops those three from resolution; the Dockerfile and README carry the --overrides invocation. Plain installs without the blinding extra are unaffected. requires-python moves 3.11 → 3.12: Smokescreen and firecrown both set a 3.12 floor, and the container base (shapepipe:develop) already runs python:3.12-slim-bookworm — this aligns pyproject with the actual runtime. ruff target-version follows. Validated: uv resolution of '.[test,glass,blinding]' with overrides (235 packages, py3.12) and the full firecrown-core + smokescreen import chain in a fresh venv on candide. The CI image build on this branch is the container-side proof. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_019bGhVMAuhTy1gF6DdUc1dQ --- Dockerfile | 13 ++++++++----- README.md | 19 +++++++++++++++++++ pyproject.toml | 27 +++++++++++++++++++++++++-- uv-overrides.txt | 21 +++++++++++++++++++++ 4 files changed, 73 insertions(+), 7 deletions(-) create mode 100644 uv-overrides.txt diff --git a/Dockerfile b/Dockerfile index c791877c..0ebb57a8 100644 --- a/Dockerfile +++ b/Dockerfile @@ -28,8 +28,11 @@ RUN uv pip install --no-cache-dir --upgrade 'cs_util>=0.2.1' WORKDIR /sp_validation COPY . /sp_validation -# Install with the test + glass extras so the image can run the unit suite in CI -# *and* the GLASS map-level mock test. `glass` (Generator for Large Scale -# Structure) ships `glass.ext.camb`; `cosmology` provides the `Cosmology` wrapper -# (`Cosmology.from_camb`) GLASS consumes. Both come in via the `[glass]` extra. -RUN uv pip install --no-cache-dir -e '.[test,glass]' +# Install with the test + glass + blinding extras so the image can run the unit +# suite in CI, the GLASS map-level mock test, *and* the SACC/Smokescreen blinding +# stack. `glass` (Generator for Large Scale Structure) ships `glass.ext.camb`; +# `cosmology` provides the `Cosmology` wrapper (`Cosmology.from_camb`) GLASS +# consumes. The `[blinding]` extra (firecrown + smokescreen) needs the override +# file: firecrown declares conda-forge-only / unused sampler connectors as hard +# deps — see uv-overrides.txt for the full story. +RUN uv pip install --no-cache-dir --overrides uv-overrides.txt -e '.[test,glass,blinding]' diff --git a/README.md b/README.md index e14a2071..640c6b20 100644 --- a/README.md +++ b/README.md @@ -87,6 +87,25 @@ docker run --rm -it ghcr.io/cosmostat/sp_validation:develop python -c "import sp We do not currently build images for Apple Silicon/arm64; however the amd64 images should work on these systems, albeit with reduced performance. +## Local Installation + +Requires Python ≥ 3.12 (the floor is set by the blinding stack; the container +already runs 3.12). With [uv](https://docs.astral.sh/uv/): + +```bash +uv venv --python 3.12 +uv pip install -e '.[test]' +``` + +To also install the data-vector blinding stack (Smokescreen + firecrown, PRD +[#241](https://github.com/CosmoStat/sp_validation/issues/241)), pass the +dependency-override file — firecrown is not pip-resolvable without it (see +`uv-overrides.txt` for why): + +```bash +uv pip install --overrides uv-overrides.txt -e '.[test,blinding]' +``` + ## Flow chart diff --git a/pyproject.toml b/pyproject.toml index bb1fe413..b282d9ad 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,7 +9,10 @@ authors = [ ] license = {text = "MIT"} readme = "README.md" -requires-python = ">=3.11" +# 3.12 floor set by Smokescreen 1.5.6 (and firecrown v1.15); the container base +# (shapepipe:develop) is already python:3.12-slim-bookworm, so this aligns +# pyproject with the actual runtime. +requires-python = ">=3.12" classifiers = [ "License :: OSI Approved :: MIT License", "Programming Language :: Python :: 3", @@ -52,6 +55,10 @@ dependencies = [ "pymaster", "regions", "reproject", + # SACC (LSST DESC's data-vector container) is the standard format for all + # data products from the tomographic round on (PRD #241); sp_validation.sacc_io + # is core library code, so sacc is a core dependency. + "sacc>=2.4,<3", # scipy 1.18 ported FITPACK from Fortran to C, changing the return shape of # RectBivariateSpline(scalar, scalar, grid=False) from 0-d `array(x)` to # shape-(1,) `array([x])`. camb's BBN Y_He predictor (bbn.py) wraps the @@ -106,6 +113,22 @@ 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). +# +# firecrown is not on PyPI and declares conda-forge-only / unused sampler +# connectors as hard deps, so installing this extra requires the dependency +# override file: `uv pip install --overrides uv-overrides.txt -e '.[blinding]'` +# (see uv-overrides.txt; the Dockerfile does this for the container). +blinding = [ + "firecrown @ git+https://github.com/LSSTDESC/firecrown.git@v1.15.1", + "smokescreen==1.5.6", + "pyccl==3.3.4", +] develop = ["sp_validation[test,docs]"] [tool.pytest.ini_options] @@ -124,7 +147,7 @@ markers = [ [tool.ruff] line-length = 88 -target-version = "py311" +target-version = "py312" # Snakemake injects a `snakemake` object into rule scripts at runtime, so ruff # can't see where it's defined. Declaring it a builtin silences the false diff --git a/uv-overrides.txt b/uv-overrides.txt new file mode 100644 index 00000000..813ad055 --- /dev/null +++ b/uv-overrides.txt @@ -0,0 +1,21 @@ +# uv dependency overrides — pass via `--overrides uv-overrides.txt` (or +# UV_OVERRIDE=uv-overrides.txt) to every `uv pip install` against this project. +# +# Why this file exists: firecrown declares its sampler *connectors* as hard +# dependencies, but we use firecrown only as the theory engine for Smokescreen +# blinding (`compute_theory_vector`); sampling stays with CosmoSIS in +# cosmo_inference. Of the three connector deps: +# +# - numcosmo-py exists only on conda-forge, so pip/uv resolution of firecrown +# is *impossible* without an override; +# - cosmosis ships sdist-only (full Fortran/C build with gsl/cfitsio) — a +# heavy, fragile compile in every CI image build, for a connector we never +# import; +# - cobaya is wheel-clean but equally unused. +# +# Each line below replaces the package's requirement (wherever it appears in +# the graph) with one gated on an always-false marker, dropping it from +# resolution. firecrown's likelihood/CCL core imports none of them. +numcosmo-py; python_version < "3" +cosmosis; python_version < "3" +cobaya; python_version < "3" From 1d4826c8abb162a473d5761822361b4023f9f14b Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 02:01:48 +0200 Subject: [PATCH 2/6] deps: firecrown pip-import fixes (numpy<2.5 cap, NumCosmo patch script) Follow-up hardening after empirical validation of the blinding stack on candide. Two facts surfaced that the resolve-only check couldn't see: - firecrown 1.15.1 subclasses npt.NDArray (DataVector); numpy 2.5 made that a non-subclassable typing alias, so firecrown breaks at import on the numpy>=2.0 resolution (2.5.1). The [blinding] extra now carries numpy>=2.2,<2.5 (2.4.3 verified against the full compiled stack + the fast suite), and the Dockerfile requests the bound explicitly to dodge uv #8410 non-movement of an already-installed numpy. - pip-installed firecrown hits NumCosmo (conda-forge-only) at import time through two paths unrelated to cosmic shear: eager re-export of LSST predefined n(z) bins in generators/__init__ (defeating upstream's own lazy __getattr__), and the cluster likelihoods -> lsstdesc-crow -> Ncm.IntegralND C-subclass at module load. scripts/patch_firecrown.py makes the bin re-export lazy, the cluster imports optional, and ships a loud numcosmo_py shim (raises on any real use). Exact-string surgery against pinned v1.15.1, idempotent, fails loudly on a version bump, ends with an import check. Verified: fresh application on a pristine install, idempotent re-run, and the full toy Smokescreen blind end-to-end after patching. CI now also smoke-tests the blinding imports in the built image - the fast suite never imports firecrown, so a broken blinding stack would otherwise ship green. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_019bGhVMAuhTy1gF6DdUc1dQ --- .github/workflows/deploy-image.yml | 6 + Dockerfile | 14 ++ README.md | 1 + pyproject.toml | 9 ++ scripts/patch_firecrown.py | 225 +++++++++++++++++++++++++++++ 5 files changed, 255 insertions(+) create mode 100644 scripts/patch_firecrown.py diff --git a/.github/workflows/deploy-image.yml b/.github/workflows/deploy-image.yml index 192e87a6..13f0fafe 100644 --- a/.github/workflows/deploy-image.yml +++ b/.github/workflows/deploy-image.yml @@ -44,6 +44,12 @@ jobs: - name: Import smoke test run: docker run --rm ${{ steps.meta.outputs.tags }} python -c "import sp_validation" + # The fast suite doesn't import the blinding stack, so a broken + # firecrown/smokescreen install would otherwise ship green. Prove the + # image can actually load it (sacc + patched firecrown + smokescreen). + - name: Blinding-stack import smoke test + run: docker run --rm ${{ steps.meta.outputs.tags }} python -c "import sacc; import firecrown.likelihood; import smokescreen" + # Run the fast test suite against the freshly-built image *before* # pushing, so a failing suite blocks publication. The image carries the # full stack and the test files (COPY . + editable install), so this diff --git a/Dockerfile b/Dockerfile index 0ebb57a8..818db67d 100644 --- a/Dockerfile +++ b/Dockerfile @@ -36,3 +36,17 @@ COPY . /sp_validation # file: firecrown declares conda-forge-only / unused sampler connectors as hard # deps — see uv-overrides.txt for the full story. RUN uv pip install --no-cache-dir --overrides uv-overrides.txt -e '.[test,glass,blinding]' + +# Same uv gotcha as the cs_util upgrade above (astral-sh/uv #8410): if the base +# image already carries a numpy that violates the [blinding] extra's new +# `numpy<2.5` cap (firecrown 1.15.1 breaks on numpy 2.5 at import), the +# editable install won't move it. Request the bound explicitly so the image is +# deterministic either way; numpy 2.4.x is ABI-compatible with the compiled +# stack (verified: pyccl/camb/treecorr/healpy/pymaster + fast suite). +RUN uv pip install --no-cache-dir 'numpy>=2.2,<2.5' + +# firecrown is distributed for conda-forge (where NumCosmo always exists) and +# hits NumCosmo at import time in a pip env, on paths unrelated to our use. +# This patches the installed tree (surgical, pinned-version-checked, loud on +# mismatch) and verifies `import firecrown.likelihood; import smokescreen`. +RUN python scripts/patch_firecrown.py diff --git a/README.md b/README.md index 640c6b20..ef9e601b 100644 --- a/README.md +++ b/README.md @@ -104,6 +104,7 @@ dependency-override file — firecrown is not pip-resolvable without it (see ```bash uv pip install --overrides uv-overrides.txt -e '.[test,blinding]' +python scripts/patch_firecrown.py # make pip-installed firecrown importable without NumCosmo ``` diff --git a/pyproject.toml b/pyproject.toml index b282d9ad..ecdec6f1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -124,10 +124,19 @@ glass = [ # connectors as hard deps, so installing this extra requires the dependency # override file: `uv pip install --overrides uv-overrides.txt -e '.[blinding]'` # (see uv-overrides.txt; the Dockerfile does this for the container). +# After installing this extra, run `python scripts/patch_firecrown.py` — it +# makes pip-installed firecrown importable without NumCosmo (conda-forge-only); +# see that script's docstring for the full story. blinding = [ "firecrown @ git+https://github.com/LSSTDESC/firecrown.git@v1.15.1", "smokescreen==1.5.6", "pyccl==3.3.4", + # firecrown 1.15.1 subclasses npt.NDArray (DataVector); numpy 2.5 turned + # npt.NDArray into a non-subclassable typing alias, breaking firecrown at + # import. firecrown's own env caps numpy<2.4; 2.4.3 is verified against + # the full compiled stack (pyccl/camb/treecorr/healpy/pymaster) + the + # sp_validation fast suite. + "numpy>=2.2,<2.5", ] develop = ["sp_validation[test,docs]"] diff --git a/scripts/patch_firecrown.py b/scripts/patch_firecrown.py new file mode 100644 index 00000000..7b528da1 --- /dev/null +++ b/scripts/patch_firecrown.py @@ -0,0 +1,225 @@ +"""Make pip-installed firecrown importable without NumCosmo. + +Run *inside* the target environment, after installing the ``[blinding]`` extra: + + python scripts/patch_firecrown.py + +Why this exists (PRD #241, PR 1): firecrown is the theory engine for +Smokescreen blinding — only ``compute_theory_vector`` on the SACC-read +cosmic-shear path is used. Upstream distributes firecrown via conda-forge, +where NumCosmo (a GObject-introspection C library, absent from PyPI) is always +present; in a pip/uv environment, firecrown 1.15.1 hits NumCosmo at *import +time* through two paths that have nothing to do with cosmic shear: + +1. ``firecrown/generators/__init__.py`` eagerly re-exports the LSST Y1/Y10 + predefined n(z) bin constants, defeating the lazy ``__getattr__`` that + ``_inferred_galaxy_zdist`` already provides — and computing those constants + imports NumCosmo. +2. ``firecrown/likelihood/__init__.py`` eagerly imports the cluster + likelihoods, which import ``crow`` (lsstdesc-crow), which subclasses a + NumCosmo C class at module load (``class CountsIntegralND(Ncm.IntegralND)``). + +This script (a) restores laziness in ``generators``, (b) makes the cluster +imports optional, and (c) installs a *loud* ``numcosmo_py`` shim so that any +genuine NumCosmo use raises immediately instead of being silently faked. +Everything is exact-string surgery against the pinned firecrown v1.15.1: if a +target string is missing (e.g. after a version bump), the script fails loudly +so the pin and the patch get reviewed together. Idempotent — safe to re-run. + +The right long-term fix is upstream (guarded/lazy imports in firecrown); until +then this file is the entire cost of staying pip-installable. +""" + +import importlib.metadata +import importlib.util +import subprocess +import sys +from pathlib import Path + +EXPECTED_FIRECROWN = "1.15.1" + +GENERATORS_OLD = """\ + # Lazy-loaded bins (via __getattr__) + Y1_LENS_BINS, + Y1_SOURCE_BINS, + Y10_LENS_BINS, + Y10_SOURCE_BINS, + LSST_Y1_LENS_HARMONIC_BIN_COLLECTION, + LSST_Y1_SOURCE_HARMONIC_BIN_COLLECTION, + LSST_Y10_LENS_HARMONIC_BIN_COLLECTION, + LSST_Y10_SOURCE_HARMONIC_BIN_COLLECTION, +) +""" + +GENERATORS_NEW = """\ +) + +# NOTE (sp_validation patch, scripts/patch_firecrown.py): the LSST Y1/Y10 +# predefined bin constants are computed lazily in _inferred_galaxy_zdist via a +# module-level __getattr__ that imports NumCosmo. Importing them EAGERLY here +# forced NumCosmo at `import firecrown.generators` (hence at +# `import firecrown.likelihood`), which pip cannot satisfy. Re-expose them +# lazily instead; the SACC-read cosmic-shear path never touches them. +_LAZY_BIN_NAMES = frozenset( + { + "Y1_LENS_BINS", + "Y1_SOURCE_BINS", + "Y10_LENS_BINS", + "Y10_SOURCE_BINS", + "LSST_Y1_LENS_HARMONIC_BIN_COLLECTION", + "LSST_Y1_SOURCE_HARMONIC_BIN_COLLECTION", + "LSST_Y10_LENS_HARMONIC_BIN_COLLECTION", + "LSST_Y10_SOURCE_HARMONIC_BIN_COLLECTION", + } +) + + +def __getattr__(name): + if name in _LAZY_BIN_NAMES: + from . import _inferred_galaxy_zdist as _z + + return getattr(_z, name) + raise AttributeError(f"module {__name__!r} has no attribute {name!r}") + +""" + +LIKELIHOOD_OLD = """\ +# Cluster statistics +from firecrown.likelihood._binned_cluster import BinnedCluster +from firecrown.likelihood._binned_cluster_number_counts import ( + BinnedClusterNumberCounts, +) +from firecrown.likelihood._binned_cluster_number_counts_shear import ( + BinnedClusterShearProfile, +) +""" + +LIKELIHOOD_NEW = """\ +# Cluster statistics. +# NOTE (sp_validation patch, scripts/patch_firecrown.py): the cluster +# likelihoods import `crow` (lsstdesc-crow), which subclasses NumCosmo C +# classes at module load. NumCosmo is conda-forge-only, so in a pip/uv env +# these imports fail. They are NOT on the cosmic-shear (TwoPoint/WeakLensing) +# path, so they become optional: without NumCosmo the cluster classes are +# unavailable but everything else loads. +try: + from firecrown.likelihood._binned_cluster import BinnedCluster + from firecrown.likelihood._binned_cluster_number_counts import ( + BinnedClusterNumberCounts, + ) + from firecrown.likelihood._binned_cluster_number_counts_shear import ( + BinnedClusterShearProfile, + ) +except (ImportError, RuntimeError, TypeError): # pragma: no cover + BinnedCluster = None # type: ignore[assignment,misc] + BinnedClusterNumberCounts = None # type: ignore[assignment,misc] + BinnedClusterShearProfile = None # type: ignore[assignment,misc] +""" + +SHIM = '''\ +"""Minimal loud shim for numcosmo_py (installed by sp_validation). + +NumCosmo is a GObject-introspection C library available only via conda-forge. +With the companion patches to firecrown (scripts/patch_firecrown.py), the +SACC-read cosmic-shear likelihood path never imports it; this shim provides +the import-time names so the patched package loads, and any genuine numerical +use of NumCosmo raises loudly rather than being silently faked. +""" + + +class _Missing: + def __init__(self, path="numcosmo_py"): + self._p = path + + def __getattr__(self, name): + return _Missing(f"{self._p}.{name}") + + def __call__(self, *a, **k): + raise RuntimeError( + f"{self._p} was called, but NumCosmo is not installed (conda-forge " + "only, not on PyPI). It is not needed for the SACC-read " + "cosmic-shear likelihood path." + ) + + def __getitem__(self, item): + return _Missing(f"{self._p}[...]") + + +Ncm = _Missing("numcosmo_py.Ncm") +Nc = _Missing("numcosmo_py.Nc") +GObject = _Missing("numcosmo_py.GObject") + + +def dict_to_var_dict(*a, **k): + raise RuntimeError("numcosmo_py.dict_to_var_dict unavailable (no NumCosmo)") + + +def var_dict_to_dict(*a, **k): + raise RuntimeError("numcosmo_py.var_dict_to_dict unavailable (no NumCosmo)") +''' + + +def patch_file(path: Path, old: str, new: str) -> str: + text = path.read_text() + if new in text: + return "already patched" + if old not in text: + sys.exit( + f"FATAL: expected text not found in {path}.\n" + "firecrown has probably been bumped past the pinned version this " + "patch targets — review scripts/patch_firecrown.py together with " + "the [blinding] pin in pyproject.toml." + ) + path.write_text(text.replace(old, new, 1)) + return "patched" + + +def main() -> None: + spec = importlib.util.find_spec("firecrown") + if spec is None or spec.origin is None: + sys.exit("FATAL: firecrown is not installed in this environment.") + pkg = Path(spec.origin).parent + + # Metadata, not `import firecrown` — pre-patch, importing is what's broken. + version = importlib.metadata.version("firecrown") + if version != EXPECTED_FIRECROWN: + sys.exit( + f"FATAL: firecrown {version} != expected {EXPECTED_FIRECROWN}; " + "review this patch against the new version before bumping " + "EXPECTED_FIRECROWN." + ) + + print( + "generators/__init__.py:", + patch_file(pkg / "generators" / "__init__.py", GENERATORS_OLD, GENERATORS_NEW), + ) + print( + "likelihood/__init__.py:", + patch_file(pkg / "likelihood" / "__init__.py", LIKELIHOOD_OLD, LIKELIHOOD_NEW), + ) + + # Loud numcosmo_py shim — only when no real NumCosmo is present. + if importlib.util.find_spec("numcosmo_py") is None: + shim_dir = pkg.parent / "numcosmo_py" + shim_dir.mkdir(exist_ok=True) + (shim_dir / "__init__.py").write_text(SHIM) + print("numcosmo_py shim: installed") + else: + print("numcosmo_py shim: skipped (numcosmo_py importable)") + + check = subprocess.run( + [ + sys.executable, + "-c", + "import firecrown.likelihood; import smokescreen", + ], + capture_output=True, + text=True, + ) + if check.returncode != 0: + sys.exit(f"FATAL: post-patch import check failed:\n{check.stderr}") + print("post-patch import check: firecrown.likelihood + smokescreen OK") + + +if __name__ == "__main__": + main() From 6050be086a9f44194ddb1a4cc7a142f7be82ccf9 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 02:36:19 +0200 Subject: [PATCH 3/6] feat(sacc_io): SACC read/write for the standard data-product layout MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add sp_validation.sacc_io: the writer/reader layer for the two-file SACC layout that becomes the package's standard data-product format. {version}.sacc analysis vector — NZ tracers, coarse xi+/-, pseudo-Cl (EE/BB/EB) with a shared BandpowerWindow, COSEBIs, pure E/B, rho/tau PSF diagnostics; one FullCovariance assembled block-diagonally (zero cross-blocks). {version}_xi_fine COSEBIs/pure-EB integration input — same NZ tracers, fine-grid xi+/-, DiagonalCovariance from TreeCorr varxip/varxim. Covariance order is point-insertion order (SACC preserves it bitwise through FITS). Writers insert in the canonical order — xi+ then xi-, Cl (ee, bb, eb), COSEBIs (all En then all Bn), pure E/B in _EB_KEYS order (xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb, matching b_modes.calculate_eb_statistics), rho, then tau — but readers never assume global order: every getter resolves indices through s.indices(dtype, tracers, **tags). assemble_covariance validates that blocks are contiguous, ascending and tile the data vector exactly, failing loud otherwise. Custom data types (pure E/B, rho, tau) all parse under sacc.parse_data_type_name. Tag filters are plain kwargs; the tags={...} form silently selects nothing and is never used. Test suite (test_sacc_io.py, all synthetic and fast): per-writer round-trips (arrays/tags/windows/NZ bitwise), covariance block alignment and zero cross-blocks, assemble_covariance failure modes, DiagonalCovariance round-trip, extract() sub-covariance alignment, a tomographic multi-pair case, reader/writer mirroring on a mixed file, and the end-to-end two-file layout. 20 passed. Co-Authored-By: Claude Opus --- src/sp_validation/sacc_io.py | 537 ++++++++++++++++++++++++ src/sp_validation/tests/test_sacc_io.py | 491 ++++++++++++++++++++++ 2 files changed, 1028 insertions(+) create mode 100644 src/sp_validation/sacc_io.py create mode 100644 src/sp_validation/tests/test_sacc_io.py diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py new file mode 100644 index 00000000..cf03f12c --- /dev/null +++ b/src/sp_validation/sacc_io.py @@ -0,0 +1,537 @@ +"""SACC_IO. + +:Name: sacc_io.py + +:Description: Read/write the standard SACC data-product layout for the + weak-lensing validation package. Two files describe each + catalogue version: + + - ``{version}.sacc`` — the analysis vector: NZ tracers, coarse + ξ±, pseudo-Cℓ (EE/BB/EB) with bandpower windows, COSEBIs, + pure E/B, and ρ/τ PSF diagnostics, all sharing a single + ``FullCovariance`` assembled block-diagonally from the + per-statistic covariances (zero cross-blocks). + - ``{version}_xi_fine.sacc`` — the COSEBIs / pure-EB integration + input: the same NZ tracers, a fine-grid ξ±, and a + ``DiagonalCovariance`` from TreeCorr ``varxip``/``varxim``. + + The covariance order is the point-insertion order (SACC preserves + it bitwise through FITS save/load). Writers below insert in the + canonical order — ξ+ then ξ−, Cℓ (ee, bb, eb), COSEBIs (all Eₙ + then all Bₙ), pure E/B (xip_E, xim_E, xip_B, xim_B, xip_amb, + xim_amb — matching ``b_modes._EB_KEYS``), ρ, then τ — but readers + never assume global order: they resolve indices through + ``Sacc.indices(dtype, tracers, **tags)``. + + Tag filters are plain keyword arguments to ``indices`` / + ``get_data_points`` / ``get_tag``; the ``tags={...}`` form + silently selects nothing and must never be used. +""" + +import numpy as np +import sacc + +PSF_TRACER = "psf_stars" + +# Standard SACC data-type strings. +XI_PLUS = "galaxy_shear_xi_plus" +XI_MINUS = "galaxy_shear_xi_minus" +CL_EE = "galaxy_shear_cl_ee" +CL_BB = "galaxy_shear_cl_bb" +CL_EB = "galaxy_shear_cl_eb" +COSEBI_EE = "galaxy_shear_cosebi_ee" +COSEBI_BB = "galaxy_shear_cosebi_bb" + +# Custom data-type strings (all parse under sacc.parse_data_type_name). +PURE_TYPES = { + "xip_E": "galaxy_shear_xiPureE_plus", + "xim_E": "galaxy_shear_xiPureE_minus", + "xip_B": "galaxy_shear_xiPureB_plus", + "xim_B": "galaxy_shear_xiPureB_minus", + "xip_amb": "galaxy_shear_xiPureAmb_plus", + "xim_amb": "galaxy_shear_xiPureAmb_minus", +} +# Insertion order of the six pure-EB blocks — matches b_modes._EB_KEYS, whose +# order is the [xip_E; xim_E; xip_B; xim_B; xip_amb; xim_amb] layout of the +# treecorr/MC pure-EB covariance (b_modes.calculate_eb_statistics, ~L392). +PURE_KEYS = ("xip_E", "xim_E", "xip_B", "xim_B", "xip_amb", "xim_amb") + +RHO_PLUS = "psf_rho{k}_xi_plus" +RHO_MINUS = "psf_rho{k}_xi_minus" +TAU_PLUS = "galaxyPsf_tau{k}_xi_plus" +TAU_MINUS = "galaxyPsf_tau{k}_xi_minus" + + +def source_name(i): + """SACC tracer name for source redshift bin ``i`` (0-based).""" + return f"source_{i}" + + +def new_sacc(nz, metadata=None): + """Create a Sacc with the survey's NZ (and PSF) tracers. + + Parameters + ---------- + nz : dict or sequence + Redshift distributions, one per source bin. Either a mapping + ``{i: (z, nz)}`` keyed by 0-based bin index, or a sequence of + ``(z, nz)`` array pairs (bin index = position). Tracers are named + ``source_{i}``. + metadata : dict, optional + Key/value pairs stored on ``s.metadata``. + + Returns + ------- + sacc.Sacc + Sacc holding the ``source_{i}`` NZ tracers and the ``psf_stars`` + Misc tracer (needed by ρ/τ diagnostics). + """ + items = nz.items() if isinstance(nz, dict) else enumerate(nz) + s = sacc.Sacc() + for i, (z, nz_i) in items: + s.add_tracer("NZ", source_name(i), np.asarray(z), np.asarray(nz_i)) + s.add_tracer("Misc", PSF_TRACER) + for key, value in (metadata or {}).items(): + s.metadata[key] = value + return s + + +def _pair(bins): + """Resolve a ``(i, j)`` bin pair to the ``(source_i, source_j)`` names.""" + i, j = bins + return (source_name(i), source_name(j)) + + +def add_xi( + s, + bins, + theta, + xip, + xim, + *, + grid, + theta_nom=None, + npairs=None, + weight=None, +): + """Add a real-space shear 2PCF (ξ+ then ξ−) for one tracer pair. + + Parameters + ---------- + s : sacc.Sacc + Target, mutated in place. + bins : tuple of int + Source bin pair ``(i, j)``. + theta : array_like + Angular separations (arcmin) — TreeCorr ``meanr``. + xip, xim : array_like + ξ+ and ξ− at ``theta``. + grid : {'coarse', 'fine'} + Distinguishes the analysis grid from the fine integration grid; + stored as the ``grid`` tag on every point. + theta_nom : array_like, optional + Nominal bin centres — TreeCorr ``rnom`` — stored as ``theta_nom``. + npairs, weight : array_like, optional + TreeCorr pair counts and weights, stored per point. + """ + tracers = _pair(bins) + for dtype, xi in ((XI_PLUS, xip), (XI_MINUS, xim)): + for n, th in enumerate(theta): + tags = {"theta": float(th), "grid": grid} + if theta_nom is not None: + tags["theta_nom"] = float(theta_nom[n]) + if npairs is not None: + tags["npairs"] = float(npairs[n]) + if weight is not None: + tags["weight"] = float(weight[n]) + s.add_data_point(dtype, tracers, float(xi[n]), **tags) + + +def add_pseudo_cl( + s, + bins, + ell_eff, + cl_ee, + cl_bb, + cl_eb, + *, + window_ells, + window_weights, +): + """Add pseudo-Cℓ (EE, BB, EB) with a shared bandpower window. + + Parameters + ---------- + s : sacc.Sacc + Target, mutated in place. + bins : tuple of int + Source bin pair ``(i, j)``. + ell_eff : array_like + Effective multipole of each bandpower. + cl_ee, cl_bb, cl_eb : array_like + EE, BB and EB bandpowers at ``ell_eff``. + window_ells : array_like + Multipoles spanned by the bandpower window matrix (shape ``(nell,)``). + window_weights : array_like + Window matrix ``W`` of shape ``(nell, nbp)`` — one column per + bandpower — from NaMaster ``get_bandpower_windows``. One + ``sacc.BandpowerWindow`` is built and shared across EE/BB/EB. + """ + tracers = _pair(bins) + window = sacc.BandpowerWindow(np.asarray(window_ells), np.asarray(window_weights)) + for dtype, cl in ((CL_EE, cl_ee), (CL_BB, cl_bb), (CL_EB, cl_eb)): + s.add_ell_cl( + dtype, *tracers, np.asarray(ell_eff), np.asarray(cl), window=window + ) + + +def add_cosebis(s, bins, En, Bn, scale_cut): + """Add COSEBIs (all Eₙ then all Bₙ) for one scale cut. + + Parameters + ---------- + s : sacc.Sacc + Target, mutated in place. + bins : tuple of int + Source bin pair ``(i, j)``. + En, Bn : array_like + E- and B-mode COSEBI amplitudes, one per logarithmic mode ``n`` + (1-based). The ``[En; Bn]`` layout matches the COSEBI covariance. + 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. + """ + tracers = _pair(bins) + theta_min, theta_max = scale_cut + for dtype, modes in ((COSEBI_EE, En), (COSEBI_BB, Bn)): + for n, value in enumerate(modes, start=1): + s.add_data_point( + dtype, + tracers, + float(value), + n=n, + theta_min=float(theta_min), + theta_max=float(theta_max), + ) + + +def add_pure_eb(s, bins, theta, xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb): + """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, + xim_B, xip_amb, xim_amb), matching ``b_modes._EB_KEYS`` and the pure-EB + covariance layout. + + Parameters + ---------- + s : sacc.Sacc + Target, mutated in place. + bins : tuple of int + Source bin pair ``(i, j)``. + theta : array_like + 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``. + """ + tracers = _pair(bins) + values = { + "xip_E": xip_E, + "xim_E": xim_E, + "xip_B": xip_B, + "xim_B": xim_B, + "xip_amb": xip_amb, + "xim_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)) + + +def add_rho(s, k, theta, rho_p, rho_m): + """Add a ρ_k PSF statistic (ρ+ then ρ−) on the ``psf_stars`` tracer. + + Parameters + ---------- + s : sacc.Sacc + Target, mutated in place. + k : int + ρ index (0…5). + theta : array_like + Angular separations (arcmin). + rho_p, rho_m : array_like + ρ_k+ and ρ_k− at ``theta``. + """ + tracers = (PSF_TRACER, PSF_TRACER) + for dtype, arr in ((RHO_PLUS.format(k=k), rho_p), (RHO_MINUS.format(k=k), rho_m)): + for n, th in enumerate(theta): + s.add_data_point(dtype, tracers, float(arr[n]), theta=float(th)) + + +def add_tau(s, bins, k, theta, tau_p, tau_m): + """Add a τ_k PSF-leakage statistic (τ+ then τ−). + + Parameters + ---------- + s : sacc.Sacc + Target, mutated in place. + bins : tuple of int + Source bin ``i`` and PSF; the τ tracers are ``(source_i, psf_stars)``. + Only ``bins[0]`` is used. + k : int + τ index (0, 2 or 5). + theta : array_like + Angular separations (arcmin). + tau_p, tau_m : array_like + τ_k+ and τ_k− at ``theta``. + """ + tracers = (source_name(bins[0]), PSF_TRACER) + for dtype, arr in ((TAU_PLUS.format(k=k), tau_p), (TAU_MINUS.format(k=k), tau_m)): + for n, th in enumerate(theta): + s.add_data_point(dtype, tracers, float(arr[n]), theta=float(th)) + + +def assemble_covariance(s, blocks): + """Assemble a block-diagonal ``FullCovariance`` from per-statistic blocks. + + Each block is validated against the current insertion order: its indices + must be contiguous and ascending, the blocks must tile ``0…len(s.mean)`` + exactly (no gap, no overlap), and each block must be square with a size + matching its index span. Any violation raises ``ValueError`` naming the + mismatch. Cross-blocks are left zero. + + Parameters + ---------- + s : sacc.Sacc + Target, mutated in place via ``add_covariance``. + blocks : sequence + Ordered ``(selector, cov)`` pairs (or a mapping of the same). Each + ``selector`` is either an index array, or a ``(data_type, tracers)`` + / ``(data_type, tracers, tags)`` tuple resolved through + ``s.indices``; ``cov`` is the block's dense covariance. + + Returns + ------- + sacc.Sacc + ``s``, with the assembled ``FullCovariance`` attached. + """ + items = blocks.items() if isinstance(blocks, dict) else blocks + ntot = len(s.mean) + full = np.zeros((ntot, ntot)) + cursor = 0 + for selector, cov in items: + idx = _resolve_indices(s, selector) + cov = np.asarray(cov) + if not np.array_equal(idx, np.arange(idx[0], idx[0] + len(idx))): + raise ValueError( + f"covariance block {selector!r} resolves to non-contiguous " + f"or non-ascending indices {idx.tolist()}" + ) + if idx[0] != cursor: + raise ValueError( + f"covariance block {selector!r} starts at index {idx[0]} but " + f"the previous blocks cover through {cursor} — blocks must tile " + "the data vector with no gap or overlap" + ) + if cov.ndim != 2 or cov.shape[0] != cov.shape[1]: + raise ValueError( + f"covariance block {selector!r} must be square; got shape {cov.shape}" + ) + if cov.shape[0] != len(idx): + raise ValueError( + f"covariance block {selector!r} has size {cov.shape[0]} but " + f"spans {len(idx)} data points" + ) + full[np.ix_(idx, idx)] = cov + cursor = idx[-1] + 1 + if cursor != ntot: + raise ValueError( + f"covariance blocks cover {cursor} of {ntot} data points — the " + "blocks must tile the whole data vector" + ) + s.add_covariance(full) + return s + + +def _resolve_indices(s, selector): + """Resolve a covariance-block selector to a sorted index array.""" + if isinstance(selector, (np.ndarray, list, tuple, range)) and not ( + len(selector) in (2, 3) and isinstance(selector[0], str) + ): + return np.asarray(selector, dtype=int) + data_type, tracers = selector[0], selector[1] + tags = selector[2] if len(selector) == 3 else {} + return np.asarray(s.indices(data_type, tuple(tracers), **tags), dtype=int) + + +def add_diagonal_covariance(s, variances): + """Attach a ``DiagonalCovariance`` from a 1-D variance array. + + The 1-D array is passed straight to ``add_covariance`` (never + ``np.diag``), which is what makes SACC store a ``DiagonalCovariance``. + + Parameters + ---------- + s : sacc.Sacc + Target, mutated in place. + variances : array_like + Per-point variances, ``len == len(s.mean)``. + + Returns + ------- + sacc.Sacc + ``s``, with the ``DiagonalCovariance`` attached. + """ + s.add_covariance(np.asarray(variances)) + return s + + +def get_nz(s, i): + """Return ``(z, nz)`` for source bin ``i``.""" + tracer = s.tracers[source_name(i)] + return tracer.z, tracer.nz + + +def get_xi(s, bins, *, grid): + """Return ``(theta, xip, xim)`` for one tracer pair and grid.""" + tracers = _pair(bins) + return ( + _sorted_tag(s, XI_PLUS, tracers, "theta", grid=grid), + _sorted_mean(s, XI_PLUS, tracers, grid=grid), + _sorted_mean(s, XI_MINUS, tracers, grid=grid), + ) + + +def get_pseudo_cl(s, bins): + """Return ``(ell_eff, cl_ee, cl_bb, cl_eb, window)`` for one tracer pair. + + ``window`` is the shared ``sacc.BandpowerWindow`` recovered via + ``get_bandpower_windows``. + """ + tracers = _pair(bins) + ell = _sorted_tag(s, CL_EE, tracers, "ell") + window = s.get_bandpower_windows(s.indices(CL_EE, tracers)) + return ( + ell, + _sorted_mean(s, CL_EE, tracers, _sort_tag="ell"), + _sorted_mean(s, CL_BB, tracers, _sort_tag="ell"), + _sorted_mean(s, CL_EB, tracers, _sort_tag="ell"), + window, + ) + + +def get_cosebis(s, bins, scale_cut=None): + """Return ``(n, En, Bn)`` for one tracer pair. + + Parameters + ---------- + scale_cut : tuple of float, optional + ``(theta_min, theta_max)`` to select when several cuts share the file. + """ + tracers = _pair(bins) + tags = ( + {"theta_min": float(scale_cut[0]), "theta_max": float(scale_cut[1])} + if scale_cut is not None + else {} + ) + modes = _sorted_tag(s, COSEBI_EE, tracers, "n", **tags) + return ( + modes.astype(int), + _sorted_mean(s, COSEBI_EE, tracers, **tags, _sort_tag="n"), + _sorted_mean(s, COSEBI_BB, tracers, **tags, _sort_tag="n"), + ) + + +def get_pure_eb(s, bins): + """Return ``(theta, {key: array})`` for the six pure-EB blocks. + + The dict is keyed by ``PURE_KEYS`` (xip_E, xim_E, …). + """ + tracers = _pair(bins) + theta = _sorted_tag(s, PURE_TYPES["xip_E"], tracers, "theta") + arrays = {key: _sorted_mean(s, PURE_TYPES[key], tracers) for key in PURE_KEYS} + return theta, arrays + + +def get_rho(s, k): + """Return ``(theta, rho_p, rho_m)`` for ρ index ``k``.""" + tracers = (PSF_TRACER, PSF_TRACER) + dt_p, dt_m = RHO_PLUS.format(k=k), RHO_MINUS.format(k=k) + return ( + _sorted_tag(s, dt_p, tracers, "theta"), + _sorted_mean(s, dt_p, tracers), + _sorted_mean(s, dt_m, tracers), + ) + + +def get_tau(s, bins, k): + """Return ``(theta, tau_p, tau_m)`` for τ index ``k`` and source bin.""" + tracers = (source_name(bins[0]), PSF_TRACER) + dt_p, dt_m = TAU_PLUS.format(k=k), TAU_MINUS.format(k=k) + return ( + _sorted_tag(s, dt_p, tracers, "theta"), + _sorted_mean(s, dt_p, tracers), + _sorted_mean(s, dt_m, tracers), + ) + + +def _order(s, data_type, tracers, sort_tag, **tag_filters): + """Indices for a selection, ordered by ascending ``sort_tag``.""" + idx = np.asarray(s.indices(data_type, tracers, **tag_filters), dtype=int) + key = np.array([s.data[i].tags[sort_tag] for i in idx]) + return idx[np.argsort(key)] + + +def _sorted_mean(s, data_type, tracers, _sort_tag="theta", **tag_filters): + """Mean values for a selection, ordered by ``_sort_tag`` ascending.""" + idx = _order(s, data_type, tracers, _sort_tag, **tag_filters) + return s.mean[idx] + + +def _sorted_tag(s, data_type, tracers, tag, **tag_filters): + """Values of ``tag`` for a selection, ordered by that tag ascending.""" + idx = _order(s, data_type, tracers, tag, **tag_filters) + return np.array([s.data[i].tags[tag] for i in idx]) + + +def extract(s, data_type=None, tracers=None, **tag_filters): + """Extract a sub-Sacc (points + aligned covariance sub-block). + + A copy is made and everything *not* matching the selection is removed, so + the covariance sub-block comes out correctly aligned and the original is + untouched. + + Parameters + ---------- + s : sacc.Sacc + Source, left unmodified. + data_type : str, optional + Data type to keep. + tracers : tuple, optional + Tracer pair to keep. + **tag_filters + Tag filters (plain kwargs, e.g. ``grid='fine'``). + + Returns + ------- + sacc.Sacc + New Sacc holding only the selected points. + """ + sub = s.copy() + selection = {} + if tracers is not None: + selection["tracers"] = tuple(tracers) + selection.update(tag_filters) + sub.keep_selection(data_type, **selection) + return sub + + +def save(s, path): + """Write ``s`` to ``path`` (FITS), overwriting any existing file.""" + s.save_fits(path, overwrite=True) + + +def load(path): + """Load a Sacc from ``path`` (FITS).""" + return sacc.Sacc.load_fits(path) diff --git a/src/sp_validation/tests/test_sacc_io.py b/src/sp_validation/tests/test_sacc_io.py new file mode 100644 index 00000000..883a2af0 --- /dev/null +++ b/src/sp_validation/tests/test_sacc_io.py @@ -0,0 +1,491 @@ +"""Tests for :mod:`sp_validation.sacc_io`. + +All synthetic, all fast: build in memory, round-trip through ``tmp_path``, +and assert arrays/tags/windows/covariances come back bitwise-identical and +correctly aligned. No cluster paths. +""" + +import numpy as np +import pytest + +from sp_validation import sacc_io as sio + + +# --------------------------------------------------------------------------- # +# Synthetic builders +# --------------------------------------------------------------------------- # +def _nz(seed, n=50): + rng = np.random.default_rng(seed) + z = np.linspace(0.01, 2.0, n) + return z, rng.uniform(0.1, 1.0, n) + + +def _theta(nbins=6): + return np.geomspace(1.0, 100.0, nbins) + + +def _spd(n, seed): + """Symmetric positive-definite matrix of size ``n``.""" + a = np.random.default_rng(seed).normal(size=(n, n)) + return a @ a.T + n * np.eye(n) + + +def _base_sacc(nbins=1): + """A Sacc with ``nbins`` NZ source tracers plus the PSF tracer.""" + return sio.new_sacc({i: _nz(i) for i in range(nbins)}) + + +# --------------------------------------------------------------------------- # +# 1. Per-writer round-trip (arrays / tags / windows / NZ bitwise) +# --------------------------------------------------------------------------- # +def _roundtrip(s, tmp_path, name="rt"): + path = tmp_path / f"{name}.sacc" + sio.save(s, str(path)) + return sio.load(str(path)) + + +def test_nz_roundtrip(tmp_path): + z, nz = _nz(3) + s = sio.new_sacc({0: (z, nz)}, metadata={"version": "v1.4.6.3"}) + s2 = _roundtrip(s, tmp_path, "nz") + z2, nz2 = sio.get_nz(s2, 0) + assert np.array_equal(z2, z) + assert np.array_equal(nz2, nz) + assert s2.metadata["version"] == "v1.4.6.3" + assert sio.PSF_TRACER in s2.tracers + + +def test_xi_roundtrip(tmp_path): + theta = _theta() + xip, xim = np.arange(6) * 1e-5, np.arange(6) * 2e-5 + npairs, weight = np.arange(6) * 1e3, np.arange(6) * 1.5 + s = _base_sacc() + sio.add_xi( + s, + (0, 0), + theta, + xip, + xim, + grid="coarse", + theta_nom=theta * 1.01, + npairs=npairs, + weight=weight, + ) + s2 = _roundtrip(s, tmp_path, "xi") + th, p, m = sio.get_xi(s2, (0, 0), grid="coarse") + assert np.array_equal(th, theta) + assert np.array_equal(p, xip) + assert np.array_equal(m, xim) + # extra tags survive + idx = s2.indices(sio.XI_PLUS, ("source_0", "source_0"), grid="coarse") + tags = s2.data[idx[0]].tags + assert tags["grid"] == "coarse" + assert set(tags) >= {"theta", "theta_nom", "npairs", "weight", "grid"} + + +def test_pseudo_cl_roundtrip(tmp_path): + ell_eff = np.array([30.0, 120.0, 210.0, 300.0]) + nell, nbp = 50, len(ell_eff) + window_ells = np.arange(2, 2 + nell).astype(float) + W = np.random.default_rng(5).uniform(size=(nell, nbp)) + ee, bb, eb = np.arange(nbp) * 1e-9, np.arange(nbp) * 2e-9, np.arange(nbp) * 3e-9 + s = _base_sacc() + sio.add_pseudo_cl( + s, + (0, 0), + ell_eff, + ee, + bb, + eb, + window_ells=window_ells, + window_weights=W, + ) + s2 = _roundtrip(s, tmp_path, "cl") + ell, cl_ee, cl_bb, cl_eb, window = sio.get_pseudo_cl(s2, (0, 0)) + assert np.array_equal(ell, ell_eff) + assert np.array_equal(cl_ee, ee) + assert np.array_equal(cl_bb, bb) + assert np.array_equal(cl_eb, eb) + assert np.array_equal(window.weight, W) + assert np.array_equal(window.values, window_ells) + + +def test_cosebis_roundtrip(tmp_path): + En, Bn = np.arange(1, 11) * 1e-6, np.arange(1, 11) * 1e-7 + s = _base_sacc() + sio.add_cosebis(s, (0, 0), En, Bn, (1.0, 100.0)) + s2 = _roundtrip(s, tmp_path, "cosebi") + n, E, B = sio.get_cosebis(s2, (0, 0)) + assert np.array_equal(n, np.arange(1, 11)) + assert np.array_equal(E, En) + assert np.array_equal(B, Bn) + idx = s2.indices(sio.COSEBI_EE, ("source_0", "source_0")) + assert s2.data[idx[0]].tags["theta_min"] == 1.0 + assert s2.data[idx[0]].tags["theta_max"] == 100.0 + + +def test_pure_eb_roundtrip(tmp_path): + theta = _theta() + arrays = {key: np.arange(6) * (i + 1) * 1e-6 for i, key in enumerate(sio.PURE_KEYS)} + s = _base_sacc() + sio.add_pure_eb(s, (0, 0), theta, **arrays) + s2 = _roundtrip(s, tmp_path, "pureeb") + th, back = sio.get_pure_eb(s2, (0, 0)) + assert np.array_equal(th, theta) + for key in sio.PURE_KEYS: + assert np.array_equal(back[key], arrays[key]) + + +def test_rho_roundtrip(tmp_path): + theta = _theta() + s = _base_sacc() + for k in range(6): + sio.add_rho( + s, k, theta, np.arange(6) * (k + 1) * 1e-6, np.arange(6) * (k + 1) * 2e-6 + ) + s2 = _roundtrip(s, tmp_path, "rho") + for k in range(6): + th, p, m = sio.get_rho(s2, k) + assert np.array_equal(th, theta) + assert np.array_equal(p, np.arange(6) * (k + 1) * 1e-6) + assert np.array_equal(m, np.arange(6) * (k + 1) * 2e-6) + + +def test_tau_roundtrip(tmp_path): + theta = _theta() + s = _base_sacc() + for k in (0, 2, 5): + sio.add_tau( + s, + (0, 0), + k, + theta, + np.arange(6) * (k + 1) * 1e-6, + np.arange(6) * (k + 1) * 2e-6, + ) + s2 = _roundtrip(s, tmp_path, "tau") + for k in (0, 2, 5): + th, p, m = sio.get_tau(s2, (0, 0), k) + assert np.array_equal(th, theta) + assert np.array_equal(p, np.arange(6) * (k + 1) * 1e-6) + assert np.array_equal(m, np.arange(6) * (k + 1) * 2e-6) + assert s2.data[s2.indices(sio.TAU_PLUS.format(k=k))[0]].tracers == ( + "source_0", + sio.PSF_TRACER, + ) + + +# --------------------------------------------------------------------------- # +# 2. Covariance alignment +# --------------------------------------------------------------------------- # +def _multi_statistic_sacc(): + """Sacc with ξ+/ξ−, Cℓ (ee/bb/eb) and COSEBIs, ready for a covariance.""" + theta = _theta() + s = _base_sacc() + sio.add_xi( + s, (0, 0), theta, np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="coarse" + ) + ell = np.array([30.0, 120.0, 210.0]) + W = np.random.default_rng(0).uniform(size=(20, 3)) + sio.add_pseudo_cl( + s, + (0, 0), + ell, + np.arange(3) * 1e-9, + np.arange(3) * 2e-9, + np.arange(3) * 3e-9, + window_ells=np.arange(2, 22).astype(float), + window_weights=W, + ) + sio.add_cosebis( + s, (0, 0), np.arange(1, 6) * 1e-6, np.arange(1, 6) * 1e-7, (1.0, 100.0) + ) + return s + + +def test_assemble_covariance_alignment(tmp_path): + s = _multi_statistic_sacc() + tr = ("source_0", "source_0") + # Block selectors in canonical (insertion) order. + xi = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)]) + cl = np.concatenate( + [s.indices(sio.CL_EE, tr), s.indices(sio.CL_BB, tr), s.indices(sio.CL_EB, tr)] + ) + co = np.concatenate([s.indices(sio.COSEBI_EE, tr), s.indices(sio.COSEBI_BB, tr)]) + cov_xi, cov_cl, cov_co = _spd(len(xi), 1), _spd(len(cl), 2), _spd(len(co), 3) + sio.assemble_covariance(s, [(xi, cov_xi), (cl, cov_cl), (co, cov_co)]) + s2 = _roundtrip(s, tmp_path, "cov") + assert type(s2.covariance).__name__ == "FullCovariance" + dense = s2.covariance.dense + # each block's sub-covariance is exactly what went in + assert np.array_equal(dense[np.ix_(xi, xi)], cov_xi) + assert np.array_equal(dense[np.ix_(cl, cl)], cov_cl) + assert np.array_equal(dense[np.ix_(co, co)], cov_co) + # zero cross-blocks + assert np.array_equal(dense[np.ix_(xi, cl)], np.zeros((len(xi), len(cl)))) + assert np.array_equal(dense[np.ix_(xi, co)], np.zeros((len(xi), len(co)))) + assert np.array_equal(dense[np.ix_(cl, co)], np.zeros((len(cl), len(co)))) + + +def test_assemble_covariance_selector_tuples(): + """Blocks addressed by (data_type, tracers) tuples, not raw indices.""" + s = _multi_statistic_sacc() + tr = ("source_0", "source_0") + xi = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)]) + cl = np.concatenate( + [s.indices(sio.CL_EE, tr), s.indices(sio.CL_BB, tr), s.indices(sio.CL_EB, tr)] + ) + sio.assemble_covariance( + s, + [ + (xi, _spd(len(xi), 1)), + (cl, _spd(len(cl), 2)), + ((sio.COSEBI_EE, tr), _spd(len(s.indices(sio.COSEBI_EE, tr)), 3)), + ((sio.COSEBI_BB, tr), _spd(len(s.indices(sio.COSEBI_BB, tr)), 4)), + ], + ) + assert type(s.covariance).__name__ == "FullCovariance" + assert s.covariance.dense.shape == (len(s.mean), len(s.mean)) + + +# --------------------------------------------------------------------------- # +# 3. assemble_covariance failure modes +# --------------------------------------------------------------------------- # +def test_assemble_covariance_wrong_dimension(): + s = _base_sacc() + sio.add_xi( + s, (0, 0), _theta(), np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="coarse" + ) + idx = np.arange(len(s.mean)) + with pytest.raises(ValueError, match="span"): + sio.assemble_covariance(s, [(idx, _spd(len(idx) - 1, 1))]) + + +def test_assemble_covariance_non_contiguous(): + s = _base_sacc() + sio.add_xi( + s, (0, 0), _theta(), np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="coarse" + ) + idx = np.array([0, 2, 4, 6, 8, 10, 1, 3]) # not contiguous/ascending + with pytest.raises(ValueError, match="non-contiguous"): + sio.assemble_covariance(s, [(idx, _spd(len(idx), 1))]) + + +def test_assemble_covariance_missing_coverage(): + s = _multi_statistic_sacc() + tr = ("source_0", "source_0") + xi = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)]) + with pytest.raises(ValueError, match="tile"): + sio.assemble_covariance( + s, [(xi, _spd(len(xi), 1))] + ) # leaves cl+cosebi uncovered + + +def test_assemble_covariance_non_square(): + s = _base_sacc() + sio.add_xi( + s, (0, 0), _theta(), np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="coarse" + ) + idx = np.arange(len(s.mean)) + with pytest.raises(ValueError, match="square"): + sio.assemble_covariance(s, [(idx, np.ones((len(idx), len(idx) - 1)))]) + + +def test_assemble_covariance_overlap(): + s = _multi_statistic_sacc() + tr = ("source_0", "source_0") + xi = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)]) + # second block starts before the first ended -> gap/overlap error + with pytest.raises(ValueError, match="tile|gap|overlap"): + sio.assemble_covariance(s, [(xi, _spd(len(xi), 1)), (xi, _spd(len(xi), 2))]) + + +# --------------------------------------------------------------------------- # +# 4. Fine file: DiagonalCovariance from 1-D variances +# --------------------------------------------------------------------------- # +def test_diagonal_covariance_roundtrip(tmp_path): + theta = np.geomspace(0.1, 250.0, 40) + n = len(theta) + xip, xim = np.arange(n) * 1e-5, np.arange(n) * 2e-5 + varxip, varxim = np.arange(1, n + 1) * 1e-12, np.arange(1, n + 1) * 2e-12 + s = _base_sacc() + sio.add_xi(s, (0, 0), theta, xip, xim, grid="fine") + variances = np.concatenate([varxip, varxim]) # [xip; xim] order + sio.add_diagonal_covariance(s, variances) + assert type(s.covariance).__name__ == "DiagonalCovariance" + s2 = _roundtrip(s, tmp_path, "fine") + assert type(s2.covariance).__name__ == "DiagonalCovariance" + assert np.array_equal(np.diag(s2.covariance.dense), variances) + + +# --------------------------------------------------------------------------- # +# 5. extract(): sub-covariance alignment; original untouched; tag filter +# --------------------------------------------------------------------------- # +def test_extract_subblock_and_original_untouched(): + s = _multi_statistic_sacc() + tr = ("source_0", "source_0") + xi = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)]) + cl = np.concatenate( + [s.indices(sio.CL_EE, tr), s.indices(sio.CL_BB, tr), s.indices(sio.CL_EB, tr)] + ) + co = np.concatenate([s.indices(sio.COSEBI_EE, tr), s.indices(sio.COSEBI_BB, tr)]) + cov_co = _spd(len(co), 3) + sio.assemble_covariance( + s, [(xi, _spd(len(xi), 1)), (cl, _spd(len(cl), 2)), (co, cov_co)] + ) + n_before = len(s.mean) + sub = sio.extract(s, data_type=sio.COSEBI_EE, tracers=tr) + # original untouched + assert len(s.mean) == n_before + # subset covariance equals the COSEBI-EE diagonal sub-block + ee_local = np.arange(len(s.indices(sio.COSEBI_EE, tr))) + assert np.array_equal(sub.covariance.dense, cov_co[np.ix_(ee_local, ee_local)]) + + +def test_extract_tag_filter(): + theta = _theta() + s = _base_sacc() + sio.add_xi( + s, (0, 0), theta, np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="coarse" + ) + sio.add_xi(s, (0, 0), theta, np.arange(6) * 3e-5, np.arange(6) * 4e-5, grid="fine") + sub = sio.extract( + s, data_type=sio.XI_PLUS, tracers=("source_0", "source_0"), grid="fine" + ) + assert len(sub.mean) == len(theta) + assert set(sub.get_tag("grid", sio.XI_PLUS)) == {"fine"} + + +# --------------------------------------------------------------------------- # +# 6. Tomographic case: >=2 bins, >=3 pairs, per-pair selection +# --------------------------------------------------------------------------- # +def test_tomographic_per_pair_selection(tmp_path): + theta = _theta() + s = _base_sacc(nbins=2) + pairs = [(0, 0), (0, 1), (1, 1)] + for k, (i, j) in enumerate(pairs): + sio.add_xi( + s, + (i, j), + theta, + np.arange(6) * (k + 1) * 1e-5, + np.arange(6) * (k + 1) * 2e-5, + grid="coarse", + ) + s2 = _roundtrip(s, tmp_path, "tomo") + for k, (i, j) in enumerate(pairs): + th, p, m = sio.get_xi(s2, (i, j), grid="coarse") + assert np.array_equal(th, theta) + assert np.array_equal(p, np.arange(6) * (k + 1) * 1e-5) + assert np.array_equal(m, np.arange(6) * (k + 1) * 2e-5) + # selecting one pair does not bleed into another + assert len(s2.indices(sio.XI_PLUS, ("source_0", "source_1"))) == len(theta) + assert len(s2.indices(sio.XI_PLUS, ("source_1", "source_1"))) == len(theta) + + +# --------------------------------------------------------------------------- # +# 7. Readers mirror writers on a mixed file +# --------------------------------------------------------------------------- # +def test_readers_on_mixed_file(tmp_path): + theta = _theta() + s = _base_sacc() + sio.add_xi( + s, (0, 0), theta, np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="coarse" + ) + ell = np.array([30.0, 120.0, 210.0]) + W = np.random.default_rng(0).uniform(size=(20, 3)) + sio.add_pseudo_cl( + s, + (0, 0), + ell, + np.arange(3) * 1e-9, + np.arange(3) * 2e-9, + np.arange(3) * 3e-9, + window_ells=np.arange(2, 22).astype(float), + window_weights=W, + ) + sio.add_cosebis( + s, (0, 0), np.arange(1, 6) * 1e-6, np.arange(1, 6) * 1e-7, (1.0, 100.0) + ) + pure = {key: np.arange(6) * (i + 1) * 1e-6 for i, key in enumerate(sio.PURE_KEYS)} + sio.add_pure_eb(s, (0, 0), theta, **pure) + for k in range(6): + sio.add_rho( + s, k, theta, np.arange(6) * (k + 1) * 1e-7, np.arange(6) * (k + 1) * 2e-7 + ) + for k in (0, 2, 5): + sio.add_tau( + s, + (0, 0), + k, + theta, + np.arange(6) * (k + 1) * 1e-8, + np.arange(6) * (k + 1) * 2e-8, + ) + s2 = _roundtrip(s, tmp_path, "mixed") + + _, p, m = sio.get_xi(s2, (0, 0), grid="coarse") + assert np.array_equal(p, np.arange(6) * 1e-5) and np.array_equal( + m, np.arange(6) * 2e-5 + ) + ell_r, ee, bb, eb, win = sio.get_pseudo_cl(s2, (0, 0)) + assert np.array_equal(ell_r, ell) and np.array_equal(win.weight, W) + n, E, B = sio.get_cosebis(s2, (0, 0)) + assert np.array_equal(E, np.arange(1, 6) * 1e-6) + _, back = sio.get_pure_eb(s2, (0, 0)) + for key in sio.PURE_KEYS: + assert np.array_equal(back[key], pure[key]) + for k in range(6): + _, rp, rm = sio.get_rho(s2, k) + assert np.array_equal(rp, np.arange(6) * (k + 1) * 1e-7) + for k in (0, 2, 5): + _, tp, tm = sio.get_tau(s2, (0, 0), k) + assert np.array_equal(tp, np.arange(6) * (k + 1) * 1e-8) + + +# --------------------------------------------------------------------------- # +# 8. End-to-end two-file layout for a synthetic catalogue version +# --------------------------------------------------------------------------- # +def test_end_to_end_two_file_layout(tmp_path): + version = "vSYNTH" + theta_c = _theta(20) + theta_f = np.geomspace(0.1, 250.0, 200) + + # analysis file + s = _base_sacc() + sio.add_xi( + s, (0, 0), theta_c, np.arange(20) * 1e-5, np.arange(20) * 2e-5, grid="coarse" + ) + sio.add_cosebis( + s, (0, 0), np.arange(1, 11) * 1e-6, np.arange(1, 11) * 1e-7, (1.0, 100.0) + ) + tr = ("source_0", "source_0") + xi = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)]) + co = np.concatenate([s.indices(sio.COSEBI_EE, tr), s.indices(sio.COSEBI_BB, tr)]) + sio.assemble_covariance(s, [(xi, _spd(len(xi), 1)), (co, _spd(len(co), 2))]) + sio.save(s, str(tmp_path / f"{version}.sacc")) + + # fine file + sf = _base_sacc() + sio.add_xi( + sf, (0, 0), theta_f, np.arange(200) * 1e-5, np.arange(200) * 2e-5, grid="fine" + ) + variances = np.concatenate([np.arange(1, 201) * 1e-12, np.arange(1, 201) * 2e-12]) + sio.add_diagonal_covariance(sf, variances) + sio.save(sf, str(tmp_path / f"{version}_xi_fine.sacc")) + + # reload both, verify everything + a = sio.load(str(tmp_path / f"{version}.sacc")) + f = sio.load(str(tmp_path / f"{version}_xi_fine.sacc")) + + th_c, p_c, _ = sio.get_xi(a, (0, 0), grid="coarse") + assert np.array_equal(th_c, theta_c) and np.array_equal(p_c, np.arange(20) * 1e-5) + n, E, B = sio.get_cosebis(a, (0, 0)) + assert np.array_equal(n, np.arange(1, 11)) + assert type(a.covariance).__name__ == "FullCovariance" + assert a.covariance.dense.shape == (len(a.mean), len(a.mean)) + + th_f, p_f, _ = sio.get_xi(f, (0, 0), grid="fine") + assert np.array_equal(th_f, theta_f) and np.array_equal(p_f, np.arange(200) * 1e-5) + assert type(f.covariance).__name__ == "DiagonalCovariance" + assert np.array_equal(np.diag(f.covariance.dense), variances) From 66bbbc0e6ad27b901fbd3009132d011dd78d9fa2 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 02:50:52 +0200 Subject: [PATCH 4/6] fix(sacc_io): enforce ascending grids; readers in insertion order MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fresh-eyes review caught a correctness bug: readers re-sorted selections by theta/ell/n, but covariance blocks and bandpower windows stay in insertion order. On a non-ascending grid the reader output silently desynchronised from its covariance, and get_pseudo_cl returned sorted cl arrays against unsorted window columns — internally inconsistent within one return tuple. Fix by construction, not by sort: - Writers validate their grids. add_xi/add_pure_eb/add_rho/add_tau require strictly ascending theta; add_pseudo_cl requires strictly ascending ell_eff (add_cosebis is inherently safe — it enumerates the mode index). Out-of-order grids raise a loud ValueError naming the argument. - Readers drop the sort entirely and return in s.indices (insertion) order, so every getter is covariance- and window-aligned for ANY file, and ascending for canonical files. The _sorted_* helpers are replaced by plain insertion-order accessors (_mean/_tag). Also: - _pair normalises (i, j) -> sorted, so get_xi(s, (1, 0)) addresses the same symmetric shear-shear pair as (0, 1) instead of a silent empty read. - Module docstring documents the tomographic ξ covariance ordering: insertion is pair-major ([pair0 xip; pair0 xim; pair1 xip; …]), supplied to assemble_covariance as one contiguous block matching add_xi call order; type-major converters (DES 2pt-FITS) permute explicitly via s.indices. - extract() docstring states tracers takes SACC names, not integer bins. New tests (26 total, was 20): writers reject non-ascending theta/ell; (1, 0) == (0, 1) round-trip; a 3-pair tomographic ξ covariance assembled as one contiguous pair-major block with per-pair sub-blocks recovered via extract(); get_pseudo_cl window column j <-> ell_eff[j] via window_ind tags. Co-Authored-By: Claude Opus --- src/sp_validation/sacc_io.py | 114 +++++++++++++------ src/sp_validation/tests/test_sacc_io.py | 140 ++++++++++++++++++++++++ 2 files changed, 218 insertions(+), 36 deletions(-) diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index cf03f12c..4ed85d0f 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -26,6 +26,18 @@ Tag filters are plain keyword arguments to ``indices`` / ``get_data_points`` / ``get_tag``; the ``tags={...}`` form silently selects nothing and must never be used. + + Tomographic ξ ordering: each ``add_xi`` call inserts one tracer + pair as ``[xip; xim]``, so a multi-pair vector is *pair-major* + — ``[pair_0 xip; pair_0 xim; pair_1 xip; …]`` — not type-major + (``[all xip; all xim]``). A tomographic ξ covariance with + cross-pair correlations is therefore supplied to + ``assemble_covariance`` as ONE contiguous block spanning the + consecutive ``add_xi`` calls, ordered pair-by-pair to match + insertion. Writers must call ``add_xi`` in the same pair order + the covariance was built in. Converters that need a type-major + layout (e.g. the DES 2pt-FITS convention) permute explicitly via + ``s.indices`` rather than assuming global order. """ import numpy as np @@ -97,11 +109,32 @@ def new_sacc(nz, metadata=None): def _pair(bins): - """Resolve a ``(i, j)`` bin pair to the ``(source_i, source_j)`` names.""" - i, j = bins + """Resolve a ``(i, j)`` bin pair to the ``(source_i, source_j)`` names. + + The pair is normalised to ``i <= j``: shear-shear statistics are symmetric + in the tracer pair, and SACC stores each pair under one ordering, so + ``(1, 0)`` must address the same points as ``(0, 1)``. + """ + i, j = sorted(bins) return (source_name(i), source_name(j)) +def _check_ascending(name, values): + """Require ``values`` to be strictly ascending; else raise ValueError. + + Insertion order is the covariance (and bandpower-window) order, and readers + return points in insertion order, so an out-of-order grid would silently + desynchronise a data vector from its covariance. Enforce monotonicity at + write time instead. + """ + values = np.asarray(values) + if not np.all(np.diff(values) > 0): + raise ValueError( + f"{name} must be strictly ascending (insertion order is the " + f"covariance order); got {values.tolist()}" + ) + + def add_xi( s, bins, @@ -134,6 +167,7 @@ def add_xi( npairs, weight : array_like, optional TreeCorr pair counts and weights, stored per point. """ + _check_ascending("theta", theta) tracers = _pair(bins) for dtype, xi in ((XI_PLUS, xip), (XI_MINUS, xim)): for n, th in enumerate(theta): @@ -177,6 +211,7 @@ def add_pseudo_cl( bandpower — from NaMaster ``get_bandpower_windows``. One ``sacc.BandpowerWindow`` is built and shared across EE/BB/EB. """ + _check_ascending("ell_eff", ell_eff) tracers = _pair(bins) window = sacc.BandpowerWindow(np.asarray(window_ells), np.asarray(window_weights)) for dtype, cl in ((CL_EE, cl_ee), (CL_BB, cl_bb), (CL_EB, cl_eb)): @@ -234,6 +269,7 @@ def add_pure_eb(s, bins, theta, xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb): xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb : array_like The six pure E/B / ambiguous mode arrays at ``theta``. """ + _check_ascending("theta", theta) tracers = _pair(bins) values = { "xip_E": xip_E, @@ -263,6 +299,7 @@ def add_rho(s, k, theta, rho_p, rho_m): rho_p, rho_m : array_like ρ_k+ and ρ_k− at ``theta``. """ + _check_ascending("theta", theta) tracers = (PSF_TRACER, PSF_TRACER) for dtype, arr in ((RHO_PLUS.format(k=k), rho_p), (RHO_MINUS.format(k=k), rho_m)): for n, th in enumerate(theta): @@ -286,6 +323,7 @@ def add_tau(s, bins, k, theta, tau_p, tau_m): tau_p, tau_m : array_like τ_k+ and τ_k− at ``theta``. """ + _check_ascending("theta", theta) tracers = (source_name(bins[0]), PSF_TRACER) for dtype, arr in ((TAU_PLUS.format(k=k), tau_p), (TAU_MINUS.format(k=k), tau_m)): for n, th in enumerate(theta): @@ -397,9 +435,9 @@ def get_xi(s, bins, *, grid): """Return ``(theta, xip, xim)`` for one tracer pair and grid.""" tracers = _pair(bins) return ( - _sorted_tag(s, XI_PLUS, tracers, "theta", grid=grid), - _sorted_mean(s, XI_PLUS, tracers, grid=grid), - _sorted_mean(s, XI_MINUS, tracers, grid=grid), + _tag(s, XI_PLUS, tracers, "theta", grid=grid), + _mean(s, XI_PLUS, tracers, grid=grid), + _mean(s, XI_MINUS, tracers, grid=grid), ) @@ -407,16 +445,17 @@ def get_pseudo_cl(s, bins): """Return ``(ell_eff, cl_ee, cl_bb, cl_eb, window)`` for one tracer pair. ``window`` is the shared ``sacc.BandpowerWindow`` recovered via - ``get_bandpower_windows``. + ``get_bandpower_windows``; its columns are in the same insertion order as + the returned ``ell_eff``/``cl`` arrays, so window column ``j`` corresponds + to ``ell_eff[j]``. """ tracers = _pair(bins) - ell = _sorted_tag(s, CL_EE, tracers, "ell") window = s.get_bandpower_windows(s.indices(CL_EE, tracers)) return ( - ell, - _sorted_mean(s, CL_EE, tracers, _sort_tag="ell"), - _sorted_mean(s, CL_BB, tracers, _sort_tag="ell"), - _sorted_mean(s, CL_EB, tracers, _sort_tag="ell"), + _tag(s, CL_EE, tracers, "ell"), + _mean(s, CL_EE, tracers), + _mean(s, CL_BB, tracers), + _mean(s, CL_EB, tracers), window, ) @@ -435,11 +474,11 @@ def get_cosebis(s, bins, scale_cut=None): if scale_cut is not None else {} ) - modes = _sorted_tag(s, COSEBI_EE, tracers, "n", **tags) + modes = _tag(s, COSEBI_EE, tracers, "n", **tags) return ( modes.astype(int), - _sorted_mean(s, COSEBI_EE, tracers, **tags, _sort_tag="n"), - _sorted_mean(s, COSEBI_BB, tracers, **tags, _sort_tag="n"), + _mean(s, COSEBI_EE, tracers, **tags), + _mean(s, COSEBI_BB, tracers, **tags), ) @@ -449,8 +488,8 @@ def get_pure_eb(s, bins): The dict is keyed by ``PURE_KEYS`` (xip_E, xim_E, …). """ tracers = _pair(bins) - theta = _sorted_tag(s, PURE_TYPES["xip_E"], tracers, "theta") - arrays = {key: _sorted_mean(s, PURE_TYPES[key], tracers) for key in PURE_KEYS} + theta = _tag(s, PURE_TYPES["xip_E"], tracers, "theta") + arrays = {key: _mean(s, PURE_TYPES[key], tracers) for key in PURE_KEYS} return theta, arrays @@ -459,9 +498,9 @@ def get_rho(s, k): tracers = (PSF_TRACER, PSF_TRACER) dt_p, dt_m = RHO_PLUS.format(k=k), RHO_MINUS.format(k=k) return ( - _sorted_tag(s, dt_p, tracers, "theta"), - _sorted_mean(s, dt_p, tracers), - _sorted_mean(s, dt_m, tracers), + _tag(s, dt_p, tracers, "theta"), + _mean(s, dt_p, tracers), + _mean(s, dt_m, tracers), ) @@ -470,28 +509,26 @@ def get_tau(s, bins, k): tracers = (source_name(bins[0]), PSF_TRACER) dt_p, dt_m = TAU_PLUS.format(k=k), TAU_MINUS.format(k=k) return ( - _sorted_tag(s, dt_p, tracers, "theta"), - _sorted_mean(s, dt_p, tracers), - _sorted_mean(s, dt_m, tracers), + _tag(s, dt_p, tracers, "theta"), + _mean(s, dt_p, tracers), + _mean(s, dt_m, tracers), ) -def _order(s, data_type, tracers, sort_tag, **tag_filters): - """Indices for a selection, ordered by ascending ``sort_tag``.""" - idx = np.asarray(s.indices(data_type, tracers, **tag_filters), dtype=int) - key = np.array([s.data[i].tags[sort_tag] for i in idx]) - return idx[np.argsort(key)] +def _mean(s, data_type, tracers, **tag_filters): + """Mean values for a selection, in ``s.indices`` (insertion) order. - -def _sorted_mean(s, data_type, tracers, _sort_tag="theta", **tag_filters): - """Mean values for a selection, ordered by ``_sort_tag`` ascending.""" - idx = _order(s, data_type, tracers, _sort_tag, **tag_filters) - return s.mean[idx] + Never re-sort: insertion order is the covariance and bandpower-window + order, so returning in ``s.indices`` order keeps every reader aligned with + the covariance for any file (and ascending for canonically-written files, + which the writers enforce). + """ + return s.mean[s.indices(data_type, tracers, **tag_filters)] -def _sorted_tag(s, data_type, tracers, tag, **tag_filters): - """Values of ``tag`` for a selection, ordered by that tag ascending.""" - idx = _order(s, data_type, tracers, tag, **tag_filters) +def _tag(s, data_type, tracers, tag, **tag_filters): + """Values of ``tag`` for a selection, in insertion order.""" + idx = s.indices(data_type, tracers, **tag_filters) return np.array([s.data[i].tags[tag] for i in idx]) @@ -509,7 +546,12 @@ def extract(s, data_type=None, tracers=None, **tag_filters): data_type : str, optional Data type to keep. tracers : tuple, optional - Tracer pair to keep. + Tracer pair to keep, as SACC tracer **names** (e.g. + ``("source_0", "source_0")`` or ``("source_0", "psf_stars")``) — *not* + integer bin indices. This differs deliberately from the ``add_*`` / + ``get_*`` interface, whose ``bins`` argument takes integer pairs: + ``extract`` is the generic selection escape hatch, mirroring + ``Sacc.keep_selection`` and addressing non-source tracers uniformly. **tag_filters Tag filters (plain kwargs, e.g. ``grid='fine'``). diff --git a/src/sp_validation/tests/test_sacc_io.py b/src/sp_validation/tests/test_sacc_io.py index 883a2af0..28bd6627 100644 --- a/src/sp_validation/tests/test_sacc_io.py +++ b/src/sp_validation/tests/test_sacc_io.py @@ -489,3 +489,143 @@ def test_end_to_end_two_file_layout(tmp_path): assert np.array_equal(th_f, theta_f) and np.array_equal(p_f, np.arange(200) * 1e-5) assert type(f.covariance).__name__ == "DiagonalCovariance" assert np.array_equal(np.diag(f.covariance.dense), variances) + + +# --------------------------------------------------------------------------- # +# 9. Ascending-grid enforcement (writers reject out-of-order grids) +# --------------------------------------------------------------------------- # +def test_add_xi_rejects_non_ascending_theta(): + s = _base_sacc() + theta = _theta()[::-1] # descending + with pytest.raises(ValueError, match="theta must be strictly ascending"): + sio.add_xi( + s, (0, 0), theta, np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="coarse" + ) + + +def test_add_pseudo_cl_rejects_non_ascending_ell(): + s = _base_sacc() + ell = np.array([210.0, 30.0, 120.0]) # not ascending + W = np.random.default_rng(0).uniform(size=(20, 3)) + with pytest.raises(ValueError, match="ell_eff must be strictly ascending"): + sio.add_pseudo_cl( + s, + (0, 0), + ell, + np.arange(3) * 1e-9, + np.arange(3) * 2e-9, + np.arange(3) * 3e-9, + window_ells=np.arange(2, 22).astype(float), + window_weights=W, + ) + + +def test_add_pure_eb_rho_tau_reject_non_ascending_theta(): + s = _base_sacc() + theta = _theta()[::-1] + pure = {key: np.arange(6) * 1e-6 for key in sio.PURE_KEYS} + with pytest.raises(ValueError, match="theta must be strictly ascending"): + sio.add_pure_eb(s, (0, 0), theta, **pure) + with pytest.raises(ValueError, match="theta must be strictly ascending"): + sio.add_rho(s, 0, theta, np.arange(6) * 1e-6, np.arange(6) * 2e-6) + with pytest.raises(ValueError, match="theta must be strictly ascending"): + sio.add_tau(s, (0, 0), 0, theta, np.arange(6) * 1e-6, np.arange(6) * 2e-6) + + +# --------------------------------------------------------------------------- # +# 10. Bin-pair normalisation: (1, 0) addresses the same points as (0, 1) +# --------------------------------------------------------------------------- # +def test_bin_pair_normalisation(): + theta = _theta() + xip, xim = np.arange(6) * 1e-5, np.arange(6) * 2e-5 + s = _base_sacc(nbins=2) + sio.add_xi(s, (0, 1), theta, xip, xim, grid="coarse") + th01, p01, m01 = sio.get_xi(s, (0, 1), grid="coarse") + th10, p10, m10 = sio.get_xi(s, (1, 0), grid="coarse") # reversed order + assert np.array_equal(th01, th10) + assert np.array_equal(p01, p10) and np.array_equal(p10, xip) + assert np.array_equal(m01, m10) and np.array_equal(m10, xim) + # writing under (1, 0) lands in the same tracer pair, not a new one + s2 = _base_sacc(nbins=2) + sio.add_xi(s2, (1, 0), theta, xip, xim, grid="coarse") + assert len(s2.indices(sio.XI_PLUS, ("source_0", "source_1"))) == len(theta) + + +# --------------------------------------------------------------------------- # +# 11. Reader/covariance alignment holds for ANY insertion order (the core +# regression the review caught): a tomographic multi-pair ξ covariance +# assembled as ONE contiguous pair-major block, per-pair sub-blocks +# recovered via extract(). +# --------------------------------------------------------------------------- # +def test_tomographic_xi_covariance_one_contiguous_block(): + theta = _theta() + nth = len(theta) + pairs = [(0, 0), (0, 1), (1, 1)] + s = _base_sacc(nbins=2) + for k, (i, j) in enumerate(pairs): + sio.add_xi( + s, + (i, j), + theta, + np.arange(nth) * (k + 1) * 1e-5, + np.arange(nth) * (k + 1) * 2e-5, + grid="coarse", + ) + # All ξ points as one contiguous block in insertion (pair-major) order. + xi_idx = np.arange(len(s.mean)) + assert np.array_equal(xi_idx, np.arange(3 * 2 * nth)) # 3 pairs x [xip; xim] + cov = _spd(len(xi_idx), 11) # dense, cross-pair correlations + sio.assemble_covariance(s, [(xi_idx, cov)]) + + # Per-pair xip sub-block: resolve indices, extract, compare to input. + for i, j in pairs: + idx_p = s.indices( + sio.XI_PLUS, ("source_" + str(min(i, j)), "source_" + str(max(i, j))) + ) + sub = sio.extract( + s, + data_type=sio.XI_PLUS, + tracers=("source_" + str(min(i, j)), "source_" + str(max(i, j))), + ) + assert np.array_equal(sub.covariance.dense, cov[np.ix_(idx_p, idx_p)]) + # readers stay covariance-aligned: get_xi returns in the same order + th, xip, _ = sio.get_xi(s, (i, j), grid="coarse") + assert np.array_equal(th, theta) + assert np.array_equal(s.mean[idx_p], xip) + + +# --------------------------------------------------------------------------- # +# 12. get_pseudo_cl window/cl column correspondence: window column j maps to +# the returned ell_eff[j] (verified through window_ind tags). +# --------------------------------------------------------------------------- # +def test_pseudo_cl_window_column_correspondence(tmp_path): + ell_eff = np.array([30.0, 120.0, 210.0, 300.0]) + nell, nbp = 40, len(ell_eff) + window_ells = np.arange(2, 2 + nell).astype(float) + # Distinct columns so a permutation would be detectable. + W = np.zeros((nell, nbp)) + for b in range(nbp): + W[b * 5 : (b + 1) * 5, b] = 1.0 + ee, bb, eb = np.arange(nbp) * 1e-9, np.arange(nbp) * 2e-9, np.arange(nbp) * 3e-9 + s = _base_sacc() + sio.add_pseudo_cl( + s, + (0, 0), + ell_eff, + ee, + bb, + eb, + window_ells=window_ells, + window_weights=W, + ) + s2 = _roundtrip(s, tmp_path, "clwin") + ell, cl_ee, _, _, window = sio.get_pseudo_cl(s2, (0, 0)) + # returned ell array is in insertion order + assert np.array_equal(ell, ell_eff) + assert np.array_equal(cl_ee, ee) + # window_ind tag on each EE point indexes the matching window column + idx = s2.indices(sio.CL_EE, ("source_0", "source_0")) + for pos, i in enumerate(idx): + col = s2.data[i].tags["window_ind"] + assert col == pos # insertion order preserved => column j <-> ell[j] + assert np.array_equal(window.weight[:, col], W[:, pos]) From db3ef5a0e91a014d28844df9891bddfb94fd1586 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 00:09:37 +0200 Subject: [PATCH 5/6] One SACC file per catalogue version: fine xi rides as grid='fine' points MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The two-file split (analysis + {version}_xi_fine.sacc) was premised on a 10000-bin fine grid; the production operating point (Paper II B-modes) is 1000 bins, where a dense per-pair fine covariance block is ~32 MB and the CosmoCov integration-binning covariance — which feeds pure-E/B and COSEBIs error propagation — has a natural home as a BlockDiagonal block alongside the analysis blocks. Layout test replaced with the one-file end-to-end case (dense fine block, extract() sub-covariance alignment, zero cross-blocks) plus a varxi-diagonal fallback test. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_019R5eiy11Lihkgn4MKufXSp --- src/sp_validation/sacc_io.py | 24 +++++--- src/sp_validation/tests/test_sacc_io.py | 82 +++++++++++++++++++------ 2 files changed, 77 insertions(+), 29 deletions(-) diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index 4ed85d0f..22fe5aa6 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -3,17 +3,23 @@ :Name: sacc_io.py :Description: Read/write the standard SACC data-product layout for the - weak-lensing validation package. Two files describe each + weak-lensing validation package. One file describes each catalogue version: - - ``{version}.sacc`` — the analysis vector: NZ tracers, coarse - ξ±, pseudo-Cℓ (EE/BB/EB) with bandpower windows, COSEBIs, - pure E/B, and ρ/τ PSF diagnostics, all sharing a single - ``FullCovariance`` assembled block-diagonally from the - per-statistic covariances (zero cross-blocks). - - ``{version}_xi_fine.sacc`` — the COSEBIs / pure-EB integration - input: the same NZ tracers, a fine-grid ξ±, and a - ``DiagonalCovariance`` from TreeCorr ``varxip``/``varxim``. + - ``{version}.sacc`` — NZ tracers, coarse ξ±, pseudo-Cℓ + (EE/BB/EB) with bandpower windows, COSEBIs, pure E/B, ρ/τ PSF + diagnostics, and the fine-grid ξ± integration input for + COSEBIs / pure-EB (``grid='fine'`` tagged points). The + covariance is assembled block-diagonally from the + per-statistic covariances (zero cross-blocks): the analysis + blocks first, then a dense per-pair fine-ξ block (the + CosmoCov integration-binning covariance when it exists — it + feeds derived-statistic error propagation — or the TreeCorr + ``varxip``/``varxim`` diagonal as degraded fallback). At the + production fine binning (1000 θ bins) a dense fine block is + ~32 MB per pair; extreme convergence-check grids (10k bins) + degrade to the diagonal fallback rather than forking the + layout. The covariance order is the point-insertion order (SACC preserves it bitwise through FITS save/load). Writers below insert in the diff --git a/src/sp_validation/tests/test_sacc_io.py b/src/sp_validation/tests/test_sacc_io.py index 28bd6627..cb9c5a37 100644 --- a/src/sp_validation/tests/test_sacc_io.py +++ b/src/sp_validation/tests/test_sacc_io.py @@ -444,14 +444,14 @@ def test_readers_on_mixed_file(tmp_path): # --------------------------------------------------------------------------- # -# 8. End-to-end two-file layout for a synthetic catalogue version +# 8. End-to-end one-file layout for a synthetic catalogue version # --------------------------------------------------------------------------- # -def test_end_to_end_two_file_layout(tmp_path): +def test_end_to_end_one_file_layout(tmp_path): version = "vSYNTH" theta_c = _theta(20) theta_f = np.geomspace(0.1, 250.0, 200) - # analysis file + # one file: analysis products first, fine-grid integration input last s = _base_sacc() sio.add_xi( s, (0, 0), theta_c, np.arange(20) * 1e-5, np.arange(20) * 2e-5, grid="coarse" @@ -459,36 +459,78 @@ def test_end_to_end_two_file_layout(tmp_path): sio.add_cosebis( s, (0, 0), np.arange(1, 11) * 1e-6, np.arange(1, 11) * 1e-7, (1.0, 100.0) ) + sio.add_xi( + s, (0, 0), theta_f, np.arange(200) * 1e-5, np.arange(200) * 2e-5, grid="fine" + ) tr = ("source_0", "source_0") - xi = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)]) + xi_c = np.concatenate( + [ + s.indices(sio.XI_PLUS, tr, grid="coarse"), + s.indices(sio.XI_MINUS, tr, grid="coarse"), + ] + ) co = np.concatenate([s.indices(sio.COSEBI_EE, tr), s.indices(sio.COSEBI_BB, tr)]) - sio.assemble_covariance(s, [(xi, _spd(len(xi), 1)), (co, _spd(len(co), 2))]) - sio.save(s, str(tmp_path / f"{version}.sacc")) - - # fine file - sf = _base_sacc() - sio.add_xi( - sf, (0, 0), theta_f, np.arange(200) * 1e-5, np.arange(200) * 2e-5, grid="fine" + xi_f = np.concatenate( + [ + s.indices(sio.XI_PLUS, tr, grid="fine"), + s.indices(sio.XI_MINUS, tr, grid="fine"), + ] + ) + # dense fine block (CosmoCov integration covariance in production) + fine_block = _spd(len(xi_f), 3) + sio.assemble_covariance( + s, + [(xi_c, _spd(len(xi_c), 1)), (co, _spd(len(co), 2)), (xi_f, fine_block)], ) - variances = np.concatenate([np.arange(1, 201) * 1e-12, np.arange(1, 201) * 2e-12]) - sio.add_diagonal_covariance(sf, variances) - sio.save(sf, str(tmp_path / f"{version}_xi_fine.sacc")) + sio.save(s, str(tmp_path / f"{version}.sacc")) - # reload both, verify everything a = sio.load(str(tmp_path / f"{version}.sacc")) - f = sio.load(str(tmp_path / f"{version}_xi_fine.sacc")) th_c, p_c, _ = sio.get_xi(a, (0, 0), grid="coarse") assert np.array_equal(th_c, theta_c) and np.array_equal(p_c, np.arange(20) * 1e-5) n, E, B = sio.get_cosebis(a, (0, 0)) assert np.array_equal(n, np.arange(1, 11)) - assert type(a.covariance).__name__ == "FullCovariance" assert a.covariance.dense.shape == (len(a.mean), len(a.mean)) - th_f, p_f, _ = sio.get_xi(f, (0, 0), grid="fine") + th_f, p_f, _ = sio.get_xi(a, (0, 0), grid="fine") assert np.array_equal(th_f, theta_f) and np.array_equal(p_f, np.arange(200) * 1e-5) - assert type(f.covariance).__name__ == "DiagonalCovariance" - assert np.array_equal(np.diag(f.covariance.dense), variances) + + # extract() of the fine selection pulls the aligned dense sub-covariance + fine = sio.extract(a, sio.XI_PLUS, tr, grid="fine") + idx_p = a.indices(sio.XI_PLUS, tr, grid="fine") + assert np.allclose(fine.covariance.dense, a.covariance.dense[np.ix_(idx_p, idx_p)]) + # zero cross-blocks between analysis and fine points + assert np.all(a.covariance.dense[np.ix_(xi_c, xi_f)] == 0) + + +def test_one_file_layout_diagonal_fine_fallback(tmp_path): + """No CosmoCov covariance: the fine block is np.diag(varxip/varxim).""" + s = _base_sacc() + theta_f = np.geomspace(0.1, 250.0, 50) + sio.add_xi( + s, (0, 0), _theta(6), np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="coarse" + ) + sio.add_xi( + s, (0, 0), theta_f, np.arange(50) * 1e-5, np.arange(50) * 2e-5, grid="fine" + ) + tr = ("source_0", "source_0") + xi_c = np.concatenate( + [ + s.indices(sio.XI_PLUS, tr, grid="coarse"), + s.indices(sio.XI_MINUS, tr, grid="coarse"), + ] + ) + xi_f = np.concatenate( + [ + s.indices(sio.XI_PLUS, tr, grid="fine"), + s.indices(sio.XI_MINUS, tr, grid="fine"), + ] + ) + variances = np.concatenate([np.arange(1, 51) * 1e-12, np.arange(1, 51) * 2e-12]) + sio.assemble_covariance(s, [(xi_c, _spd(len(xi_c), 1)), (xi_f, np.diag(variances))]) + sio.save(s, str(tmp_path / "vDIAG.sacc")) + a = sio.load(str(tmp_path / "vDIAG.sacc")) + assert np.array_equal(np.diag(a.covariance.dense[np.ix_(xi_f, xi_f)]), variances) # --------------------------------------------------------------------------- # From 873ffd311faef59d131add20b97e1617b419cc87 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 11:55:55 +0200 Subject: [PATCH 6/6] Word fine-block covariance source tool-agnostically (OneCovariance go-forward) Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_019R5eiy11Lihkgn4MKufXSp --- src/sp_validation/sacc_io.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index 22fe5aa6..157f409e 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -13,7 +13,7 @@ covariance is assembled block-diagonally from the per-statistic covariances (zero cross-blocks): the analysis blocks first, then a dense per-pair fine-ξ block (the - CosmoCov integration-binning covariance when it exists — it + analytic integration-binning covariance when it exists — it feeds derived-statistic error propagation — or the TreeCorr ``varxip``/``varxim`` diagonal as degraded fallback). At the production fine binning (1000 θ bins) a dense fine block is