diff --git a/.gitignore b/.gitignore index 4af80f09..7386792b 100644 --- a/.gitignore +++ b/.gitignore @@ -213,3 +213,6 @@ residual.md # Generated by docs/source/_ext/catalogue_columns.py from config/columns/ docs/source/catalogue_columns.md + +# Blind records hold a secret seed and never enter git. +*.blind.json diff --git a/CLAUDE.md b/CLAUDE.md index 9ae25ba6..745f590f 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -42,6 +42,7 @@ is the container (full scientific stack pre-built). For a local dev environment: - `cat.py`: Catalogue handling and manipulation - `cosmo_val.py`: Cosmology validation routines - `cosmology.py`: Cosmological calculations and theory +- `blinding.py`: Blinds, the theory `theory(params, s)` they shift by, and `conceal`, which shifts a SACC by one - `galaxy.py`: Galaxy-specific processing - `grammar.py`: Column grammars: ShapePipe v1 -> v2 adapter and `column_map` renames (`adapt`) - `io.py`: Input/output; the catalogue reader (`read_catalogue`, `Catalogue`), which detects FITS/HDF5 layouts @@ -80,6 +81,12 @@ Main configuration in `scripts/calibration/params.py` with parameters: - emcee for MCMC sampling - pyccl for cosmological calculations +## Blinded catalogues +Each `cosmo_val/cat_config.yaml` entry declares `blind: none` or a blind's name. +- For an entry with `blind: `, signal (ξ±, Cℓ, any cosmological statistic) leaves the function that measures it only concealed: compute it and `sacc_io.save(..., blind=)` it in that function and return the saved part, as `CosmologyValidation` and the workflow do. Never measure it from the file and keep the raw values. +- Never set `blind: none` to get a run through. +- Never print, paste or commit a `.blind.json`. + ## Container Usage Nothing is hand-built. CI publishes `ghcr.io/cosmostat/sp_validation:` on every push, and each person keeps their own copy at diff --git a/cosmo_val/cat_config.yaml b/cosmo_val/cat_config.yaml index 5dc90cdd..41eb3281 100644 --- a/cosmo_val/cat_config.yaml +++ b/cosmo_val/cat_config.yaml @@ -1,4 +1,7 @@ +# One entry per catalogue. `blind:` is `none` (public) or the name of the blind +# its signal is concealed under (sp_validation.blinding). DES: + blind: none subdir: /n17data/mkilbing/astro/data/DES pipeline: SP colour: orange @@ -40,6 +43,7 @@ DES: path: psf_y3a1-v29.fits patch_number: 120 SP_axel_v0.0: + blind: none subdir: /n17data/mkilbing/astro/data/CFIS/v0.0 pipeline: SP colour: cyan @@ -82,6 +86,7 @@ SP_axel_v0.0: path: star_cat.fits patch_number: 120 SP_v0.1.1: + blind: none subdir: /n17data/mkilbing/astro/data/CFIS/v0.0 pipeline: SP colour: cyan @@ -116,6 +121,7 @@ SP_v0.1.1: path: star_cat.fits patch_number: 150 SP_v1.3: + blind: none subdir: /n17data/mkilbing/astro/data/CFIS/v1.0/ShapePipe pipeline: SP colour: green @@ -152,6 +158,7 @@ SP_v1.3: path: unions_shapepipe_star_2022_v1.0.3.fits patch_number: 150 SP_v1.3.6: + blind: none subdir: /n17data/UNIONS/WL/v1.3.x pipeline: SP colour: coral @@ -196,6 +203,7 @@ SP_v1.3.6: path: unions_shapepipe_star_2022_v1.0.3.fits patch_number: 150 SP_v1.4.5: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: dodgerblue @@ -240,6 +248,7 @@ SP_v1.4.5: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.5.A: + blind: none subdir: /n17data/guinot/CFIS_3500_cat/catalogues_SPv1_v1.4.5 pipeline: SP colour: red @@ -283,6 +292,7 @@ SP_v1.4.5.A: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.5_bright: + blind: none subdir: /n17data/murray/unions_cats pipeline: SP colour: mediumblue @@ -324,6 +334,7 @@ SP_v1.4.5_bright: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.5_faint: + blind: none subdir: /n17data/murray/unions_cats pipeline: SP colour: forestgreen @@ -365,6 +376,7 @@ SP_v1.4.5_faint: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.6_glass_mock: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: darkgreen @@ -408,6 +420,7 @@ SP_v1.4.6_glass_mock: path: unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.5_intermediate: + blind: none subdir: /n17data/murray/unions_cats pipeline: SP colour: orange @@ -449,6 +462,7 @@ SP_v1.4.5_intermediate: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.6: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: green @@ -493,6 +507,7 @@ SP_v1.4.6: path: unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.6.3: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: orange @@ -537,6 +552,7 @@ SP_v1.4.6.3: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.6.3_B: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: green @@ -581,6 +597,7 @@ SP_v1.4.6.3_B: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.6.3_C: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: green @@ -625,6 +642,7 @@ SP_v1.4.6.3_C: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.6.3_A: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: green @@ -669,6 +687,7 @@ SP_v1.4.6.3_A: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.6_ecut07: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: darkgreen @@ -713,6 +732,7 @@ SP_v1.4.6_ecut07: path: unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.6.1: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: red @@ -757,6 +777,7 @@ SP_v1.4.6.1: path: unions_shapepipe_star_2024_v1.4.a.fits patch_number: 100 SP_v1.4.7: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: darkorchid @@ -803,6 +824,7 @@ SP_v1.4.7: path: unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.8: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: darkorchid @@ -847,6 +869,7 @@ SP_v1.4.8: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.11.2: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: orange @@ -891,6 +914,7 @@ SP_v1.4.11.2: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.11.3: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: red @@ -935,6 +959,7 @@ SP_v1.4.11.3: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.12.3: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: lightblue @@ -979,6 +1004,7 @@ SP_v1.4.12.3: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.13.3: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: cyan @@ -1023,6 +1049,7 @@ SP_v1.4.13.3: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.11.3_ecut07: + blind: none subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP colour: darkblue @@ -1067,6 +1094,7 @@ SP_v1.4.11.3_ecut07: path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits patch_number: 150 SP_v1.4.6.3_uncal: + blind: none pipeline: SP subdir: /n17data/UNIONS/WL/v1.4.x shear: @@ -1082,6 +1110,7 @@ SP_v1.4.6.3_uncal: hdu: 1 patch_number: 100 SP_v1.4.6.3_uncal_w_iv: + blind: none pipeline: SP subdir: /n17data/UNIONS/WL/v1.4.x shear: @@ -1097,6 +1126,7 @@ SP_v1.4.6.3_uncal_w_iv: hdu: 1 patch_number: 100 SP_v1.4.6.3_uncal_w_1: + blind: none pipeline: SP subdir: /n17data/UNIONS/WL/v1.4.x shear: @@ -1112,6 +1142,7 @@ SP_v1.4.6.3_uncal_w_1: hdu: 1 patch_number: 100 SP_v1.4.5_uncal: + blind: none pipeline: SP subdir: /n17data/UNIONS/WL/v1.4.x shear: @@ -1127,6 +1158,7 @@ SP_v1.4.5_uncal: hdu: 1 patch_number: 100 SP_v1.4.7_uncal: + blind: none pipeline: SP subdir: /n17data/UNIONS/WL/v1.4.x shear: @@ -1142,6 +1174,7 @@ SP_v1.4.7_uncal: hdu: 1 patch_number: 100 SP_v1.4.8_uncal: + blind: none pipeline: SP subdir: /n17data/UNIONS/WL/v1.4.x shear: @@ -1163,7 +1196,10 @@ nz: path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_{pipeline}_A.txt paths: output: ./output + # Where the .blind.json records live. + blinds: /n17data/UNIONS/WL/blinds SP_v1.6.6: + blind: none subdir: /n17data/UNIONS/WL/v1.6.x pipeline: SP colour: violet @@ -1210,6 +1246,7 @@ SP_v1.6.6: patch_number: 150 GLASS_mock_validation: + blind: none subdir: /n09data/guerrini/glass_mock_test/results/ pipeline: SP colour: violet diff --git a/papers/bmodes/rules/figures.smk b/papers/bmodes/rules/figures.smk index b9095aaf..9506b81e 100644 --- a/papers/bmodes/rules/figures.smk +++ b/papers/bmodes/rules/figures.smk @@ -82,10 +82,10 @@ def _reporting_cov_path(version): def _xi_integration_path(version): - """Path to fine-binned 2PCF integration file. Unpatched: values only, no covariance.""" + """Path to the fine-binned ξ± integration part (unpatched).""" return ( - f"{COSMO_VAL_OUTPUT}/xi_{version}_tomo_bin_all_minsep={FIDUCIAL['min_sep_int']}" - f"_maxsep={FIDUCIAL['max_sep_int']}_nbins={FIDUCIAL['nbins_int']}_npatch=1.txt" + f"{COSMO_VAL_OUTPUT}/{version}_xi_minsep={FIDUCIAL['min_sep_int']}" + f"_maxsep={FIDUCIAL['max_sep_int']}_nbins={FIDUCIAL['nbins_int']}_npatch=1.sacc" ) diff --git a/papers/bmodes/scripts/bb_covariance_nz_independence.py b/papers/bmodes/scripts/bb_covariance_nz_independence.py index a11f9916..1e4cdb52 100644 --- a/papers/bmodes/scripts/bb_covariance_nz_independence.py +++ b/papers/bmodes/scripts/bb_covariance_nz_independence.py @@ -20,12 +20,12 @@ import matplotlib.pyplot as plt import numpy as np -import treecorr import yaml from astropy.io import fits from plotting_utils import PAPER_MPLSTYLE from pseudo_cl_io import load_pseudo_cl_data +from sp_validation import sacc_io from sp_validation.b_modes import calculate_cosebis plt.style.use(PAPER_MPLSTYLE) @@ -82,10 +82,7 @@ def load_cosebis_diagonals( Returns E_n and B_n covariance diagonals. """ # Load fine-binned 2PCF (need the binning info for COSEBIS calculation) - gg = treecorr.GGCorrelation( - min_sep=min_sep_int, max_sep=max_sep_int, nbins=nbins_int, sep_units="arcmin" - ) - gg.read(xi_integration_path) + gg = sacc_io.xi_correlation(sacc_io.load(xi_integration_path)) # Compute COSEBIS with this realisation's covariance results = calculate_cosebis( @@ -511,8 +508,8 @@ def _from_cli(argv=None): } xi_integration_path = os.path.join( a.cosmo_val_dir, - f"xi_{version}_tomo_bin_all_minsep={min_sep_int}_maxsep={max_sep_int}" - f"_nbins={nbins_int}_npatch={npatch}.txt", + f"{version}_xi_minsep={min_sep_int}_maxsep={max_sep_int}" + f"_nbins={nbins_int}_npatch={npatch}.sacc", ) pseudo_cl_path = os.path.join( a.cosmo_val_dir, f"pseudo_cl_{version}_powspace_nbins=32.sacc" diff --git a/papers/bmodes/scripts/compute_cosebis_pte_single.py b/papers/bmodes/scripts/compute_cosebis_pte_single.py index 6fc5bf92..d80f9eee 100644 --- a/papers/bmodes/scripts/compute_cosebis_pte_single.py +++ b/papers/bmodes/scripts/compute_cosebis_pte_single.py @@ -15,9 +15,9 @@ from pathlib import Path import numpy as np -import treecorr from plotting_utils import compute_chi2_pte +from sp_validation import sacc_io from sp_validation.b_modes import calculate_cosebis @@ -89,14 +89,7 @@ def main(config, xi_integration, cov_integration, out_dir, version=None): version = version if version is not None else fid["version"] nmodes = int(fid["nmodes"]) # 20 for full computation - min_sep_int = fid["min_sep_int"] - max_sep_int = fid["max_sep_int"] - nbins_int = fid["nbins_int"] - - gg = treecorr.GGCorrelation( - min_sep=min_sep_int, max_sep=max_sep_int, nbins=nbins_int, sep_units="arcmin" - ) - gg.read(xi_integration) + gg = sacc_io.xi_correlation(sacc_io.load(xi_integration)) # Reporting theta grid (nbins+1 = 21 edges): geomspace(1', 250', 21) theta_grid = np.geomspace(fid["min_sep"], fid["max_sep"], fid["nbins"] + 1) @@ -183,7 +176,7 @@ def _from_cli(argv=None): ap.add_argument( "--xi-integration", required=True, - help="Fiducial 1000-bin integration-grid TreeCorr xi_pm .txt dump", + help="Fiducial 1000-bin integration-grid ξ± SACC part", ) ap.add_argument( "--cov-integration", diff --git a/papers/bmodes/scripts/cosebis_data_vector.py b/papers/bmodes/scripts/cosebis_data_vector.py index 89631f83..ea431569 100644 --- a/papers/bmodes/scripts/cosebis_data_vector.py +++ b/papers/bmodes/scripts/cosebis_data_vector.py @@ -17,12 +17,12 @@ import matplotlib.pyplot as plt import numpy as np import seaborn as sns -import treecorr from plotting_utils import ( FIG_WIDTH_SINGLE, PAPER_MPLSTYLE, ) +from sp_validation import sacc_io from sp_validation.b_modes import calculate_cosebis plt.style.use(PAPER_MPLSTYLE) @@ -153,21 +153,11 @@ def main(config, xi_integration, cov_integration, out_dir): "full": full_scale_cut, } - min_sep_int = float(config["fiducial"]["min_sep_int"]) - max_sep_int = float(config["fiducial"]["max_sep_int"]) - nbins_int = int(config["fiducial"]["nbins_int"]) - out_dir = Path(out_dir) out_dir.mkdir(parents=True, exist_ok=True) # Load the fiducial fine-binned 2PCF once (integration grid) - gg = treecorr.GGCorrelation( - min_sep=min_sep_int, - max_sep=max_sep_int, - nbins=nbins_int, - sep_units="arcmin", - ) - gg.read(xi_integration) + gg = sacc_io.xi_correlation(sacc_io.load(xi_integration)) # Compute the E_n / B_n data vectors + T^T C_xi T mode covariance at both # scale cuts. Same calculate_cosebis call as the original per-version loop. @@ -235,7 +225,7 @@ def _from_cli(argv=None): ap.add_argument( "--xi-integration", required=True, - help="Fiducial fine-binned (1000-bin) TreeCorr xi_pm .txt dump", + help="Fiducial fine-binned (1000-bin) ξ± SACC part", ) ap.add_argument( "--cov-integration", diff --git a/papers/bmodes/scripts/cosebis_version_comparison.py b/papers/bmodes/scripts/cosebis_version_comparison.py index 51fa9d53..eca3deaa 100644 --- a/papers/bmodes/scripts/cosebis_version_comparison.py +++ b/papers/bmodes/scripts/cosebis_version_comparison.py @@ -13,7 +13,6 @@ import matplotlib.pyplot as plt import numpy as np import seaborn as sns -import treecorr from plotting_utils import ( ERRORBAR_DEFAULTS, FIG_WIDTH_FULL, @@ -26,6 +25,7 @@ version_label, ) +from sp_validation import sacc_io from sp_validation.b_modes import calculate_cosebis plt.style.use(PAPER_MPLSTYLE) @@ -148,7 +148,7 @@ def _create_stacked_bmode_figure( def _xi_integration(results_dir, ver): - return f"{results_dir}/xi_{ver}_tomo_bin_all_minsep=0.5_maxsep=300.0_nbins=1000_npatch=1.txt" + return f"{results_dir}/{ver}_xi_minsep=0.5_maxsep=300.0_nbins=1000_npatch=1.sacc" def _cov_integration(cov_dir, ver): @@ -212,10 +212,6 @@ def main( "full": full_scale_cut, } - min_sep_int = float(config["fiducial"]["min_sep_int"]) - max_sep_int = float(config["fiducial"]["max_sep_int"]) - nbins_int = int(config["fiducial"]["nbins_int"]) - # Color/marker assignment: pair by parent version if ecut versions present has_ecut = any("_ecut" in v for v in versions) if has_ecut: @@ -251,13 +247,7 @@ def main( marker = MARKER_STYLES[i] if i < len(MARKER_STYLES) else "o" fillstyle = "full" - gg = treecorr.GGCorrelation( - min_sep=min_sep_int, - max_sep=max_sep_int, - nbins=nbins_int, - sep_units="arcmin", - ) - gg.read(xi_path) + gg = sacc_io.xi_correlation(sacc_io.load(xi_path)) results = calculate_cosebis( gg, @@ -368,7 +358,7 @@ def _from_cli(argv=None): ap.add_argument( "--results-dir", required=True, - help="COSMO_VAL output dir with per-version integration-grid xi_pm .txt dumps", + help="COSMO_VAL output dir with per-version integration-grid ξ± SACC parts", ) ap.add_argument( "--cov-dir", diff --git a/papers/bmodes/scripts/harmonic_config_cosebis_comparison.py b/papers/bmodes/scripts/harmonic_config_cosebis_comparison.py index 2fc52dae..39cadec3 100644 --- a/papers/bmodes/scripts/harmonic_config_cosebis_comparison.py +++ b/papers/bmodes/scripts/harmonic_config_cosebis_comparison.py @@ -20,7 +20,6 @@ import matplotlib.pyplot as plt import numpy as np import seaborn as sns -import treecorr import yaml from astropy.io import fits from cosmo_numba.B_modes.cosebis import COSEBIS @@ -34,6 +33,7 @@ ) from pseudo_cl_io import load_pseudo_cl_data +from sp_validation import sacc_io from sp_validation.b_modes import calculate_cosebis plt.style.use(PAPER_MPLSTYLE) @@ -139,17 +139,7 @@ def _compute_config_cosebis(xi_path, cov_path, nmodes, scale_cut, config): Returns (En, Bn, cov). """ - min_sep_int = float(config["fiducial"]["min_sep_int"]) - max_sep_int = float(config["fiducial"]["max_sep_int"]) - nbins_int = int(config["fiducial"]["nbins_int"]) - - gg = treecorr.GGCorrelation( - min_sep=min_sep_int, - max_sep=max_sep_int, - nbins=nbins_int, - sep_units="arcmin", - ) - gg.read(xi_path) + gg = sacc_io.xi_correlation(sacc_io.load(xi_path)) results = calculate_cosebis( gg, nmodes=nmodes, scale_cuts=[scale_cut], cov_path=cov_path @@ -854,7 +844,7 @@ def _from_cli(argv=None): required=True, help=( "COSMO_VAL output dir (pseudo_cl SACC parts, pseudo_cl_cov FITS " - "+ xi_integration txt)" + "+ xi_integration SACC part)" ), ) ap.add_argument( @@ -888,9 +878,9 @@ def _from_cli(argv=None): "--fiducial-xi-path", default=None, help=( - "Explicit path to the fiducial 1000-bin integration xi_pm text file " - "reproduced by lc (e.g. xi_SP_v1.4.6.3_leak_corr_tomo_bin_all_minsep=0.5_maxsep=300.0_" - "nbins=1000_npatch=1.txt), overriding the --cosmo-val-dir pattern lookup " + "Explicit path to the fiducial 1000-bin integration ξ± SACC part " + "reproduced by lc (e.g. SP_v1.4.6.3_leak_corr_xi_minsep=0.5_maxsep=300.0_" + "nbins=1000_npatch=1.sacc), overriding the --cosmo-val-dir pattern lookup " "for --fiducial-version." ), ) @@ -946,8 +936,8 @@ def _from_cli(argv=None): if is_fiducial and a.fiducial_xi_path else os.path.join( a.cosmo_val_dir, - f"xi_{ver}_tomo_bin_all_minsep={min_sep_int}_maxsep={max_sep_int}" - f"_nbins={nbins_int}_npatch={npatch}.txt", + f"{ver}_xi_minsep={min_sep_int}_maxsep={max_sep_int}" + f"_nbins={nbins_int}_npatch={npatch}.sacc", ) ) inputs[f"cov_{ver}"] = ( diff --git a/papers/bmodes/scripts/pseudo_cl_io.py b/papers/bmodes/scripts/pseudo_cl_io.py index b7b1c796..b9497105 100644 --- a/papers/bmodes/scripts/pseudo_cl_io.py +++ b/papers/bmodes/scripts/pseudo_cl_io.py @@ -7,8 +7,6 @@ def load_pseudo_cl_data(path): """Return a pseudo-Cℓ part's spectra as ``{"ELL", "EE", "EB", "BB"}`` arrays. The bandpower covariance is a separate FITS product (``pseudo_cl_cov``). - Loading goes through ``sacc_io.load``, so an unblinded real-data part is - refused. """ ell, ee, bb, eb, _window = sacc_io.get_pseudo_cl(sacc_io.load(str(path)), (0, 0)) missing = [name for name, values in (("BB", bb), ("EB", eb)) if values is None] diff --git a/papers/bmodes/scripts/pure_eb_modes.py b/papers/bmodes/scripts/pure_eb_modes.py index 643fe12c..a9c0898d 100644 --- a/papers/bmodes/scripts/pure_eb_modes.py +++ b/papers/bmodes/scripts/pure_eb_modes.py @@ -1,6 +1,6 @@ """Pure E/B modes and their exact covariance for one version. -Reads the fine-grid ξ± (TreeCorr text dump, with its pair weights) and the +Reads the fine-grid ξ± SACC part (with its TreeCorr pair weights) and the Gaussian ξ± covariance on the same grid, and applies ``sp_validation.b_modes.calculate_pure_eb_correlation``: the fixed-operator pure-E/B estimator averaged into the reporting bins, with covariance @@ -8,7 +8,7 @@ pure-mode plot / PTE. python pure_eb_modes.py \ - --xi-integration \ + --xi-integration \ --cov-integration \ --min-sep 1.0 --max-sep 250.0 --nbins 20 \ --min-sep-int 0.5 --max-sep-int 300.0 --nbins-int 1000 \ @@ -20,18 +20,15 @@ import numpy as np +from sp_validation import sacc_io from sp_validation.b_modes import calculate_pure_eb_correlation from sp_validation.sacc_io import PURE_KEYS -def _load_xi(path, nbins): - """``(meanr, xip, xim, weight)`` from a TreeCorr text dump. - - TreeCorr's ASCII header is - r_nom meanr meanlogr xip xim xip_im xim_im sigma_xip sigma_xim weight npairs. - """ - data = np.loadtxt(path, comments="#", max_rows=nbins) - return data[:, 1], data[:, 3], data[:, 4], data[:, 9] +def _load_xi(path): + """``(meanr, xip, xim, weight)`` from a ξ± SACC part.""" + gg = sacc_io.xi_correlation(sacc_io.load(path)) + return gg.meanr, gg.xip, gg.xim, gg.weight def pure_eb_modes( @@ -46,7 +43,7 @@ def pure_eb_modes( out_path, ): results = calculate_pure_eb_correlation( - *_load_xi(xi_integration, nbins_int), + *_load_xi(xi_integration), np.geomspace(min_sep_int, max_sep_int, nbins_int + 1), np.loadtxt(cov_integration), np.geomspace(min_sep, max_sep, nbins + 1), diff --git a/papers/catalog/2025_10_02_xi_sys_rho_tau.py b/papers/catalog/2025_10_02_xi_sys_rho_tau.py index ab95d69c..3a73c3dd 100644 --- a/papers/catalog/2025_10_02_xi_sys_rho_tau.py +++ b/papers/catalog/2025_10_02_xi_sys_rho_tau.py @@ -13,6 +13,8 @@ from getdist import MCSamples, plots from shear_psf_leakage.rho_tau_stat import PSFErrorFit, RhoStat, TauStat +from sp_validation import sacc_io + plt.style.use("./matplotlib_config/paper.mplstyle") plt.rcParams["text.usetex"] = True @@ -127,15 +129,12 @@ xi_psf_sys_quant = np.quantile(xi_psf_sys_samples, quantiles, axis=0) theta = psf_fitter.rho_stat_handler.rho_stats["theta"] - if os.path.exists( + xi_part = ( base_dir - + f"/../{ver}_xi_minsep={theta_min}_maxsep={theta_max}_nbins={nbins}_npatch=1.txt" - ): - xip = np.loadtxt( - base_dir - + f"/../{ver}_xi_minsep={theta_min}_maxsep={theta_max}_nbins={nbins}_npatch=1.txt", - usecols=3, - )[:nbins] ## This will need modifications as the pipeline evolved + + f"/../{ver}_xi_minsep={theta_min}_maxsep={theta_max}_nbins={nbins}_npatch=1.sacc" + ) + if os.path.exists(xi_part): + xip = sacc_io.xi_correlation(sacc_io.load(xi_part)).xip ratio_mean = xi_psf_sys_mean / xip ratio_quant = xi_psf_sys_quant / xip @@ -275,15 +274,12 @@ xi_psf_sys_quant = np.quantile(xi_psf_sys_samples, quantiles, axis=0) theta = psf_fitter.rho_stat_handler.rho_stats["theta"] - if os.path.exists( + xi_part = ( base_dir - + f"/../{ver}_xi_minsep={theta_min}_maxsep={theta_max}_nbins={nbins}_npatch=1.txt" - ): - xip = np.loadtxt( - base_dir - + f"/../{ver}_xi_minsep={theta_min}_maxsep={theta_max}_nbins={nbins}_npatch=1.txt", - usecols=3, - )[:nbins] ## This will need modifications as the pipeline evolved + + f"/../{ver}_xi_minsep={theta_min}_maxsep={theta_max}_nbins={nbins}_npatch=1.sacc" + ) + if os.path.exists(xi_part): + xip = sacc_io.xi_correlation(sacc_io.load(xi_part)).xip ratio_mean = xi_psf_sys_mean / xip ratio_quant = xi_psf_sys_quant / xip diff --git a/pyproject.toml b/pyproject.toml index 13cbadac..7d35d8b5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,11 +27,6 @@ dependencies = [ "camb>=2.0", "clmm", "colorama", - # Blinding closure (core, no extra): cryptography + sacc here, plus the - # Smokescreen fork pin below; pyccl (the theory backend) is already core. - # cryptography and sacc are declared explicitly so the core runtime - # closure is self-documenting and independent of fork-metadata drift. - "cryptography", # Track cs_util's develop branch directly (git dependency) rather than a # PyPI pin: the two repos are iterating together heavily and cs_util # releases are infrequent. This PR's cosmology repoint needs get_cosmo / diff --git a/scripts/map2.py b/scripts/map2.py deleted file mode 100755 index 361fe48e..00000000 --- a/scripts/map2.py +++ /dev/null @@ -1,151 +0,0 @@ -#!/usr/bin/env python3 - -"""map2.py - -Compute aperture-mass dispersion using ``treecorr``. -Requires input file with xi+ and xi- previously -computed with ``treecorr``. - -:Author: Martin Kilbinger - -""" - -import sys -from optparse import OptionParser - -import numpy as np -import treecorr -from cs_util import logging - - -def params_default(): - - params = { - "input_path": "xip_xim.txt", - "theta_min": 0.5, - "theta_max": 200, - "n_theta": 20, - "n_Theta": 10, - "output_path": "./map2.txt", - } - - short_options = { - "input_path": "-i", - "output_path": "-o", - } - - types = { - "theta_min": "float", - "theta_max": "float", - "n_theta": "int", - "n_Theta": "int", - } - - help_strings = { - "input_path": "input file containing xi+ and xi-, default={}", - "theta_min": "mininum angular scale [arcmin], default={}", - "theta_max": "maximum angular scale [arcmin], default={}", - "n_theta": "number of angular scales on input, default={}", - "n_Theta": "number of angular scales on output, default={}", - "output_path": "output path, default={}", - } - - return params, short_options, types, help_strings - - -def parse_options(p_def, short_options, types, help_strings): - """Parse command line options. - - Parameters - ---------- - p_def : dict - default parameter values - help_strings : dict - help strings for options - - Returns - ------- - options: tuple - Command line options - """ - - usage = "%prog [OPTIONS]" - parser = OptionParser(usage=usage) - - for key in p_def: - if key in help_strings: - if key in short_options: - short = short_options[key] - else: - short = "" - - if key in types: - typ = types[key] - else: - typ = "string" - - parser.add_option( - short, - f"--{key}", - dest=key, - type=typ, - default=p_def[key], - help=help_strings[key].format(p_def[key]), - ) - - parser.add_option( - "-v", "--verbose", dest="verbose", action="store_true", help="verbose output" - ) - - options, args = parser.parse_args() - - return options - - -def main(argv=None): - - params, short_options, types, help_strings = params_default() - - options = parse_options(params, short_options, types, help_strings) - - # Update parameter values - for key in vars(options): - params[key] = getattr(options, key) - - # Save calling command - logging.log_command(argv) - - # Initizlies correlation object - gg = treecorr.GGCorrelation( - min_sep=params["theta_min"], - max_sep=params["theta_max"], - bin_size=params["n_theta"], - ) - - # Open input catalogue - if params["verbose"]: - print(f"Reading xi+ and xi- input file {params['input_path']}...") - gg.read(params["input_path"]) - - # Set up angular smoothing scales on output - R = np.geomspace( - params["theta_min"] * 5, - params["theta_max"] / 2, - params["n_Theta"], - ) - - # Compute correlation - if params["verbose"]: - print("Computing aperture mass dispersion...") - gg.calculateMapSq(R, m2_uform="Schneider") - - # Write to file - if params["verbose"]: - print(f"Writing output file {params['output_path']}") - gg.writeMapSq(params["output_path"], R=R, m2_uform="Schneider") - - return 0 - - -if __name__ == "__main__": - sys.exit(main(sys.argv)) diff --git a/scripts/xip_xim.py b/scripts/xip_xim.py deleted file mode 100755 index 763752ce..00000000 --- a/scripts/xip_xim.py +++ /dev/null @@ -1,217 +0,0 @@ -#!/usr/bin/env python3 - -"""xip_xim.py - -Compute shear correlation functions using ``treecorr``. - -:Author: Martin Kilbinger - -""" - -import sys -from optparse import OptionParser - -import treecorr -from astropy.io import fits -from cs_util import logging - - -def params_default(): - - params = { - "input_path": "shape_catalog_ngmix.fits", - "tomo_bin1": None, - "tomo_bin2": None, - "key_ra": "RA", - "key_dec": "DEC", - "key_e1": "e1", - "key_e2": "e2", - "sign_e1": +1, - "sign_e2": +1, - "theta_min": 0.5, - "theta_max": 200, - "n_theta": 20, - "output_path": "./xip_xim.txt", - "key_tomo_bin_col": "tomo_bin_id", - } - - short_options = { - "input_path": "-i", - "output_path": "-o", - "tomo_bin1": "-b1", - "tomo_bin2": "-b2", - } - - types = { - "sign_e1": "int", - "sign_e2": "int", - "tomo_bin1": "int", - "tomo_bin2": "int", - } - - help_strings = { - "input_path": "shear catalogue input path, default={}", - "tomo_bin1": "First tomographic bin, default={} (non-tomographic)", - "tomo_bin2": "Second tomographic bin, default={} (non-tomographic)", - "key_ra": "column name for right ascension, default={}", - "key_dec": "column name for declination, default={}", - "key_e1": "column name for ellipticity component 1, default={}", - "key_e2": "column name for ellipticity component 2, default={}", - "sign_e1": "sign for ellipticity component 1, default={}", - "sign_e2": "sign for ellipticity component 2, default={}", - "theta_min": "mininum angular scale [arcmin], default={}", - "theta_max": "maximum angular scale [arcmin], default={}", - "n_theta": "number of angular scales, default={}", - "output_path": "output path, default={}", - "key_tomo_bin_col": "column name for tomography bin, default={}", - } - - return params, short_options, types, help_strings - - -def parse_options(p_def, short_options, types, help_strings): - """Parse command line options. - - Parameters - ---------- - p_def : dict - default parameter values - help_strings : dict - help strings for options - - Returns - ------- - options: tuple - Command line options - """ - - usage = "%prog [OPTIONS]" - parser = OptionParser(usage=usage) - - for key in p_def: - if key in help_strings: - if key in short_options: - short = short_options[key] - else: - short = "" - - if key in types: - typ = types[key] - else: - typ = "string" - - parser.add_option( - short, - f"--{key}", - dest=key, - type=typ, - default=p_def[key], - help=help_strings[key].format(p_def[key]), - ) - - parser.add_option( - "-v", "--verbose", dest="verbose", action="store_true", help="verbose output" - ) - - options, args = parser.parse_args() - - return options - - -def main(argv=None): - - params, short_options, types, help_strings = params_default() - - options = parse_options(params, short_options, types, help_strings) - - # Update parameter values - for key in vars(options): - params[key] = getattr(options, key) - - # Save calling command - logging.log_command(argv) - - # Open input catalogue - if params["verbose"]: - print(f"Reading catalogue {params['input_path']}...") - data = fits.getdata(params["input_path"]) - - coord_units = "degrees" - if params["verbose"]: - print( - "Signs for ellipticity components =" - + f" ({params['sign_e1']:+d}, {params['sign_e2']:+d})" - ) - tomographic = params["tomo_bin1"] is not None and params["tomo_bin2"] is not None - cat2 = None - if tomographic: - if params["verbose"]: - print( - f"Calculating 2PCF for tomographic bins {params['tomo_bin1']} and {params['tomo_bin2']}" - ) - tomo_bin1_idx = data[params["key_tomo_bin_col"]] == params["tomo_bin1"] - tomo_bin2_idx = data[params["key_tomo_bin_col"]] == params["tomo_bin2"] - - g1 = data[params["key_e1"]] * params["sign_e1"] - g2 = data[params["key_e2"]] * params["sign_e2"] - - cat1 = treecorr.Catalog( - ra=data[params["key_ra"]][tomo_bin1_idx], - dec=data[params["key_dec"]][tomo_bin1_idx], - g1=g1[tomo_bin1_idx], - g2=g2[tomo_bin1_idx], - w=data["w"][tomo_bin1_idx], - ra_units=coord_units, - dec_units=coord_units, - ) - if params["tomo_bin1"] != params["tomo_bin2"]: - cat2 = treecorr.Catalog( - ra=data[params["key_ra"]][tomo_bin2_idx], - dec=data[params["key_dec"]][tomo_bin2_idx], - g1=g1[tomo_bin2_idx], - g2=g2[tomo_bin2_idx], - w=data["w"][tomo_bin2_idx], - ra_units=coord_units, - dec_units=coord_units, - ) - else: - if params["verbose"]: - print("Calculating non-tomographic 2PCF") - g1 = data[params["key_e1"]] * params["sign_e1"] - g2 = data[params["key_e2"]] * params["sign_e2"] - cat1 = treecorr.Catalog( - ra=data[params["key_ra"]], - dec=data[params["key_dec"]], - g1=g1, - g2=g2, - w=data["w"], - ra_units=coord_units, - dec_units=coord_units, - ) - # Set treecorr config info for correlation - sep_units = "arcmin" - TreeCorrConfig = { - "ra_units": coord_units, - "dec_units": coord_units, - "min_sep": params["theta_min"], - "max_sep": params["theta_max"], - "sep_units": sep_units, - "nbins": params["n_theta"], - } - gg = treecorr.GGCorrelation(TreeCorrConfig) - - # Compute correlation - if params["verbose"]: - print("Correlating...") - gg.process(cat1, cat2=cat2) - - # Write to file - if params["verbose"]: - print(f"Writing output file {params['output_path']}") - gg.write(params["output_path"]) - - return 0 - - -if __name__ == "__main__": - sys.exit(main(sys.argv)) diff --git a/src/sp_validation/blinding.py b/src/sp_validation/blinding.py new file mode 100644 index 00000000..09c49718 --- /dev/null +++ b/src/sp_validation/blinding.py @@ -0,0 +1,350 @@ +"""Blinds, the theory a blind shifts by, and the call that conceals a SACC. + +A blind is a named secret, ``/.blind.json`` holding a +seed, the envelope the hidden point is drawn in and the fiducial point. It is +drawn once, ``python -m sp_validation.blinding init ``, and never +committed (``*.blind.json`` is gitignored). Each catalogue entry of the +catalogue config declares ``blind: none`` (public) or ``blind: `` +(:func:`blind_of`). + +:func:`conceal` adds the blind's shift, t(hidden) − t(fiducial), to every row +of a SACC; Smokescreen draws the hidden point from the seed and evaluates the +theory t at both points. + +A theory is a function ``theory(params, s)``: + +- ``params`` is a dict of cosmological parameters with the keys of + :func:`fiducial` (``S8``, ``Omega_m``, ``Omega_b``, ``h``, ``n_s``, + ``m_nu``, ``w0``, ``wa``, ``logT_AGN``, ``A_IA``). The fiducial and hidden + points differ only in ``S8`` and ``Omega_m``. +- ``s`` is the ``sacc.Sacc`` being concealed: its rows carry a data type, a + tracer pair and a θ or ℓ tag, its tracers the n(z). +- It returns an array shaped like the data, ``len(s.mean)`` values in the + order of ``s.data``: the prediction for each row, zero where the cosmology + has no effect. + +For Cℓ_EE/BB/EB bandpowers, for example:: + + def my_theory(params, s): + out = np.zeros(len(s.mean)) + ee = s.indices(sacc_io.CL_EE) + window = s.get_bandpower_windows(ee) + out[ee] = window.weight.T @ my_cl_ee(params, window.values) + return out # Cℓ_BB and Cℓ_EB stay zero + +:func:`shear` is the default theory; :func:`no_signal` is that of a statistic +with no cosmological signal (ρ/τ). + +The host Snakemake loads this module by path, so at import it needs only the +standard library. +""" + +import argparse +import dataclasses +import functools +import json +import math +import os +import re +import secrets +from pathlib import Path + +ENVELOPE = {"S8": 0.075, "Omega_m": 0.1} +REPO_CAT_CONFIG = Path(__file__).parents[2] / "cosmo_val" / "cat_config.yaml" +_NAME = re.compile(r"[a-z0-9][a-z0-9_.-]*") + + +class BlindingError(RuntimeError): + """A blind cannot be opened or applied; the message carries no hidden value.""" + + +@dataclasses.dataclass(frozen=True) +class Blind: + """A blind's name and the path of its record; ``none`` has no record.""" + + name: str + path: Path | None = None + + def record(self): + """The blind's ``{seed, envelope, fiducial}``.""" + return json.loads(self.path.read_text()) + + +NONE = Blind("none") + + +def _shear_file(entry): + shear = entry.get("shear") if isinstance(entry, dict) else None + if not isinstance(shear, dict) or "path" not in shear: + return None + return os.path.normpath(os.path.join(str(entry.get("subdir", "")), shear["path"])) + + +def blind_of(catalogues, entry): + """The ``blind:`` of catalogue-config entry ``entry``. + + Refuses an entry that declares none, and entries reading one shear file + under different blinds. + """ + blind = catalogues[entry].get("blind") + if not isinstance(blind, str): + raise ValueError( + f"{entry} declares no `blind:`; declare `blind: none` (public) or " + "`blind: ` in its catalogue config entry" + ) + path = _shear_file(catalogues[entry]) + readers = { + name: other.get("blind") + for name, other in catalogues.items() + if path and _shear_file(other) == path + } + if len(set(readers.values())) > 1: + raise ValueError(f"{path} is read under different blinds: {readers}") + return blind + + +def registry(catalogues): + """The directory a catalogue config keeps its blinds in (``paths.blinds``).""" + return Path(os.path.expanduser(catalogues["paths"]["blinds"])) + + +def open_blind(name, catalogues): + """Blind ``name`` of the catalogue config ``catalogues``; its record must exist.""" + if name == NONE.name: + return NONE + path = registry(catalogues) / f"{name}.blind.json" + if not os.access(path, os.R_OK): + raise BlindingError( + f"cannot read blind {name} at {path}; `python -m sp_validation.blinding " + f"init {name}` draws it, if it is meant to be new" + ) + return Blind(name, path) + + +def conceal(s, blind, theory=None): + """A copy of SACC ``s`` with ``blind``'s shift added to every row. + + The shift is t(hidden) − t(fiducial), t being ``theory`` (default + :func:`shear`). A failure raises :class:`BlindingError` naming only the + exception type, so no message or traceback shows the hidden point. + """ + import numpy as np + from smokescreen.datavector import concealing_factor + + theory = theory or shear + record, failure = blind.record(), None + try: + shift = concealing_factor( + record["fiducial"], + record["envelope"], + seed=record["seed"], + theory_fn=lambda params: theory(params, s), + ) + except Exception as err: + failure = type(err).__name__ + if failure: # raised outside the except, so nothing is chained + raise BlindingError(f"the theory failed under blind {blind.name} ({failure})") + shift = np.asarray(shift, float) + if shift.shape != (len(s.mean),): + raise BlindingError( + f"the theory returned {shift.shape[0] if shift.ndim else 0} values for " + f"{len(s.mean)} rows" + ) + out = s.copy() + for dp, d in zip(out.data, shift): + dp.value += float(d) + return out + + +def init(name, catalogues): + """Draw blind ``name`` into the catalogue config's registry, once.""" + if not _NAME.fullmatch(name) or name == NONE.name: + raise BlindingError(f"blind name {name!r}: use a-z, 0-9, - _ .") + where = registry(catalogues) + where.mkdir(mode=0o700, parents=True, exist_ok=True) + record = { + "seed": secrets.token_hex(32), + "envelope": ENVELOPE, + "fiducial": fiducial(), + } + path = where / f"{name}.blind.json" + try: + fd = os.open(path, os.O_WRONLY | os.O_CREAT | os.O_EXCL, 0o440) + except FileExistsError: + raise BlindingError(f"blind {name} exists at {path}") from None + with os.fdopen(fd, "w") as f: + json.dump(record, f, indent=1) + print(f"[blinding] drew blind {name} at {path}") + return Blind(name, path) + + +# --------------------------------------------------------------------------- # +# Theories +# --------------------------------------------------------------------------- # +_COSMOLOGY = ("S8", "Omega_m", "Omega_b", "h", "n_s", "m_nu", "w0", "wa", "logT_AGN") + + +def fiducial(): + """cs_util's Planck 2018 point, with feedback and no IA.""" + from cs_util.cosmo import PLANCK18 as p + + return { + "S8": float(p["sigma_8"] * math.sqrt(p["Omega_m"] / 0.3)), + **{ + k: float(p[k]) + for k in ("Omega_m", "Omega_b", "h", "n_s", "m_nu", "w0", "wa") + }, + "logT_AGN": 7.5, + "A_IA": 0.0, + } + + +def no_signal(params, s): + """Zeros: the theory of a statistic with no cosmological signal (ρ/τ).""" + import numpy as np + + return np.zeros(len(s.mean)) + + +def shear(params, s): + """Cosmic-shear ξ± and Cℓ_EE from pyccl; zeros for Cℓ_BB and Cℓ_EB. + + ξ± is evaluated at each row's stored θ (TreeCorr's ``meanr``) and Cℓ_EE + through each row's bandpower window. Any other data type raises: pass your + own theory. + """ + import numpy as np + + from .sacc_io import CL_BB, CL_EB, CL_EE, XI_MINUS, XI_PLUS + + out = np.zeros(len(s.mean)) + groups = {} + for i, dp in enumerate(s.data): + if dp.data_type in (XI_PLUS, XI_MINUS): + key = (_xi, dp.tracers) + elif dp.data_type == CL_EE: + key = (_cl, dp.tracers, id(dp.tags.get("window"))) + elif dp.data_type in (CL_BB, CL_EB): + continue + else: + raise ValueError( + f"blinding.shear has no prediction for {dp.data_type}; pass your " + "own theory" + ) + groups.setdefault(key, []).append(i) + # One call per tracer pair (and window), so ξ+ and ξ− share their Cℓ. + for (function, *_), rows in groups.items(): + out[rows] = function(params, s, np.asarray(rows)) + return out + + +def cosmology(params): + """The ``pyccl.Cosmology`` at ``params``, built once per point. + + Built through ``cs_util.cosmo.get_cosmo`` on the CAMB HMCode2020-feedback + route the CosmoSIS inference runs, with σ8 = S8/√(Ωm/0.3); ``get_cosmo`` + subtracts CAMB's Ω_ν from ``Omega_m`` for Ω_c. + """ + return _cosmology(tuple(float(params[k]) for k in _COSMOLOGY)) + + +@functools.lru_cache(maxsize=4) +def _cosmology(point): + from cs_util.cosmo import get_cosmo + + p = dict(zip(_COSMOLOGY, point)) + cosmo = get_cosmo( + Omega_m=p["Omega_m"], + Omega_b=p["Omega_b"], + h=p["h"], + sig8=p["S8"] / math.sqrt(p["Omega_m"] / 0.3), + ns=p["n_s"], + w0=p["w0"], + wa=p["wa"], + mnu=p["m_nu"], + # get_cosmo's default "halofit" ignores extra_params. + matter_power_spectrum="camb", + extra_params={ + "camb": { + "halofit_version": "mead2020_feedback", + "HMCode_logT_AGN": p["logT_AGN"], + } + }, + ) + # CCL's Hankel transform extrapolates C_ℓ to ELL_MAX_CORR; its default + # (6·10⁴) rings in ξ− below a few arcmin. + cosmo.cosmo.spline_params.ELL_MAX_CORR = 10_000_000 + cosmo.cosmo.spline_params.N_ELL_CORR = 5_000 + return cosmo + + +@functools.cache +def _ell(): + """The multipoles the ξ± Hankel transform integrates over: every ℓ below 50, + then 200 log-spaced up to 6·10⁴.""" + import numpy as np + + return np.unique(np.concatenate([np.arange(2, 50), np.geomspace(50, 6e4, 200)])) + + +def _lensing(cosmo, params, s, name): + import numpy as np + import pyccl as ccl + + tracer = s.tracers[name] + if getattr(tracer, "nz", None) is None: + raise ValueError(f"tracer {name} has no n(z)") + z, n = np.asarray(tracer.z, float), np.asarray(tracer.nz, float) + ia = None if params["A_IA"] == 0 else (z, np.full_like(z, params["A_IA"])) + return ccl.WeakLensingTracer(cosmo, dndz=(z, n), ia_bias=ia) + + +def _shear_cl(params, s, rows, ell): + import pyccl as ccl + + cosmo = cosmology(params) + a, b = (_lensing(cosmo, params, s, t) for t in s.data[rows[0]].tracers) + return ccl.angular_cl(cosmo, a, b, ell) + + +def _xi(params, s, rows): + """ξ± of one shear pair: Limber Cℓ on :func:`_ell`, then CCL's Hankel transform.""" + import numpy as np + import pyccl as ccl + + from .sacc_io import XI_PLUS + + ell = _ell() + cosmo, cl = cosmology(params), _shear_cl(params, s, rows, ell) + theta = np.array([s.data[i].tags["theta"] for i in rows]) / 60.0 + xip, xim = ( + ccl.correlation(cosmo, ell=ell, C_ell=cl, theta=theta, type=t) + for t in ("GG+", "GG-") + ) + plus = np.array([s.data[i].data_type == XI_PLUS for i in rows]) + return np.where(plus, xip, xim) + + +def _cl(params, s, rows): + """Cℓ_EE of one shear pair through its bandpower window.""" + import numpy as np + + window = s.get_bandpower_windows(rows) + cl = _shear_cl(params, s, rows, np.asarray(window.values, float)) + return np.asarray(window.weight).T @ cl + + +def main(argv=None): + import yaml + + parser = argparse.ArgumentParser(prog="python -m sp_validation.blinding") + parser.add_argument("command", choices=("init",)) + parser.add_argument("blind") + parser.add_argument("--cat-config", default=str(REPO_CAT_CONFIG)) + a = parser.parse_args(argv) + init(a.blind, yaml.safe_load(Path(a.cat_config).read_text())) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py index 4fbd2a89..8159a1a0 100644 --- a/src/sp_validation/cosmo_val/core.py +++ b/src/sp_validation/cosmo_val/core.py @@ -2,6 +2,7 @@ import copy import itertools import os +import re from pathlib import Path import colorama @@ -10,7 +11,7 @@ from cs_util.cosmo import get_cosmo from shear_psf_leakage import leakage, run_object, run_scale -from .. import io +from .. import blinding, io from ..b_modes import ( _get_pte_from_scale_cut, covariance_label, @@ -386,6 +387,9 @@ def __init__( self.catalog_config_path = Path(catalog_config) with self.catalog_config_path.open("r") as file: self.cc = cc = yaml.load(file, Loader=yaml.FullLoader) + # The config as written, before virtual versions and resolved paths: + # what each catalogue's blind is checked against. + self._config_as_read = copy.deepcopy(cc) def resolve_paths_for_version(ver): """Resolve relative paths for a version using its subdir.""" @@ -467,6 +471,8 @@ def ensure_version_exists(ver): # B-mode results storage for summarize_bmodes() self._pure_eb_results = {} self._cosebis_results = {} + # The sealed ξ± parts calculate_2pcf_version made, by (version, grid). + self.xi_parts = {} def _output_path(self, *parts): """Absolute path under the catalog config's output directory. @@ -578,6 +584,36 @@ def results_objectwise(self): self._results_objectwise = self.init_results(objectwise=True) return self._results_objectwise + def blind(self, version): + """The blind every SACC this object writes for ``version`` is saved + under, as the catalogue config declares it.""" + entry = version + while entry not in self._config_as_read: + base = re.sub(r"_leak_corr$", "", entry) + if base == entry: + base = self._split_seed_variant(entry)[0] or entry + if base == entry: + raise ValueError(f"no catalogue {version} in the catalogue config") + entry = base + catalogues = self._config_as_read + return blinding.open_blind(blinding.blind_of(catalogues, entry), catalogues) + + def _refuse_blinded_tomography(self, version): + """Refuse to measure tomographic shear signal on a blinded catalogue. + + Only the ``("all", "all")`` pair is born as a SACC part, whose tracer + carries the n(z) the blind's theory shifts it by. The tomographic ξ±, + ⟨M_ap²⟩ and pseudo-Cℓ are TreeCorr objects and FITS tables with no + per-bin tracer, so on a blinded catalogue they would hold true values. + """ + blind = self.blind(version) + if blind.name != blinding.NONE.name: + raise blinding.BlindingError( + f"{version} is blinded under {blind.name}, and blinding covers " + "its non-tomographic pair only; tomographic shear statistics are " + "measured on public (`blind: none`) catalogues" + ) + def basename( self, version, @@ -634,8 +670,7 @@ def _calibrated_g(self, ver): Applies additive-bias subtraction and the multiplicative response: ``g = (e − c) / R``. For DES the response is the catalog-averaged per-component ``R11``/``R22`` (column names in the config); for every - other version it is the scalar ``R`` from the config. Used identically - by :meth:`calculate_2pcf` and :meth:`calculate_aperture_mass_dispersion`. + other version it is the scalar ``R`` from the config. Must be called inside a ``self.results[ver].temporarily_read_data()`` context, since it reads ``dat_shear`` columns. diff --git a/src/sp_validation/cosmo_val/cosebis.py b/src/sp_validation/cosmo_val/cosebis.py index 079720bf..bea1df16 100644 --- a/src/sp_validation/cosmo_val/cosebis.py +++ b/src/sp_validation/cosmo_val/cosebis.py @@ -10,6 +10,7 @@ from ..b_modes import ( calculate_cosebis, find_conservative_scale_cut_key, + log_bin_edges, plot_cosebis_covariance_matrix, plot_cosebis_modes, plot_cosebis_scale_cut_heatmap, @@ -98,13 +99,8 @@ def calculate_cosebis( "tomographic bin pair." ) - # Set up parameters with defaults - npatch = npatch or self.npatch - - # Always use integration binning for COSEBIs calculation (fine binning) - treecorr_config = self._binning(min_sep_int, max_sep_int, nbins_int) - - # Calculate single fine-binned correlation function for COSEBIs + # Calculate single fine-binned correlation function for COSEBIs, as the + # sealed part: on a blinded catalogue the COSEBIs are concealed with it. print( f"Computing fine-binned 2PCF with {nbins_int} bins from {min_sep_int} to " f"{max_sep_int} arcmin" @@ -114,7 +110,10 @@ def calculate_cosebis( version, npatch=npatch, compute_tomography=compute_tomography, - **treecorr_config, + grid="integration", + min_sep=min_sep_int, + max_sep=max_sep_int, + nbins=nbins_int, ) results = { bin_key: None for bin_key in ggs.keys() @@ -270,7 +269,6 @@ def plot_cosebis( else: results_tomo = {"tomo_bin_all_tomo_bin_all": results} - ggs_temp = None for bin_key, bin_results in results_tomo.items(): bin_stub = ( out_stub @@ -300,18 +298,12 @@ def plot_cosebis( ) if multiple_scale_cuts and len(bin_results) > 1: - if ggs_temp is None: - treecorr_config_temp = self._binning(min_sep, max_sep, nbins) - ggs_temp = self.calculate_2pcf_version( - version, - npatch=npatch, - compute_tomography=compute_tomography, - **treecorr_config_temp, - ) - gg_temp = ggs_temp[bin_key] + binning = self._binning(min_sep, max_sep, nbins) plot_cosebis_scale_cut_heatmap( bin_results, - (gg_temp.left_edges, gg_temp.right_edges), + log_bin_edges( + binning["min_sep"], binning["max_sep"], binning["nbins"] + ), version, bin_stub + "_scalecut_ptes.png", fiducial_scale_cut=fiducial_scale_cut, diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index 686cda01..a728d3eb 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -84,6 +84,12 @@ def plot_pseudo_cl_spectrum(datasets, spectrum, output_path): plt.close(fig) +def _spectra(part): + """The ELL/EE/EB/BE/BB dict of a single-field pseudo-Cl part, whose BE is EB.""" + ell, ee, bb, eb, _window = sacc_io.get_pseudo_cl(part, SACC_BIN) + return {"ELL": ell, "EE": ee, "EB": eb, "BE": eb.copy(), "BB": bb} + + class PseudoClMixin: # ---------------- Pseudo-Cl properties ---------------- # @property @@ -114,10 +120,14 @@ def calculate_pseudo_cl(self, compute_tomography=True, out_path=None): Compute the pseudo-Cl of a `CosmologyValidation` inputs with tomography. With ``compute_tomography=False`` the single ``("all", "all")`` pair is - computed and written as a SACC part, by default to - ``pseudo_cl_{ver}.sacc``; ``out_path`` overrides that destination (one - version only). With ``compute_tomography=True`` every tomographic bin - pair is written to its own FITS file and ``out_path`` is not accepted. + computed and written as a SACC part sealed under the version's blind, + by default to ``pseudo_cl_{ver}.sacc``; ``out_path`` overrides that + destination (one version only). An existing part is reused only if it + was sealed under the blind declared now. With + ``compute_tomography=True`` every tomographic bin pair is written to its + own FITS file and ``out_path`` is not accepted; tomographic pairs are + measured only on a public catalogue + (:meth:`_refuse_blinded_tomography`). """ if out_path is not None: if compute_tomography: @@ -149,6 +159,7 @@ def calculate_pseudo_cl(self, compute_tomography=True, out_path=None): self._pseudo_cls[ver] = {} if compute_tomography: + self._refuse_blinded_tomography(ver) tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver) if tomo_bin_ids is None or tomo_bin_pairs is None: @@ -174,13 +185,18 @@ def calculate_pseudo_cl(self, compute_tomography=True, out_path=None): pair_out_path = out_path or self._output_path_pseudo_cl( ver, tomo_bin_pair=(bin_key1, bin_key2) ) - if os.path.exists(pair_out_path) and not self.force_run: + held = ( + self._load_pseudo_cl(ver, pair_out_path, (bin_key1, bin_key2)) + if os.path.exists(pair_out_path) and not self.force_run + else None + ) + if held is not None: self.print_done( f"Skipping Pseudo-Cl's calculation, {pair_out_path} exists" ) self._pseudo_cls[ver][f"tomo_bin_{bin_key1}_tomo_bin_{bin_key2}"][ "pseudo_cl" - ] = self._load_pseudo_cl(pair_out_path, (bin_key1, bin_key2)) + ] = held continue if self.cell_method == "map": @@ -289,13 +305,11 @@ def calculate_pseudo_cl_map(self, ver, nside, out_path, tomo_bin_a, tomo_bin_b): self.print_cyan("Saving pseudo-Cl's...") tomo_bin_pair = (tomo_bin_a, tomo_bin_b) - self._save_pseudo_cl( - ver, out_path, tomo_bin_pair, ell_eff, cl_shear, wsp, nside=nside - ) - self._pseudo_cls[ver][f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ "pseudo_cl" - ] = self._load_pseudo_cl(out_path, tomo_bin_pair) + ] = self._save_pseudo_cl( + ver, out_path, tomo_bin_pair, ell_eff, cl_shear, wsp, nside=nside + ) def calculate_pseudo_cl_catalog(self, ver, out_path, tomo_bin_a, tomo_bin_b): assert (tomo_bin_a == "all" and tomo_bin_b == "all") or ( @@ -314,11 +328,9 @@ def calculate_pseudo_cl_catalog(self, ver, out_path, tomo_bin_a, tomo_bin_b): self.print_cyan("Saving pseudo-Cl's...") tomo_bin_pair = (tomo_bin_a, tomo_bin_b) - self._save_pseudo_cl(ver, out_path, tomo_bin_pair, ell_eff, cl_shear, wsp) - self._pseudo_cls[ver][f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ "pseudo_cl" - ] = self._load_pseudo_cl(out_path, tomo_bin_pair) + ] = self._save_pseudo_cl(ver, out_path, tomo_bin_pair, ell_eff, cl_shear, wsp) def calculate_pseudo_cl_inka_cov( self, compute_tomography=True, load_all_block=False @@ -1206,28 +1218,32 @@ def _output_path_iNKA_block_cov(self, ver, tomo_bin_quad): def _save_pseudo_cl( self, ver, out_path, tomo_bin_pair, ell_eff, cl_all, wsp, nside=None ): - """Write one pair as SACC or FITS, forwarding ``nside`` for map spectra.""" - if tuple(tomo_bin_pair) == ("all", "all"): - self.pseudo_cl_to_sacc_part( - ver, out_path, ell_eff, cl_all, wsp, nside=nside - ) - else: - self.save_pseudo_cl(ell_eff, cl_all, out_path) + """Write one pair as SACC or FITS and return its spectra as written. - def _load_pseudo_cl(self, out_path, tomo_bin_pair): - """Read one bin pair's pseudo-Cl product written by ``_save_pseudo_cl``.""" + Forwards ``nside`` for map spectra. The ``("all", "all")`` spectra are + those of the part as sealed under the version's blind. + """ if tuple(tomo_bin_pair) == ("all", "all"): - return self._load_pseudo_cl_sacc(out_path) + return _spectra( + self.pseudo_cl_to_sacc_part( + ver, out_path, ell_eff, cl_all, wsp, nside=nside + ) + ) + self.save_pseudo_cl(ell_eff, cl_all, out_path) return fits.getdata(out_path) - @staticmethod - def _load_pseudo_cl_sacc(out_path): - """Read the single-field auto-spectrum, whose BE equals EB.""" - # Readback of a part this producer just wrote — a legitimate pre-blind - # consumer, so the fail-closed load is opted out of. - s = sacc_io.load(out_path, allow_unblinded=True) - ell, ee, bb, eb, _window = sacc_io.get_pseudo_cl(s, SACC_BIN) - return {"ELL": ell, "EE": ee, "EB": eb, "BE": eb.copy(), "BB": bb} + def _load_pseudo_cl(self, ver, out_path, tomo_bin_pair): + """One bin pair's pseudo-Cl product written by ``_save_pseudo_cl``. + + ``None`` for an ``("all", "all")`` part sealed under a blind other than + the one ``ver`` declares now, which must be measured again. + """ + if tuple(tomo_bin_pair) != ("all", "all"): + return fits.getdata(out_path) + part = sacc_io.load(out_path) + if sacc_io.stamp(part) != self.blind(ver).name: + return None + return _spectra(part) def pseudo_cl_to_sacc_part( self, version, out_path, ell_eff, cl_all, wsp, nside=None @@ -1237,7 +1253,8 @@ def pseudo_cl_to_sacc_part( ``cl_all`` is NaMaster's decoupled ``(4, nbp)`` array (EE, EB, BE, BB). Map-based callers pass ``nside`` so the window includes ``pw²(ℓ)``; catalogue-based callers leave it unset because their spectra have no - HEALPix pixel window. No covariance is attached here. + HEALPix pixel window. No covariance is attached here. Returns the part + as sealed under the version's blind. """ s = pseudo_cl_to_sacc( self.sacc_nz(version), @@ -1247,7 +1264,7 @@ def pseudo_cl_to_sacc_part( wsp, nside=nside, ) - sacc_io.save(s, out_path, type="data") + return sacc_io.save(s, out_path, blind=self.blind(version)) def save_pseudo_cl(self, ell_eff, pseudo_cl, out_path): """ diff --git a/src/sp_validation/cosmo_val/psf_systematics.py b/src/sp_validation/cosmo_val/psf_systematics.py index b189e622..991567fb 100644 --- a/src/sp_validation/cosmo_val/psf_systematics.py +++ b/src/sp_validation/cosmo_val/psf_systematics.py @@ -18,7 +18,7 @@ from shear_psf_leakage.rho_tau_stat import PSFErrorFit from uncertainties import ufloat -from .. import sacc_io +from .. import blinding, sacc_io from ..io import open_entry from ..rho_tau import ( get_rho_tau_w_cov, @@ -134,7 +134,7 @@ def calculate_rho_tau_stats(self, tomography=True): def rho_tau_to_sacc_part( self, version, out_dir, base, rho_stat_handler, tau_stat_handler ): - """Write the ρ/τ SACC part for one version. + """Write the ρ/τ SACC part for one version, under the version's blind. ρ_0…ρ_5 autos and τ_0/τ_2/τ_5 leakage from the handler tables. The ``CovTauTh`` theory covariance ``cov_tau_{base}_th.npy`` is passed as @@ -157,7 +157,7 @@ def rho_tau_to_sacc_part( tau_cov_th=tau_cov_th, ) out_path = os.path.join(out_dir, f"rho_tau_{base}.sacc") - sacc_io.save(s, out_path, type="data") + sacc_io.save(s, out_path, blind=self.blind(version), theory=blinding.no_signal) def calculate_rho_tau_fits(self, tomography=True, track_result=True): assert self.rho_tau_method != "none" diff --git a/src/sp_validation/cosmo_val/pure_eb.py b/src/sp_validation/cosmo_val/pure_eb.py index 503b02e8..793fa826 100644 --- a/src/sp_validation/cosmo_val/pure_eb.py +++ b/src/sp_validation/cosmo_val/pure_eb.py @@ -36,10 +36,12 @@ def calculate_pure_eb( """ Calculate the pure E/B modes for the given catalog version. - ξ± is measured on the fine integration grid only; the reporting - binning (the instance's treecorr_config unless overridden) enters as - bin edges, snapped onto the fine edges, into which - :func:`~sp_validation.b_modes.pure_eb_operator` averages the modes. + ξ± is measured on the fine integration grid only, as the version's + sealed part (:meth:`calculate_2pcf_version`), so a blinded catalogue's + modes are concealed. The reporting binning (the instance's treecorr_config + unless overridden) enters as bin edges, snapped onto the fine edges, + into which :func:`~sp_validation.b_modes.pure_eb_operator` averages + the modes. Parameters ---------- @@ -57,8 +59,8 @@ def calculate_pure_eb( instead of the non-tomographic pair. Defaults to False. cov_path_int : str, optional Analytic ξ± covariance on the integration grid. Without it the - covariance is the jackknife of the integration-grid ξ±. One file - holds one pair's covariance, so it cannot be combined with + covariance is the jackknife the integration-grid part carries. One + file holds one pair's covariance, so it cannot be combined with ``compute_tomography``. Returns @@ -79,6 +81,10 @@ def calculate_pure_eb( "tomographic bin pair." ) + integration = self._binning(min_sep_int, max_sep_int, nbins_int) + integration_edges = np.geomspace( + integration["min_sep"], integration["max_sep"], integration["nbins"] + 1 + ) reporting = self._binning(min_sep, max_sep, nbins) reporting_edges = np.geomspace( reporting["min_sep"], reporting["max_sep"], reporting["nbins"] + 1 @@ -87,7 +93,10 @@ def calculate_pure_eb( version, npatch=npatch, compute_tomography=compute_tomography, - **self._binning(min_sep_int, max_sep_int, nbins_int), + grid="integration", + min_sep=integration["min_sep"], + max_sep=integration["max_sep"], + nbins=integration["nbins"], ) results = {} @@ -95,14 +104,13 @@ def calculate_pure_eb( if cov_path_int is not None: cov_xi, cov_npatch = np.loadtxt(cov_path_int), None else: - cov_xi = gg_int.estimate_cov("jackknife", cross_patch_weight="match") - cov_npatch = gg_int.npatch1 + cov_xi, cov_npatch = gg_int.cov, gg_int.npatch1 results[bin_key] = calculate_pure_eb_correlation( gg_int.meanr, gg_int.xip, gg_int.xim, gg_int.weight, - np.append(gg_int.left_edges, gg_int.right_edges[-1]), + integration_edges, cov_xi, reporting_edges, npatch=cov_npatch, diff --git a/src/sp_validation/cosmo_val/real_space.py b/src/sp_validation/cosmo_val/real_space.py index 7dea3542..8678b828 100644 --- a/src/sp_validation/cosmo_val/real_space.py +++ b/src/sp_validation/cosmo_val/real_space.py @@ -13,6 +13,11 @@ from sp_validation.statistics import jackknife_patch_centers +from .. import sacc_io +from .sacc_writers import xi_to_sacc + +ALL_PAIR = "tomo_bin_all_tomo_bin_all" + class RealSpaceMixin: def calculate_2pcf_version( @@ -20,52 +25,49 @@ def calculate_2pcf_version( ver, npatch=None, compute_tomography=False, + *, + grid="reporting", + out=None, **treecorr_config, ): - """ - Calculate the two-point correlation function (2PCF) ξ± for a single catalog - version with TreeCorr. - - This is the per-version child function. Use :meth:`calculate_2pcf` to run over every - version in ``self.versions`` in one call. - - By default the class instance's `npatch` and `treecorr_config` entries are - used to initialize the TreeCorr Catalog and GGCorrelation objects, but may - be overridden by passing keyword arguments. + """ξ± of one catalogue version on one binning, for each bin pair. + + The non-tomographic ``("all", "all")`` pair is born as its SACC part, + sealed under the version's blind (:func:`sp_validation.sacc_io.seal`) + before it is kept in ``self.xi_parts[ver, grid]``, returned or written, + so on a blinded catalogue its true values never leave this method. A + part already held for ``(ver, grid)`` is reused when no ``npatch``, + binning or ``out`` is asked for, so the figure rules draw the parts they + are handed. Tomographic pairs are measured only on a public catalogue + (:meth:`_refuse_blinded_tomography`). Use :meth:`calculate_2pcf` to run + over every version in ``self.versions``. Parameters: - ver (str): The catalog version to process. - - npatch (int, optional): The number of patches to use for the - calculation. Defaults to the instance's `npatch` attribute. - - compute_tomography (bool, optional): Whether to compute tomographic - correlations. Defaults to False. - - **treecorr_config: Additional TreeCorr configuration parameters that - will override the instance's default `treecorr_config`. For example, - `min_sep=1`. + ver (str): The catalogue version to measure. + npatch (int, optional): Jackknife patches; the instance's + ``npatch`` by default. With patches the part carries the + jackknife covariance, without them the shot-noise diagonal. + Seeded patch centres come from the full catalogue and are + shared by every bin pair. + compute_tomography (bool, optional): Measure the tomographic bin + pairs instead of ``("all", "all")``. + grid (str): The grid tag the part's points carry. + out (str, optional): Where to write the part as well. + **treecorr_config: Overrides of the instance's ``treecorr_config``, + e.g. ``min_sep=1``. Returns: - dict: Mapping of ``"tomo_bin_{b1}_tomo_bin_{b2}"`` to the corresponding - treecorr.GGCorrelation object. For the non-tomographic case the single - key is ``"tomo_bin_all_tomo_bin_all"``. - - Notes: - - The non-tomographic pair is written to the columns-only TreeCorr - dump ``xi_{basename}.txt``. If that file already exists, the pair is - read back from it instead of being recomputed. - - Seeded patch centres are computed once from the full catalogue and - shared by every tomographic bin pair. + dict: ``"tomo_bin_{b1}_tomo_bin_{b2}"`` → ξ± shaped like a TreeCorr + ``GGCorrelation``: the part's + :func:`~sp_validation.sacc_io.xi_correlation` view for + ``("all", "all")``, the TreeCorr object for a tomographic pair. """ - - npatch = npatch or self.npatch - treecorr_config = { - **self._binning(**treecorr_config), - "var_method": "jackknife" if int(npatch) > 1 else "shot", - } + reuse = npatch is None and out is None and not treecorr_config + npatch = int(npatch or self.npatch) + blind = self.blind(ver) if compute_tomography: + self._refuse_blinded_tomography(ver) tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver) if tomo_bin_ids is None or tomo_bin_pairs is None: raise ValueError(f"Version {ver} does not have tomography information.") @@ -77,24 +79,26 @@ def calculate_2pcf_version( tomo_bin_pairs = [("all", "all")] ggs = {f"tomo_bin_{b1}_tomo_bin_{b2}": None for b1, b2 in tomo_bin_pairs} - to_compute = [] - for bin1, bin2 in tomo_bin_pairs: - if (bin1, bin2) == ("all", "all"): - out_fname = self._xi_txt_path(ver, treecorr_config, npatch) - if os.path.exists(out_fname): - self.print_done(f"Skipping 2PCF calculation, {out_fname} exists") - gg = treecorr.GGCorrelation(treecorr_config) - gg.read(out_fname) - ggs["tomo_bin_all_tomo_bin_all"] = gg - continue - to_compute.append((bin1, bin2)) + held = self.xi_parts.get((ver, grid)) if reuse else None + if ALL_PAIR in ggs and held is not None: + ggs[ALL_PAIR] = sacc_io.xi_correlation(held) + to_compute = [ + (b1, b2) + for b1, b2 in tomo_bin_pairs + if ggs[f"tomo_bin_{b1}_tomo_bin_{b2}"] is None + ] if to_compute: + jackknife = npatch > 1 + gg_config = { + **self._binning(**treecorr_config), + "var_method": "jackknife" if jackknife else "shot", + } cols = self._shear_columns(ver, compute_tomography) patch_centers = self._patch_centers(cols, npatch) for bin1, bin2 in to_compute: - gg = treecorr.GGCorrelation(treecorr_config) + gg = treecorr.GGCorrelation(gg_config) cat_gal1 = self._bin_catalog(cols, bin1, npatch, patch_centers) cat_gal2 = ( @@ -105,29 +109,35 @@ def calculate_2pcf_version( gg.process(cat_gal1, cat2=cat_gal2) - if (bin1, bin2) == ("all", "all"): - # Columns only. The covariance matrix lives in the SACC part; - # a per-patch ξ± realisation is an unblinded data vector - # nothing reads; and TreeCorr cannot read back a text file - # carrying the matrix without the per-patch results. - gg.write( - self._xi_txt_path(ver, treecorr_config, npatch), - write_patch_results=False, - write_cov=False, - ) - - ggs[f"tomo_bin_{bin1}_tomo_bin_{bin2}"] = gg + if (bin1, bin2) != ("all", "all"): + ggs[f"tomo_bin_{bin1}_tomo_bin_{bin2}"] = gg + continue + s = xi_to_sacc( + self.sacc_nz(ver), + {**self.sacc_metadata(ver), "npatch": npatch}, + gg.meanr, + gg.xip, + gg.xim, + grid=grid, + theta_nom=gg.rnom, + npairs=gg.npairs, + weight=gg.weight, + covariance=gg.cov if jackknife else None, + variances=( + None if jackknife else np.concatenate([gg.varxip, gg.varxim]) + ), + ) + if out: + part = sacc_io.save(s, out, blind=blind) + else: + part = sacc_io.seal(s, blind) + self.xi_parts[ver, grid] = part + ggs[ALL_PAIR] = sacc_io.xi_correlation(part) self.print_done(f"Done 2PCF for {ver}.") return ggs - def _xi_txt_path(self, ver, treecorr_config, npatch): - """Path of the non-tomographic ξ± TreeCorr dump for a version.""" - return self._output_path( - f"xi_{self.basename(ver, treecorr_config=treecorr_config, npatch=npatch)}.txt" - ) - def _shear_columns(self, ver, compute_tomography): """Positions, calibrated shears, weights and bin labels of a version. @@ -204,7 +214,7 @@ def calculate_2pcf( Returns: dict: ``self.cat_ggs``, mapping each version to its - ``{"tomo_bin_{b1}_tomo_bin_{b2}": treecorr.GGCorrelation}`` dict. + ``{"tomo_bin_{b1}_tomo_bin_{b2}": ξ±}`` dict. """ self.cat_ggs = {} for ver in self.versions: @@ -223,56 +233,41 @@ def calculate_aperture_mass_dispersion( theta_max=200, nbins=500, nbins_map=15, - npatch=25, + npatch=None, compute_tomography=False, ): - self._map2 = {} - theta_map = np.geomspace(theta_min * 5, theta_max / 2, nbins_map) - self._map2["theta_map"] = theta_map + """⟨M_ap²⟩ and ⟨M_ײ⟩ of every version and bin pair, from its ξ±. - treecorr_config = self._binning(theta_min, theta_max, nbins) + Both are linear in ξ± (TreeCorr's ``calculateMapSq`` sum, Schneider et + al. 2002 filter), so they and their covariance T·C·Tᵀ come from the ξ± + of :meth:`calculate_2pcf_version` on a fine grid: for ``("all", + "all")`` its sealed part, concealed on a blinded catalogue, whose shift, + the same in every patch, leaves C unchanged. + """ + theta_map = np.geomspace(theta_min * 5, theta_max / 2, nbins_map) + self._map2 = {"theta_map": theta_map} + bin_size = np.log(theta_max / theta_min) / nbins for ver in self.versions: - if compute_tomography: - tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver) - if tomo_bin_ids is None or tomo_bin_pairs is None: - raise ValueError( - f"Version {ver} does not have tomography information." - ) - self.print_magenta( - f"Computing MAP for {ver} with {len(tomo_bin_pairs)} bins." - ) - else: - self.print_magenta(f"Computing non-tomographic MAP for {ver}.") - - tomo_bin_pairs = [("all", "all")] - - self._map2.setdefault(ver, {}) - cols = self._shear_columns(ver, compute_tomography) - patch_centers = self._patch_centers(cols, npatch) - - for bin1, bin2 in tomo_bin_pairs: - gg = treecorr.GGCorrelation(treecorr_config) - - cat_gal1 = self._bin_catalog(cols, bin1, npatch, patch_centers) - cat_gal2 = ( - self._bin_catalog(cols, bin2, npatch, patch_centers) - if bin1 != bin2 - else None - ) - - gg.process(cat_gal1, cat2=cat_gal2) - - mapsq, mapsq_im, mxsq, mxsq_im, varmapsq = gg.calculateMapSq( - R=theta_map, - m2_uform="Schneider", - ) - self._map2[ver][f"tomo_bin_{bin1}_tomo_bin_{bin2}"] = { + ggs = self.calculate_2pcf_version( + ver, + npatch=npatch, + compute_tomography=compute_tomography, + grid="aperture_mass", + min_sep=theta_min, + max_sep=theta_max, + nbins=nbins, + ) + self._map2[ver] = {} + for key, gg in ggs.items(): + transform = _map2_transform(theta_map, gg.meanr, bin_size) + mapsq, mxsq = np.split(transform @ np.concatenate([gg.xip, gg.xim]), 2) + variances = np.split(np.diag(transform @ gg.cov @ transform.T), 2) + self._map2[ver][key] = { "mapsq": mapsq, - "mapsq_im": mapsq_im, "mxsq": mxsq, - "mxsq_im": mxsq_im, - "varmapsq": varmapsq, + "varmapsq": variances[0], + "varmxsq": variances[1], } self.print_done(f"Done aperture-mass dispersion for {ver}.") @@ -712,13 +707,11 @@ def _mapsq_mxsq_sample_x_y_plot_function( y_plus = map2["mapsq"] * scale y_minus = map2["mxsq"] * scale - # Both E- and B-mode share the same variance estimate. - yerr = np.sqrt(map2["varmapsq"]) * scale ax_plus.errorbar( jittered_theta, y_plus, - yerr=yerr, + yerr=np.sqrt(map2["varmapsq"]) * scale, color=color, alpha=alpha, fmt="o", @@ -729,7 +722,7 @@ def _mapsq_mxsq_sample_x_y_plot_function( ax_minus.errorbar( jittered_theta, y_minus, - yerr=yerr, + yerr=np.sqrt(map2["varmxsq"]) * scale, color=color, alpha=alpha, fmt="o", @@ -781,3 +774,23 @@ def _ratio_xi_sys_xi_x_y_plot_function( ) if idx == 0: ax.axhspan(-threshold, threshold, color="black", alpha=0.1) + + +def _map2_transform(radii, theta, bin_size): + """The matrix taking [ξ+, ξ−] on a log grid to [⟨M_ap²⟩, ⟨M_ײ⟩] at ``radii``. + + TreeCorr's ``calculateMapSq`` sum with the Schneider et al. (2002) filter: + ⟨M_ap²⟩, ⟨M_ײ⟩ = Σ s² (T+ ξ+ ± T− ξ−) dlnθ / 2, s = θ/R, T± zero for s ≥ 2. + """ + s = np.minimum(np.outer(1.0 / radii, theta), 2.0) + ssq = s * s + tp = 12.0 / (5.0 * np.pi) * (2.0 - 15.0 * ssq) * np.arccos(s / 2.0) + tp += ( + s + * np.sqrt(4.0 - ssq) + * (120.0 + ssq * (2320.0 + ssq * (-754.0 + ssq * (132.0 - 9.0 * ssq)))) + / (100.0 * np.pi) + ) + tm = 3.0 / (70.0 * np.pi) * s * ssq * (4.0 - ssq) ** 3.5 + tp, tm = (x * ssq * 0.5 * bin_size for x in (tp, tm)) + return np.block([[tp, tm], [tp, -tm]]) diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index ff11db0b..2c2e7919 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -85,7 +85,9 @@ itself is not a dependency. """ +import functools import os +import types import numpy as np import sacc @@ -580,9 +582,69 @@ def get_xi(s, bins, *, grid): return _get_pm(s, XI_PLUS, XI_MINUS, _pair(bins), grid=grid) -def get_xi_weight(s, bins, *, grid): - """Return the TreeCorr pair weights stored with :func:`add_xi`'s ξ+ points.""" - return _tag(s, XI_PLUS, _pair(bins), "weight", grid=grid) +class _XiView(types.SimpleNamespace): + """What :func:`xi_correlation` returns; ``cov`` is sliced out when first read.""" + + @functools.cached_property + def cov(self): + return ( + None if self._covariance is None else self._covariance.get_block(self._rows) + ) + + +def xi_correlation(s, bins=(0, 0), grid=None): + """One ξ± series of ``s``, shaped like the TreeCorr ``GGCorrelation`` it came from. + + Carries ``meanr``, ``rnom``, ``xip``, ``xim``, ``varxip``, ``varxim``, + ``cov``, ``npairs``, ``weight``, ``left_edges``, ``right_edges`` and + ``npatch1``, so code written against a measurement reads a part unchanged. + ``grid`` may be left out when the pair's ξ± lies on one grid. The edges are + the log-binning edges about the nominal centres (``theta_nom``), and ``cov`` + is the series' block of the part's covariance. + """ + tracers = _pair(bins) + if grid is None: + grids = {s.data[i].tags.get("grid") for i in _indices(s, XI_PLUS, tracers)} + if len(grids) > 1: + raise ValueError( + f"ξ± of {tracers} lies on grids {sorted(grids)}; pass grid=" + ) + (grid,) = grids + tags = {} if grid is None else {"grid": grid} + plus, minus = (_indices(s, t, tracers, **tags) for t in (XI_PLUS, XI_MINUS)) + rows = np.concatenate([plus, minus]) + + def tag(name): + values = [s.data[i].tags.get(name) for i in plus] + return None if None in values else np.array(values, float) + + covariance = s.covariance + if covariance is None: + variances = np.full(len(rows), np.nan) + elif isinstance(covariance, sacc.covariance.DiagonalCovariance): + variances = np.asarray(covariance.diag)[rows] + else: + variances = np.diagonal(covariance.dense)[rows] + rnom = tag("theta_nom") + edges = (None, None) + if rnom is not None and len(rnom) > 1: + half_bin = np.log(rnom[-1] / rnom[0]) / (len(rnom) - 1) / 2 + edges = (rnom * np.exp(-half_bin), rnom * np.exp(half_bin)) + return _XiView( + meanr=tag("theta"), + rnom=rnom, + xip=s.mean[plus], + xim=s.mean[minus], + varxip=variances[: len(plus)], + varxim=variances[len(plus) :], + npairs=tag("npairs"), + weight=tag("weight"), + left_edges=edges[0], + right_edges=edges[1], + npatch1=int(s.metadata.get("npatch", 1)), + _covariance=covariance, + _rows=rows, + ) def get_pseudo_cl(s, bins): @@ -766,9 +828,8 @@ def merge(saccs): block-diagonal is out of scope here; see ``assemble_covariance``). Metadata must be consistent: keys present in several inputs must carry - equal values (a ``type: data`` file cannot merge with a ``type: mock`` - file), and the union lands on the result. This deliberately replaces the - library's clash behaviour, which mangles clashing keys by appending + equal values, and the union lands on the result. This deliberately replaces + the library's clash behaviour, which mangles clashing keys by appending labels. Grid consistency follows tagging semantics: the ``grid`` tag declares @@ -910,106 +971,76 @@ def _check_grid_consistency(s, angle): ) -def update_statistic(s, sub): - """Overwrite the values of ``s``'s points that match ``sub``'s, in place. +# --------------------------------------------------------------------------- # +# Saving: a birth is concealed under its catalogue's blind, a derivation +# inherits its inputs' stamp +# --------------------------------------------------------------------------- # +STAMP_KEY = "blind" - The merge-back half of the extract → conceal → merge blinding flow - (PRD #241 §4): each point of ``sub`` is matched to exactly one point of - ``s`` by ``(data_type, tracers, tags)``, and that point's *value* is - replaced. Nothing else changes — insertion order, tags, windows and the - covariance are untouched (blinding shifts the mean only), so ``sub``'s - own covariance (e.g. the sub-block ``extract`` attaches) is deliberately - not consulted. A ``sub`` point with no match, or with several, raises - ``ValueError``. - - Parameters - ---------- - s : sacc.Sacc - Target, mutated in place. - sub : sacc.Sacc - The replacement block, e.g. ``extract(s, ...)`` after concealment. - """ - claimed = set() - for point in sub.data: - idx = s.indices(point.data_type, point.tracers, **point.tags) - if len(idx) != 1: - raise ValueError( - f"update_statistic: {len(idx)} points in the target match " - f"({point.data_type}, {point.tracers}, {point.tags}) — need " - "exactly one" - ) - if idx[0] in claimed: - raise ValueError( - f"update_statistic: two sub points match the same target " - f"point ({point.data_type}, {point.tracers}, {point.tags})" - ) - claimed.add(idx[0]) - s.data[idx[0]].value = point.value +def stamp(s): + """The blind ``s`` was saved under: its ``blind`` metadata, ``none`` if absent.""" + return s.metadata.get(STAMP_KEY, "none") -def save(s, path, *, type): - """Write ``s`` to ``path`` (FITS), overwriting any existing file. - Parameters - ---------- - s : sacc.Sacc - Data set to write; its metadata is stamped in place. - type : {'data', 'mock'} - Provenance of the underlying catalogue, stored as the required - ``type`` metadata tag (PRD #241 §4, "Mocks vs data"). The caller — - the pipeline computing the data vector — knows whether its input - catalogue is a mock; there is deliberately no default. ``load`` - refuses ``type='data'`` files that are not blinded. - """ - if type not in ("data", "mock"): - raise ValueError(f"type must be 'data' or 'mock'; got {type!r}") - if s.metadata.get("type", type) != type: +def seal(s, blind, theory=None): + """What :func:`save` writes for a birth, kept in memory: a copy of ``s`` + concealed under ``blind`` with ``theory`` and stamped with its name.""" + if STAMP_KEY in s.metadata: raise ValueError( - f"Sacc metadata already carries type={s.metadata['type']!r}; " - f"refusing to re-stamp as {type!r}" + "this SACC already carries a blind stamp, so it is not a birth; " + "save it with derived_from=[its input parts]" ) - s.metadata["type"] = type - s.save_fits(path, overwrite=True) + if blind.name == "none": + out = s.copy() + else: + from .blinding import conceal + out = conceal(s, blind, theory) + out.metadata[STAMP_KEY] = blind.name + return out -def load(path, *, allow_unblinded=False): - """Load a Sacc from ``path`` (FITS), failing closed on unblinded data. - Every sacc_io file carries a ``type: data|mock`` metadata tag (stamped by - ``save``); blinded files are additionally stamped ``concealed=True`` by - Smokescreen. A ``type='data'`` file without that stamp is real, unblinded - data, and loading it raises — skipping the blind can never silently - expose the measured vector (PRD #241 §4). Mocks load freely, blinded or - not. +def save(s, path, *, blind=None, theory=None, derived_from=None): + """Write ``s`` to ``path`` (FITS) and return what was written. - Parameters - ---------- - path : str - File to load. - allow_unblinded : bool, optional - Escape hatch for the two legitimate consumers of unblinded data: - the blinding step itself (which must read the true vector to conceal - it) and the unblinding/verification tooling. Nothing else — no - analysis, plotting or inference code — may pass ``True``. + - ``save(s, path, blind=b)``: a birth. ``b`` is a + :class:`sp_validation.blinding.Blind`; unless it is ``none``, every row + of ``s`` is shifted by t(hidden) − t(fiducial), t being ``theory`` + (:mod:`sp_validation.blinding` describes a theory; the default is + :func:`sp_validation.blinding.shear`). The output is stamped with + ``b``'s name. + - ``save(s, path, derived_from=parts)``: a derivation (COSEBIs, pure-E/B, + an assembly); the parts must share one :func:`stamp`, which ``s`` takes, + unshifted. With ``blind=b`` too, that stamp must be ``b``'s name. - Returns - ------- - sacc.Sacc - The loaded data set. + A measurement computes its signal and saves (or seals) it in the same + function, returning the sealed part, so on a blinded catalogue its raw + signal never leaves that function. """ - s = sacc.Sacc.load_fits(path) - if ( - s.metadata["type"] == "data" - and not s.metadata.get("concealed", False) - and not allow_unblinded - ): - raise ValueError( - f"{path} holds real data (type='data') without the " - "concealed=True blinding stamp — refusing to load an unblinded " - "data vector. Only the blinding/unblinding tooling may pass " - "allow_unblinded=True." - ) - return s + if derived_from is not None: + stamps = {stamp(p) for p in derived_from} + if len(stamps) != 1: + raise ValueError(f"input parts carry different blinds: {sorted(stamps)}") + (name,) = stamps + if blind is not None and blind.name != name: + raise ValueError( + f"parts are stamped {name}, but the catalogue is declared under " + f"{blind.name}" + ) + out = s.copy() + out.metadata[STAMP_KEY] = name + elif blind is not None: + out = seal(s, blind, theory) + else: + raise ValueError("save needs blind= (a birth) or derived_from= (a derivation)") + out.save_fits(str(path), overwrite=True) + return out + + +def load(path): + """Load the SACC at ``path``.""" + return sacc.Sacc.load_fits(str(path)) # ============================================================================= diff --git a/src/sp_validation/tests/regression/test_tomo_drops_min_top_pin.py b/src/sp_validation/tests/regression/test_tomo_drops_min_top_pin.py index 09cea030..dc356c16 100644 --- a/src/sp_validation/tests/regression/test_tomo_drops_min_top_pin.py +++ b/src/sp_validation/tests/regression/test_tomo_drops_min_top_pin.py @@ -24,7 +24,11 @@ def _cosmology_validation(tmp_path): cfg = { "paths": {"output": str(tmp_path / "out")}, "nz": {"subdir": str(tmp_path), "dndz": {"path": "nz.txt"}}, - "synthetic": {"subdir": str(tmp_path), "shear": {"path": "cat.fits"}}, + "synthetic": { + "blind": "none", + "subdir": str(tmp_path), + "shear": {"path": "cat.fits"}, + }, } cfg_path = tmp_path / "cat_config.yaml" cfg_path.write_text(yaml.safe_dump(cfg)) diff --git a/src/sp_validation/tests/regression/test_tomo_plot_cosebis_ignores_tomography.py b/src/sp_validation/tests/regression/test_tomo_plot_cosebis_ignores_tomography.py index 63de87b7..645a0bb6 100644 --- a/src/sp_validation/tests/regression/test_tomo_plot_cosebis_ignores_tomography.py +++ b/src/sp_validation/tests/regression/test_tomo_plot_cosebis_ignores_tomography.py @@ -47,6 +47,7 @@ def _make_cv(tmp_path): out.mkdir() cfg = { VER: { + "blind": "none", "subdir": str(tmp_path), "shear": { "path": str(cat), diff --git a/src/sp_validation/tests/test_assemble_sacc.py b/src/sp_validation/tests/test_assemble_sacc.py index f17d5468..5c340e24 100644 --- a/src/sp_validation/tests/test_assemble_sacc.py +++ b/src/sp_validation/tests/test_assemble_sacc.py @@ -18,6 +18,7 @@ import pytest from sp_validation import sacc_io as sio +from sp_validation.blinding import NONE from sp_validation.cosmo_val import sacc_writers as sw @@ -157,7 +158,7 @@ def get_bandpower_windows(self): paths = {} for name, part in parts.items(): p = tmp_path / f"{name}.sacc" - sio.save(part, str(p), type="mock") + sio.save(part, str(p), blind=NONE) paths[name] = str(p) return paths @@ -170,7 +171,12 @@ def test_assemble_sacc_canonical_order(tmp_path): cl_cov_path, _blocks = _pseudo_cl_cov_fits(tmp_path) out = tmp_path / "vSYNTH.sacc" s = asm.assemble_sacc( - "vSYNTH", paths, str(out), xi_cov=cov_path, pseudo_cl_cov=cl_cov_path + "vSYNTH", + paths, + str(out), + blind="none", + xi_cov=cov_path, + pseudo_cl_cov=cl_cov_path, ) assert out.exists() assert type(s.covariance).__name__ == "BlockDiagonalCovariance" @@ -209,12 +215,17 @@ def test_injected_xi_covariance_replaces_the_parts_own(tmp_path): cov_path, xi_cov = _xi_cov_txt(tmp_path) cl_cov_path, _blocks = _pseudo_cl_cov_fits(tmp_path) - born = sio.load(paths["xi_reporting"], allow_unblinded=True).covariance.dense + born = sio.load(paths["xi_reporting"]).covariance.dense assert not np.allclose(born, xi_cov) # the two are distinguishable out = tmp_path / "vSYNTH.sacc" s = asm.assemble_sacc( - "vSYNTH", paths, str(out), xi_cov=cov_path, pseudo_cl_cov=cl_cov_path + "vSYNTH", + paths, + str(out), + blind="none", + xi_cov=cov_path, + pseudo_cl_cov=cl_cov_path, ) tr = ("source_0", "source_0") xi_idx = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)]) @@ -232,7 +243,7 @@ def test_assemble_sacc_injects_pseudo_cl_covariance(tmp_path): out = tmp_path / "vSYNTH.sacc" s = asm.assemble_sacc( - "vSYNTH", paths, str(out), xi_cov=cov_path, pseudo_cl_cov=cov_fits + "vSYNTH", paths, str(out), blind="none", xi_cov=cov_path, pseudo_cl_cov=cov_fits ) tr = ("source_0", "source_0") cl_idx = np.concatenate( @@ -253,7 +264,7 @@ def test_missing_injected_covariance_raises(tmp_path): paths = _write_parts(tmp_path, cov_less=()) # every part born with a block out = tmp_path / "vSYNTH.sacc" with pytest.raises(ValueError, match="takes its analysis covariance from"): - asm.assemble_sacc("vSYNTH", paths, str(out)) + asm.assemble_sacc("vSYNTH", paths, str(out), blind="none") def test_assemble_sacc_respects_pseudo_cl_toggle(tmp_path): @@ -262,7 +273,7 @@ def test_assemble_sacc_respects_pseudo_cl_toggle(tmp_path): assert "pseudo_cl" not in paths cov_path, _xi_cov = _xi_cov_txt(tmp_path) out = tmp_path / "vSYNTH.sacc" - s = asm.assemble_sacc("vSYNTH", paths, str(out), xi_cov=cov_path) + s = asm.assemble_sacc("vSYNTH", paths, str(out), blind="none", xi_cov=cov_path) tr = ("source_0", "source_0") assert len(s.indices(sio.CL_EE, tr)) == 0 # Round-trips as a valid BlockDiagonalCovariance over the remaining points. @@ -284,6 +295,7 @@ def test_assemble_sacc_expected_part_missing_raises(tmp_path): "vSYNTH", paths, str(out), + blind="none", expected=["xi_reporting", "pseudo_cl", "cosebis", "pure_eb", "rho_tau"], xi_cov=cov_path, ) @@ -296,5 +308,10 @@ def test_assemble_sacc_expected_rejects_unknown_name(tmp_path): out = tmp_path / "vSYNTH.sacc" with pytest.raises(ValueError, match="not assemblable statistics"): asm.assemble_sacc( - "vSYNTH", paths, str(out), expected=["cosebi"], xi_cov=cov_path + "vSYNTH", + paths, + str(out), + blind="none", + expected=["cosebi"], + xi_cov=cov_path, ) diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py new file mode 100644 index 00000000..d148e88b --- /dev/null +++ b/src/sp_validation/tests/test_blinding.py @@ -0,0 +1,104 @@ +"""What a blind shifts, and what the catalogue config may declare. + +The shifts run the default theory, :func:`sp_validation.blinding.shear`, under +a blind whose hidden point raises S8 above the fiducial. +""" + +import json + +import numpy as np +import pytest + +from sp_validation import b_modes +from sp_validation import blinding as bd +from sp_validation import sacc_io as sio + + +@pytest.fixture(scope="module") +def up(tmp_path_factory): + path = tmp_path_factory.mktemp("blinds") / "up.blind.json" + path.write_text( + json.dumps( + {"seed": "s", "envelope": {"S8": [0.05, 0.06]}, "fiducial": bd.fiducial()} + ) + ) + return bd.Blind("up", path) + + +def _sacc(): + z = np.linspace(0.0, 3.0, 200) + return sio.new_sacc({0: (z, np.exp(-0.5 * ((z - 0.7) / 0.2) ** 2))}) + + +@pytest.mark.slow +def test_a_blind_shifts_xi_and_cl_ee_and_leaves_bb_eb(up): + s = _sacc() + theta = np.geomspace(2.0, 200.0, 6) + sio.add_xi(s, (0, 0), theta, 0 * theta, 0 * theta, grid="reporting") + ell = np.array([30.0, 80.0, 150.0, 280.0, 450.0]) + w_ell = np.arange(2, 501).astype(float) + window = np.exp(-0.5 * ((w_ell[:, None] - ell) / 40.0) ** 2) + zeros = 0 * ell + sio.add_pseudo_cl( + s, + (0, 0), + ell, + zeros, + zeros, + zeros, + window_ells=w_ell, + window_weights=window / window.sum(axis=0), + ) + + sealed = sio.seal(s, up) + delta = sealed.mean - s.mean + for data_type in (sio.XI_PLUS, sio.XI_MINUS, sio.CL_EE): + assert np.all(delta[s.indices(data_type)] > 0), data_type + for data_type in (sio.CL_BB, sio.CL_EB): + assert np.all(delta[s.indices(data_type)] == 0), data_type + assert sio.stamp(sealed) == "up" + + +@pytest.mark.slow +def test_pure_eb_under_a_blind_shifts_e_and_not_b(up): + """The pure-E/B operator on the integration-grid ξ± is linear, so a blind + moves the modes by its response to the shift: E by several σ, B by under + a tenth of σ (σ: UNIONS shape noise through the operator).""" + edges_int = np.geomspace(0.08, 300.0, 1001) + theta_int = np.sqrt(edges_int[:-1] * edges_int[1:]) + s = _sacc() + sio.add_xi(s, (0, 0), theta_int, 0 * theta_int, 0 * theta_int, grid="integration") + shifted = sio.xi_correlation(sio.seal(s, up)) + delta = np.r_[shifted.xip, shifted.xim] + + annuli = np.diff(edges_int**2) # ∝ pair counts on a uniform field + operator, _, _ = b_modes.pure_eb_operator( + annuli, edges_int, np.geomspace(1.0, 250.0, 21) + ) + # Var ξ± = σ_e⁴ / (2 N_pairs), SP_v1.4.6.3's area, n_eff and σ_e. + n_eff, area, sigma_e = 4.96, 2894.0 * 3600, 0.378 + variance = np.tile(sigma_e**4 / (n_eff**2 * area * np.pi * annuli), 2) + sigma = np.sqrt(np.einsum("ij,j,ij->i", operator, variance, operator)) + response = operator @ delta / sigma + + n = len(response) // 6 # blocks in _EB_KEYS order: E+, E−, B+, B−, amb + assert np.max(np.abs(response[: 2 * n])) > 1 + assert np.max(np.abs(response[2 * n : 4 * n])) < 0.1 + + +def test_seal_refuses_a_sacc_already_stamped(): + sealed = sio.seal(_sacc(), bd.NONE) + with pytest.raises(ValueError, match="already carries a blind stamp"): + sio.seal(sealed, bd.NONE) + + +def test_entries_reading_one_file_declare_one_blind(): + catalogues = { + "paths": {"blinds": "/blinds"}, + "Y3": {"blind": "y3", "subdir": "/d", "shear": {"path": "y3.fits"}}, + "TWIN": {"blind": "none", "subdir": "/d", "shear": {"path": "twin.fits"}}, + } + assert bd.blind_of(catalogues, "TWIN") == "none" + catalogues["TWIN"]["shear"]["path"] = "/d/sub/../y3.fits" + with pytest.raises(ValueError, match="different blinds"): + bd.blind_of(catalogues, "TWIN") diff --git a/src/sp_validation/tests/test_column_schema.py b/src/sp_validation/tests/test_column_schema.py index d4ad4cfd..da073b80 100644 --- a/src/sp_validation/tests/test_column_schema.py +++ b/src/sp_validation/tests/test_column_schema.py @@ -64,6 +64,8 @@ # NaMaster field/workspace dicts keyed by bin; cosmology parameter dicts "W{}", "H0", + "A_IA", + "S8", # GLASS mock catalogues, which are not ShapePipe products "TOM_BIN_ID", } diff --git a/src/sp_validation/tests/test_config_paths_exist.py b/src/sp_validation/tests/test_config_paths_exist.py index 52646c9b..b5f593e9 100644 --- a/src/sp_validation/tests/test_config_paths_exist.py +++ b/src/sp_validation/tests/test_config_paths_exist.py @@ -177,12 +177,13 @@ def _located_elsewhere(source: Path, key: str, value: str) -> bool: """Whether a path-shaped value names nothing the config itself locates. ``paths.output`` in the catalogue config is where cosmo_val writes, created - by the run. A calibration ``params.input_path`` with no directory is opened + by the run, and ``paths.blinds`` the blind registry ``blinding init`` + creates. A calibration ``params.input_path`` with no directory is opened relative to its consumer's run directory (the image-sim run, or the catalogue directory ``scripts/masking.py`` prefixes). """ if source.name == "cat_config.yaml": - return key == "paths.output" + return key in ("paths.output", "paths.blinds") return ( source.parent.name == "calibration" and key == "params.input_path" diff --git a/src/sp_validation/tests/test_cosmo_val.py b/src/sp_validation/tests/test_cosmo_val.py index 7029cdcd..25cf0865 100644 --- a/src/sp_validation/tests/test_cosmo_val.py +++ b/src/sp_validation/tests/test_cosmo_val.py @@ -15,6 +15,7 @@ import pytest import yaml +from sp_validation import sacc_io from sp_validation.cosmo_val import CosmologyValidation # These tests load real UNIONS catalogues from the cluster filesystem. Skip them @@ -77,6 +78,7 @@ def _make_seed_config(tmp_path, shear_filename): }, "paths": {"output": str(output_dir)}, base_version: { + "blind": "none", "subdir": str(base_dir), "pipeline": "SP", "shear": { @@ -193,6 +195,24 @@ def test_seed_variant_updates_shear_path(self, tmp_path): assert seed_version in cv.cc assert cv.cc[seed_version]["shear"]["path"].endswith("shear_seed_007.fits") + def test_blind_is_checked_against_the_config_as_written(self, tmp_path): + """A twin entry reading the same shear file under another blind is + refused, though the selected entry's paths were resolved; a seeded + leak-corrected variant takes its base entry's blind.""" + params, base_version = self._make_seed_config( + tmp_path, shear_filename="shear_seed_1234.fits" + ) + config = yaml.safe_load(open(params["catalog_config"])) + cv = CosmologyValidation(versions=[f"{base_version}_seed7_leak_corr"], **params) + assert cv.blind(f"{base_version}_seed7_leak_corr").name == "none" + + config["Twin"] = {**config[base_version], "blind": "y3"} + with open(params["catalog_config"], "w") as f: + yaml.dump(config, f, sort_keys=False) + cv = CosmologyValidation(versions=[base_version], **params) + with pytest.raises(ValueError, match="different blinds"): + cv.blind(base_version) + def test_seed_leak_corr_materializes_seed_first(self, tmp_path): """_seed_leak_corr should clone the seed variant before leak fixes.""" params, base_version = self._make_seed_config( @@ -376,6 +396,7 @@ def _write_synthetic_catalogs( "e2_col": "HSM_G2_PSF", } version_cfg = { + "blind": "none", "subdir": str(cat_dir), "pipeline": "SP", "colour": "tab:blue", @@ -416,13 +437,15 @@ def _write_synthetic_catalogs( return params, version def test_calculate_2pcf_runs_on_synthetic_catalog(self, tmp_path): - """calculate_2pcf wires catalog+config into treecorr GGCorrelation. + """calculate_2pcf_version wires catalog+config into a ξ± part. - Smoke-integration: assert the xi+/- data vector is computed with the - configured number of angular bins and is finite. Numerical values are - deliberately not asserted (could be tightened to allclose vs. a - committed reference later for value-drift coverage). + Smoke-integration: the part's ξ± has the configured number of angular + bins and is finite, and the part reads back with the binning's edges, + by which scale cuts select bins. Numerical values are deliberately not + asserted. """ + from sp_validation.b_modes import log_bin_edges + pytest.importorskip("treecorr") params, version = self._write_synthetic_catalogs(tmp_path) @@ -436,14 +459,16 @@ def test_calculate_2pcf_runs_on_synthetic_catalog(self, tmp_path): **params, ) - ggs = cv.calculate_2pcf_version(version) + cv.calculate_2pcf_version(version) + gg = sacc_io.xi_correlation(cv.xi_parts[version, "reporting"]) - # treecorr GGCorrelation with xi+/- on the configured angular grid - # LG:Hardcoded to be non-tomographic for now - assert ggs["tomo_bin_all_tomo_bin_all"].xip.shape == (nbins,) - assert ggs["tomo_bin_all_tomo_bin_all"].xim.shape == (nbins,) - assert np.all(np.isfinite(ggs["tomo_bin_all_tomo_bin_all"].xip)) - assert np.all(np.isfinite(ggs["tomo_bin_all_tomo_bin_all"].xim)) + edges = log_bin_edges(5.0, 100.0, nbins) + np.testing.assert_allclose(gg.left_edges, edges[0], rtol=1e-12) + np.testing.assert_allclose(gg.right_edges, edges[1], rtol=1e-12) + assert gg.xip.shape == (nbins,) + assert gg.xim.shape == (nbins,) + assert np.all(np.isfinite(gg.xip)) + assert np.all(np.isfinite(gg.xim)) # The additive-bias subtraction in the pipeline must have run. assert version in cv.c1 and version in cv.c2 @@ -460,8 +485,6 @@ def test_xi_part_carries_the_covariance_the_measurement_estimated( import sacc - from sp_validation import sacc_io - script = Path(__file__).resolve().parents[3] / "workflow/scripts/run_2pcf.py" spec = importlib.util.spec_from_file_location("run_2pcf_part", script) run_2pcf = importlib.util.module_from_spec(spec) @@ -469,7 +492,7 @@ def test_xi_part_carries_the_covariance_the_measurement_estimated( params, version = self._write_synthetic_catalogs(tmp_path) part = tmp_path / "part.sacc" - gg = run_2pcf.run_2pcf( + written = run_2pcf.run_2pcf( ver=version, min_sep=5.0, max_sep=100.0, @@ -480,37 +503,33 @@ def test_xi_part_carries_the_covariance_the_measurement_estimated( sacc_out=str(part), ) - cov = sacc_io.load(str(part), allow_unblinded=True).covariance - if npatch > 1: - assert isinstance(cov, sacc.covariance.FullCovariance) - np.testing.assert_array_equal(cov.dense, gg.cov) - else: - assert isinstance(cov, sacc.covariance.DiagonalCovariance) - np.testing.assert_array_equal( - cov.diag, np.concatenate([gg.varxip, gg.varxim]) - ) - - def test_a_patched_xi_dump_reads_back(self, tmp_path): - """calculate_2pcf_version reads back the dump a patched measurement wrote. + cov = sacc_io.load(str(part)).covariance + assert isinstance( + cov, + sacc.covariance.FullCovariance + if npatch > 1 + else sacc.covariance.DiagonalCovariance, + ) + np.testing.assert_array_equal(cov.dense, written.covariance.dense) + assert np.all(np.diag(cov.dense) > 0) - The ξ± figure rules re-enter calculate_2pcf on the reporting grid, which - has patches, and are handed the dump rule xi wrote (to its precision). - """ + def test_a_held_part_is_drawn_without_measuring(self, tmp_path, monkeypatch): + """The ξ± figure rules hand calculate_2pcf_version the reporting parts + rule xi wrote; it draws them and never reads the catalogue.""" params, version = self._write_synthetic_catalogs(tmp_path) - binning = dict(npatch=4, min_sep=5.0, max_sep=100.0, nbins=6) - pair = "tomo_bin_all_tomo_bin_all" - measured = CosmologyValidation( - versions=[version], **params - ).calculate_2pcf_version(version, **binning)[pair] - dumps = list(Path(params["output_dir"]).glob(f"xi_{version}_*.txt")) - assert len(dumps) == 1 - read = CosmologyValidation(versions=[version], **params).calculate_2pcf_version( - version, **binning - )[pair] - for column in ("meanr", "npairs", "xip", "xim", "varxip", "varxim"): - np.testing.assert_allclose( - getattr(read, column), getattr(measured, column), rtol=1e-4 - ) + measuring = CosmologyValidation(versions=[version], npatch=4, **params) + measuring.calculate_2pcf_version(version) + held = measuring.xi_parts[version, "reporting"] + + cv = CosmologyValidation(versions=[version], npatch=4, **params) + cv.xi_parts[version, "reporting"] = held + + def no_catalogue(*args, **kwargs): + raise AssertionError("measured instead of drawing the held part") + + monkeypatch.setattr(cv, "_shear_columns", no_catalogue) + drawn = cv.calculate_2pcf()[version]["tomo_bin_all_tomo_bin_all"] + np.testing.assert_array_equal(drawn.xip, sacc_io.xi_correlation(held).xip) def test_calculate_2pcf_is_reproducible_across_machines( self, tmp_path, monkeypatch @@ -519,11 +538,12 @@ def test_calculate_2pcf_is_reproducible_across_machines( import treecorr import treecorr.field - patches = [] + patches, measured = [], [] process = treecorr.GGCorrelation.process def recording_process(gg, cat, *args, **kwargs): patches.append(np.array(cat.patch)) + measured.append(gg) return process(gg, cat, *args, **kwargs) monkeypatch.setattr(treecorr.GGCorrelation, "process", recording_process) @@ -552,7 +572,9 @@ def recording_process(gg, cat, *args, **kwargs): ] xi[tree] = np.concatenate([gg.xip, gg.xim]) var[tree] = np.concatenate([gg.varxip, gg.varxim]) - counts[tree] = {key: result.npairs for key, result in gg.results.items()} + counts[tree] = { + key: result.npairs for key, result in measured[-1].results.items() + } np.testing.assert_array_equal(patches[0], patches[1]) assert counts["a"].keys() == counts["b"].keys() @@ -607,6 +629,27 @@ def recording_bin_catalog(cols, bin_id, npatch, patch_centers=None): cross_pair_catalogs[1][2], expected_centers, rtol=0, atol=1e-14 ) + def test_a_blinded_catalogue_refuses_tomographic_signal(self, tmp_path): + """Blinding shifts the ("all", "all") parts only, so a blinded + catalogue's tomographic ξ± and pseudo-Cℓ are refused before measuring.""" + from sp_validation import blinding + + params, version = self._write_synthetic_catalogs( + tmp_path, n_gal=400, with_tomography=True + ) + config = yaml.safe_load(open(params["catalog_config"])) + config["paths"]["blinds"] = str(tmp_path / "blinds") + config[version]["blind"] = "toy" + with open(params["catalog_config"], "w") as f: + yaml.dump(config, f, sort_keys=False) + blinding.init("toy", config) + cv = CosmologyValidation(versions=[version], npatch=1, **params) + + with pytest.raises(blinding.BlindingError, match="tomographic"): + cv.calculate_2pcf_version(version, compute_tomography=True) + with pytest.raises(blinding.BlindingError, match="tomographic"): + cv.calculate_pseudo_cl(compute_tomography=True) + def test_treecorr_runs_on_the_cpus_the_process_holds(self, tmp_path): """By default TreeCorr takes the process's CPU affinity, not the node's count.""" import treecorr @@ -657,17 +700,14 @@ def test_calculate_pure_eb_runs_on_synthetic_catalog(self, tmp_path, pure_eb_xi) The ξ± it measures equal the committed ``pure_eb_xi`` (``test_b_modes`` pins the operator on the same ξ±, so a failure names the step that - moved), and the jackknife covariance of the modes is the jackknife ξ± - covariance through the operator: for a linear estimator that equals - TreeCorr's per-patch jackknife of the modes themselves. + moved), and the covariance of the modes is the jackknife ξ± covariance + of its integration-grid part through the operator. ξ±: exact binning (bin_slop = angle_slop = 0) makes ξ± a plain pair sum, - independent of the tree and so of the jackknife patches, whose k-means - centres this test does not fix. + independent of the tree and so of the jackknife patches. """ pytest.importorskip("treecorr") pytest.importorskip("cosmo_numba") - import treecorr from sp_validation import b_modes @@ -720,25 +760,9 @@ def test_calculate_pure_eb_runs_on_synthetic_catalog(self, tmp_path, pure_eb_xi) cov = np.asarray(results["cov"]) assert cov.shape == (6 * nbins, 6 * nbins) assert results["npatch"] == npatch - # The integration-grid ξ± again, measured on the jackknife patches the - # first measurement wrote (its columns-only dump carries no patches). - for dump in Path(params["output_dir"]).glob(f"xi_{version}_*.txt"): - dump.unlink() - gg_int = cv.calculate_2pcf_version( - version, - npatch=npatch, - min_sep=integration["min_sep_int"], - max_sep=integration["max_sep_int"], - nbins=integration["nbins_int"], - )["tomo_bin_all_tomo_bin_all"] - jackknife_of_modes = treecorr.estimate_multi_cov( - [gg_int], - "jackknife", - func=lambda corrs: operator @ np.concatenate([corrs[0].xip, corrs[0].xim]), - cross_patch_weight="match", - ) + cov_xi = sacc_io.xi_correlation(cv.xi_parts[version, "integration"]).cov np.testing.assert_allclose( - cov, jackknife_of_modes, rtol=0, atol=1e-10 * np.abs(cov).max() + cov, operator @ cov_xi @ operator.T, rtol=0, atol=1e-10 * np.abs(cov).max() ) def test_plot_2pcf_writes_the_non_tomographic_figures(self, tmp_path): @@ -867,3 +891,31 @@ def test_pol_factor_must_be_plus_or_minus_one(self, tmp_path, pol_factor): params, version = self._write_synthetic_catalogs(tmp_path) with pytest.raises(ValueError, match="pol_factor"): CosmologyValidation(versions=[version], pol_factor=pol_factor, **params) + + +def test_map2_transform_is_treecorrs_calculate_map_sq(): + """⟨M_ap²⟩, ⟨M_ײ⟩ from the ξ± transform equal TreeCorr's own sum.""" + import treecorr + + from sp_validation.cosmo_val.real_space import _map2_transform + + rng = np.random.default_rng(1) + n = 5000 + cat = treecorr.Catalog( + ra=rng.uniform(0, 5, n), + dec=rng.uniform(0, 5, n), + g1=rng.normal(0, 0.3, n), + g2=rng.normal(0, 0.3, n), + ra_units="deg", + dec_units="deg", + ) + gg = treecorr.GGCorrelation( + min_sep=0.5, max_sep=200, nbins=200, sep_units="arcmin", bin_slop=0 + ) + gg.process(cat) + radii = np.geomspace(1, 60, 8) + mapsq, _, mxsq, _, _ = gg.calculateMapSq(R=radii, m2_uform="Schneider") + ours = _map2_transform(radii, gg.meanr, gg.bin_size) @ np.r_[gg.xip, gg.xim] + treecorrs = np.r_[mapsq, mxsq] + atol = 1e-12 * np.max(np.abs(treecorrs)) + np.testing.assert_allclose(ours, treecorrs, rtol=1e-12, atol=atol) diff --git a/src/sp_validation/tests/test_pseudo_cl.py b/src/sp_validation/tests/test_pseudo_cl.py index e3dae689..5e4724f0 100644 --- a/src/sp_validation/tests/test_pseudo_cl.py +++ b/src/sp_validation/tests/test_pseudo_cl.py @@ -160,6 +160,7 @@ def _write_synthetic_config(tmp_path): "star": {**psf_cfg}, "psf": psf_cfg, "patch_number": 150, + "blind": "none", }, } config_path = tmp_path / "config.yaml" @@ -694,10 +695,10 @@ def capture_workspace( cv.calculate_pseudo_cl_catalog(ver, out_path, tomo_bin_a="all", tomo_bin_b="all") assert os.path.exists(out_path) - s = sacc_io.load(out_path, allow_unblinded=True) + s = sacc_io.load(out_path) ell, ee, bb, eb, window = sacc_io.get_pseudo_cl(s, SACC_BIN) assert window is not None # the shared BandpowerWindow rides the part - readback = cv._load_pseudo_cl(out_path, ("all", "all")) + readback = cv._load_pseudo_cl(ver, out_path, ("all", "all")) npt.assert_array_equal(readback["BE"], eb) assert not np.shares_memory(readback["BE"], readback["EB"]) assert saved["nside"] is None @@ -796,7 +797,7 @@ def capture_workspace( out_path = cv._output_path(f"pseudo_cl_map_{ver}.sacc") cv.calculate_pseudo_cl_map(ver, NSIDE, out_path, "all", "all") - s = sacc_io.load(out_path, allow_unblinded=True) + s = sacc_io.load(out_path) _ell, _ee, _bb, _eb, window = sacc_io.get_pseudo_cl(s, SACC_BIN) assert saved["nside"] == NSIDE window_ells, unwindowed_weights = bandpower_window_from_workspace( diff --git a/src/sp_validation/tests/test_sacc_io.py b/src/sp_validation/tests/test_sacc_io.py index 872a8198..d2279969 100644 --- a/src/sp_validation/tests/test_sacc_io.py +++ b/src/sp_validation/tests/test_sacc_io.py @@ -9,6 +9,7 @@ import pytest from sp_validation import sacc_io as sio +from sp_validation.blinding import NONE # --------------------------------------------------------------------------- # @@ -64,7 +65,7 @@ def _cosebi_block(s, tr): # --------------------------------------------------------------------------- # def _roundtrip(s, tmp_path, name="rt"): path = tmp_path / f"{name}.sacc" - sio.save(s, str(path), type="mock") + sio.save(s, str(path), blind=NONE) return sio.load(str(path)) @@ -100,7 +101,7 @@ def test_xi_roundtrip(tmp_path): assert np.array_equal(th, theta) assert np.array_equal(p, xip) assert np.array_equal(m, xim) - assert np.array_equal(sio.get_xi_weight(s2, (0, 0), grid="reporting"), weight) + assert np.array_equal(sio.xi_correlation(s2, grid="reporting").weight, weight) # extra tags survive idx = s2.indices(sio.XI_PLUS, ("source_0", "source_0"), grid="reporting") tags = s2.data[idx[0]].tags @@ -466,7 +467,7 @@ def test_end_to_end_one_file_layout(tmp_path): s, [(xi_c, _spd(len(xi_c), 1)), (co, _spd(len(co), 2)), (xi_f, fine_block)], ) - sio.save(s, str(tmp_path / f"{version}.sacc"), type="mock") + sio.save(s, str(tmp_path / f"{version}.sacc"), blind=NONE) a = sio.load(str(tmp_path / f"{version}.sacc")) @@ -505,7 +506,7 @@ def test_one_file_layout_diagonal_fine_fallback(tmp_path): xi_f = _xi_block(s, tr, grid="integration") 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"), type="mock") + sio.save(s, str(tmp_path / "vDIAG.sacc"), blind=NONE) a = sio.load(str(tmp_path / "vDIAG.sacc")) assert np.array_equal(np.diag(a.covariance.dense[np.ix_(xi_f, xi_f)]), variances) @@ -613,80 +614,10 @@ def test_tomographic_xi_covariance_one_contiguous_block(): assert np.array_equal(s.mean[idx_p], xip) -# --------------------------------------------------------------------------- # -# 12. type stamping + fail-closed load (data/mock x concealed/not; escape -# hatch for the blinding/unblinding tooling) -# --------------------------------------------------------------------------- # -def _saved(tmp_path, name, *, type, concealed=None): - s = _base_sacc() - _add_xi(s) - if concealed is not None: - s.metadata["concealed"] = concealed - path = str(tmp_path / f"{name}.sacc") - sio.save(s, path, type=type) - return path - - -def test_load_mock_unconcealed(tmp_path): - s = sio.load(_saved(tmp_path, "m0", type="mock")) - assert s.metadata["type"] == "mock" - - -def test_load_mock_concealed(tmp_path): - s = sio.load(_saved(tmp_path, "m1", type="mock", concealed=True)) - assert s.metadata["concealed"] - - -def test_load_data_concealed(tmp_path): - s = sio.load(_saved(tmp_path, "d1", type="data", concealed=True)) - assert s.metadata["type"] == "data" - - -def test_load_data_unconcealed_fails_closed(tmp_path): - path = _saved(tmp_path, "d0", type="data") - with pytest.raises(ValueError, match="unblinded"): - sio.load(path) - # concealed=False is as unblinded as no stamp at all - path_f = _saved(tmp_path, "d0f", type="data", concealed=False) - with pytest.raises(ValueError, match="unblinded"): - sio.load(path_f) - - -def test_load_data_unconcealed_escape_hatch(tmp_path): - s = sio.load(_saved(tmp_path, "d0h", type="data"), allow_unblinded=True) - assert s.metadata["type"] == "data" - - -def test_save_requires_valid_type(tmp_path): - s = _base_sacc() - with pytest.raises(TypeError): - sio.save(s, str(tmp_path / "x.sacc")) # type is required - with pytest.raises(ValueError, match="'data' or 'mock'"): - sio.save(s, str(tmp_path / "x.sacc"), type="simulation") - - -def test_save_refuses_type_restamp(tmp_path): - s = _base_sacc() - sio.save(s, str(tmp_path / "x.sacc"), type="mock") - with pytest.raises(ValueError, match="re-stamp"): - sio.save(s, str(tmp_path / "x.sacc"), type="data") - - -def test_load_requires_type_tag(tmp_path): - import sacc as sacc_lib - - s = _base_sacc() # never stamped - path = str(tmp_path / "untyped.sacc") - s.save_fits(path, overwrite=True) - with pytest.raises(KeyError): - sio.load(path) - assert sacc_lib.Sacc.load_fits(path) is not None # raw loader still works - - # --------------------------------------------------------------------------- # # 13. merge(): per-statistic files combine into one; shared tracers stored # once; covariance block-diagonal (all-or-none); metadata union with -# loud conflicts. update_statistic(): value-only merge-back. +# loud conflicts. # --------------------------------------------------------------------------- # def _xi_sacc(metadata=None): s = sio.new_sacc({0: _nz(0)}, metadata=metadata) @@ -703,18 +634,18 @@ def _cosebi_sacc(metadata=None): def test_merge_per_statistic_files(tmp_path): - meta = {"version": "vM", "type": "mock"} + meta = {"version": "vM", "survey": "UNIONS"} s_xi, s_co = _xi_sacc(meta), _cosebi_sacc(meta) merged = sio.merge([s_xi, s_co]) # shared tracers stored once; all points present, xi first assert set(merged.tracers) == {"source_0", sio.PSF_TRACER} assert len(merged.mean) == len(s_xi.mean) + len(s_co.mean) assert np.array_equal(merged.mean, np.concatenate([s_xi.mean, s_co.mean])) - assert merged.metadata["version"] == "vM" and merged.metadata["type"] == "mock" + assert merged.metadata["version"] == "vM" and merged.metadata["survey"] == "UNIONS" # inputs untouched assert s_xi.metadata["version"] == "vM" # readers work on the merged file after a round-trip - sio.save(merged, str(tmp_path / "vM.sacc"), type="mock") + sio.save(merged, str(tmp_path / "vM.sacc"), blind=NONE) merged_rt = sio.load(str(tmp_path / "vM.sacc")) _, p, _ = sio.get_xi(merged_rt, (0, 0), grid="reporting") assert np.array_equal(p, np.arange(6) * 1e-5) @@ -749,7 +680,7 @@ def test_merge_block_diagonal_covariance_stays_block_diagonal(tmp_path): assert type(s_xi.covariance).__name__ == "BlockDiagonalCovariance" merged = sio.merge([s_xi, s_co]) assert type(merged.covariance).__name__ == "BlockDiagonalCovariance" - sio.save(merged, str(tmp_path / "vBLK.sacc"), type="mock") + sio.save(merged, str(tmp_path / "vBLK.sacc"), blind=NONE) merged_rt = sio.load(str(tmp_path / "vBLK.sacc")) assert type(merged_rt.covariance).__name__ == "BlockDiagonalCovariance" n_xi = len(s_xi.mean) @@ -769,52 +700,12 @@ def test_merge_mixed_covariance_fails(): def test_merge_conflicting_metadata_fails(): - s_xi = _xi_sacc({"type": "data"}) - s_co = _cosebi_sacc({"type": "mock"}) + s_xi = _xi_sacc({"survey": "UNIONS"}) + s_co = _cosebi_sacc({"survey": "KiDS"}) with pytest.raises(ValueError, match="conflicting metadata"): sio.merge([s_xi, s_co]) -def test_update_statistic_values_only(): - s = _multi_statistic_sacc() - tr = ("source_0", "source_0") - xi, cl, co = _xi_block(s, tr), _cl_block(s, tr), _cosebi_block(s, tr) - cov = [(xi, _spd(len(xi), 1)), (cl, _spd(len(cl), 2)), (co, _spd(len(co), 3))] - sio.assemble_covariance(s, cov) - dense_before = s.covariance.dense.copy() - mean_before = s.mean.copy() - - # extract -> shift (a stand-in for conceal) -> merge back - sub = sio.extract(s, data_type=sio.XI_PLUS, tracers=tr) - for point in sub.data: - point.value += 1e-4 - sio.update_statistic(s, sub) - - idx_p = s.indices(sio.XI_PLUS, tr) - assert np.allclose(s.mean[idx_p], mean_before[idx_p] + 1e-4) - untouched = np.setdiff1d(np.arange(len(s.mean)), idx_p) - assert np.array_equal(s.mean[untouched], mean_before[untouched]) - # covariance and insertion order untouched - assert np.array_equal(s.covariance.dense, dense_before) - - -def test_update_statistic_requires_unique_match(): - s = _xi_sacc() - sub = sio.extract(s, data_type=sio.XI_PLUS, tracers=("source_0", "source_0")) - missing = sub.copy() - missing.data[0].tags["theta"] = 999.0 # matches nothing in s - with pytest.raises(ValueError, match="0 points"): - sio.update_statistic(s, missing) - - -def test_update_statistic_rejects_duplicate_sub_points(): - s = _xi_sacc() - sub = sio.extract(s, data_type=sio.XI_PLUS, tracers=("source_0", "source_0")) - sub.data[1].tags = dict(sub.data[0].tags) # two sub points -> one target - with pytest.raises(ValueError, match="same target"): - sio.update_statistic(s, sub) - - def test_merge_rejects_divergent_shared_tracer(): s_xi, s_co = _xi_sacc(), _cosebi_sacc() s_co.tracers["source_0"].nz = s_co.tracers["source_0"].nz * 2.0 diff --git a/src/sp_validation/tests/test_sacc_writers.py b/src/sp_validation/tests/test_sacc_writers.py index 74181190..a3eb6618 100644 --- a/src/sp_validation/tests/test_sacc_writers.py +++ b/src/sp_validation/tests/test_sacc_writers.py @@ -15,6 +15,7 @@ import pytest from sp_validation import sacc_io as sio +from sp_validation.blinding import NONE from sp_validation.cosmo_val import sacc_writers as sw @@ -34,7 +35,7 @@ def _theta(n=6): def _roundtrip(s, tmp_path, name): p = tmp_path / f"{name}.sacc" - sio.save(s, str(p), type="mock") + sio.save(s, str(p), blind=NONE) return sio.load(str(p)) @@ -143,7 +144,7 @@ def get_bandpower_windows(self): part = sw.pseudo_cl_to_sacc({0: _nz()}, META, ell, cl_all, _Workspace()) path = tmp_path / "pseudo_cl.sacc" - sio.save(part, str(path), type="mock") + sio.save(part, str(path), blind=NONE) data = _paper_pseudo_cl_reader()(path) assert np.array_equal(data["ELL"], ell) @@ -361,7 +362,7 @@ def test_assemble_from_reloaded_parts(tmp_path): parts = _make_parts(nz) reloaded = [] for i, part in enumerate(parts): - sio.save(part, str(tmp_path / f"part{i}.sacc"), type="mock") + sio.save(part, str(tmp_path / f"part{i}.sacc"), blind=NONE) reloaded.append(sio.load(str(tmp_path / f"part{i}.sacc"))) s = sw.assemble_analysis_sacc(reloaded) assert type(s.covariance).__name__ == "BlockDiagonalCovariance" diff --git a/uv.lock b/uv.lock index 7d11c489..cc421f6f 100644 --- a/uv.lock +++ b/uv.lock @@ -3885,7 +3885,6 @@ dependencies = [ { name = "clmm" }, { name = "colorama" }, { name = "cosmo-numba" }, - { name = "cryptography" }, { name = "cs-util" }, { name = "emcee" }, { name = "getdist" }, @@ -3968,7 +3967,6 @@ requires-dist = [ { name = "cosmo-numba", git = "https://github.com/cailmdaley/cosmo-numba.git?rev=a64cb2ed13e595a16ad056a128ef24594b811504" }, { name = "cosmology-compat-camb", marker = "extra == 'glass'", specifier = "==0.2.0" }, { name = "cosmosis", marker = "extra == 'workflow'", specifier = ">=3.25" }, - { name = "cryptography" }, { name = "cs-util", git = "https://github.com/CosmoStat/cs_util.git?rev=develop" }, { name = "emcee" }, { name = "fast-pt", marker = "extra == 'workflow'", specifier = ">=3.2,<4" }, diff --git a/workflow/README.md b/workflow/README.md index eda38c78..8af65cab 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -25,6 +25,35 @@ Runs stay modular, not monolithic: a paper or run composes these rules with Snakemake's `module` directive under its own config and an output `prefix`, so each namespaces cleanly under `results//`. +## Blinding + +Every entry of `cosmo_val/cat_config.yaml` declares `blind: none` or +`blind: `; an entry without one stops the launch. `_leak_corr` versions +take their entry's blind, and entries reading one shear file declare one blind. + +A blind is a secret record, `/.blind.json` (seed, envelope, +fiducial), drawn once with +`spv-container exec python -m sp_validation.blinding init ` and never +committed. On a blinded catalogue, ξ± and Cℓ_EE are shifted by +t(hidden) − t(fiducial) before they are first written: Smokescreen draws the +hidden point (S8 and Ωm) from the seed, and the theory `t` is +`sp_validation.blinding.shear`. Cℓ_BB and Cℓ_EB get zero shift; ρ/τ are saved +with `blinding.no_signal`. Another statistic passes its own theory, +`theory(params, s) -> array the length of s.mean` (see +`src/sp_validation/blinding.py`). + +A measurement calculates, then saves: one function computes the signal, saves +it with `sacc_io.save(s, path, blind=...)` (or `sacc_io.seal`), and returns the +concealed part, so raw signal never leaves it +(`CosmologyValidation.calculate_2pcf_version` is the pattern). Statistics computed from +parts (COSEBIs, pure-E/B, the assembled `{version}.sacc`) are saved with +`derived_from=parts` and carry their inputs' stamp, `s.metadata["blind"]`. + +Producer rules carry the catalogue's blind in `params.blind`, so unblinding is +flipping the declaration to `blind: none` and rerunning: Snakemake reruns +exactly the parts it touches, and assembly refuses a part whose stamp is not the +catalogue's declared blind. + ## Running on the cluster — the candide profile `profiles/candide/config.yaml` is the committed SLURM profile: it hands diff --git a/workflow/common.py b/workflow/common.py index 111da5f5..9290a83e 100644 --- a/workflow/common.py +++ b/workflow/common.py @@ -33,21 +33,33 @@ def _plain(path): REPO_ROOT = _plain(__file__).parents[1] REPO_SRC = REPO_ROOT / "src" -# The container model lives in the package (``sp_validation/container.py``). -# Taken from *this checkout's* src/, so the workflow and the ``spv-container`` -# CLI can never disagree about which image to run. -_container = importlib.util.module_from_spec( - importlib.util.spec_from_file_location( - "_spv_container", REPO_SRC / "sp_validation" / "container.py" + +def _load_checkout_module(name): + """Load ``sp_validation/.py``, a stdlib-only module, from this checkout. + + By path, since the host Snakemake has no sp_validation installed; from + *this checkout's* src/, so the workflow and the package never disagree. + """ + alias = f"_spv_{name}" + module = importlib.util.module_from_spec( + importlib.util.spec_from_file_location( + alias, REPO_SRC / "sp_validation" / f"{name}.py" + ) ) -) -sys.modules["_spv_container"] = _container -_container.__loader__.exec_module(_container) + sys.modules[alias] = module + module.__loader__.exec_module(module) + return module + +# The container model, shared with the ``spv-container`` CLI. +_container = _load_checkout_module("container") compare_revision = _container.compare_revision image_revision = _container.image_revision resolve_image = _container.resolve_image +# Each catalogue's declared blind (``blinding.blind_of``). +_blinding = _load_checkout_module("blinding") + # Every job inherits this launch's environment (the slurm executor submits with # --export=ALL), and the Snakemake each job step starts keeps its source cache # under XDG_CACHE_HOME, which a login shell may point at node-local storage. @@ -193,6 +205,16 @@ def configure(workflow_config): ) with open(COSMOLOGY_PARAMS) as f: PLANCK18 = json.load(f) + # A run catalogue that declares no blind stops the launch before any job. + fiducial = (FIDUCIAL.get(k) for k in ("version", "mock_version")) + for version in [*workflow_config.get("versions", []), *filter(None, fiducial)]: + blind_of(version) + + +def blind_of(version): + """``version``'s declared blind: a producer's ``params.blind``, so flipping + it reruns exactly the jobs it touches.""" + return _blinding.blind_of(CATALOG_CONFIG, catalogue_name(version)) def fiducial_binning_suffix(fiducial=None): @@ -260,12 +282,17 @@ def base_version(version): return re.sub(r"_ecut\d+", "", re.sub(r"_leak_corr$", "", version)) -def catalogue_entry(version): - """The catalogue-config entry describing ``version`` (its own, or the one - its ``_leak_corr`` variant derives from).""" +def catalogue_name(version): + """The name of the catalogue-config entry describing ``version`` (its own, + or the one its ``_leak_corr`` / ``_seed`` variant derives from).""" if version in CATALOG_CONFIG: - return CATALOG_CONFIG[version] - return CATALOG_CONFIG[re.sub(r"_leak_corr$", "", version)] + return version + return re.sub(r"_seed\d+$", "", re.sub(r"_leak_corr$", "", version)) + + +def catalogue_entry(version): + """The catalogue-config entry describing ``version``.""" + return CATALOG_CONFIG[catalogue_name(version)] def redshift_path(version): diff --git a/workflow/rules/cosmo_val.smk b/workflow/rules/cosmo_val.smk index f4eb0e2a..81d305d1 100644 --- a/workflow/rules/cosmo_val.smk +++ b/workflow/rules/cosmo_val.smk @@ -9,18 +9,18 @@ # catalogue ──→ xi (one job per grid: reporting, integration) # ├─ integration part ─┬──→ pure_eb (part, npz, figures) # │ └──→ cosebis (part, npz, figures) -# └─ reporting .txt ──→ 2pcf plot, ratio_xi_sys_xi +# └─ reporting part ──→ 2pcf plot, ratio_xi_sys_xi # CosmoCov ξ±, integration grid ──→ pure_eb, cosebis (their covariances) # catalogue ──→ pseudo_cl (part) ──┬─→ pseudo-Cl figures # CosmoCov ξ±, reporting grid ─────┤ # rho/tau (part + FITS) ───────────┼─→ summarize_bmodes (reads the products) # └─→ assemble_sacc ──→ {version}.sacc # -# The B-mode rules are ingests: pure_eb, cosebis, the pseudo-Cl figures and the -# summary all work from the parts and the covariance inputs, never from a -# catalogue, so a blinded part keeps everything downstream blinded. The -# analytic covariances (CosmoCov ξ±, NaMaster pseudo-Cℓ) are what assembly puts -# in the terminal file, replacing the estimates a part was born with. +# The B-mode rules and the ξ± figures are ingests: they work from the parts and +# the covariance inputs, never from a catalogue, so a blinded part keeps +# everything downstream blinded. The analytic covariances (CosmoCov ξ±, +# NaMaster pseudo-Cℓ) are what assembly puts in the terminal file, replacing +# the estimates a part was born with. # # Methods that only emit figures, or whose figure paths derive from internal # handler state (rho/tau plots, rho_tau_fits, objectwise leakage, 2pcf @@ -49,15 +49,6 @@ CV_BINNING = ( ) -def cv_xi_txt(version): - """Path to the ξ± TreeCorr dump calculate_2pcf_version writes for a version. - - Mirrors RealSpaceMixin._xi_txt_path for the ("all", "all") pair: - xi_{ver}_tomo_bin_all_minsep=..._maxsep=..._nbins=..._npatch=...txt - """ - return str(COSMO_VAL / f"xi_{version}_tomo_bin_all_{xi_binning('reporting')}.txt") - - def _pure_eb_stub(version): """Shared stem of the pure-E/B diagnostic products (npz + figures).""" eb = XI_GRIDS["integration"] @@ -197,12 +188,11 @@ CV_INIT = cv_init_params(config) # --------------------------------------------------------------------------- # PSF diagnostics: rho/tau statistics and the PSF-error fit # --------------------------------------------------------------------------- -# The rho/tau FITS and the xi data vector are produced by the generic compute -# rules `rho_tau_stats` and `xi` already defined in workflow/rules/twopoint.smk -# (same filename convention, same CosmologyValidation methods). The cosmo_val -# suite consumes those outputs rather than re-declaring colliding rules. The -# helper paths cv_rho_stats / cv_tau_stats / cv_xi_txt resolve to exactly the -# files those rules write, so requesting them triggers the existing compute. +# The rho/tau FITS and the ξ± parts are produced by the generic compute rules +# `rho_tau_stats` and `xi` in workflow/rules/twopoint.smk. The cosmo_val suite +# consumes those outputs rather than re-declaring colliding rules: the helper +# paths cv_rho_stats / cv_tau_stats / cv_xi_sacc resolve to exactly the files +# those rules write, so requesting them triggers the existing compute. rule cv_plot_rho_stats: @@ -319,9 +309,9 @@ rule cv_additive_bias: rule cv_plot_2pcf: - """n_pairs / xi± overlay across versions.""" + """n_pairs / xi± overlay across versions, from the reporting parts.""" input: - xi=[cv_xi_txt(v) for v in CV_VERSIONS], + xi=[cv_xi_sacc(v, "reporting") for v in CV_VERSIONS], output: sentinel=str(CV_SENTINELS / "plot_2pcf.done"), params: @@ -335,7 +325,7 @@ rule cv_plot_2pcf: rule cv_ratio_xi_sys_xi: """Ratio of PSF systematics (xi_psf_sys) to the cosmic-shear signal (xi+).""" input: - xi=[cv_xi_txt(v) for v in CV_VERSIONS], + xi=[cv_xi_sacc(v, "reporting") for v in CV_VERSIONS], rho=[cv_rho_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], tau=[cv_tau_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], output: @@ -502,7 +492,7 @@ rule assemble_sacc: sacc=cv_analysis_sacc("{version}"), params: version="{version}", - type=CV.get("type", "data"), + blind=lambda w: blind_of(w.version), # The statistics this rule wired, so a typo'd input keyword cannot # silently drop one. expected=lambda w: [ diff --git a/workflow/rules/twopoint.smk b/workflow/rules/twopoint.smk index 93abebe4..5d6450e6 100644 --- a/workflow/rules/twopoint.smk +++ b/workflow/rules/twopoint.smk @@ -14,7 +14,7 @@ def xi_binning(grid): rule xi: - """TreeCorr ξ±(θ) for one version on one angular grid. + """TreeCorr ξ±(θ) for one version on one angular grid, as its sealed SACC part. One rule for every grid: outputs are named by their binning, so a request binds the wildcards and the grid label resolves from them. @@ -22,7 +22,6 @@ rule xi: input: catalog=get_shear_catalog, output: - txt=str(COSMO_VAL / "xi_{version}_tomo_bin_all_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.txt"), sacc=str(COSMO_VAL / "{version}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.sacc"), threads: 24 params: @@ -34,6 +33,7 @@ rule xi: cat_config=CAT_CONFIG, output_dir=str(COSMO_VAL), grid=lambda w: grid_of(XI_GRIDS, w), + blind=lambda w: blind_of(w.version), resources: # The fine integration grid needs more memory and wall time than the # ~20-bin reporting one; scale on nbins rather than splitting the rule. @@ -59,6 +59,7 @@ rule rho_tau_stats: npatch="{npatch}", cat_config=CAT_CONFIG, output_dir=str(COSMO_VAL), + blind=lambda w: blind_of(w.version), resources: mem_mb=64000, disk_mb=20000, @@ -86,6 +87,7 @@ rule pseudo_cl: binning="{binning}", nbins=lambda w: int(w.nbins), power=0.5, + blind=lambda w: blind_of(w.version), resources: mem_mb=32000, runtime=120, diff --git a/workflow/scripts/assemble_sacc.py b/workflow/scripts/assemble_sacc.py index 50d1d62e..d756bce2 100644 --- a/workflow/scripts/assemble_sacc.py +++ b/workflow/scripts/assemble_sacc.py @@ -18,6 +18,7 @@ import numpy as np from sp_validation import sacc_io +from sp_validation.blinding import Blind from sp_validation.cosmo_val.sacc_writers import assemble_analysis_sacc # NaMaster iNKA covariance FITS: per-spectrum HDU names, in SACC insertion order. @@ -91,13 +92,16 @@ def assemble_sacc( part_paths, out_path, *, + blind, expected=None, xi_cov=None, pseudo_cl_cov=None, - allow_unblinded=False, ): """Assemble ``{version}.sacc`` from the per-statistic ``part_paths`` mapping. + Every part must be stamped with ``blind``, the catalogue's declared blind, + so a stale part from before a flip is refused. + Parameters ---------- version : str @@ -108,10 +112,10 @@ def assemble_sacc( expected : sequence of str, optional Statistics that must be present, from the caller's config toggles. A typo'd input keyword would otherwise silently drop a statistic. + blind : str + The name of the blind ``version`` is declared under. xi_cov, pseudo_cl_cov Covariance sourcing — see the module docstring. - allow_unblinded : bool, optional - Passed to :func:`sacc_io.load` for every part; ``True`` only for mocks. """ if expected is not None: unknown = [name for name in expected if name not in CANONICAL] @@ -132,12 +136,12 @@ def assemble_sacc( path = part_paths.get(name) if path is None: continue - part = sacc_io.load(path, allow_unblinded=allow_unblinded) - parts.append(_attach_cov(part, name, xi_cov, pseudo_cl_cov)) + parts.append(_attach_cov(sacc_io.load(path), name, xi_cov, pseudo_cl_cov)) if not parts: raise ValueError(f"no parts found for {version}: {part_paths}") - s = assemble_analysis_sacc(parts) - sacc_io.save(s, out_path, type=s.metadata["type"]) + s = sacc_io.save( + assemble_analysis_sacc(parts), out_path, derived_from=parts, blind=Blind(blind) + ) print(f"Assembled {len(parts)} parts -> {out_path}") return s @@ -154,10 +158,10 @@ def _from_snakemake(smk): version=p["version"], part_paths=part_paths, out_path=str(smk.output[0]), + blind=p["blind"], expected=list(p["expected"]), xi_cov=getattr(inp, "xi_cov", None), pseudo_cl_cov=getattr(inp, "pseudo_cl_cov", None), - allow_unblinded=(p.get("type", "data") == "mock"), ) @@ -168,11 +172,7 @@ def _from_cli(argv=None): ap.add_argument("--version", required=True, help="Catalogue version") ap.add_argument("--out", required=True, help="Output {version}.sacc path") ap.add_argument( - "--type", - choices=("data", "mock"), - default="data", - help="Run type. 'mock' reads parts freely; 'data' fails closed on " - "unblinded parts (only concealed/blinded parts load).", + "--blind", required=True, help="The blind the catalogue is declared under" ) for name in CANONICAL: ap.add_argument( @@ -188,9 +188,9 @@ def _from_cli(argv=None): version=a.version, part_paths=part_paths, out_path=a.out, + blind=a.blind, xi_cov=a.xi_cov, pseudo_cl_cov=a.pseudo_cl_cov, - allow_unblinded=(a.type == "mock"), ) diff --git a/workflow/scripts/cv_cosebis.py b/workflow/scripts/cv_cosebis.py index d431b470..ee9d8aef 100644 --- a/workflow/scripts/cv_cosebis.py +++ b/workflow/scripts/cv_cosebis.py @@ -3,7 +3,8 @@ A consumer of the integration-grid ξ± part and the CosmoCov ξ± covariance on the same grid — nothing here touches a catalogue. The covariance goes through the same linear kernel as the modes to give the COSEBIs covariance; it is -analytic, so no Hartlap debiasing applies. +analytic, so no Hartlap debiasing applies. The COSEBIs part is a derivation of +the ξ± part and carries its blind stamp. """ import numpy as np @@ -62,9 +63,7 @@ save_cosebis_results(results, snakemake.output["npz"], fiducial_scale_cut) -# The part inherits the ξ± part's provenance; `type` is re-stamped on save. -metadata = {k: v for k, v in part.metadata.items() if k != "type"} -s = cosebis_to_sacc({0: sacc_io.get_nz(part, 0)}, metadata, fiducial, fiducial_key) -sacc_io.save(s, snakemake.output["sacc"], type="data") +s = cosebis_to_sacc({0: sacc_io.get_nz(part, 0)}, part.metadata, fiducial, fiducial_key) +sacc_io.save(s, snakemake.output["sacc"], derived_from=[part]) verify_outputs(snakemake) diff --git a/workflow/scripts/cv_plot_2pcf.py b/workflow/scripts/cv_plot_2pcf.py index c589edb3..5bb3ac32 100644 --- a/workflow/scripts/cv_plot_2pcf.py +++ b/workflow/scripts/cv_plot_2pcf.py @@ -1,14 +1,15 @@ """Rule cv_plot_2pcf: xi± overlay across versions. -Reads each version's xi txt (declared inputs, produced by the xi rule); calls -plot_2pcf, which reads the existing txt files back rather than recomputing, for -the non-tomographic ("all", "all") pair. Writes figures under the output dir. -Sentinel-tracked: plot_2pcf emits several figures whose names are internal. +Draws each version's reporting ξ± part (the declared inputs) for the +non-tomographic ("all", "all") pair; plot_2pcf draws the parts it is handed +rather than measuring. Writes figures under the output dir. Sentinel-tracked: +plot_2pcf emits several figures whose names are internal. """ -from cv_runner import _unbuffer_streams, make_cv, touch_sentinels +from cv_runner import _unbuffer_streams, make_cv, take_reporting_parts, touch_sentinels _unbuffer_streams() cv = make_cv(snakemake) +take_reporting_parts(cv, snakemake.input["xi"]) cv.plot_2pcf(tomography=False, show=False) touch_sentinels(snakemake) diff --git a/workflow/scripts/cv_pure_eb.py b/workflow/scripts/cv_pure_eb.py index 3f978d3a..3b582699 100644 --- a/workflow/scripts/cv_pure_eb.py +++ b/workflow/scripts/cv_pure_eb.py @@ -2,10 +2,10 @@ A consumer of the integration-grid ξ± part plus a ξ± covariance on that grid — nothing here touches a catalogue. The estimator is one fixed linear operator -on the fine ξ± (b_modes.pure_eb_operator), averaged into the reporting bins with -the part's TreeCorr pair weights into reporting bins snapped onto the fine -edges, so its covariance is the supplied ξ± covariance pushed -exactly through that operator. +on the fine ξ± (b_modes.pure_eb_operator), averaged with the part's TreeCorr +pair weights into reporting bins snapped onto the fine edges, so its covariance +is the supplied ξ± covariance pushed exactly through that operator. The +pure-E/B part is a derivation of the ξ± part and carries its blind stamp. """ import numpy as np @@ -30,17 +30,16 @@ fiducial_scale_cut = tuple(p["fiducial_scale_cut"]) part = sacc_io.load(snakemake.input["xi_integration"]) -theta_int, xip_int, xim_int = sacc_io.get_xi(part, (0, 0), grid="integration") -weight_int = sacc_io.get_xi_weight(part, (0, 0), grid="integration") +gg_int = sacc_io.xi_correlation(part, grid="integration") # A part stores bin centres; the edges come from the grid it was measured on. grid = p["integration"] edges_int = np.geomspace(grid["min_sep"], grid["max_sep"], grid["nbins"] + 1) results = calculate_pure_eb_correlation( - theta_int, - xip_int, - xim_int, - weight_int, + gg_int.meanr, + gg_int.xip, + gg_int.xim, + gg_int.weight, edges_int, np.loadtxt(snakemake.input["cov_integration"]), np.geomspace(p["min_sep"], p["max_sep"], p["nbins"] + 1), @@ -73,15 +72,16 @@ save_pure_eb_results(results, snakemake.output["npz"]) -# The part inherits the ξ± part's provenance; `type` is re-stamped on save. -metadata = {k: v for k, v in part.metadata.items() if k != "type"} -s = pure_eb_to_sacc( - {0: sacc_io.get_nz(part, 0)}, - metadata, - results["theta"], - {key: results[key] for key in sacc_io.PURE_KEYS}, - covariance=results["cov"], +sacc_io.save( + pure_eb_to_sacc( + {0: sacc_io.get_nz(part, 0)}, + part.metadata, + results["theta"], + {key: results[key] for key in sacc_io.PURE_KEYS}, + covariance=results["cov"], + ), + snakemake.output["sacc"], + derived_from=[part], ) -sacc_io.save(s, snakemake.output["sacc"], type="data") verify_outputs(snakemake) diff --git a/workflow/scripts/cv_ratio_xi_sys_xi.py b/workflow/scripts/cv_ratio_xi_sys_xi.py index 3aebd8b8..9e1e4396 100644 --- a/workflow/scripts/cv_ratio_xi_sys_xi.py +++ b/workflow/scripts/cv_ratio_xi_sys_xi.py @@ -1,17 +1,17 @@ """Rule cv_ratio_xi_sys_xi: ratio of PSF systematics to cosmic-shear signal. -Joins two upstream chains for the non-tomographic ("all", "all") pair: the -2pcf data vector (xi txt, read back by calculate_2pcf) and xi_psf_sys -(recomputed in memory from the rho/tau FITS via the lazy cv.xi_psf_sys -property — the PSF-error fit is not persisted). Both are declared as inputs so -the DAG shows the join. Writes ratio_xi_sys_xi.png at a fixed path (declared -output). +Joins two upstream chains for the non-tomographic ("all", "all") pair: each +version's reporting ξ± part and xi_psf_sys (recomputed in memory from the +rho/tau FITS via the lazy cv.xi_psf_sys property — the PSF-error fit is not +persisted). Both are declared as inputs so the DAG shows the join. Writes +ratio_xi_sys_xi.png at a fixed path (declared output). """ -from cv_runner import _unbuffer_streams, make_cv, verify_outputs +from cv_runner import _unbuffer_streams, make_cv, take_reporting_parts, verify_outputs _unbuffer_streams() cv = make_cv(snakemake) +take_reporting_parts(cv, snakemake.input["xi"]) cv.plot_ratio_xi_sys_xi( tomography=False, offset=snakemake.params.get("offset", 0.1), show=False ) diff --git a/workflow/scripts/cv_runner.py b/workflow/scripts/cv_runner.py index 601f2e14..bb67cdfb 100644 --- a/workflow/scripts/cv_runner.py +++ b/workflow/scripts/cv_runner.py @@ -7,7 +7,7 @@ common machinery so each rule script stays a few readable lines. The in-memory ``cv`` object is re-instantiated per rule. The DAG links rules -through the *real* data products each method writes (rho/tau FITS, xi text, +through the *real* data products each method writes (rho/tau FITS, ξ± parts, pseudo-Cl FITS, E/B npz, leakage pkl); the lazy ``cv`` properties that aren't persisted (``c1``/``c2``, ``xi_psf_sys``, ``rho_tau_fits``) are recomputed in whichever rule needs them, which is cheap next to the science compute they @@ -36,6 +36,14 @@ def make_cv(snakemake): return CosmologyValidation(**dict(snakemake.params["cv_init"])) +def take_reporting_parts(cv, paths): + """Hand ``cv`` the reporting ξ± parts at ``paths``, one per version, in order.""" + from sp_validation import sacc_io + + for version, path in zip(cv.versions, paths, strict=True): + cv.xi_parts[version, "reporting"] = sacc_io.load(path) + + def verify_outputs(snakemake): """Fail loudly if any declared output is missing after the method ran.""" missing = [str(o) for o in snakemake.output if not Path(o).exists()] diff --git a/workflow/scripts/generate_pseudo_cl.py b/workflow/scripts/generate_pseudo_cl.py index 53679944..1decf091 100644 --- a/workflow/scripts/generate_pseudo_cl.py +++ b/workflow/scripts/generate_pseudo_cl.py @@ -24,7 +24,6 @@ import json import os -from sp_validation import sacc_io from sp_validation.cosmo_val import CosmologyValidation @@ -47,7 +46,7 @@ def generate_pseudo_cl( Catalog version (e.g., "SP_v1.4.6_leak_corr") out_path : str Exact destination the SACC part is born at — its final (possibly tagged) - name. Skip-if-exists keys on it, so no two rules share a basename. + name. cat_config : str Path to catalog configuration YAML nside : int @@ -128,12 +127,9 @@ def generate_pseudo_cl( # Pseudo-Cls only (no covariance), born directly at the final out_path. cv.calculate_pseudo_cl(compute_tomography=False, out_path=out_path) - if os.path.exists(out_path): - # Readback of the part just written — a legitimate pre-blind consumer. - s = sacc_io.load(out_path, allow_unblinded=True) - ell = sacc_io.get_pseudo_cl(s, (0, 0))[0] - print(f"Generated pseudo-Cl with {len(ell)} ell bins") - print(f"ell range: [{ell.min():.1f}, {ell.max():.1f}]") + ell = cv.pseudo_cls[version]["tomo_bin_all_tomo_bin_all"]["pseudo_cl"]["ELL"] + print(f"Generated pseudo-Cl with {len(ell)} ell bins") + print(f"ell range: [{ell.min():.1f}, {ell.max():.1f}]") return out_path diff --git a/workflow/scripts/run_2pcf.py b/workflow/scripts/run_2pcf.py index 1fcd9a0c..980eaf34 100644 --- a/workflow/scripts/run_2pcf.py +++ b/workflow/scripts/run_2pcf.py @@ -13,13 +13,12 @@ --out The measurement is binning-agnostic: the reporting and the fine integration -grids are the same compute with different ``--min-sep/--max-sep/--nbins``. -``CosmologyValidation.calculate_2pcf_version`` measures the non-tomographic -``("all", "all")`` pair and writes its ``xi_{basename}.txt`` dump (a raw -byproduct); the ξ± data product is born as SACC here, a *part* named by its -binning and tagged with its ``--grid``. The part carries the covariance the -measurement estimated: the dense jackknife covariance when it had patches, the -shot-noise ``varxip``/``varxim`` diagonal when it had none. +grids are the same compute with different ``--min-sep/--max-sep/--nbins``. The +non-tomographic ``("all", "all")`` ξ± is born as a SACC part, named by its +binning, tagged with its ``--grid`` and sealed under the catalogue's blind by +``CosmologyValidation.calculate_2pcf_version``. The part carries the covariance +the measurement estimated: the dense jackknife covariance when it had patches, +the shot-noise ``varxip``/``varxim`` diagonal when it had none. ``output_dir`` is passed explicitly so lc can point each run at its own ``{output}`` tree. @@ -28,11 +27,7 @@ import argparse import os -import numpy as np - -from sp_validation import sacc_io from sp_validation.cosmo_val import CosmologyValidation -from sp_validation.cosmo_val.sacc_writers import xi_to_sacc def run_2pcf( @@ -46,60 +41,42 @@ def run_2pcf( sacc_out=None, grid="reporting", ): - """Measure ξ±(θ) for ``ver`` and write its reporting SACC part. + """Measure ξ±(θ) for ``ver`` and write its sealed SACC part. Parameters mirror the TreeCorr reporting/integration grids: ``min_sep`` / ``max_sep`` in arcmin, ``nbins`` logarithmic bins, ``npatch`` spatial patches (1 for the paper fiducial). ``cat_config`` is an absolute path to the catalog configuration; ``output_dir`` overrides - ``cat_config['paths']['output']`` so the ``.txt`` byproduct lands where lc - expects. ``sacc_out`` is the exact destination for the SACC part (the - Snakemake-declared output); it defaults to a binning-derived name under - the resolved output directory for the CLI path. + ``cat_config['paths']['output']``. ``sacc_out`` is the exact destination for + the part (the Snakemake-declared output); it defaults to a binning-derived + name under the resolved output directory for the CLI path. Returns ------- - treecorr.GGCorrelation - The measured correlation object (also the source of the SACC part). + sacc.Sacc + The part as written. """ cv = CosmologyValidation( versions=[ver], catalog_config=cat_config, output_dir=output_dir, - # so the SACC provenance metadata stamps the npatch actually measured - npatch=npatch, ) - gg = cv.calculate_2pcf_version( + out_path = sacc_out or os.path.join( + output_dir or cv.cc["paths"]["output"], + f"{ver}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.sacc", + ) + cv.calculate_2pcf_version( ver, npatch=npatch, - compute_tomography=False, + grid=grid, + out=out_path, min_sep=min_sep, max_sep=max_sep, nbins=nbins, - )["tomo_bin_all_tomo_bin_all"] - - # Born-as-SACC ξ± part. theta = meanr; theta_nom = rnom. - jackknife = gg.var_method == "jackknife" - s = xi_to_sacc( - cv.sacc_nz(ver), - cv.sacc_metadata(ver), - gg.meanr, - gg.xip, - gg.xim, - grid=grid, - theta_nom=gg.rnom, - npairs=gg.npairs, - weight=gg.weight, - covariance=gg.cov if jackknife else None, - variances=None if jackknife else np.concatenate([gg.varxip, gg.varxim]), - ) - out_path = sacc_out or os.path.join( - output_dir or cv.cc["paths"]["output"], - f"{ver}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.sacc", ) - sacc_io.save(s, out_path, type="data") + part = cv.xi_parts[ver, grid] print(f"Wrote {grid} ξ± SACC part: {out_path}") - return gg + return part def _from_snakemake(smk): @@ -113,8 +90,6 @@ def _from_snakemake(smk): cat_config=p["cat_config"], output_dir=p["output_dir"], grid=p.get("grid", "reporting"), - # The SACC part goes exactly where the rule declares it; the .txt - # byproduct still lands under the resolved output dir. sacc_out=smk.output["sacc"], ) diff --git a/workflow/tests/conftest.py b/workflow/tests/conftest.py index 735d5ebf..202e75f3 100644 --- a/workflow/tests/conftest.py +++ b/workflow/tests/conftest.py @@ -10,10 +10,10 @@ Snakemake is in. The ``candide`` tests need candide itself; CI deselects them. The ``toy`` fixture is a disposable checkout: copies of ``workflow/`` and -``papers/cosmo_val/``, this checkout's ``src/`` symlinked in, a one-catalogue -``cosmo_val/cat_config.yaml``, a touched catalogue file, the processed CosmoCov -covariances already in place (their inputs live on candide), and both output -roots in tmp. +``papers/cosmo_val/``, this checkout's ``src/`` symlinked in, a blinded and an +undeclared catalogue, stand-ins for the processed CosmoCov covariances (their +inputs live on candide), and both output roots in tmp. The host never opens a +blind, so there is no blind record. """ import dataclasses @@ -30,8 +30,10 @@ REPO = Path(__file__).resolve().parents[2] -# The toy catalogue and its leakage-corrected variant. +# The toy catalogue and its leakage-corrected variant: blinded under `toy`. VERSIONS = ("SP_v0.1", "SP_v0.1_leak_corr") +# Declares no blind. +UNDECLARED = "SP_v0.5" @dataclasses.dataclass @@ -103,23 +105,27 @@ class Toy: common: object covariances: dict # (version, "g" | "ng") -> the processed CosmoCov file - def snakemake(self, *args, cwd=None, env=None, timeout=300): + def snakemake(self, *args, config=(), cwd=None, env=None, timeout=300): """Run the host Snakemake in the toy's paper directory, or in ``cwd``. - ``env`` replaces the toy's environment. The workflow profile is off: it - enables apptainer, which CI runners lack. + ``config`` adds ``KEY=VALUE`` overrides. ``env`` replaces the toy's + environment. The workflow profile is off: it enables apptainer, which CI + runners lack. """ + cmd = [ + sys.executable, + "-m", + "snakemake", + "--cores", + "1", + "--workflow-profile", + "none", + *args, + ] + if config: + cmd += ["--config", *config] return subprocess.run( - [ - sys.executable, - "-m", - "snakemake", - "--cores", - "1", - "--workflow-profile", - "none", - *args, - ], + cmd, cwd=cwd or self.rundir, env=env or self.env, text=True, @@ -130,24 +136,35 @@ def snakemake(self, *args, cwd=None, env=None, timeout=300): ) -def _cat_config(catalogue): - entry = { - "subdir": str(catalogue.parent), - "pipeline": "SP", - "colour": "orange", - "marker": "^", - "cov_th": {"A": 100.0, "n_e": 5.0, "sigma_e": 0.3}, - "shear": { - "path": str(catalogue), - "redshift_path": str(catalogue.parent / "nz_SP_v0.1_A.txt"), - "w_col": "w", - "e1_col": "e1", - "e2_col": "e2", - "e1_col_corrected": "e1_leak_corrected", - "e2_col_corrected": "e2_leak_corrected", - }, +def _cat_config(data): + """A blinded and an undeclared catalogue, each reading its own touched file.""" + + def entry(name, **declaration): + catalogue = data / f"{name}.fits" + catalogue.touch() + return { + "subdir": str(data), + "pipeline": "SP", + "colour": "orange", + "marker": "^", + "cov_th": {"A": 100.0, "n_e": 5.0, "sigma_e": 0.3}, + "shear": { + "path": str(catalogue), + "redshift_path": str(data / "nz_SP_v0.1_A.txt"), + "w_col": "w", + "e1_col": "e1", + "e2_col": "e2", + "e1_col_corrected": "e1_leak_corrected", + "e2_col_corrected": "e2_leak_corrected", + }, + **declaration, + } + + return { + VERSIONS[0]: entry(VERSIONS[0], blind="toy"), + UNDECLARED: entry(UNDECLARED), + "paths": {"output": "./output"}, } - return {VERSIONS[0]: entry, "paths": {"output": "./output"}} @pytest.fixture(scope="session") @@ -160,12 +177,10 @@ def toy(tmp_path_factory): ) (root / "src").symlink_to(REPO / "src") - catalogue = root / "data" / "toy_shear.fits" - catalogue.parent.mkdir() - catalogue.touch() + (root / "data").mkdir() (root / "cosmo_val").mkdir() (root / "cosmo_val" / "cat_config.yaml").write_text( - yaml.safe_dump(_cat_config(catalogue)) + yaml.safe_dump(_cat_config(root / "data")) ) rundir = root / "papers" / "cosmo_val" diff --git a/workflow/tests/test_dag.py b/workflow/tests/test_dag.py index 438f91a1..7b1e4e3d 100644 --- a/workflow/tests/test_dag.py +++ b/workflow/tests/test_dag.py @@ -1,12 +1,27 @@ """DAG properties, checked through the host launcher (see conftest.py).""" import os +import re import subprocess import sys from pathlib import Path import pytest -from conftest import REPO, VERSIONS, _load_module, on_candide, parse_jobs +import yaml +from conftest import REPO, UNDECLARED, VERSIONS, _load_module, on_candide, parse_jobs + + +@pytest.fixture(scope="module") +def forced(toy): + """``snakemake -F -n all`` on the toy: its output and every job it schedules.""" + result = toy.snakemake("-F", "-n", "all") + assert result.returncode == 0, result.stdout + return result.stdout, parse_jobs(result.stdout) + + +@pytest.fixture(scope="module") +def grids(toy): + return toy.common.xi_grids(toy.config, toy.config["fiducial"]) def test_assemble_resolves(toy): @@ -37,13 +52,24 @@ def test_assemble_resolves(toy): }, job.input -def test_one_integration_grid(toy): - """COSEBIs and pure-E/B share the integration-grid part and its covariance.""" - result = toy.snakemake("-n", "assemble_sacc_all") - assert result.returncode == 0, result.stdout - jobs = parse_jobs(result.stdout) - grids = toy.common.xi_grids(toy.config, toy.config["fiducial"]) +def test_xi_leaves_a_measurement_only_as_a_part(toy, forced, grids): + """No job reads or writes a ξ± text dump; the ξ± figures draw the parts.""" + jobs = forced[1] + dumps = [ + f for j in jobs for f in j.input + j.output if re.search(r"_xi_.*\.txt$", f) + ] + assert not dumps, dumps + reporting = toy.common.grid_binning(grids["reporting"]) + for rule in ("cv_plot_2pcf", "cv_ratio_xi_sys_xi"): + (job,) = [j for j in jobs if j.rule == rule] + assert {Path(f).name for f in job.input if "_xi_" in Path(f).name} == { + f"{v}_xi_{reporting}.sacc" for v in VERSIONS + }, job.input + +def test_one_integration_grid(toy, forced, grids): + """COSEBIs and pure-E/B share the integration-grid part and its covariance.""" + jobs = forced[1] tag = toy.common.grid_binning(grids["integration"]) for version in VERSIONS: by_rule = { @@ -55,6 +81,47 @@ def test_one_integration_grid(toy): assert by_rule["cv_pure_eb"] == {part, covariance}, by_rule["cv_pure_eb"] +def test_a_blind_flip_reruns_the_catalogues_parts(toy, grids, tmp_path): + """A part's params carry its catalogue's blind, so declaring the catalogue + public reruns the part and nothing else does.""" + env = toy.env | {"COSMO_VAL": str(tmp_path / "cosmo_val")} + reporting = toy.common.grid_binning(grids["reporting"]) + part = tmp_path / "cosmo_val" / f"{VERSIONS[0]}_xi_{reporting}.sacc" + # Records the job's params; the output itself must exist to be touched. + touched = toy.snakemake("--touch", str(part), env=env) + assert touched.returncode == 0, touched.stdout + part.parent.mkdir(exist_ok=True) + part.touch() + + def scheduled(): + result = toy.snakemake("-n", str(part), env=env) + assert result.returncode == 0, result.stdout + return [j.rule for j in parse_jobs(result.stdout)], result.stdout + + rules, output = scheduled() + assert rules == [], output + cat_config = toy.root / "cosmo_val" / "cat_config.yaml" + declared = cat_config.read_text() + flipped = yaml.safe_load(declared) + flipped[VERSIONS[0]]["blind"] = "none" + cat_config.write_text(yaml.safe_dump(flipped)) + try: + rules, output = scheduled() + finally: + cat_config.write_text(declared) + assert rules == ["xi"], output + assert "params have changed" in output.lower(), output + + +def test_a_launch_stops_on_a_catalogue_declaring_no_blind(toy): + result = toy.snakemake( + "-n", "assemble_sacc_all", config=[f'versions=["{UNDECLARED}"]'] + ) + assert result.returncode != 0, result.stdout + assert "declares no `blind:`" in result.stdout, result.stdout + assert "rule assemble_sacc" not in result.stdout + + @pytest.mark.parametrize("named", [True, False], ids=["named", "unnamed"]) def test_outputs_stay_in_the_output_roots(toy, named): """Nothing the suite declares lands outside the configured output roots.