sacc_io: SACC read/write layer for the standard data-product layout

[cailmdaley]

Closes #244. Row 2 of the PRD (#241 §8): the sacc_io module — every data product written to and read from SACC, with round-trip tests. This is the layout contract PRs 3 (converter), 4 (migration), and 6 (blinding) build on. The PR touches exactly two files: sacc_io.py and its test file (most of the 1500 lines are tests).

SACC in five sentences

A SACC file is one object holding tracers (here: one n(z) per tomographic bin, named source_0, source_1, …), a flat ordered list of data points, and optionally a covariance. Each data point is a data type string (galaxy_shear_xi_plus, …), a tuple of tracer names, and one value — plus a free-form dict of point tags: per-point metadata whose keys are strings and whose values are whatever FITS can store (floats, ints, strings). Tags do double duty: they carry the point's coordinates (theta, ell, n) and bookkeeping (npairs, weight), and they are the query keys readers select points by (e.g. grid='integration'). The covariance has no linkage to the points other than order: row/column i of the matrix is the i-th inserted point, full stop — which is why insertion order is load-bearing everywhere below. Everything, order included, survives FITS save/load bitwise.

What's in the file

Product	SACC data type	Tracers	Point tags
ξ± (reporting + integration grids)	galaxy_shear_xi_plus/minus (standard)	(source_i, source_j)	theta (arcmin, TreeCorr meanr), theta_nom, npairs, weight, grid
pseudo-Cℓ EE/BB/EB	galaxy_shear_cl_ee/bb/eb (standard)	(source_i, source_j)	ell + shared BandpowerWindow
COSEBIs	galaxy_shear_cosebi_ee/bb (standard)	(source_i, source_j)	n, theta_min, theta_max
pure E/B	galaxy_shear_xiPure{E,B,Amb}_{plus,minus} (custom)	(source_i, source_j)	theta
ρ₀…ρ₅	psf_rho{k}_xi_{plus,minus} (custom)	(psf_stars, psf_stars)	theta
τ₀,τ₂,τ₅	galaxyPsf_tau{k}_xi_{plus,minus} (custom)	(source_i, psf_stars)	theta
n(z)	NZ tracers source_0, source_1, …

Tomography-native throughout: every writer/reader takes a bin pair (i, j); the current single-bin analysis is the (0, 0) case. ξ± lives on two θ grids — the reporting grid (the coarse analysis binning that goes to inference) and the integration grid (the dense input the COSEBIs / pure-E/B integrals consume). Both share data type and tracer pair, so the grid tag is the sole disambiguator; the grid names follow the convention already used throughout pure_eb.py / cosebis.py / b_modes.py.

One file per catalogue version

