diff --git a/Dockerfile b/Dockerfile index f9272fd3d..c6515bc0d 100644 --- a/Dockerfile +++ b/Dockerfile @@ -124,8 +124,9 @@ RUN bash -c 'set -eo pipefail; \ cmds=$(python -m cosmosis.configure); \ eval "$cmds"; \ export GSL_INC=/usr/include GSL_LIB=/usr/lib/x86_64-linux-gnu; \ + export FFTW_INCLUDE_DIR=/usr/include FFTW_LIBRARY=/usr/lib/x86_64-linux-gnu; \ make -C shear; \ - test -f shear/cl_to_xi_nicaea/nicaea_interface.so' + python -c "import ctypes; ctypes.CDLL(\"shear/cl_to_xi_nicaea/nicaea_interface.so\")"' # Install sp_validation itself (editable) into the same venv; deps are already # satisfied by the sync above. diff --git a/cosmo_inference/README.md b/cosmo_inference/README.md index d45eb148b..5c9dd6e1b 100644 --- a/cosmo_inference/README.md +++ b/cosmo_inference/README.md @@ -1,76 +1,88 @@ -# UNIONS Cosmological Inference Pipeline -by Lisa Goh and Sacha Guerrini, CEA Paris-Saclay +# SACC → CosmoSIS inference -This folder contains the files neccessary to run the cosmological inference pipeline on the UNIONS galaxy catalogues. +The shared workflow converts an analysis SACC to CosmoSIS 2pt-FITS, renders a pipeline INI from explicit template/values/priors inputs, samples under MPI, and writes a weighted JSON summary. +It never reads a catalogue or a blind record. +A concealed SACC stays concealed; its `blind` stamp is recorded in the input provenance. -### Requirements -Everything the pipeline needs ships in the container: nothing to install, and no -paths to edit before a run. +```bash +snakemake --profile workflow/profiles/default -j 32 \ + -s workflow/inference/Snakefile inference_fiducial --configfile run.yaml +``` -[CosmoSIS](https://cosmosis.readthedocs.io/en/latest/) comes in via the -`workflow` extra, built with MPI support. The CosmoSIS Standard Library — the -tree of modules the `.ini` pipelines name — is built into the image at -`/opt/cosmosis-standard-library`, with `CSL_DIR` pointing there. CosmoSIS reads -environment variables into an `.ini`'s `[DEFAULT]` section, so the templates' -`COSMOSIS_DIR = %(CSL_DIR)s` resolves to it. Outside the container, export -`CSL_DIR` at a build of your own and the same templates work unchanged. +Run Snakemake on the host, from the real checkout path. +Its jobs run in the CI dependency image and import the launched checkout. +The standalone entry point needs only the config below; `workflow/Snakefile` includes the same rules for composition with measurement workflows. +`inference_fiducial` gathers all named runs, not just preparation products. +A single run can be requested as `//summary.json`. -CSL is pinned to the **UNIONS-WL org fork** -([UNIONS-WL/cosmosis-standard-library](https://github.com/UNIONS-WL/cosmosis-standard-library/)) -at `b7b1552a`, the fork's main branch: Sacha's four UNIONS commits are -reapplied on current upstream, including the scipy `lpn` fix. +```yaml +inference: + output_dir: /path/to/run/inference + runs: + mock_00001: + sacc: /path/to/MOCK_00001_clean_tomo.sacc + covariance: /path/to/current_reference_tomo.sacc # optional replacement + n_bins: 4 + template: cosmo_inference/cosmosis_config/templates/cosmosis_pipeline_A_ia.ini + values: /path/to/mock_values.ini + priors: /path/to/mock_priors.ini + threads: 32 + sampler: nautilus + sampler_settings: + n_live: 300 + n_eff: 1500 + resume: T + scale_cuts: + XI_PLUS: [12.0, 83.0] + XI_MINUS: [12.0, 83.0] + overrides: + pipeline: + debug: F + timing: F +``` -Launch sampling under MPI (`mpiexec -n N cosmosis --mpi ...`), not `--smp`: -CosmoSIS's shared-memory pool is unmaintained and still crashes after sampling -completes (`Pool` has no attribute `data`, `runtime/process_pool.py`) as of -3.25.2. +Each run produces `data.fits`, `input.json`, `pipeline.ini`, `chain.txt`, `cosmosis.log`, and `summary.json`. +Keep values and priors outside the code: every sampled bin needs its own m and redshift-shift entries. +The example sampler settings are economical mock-validation settings, not a convergence recommendation for a final analysis. +The summary uses Nautilus log weights (or equal weights for unweighted chains) and reports means, standard deviations and central 68% intervals. +MCMC burn-in must be handled by the sampler configuration before summarization. -### To Run -The inference pipeline is orchestrated through Snakemake. On the candide -cluster, drive it with the committed profile — see -[`workflow/README.md`](../workflow/README.md) for the one-time -`uv tool install` setup and the full explanation. From the repository root: +## Inputs and covariance -```bash -snakemake --profile workflow/profiles/candide \ - -s workflow/Snakefile \ - inference_fiducial --configfile -``` +Conversion writes all ξ+ pairs, then all ξ− pairs, in numeric upper-triangle order, with 1-indexed BIN1/BIN2 and one NZ column per source bin. +The covariance follows exactly the same SACC-index permutation, including cross-pair terms. +Integration-grid and unrelated diagnostic points are excluded. +A covariance reference contributes only its covariance, never its means or n(z). +Reference rows are matched by statistic, tracer pair, grid and nominal angle (or measured angle if no nominal bin is stored). +This allows a corrected data covariance to be used on a clean mock with identical nominal bins. -Off-cluster, drop `--profile` and add `-j ` instead. Each job runs -inside the sp_validation container automatically — no `apptainer shell` or -`apptainer exec` needed by hand. +The direct `sacc_to_twopoint_fits` converter retains the historical mixed ξ/Cℓ file layout and single-bin bytes. +CSL's legacy reader cannot associate two separate covariances with such a mixed file, so inference preparation writes only the requested likelihood statistic. +The default is `statistic: xi`. +`statistic: cell` selects a CELL_EE-only likelihood product, but the SAMPLE representation does not propagate SACC bandpower windows: a harmonic inference run must validate its window treatment separately. -This will automatically execute all steps: -1. Calculate 2PCF ($\xi_{pm}$) via `cosmo_val.py` -2. Compute covariance matrices using CosmoCov -3. Prepare CosmoSIS data (FITS) via `cosmosis_fitting.py` -4. Run CosmoSIS inference +For the legacy single-bin PSF-leakage pipeline, add `rho_stats` and `tau_stats` sidecar FITS paths and select `cosmosis_pipeline_A_psf.ini`. +Both sidecars are required; RHO_STATS carries variances absent from the SACC. +The τ covariance comes from the SACC, including τ0–τ2 correlations. +Tomographic PSF sidecars are not supported by that legacy model. +The old ξ-file preparation CLI `scripts/cosmosis_fitting.py` is retired; its FITS builders survive only as a frozen byte-compatibility test oracle. -For standalone FITS data preparation (real-space inputs plus optional pseudo-$C_\ell$ data), you can also use the Python script directly: +## Runtime -```bash -python scripts/cosmosis_fitting.py \ - --cosmosis-root "catalog_version_config" \ - --data-dir "/path/to/output/chains" \ - --nz-file "/path/to/nz_file.txt" \ - --out-file "/path/to/output.fits" \ - --xi "/path/to/xi_plus.fits" "/path/to/xi_minus.fits" \ - --cov-xi "/path/to/covariance.txt" \ - --use-rho-tau \ - --rho-stats "/path/to/rho_stats.fits" \ - --tau-stats "/path/to/tau_stats.fits" \ - --cov-tau "/path/to/cov_tau.npy" \ - --cl-file "/path/to/pseudo_cl.fits" \ - --cov-cl "/path/to/pseudo_cl_cov.fits" -``` - -You can view all available options with: -```bash -python scripts/cosmosis_fitting.py --help -``` +CosmoSIS and the UNIONS-WL CSL fork (`b7b1552a`) ship in the CI image at `/opt/cosmosis-standard-library`. +The renderer resolves `CSL_DIR` explicitly, and the chain rule evaluates `python -m cosmosis.configure` before launch so compiled modules can find `libcosmosis.so`. +Run sampling as `mpiexec -n N cosmosis --mpi`, never `--smp`. +PolyChord requires its optional sampler package; Nautilus is available in the current image. +Use resumable settings for long runs. -Ensure the pseudo-$C_\ell$ spectra (`pseudo_cl_*.fits`) and their covariance (`pseudo_cl_cov_*.fits`) produced by `cosmo_val.py` exist for the requested catalog version (or mock seed) before running the standalone command. +Each chain rule runs on one allocated node and disables Open MPI's attempt to rediscover Slurm, because the image has no `srun` binary. +For manual single-node launches, use the same environment settings: +```bash +export APPTAINERENV_PRTE_MCA_ras='^slurm' +export APPTAINERENV_PRTE_MCA_plm='^slurm' +``` -This is the pipeline used to derive cosmological constraints with cosmic shear data from the UNIONS v1.4 catalogue. +Set `inference.mpi_args: "--bind-to none"` when running several chains concurrently inside the allocation. +Never use this local-launch setting to span unallocated nodes. +The allocation, its CPU limit, and Snakemake's `--cores` remain the resource limits. diff --git a/cosmo_inference/cosmosis_config/templates/cosmosis_pipeline_A_psf.ini b/cosmo_inference/cosmosis_config/templates/cosmosis_pipeline_A_psf.ini index 341baf214..dc88682e3 100644 --- a/cosmo_inference/cosmosis_config/templates/cosmosis_pipeline_A_psf.ini +++ b/cosmo_inference/cosmosis_config/templates/cosmosis_pipeline_A_psf.ini @@ -107,6 +107,8 @@ file = %(COSMOSIS_DIR)s/shear/xi_sys/tau_from_rho.py data_file=%(FITS_FILE)s [2pt_like] +file = %(COSMOSIS_DIR)s/likelihood/2pt/2pt_like_xi_sys.py +add_xi_sys = T data_file=%(FITS_FILE)s gaussian_covariance=F covmat_name=COVMAT diff --git a/cosmo_inference/scripts/cosmosis_fitting.py b/cosmo_inference/scripts/cosmosis_fitting.py deleted file mode 100644 index 9832984c2..000000000 --- a/cosmo_inference/scripts/cosmosis_fitting.py +++ /dev/null @@ -1,789 +0,0 @@ -#!/usr/bin/env python -# coding: utf-8 -"""Prepare CosmoSIS inputs from UNIONS validation outputs. - -The script lives in ``cosmo_inference/scripts``. By default it reads templates -from ``cosmo_inference/cosmosis_config/templates`` and writes data products beneath -``cosmo_inference/data`` and ``cosmo_inference/cosmosis_config/output``. Override -``--template-dir`` or ``--output-root`` to use alternative locations. -""" - -import argparse -import os -import re -import sys -from pathlib import Path - -import numpy as np -from astropy.io import fits - -# ============================================================================ -# SECTION 1: FITS FORMAT CONVERTERS (data → FITS HDU) -# ============================================================================ - - -def _create_2pt_hdu(values, theta, name, quant1, quant2): - """Create standardized 2-point correlation FITS HDU.""" - nbins = len(values) - lst = np.arange(1, nbins + 1) - - col1 = fits.Column(name="BIN1", format="K", array=np.ones(nbins)) - col2 = fits.Column(name="BIN2", format="K", array=np.ones(nbins)) - col3 = fits.Column(name="ANGBIN", format="K", array=lst) - col4 = fits.Column(name="VALUE", format="D", array=values) - col5 = fits.Column(name="ANG", format="D", unit="arcmin", array=theta) - - coldefs = fits.ColDefs([col1, col2, col3, col4, col5]) - hdu = fits.BinTableHDU.from_columns(coldefs, name=name) - - hdu_dict = { - "2PTDATA": "T", - "QUANT1": quant1, - "QUANT2": quant2, - "KERNEL_1": "NZ_SOURCE", - "KERNEL_2": "NZ_SOURCE", - "WINDOWS": "SAMPLE", - } - for key, value in hdu_dict.items(): - hdu.header[key] = value - - return hdu - - -def nz_to_fits(filename): - """Convert n(z) text file to FITS format.""" - line = np.loadtxt(filename, max_rows=1) - nbins = len(line) - 1 - - z_mid = np.loadtxt(filename, usecols=0) - nstep = z_mid[1] - z_mid[0] - - z_low = z_mid - nstep / 2 - z_high = z_mid + nstep / 2 - - col1 = fits.Column(name="Z_LOW", format="D", array=z_low) - col2 = fits.Column(name="Z_MID", format="D", array=z_mid) - col3 = fits.Column(name="Z_HIGH", format="D", array=z_high) - cols = [col1, col2, col3] - - for i in range(nbins): - bin_col = np.loadtxt(filename, usecols=i + 1) - hdu_col = fits.Column(name="BIN%d" % (i + 1), format="D", array=bin_col) - cols.append(hdu_col) - - coldefs = fits.ColDefs(cols) - nz_hdu = fits.BinTableHDU.from_columns(coldefs, name="NZDATA") - - nz_lens_dict = { - "NZDATA": "T ", - "EXTNAME": "NZ_SOURCE", - "NBIN": nbins, - "NZ": len(z_low), - } - - for key, value in nz_lens_dict.items(): - nz_hdu.header[key] = value - - return nz_hdu - - -def treecorr_to_fits(filename1, filename2): - """Load xi+ and xi- from separate TreeCorr FITS files.""" - xiplus_hdu = fits.open(filename1) - ximinus_hdu = fits.open(filename2) - return xiplus_hdu[1], ximinus_hdu[1] - - -def parse_combined_xi_fits(filepath): - """Parse combined FITS file with xip/xim columns (mock format).""" - with fits.open(filepath) as hdul: - xi_data = hdul[1].data - xip_vals = xi_data["xip"] - xim_vals = xi_data["xim"] - xi_theta = xi_data["ANG"] if "ANG" in xi_data.names else xi_data["meanr"] - - xip_hdu = _create_2pt_hdu(xip_vals, xi_theta, "XI_PLUS", "G+R", "G+R") - xim_hdu = _create_2pt_hdu(xim_vals, xi_theta, "XI_MINUS", "G-R", "G-R") - - return xip_hdu, xim_hdu - - -def load_pseudo_cl(cl_file): - """Load pseudo-C_ell spectra from .npy or FITS formats.""" - if cl_file.endswith(".npy"): - cl_block = np.load(cl_file) - if cl_block.shape[0] < 5: - raise ValueError( - f"Unexpected C_ell array shape {cl_block.shape} for {cl_file}" - ) - ell = np.asarray(cl_block[0], dtype=np.float64) - cl_ee = np.asarray(cl_block[1], dtype=np.float64) - cl_bb = np.asarray(cl_block[4], dtype=np.float64) - return ell, cl_ee, cl_bb - - if cl_file.endswith(".fits"): - with fits.open(cl_file) as hdul: - data = hdul[1].data - ell = np.asarray(data["ELL"], dtype=np.float64) - cl_ee = np.asarray(data["EE"], dtype=np.float64) - cl_bb = np.asarray(data["BB"], dtype=np.float64) - return ell, cl_ee, cl_bb - - raise NotImplementedError(f"Cl format not supported: {cl_file}") - - -def cl_to_fits(ell, cl_ee, cl_bb): - """Convert C_ell arrays to CosmoSIS FITS HDUs for EE and BB.""" - nbins = len(ell) - lst = np.arange(1, nbins + 1) - - col1 = fits.Column(name="BIN1", format="K", array=np.ones(nbins)) - col2 = fits.Column(name="BIN2", format="K", array=np.ones(nbins)) - col3 = fits.Column(name="ANGBIN", format="K", array=lst) - col5 = fits.Column(name="ANG", format="D", array=ell) - - col4_ee = fits.Column(name="VALUE", format="D", array=cl_ee) - coldefs_ee = fits.ColDefs([col1, col2, col3, col4_ee, col5]) - cl_ee_hdu = fits.BinTableHDU.from_columns(coldefs_ee, name="CELL_EE") - - cl_ee_dict = { - "2PTDATA": "T", - "QUANT1": "GEF", - "QUANT2": "GEF", - "KERNEL_1": "NZ_SOURCE", - "KERNEL_2": "NZ_SOURCE", - "WINDOWS": "SAMPLE", - } - for key, value in cl_ee_dict.items(): - cl_ee_hdu.header[key] = value - - col4_bb = fits.Column(name="VALUE", format="D", array=cl_bb) - coldefs_bb = fits.ColDefs([col1, col2, col3, col4_bb, col5]) - cl_bb_hdu = fits.BinTableHDU.from_columns(coldefs_bb, name="CELL_BB") - - cl_bb_dict = { - "2PTDATA": "T", - "QUANT1": "GBF", - "QUANT2": "GBF", - "KERNEL_1": "NZ_SOURCE", - "KERNEL_2": "NZ_SOURCE", - "WINDOWS": "SAMPLE", - } - for key, value in cl_bb_dict.items(): - cl_bb_hdu.header[key] = value - - return cl_ee_hdu, cl_bb_hdu - - -def cov_cl_to_fits(cov_file, cov_hdu="COVAR_FULL"): - """Convert pseudo-C_ell covariance to a CosmoSIS ImageHDU.""" - if cov_file.endswith(".fits"): - with fits.open(cov_file) as hdul: - cov_data = np.asarray(hdul[cov_hdu].data, dtype=np.float64) - elif cov_file.endswith(".npy"): - cov_data = np.load(cov_file) - else: - raise NotImplementedError( - f"Unsupported pseudo-Cl covariance format: {cov_file}" - ) - - if cov_data.shape[0] != cov_data.shape[1]: - raise ValueError("Pseudo-Cl covariance matrix must be square") - - cov_hdu = fits.ImageHDU(cov_data, name="COVMAT_CELL") - cov_dict = { - "COVDATA": "True", - "EXTNAME": "COVMAT_CELL", - "NAME_0": "CELL_EE", - "STRT_0": 0, - } - - for key, value in cov_dict.items(): - cov_hdu.header[key] = value - - return cov_hdu - - -def tau_to_fits(filename, theta=None): - """ - Convert tau statistics to FITS format. - If theta provided, override original values for consistency with xi. - """ - tau_stats = fits.getdata(filename) - - if theta is not None: - ang = theta - print("Using provided theta values for tau statistics (forcing consistency)") - else: - ang = tau_stats["theta"] - print("Using original tau theta values") - - tau_0_p_hdu = _create_2pt_hdu(tau_stats["tau_0_p"], ang, "TAU_0_PLUS", "G+R", "P+R") - tau_2_p_hdu = _create_2pt_hdu( - tau_stats["tau_2_p"], ang, "TAU_2_PLUS", "G+R", "SR+R" - ) - - return tau_0_p_hdu, tau_2_p_hdu - - -def rho_to_fits(filename, theta=None): - """ - Convert rho statistics to FITS format. - If theta provided, override original values for consistency with xi. - """ - rho_stat_hdul = fits.open(filename) - rho_stat_hdu = rho_stat_hdul[1].copy() - rho_stat_hdu.name = "RHO_STATS" - - if theta is not None: - print( - "Forcing rho statistics to use provided theta values (forcing consistency)" - ) - rho_stat_hdu.data = rho_stat_hdu.data.copy() - rho_stat_hdu.data["theta"] = theta - else: - print("Using original rho theta values") - - rho_stat_hdul.close() - return rho_stat_hdu - - -def covdat_to_fits( - filename_cov_xi, filename_cov_tau=None, filename_cov_cl=None, cov_hdu=None -): - """ - Convert CosmoCov covariance matrix to FITS format. - - If tau covariance provided, block with xi covariance. - """ - covmat_xi = np.loadtxt(filename_cov_xi) - - if filename_cov_tau is not None: - covmat_tau = np.load(filename_cov_tau) - nbins = int(len(covmat_tau) / 3) - covmat_tau = covmat_tau[: 2 * nbins, : 2 * nbins] - covmat = np.block( - [ - [covmat_xi, np.zeros((len(covmat_xi), len(covmat_tau)))], - [np.zeros((len(covmat_tau), len(covmat_xi))), covmat_tau], - ] - ) - else: - covmat = covmat_xi - - if filename_cov_cl is not None: - if filename_cov_cl.endswith(".fits"): - with fits.open(filename_cov_cl) as hdul: - cov_data = np.asarray(hdul[cov_hdu].data, dtype=np.float64) - elif filename_cov_cl.endswith(".npy"): - cov_data = np.load(filename_cov_cl) - else: - raise NotImplementedError( - f"Unsupported pseudo-Cl covariance format: {filename_cov_cl}" - ) - covmat = np.block( - [ - [covmat, np.zeros((len(covmat), len(cov_data)))], - [np.zeros((len(cov_data), len(covmat))), cov_data], - ] - ) - - if len(covmat) != len(covmat[0]): - raise RuntimeError("Covariance matrix is not square") - - cov_hdu = fits.ImageHDU(covmat) - - cov_dict = { - "COVDATA": "True", - "EXTNAME": "COVMAT", - "NAME_0": "XI_PLUS", - "STRT_0": 0, - "NAME_1": "XI_MINUS", - "STRT_1": int(len(covmat_xi) / 2), - } - - if filename_cov_tau: - cov_dict.update( - { - "NAME_2": "TAU_0_PLUS", - "STRT_2": len(covmat_xi), - "NAME_3": "TAU_2_PLUS", - "STRT_3": len(covmat_xi) + int(len(covmat_tau) / 2), - } - ) - - filename_cov_cl and cov_dict.update( - { - "NAME_4": "CELL_EE", - "STRT_4": len(covmat_xi) + (len(covmat_tau) if filename_cov_tau else 0), - } - ) - - for key, value in cov_dict.items(): - cov_hdu.header[key] = value - - return cov_hdu - - -# ============================================================================ -# SECTION 2: VALIDATION & CONSISTENCY CHECKS -# ============================================================================ - - -def check_meanr_consistency(xi_theta, tau_theta, rho_theta, threshold=1.0): - """Check and report theta consistency across xi, tau, and rho statistics.""" - print("=" * 60) - print("MEANR CONSISTENCY CHECK") - print("=" * 60) - - tau_diff = np.abs((tau_theta - xi_theta) / xi_theta) * 100 - rho_diff = np.abs((rho_theta - xi_theta) / xi_theta) * 100 - - print(f"Xi theta range: {xi_theta.min():.6f} - {xi_theta.max():.6f} arcmin") - print(f"Tau theta range: {tau_theta.min():.6f} - {tau_theta.max():.6f} arcmin") - print(f"Rho theta range: {rho_theta.min():.6f} - {rho_theta.max():.6f} arcmin") - print() - print(f"Max tau-xi relative difference: {tau_diff.max():.3f}%") - print(f"Mean tau-xi relative difference: {tau_diff.mean():.3f}%") - print(f"Max rho-xi relative difference: {rho_diff.max():.3f}%") - print(f"Mean rho-xi relative difference: {rho_diff.mean():.3f}%") - - if tau_diff.max() > threshold or rho_diff.max() > threshold: - print(f"\nWARNING: Meanr differences exceed {threshold}% threshold") - if tau_diff.max() > threshold: - print(f" Tau-xi: max {tau_diff.max():.3f}%") - if rho_diff.max() > threshold: - print(f" Rho-xi: max {rho_diff.max():.3f}%") - else: - print(f"✓ All differences below {threshold}% threshold") - - print("=" * 60) - print() - - -# ============================================================================ -# SECTION 3: COSMOSIS CONFIGURATION GENERATION -# ============================================================================ - - -def _generate_ini_file( - args, - template_base, - priors_file, - values_file, - suffix="", - is_harmonic=False, -): - """Generate a CosmoSIS INI configuration file from template with modifications.""" - template_path = Path(args.template_dir) / template_base - output_path = ( - Path(args.output_config_dir) - / f"cosmosis_pipeline_{args.config_name_base}{suffix}.ini" - ) - - with open(template_path, "r") as f: - config_content = f.read() - - modifications = [] - - relative_fits_file = args.config_relative_fits - default_section = ( - f"[DEFAULT]\nSCRATCH = {args.data_dir}\nFITS_FILE = {relative_fits_file}" - ) - modifications.append((r"^\[DEFAULT\]", default_section)) - - output_section = ( - f"[output]\nfilename = %(SCRATCH)s/{args.cosmosis_root}/samples_" - f"{args.cosmosis_root}{suffix}.txt" - ) - modifications.append((r"^\[output\]", output_section)) - - pipeline_section = ( - f"[pipeline]\nvalues = cosmosis_config/templates/{values_file}\npriors = " - f"cosmosis_config/templates/{priors_file}" - ) - modifications.append((r"^\[pipeline\]", pipeline_section)) - - if not is_harmonic: - if args.use_rho_tau: - like_section = ( - "[2pt_like]\nfile = %(COSMOSIS_DIR)s/likelihood/2pt/2pt_like_xi_sys.py" - "\ndata_sets=XI_PLUS XI_MINUS TAU_0_PLUS TAU_2_PLUS\nadd_xi_sys=T" - ) - else: - like_section = ( - "[2pt_like]\nfile = %(COSMOSIS_DIR)s/likelihood/2pt/2pt_like.py" - "\ndata_sets=XI_PLUS XI_MINUS" - ) - else: - like_section = ( - "[2pt_like]\nfile = %(COSMOSIS_DIR)s/likelihood/2pt/2pt_like.py" - "\ndata_sets=CELL_EE" - ) - - covmat_line = "covmat_name=COVMAT" - - modifications.append((r"^\[2pt_like\]", like_section)) - modifications.append((r"^covmat_name=.*", covmat_line)) - - poly_section = f"[polychord]\npolychord_outfile_root = {args.cosmosis_root}{suffix}" - modifications.append((r"^\[polychord\]", poly_section)) - - test_section = ( - f"[test]\nsave_dir = %(SCRATCH)s/best_fit/{args.cosmosis_root}{suffix}" - ) - modifications.append((r"^\[test\]", test_section)) - - for pattern, replacement in modifications: - config_content = re.sub( - pattern, replacement, config_content, flags=re.MULTILINE - ) - - output_path.parent.mkdir(parents=True, exist_ok=True) - - with open(output_path, "w") as f: - f.write(config_content) - - config_type = "harmonic-space" if is_harmonic else "real-space" - print(f"Prepared CosmoSIS {config_type} configuration file in {output_path}") - print(f"You can now run inference with the command: cosmosis {output_path}") - - -def generate_cosmosis_config(args): - """Generate CosmoSIS INI files (real-space and optional harmonic-space).""" - if args.use_rho_tau: - template_base_realspace = "cosmosis_pipeline_A_psf.ini" - values_file = "values_psf.ini" - else: - template_base_realspace = "cosmosis_pipeline_A_ia.ini" - values_file = "values_ia.ini" - - if args.mock: - priors_file = "priors_mock.ini" - elif args.use_rho_tau: - priors_file = "priors_psf.ini" - else: - priors_file = "priors.ini" - - os.makedirs(args.output_config_dir, exist_ok=True) - - if args.xi: - _generate_ini_file( - args, - template_base_realspace, - priors_file, - values_file, - suffix="", - is_harmonic=False, - ) - - if args.cl_file: - template_base_harmonic = "cosmosis_pipeline_A_ia_cell.ini" - _generate_ini_file( - args, - template_base_harmonic, - priors_file, - values_file, - suffix="_cell", - is_harmonic=True, - ) - - -# ============================================================================ -# SECTION 4: MAIN WORKFLOW -# ============================================================================ - - -def parse_args(): - """Parse command-line arguments for the CosmoSIS preparation workflow. - - The script detects its home directory (``cosmo_inference``) automatically - and uses that as the default anchor for CosmoSIS templates as well as the - generated outputs. Override ``--template-dir`` or ``--output-root`` when - running from a different checkout or writing to scratch space. - """ - cosmo_inference_root = Path(__file__).resolve().parent.parent - - parser = argparse.ArgumentParser( - description="Prepare CosmoSIS inference FITS files from real or mock data. " - "Supports multiple xi input formats (separate files, combined FITS).", - epilog=""" -Example for SP_v1.4.6_leak_corr (real data): - python /n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_inference/scripts/cosmosis_fitting.py \\ - --cosmosis-root "SP_v1.4.6_leak_corr_A_minsep=1.0_maxsep=250.0_nbins=20_npatch=1" \\ - --data-dir "/n09data/guerrini/output_chains/SP_v1.4.6_leak_corr_A_minsep=1.0_maxsep=250.0_nbins=20_npatch=1" \\ - --nz-file "/n17data/sguerrini/UNIONS/WL/nz/v1.4.6/nz_SP_v1.4.6_A.txt" \\ - --output-root "/home/guerrini/sp_validation/cosmo_inference" \\ - --xi "/n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_val/output/xi_plus_SP_v1.4.6_leak_corr_minsep=1.0_maxsep=250.0_nbins=20_npatch=1.fits" \\ - "/n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_val/output/xi_minus_SP_v1.4.6_leak_corr_minsep=1.0_maxsep=250.0_nbins=20_npatch=1.fits" \\ - --cov-xi "/n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_inference/data/covariance/covariance_SP_v1.4.6_leak_corr_A_ng_minsep=1.0_maxsep=250.0_nbins=20_masked/covariance_SP_v1.4.6_leak_corr_A_ng_minsep=1.0_maxsep=250.0_nbins=20_masked_processed.txt" \\ - --use-rho-tau \\ - --rho-stats "/n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_val/output/rho_tau_stats/rho_stats_SP_v1.4.6_leak_corr_minsep=1.0_maxsep=250.0_nbins=20_npatch=1.fits" \\ - --tau-stats "/n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_val/output/rho_tau_stats/tau_stats_SP_v1.4.6_leak_corr_minsep=1.0_maxsep=250.0_nbins=20_npatch=1.fits" \\ - --cov-tau "/n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_val/output/rho_tau_stats/cov_tau_SP_v1.4.6_leak_corr_minsep=1.0_maxsep=250.0_nbins=20_npatch=1_th.npy" \\ - --cl-file "/home/guerrini/sp_validation/cosmo_val/output/pseudo_cl_SP_v1.4.6_leak_corr.fits" \\ - --cov-cl "/home/guerrini/sp_validation/cosmo_val/output/pseudo_cl_cov_g_ng_iNKA_SP_v1.4.6_leak_corr.fits" - -Example for glass mock v0 (mock data): - python /n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_inference/scripts/cosmosis_fitting.py \\ - --mock \\ - --cosmosis-root "glass_mock_v0_00001" \\ - --data-dir "/n09data/guerrini/glass_mock_chains/glass_mock_v1.4.6.3_v2/glass_mock_v0_00001" \\ - --nz-file "/n17data/sguerrini/UNIONS/WL/nz/v1.4.6/nz_SP_v1.4.6_A.txt" \\ - --output-root "/home/guerrini/sp_validation/cosmo_inference" \\ - --output-basename "glass_mocks/glass_mock_v1.4.6.3_v2/v0/glass_mock_00001" \\ - --xi "/n09data/guerrini/glass_mock_v1.4.6.3_v2/results/xi_glass_mock_00001_4096_nbins=20.fits" \\ - --cov-xi "/n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_inference/data/covariance/covariance_SP_v1.4.6_A_ng_minsep=1.0_maxsep=250.0_nbins=20_masked/covariance_SP_v1.4.6_A_ng_minsep=1.0_maxsep=250.0_nbins=20_masked_processed.txt" \\ - --use-rho-tau \\ - --rho-stats "/n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_val/output/rho_tau_stats/rho_stats_SP_v1.4.6_minsep=1.0_maxsep=250.0_nbins=20_npatch=1.fits" \\ - --tau-stats "/n17data/cdaley/unions/pure_eb/results/glass_mock_rhotau_samples/00001/tau_stats_sampled.fits" \\ - --cov-tau "/n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_val/output/rho_tau_stats/cov_tau_SP_v1.4.6_minsep=1.0_maxsep=250.0_nbins=20_npatch=1_th.npy" \\ - --cl-file "/n09data/guerrini/glass_mock_v1.4.6.3_v2/results/cl_glass_mock_00001_4096.npy" \\ - --cov-cl "/home/guerrini/sp_validation/cosmo_val/output/pseudo_cl_cov_g_ng_iNKA_SP_v1.4.6.fits" - """, - formatter_class=argparse.RawDescriptionHelpFormatter, - ) - - parser.add_argument( - "--cosmosis-root", type=str, required=True, help="CosmoSIS root identifier" - ) - parser.add_argument( - "--data-dir", type=str, required=True, help="Output MCMC chain directory" - ) - parser.add_argument("--nz-file", type=str, required=True, help="Path to n(z) file") - parser.add_argument( - "--xi", - nargs="+", - required=False, - help="Xi files: 1 (mock FITS) or 2 (data: xi_plus.fits, xi_minus.fits)", - ) - parser.add_argument( - "--cov-xi", type=str, required=False, help="Xi covariance matrix file" - ) - parser.add_argument( - "--use-rho-tau", - action="store_true", - help="Include PSF systematics (requires --rho-stats, --tau-stats, --cov-tau)", - ) - parser.add_argument( - "--rho-stats", - type=str, - required=False, - help="Path to rho statistics FITS file (required if --use-rho-tau)", - ) - parser.add_argument( - "--tau-stats", - type=str, - required=False, - help="Path to tau statistics FITS file (required if --use-rho-tau)", - ) - parser.add_argument( - "--cov-tau", - type=str, - required=False, - help="Path to tau covariance matrix (required if --use-rho-tau)", - ) - parser.add_argument( - "--cl-file", - type=str, - required=False, - help="Path to C_ell data file (.npy, optional for data and mock)", - ) - parser.add_argument( - "--cov-cl", - type=str, - required=False, - help="Path to pseudo-C_ell covariance matrix (required if --cl-file)", - ) - parser.add_argument("--mock", action="store_true", help="Mock data mode") - parser.add_argument( - "--output-root", - type=str, - default=str(cosmo_inference_root), - help=( - "Base directory for generated outputs (defaults to the cosmo_inference " - "folder containing this script)." - ), - ) - parser.add_argument( - "--output-basename", - type=str, - required=False, - help=( - "Optional override for the output sub-directory/basename inside --output-root. " - "Defaults to --cosmosis-root." - ), - ) - parser.add_argument( - "--template-dir", - type=str, - default=str(cosmo_inference_root / "cosmosis_config" / "templates"), - help=( - "Directory containing CosmoSIS template INI files (defaults to the " - "cosmosis_config folder next to this script)." - ), - ) - - return parser.parse_args() - - -if __name__ == "__main__": - args = parse_args() - - try: - output_basename = args.output_basename or args.cosmosis_root - output_root_path = Path(args.output_root).expanduser().resolve() - template_dir_path = Path(args.template_dir).expanduser().resolve() - output_basename_path = Path(output_basename) - data_dir_root = output_root_path / "data" / output_basename_path - config_dir_root = output_root_path / "cosmosis_config" / "output" - data_dir_root.mkdir(parents=True, exist_ok=True) - config_dir_root.mkdir(parents=True, exist_ok=True) - out_file_path = data_dir_root / f"cosmosis_{args.cosmosis_root}.fits" - args.output_root = str(output_root_path) - args.output_config_dir = str(config_dir_root) - args.output_basename = str(output_basename_path) - args.template_dir = str(template_dir_path) - args.config_name_base = str(output_basename_path).replace("/", "_") - args.config_relative_fits = str( - Path("data") / output_basename_path / f"cosmosis_{args.cosmosis_root}.fits" - ) - - print("=" * 60) - print("COSMOSIS_FITTING.PY") - print("=" * 60) - print(f"cosmosis_root: {args.cosmosis_root}") - print(f"data_dir: {args.data_dir}") - print(f"nz_file: {args.nz_file}") - print(f"output_root: {args.output_root}") - print(f"output_basename: {args.output_basename}") - print(f"out_file: {out_file_path}") - if args.xi: - print(f"xi files: {args.xi}") - if args.cov_xi: - print(f"cov_xi: {args.cov_xi}") - print(f"use_rho_tau: {args.use_rho_tau}") - if args.use_rho_tau: - print(f"rho_stats: {args.rho_stats}") - print(f"tau_stats: {args.tau_stats}") - print(f"cov_tau: {args.cov_tau}") - if args.cl_file: - print(f"cl_file: {args.cl_file}") - if args.cov_cl: - print(f"cov_cl: {args.cov_cl}") - print("=" * 60) - print() - - if args.xi and not args.cov_xi: - raise ValueError("--cov-xi is required when --xi files are provided") - if args.cov_xi and not args.xi: - raise ValueError("--xi files are required when --cov-xi is provided") - if args.use_rho_tau and not args.xi: - raise ValueError("--xi files are required when --use-rho-tau is set") - if args.use_rho_tau and not all([args.rho_stats, args.tau_stats, args.cov_tau]): - raise ValueError( - "--use-rho-tau requires: --rho-stats, --tau-stats, --cov-tau" - ) - if args.cl_file and not args.cov_cl: - raise ValueError("--cov-cl is required when --cl-file is provided") - if args.cov_cl and not args.cl_file: - raise ValueError("--cl-file is required when --cov-cl is provided") - - os.makedirs(args.data_dir, exist_ok=True) - cov_hdu = None - xip_hdu = None - xim_hdu = None - if args.xi: - print("Loading xi correlation functions...") - if args.mock: - # Mock mode expects a single combined FITS file - xip_hdu, xim_hdu = parse_combined_xi_fits(args.xi[0]) - else: - # Data mode expects two TreeCorr files (xi_plus, xi_minus) - if len(args.xi) != 2: - raise ValueError( - f"Data mode requires exactly 2 xi files, got {len(args.xi)}" - ) - xip_hdu, xim_hdu = treecorr_to_fits(args.xi[0], args.xi[1]) - - xi_theta = xip_hdu.data["ANG"] - print(f"Loaded xi: {len(xip_hdu.data)} bins") - - print("Loading covariance matrix...") - cov_hdu = covdat_to_fits(args.cov_xi, filename_cov_tau=None) - print(f"Loaded covariance: shape {cov_hdu.data.shape}") - - print("Loading n(z)...") - nz_hdu = nz_to_fits(args.nz_file) - print("Loaded n(z)") - - cl_ee_hdu = None - cl_bb_hdu = None - cov_cl_hdu = None - if args.cl_file: - print("Loading Cl data...") - ell, cl_ee, cl_bb = load_pseudo_cl(args.cl_file) - print(f"Loaded Cl: {len(ell)} multipoles") - cl_ee_hdu, cl_bb_hdu = cl_to_fits(ell, cl_ee, cl_bb) - cov_cl_hdu = cov_cl_to_fits(args.cov_cl, cov_hdu="COVAR_FULL") - print("Loaded pseudo-Cl covariance") - - rho_hdu = None - tau_0_p_hdu = None - tau_2_p_hdu = None - if args.use_rho_tau: - print("Loading rho/tau statistics...") - tau_stats = fits.getdata(args.tau_stats) - rho_stats = fits.getdata(args.rho_stats) - tau_theta = tau_stats["theta"] - rho_theta = rho_stats["theta"] - - check_meanr_consistency(xi_theta, tau_theta, rho_theta, threshold=5.0) - print("✓ Forcing rho and tau to use xi meanr values for consistency") - - rho_hdu = rho_to_fits(args.rho_stats, theta=xi_theta) - tau_0_p_hdu, tau_2_p_hdu = tau_to_fits(args.tau_stats, theta=xi_theta) - print("Loaded rho/tau statistics") - - if args.cl_file: - cov_hdu = covdat_to_fits( - args.cov_xi, - filename_cov_tau=args.cov_tau, - filename_cov_cl=args.cov_cl, - cov_hdu="COVAR_FULL", - ) - else: - cov_hdu = covdat_to_fits(args.cov_xi, filename_cov_tau=args.cov_tau) - print("Loaded combined covariance with tau") - - pri_hdr = fits.Header() - pri_hdu = fits.PrimaryHDU(header=pri_hdr) - - print("Assembling FITS file...") - hdu_list = [pri_hdu, nz_hdu] - if cov_hdu is not None: - hdu_list.append(cov_hdu) - else: - cov_cl_hdu.header["EXTNAME"] = ( - "COVMAT" # If no xi use covmat for the extname of the cl cov - ) - if cov_cl_hdu is not None: - hdu_list.append(cov_cl_hdu) - if args.xi: - hdu_list.extend([xip_hdu, xim_hdu]) - if args.cl_file: - hdu_list.extend([cl_ee_hdu]) - - if args.use_rho_tau: - hdu_list.extend([tau_0_p_hdu, tau_2_p_hdu, rho_hdu]) - - hdul = fits.HDUList(hdu_list) - hdul.writeto(str(out_file_path), overwrite=True) - print(f"✓ FITS file written to {out_file_path}") - print() - - generate_cosmosis_config(args) - - except Exception as e: - print(f"ERROR: {e}", file=sys.stderr) - sys.exit(1) diff --git a/papers/bmodes/config/config.yaml b/papers/bmodes/config/config.yaml index 5accc706e..b8b3280ec 100644 --- a/papers/bmodes/config/config.yaml +++ b/papers/bmodes/config/config.yaml @@ -83,15 +83,6 @@ harmonic: binning: powspace nbins: 32 -# Cosmological inference data-product locations (dormant subsystem). Output -# products land beside the code (workflow/common.py COSMO_INFERENCE); only the -# external chain directories are configured here. -inference: - # CosmoSIS chain output root for real-data runs. - chains_dir: "/n09data/guerrini/output_chains" - # CosmoSIS chain output root for GLASS mock runs. - glass_mock_chains_dir: "/n09data/guerrini/glass_mock_chains" - # COSEBIS calculation parameters cosebis: theta_min: 1.0 # arcmin, angular scale range diff --git a/papers/cosmo_val/config/config.yaml b/papers/cosmo_val/config/config.yaml index 99f553c98..46fea5b0d 100644 --- a/papers/cosmo_val/config/config.yaml +++ b/papers/cosmo_val/config/config.yaml @@ -128,11 +128,6 @@ harmonic: binning: powspace nbins: 32 -# Cosmological inference data-product locations (dormant subsystem). -inference: - chains_dir: "/n09data/guerrini/output_chains" - glass_mock_chains_dir: "/n09data/guerrini/glass_mock_chains" - cosebis: theta_min: 1.0 theta_max: 250.0 diff --git a/src/sp_validation/inference.py b/src/sp_validation/inference.py new file mode 100644 index 000000000..c53823731 --- /dev/null +++ b/src/sp_validation/inference.py @@ -0,0 +1,172 @@ +"""SACC inputs, CosmoSIS configuration and weighted chain summaries. + +The workflow owns execution; values, priors and sampler settings are run inputs. +No catalogue or blind record is read here. +""" + +import configparser +import json +import os +from pathlib import Path + +import numpy as np +from astropy.io import fits + +from . import sacc_io as sio + + +def prepare( + sacc_path, + output, + *, + n_bins=1, + statistic="xi", + covariance=None, + rho_stats=None, + tau_stats=None, +): + """Convert the selected reporting statistic, optionally replacing its covariance. + + A covariance reference contributes only its covariance, never its signal or + n(z). Selection keys enforce the same point order in the two products. + A ξ-only file is intentional: CSL's twopoint reader cannot attach separate + COVMAT and COVMAT_CELL matrices to a mixed ξ/Cℓ file. + """ + s = sio.load(sacc_path) + types = {sio.XI_PLUS, sio.XI_MINUS} + if statistic == "cell": + types.add(sio.CL_EE) + if rho_stats is not None: + types.update(sio.TAU_PLUS.format(k=k) for k in (0, 2)) + rows = [ + i + for i, p in enumerate(s.data) + if p.data_type in types and p.tags.get("grid") == "reporting" + ] + s.keep_indices(rows) + if covariance: + ref = sio.load(covariance) + + def key(p): + return ( + p.data_type, + p.tracers, + p.tags.get("grid"), + p.tags.get("theta_nom", p.tags.get("theta")), + p.tags.get("ell"), + ) + + lookup = {key(p): i for i, p in enumerate(ref.data)} + indices = [lookup[key(p)] for p in s.data] + s.add_covariance(ref.covariance.get_block(indices), overwrite=True) + Path(output).parent.mkdir(parents=True, exist_ok=True) + sidecars = {} + for name, path in (("rho_stats_hdu", rho_stats), ("tau_stats_hdu", tau_stats)): + if path: + with fits.open(path) as hdus: + sidecars[name] = hdus[1].copy() + hdus = sio.sacc_to_twopoint_fits(s, output, n_bins=n_bins, **sidecars) + if statistic == "cell": + # Separate likelihood files avoid the legacy reader's mixed-covariance + # limitation; conversion itself retains the historical single-bin bytes. + del hdus["XI_PLUS"] + del hdus["XI_MINUS"] + del hdus["COVMAT"] + hdus.writeto(output, overwrite=True) + return { + "sacc": str(sacc_path), + "covariance": str(covariance or sacc_path), + "blind": sio.stamp(s), + "n_bins": n_bins, + "statistic": statistic, + } + + +def render( + template, + output, + *, + fits_path, + values, + priors, + chain, + n_bins=1, + sampler="nautilus", + settings=None, + scale_cuts=None, + overrides=None, +): + """Render a small run INI without changing the supplied values or priors.""" + ini = configparser.ConfigParser(interpolation=None) + ini.read(template) + ini["DEFAULT"]["CSL_DIR"] = os.environ.get( + "CSL_DIR", "/opt/cosmosis-standard-library" + ) + ini["DEFAULT"]["SCRATCH"] = str(Path(chain).resolve().parent) + ini["DEFAULT"]["FITS_FILE"] = str(Path(fits_path).resolve()) + ini["pipeline"]["values"] = str(Path(values).resolve()) + ini["pipeline"]["priors"] = str(Path(priors).resolve()) + ini["runtime"]["sampler"] = sampler + ini["output"]["filename"] = str(Path(chain).resolve()) + if sampler not in ini: + ini[sampler] = {} + for key, value in (settings or {}).items(): + ini[sampler][key] = str(value) + # Remove the Paper-III single-bin cuts; every run supplies its own cuts. + like = ini["2pt_like"] + for key in list(like): + if key.startswith("angle_range_"): + del like[key] + for kind, limits in (scale_cuts or {}).items(): + for i in range(1, n_bins + 1): + for j in range(i, n_bins + 1): + like[f"angle_range_{kind}_{i}_{j}"] = " ".join( + str(float(x)) for x in limits + ) + for section, options in (overrides or {}).items(): + if section not in ini: + ini[section] = {} + for key, value in options.items(): + ini[section][key] = str(value) + Path(output).parent.mkdir(parents=True, exist_ok=True) + Path(chain).parent.mkdir(parents=True, exist_ok=True) + with open(output, "w") as stream: + ini.write(stream) + + +def summarize(chain, output): + """Weighted means, standard deviations and central 68% intervals. + + Nautilus writes log weights; emcee and other equal-weight chains do not. + Burn-in for an MCMC run must be handled by its sampler configuration. + """ + with open(chain) as stream: + names = stream.readline().lstrip("#").lower().split() + data = np.loadtxt(chain, ndmin=2) + columns = dict(zip(names, data.T)) + if "log_weight" in columns: + logw = columns["log_weight"] + weights = np.exp(logw - logw.max()) + else: + weights = columns.get("weight", np.ones(len(data))) + weights = weights / weights.sum() + result = { + "samples": len(data), + "effective_samples": float(1 / (weights @ weights)), + "parameters": {}, + } + for name, x in columns.items(): + if "--" not in name: + continue + mean = weights @ x + order = np.argsort(x) + q = np.interp([0.15865, 0.5, 0.84135], np.cumsum(weights[order]), x[order]) + result["parameters"][name] = { + "mean": float(mean), + "sigma": float(np.sqrt(weights @ (x - mean) ** 2)), + "q16": float(q[0]), + "median": float(q[1]), + "q84": float(q[2]), + } + Path(output).write_text(json.dumps(result, indent=2) + "\n") + return result diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index 2ae123220..eb2536940 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -76,11 +76,11 @@ the CosmoSIS "2pt FITS" (``sacc_to_twopoint_fits``) and OneCovariance's redshift / covariance files (``write_nz``, ``nz_config_stanza``, ``read_nz``, ``covariance_blocks``). The - 2pt-FITS converter reproduces today's hand-assembled product from - ``cosmo_inference/scripts/cosmosis_fitting.py`` HDU-for-HDU and - byte-for-byte (verified against that writer), and is single-bin - only today — it fails fast on a multi-bin SACC (tomographic - emission lands with the tomographic round). The OneCovariance + 2pt-FITS converter reproduces the historical hand-assembled product + from the retired ``cosmosis_fitting.py`` HDU-for-HDU and + byte-for-byte for one bin (verified against that writer), and + permutes tomographic products to type-major, pair-ordered layout. + The OneCovariance converters are coupled to SACC by file format only; OneCovariance itself is not a dependency. """ @@ -1091,7 +1091,9 @@ def load(path, blind=None): } -def _twopoint_hdu(name, values, ang, *, ang_unit=None): +def _twopoint_hdu( + name, values, ang, *, ang_unit=None, bin1=None, bin2=None, angbin=None +): """Build one 2pt BinTableHDU (BIN1/BIN2/ANGBIN/VALUE/ANG). Reproduces ``cosmosis_fitting.py._create_2pt_hdu`` /``cl_to_fits`` exactly: @@ -1100,10 +1102,15 @@ def _twopoint_hdu(name, values, ang, *, ang_unit=None): column ("arcmin" for real-space ξ/τ; unset for Cℓ, whose ANG is ℓ). """ nbins = len(values) - angbin = np.arange(1, nbins + 1) + if angbin is None: + angbin = np.arange(1, nbins + 1) columns = [ - fits.Column(name="BIN1", format="K", array=np.ones(nbins)), - fits.Column(name="BIN2", format="K", array=np.ones(nbins)), + fits.Column( + name="BIN1", format="K", array=np.ones(nbins) if bin1 is None else bin1 + ), + fits.Column( + name="BIN2", format="K", array=np.ones(nbins) if bin2 is None else bin2 + ), fits.Column(name="ANGBIN", format="K", array=angbin), fits.Column(name="VALUE", format="D", array=values), fits.Column(name="ANG", format="D", unit=ang_unit, array=ang), @@ -1184,41 +1191,38 @@ def _cov_hdu(matrix, block_names, block_starts, extname="COVMAT", name_in_ctor=F return hdu -def _type_major_xi(s, bins): - """Return ``(theta, xip, xim)`` for one bin pair from the SACC reporting grid. - - ``get_xi`` already returns each statistic in insertion (= ascending - θ) order; the type-major split (all ξ+, then all ξ−) is exactly the two - arrays it hands back, so no further permutation is needed for a single pair. - """ - return get_xi(s, bins, grid="reporting") - +def _twopoint_series(s, name, dtype, pairs, *, angle, ang_unit=None): + """One type-major table and its SACC indices, in numeric bin-pair order.""" + indices, bin1, bin2, angbin = [], [], [], [] + for i, j in pairs: + idx = _indices(s, dtype, _pair((i, j)), grid="reporting") + indices.extend(idx) + bin1.extend([i + 1] * len(idx)) + bin2.extend([j + 1] * len(idx)) + angbin.extend(range(1, len(idx) + 1)) + indices = np.asarray(indices, dtype=int) + hdu = _twopoint_hdu( + name, + s.mean[indices], + [s.data[k].tags[angle] for k in indices], + ang_unit=ang_unit, + bin1=bin1, + bin2=bin2, + angbin=angbin, + ) + return hdu, indices -def _require_single_bin(s, n_bins): - """Fail fast unless the SACC is a valid single-bin ξ product. - The converter emits the single-bin 2pt-FITS today's CosmoSIS pipeline reads - (BIN1/BIN2 all 1, one NZ column). A tomographic SACC would otherwise slip - through silently — ``n_bins`` alone drives the NZDATA column count while the - ξ/covariance are read from bin ``(0, 0)`` only, so a 2-bin file would emit a - ``NBIN=2`` n(z) beside a data vector holding just the ``(0, 0)`` pair. - Guards both the empty-ξ case and the single-bin contract; tomographic - emission lands with the tomographic round. - """ - pairs = s.get_tracer_combinations(XI_PLUS) - if not pairs: - raise ValueError( - f"SACC has no {XI_PLUS} points — nothing to convert; the " - "2pt-FITS data vector is built from the ξ± statistics" - ) - expected = (source_name(0), source_name(0)) - if n_bins != 1 or set(pairs) != {expected}: - raise ValueError( - f"converter is single-bin only (n_bins=1, ξ pairs == {{{expected}}}); " - f"got n_bins={n_bins} and ξ pairs {sorted(pairs)}. Tomographic " - "emission (multiple bin pairs, per-pair BIN1/BIN2, one NZ column per " - "bin) lands with the tomographic round." - ) +def _source_pairs(s, n_bins): + """Require the complete upper triangle of reporting ξ source pairs.""" + indices = s.indices(XI_PLUS, grid="reporting") + if len(indices) == 0: + raise ValueError(f"SACC has no reporting {XI_PLUS} — nothing to convert") + pairs = [(i, j) for i in range(n_bins) for j in range(i, n_bins)] + actual = {s.data[k].tracers for k in indices} + if actual != {_pair(p) for p in pairs}: + raise ValueError(f"n_bins={n_bins} does not match reporting ξ pairs {actual}") + return pairs def sacc_to_twopoint_fits( @@ -1231,7 +1235,7 @@ def sacc_to_twopoint_fits( ): """Convert an analysis SACC to a CosmoSIS 2pt-FITS file. - The assembled ``HDUList`` matches today's ``cosmosis_fitting.py`` product + The assembled ``HDUList`` matches the legacy ``cosmosis_fitting.py`` product for the configuration the SACC describes: PRIMARY, NZ_SOURCE, COVMAT, then (if present) COVMAT_CELL, XI_PLUS, XI_MINUS, (if present) CELL_EE / CELL_BB, and (if the rho/tau sidecars are supplied) TAU_0_PLUS, TAU_2_PLUS, @@ -1252,11 +1256,9 @@ def sacc_to_twopoint_fits( alone cannot rebuild the ``varrho_*`` columns Sacha's fork reads. When omitted, a pure ξ (± Cℓ) product is written. n_bins : int, optional - Number of source tomographic bins. Must be ``1``: this converter emits - the single-bin 2pt-FITS today's CosmoSIS pipeline consumes. Tomographic - emission (multiple bin pairs, per-pair BIN1/BIN2, one NZ column per bin) - lands with the tomographic round; the converter fails fast on anything - else rather than silently truncating to bin ``(0, 0)``. + Number of source bins (default 1). Emits the complete upper triangle + in numeric pair order within each type, with 1-indexed BIN1/BIN2 and + one NZ column per bin. The legacy ρ/τ sidecars support one bin only. Returns ------- @@ -1266,27 +1268,29 @@ def sacc_to_twopoint_fits( Raises ------ ValueError - If the SACC has no ξ points; if ``n_bins != 1`` or the SACC's ξ tracer - pairs are anything other than exactly ``{(source_0, source_0)}`` (the - single-bin contract); or if exactly one of the ρ/τ sidecars is supplied. + If reporting ξ pairs do not match ``n_bins``, or if the legacy ρ/τ + sidecars are supplied incompletely or for multiple source bins. """ if (rho_stats_hdu is None) != (tau_stats_hdu is None): raise ValueError( "rho_stats_hdu and tau_stats_hdu must be supplied together " "(the ρ/τ product needs both, or neither for a pure-ξ product)" ) - _require_single_bin(s, n_bins) + pairs = _source_pairs(s, n_bins) use_rho_tau = rho_stats_hdu is not None - bins = (0, 0) + if use_rho_tau and n_bins != 1: + raise ValueError("the legacy ρ/τ sidecars support one source bin only") nz_hdu = _nz_hdu(s, n_bins) - theta, xip, xim = _type_major_xi(s, bins) - xip_hdu = _twopoint_hdu("XI_PLUS", xip, theta, ang_unit="arcmin") - xim_hdu = _twopoint_hdu("XI_MINUS", xim, theta, ang_unit="arcmin") - - cell_hdu, cov_cell_hdu = _build_cell(s, bins) - - cov_hdu = _build_covmat(s, bins, use_rho_tau=use_rho_tau) + xip_hdu, idx_p = _twopoint_series( + s, "XI_PLUS", XI_PLUS, pairs, angle="theta", ang_unit="arcmin" + ) + xim_hdu, idx_m = _twopoint_series( + s, "XI_MINUS", XI_MINUS, pairs, angle="theta", ang_unit="arcmin" + ) + theta = xip_hdu.data["ANG"] + cell_hdu, cov_cell_hdu = _build_cell(s, pairs) + cov_hdu = _build_covmat(s, idx_p, idx_m, use_rho_tau=use_rho_tau) tau_hdus, rho_hdu = _build_rho_tau(rho_stats_hdu, tau_stats_hdu, theta, use_rho_tau) @@ -1306,7 +1310,7 @@ def sacc_to_twopoint_fits( return hdul -def _build_cell(s, bins): +def _build_cell(s, pairs): """Build the CELL_EE 2pt HDU plus the COVMAT_CELL HDU from the SACC pseudo-Cℓ. Returns ``(None, None)`` when the SACC has no pseudo-Cℓ. Only CELL_EE is @@ -1320,17 +1324,15 @@ def _build_cell(s, bins): if CL_EE not in s.get_data_types(): return None, None - ell, cl_ee, _cl_bb, _cl_eb, _window = get_pseudo_cl(s, bins) - cell_hdu = _twopoint_hdu("CELL_EE", cl_ee, ell) - cell_idx = _indices(s, CL_EE, _pair(bins), grid="reporting") - cov_cell = s.covariance.dense[np.ix_(cell_idx, cell_idx)] + cell_hdu, cell_idx = _twopoint_series(s, "CELL_EE", CL_EE, pairs, angle="ell") + cov_cell = s.covariance.get_block(cell_idx) cov_cell_hdu = _cov_hdu( cov_cell, ["CELL_EE"], [0], extname="COVMAT_CELL", name_in_ctor=True ) return cell_hdu, cov_cell_hdu -def _build_covmat(s, bins, *, use_rho_tau): +def _build_covmat(s, idx_p, idx_m, *, use_rho_tau): """Assemble the blocked COVMAT (ξ± type-major, then the τ blocks). The ξ covariance is pulled from the SACC as the contiguous ξ+/ξ− block for @@ -1338,12 +1340,9 @@ def _build_covmat(s, bins, *, use_rho_tau): ``use_rho_tau`` the τ_0/τ_2 covariance blocks are appended block-diagonally with zero ξ↔τ cross-blocks, exactly as ``covdat_to_fits`` builds them. """ - pair = _pair(bins) - idx_p = _indices(s, XI_PLUS, pair, grid="reporting") - idx_m = _indices(s, XI_MINUS, pair, grid="reporting") n_theta = len(idx_p) xi_idx = np.concatenate([idx_p, idx_m]) # type-major permutation - xi_cov = s.covariance.dense[np.ix_(xi_idx, xi_idx)] + xi_cov = s.covariance.get_block(xi_idx) names = ["XI_PLUS", "XI_MINUS"] starts = [0, n_theta] @@ -1359,7 +1358,7 @@ def _build_covmat(s, bins, *, use_rho_tau): idx_tau0 = _indices(s, TAU_PLUS.format(k=0), tau_pair, grid="reporting") idx_tau2 = _indices(s, TAU_PLUS.format(k=2), tau_pair, grid="reporting") tau_idx = np.concatenate([idx_tau0, idx_tau2]) - tau_cov = s.covariance.dense[np.ix_(tau_idx, tau_idx)] + tau_cov = s.covariance.get_block(tau_idx) matrix = _block_diag(matrix, tau_cov) names += ["TAU_0_PLUS", "TAU_2_PLUS"] starts += [2 * n_theta, 2 * n_theta + len(idx_tau0)] diff --git a/src/sp_validation/tests/legacy_twopoint_reference.py b/src/sp_validation/tests/legacy_twopoint_reference.py new file mode 100644 index 000000000..9ed3efe21 --- /dev/null +++ b/src/sp_validation/tests/legacy_twopoint_reference.py @@ -0,0 +1,286 @@ +"""Frozen pre-SACC FITS builders; byte-compatibility oracle, never production code.""" + +import numpy as np +from astropy.io import fits + + +def _create_2pt_hdu(values, theta, name, quant1, quant2): + """Create standardized 2-point correlation FITS HDU.""" + nbins = len(values) + lst = np.arange(1, nbins + 1) + + col1 = fits.Column(name="BIN1", format="K", array=np.ones(nbins)) + col2 = fits.Column(name="BIN2", format="K", array=np.ones(nbins)) + col3 = fits.Column(name="ANGBIN", format="K", array=lst) + col4 = fits.Column(name="VALUE", format="D", array=values) + col5 = fits.Column(name="ANG", format="D", unit="arcmin", array=theta) + + coldefs = fits.ColDefs([col1, col2, col3, col4, col5]) + hdu = fits.BinTableHDU.from_columns(coldefs, name=name) + + hdu_dict = { + "2PTDATA": "T", + "QUANT1": quant1, + "QUANT2": quant2, + "KERNEL_1": "NZ_SOURCE", + "KERNEL_2": "NZ_SOURCE", + "WINDOWS": "SAMPLE", + } + for key, value in hdu_dict.items(): + hdu.header[key] = value + + return hdu + + +def nz_to_fits(filename): + """Convert n(z) text file to FITS format.""" + line = np.loadtxt(filename, max_rows=1) + nbins = len(line) - 1 + + z_mid = np.loadtxt(filename, usecols=0) + nstep = z_mid[1] - z_mid[0] + + z_low = z_mid - nstep / 2 + z_high = z_mid + nstep / 2 + + col1 = fits.Column(name="Z_LOW", format="D", array=z_low) + col2 = fits.Column(name="Z_MID", format="D", array=z_mid) + col3 = fits.Column(name="Z_HIGH", format="D", array=z_high) + cols = [col1, col2, col3] + + for i in range(nbins): + bin_col = np.loadtxt(filename, usecols=i + 1) + hdu_col = fits.Column(name="BIN%d" % (i + 1), format="D", array=bin_col) + cols.append(hdu_col) + + coldefs = fits.ColDefs(cols) + nz_hdu = fits.BinTableHDU.from_columns(coldefs, name="NZDATA") + + nz_lens_dict = { + "NZDATA": "T ", + "EXTNAME": "NZ_SOURCE", + "NBIN": nbins, + "NZ": len(z_low), + } + + for key, value in nz_lens_dict.items(): + nz_hdu.header[key] = value + + return nz_hdu + + +def load_pseudo_cl(cl_file): + """Load pseudo-C_ell spectra from .npy or FITS formats.""" + if cl_file.endswith(".npy"): + cl_block = np.load(cl_file) + if cl_block.shape[0] < 5: + raise ValueError( + f"Unexpected C_ell array shape {cl_block.shape} for {cl_file}" + ) + ell = np.asarray(cl_block[0], dtype=np.float64) + cl_ee = np.asarray(cl_block[1], dtype=np.float64) + cl_bb = np.asarray(cl_block[4], dtype=np.float64) + return ell, cl_ee, cl_bb + + if cl_file.endswith(".fits"): + with fits.open(cl_file) as hdul: + data = hdul[1].data + ell = np.asarray(data["ELL"], dtype=np.float64) + cl_ee = np.asarray(data["EE"], dtype=np.float64) + cl_bb = np.asarray(data["BB"], dtype=np.float64) + return ell, cl_ee, cl_bb + + raise NotImplementedError(f"Cl format not supported: {cl_file}") + + +def cl_to_fits(ell, cl_ee, cl_bb): + """Convert C_ell arrays to CosmoSIS FITS HDUs for EE and BB.""" + nbins = len(ell) + lst = np.arange(1, nbins + 1) + + col1 = fits.Column(name="BIN1", format="K", array=np.ones(nbins)) + col2 = fits.Column(name="BIN2", format="K", array=np.ones(nbins)) + col3 = fits.Column(name="ANGBIN", format="K", array=lst) + col5 = fits.Column(name="ANG", format="D", array=ell) + + col4_ee = fits.Column(name="VALUE", format="D", array=cl_ee) + coldefs_ee = fits.ColDefs([col1, col2, col3, col4_ee, col5]) + cl_ee_hdu = fits.BinTableHDU.from_columns(coldefs_ee, name="CELL_EE") + + cl_ee_dict = { + "2PTDATA": "T", + "QUANT1": "GEF", + "QUANT2": "GEF", + "KERNEL_1": "NZ_SOURCE", + "KERNEL_2": "NZ_SOURCE", + "WINDOWS": "SAMPLE", + } + for key, value in cl_ee_dict.items(): + cl_ee_hdu.header[key] = value + + col4_bb = fits.Column(name="VALUE", format="D", array=cl_bb) + coldefs_bb = fits.ColDefs([col1, col2, col3, col4_bb, col5]) + cl_bb_hdu = fits.BinTableHDU.from_columns(coldefs_bb, name="CELL_BB") + + cl_bb_dict = { + "2PTDATA": "T", + "QUANT1": "GBF", + "QUANT2": "GBF", + "KERNEL_1": "NZ_SOURCE", + "KERNEL_2": "NZ_SOURCE", + "WINDOWS": "SAMPLE", + } + for key, value in cl_bb_dict.items(): + cl_bb_hdu.header[key] = value + + return cl_ee_hdu, cl_bb_hdu + + +def cov_cl_to_fits(cov_file, cov_hdu="COVAR_FULL"): + """Convert pseudo-C_ell covariance to a CosmoSIS ImageHDU.""" + if cov_file.endswith(".fits"): + with fits.open(cov_file) as hdul: + cov_data = np.asarray(hdul[cov_hdu].data, dtype=np.float64) + elif cov_file.endswith(".npy"): + cov_data = np.load(cov_file) + else: + raise NotImplementedError( + f"Unsupported pseudo-Cl covariance format: {cov_file}" + ) + + if cov_data.shape[0] != cov_data.shape[1]: + raise ValueError("Pseudo-Cl covariance matrix must be square") + + cov_hdu = fits.ImageHDU(cov_data, name="COVMAT_CELL") + cov_dict = { + "COVDATA": "True", + "EXTNAME": "COVMAT_CELL", + "NAME_0": "CELL_EE", + "STRT_0": 0, + } + + for key, value in cov_dict.items(): + cov_hdu.header[key] = value + + return cov_hdu + + +def tau_to_fits(filename, theta=None): + """ + Convert tau statistics to FITS format. + If theta provided, override original values for consistency with xi. + """ + tau_stats = fits.getdata(filename) + + if theta is not None: + ang = theta + print("Using provided theta values for tau statistics (forcing consistency)") + else: + ang = tau_stats["theta"] + print("Using original tau theta values") + + tau_0_p_hdu = _create_2pt_hdu(tau_stats["tau_0_p"], ang, "TAU_0_PLUS", "G+R", "P+R") + tau_2_p_hdu = _create_2pt_hdu( + tau_stats["tau_2_p"], ang, "TAU_2_PLUS", "G+R", "SR+R" + ) + + return tau_0_p_hdu, tau_2_p_hdu + + +def rho_to_fits(filename, theta=None): + """ + Convert rho statistics to FITS format. + If theta provided, override original values for consistency with xi. + """ + rho_stat_hdul = fits.open(filename) + rho_stat_hdu = rho_stat_hdul[1].copy() + rho_stat_hdu.name = "RHO_STATS" + + if theta is not None: + print( + "Forcing rho statistics to use provided theta values (forcing consistency)" + ) + rho_stat_hdu.data = rho_stat_hdu.data.copy() + rho_stat_hdu.data["theta"] = theta + else: + print("Using original rho theta values") + + rho_stat_hdul.close() + return rho_stat_hdu + + +def covdat_to_fits( + filename_cov_xi, filename_cov_tau=None, filename_cov_cl=None, cov_hdu=None +): + """ + Convert CosmoCov covariance matrix to FITS format. + + If tau covariance provided, block with xi covariance. + """ + covmat_xi = np.loadtxt(filename_cov_xi) + + if filename_cov_tau is not None: + covmat_tau = np.load(filename_cov_tau) + nbins = int(len(covmat_tau) / 3) + covmat_tau = covmat_tau[: 2 * nbins, : 2 * nbins] + covmat = np.block( + [ + [covmat_xi, np.zeros((len(covmat_xi), len(covmat_tau)))], + [np.zeros((len(covmat_tau), len(covmat_xi))), covmat_tau], + ] + ) + else: + covmat = covmat_xi + + if filename_cov_cl is not None: + if filename_cov_cl.endswith(".fits"): + with fits.open(filename_cov_cl) as hdul: + cov_data = np.asarray(hdul[cov_hdu].data, dtype=np.float64) + elif filename_cov_cl.endswith(".npy"): + cov_data = np.load(filename_cov_cl) + else: + raise NotImplementedError( + f"Unsupported pseudo-Cl covariance format: {filename_cov_cl}" + ) + covmat = np.block( + [ + [covmat, np.zeros((len(covmat), len(cov_data)))], + [np.zeros((len(cov_data), len(covmat))), cov_data], + ] + ) + + if len(covmat) != len(covmat[0]): + raise RuntimeError("Covariance matrix is not square") + + cov_hdu = fits.ImageHDU(covmat) + + cov_dict = { + "COVDATA": "True", + "EXTNAME": "COVMAT", + "NAME_0": "XI_PLUS", + "STRT_0": 0, + "NAME_1": "XI_MINUS", + "STRT_1": int(len(covmat_xi) / 2), + } + + if filename_cov_tau: + cov_dict.update( + { + "NAME_2": "TAU_0_PLUS", + "STRT_2": len(covmat_xi), + "NAME_3": "TAU_2_PLUS", + "STRT_3": len(covmat_xi) + int(len(covmat_tau) / 2), + } + ) + + filename_cov_cl and cov_dict.update( + { + "NAME_4": "CELL_EE", + "STRT_4": len(covmat_xi) + (len(covmat_tau) if filename_cov_tau else 0), + } + ) + + for key, value in cov_dict.items(): + cov_hdu.header[key] = value + + return cov_hdu diff --git a/src/sp_validation/tests/test_column_schema.py b/src/sp_validation/tests/test_column_schema.py index 6181c886e..1f4527990 100644 --- a/src/sp_validation/tests/test_column_schema.py +++ b/src/sp_validation/tests/test_column_schema.py @@ -39,6 +39,10 @@ "R", "R11", "R22", + "CSL_DIR", + "DEFAULT", + "FITS_FILE", + "SCRATCH", # FITS header keywords "EXTNAME", "R_S11", @@ -62,6 +66,11 @@ "ELL", "NAME_{}", "STRT_{}", + # CosmoSIS/twopoint FITS products + "ANG", + "COVMAT", + "XI_MINUS", + "XI_PLUS", # in-memory tables sp_validation builds itself "C11", "C22", diff --git a/src/sp_validation/tests/test_cosmosis_fitting.py b/src/sp_validation/tests/test_cosmosis_fitting.py deleted file mode 100644 index b9625c746..000000000 --- a/src/sp_validation/tests/test_cosmosis_fitting.py +++ /dev/null @@ -1,529 +0,0 @@ -"""Integration tests for the CosmoSIS data-FITS assembly. - -This pins the riskiest, least-tested step in the inference data-product path: -``cosmo_inference/scripts/cosmosis_fitting.py`` turning UNIONS validation -outputs (xi+/-, n(z), a blocked covariance, and -- with ``--use-rho-tau`` -- -rho/tau statistics and their covariance) into the single CosmoSIS data FITS that -a chain consumes. We stop short of the chain itself: no cosmosis, no MCMC. - -Two danger surfaces drive the design: - -1. **The HDU set + ordering.** The assembled ``HDUList`` is built by hand in the - script's ``__main__`` block (not a reusable function), so the only way it can - silently change shape is for nobody to be watching it. The CLI tests below - watch it: they assert the exact extension-name set and the data-vector - ordering (xi+ then xi-), for both the plain-xi product and the rho/tau one. - -2. **The blocked covariance offsets.** ``covdat_to_fits`` lays the xi+ block, - the xi- block, and -- under ``--use-rho-tau`` -- the tau_0/tau_2 blocks onto - a single matrix, recording each block's start index in the FITS header - (``STRT_0..STRT_3``). The tau covariance is *truncated* from 3*nbins to - 2*nbins before blocking; a wrong offset or a wrong truncation lands the data - on the wrong rows and is invisible until a chain quietly mis-fits. These are - pinned at the function level with deterministic synthetic inputs whose offsets - we work out by hand, plus teeth that flip red on a perturbed input or a - wrong-offset expectation. - -The script imports only numpy + astropy.io.fits (no cosmosis), so the whole -assembly path is exercisable in the production container. The end-to-end CLI run -also triggers ``generate_cosmosis_config``; we feed it minimal synthetic template -INIs in ``tmp_path`` so the run completes without touching the repo's real -templates or any cluster data. - -:Author: cdaley - -""" - -import importlib.util -import subprocess -import sys -from pathlib import Path - -import numpy as np -import pytest -from astropy.io import fits - -# --- load cosmosis_fitting.py as a module (it is a script, not a package) ---- - -_SCRIPT = ( - Path(__file__).resolve().parents[3] - / "cosmo_inference" - / "scripts" - / "cosmosis_fitting.py" -) - - -def _load_module(): - """Import cosmosis_fitting.py by path; skip cleanly if a dep is missing.""" - if not _SCRIPT.exists(): - pytest.skip(f"cosmosis_fitting.py not found at {_SCRIPT}") - spec = importlib.util.spec_from_file_location("cosmosis_fitting", _SCRIPT) - module = importlib.util.module_from_spec(spec) - try: - spec.loader.exec_module(module) - except ImportError as exc: # pragma: no cover - container has numpy/astropy - pytest.importorskip(getattr(exc, "name", "") or "cosmosis_fitting_dependency") - raise - return module - - -cf = _load_module() - - -# --- deterministic synthetic dimensions -------------------------------------- -# -# N_ANG angular bins. The xi data vector is xi+ (N_ANG) then xi- (N_ANG); the xi -# covariance is therefore (2*N_ANG) square. The tau covariance is stored as -# 3 statistics x N_ANG bins = (3*N_ANG) square, and covdat_to_fits keeps only -# the first 2*N_ANG rows/cols (tau_0, tau_2) before blocking it onto the xi -# covariance. With N_ANG=3 every offset is small enough to verify by eye. - -N_ANG = 3 - - -def _xi_cov_matrix(): - """A (2*N_ANG) covariance with unique, recognizable entries.""" - n = 2 * N_ANG - return np.arange(1, n * n + 1, dtype=float).reshape(n, n) - - -def _tau_cov_full(): - """A (3*N_ANG) tau covariance, offset by 1000 so its block is identifiable.""" - n = 3 * N_ANG - return np.arange(1, n * n + 1, dtype=float).reshape(n, n) + 1000.0 - - -def _write_xi_cov(path): - np.savetxt(path, _xi_cov_matrix()) - return path - - -def _write_tau_cov(path): - np.save(path, _tau_cov_full()) - return path - - -# ============================================================================= -# Function-level: the blocked covariance (covdat_to_fits) -- the danger -# ============================================================================= - - -def test_cov_xi_only_offsets(tmp_path): - """xi-only covariance: shape (2*N_ANG), XI_MINUS block starts at N_ANG. - - With only xi the covariance is passed through verbatim and the header marks - the xi+ block at row 0 and the xi- block at row N_ANG (= half the matrix). - """ - cov_file = _write_xi_cov(tmp_path / "cov_xi.txt") - hdu = cf.covdat_to_fits(str(cov_file), filename_cov_tau=None) - - assert hdu.data.shape == (2 * N_ANG, 2 * N_ANG) - assert hdu.header["EXTNAME"] == "COVMAT" - assert hdu.header["NAME_0"] == "XI_PLUS" - assert hdu.header["STRT_0"] == 0 - assert hdu.header["NAME_1"] == "XI_MINUS" - assert hdu.header["STRT_1"] == N_ANG - # No tau augmentation when filename_cov_tau is None. - assert "NAME_2" not in hdu.header - # The matrix is the input, unmodified. - assert np.array_equal(hdu.data, _xi_cov_matrix()) - - -def test_cov_with_tau_block_layout(tmp_path): - """xi+tau covariance: the four blocks land at the offsets we predict by hand. - - Layout for N_ANG=3: - rows 0..2 -> xi+ (STRT_0 = 0) - rows 3..5 -> xi- (STRT_1 = 3) - rows 6..8 -> tau_0 (STRT_2 = 6 = len(covmat_xi)) - rows 9..11 -> tau_2 (STRT_3 = 9 = len(covmat_xi) + len(covmat_tau)/2) - The tau covariance is truncated from 3*N_ANG to 2*N_ANG before blocking, and - the cross blocks between xi and tau are zero. - """ - cov_xi = _write_xi_cov(tmp_path / "cov_xi.txt") - cov_tau = _write_tau_cov(tmp_path / "cov_tau.npy") - - hdu = cf.covdat_to_fits(str(cov_xi), filename_cov_tau=str(cov_tau)) - - # Combined size: xi (2*N_ANG) + truncated tau (2*N_ANG). - assert hdu.data.shape == (4 * N_ANG, 4 * N_ANG) - - # Block-start offsets, pinned to concrete integers. - assert hdu.header["STRT_0"] == 0 - assert hdu.header["STRT_1"] == N_ANG - assert hdu.header["STRT_2"] == 2 * N_ANG - assert hdu.header["STRT_3"] == 3 * N_ANG - assert hdu.header["NAME_0"] == "XI_PLUS" - assert hdu.header["NAME_1"] == "XI_MINUS" - assert hdu.header["NAME_2"] == "TAU_0_PLUS" - assert hdu.header["NAME_3"] == "TAU_2_PLUS" - - cov = hdu.data - xi = _xi_cov_matrix() - tau_trunc = _tau_cov_full()[: 2 * N_ANG, : 2 * N_ANG] - - # xi block sits top-left; tau block sits bottom-right (truncated). - assert np.array_equal(cov[: 2 * N_ANG, : 2 * N_ANG], xi) - assert np.array_equal(cov[2 * N_ANG :, 2 * N_ANG :], tau_trunc) - # The tau block is the *truncated* matrix: its [0,0] is the full matrix's - # [0,0], and the dropped tau_5 rows (>= 2*N_ANG) never appear. - assert cov[2 * N_ANG, 2 * N_ANG] == _tau_cov_full()[0, 0] - assert _tau_cov_full()[2 * N_ANG, 2 * N_ANG] not in cov - - # Cross blocks between xi and tau are exactly zero. - assert np.all(cov[: 2 * N_ANG, 2 * N_ANG :] == 0.0) - assert np.all(cov[2 * N_ANG :, : 2 * N_ANG] == 0.0) - - -def test_cov_with_tau_offsets_have_teeth(tmp_path): - """Teeth: the tau block does NOT live at the xi- offset. - - A common wrong expectation is that the tau block starts where xi- starts - (N_ANG) rather than after all of xi (2*N_ANG). Pin that the truncated tau - block is genuinely at 2*N_ANG, so an off-by-a-block error would be caught. - """ - cov_xi = _write_xi_cov(tmp_path / "cov_xi.txt") - cov_tau = _write_tau_cov(tmp_path / "cov_tau.npy") - hdu = cf.covdat_to_fits(str(cov_xi), filename_cov_tau=str(cov_tau)) - - tau_trunc = _tau_cov_full()[: 2 * N_ANG, : 2 * N_ANG] - # Correct offset reproduces the tau block... - assert np.array_equal(hdu.data[2 * N_ANG :, 2 * N_ANG :], tau_trunc) - # ...and a wrong (xi-) offset does not. - wrong = hdu.data[N_ANG : N_ANG + 2 * N_ANG, N_ANG : N_ANG + 2 * N_ANG] - assert not np.array_equal(wrong, tau_trunc) - - -def test_cov_non_square_raises(tmp_path): - """A non-square xi covariance must fail loudly, not assemble silently.""" - bad = tmp_path / "cov_bad.txt" - np.savetxt(bad, np.arange(6, dtype=float).reshape(2, 3)) - with pytest.raises(Exception): - cf.covdat_to_fits(str(bad), filename_cov_tau=None) - - -# ============================================================================= -# Function-level: the 2pt and n(z) HDUs feeding the data vector -# ============================================================================= - - -def test_2pt_hdu_columns_and_header(): - """A 2pt HDU carries VALUE/ANG and the CosmoSIS 2PTDATA header marker.""" - values = np.array([1.0, 2.0, 3.0]) - theta = np.array([5.0, 50.0, 500.0]) - hdu = cf._create_2pt_hdu(values, theta, "XI_PLUS", "G+R", "G+R") - - assert hdu.name == "XI_PLUS" - assert set(["BIN1", "BIN2", "ANGBIN", "VALUE", "ANG"]).issubset(set(hdu.data.names)) - assert np.array_equal(hdu.data["VALUE"], values) - assert np.array_equal(hdu.data["ANG"], theta) - assert np.array_equal(hdu.data["ANGBIN"], np.arange(1, len(values) + 1)) - assert hdu.header["2PTDATA"] == "T" - assert hdu.header["QUANT1"] == "G+R" - - -def test_nz_to_fits(tmp_path): - """n(z) text -> NZDATA HDU named NZ_SOURCE with correct bin edges/count.""" - z_mid = np.array([0.1, 0.3, 0.5]) - bin1 = np.array([0.2, 0.5, 0.3]) - nz_path = tmp_path / "nz.txt" - np.savetxt(nz_path, np.column_stack([z_mid, bin1])) - - hdu = cf.nz_to_fits(str(nz_path)) - - assert hdu.name == "NZ_SOURCE" - assert hdu.header["NZDATA"].strip() == "T" - assert hdu.header["NBIN"] == 1 - assert hdu.header["NZ"] == len(z_mid) - assert np.array_equal(hdu.data["Z_MID"], z_mid) - # Edges are midpoints +/- half the (uniform) step. - step = z_mid[1] - z_mid[0] - assert np.allclose(hdu.data["Z_LOW"], z_mid - step / 2) - assert np.allclose(hdu.data["Z_HIGH"], z_mid + step / 2) - assert np.array_equal(hdu.data["BIN1"], bin1) - - -# ============================================================================= -# End-to-end CLI: the assembled HDU set + ordering (the __main__ assembly) -# ============================================================================= -# -# The full HDUList is built in __main__, reachable only through the CLI. We run -# the script as a subprocess against synthetic inputs in tmp_path -- including -# minimal template INIs so generate_cosmosis_config completes -- and inspect the -# written FITS. Nothing here touches the repo's real templates or cluster data. - -# Minimal template containing exactly the section headers _generate_ini_file's -# regexes rewrite. Anything else would just be passed through untouched. -_TEMPLATE_INI = "\n".join( - [ - "[DEFAULT]", - "[runtime]", - "[output]", - "[pipeline]", - "[2pt_like]", - "covmat_name=PLACEHOLDER", - "[test]", - "[polychord]", - "", - ] -) - -# Theta grid shared by xi (and forced onto rho/tau for consistency). -_THETA = np.array([5.0, 50.0, 500.0]) -_XIP = np.array([10.0, 11.0, 12.0]) # xi+ data vector -_XIM = np.array([20.0, 21.0, 22.0]) # xi- data vector - - -def _write_treecorr_xi(path, values, extname): - """A TreeCorr-style xi FITS: data in HDU 1 with ANG + VALUE columns. - - Data mode's ``treecorr_to_fits`` returns ``hdul[1]`` verbatim, so the - extension name carried by the real TreeCorr file is what lands in the - assembled product. CosmoSIS looks the blocks up by name - (``data_sets=XI_PLUS XI_MINUS``), so the inputs must carry those EXTNAMEs; - we set them here to be faithful to the real files. - """ - cols = fits.ColDefs( - [ - fits.Column(name="ANG", format="D", unit="arcmin", array=_THETA), - fits.Column(name="VALUE", format="D", array=values), - fits.Column(name="meanr", format="D", array=_THETA), - fits.Column(name="xip", format="D", array=values), - fits.Column(name="xim", format="D", array=values), - ] - ) - hdul = fits.HDUList( - [fits.PrimaryHDU(), fits.BinTableHDU.from_columns(cols, name=extname)] - ) - hdul.writeto(path, overwrite=True) - return path - - -def _write_rho_stats(path): - """rho-stats FITS: HDU 1 carries a theta column (values irrelevant here).""" - cols = fits.ColDefs( - [ - fits.Column(name="theta", format="D", array=_THETA), - fits.Column(name="rho_0", format="D", array=np.ones(N_ANG)), - ] - ) - hdul = fits.HDUList([fits.PrimaryHDU(), fits.BinTableHDU.from_columns(cols)]) - hdul.writeto(path, overwrite=True) - return path - - -def _write_tau_stats(path): - """tau-stats FITS: theta + tau_0_p/tau_2_p columns (read by tau_to_fits).""" - cols = fits.ColDefs( - [ - fits.Column(name="theta", format="D", array=_THETA), - fits.Column(name="tau_0_p", format="D", array=np.array([1.0, 2.0, 3.0])), - fits.Column(name="tau_2_p", format="D", array=np.array([4.0, 5.0, 6.0])), - ] - ) - hdul = fits.HDUList([fits.PrimaryHDU(), fits.BinTableHDU.from_columns(cols)]) - hdul.writeto(path, overwrite=True) - return path - - -def _write_nz(path): - z_mid = np.array([0.1, 0.3, 0.5]) - np.savetxt(path, np.column_stack([z_mid, np.array([0.2, 0.5, 0.3])])) - return path - - -def _run_cli(tmp_path, *, use_rho_tau): - """Run cosmosis_fitting.py against synthetic tmp_path inputs; return FITS path. - - Builds a self-contained output root (with a cosmosis_config holding the - minimal template + value/prior INIs config generation expects) so the run - completes end-to-end. Returns the path to the assembled data FITS. - """ - out_root = tmp_path / "out_root" - template_dir = out_root / "cosmosis_config" - template_dir.mkdir(parents=True) - - # Templates + the value/prior INIs generate_cosmosis_config reads/writes. - for name in ( - "cosmosis_pipeline_A_ia.ini", - "cosmosis_pipeline_A_psf.ini", - ): - (template_dir / name).write_text(_TEMPLATE_INI) - for name in ("values_ia.ini", "values_psf.ini", "priors.ini", "priors_psf.ini"): - (template_dir / name).write_text("[parameters]\n") - - nz = _write_nz(tmp_path / "nz.txt") - xip = _write_treecorr_xi(tmp_path / "xi_plus.fits", _XIP, "XI_PLUS") - xim = _write_treecorr_xi(tmp_path / "xi_minus.fits", _XIM, "XI_MINUS") - cov_xi = _write_xi_cov(tmp_path / "cov_xi.txt") - - root = "TESTROOT" - cmd = [ - sys.executable, - str(_SCRIPT), - "--cosmosis-root", - root, - "--data-dir", - str(tmp_path / "chains"), - "--nz-file", - str(nz), - "--output-root", - str(out_root), - "--template-dir", - str(template_dir), - "--xi", - str(xip), - str(xim), - "--cov-xi", - str(cov_xi), - ] - if use_rho_tau: - cmd += [ - "--use-rho-tau", - "--rho-stats", - str(_write_rho_stats(tmp_path / "rho.fits")), - "--tau-stats", - str(_write_tau_stats(tmp_path / "tau.fits")), - "--cov-tau", - str(_write_tau_cov(tmp_path / "cov_tau.npy")), - ] - - proc = subprocess.run(cmd, capture_output=True, text=True) - assert proc.returncode == 0, ( - f"cosmosis_fitting.py failed:\nSTDOUT:\n{proc.stdout}\nSTDERR:\n{proc.stderr}" - ) - - out_fits = out_root / "data" / root / f"cosmosis_{root}.fits" - assert out_fits.exists(), f"expected FITS not written: {out_fits}" - return out_fits - - -def test_cli_plain_xi_hdu_set_and_ordering(tmp_path): - """Plain-xi product: exact HDU set, xi+ -> xi- data ordering, cov shape.""" - out_fits = _run_cli(tmp_path, use_rho_tau=False) - - with fits.open(out_fits) as hdul: - names = [h.name for h in hdul] - # Order from __main__: PRIMARY, NZ_SOURCE, COVMAT, XI_PLUS, XI_MINUS. - assert names == ["PRIMARY", "NZ_SOURCE", "COVMAT", "XI_PLUS", "XI_MINUS"] - - xip = hdul["XI_PLUS"].data["VALUE"] - xim = hdul["XI_MINUS"].data["VALUE"] - assert np.array_equal(xip, _XIP) - assert np.array_equal(xim, _XIM) - - # The data vector is xi+ THEN xi- -- ordering is load-bearing for the - # covariance blocks below. - data_vector = np.concatenate([xip, xim]) - assert np.array_equal(data_vector, np.concatenate([_XIP, _XIM])) - assert len(data_vector) == 2 * N_ANG - - cov = hdul["COVMAT"].data - assert cov.shape == (2 * N_ANG, 2 * N_ANG) - assert hdul["COVMAT"].header["STRT_1"] == N_ANG - - -def test_cli_rho_tau_hdu_set_and_cov_block(tmp_path): - """--use-rho-tau product: rho/tau HDUs appended, covariance augmented. - - The HDU set gains TAU_0_PLUS, TAU_2_PLUS, RHO_STATS, and the COVMAT grows to - (4*N_ANG) with the tau blocks at offsets 2*N_ANG and 3*N_ANG. - """ - out_fits = _run_cli(tmp_path, use_rho_tau=True) - - with fits.open(out_fits) as hdul: - names = [h.name for h in hdul] - # Order from __main__: PRIMARY, NZ_SOURCE, COVMAT, XI_PLUS, XI_MINUS, - # then TAU_0_PLUS, TAU_2_PLUS, RHO_STATS. - assert names == [ - "PRIMARY", - "NZ_SOURCE", - "COVMAT", - "XI_PLUS", - "XI_MINUS", - "TAU_0_PLUS", - "TAU_2_PLUS", - "RHO_STATS", - ] - - # xi data vector unchanged and still xi+ then xi-. - assert np.array_equal(hdul["XI_PLUS"].data["VALUE"], _XIP) - assert np.array_equal(hdul["XI_MINUS"].data["VALUE"], _XIM) - - # Covariance now carries the augmented block layout. - cov_hdu = hdul["COVMAT"] - assert cov_hdu.data.shape == (4 * N_ANG, 4 * N_ANG) - assert cov_hdu.header["STRT_0"] == 0 - assert cov_hdu.header["STRT_1"] == N_ANG - assert cov_hdu.header["STRT_2"] == 2 * N_ANG - assert cov_hdu.header["STRT_3"] == 3 * N_ANG - assert cov_hdu.header["NAME_2"] == "TAU_0_PLUS" - assert cov_hdu.header["NAME_3"] == "TAU_2_PLUS" - - # The tau block is the truncated tau covariance at the bottom-right. - tau_trunc = _tau_cov_full()[: 2 * N_ANG, : 2 * N_ANG] - assert np.array_equal(cov_hdu.data[2 * N_ANG :, 2 * N_ANG :], tau_trunc) - - # rho/tau theta is forced onto the xi theta grid (consistency step). - assert np.allclose(hdul["RHO_STATS"].data["theta"], _THETA) - assert np.allclose(hdul["TAU_0_PLUS"].data["ANG"], _THETA) - - -def test_cli_perturbed_xi_changes_data_vector(tmp_path): - """Teeth: a changed xi+ input must change the written data vector. - - If the assembly silently dropped or reordered the xi+ values, this would not - move -- so a real input perturbation flowing through to the FITS is the - end-to-end guarantee. - """ - baseline = _run_cli(tmp_path, use_rho_tau=False) - with fits.open(baseline) as hdul: - base_xip = hdul["XI_PLUS"].data["VALUE"].copy() - - # Perturb xi+ and rerun in a fresh output area. - tmp2 = tmp_path / "perturbed" - tmp2.mkdir() - perturbed_xip = _XIP + 100.0 - # Reuse the runner but with a different xi+ vector by monkey-free rewrite. - out_root = tmp2 / "out_root" - template_dir = out_root / "cosmosis_config" - template_dir.mkdir(parents=True) - (template_dir / "cosmosis_pipeline_A_ia.ini").write_text(_TEMPLATE_INI) - (template_dir / "values_ia.ini").write_text("[parameters]\n") - (template_dir / "priors.ini").write_text("[parameters]\n") - - nz = _write_nz(tmp2 / "nz.txt") - xip = _write_treecorr_xi(tmp2 / "xi_plus.fits", perturbed_xip, "XI_PLUS") - xim = _write_treecorr_xi(tmp2 / "xi_minus.fits", _XIM, "XI_MINUS") - cov_xi = _write_xi_cov(tmp2 / "cov_xi.txt") - root = "TESTROOT2" - proc = subprocess.run( - [ - sys.executable, - str(_SCRIPT), - "--cosmosis-root", - root, - "--data-dir", - str(tmp2 / "chains"), - "--nz-file", - str(nz), - "--output-root", - str(out_root), - "--template-dir", - str(template_dir), - "--xi", - str(xip), - str(xim), - "--cov-xi", - str(cov_xi), - ], - capture_output=True, - text=True, - ) - assert proc.returncode == 0, proc.stderr - out_fits = out_root / "data" / root / f"cosmosis_{root}.fits" - with fits.open(out_fits) as hdul: - new_xip = hdul["XI_PLUS"].data["VALUE"] - - assert not np.array_equal(base_xip, new_xip) - assert np.array_equal(new_xip, perturbed_xip) diff --git a/src/sp_validation/tests/test_inference.py b/src/sp_validation/tests/test_inference.py new file mode 100644 index 000000000..428bae295 --- /dev/null +++ b/src/sp_validation/tests/test_inference.py @@ -0,0 +1,101 @@ +"""Synthetic tests of the SACC-to-CosmoSIS workflow boundary.""" + +import configparser +import json +from pathlib import Path + +import numpy as np +from astropy.io import fits + +from sp_validation import sacc_io as sio +from sp_validation.inference import prepare, render, summarize + + +def test_prepare_uses_reference_covariance_not_reference_signal(tmp_path): + z = np.linspace(0.01, 2.01, 20) + s = sio.new_sacc([(z, np.ones(20))], metadata={"blind": "mock"}) + theta = np.array([2.0, 4.0, 8.0]) + sio.add_xi(s, (0, 0), theta, theta, -theta, grid="reporting", theta_nom=theta) + s.add_covariance(np.eye(6)) + original = tmp_path / "mock.sacc" + s.save_fits(original) + ref = s.copy() + # Reference means and measured angles differ, nominal bins identify rows. + for p in ref.data: + p.value = 999.0 + p.tags["theta"] *= 1.001 + ref.add_covariance(np.eye(6) * 7, overwrite=True) + reference = tmp_path / "covariance.sacc" + ref.save_fits(reference) + out = tmp_path / "out.fits" + provenance = prepare(original, out, covariance=reference) + assert provenance["blind"] == "mock" + with fits.open(out) as hdus: + np.testing.assert_array_equal(hdus["XI_PLUS"].data["VALUE"], theta) + np.testing.assert_array_equal(hdus["XI_PLUS"].data["ANG"], theta) + np.testing.assert_array_equal(hdus["COVMAT"].data, np.eye(6) * 7) + np.testing.assert_array_equal(hdus["NZ_SOURCE"].data["BIN1"], np.ones(20)) + + +def test_render_explicit_inputs_and_all_bin_cuts(tmp_path): + template = tmp_path / "template.ini" + template.write_text( + "[pipeline]\nvalues = old\n[runtime]\nsampler = polychord\n" + "[output]\nformat = text\n[2pt_like]\n" + "angle_range_XI_PLUS_1_1 = 10.0 200.0\n" + ) + out = tmp_path / "pipeline.ini" + render( + template, + out, + fits_path=tmp_path / "data.fits", + values=tmp_path / "v.ini", + priors=tmp_path / "p.ini", + chain=tmp_path / "chain.txt", + n_bins=4, + settings={"n_live": 300}, + scale_cuts={"XI_PLUS": [12, 83]}, + ) + ini = configparser.ConfigParser(interpolation=None) + ini.read(out) + assert ini["pipeline"]["values"] == str(tmp_path / "v.ini") + assert ini["pipeline"]["priors"] == str(tmp_path / "p.ini") + assert ini["runtime"]["sampler"] == "nautilus" + assert ini["nautilus"]["n_live"] == "300" + assert ini["2pt_like"]["angle_range_XI_PLUS_3_4"] == "12.0 83.0" + assert len([k for k in ini["2pt_like"] if k.startswith("angle_range")]) == 10 + + +def test_psf_template_uses_likelihood_that_consumes_xi_sys(tmp_path): + root = Path(__file__).resolve().parents[3] + template = ( + root / "cosmo_inference/cosmosis_config/templates/cosmosis_pipeline_A_psf.ini" + ) + out = tmp_path / "pipeline.ini" + render( + template, + out, + fits_path=tmp_path / "data.fits", + values=tmp_path / "v.ini", + priors=tmp_path / "p.ini", + chain=tmp_path / "chain.txt", + ) + ini = configparser.ConfigParser() + ini.read(out) + assert ini["2pt_like"]["file"].endswith("/likelihood/2pt/2pt_like_xi_sys.py") + assert ini["2pt_like"]["add_xi_sys"] == "T" + + +def test_weighted_summary(tmp_path): + chain = tmp_path / "chain.txt" + chain.write_text( + "#cosmological_parameters--s_8 log_weight prior post\n" + "0.7 0 0 0\n0.9 1.0986122886681098 0 0\n" + ) + out = tmp_path / "summary.json" + summarize(chain, out) + result = json.loads(out.read_text()) + np.testing.assert_allclose( + result["parameters"]["cosmological_parameters--s_8"]["mean"], 0.85 + ) + np.testing.assert_allclose(result["effective_samples"], 1.6) diff --git a/src/sp_validation/tests/test_sacc_io_realdata.py b/src/sp_validation/tests/test_sacc_io_realdata.py index abe1f230b..e839cdbc7 100644 --- a/src/sp_validation/tests/test_sacc_io_realdata.py +++ b/src/sp_validation/tests/test_sacc_io_realdata.py @@ -38,12 +38,7 @@ "glass_mock_00001": _DATA / "glass_mock_00001" / "cosmosis_glass_mock_00001.fits", } -_SCRIPT = ( - Path(__file__).resolve().parents[3] - / "cosmo_inference" - / "scripts" - / "cosmosis_fitting.py" -) +_SCRIPT = Path(__file__).with_name("legacy_twopoint_reference.py") def _load_cf(): diff --git a/src/sp_validation/tests/test_sacc_io_twopoint.py b/src/sp_validation/tests/test_sacc_io_twopoint.py index 78f6809da..ed755a39a 100644 --- a/src/sp_validation/tests/test_sacc_io_twopoint.py +++ b/src/sp_validation/tests/test_sacc_io_twopoint.py @@ -1,7 +1,7 @@ """Byte-compare tests for the SACC -> 2pt-FITS converter. The converter (:mod:`sp_validation.sacc_io`) must reproduce the CosmoSIS -2pt-FITS that ``cosmo_inference/scripts/cosmosis_fitting.py`` assembles today, +2pt-FITS that the retired ``cosmosis_fitting.py`` assembled, so the inference chain (``2pt_like`` and Sacha Guerrini's rho/tau ``2pt_like_xi_sys`` fork) runs untouched behind it. The strongest possible check is *byte* equality, and astropy writes FITS deterministically, so that is what we @@ -24,9 +24,8 @@ converter copies them from the sidecar exactly as today's assembly does. A teeth test pins that a perturbed input moves the output. -The reference builder imports ``cosmosis_fitting.py`` by path (it is a script, -not a package module), skipping cleanly if a dependency is missing -- the same -loader pattern as ``test_cosmosis_fitting.py``. +The reference builder imports the frozen legacy builders from +``legacy_twopoint_reference.py``; no obsolete preparation CLI is retained. """ import importlib.util @@ -38,18 +37,13 @@ from sp_validation import sacc_io -_SCRIPT = ( - Path(__file__).resolve().parents[3] - / "cosmo_inference" - / "scripts" - / "cosmosis_fitting.py" -) +_SCRIPT = Path(__file__).with_name("legacy_twopoint_reference.py") def _load_cf(): - """Import cosmosis_fitting.py by path; skip cleanly if a dep is missing.""" + """Import the frozen byte-compatibility oracle by path.""" if not _SCRIPT.exists(): - pytest.skip(f"cosmosis_fitting.py not found at {_SCRIPT}") + pytest.fail(f"legacy byte-compatibility oracle not found at {_SCRIPT}") spec = importlib.util.spec_from_file_location("cosmosis_fitting", _SCRIPT) module = importlib.util.module_from_spec(spec) try: @@ -475,23 +469,15 @@ def test_rho_tau_sidecars_required_together(tmp_path): # ============================================================================= -def test_tomographic_sacc_raises(tmp_path): - """A multi-bin SACC fails fast instead of silently truncating to (0, 0). - - Review finding (HIGH): ``n_bins`` alone drove the NZDATA column count while - the data vector and covariance were read from bin ``(0, 0)`` only, so a - 2-bin SACC + ``n_bins=2`` emitted a plausible-looking FITS carrying 1/3 of - the data. Both the ``n_bins`` and the tracer-pair mismatch must raise. - """ +def test_wrong_bin_count_raises(tmp_path): + """A wrong n_bins must not silently truncate a tomographic product.""" inp = _inputs(seed=30) s = sacc_io.new_sacc({0: (inp["z"], inp["nz"]), 1: (inp["z"], inp["nz"])}) for pair in [(0, 0), (0, 1), (1, 1)]: sacc_io.add_xi(s, pair, inp["theta"], inp["xip"], inp["xim"], grid="reporting") s.add_covariance(np.eye(len(s.mean))) - with pytest.raises(ValueError, match="single-bin only"): - sacc_io.sacc_to_twopoint_fits(s, str(tmp_path / "x.fits"), n_bins=2) - with pytest.raises(ValueError, match="single-bin only"): + with pytest.raises(ValueError, match="does not match"): sacc_io.sacc_to_twopoint_fits(s, str(tmp_path / "x.fits"), n_bins=1) assert not (tmp_path / "x.fits").exists() diff --git a/src/sp_validation/tests/test_twopoint_tomography.py b/src/sp_validation/tests/test_twopoint_tomography.py new file mode 100644 index 000000000..e8cc6bf53 --- /dev/null +++ b/src/sp_validation/tests/test_twopoint_tomography.py @@ -0,0 +1,91 @@ +"""Round-trip through the same twopoint reader used by CSL's 2pt_like.""" + +import importlib.util +import os +from pathlib import Path + +import numpy as np +import pytest + +from sp_validation import sacc_io as sio + + +def _reader(): + path = Path(os.environ.get("CSL_DIR", "/opt/cosmosis-standard-library")) + path /= "likelihood/2pt/twopoint.py" + if not path.exists(): + pytest.skip("CSL twopoint reader is not installed") + spec = importlib.util.spec_from_file_location("csl_twopoint", path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +@pytest.mark.parametrize("n_bins", [1, 4]) +def test_twopoint_roundtrip(tmp_path, n_bins): + z = np.linspace(0.01, 2.01, 101) + nz = {i: (z, np.exp(-(((z - 0.3 * (i + 1)) / 0.2) ** 2))) for i in range(n_bins)} + s = sio.new_sacc(nz) + pairs = [(i, j) for i in range(n_bins) for j in range(i, n_bins)] + theta, ell = np.geomspace(1, 250, 5), np.array([100.0, 300.0, 900.0]) + # Reverse insertion order so a converter assuming canonical SACC order fails. + for k, pair in enumerate(reversed(pairs)): + sio.add_xi(s, pair, theta, k + theta, -k - theta, grid="reporting") + sio.add_pseudo_cl( + s, pair, ell, k + ell, window_ells=ell, window_weights=np.eye(len(ell)) + ) + sio.add_xi(s, pair, theta * 0.01, theta * 99, theta * 99, grid="integration") + n = len(s.mean) + # Asymmetric row/column encoding detects every wrong permutation. + cov = np.arange(n * n, dtype=float).reshape(n, n) + s.add_covariance(cov) + path = tmp_path / "two.fits" + sio.sacc_to_twopoint_fits(s, path, n_bins=n_bins) + reader = _reader() + # The legacy file has separate ξ and Cℓ covariances. Read spectra without + # assigning errors, and read each covariance with CSL's own covariance reader. + two = reader.TwoPointFile.from_fits(path, covmat_name=None) + harmonic = two + from astropy.io import fits + + with fits.open(path) as hdus: + xi_cov = reader.CovarianceMatrixInfo.from_fits(hdus["COVMAT"]) + cl_cov = reader.CovarianceMatrixInfo.from_fits(hdus["COVMAT_CELL"]) + + def rows(dtype): + return np.array( + [ + k + for pair in pairs + for k, p in enumerate(s.data) + if p.data_type == dtype + and p.tracers == sio._pair(pair) + and p.tags["grid"] == "reporting" + ] + ) + + for name, dtype, product, angle in [ + ("XI_PLUS", sio.XI_PLUS, two, theta), + ("XI_MINUS", sio.XI_MINUS, two, theta), + ("CELL_EE", sio.CL_EE, harmonic, ell), + ]: + spectrum = product.get_spectrum(name) + np.testing.assert_array_equal(spectrum.value, s.mean[rows(dtype)]) + np.testing.assert_array_equal( + spectrum.bin1, np.repeat([i + 1 for i, j in pairs], len(angle)) + ) + np.testing.assert_array_equal( + spectrum.bin2, np.repeat([j + 1 for i, j in pairs], len(angle)) + ) + np.testing.assert_array_equal( + spectrum.angular_bin, np.tile(np.arange(1, len(angle) + 1), len(pairs)) + ) + np.testing.assert_array_equal(spectrum.angle, np.tile(angle, len(pairs))) + xi_rows = np.concatenate([rows(sio.XI_PLUS), rows(sio.XI_MINUS)]) + np.testing.assert_array_equal(xi_cov.covmat, cov[np.ix_(xi_rows, xi_rows)]) + cl_rows = rows(sio.CL_EE) + np.testing.assert_array_equal(cl_cov.covmat, cov[np.ix_(cl_rows, cl_rows)]) + kernel = two.get_kernel("NZ_SOURCE") + np.testing.assert_array_equal(kernel.z, z) + for i in range(n_bins): + np.testing.assert_array_equal(kernel.nzs[i], nz[i][1]) diff --git a/workflow/README.md b/workflow/README.md index fa84a7b4b..b5b5d23b4 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -122,15 +122,16 @@ of it: expect to edit the default profile's `--bind` list for your machine. Set `glass_mocks.data_dir` in the paper config to the suite's `results` directory. Both paper configs use `/n09data/guerrini/glass_mock_v1.4.6.3_v2/results`, the suite supplying Paper II's 350 mocks. -Catalogue measurements, mock covariance, and inference preparation all read this key. +Catalogue measurements and mock covariance read this key. Generated mock products carry the parent directory's suite name (`glass_mock_v1.4.6.3_v2`) in their output paths, so suites don't share products. -`glass_mocks.version` labels the inference product version within a suite; `seed_range` selects the covariance and inference seeds. +`seed_range` selects the covariance seeds. +Inference consumes named SACC inputs under `inference.runs`; see [`cosmo_inference/README.md`](../cosmo_inference/README.md). The fine-grid measurement aggregators use the first 100 mocks. Sampled tau statistics depend on the real catalogue, not the mock suite, and retain their catalogue-derived inputs and output paths. ```bash snakemake -n --profile workflow/profiles/default -s workflow/Snakefile \ - covariance_glass_mock glass_mock_all_xi inference_glass_mocks \ + covariance_glass_mock glass_mock_all_xi \ --configfile papers/bmodes/config/config.yaml ``` diff --git a/workflow/inference/Snakefile b/workflow/inference/Snakefile new file mode 100644 index 000000000..a407a56ec --- /dev/null +++ b/workflow/inference/Snakefile @@ -0,0 +1,11 @@ +# Standalone entry point: needs only inference.runs, not measurement configs. +import os +import sys +from pathlib import Path + +sys.path.insert(0, os.path.realpath(str(Path(workflow.basedir).parent))) +import common + +common.inject_checkout_pythonpath(config) +container: common.resolve_container(config.get("container")) +include: "../rules/inference.smk" diff --git a/workflow/rules/inference.smk b/workflow/rules/inference.smk index 154fa9ac1..8d06937ee 100644 --- a/workflow/rules/inference.smk +++ b/workflow/rules/inference.smk @@ -1,178 +1,90 @@ -# Imports from Snakefile: FIDUCIAL, COSMO_INFERENCE, COSMO_VAL, covariance_path, redshift_path, fiducial_binning_suffix, pseudo_cl_tag -# NOTE: dormant subsystem. The file-name plumbing (config-driven paths + the -# producer-tagged pseudo-Cl names) is fixed and the DAG is valid, but it has not -# been run end-to-end. Reviving it still needs the FITS-CONTENT plumbing -# reconciled: cosmosis_fitting.py reads ELL/EE/BB + COVAR_FULL, while the -# producers write PSEUDO_CELL/ELL + COVAR_BB_BB. +# SACC → likelihood FITS → run INI → MPI chain → weighted summary. +# Each named run owns its inputs and priors; no catalogue is remeasured here. +INFERENCE = config.get("inference", {}) +INFERENCE_RUNS = INFERENCE.get("runs", {}) +INFERENCE_OUT = Path(INFERENCE.get("output_dir", "results/inference")).resolve() +INFERENCE_PATTERN = str(INFERENCE_OUT / "{run}") -# Output root for CosmoSIS data products + configs. COSMO_INFERENCE (common.py) -# already resolves to THIS repo's cosmo_inference dir, so the products land -# beside the code that builds them rather than in a contributor's home. -COSMO_INFERENCE_PROD = COSMO_INFERENCE -# Working directory for the cosmosis_fitting.py invocation — the same repo dir. -COSMO_INFERENCE_RUNDIR = str(COSMO_INFERENCE) -# External chain/mock locations are deployment-specific, so they live in config. -INFERENCE = config["inference"] -CHAINS_DIR = INFERENCE["chains_dir"] # CosmoSIS chain output root (real data) -GLASS_MOCK_CHAINS_DIR = INFERENCE["glass_mock_chains_dir"] # mock chain output root +def inference_run(w): + return INFERENCE_RUNS[w.run] -PSEUDO_CL_DIR = COSMO_VAL # producer (twopoint.smk) writes pseudo_cl* here -GLASS_MOCK_VERSION = config["glass_mocks"].get("version", "v0") -GLASS_MOCK_SEED_RANGE = config["glass_mocks"]["seed_range"] -GLASS_MOCK_FITS_PATTERN = str( - COSMO_INFERENCE_PROD - / f"data/glass_mocks/{GLASS_MOCK_SUITE}/{GLASS_MOCK_VERSION}/glass_mock_{{mock_id}}" - / f"cosmosis_glass_mock_{GLASS_MOCK_VERSION}_{{mock_id}}.fits" -) -GLASS_MOCK_CONFIG_PATTERN = str( - COSMO_INFERENCE_PROD - / f"cosmosis_config/output/{GLASS_MOCK_SUITE}/cosmosis_pipeline_glass_mocks_{GLASS_MOCK_VERSION}_glass_mock_{{mock_id}}.ini" -) - -# Fiducial harmonic-binning tag the pseudo-Cl producer (twopoint.smk) stamps -# into the filename. These are NOT inference_prep wildcards, so the consumer -# reads them from config to reconstruct the exact name the producer emits. -PSEUDO_CL_TAG = pseudo_cl_tag(config) - - -def pseudo_cl_assets(version): - """Return pseudo-Cl and covariance paths for the requested catalog version. +rule inference_fiducial: + input: + expand(INFERENCE_PATTERN + "/summary.json", run=list(INFERENCE_RUNS)) - The producer (twopoint.smk rules pseudo_cl / pseudo_cl_cov) writes - wildcard-tagged names; the consumer reconstructs them from the fiducial - harmonic-binning config so the requested path matches byte-for-byte. - """ - cl_path = PSEUDO_CL_DIR / f"pseudo_cl_{version}_{PSEUDO_CL_TAG}.fits" - cov_path = PSEUDO_CL_DIR / f"pseudo_cl_cov_{version}_{PSEUDO_CL_TAG}.fits" - return str(cl_path), str(cov_path) rule inference_prep: input: - # Processed covariance matrix - use centralized covariance_path() - cov_matrix=lambda w: covariance_path(w.version, min_sep=w.min_sep, max_sep=w.max_sep, nbins=w.nbins), - # Xi FITS files - xi_plus=str(COSMO_VAL / "xi_plus_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), - xi_minus=str(COSMO_VAL / "xi_minus_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), - # n(z) file: the catalogue entry's - nz_file=lambda w: redshift_path(w.version), - # rho/tau stats - rho_stats=lambda w: cv_rho_stats(w.version, w), - tau_stats=lambda w: cv_tau_stats(w.version, w), - # tau covariance (tracked as dependency) - tau_cov=lambda w: cv_cov_tau(w.version, w), - pseudo_cl=lambda w: pseudo_cl_assets(w.version)[0], - pseudo_cl_cov=lambda w: pseudo_cl_assets(w.version)[1], + sacc=lambda w: inference_run(w)["sacc"], + covariance=lambda w: [inference_run(w)["covariance"]] if inference_run(w).get("covariance") else [], + rho=lambda w: [inference_run(w)["rho_stats"]] if inference_run(w).get("rho_stats") else [], + tau=lambda w: [inference_run(w)["tau_stats"]] if inference_run(w).get("tau_stats") else [], output: - fits_file=str( - COSMO_INFERENCE_PROD - / "data/{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}/cosmosis_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits" - ), - config_file=str( - COSMO_INFERENCE_PROD - / "cosmosis_config/output/cosmosis_pipeline_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.ini" - ) + fits=INFERENCE_PATTERN + "/data.fits", + provenance=INFERENCE_PATTERN + "/input.json", params: - cosmosis_root="{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}", - data_dir=f"{CHAINS_DIR}/{{version}}_minsep={{min_sep}}_maxsep={{max_sep}}_nbins={{nbins}}_npatch={{npatch}}", - output_root=str(COSMO_INFERENCE_PROD), + run=inference_run, threads: 1 resources: mem_mb=8000, runtime=10, - shell: - """ - cd {COSMO_INFERENCE_RUNDIR} + script: + "../scripts/inference_prepare.py" - # Run inference preparation step with cosmosis_fitting.py - python scripts/cosmosis_fitting.py \ - --cosmosis-root {params.cosmosis_root} \ - --nz-file {input.nz_file} \ - --data-dir {params.data_dir} \ - --output-root {params.output_root} \ - --xi {input.xi_plus} {input.xi_minus} \ - --cov-xi {input.cov_matrix} \ - --use-rho-tau \ - --rho-stats {input.rho_stats} \ - --tau-stats {input.tau_stats} \ - --cov-tau {input.tau_cov} \ - --cl-file {input.pseudo_cl} \ - --cov-cl {input.pseudo_cl_cov} - """ - -rule inference_fiducial: - input: - # Use the same output patterns as inference_prep with FIDUCIAL params - rules.inference_prep.output.fits_file.format( - version=FIDUCIAL["version"], - min_sep=FIDUCIAL["min_sep"], max_sep=FIDUCIAL["max_sep"], - nbins=FIDUCIAL["nbins"], npatch=FIDUCIAL["npatch"] - ), - rules.inference_prep.output.config_file.format( - version=FIDUCIAL["version"], - min_sep=FIDUCIAL["min_sep"], max_sep=FIDUCIAL["max_sep"], - nbins=FIDUCIAL["nbins"], npatch=FIDUCIAL["npatch"] - ) - - -rule inference_glass_mocks: +rule inference_config: input: - expand(GLASS_MOCK_FITS_PATTERN, mock_id=[f"{i:05d}" for i in range(GLASS_MOCK_SEED_RANGE[0], GLASS_MOCK_SEED_RANGE[1] + 1)]) + fits=rules.inference_prep.output.fits, + template=lambda w: inference_run(w)["template"], + values=lambda w: inference_run(w)["values"], + priors=lambda w: inference_run(w)["priors"], + output: + ini=INFERENCE_PATTERN + "/pipeline.ini", + params: + run=inference_run, + chain=lambda w: INFERENCE_PATTERN.format(run=w.run) + "/chain.txt", + threads: 1 + resources: + mem_mb=1000, + runtime=5, + script: + "../scripts/inference_config.py" -rule inference_prep_glass_mock: +rule inference_chain: input: - xi=f"{GLASS_MOCK_DIR}/xi_glass_mock_{{mock_id}}_4096_nbins=20.fits", - # Use centralized covariance_path() with fiducial mock version - cov_matrix=covariance_path(FIDUCIAL["mock_version"]), - # n(z) file - nz_file=redshift_path(FIDUCIAL["mock_version"]), - # Rho/tau stats: rho from real data, tau sampled - rho_stats=cv_rho_stats(FIDUCIAL["mock_version"]), - tau_stats="results/glass_mock_rhotau_samples/{mock_id}/tau_stats_sampled.fits", - # Tau covariance (real data) - tau_cov=cv_cov_tau(FIDUCIAL["mock_version"]), - # C_ell data for dual config generation - cl_file=f"{GLASS_MOCK_DIR}/cl_glass_mock_{{mock_id}}_4096.npy", - cl_cov=pseudo_cl_assets(FIDUCIAL["mock_version"])[1], + ini=rules.inference_config.output.ini, + fits=rules.inference_prep.output.fits, + values=lambda w: inference_run(w)["values"], + priors=lambda w: inference_run(w)["priors"], output: - fits_file=GLASS_MOCK_FITS_PATTERN, - config_file=GLASS_MOCK_CONFIG_PATTERN, + chain=INFERENCE_PATTERN + "/chain.txt", + log: + INFERENCE_PATTERN + "/cosmosis.log", params: - cosmosis_root=f"glass_mock_{GLASS_MOCK_VERSION}_{{mock_id}}", - data_dir=f"{GLASS_MOCK_CHAINS_DIR}/{GLASS_MOCK_SUITE}/glass_mock_{GLASS_MOCK_VERSION}_{{mock_id}}", - output_root=str(COSMO_INFERENCE_PROD), - output_basename=f"glass_mocks/{GLASS_MOCK_SUITE}/{GLASS_MOCK_VERSION}/glass_mock_{{mock_id}}", - threads: 1 + mpi_args=INFERENCE.get("mpi_args", ""), + threads: lambda w: inference_run(w).get("threads", 8) resources: - mem_mb=8000, - runtime=10, + mem_mb=lambda w: inference_run(w).get("mem_mb", 16000), + runtime=lambda w: inference_run(w).get("runtime", 720), shell: - """ - cd {COSMO_INFERENCE_RUNDIR} + "eval \"$(python -m cosmosis.configure)\"; " + # Each rule is one allocated node; the image has no srun executable. + "PRTE_MCA_ras=^slurm PRTE_MCA_plm=^slurm " + "OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 mpiexec {params.mpi_args} -n {threads} " + "cosmosis --mpi {input.ini:q} > {log:q} 2>&1" - mkdir -p {params.data_dir} - python scripts/cosmosis_fitting.py --mock \ - --cosmosis-root {params.cosmosis_root} \ - --nz-file {input.nz_file} \ - --data-dir {params.data_dir} \ - --output-root {params.output_root} \ - --output-basename {params.output_basename} \ - --xi {input.xi} \ - --cov-xi {input.cov_matrix} \ - --use-rho-tau \ - --rho-stats {input.rho_stats} \ - --tau-stats {input.tau_stats} \ - --cov-tau {input.tau_cov} \ - --cl-file {input.cl_file} \ - --cov-cl {input.cl_cov} - """ - -localrules: - inference_prep, - inference_prep_glass_mock, - inference_fiducial, - inference_glass_mocks, +rule inference_summary: + input: + chain=rules.inference_chain.output.chain, + output: + summary=INFERENCE_PATTERN + "/summary.json", + threads: 1 + resources: + mem_mb=2000, + runtime=5, + script: + "../scripts/inference_summary.py" diff --git a/workflow/scripts/inference_config.py b/workflow/scripts/inference_config.py new file mode 100644 index 000000000..22615514f --- /dev/null +++ b/workflow/scripts/inference_config.py @@ -0,0 +1,18 @@ +"""Render a run INI from an explicit template, values and priors.""" + +from sp_validation.inference import render + +run = snakemake.params.run +render( + snakemake.input.template, + snakemake.output.ini, + fits_path=snakemake.input.fits, + values=snakemake.input.values, + priors=snakemake.input.priors, + chain=snakemake.params.chain, + n_bins=run.get("n_bins", 1), + sampler=run.get("sampler", "nautilus"), + settings=run.get("sampler_settings"), + scale_cuts=run.get("scale_cuts"), + overrides=run.get("overrides"), +) diff --git a/workflow/scripts/inference_prepare.py b/workflow/scripts/inference_prepare.py new file mode 100644 index 000000000..68574e963 --- /dev/null +++ b/workflow/scripts/inference_prepare.py @@ -0,0 +1,18 @@ +"""Convert a reporting SACC and record the likelihood's input provenance.""" + +import json +from pathlib import Path + +from sp_validation.inference import prepare + +run = snakemake.params.run +provenance = prepare( + snakemake.input.sacc, + snakemake.output.fits, + n_bins=run.get("n_bins", 1), + statistic=run.get("statistic", "xi"), + covariance=run.get("covariance"), + rho_stats=run.get("rho_stats"), + tau_stats=run.get("tau_stats"), +) +Path(snakemake.output.provenance).write_text(json.dumps(provenance, indent=2) + "\n") diff --git a/workflow/scripts/inference_summary.py b/workflow/scripts/inference_summary.py new file mode 100644 index 000000000..a7432fe9a --- /dev/null +++ b/workflow/scripts/inference_summary.py @@ -0,0 +1,5 @@ +"""Summarize the completed weighted CosmoSIS chain.""" + +from sp_validation.inference import summarize + +summarize(snakemake.input.chain, snakemake.output.summary) diff --git a/workflow/tests/test_glass_suite.py b/workflow/tests/test_glass_suite.py index 50e4e84b5..bd2bad1d8 100644 --- a/workflow/tests/test_glass_suite.py +++ b/workflow/tests/test_glass_suite.py @@ -16,7 +16,6 @@ def test_paper_glass_suite(paper): assert config["glass_mocks"]["data_dir"] == ( "/n09data/guerrini/glass_mock_v1.4.6.3_v2/results" ) - assert "glass_mock_data_dir" not in config["inference"] assert config["glass_mocks"]["seed_range"] == [1, 350] @@ -64,7 +63,6 @@ def test_glass_rules_share_suite(toy, tmp_path, suite): "covariance_glass_mock", "glass_mock_all_xi", "glass_mock_all_pseudo_cl", - "inference_glass_mocks", "mock_cosebis_bias_test", "--configfile", str(config_path), @@ -75,7 +73,6 @@ def test_glass_rules_share_suite(toy, tmp_path, suite): ("covariance_glass_mock", 1, 4), ("glass_mock_xi_fine", 100, 1), ("glass_mock_pseudo_cl", 100, 1), - ("inference_prep_glass_mock", 2, 2), ): selected = [job for job in jobs if job.rule == rule] assert len(selected) == count diff --git a/workflow/tests/test_inference_dag.py b/workflow/tests/test_inference_dag.py new file mode 100644 index 000000000..1fcfa5bb9 --- /dev/null +++ b/workflow/tests/test_inference_dag.py @@ -0,0 +1,55 @@ +"""Inference starts from declared SACCs, not obsolete ξ-file producers.""" + +import yaml +from conftest import parse_jobs + + +def test_sacc_to_summary_dag(toy, tmp_path): + inputs = {} + for key in ( + "sacc", + "covariance", + "template", + "values", + "priors", + "rho_stats", + "tau_stats", + ): + path = tmp_path / key + path.touch() + inputs[key] = str(path) + inputs.update(n_bins=1, threads=8) + out = tmp_path / "out" + config = tmp_path / "inference.yaml" + config.write_text( + yaml.safe_dump( + {"inference": {"output_dir": str(out), "runs": {"mock": inputs}}} + ) + ) + result = toy.snakemake( + "-n", + "-s", + str(toy.root / "workflow/inference/Snakefile"), + "inference_fiducial", + "--configfile", + str(config), + ) + assert result.returncode == 0, result.stdout + jobs = {j.rule: j for j in parse_jobs(result.stdout)} + assert set(jobs) == { + "inference_prep", + "inference_config", + "inference_chain", + "inference_summary", + "inference_fiducial", + } + assert set(jobs["inference_prep"].input) == { + inputs[k] for k in ("sacc", "covariance", "rho_stats", "tau_stats") + } + assert set(jobs["inference_config"].input) == { + str(out / "mock/data.fits"), + *[inputs[k] for k in ("template", "values", "priors")], + } + assert str(out / "mock/chain.txt") in jobs["inference_summary"].input + assert str(out / "mock/summary.json") in jobs["inference_fiducial"].input + assert all("xi_plus_" not in p for job in jobs.values() for p in job.input)