{version}.sacc carries everything in the table: the analysis vector and the integration-grid ξ±. The covariance is assembled block-diagonally (zero cross-blocks — the same independence assumption today's cosmosis_fitting.py assembly makes): the per-statistic analysis blocks first, then a dense per-pair integration-ξ block — the analytic covariance at the integration binning when it exists (it is the input to the pure-E/B MC error propagation and the COSEBIs covariance transform), with the TreeCorr varxip/varxim diagonal as the explicit degraded fallback.

Two things this file is not:

• Not where blinding happens. Per PRD §4, each blindable part (reporting ξ±, integration ξ±, pseudo-Cℓ) is blinded at birth, the moment it is computed (PR 6) — nobody waits on the terminal file for a blinded product they care about. This file is a gather of already-blinded parts (PR 4); it never sees a plaintext real data vector and is never itself a blinding step.
• Not necessarily the archival product, in full. The integration-grid block is an intermediate — consumed by the COSEBIs / pure-E/B derivation before this file is even assembled — and it dominates the size: at the production binning (1000 θ bins, 0.5–300′) the covariance blocks alone are ~32 MB per bin pair, by far the largest thing in a tomographic file. Open question for review: keep it in {version}.sacc for provenance, or drop it from the terminal gather and let it persist only as its per-part file? The writers support both; it's a PR-4 assembly choice, not a format question.

Beyond the PRD's row-2 list, pure-E/B writers are included: §4 derives pure-E/B from blinded ξ± and PR 6's acceptance test measures them, so they need a home in the layout.

Design notes

• Optional components are genuinely optional. A file carries whichever components of a statistic were actually computed: add_pseudo_cl requires EE (the bandpower window's reference series) but BB/EB are optional; COSEBIs Bₙ and the pure-E/B B/ambiguous pairs likewise (both-or-neither per pair). Composite readers return None for an absent optional component; explicit selections still fail loud.
• Covariances are stored block-diagonally: assemble_covariance hands sacc a list of blocks → BlockDiagonalCovariance, one FITS table per block, Σ block² on disk instead of dense N² (~tens of MB for a 5-bin analysis vector instead of hundreds); merge preserves the blocks.
• Insertion order is the covariance order (see the primer) — alignment holds by construction: writers enforce strictly ascending θ/ℓ grids (loud ValueError otherwise) and readers return insertion order, so reader output, covariance sub-blocks, and bandpower-window columns can never silently diverge. Bin pairs normalize to i≤j. Insertion is pair-major ([pair₀ ξ+; pair₀ ξ−; pair₁ ξ+; …]); type-major consumers (DES 2pt-FITS) permute explicitly via Sacc.indices. Analysis points insert before integration points so covariance blocks tile the mean contiguously.
• Unmatched selections fail loud. sacc.Sacc.indices returns an empty array (warning only) when a selection matches nothing — a typo'd tag or a float scale cut that isn't bitwise-identical to the stored one would silently produce empty data vectors. Every reader, extract(), and covariance-block selector funnels through a shared guard that raises instead.
• assemble_covariance(s, blocks) places per-statistic covariance blocks with zero cross-blocks and fails loud on any gap, overlap, or size mismatch.
• extract() (copy + keep_selection) pulls a subset with its correctly-aligned covariance sub-block, and update_statistic(s, sub) is the value-only merge-back: each sub point matches exactly one target point by (data_type, tracers, tags) and only its value is replaced — insertion order and covariance untouched. Together they are PR 6's extract → conceal → merge-back primitive (applied per part, at birth) and the selection layer PR 3's converter builds on.
• The pseudo-Cℓ writer takes the bandpower window matrix as explicit input: NaMaster computes it today but nothing saves it; PR 4 wires get_bandpower_windows() through.
• The psf_stars Misc tracer sits beside the NZ tracers because a Sacc has one flat tracer namespace — bookkeeping, not physics; it exists solely so ρ/τ points have a referent.

Blinding hygiene (PRD §4, "Mocks vs data")

save(s, path, type=...) requires type='data'|'mock' (no default) and stamps it into metadata; load(path) fails closed — a type='data' file without the concealed=True Smokescreen stamp raises, so an unblinded real data vector can never be loaded silently. allow_unblinded=True is the loudly-documented escape hatch reserved for the blinding/unblinding tooling itself. Mocks load freely, blinded or not.

Merge

merge(saccs) thinly wraps sacc.concatenate_data_sets: per-statistic files combine in insertion order, shared tracers are stored once, covariances combine block-diagonally with the library's all-or-none rule (anything richer is #247), and metadata is unioned with loud conflicts (a type='data' input cannot merge with a type='mock' one). It also guards θ consistency by tag semantics: all θ-carrying writers (ξ±, ρ, τ, pure-E/B) tag grid="reporting" by default, and same-length θ arrays sharing a grid tag value must be bitwise identical — a near-miss grid (the historical ρ/τ-vs-ξ CosmoSIS failure) raises a loud ValueError; grids under different tags (e.g. integration) differ by design and are unconstrained. The same rule covers the harmonic side: pseudo-Cℓ points tag grid="reporting" too, effective-ℓ arrays sharing a tag value must be bitwise identical, and series on the same ℓ grid must carry identical bandpower windows — the windows, not the effective ℓ, are what theory integrates against.

Testing

68 tests, all synthetic and fast: per-writer round-trips (values, tags, windows, n(z) bitwise), covariance block alignment incl. zero cross-blocks, assemble_covariance failure modes, extract sub-covariance alignment, a 2-bin/3-pair tomographic case, reader/writer mirroring on a mixed file, the end-to-end one-file layout, the varxi-diagonal fallback, the full data/mock × concealed/unconcealed load matrix plus the escape hatch, and merge/update_statistic behaviour (point order, shared tracers, block-diagonal covariance, metadata conflicts, values-only substitution).

Two adversarial review passes have run (8 finder angles, every finding independently verified against the sacc 2.4 source). Pass 1 caught reader-side sorting (the silent covariance-desync class this module is built to prevent); pass 2 caught the silent-empty-selection family now closed by the loud guard, plus the merge/update edge cases above — all fixed, with regression tests.

Not in scope here: cosmo_val mixin migration (PR 4), SACC→2pt-FITS converter and OneCovariance glue (PR 3).

— Claude (Fable) on behalf of Cail.

[cailmdaley]

@sachaguer @LisaGoh this should be ready for review! the description should convey the high-level design, let me know if you have thoughts on that or the implementation details. i'm sure we'll make some changes as we run into issues integrating these functions into the codebase in the subsequent PRs, but it'll be good to start from a well thought out foundation

[LisaGoh]

Thank you for this! Just so I understand correctly, this script is meant to be called as part of the Snakemake workflow where each statistic is calculated in cosmo_val, then transformed into separate sacc files, blinded, then once all the components (statistics, covariances, nz) are transformed and calculated, this script is called again to merge them into a single sacc file to rule them all?

Also, should we impose a catch to check that all the theta values are the same? (unless it's already done and I missed it) Because I remember the last time the theta values were not exactly the same for the rho/tau stats which led to CosmoSIS complaining, and we had to manually flush the values.

Finally, can we have the option to save different covariances (eg. one calculated with OneCov, NaMaster etc.) for validation purposes?

[cailmdaley]

Thanks Lisa — good questions, all three.

Workflow: exactly right, with one nuance: this PR is pure library — no scripts and no integration with existing functionality yet. These functions get wired into anything that writes data vectors or n(z)s in the follow-up PRs; the flow you describe (per-statistic parts in cosmo_val → blind → final gather into one file) is precisely the design of the Snakemake migration PR (#247), which calls merge() for that last step.

θ consistency: you didn't miss it — there was no such check, and sacc itself enforces nothing (θ is an opaque tag to the library). Given the ρ/τ history we've now added a guard to merge(): points carrying the same grid tag are declared to live on the same binning, so their θ arrays must be exactly identical — anything less raises a loud ValueError naming the offending statistics, so mismatches get harmonized upstream instead of surfacing as a CosmoSIS complaint. All θ-carrying writers (ξ±, ρ, τ, pure-E/B) tag grid="reporting" by default, so they all sit in one consistency group; grids under different tags (e.g. the integration binning) differ by design and aren't constrained.

Multiple covariances: yes — a SACC file holds exactly one covariance, so we'll simply produce multiple merged files: a fiducial {version}.sacc plus any number of tagged variants. Since the final file is a cheap gather over per-statistic parts, a variant covariance reuses the whole data vector. How variants are declared and named is part of #247.

— Claude (Fable) on behalf of Cail

[sachaguer]

The overall spirit seems good to me. What I don't have enough understanding yet, is whether we really want our final version sacc file to contain everything. I am tempted to say no and that we would tipycally want to remove:

• CL_EB and CL_BB part of the power spectra
• The finely binned 2PCF
• Rho/Tau-statistics (you never know maybe leakage will vanish in the next iteration)
• Their related covariance.

For some of them the existing code seems quite modular to remove these tracers. The only doubt I have is for the null test part of the angular power spectra that do not have a route where it is not saved if it is not given as input or if a certain boolean is True. Might be worth thinking about this at that stage.

[cailmdaley]

thanks @sachaguer. what's the argument for not including ~everything in the final sacc file? my thinking was that including everything gives confidence that the data products passed to inference are the same ones that passed all null tests & validation. for example:

(you never know maybe leakage will vanish in the next iteration)

yes! we should be absolutely sure we're using the leakage-free iteration for inference

it should be cheap to do this, since it's just a merge of existing products, and we can still allow the creation of inference FITS from sacc files that don't have all of these components.

i suppose the counter-argument is that if the leakage vanishes, that means something changed upstream and snakemake will also know to recompute the inference data vectors. but still i like the idea of having a single source of truth. i'm happy to let go of this if there's a drawback i haven't considered though

[sachaguer]

I think the argument is what seems to me an unnecessary headache and use of storage. Typically, for the Cl, we have the cross-covariance for all polarisations across all tomographic bins, so it's huge, and we probably don't need to save this in the input for the inference.

I agree with the idea of a single source of truth, but if this is a merged file, that is already guaranteed by the upstream files used to generate it. As such, our snakemake workflow should be sufficient to determine, e.g., how a Cl_EB file was produced and whether it matches the merged file (which is essentially your final argument).

Essentially, I am not at ease using several times more storage than needed if we are in any case able to trace the source files (with a header for example). That being said, I have not computed how much storage it would cost to store Cl_eb/bb, or not. I just think, for this PR, it would be nice to have the option to include it in the sacc file as an additional degree of freedom if needed later.

[cailmdaley]

Sacha and I just chatted in person; we concluded it would be overkill to include the integration grids and their covariances in the final sacc. However if we take advantage of the block covariance classes in sacc, storing all the reporting-grid data products and covariances in a single file would only be 10s of MBs which should be fine. Nevertheless the storage should be flexible.

This is largely moot since this PR is just I/O, so after adding the block covariance functionality and making the save functions a little more flexible (e.g. not requiring EB data vectors) this is ready to be merged.