diff --git a/Dockerfile b/Dockerfile index c791877c..ff0ea389 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,5 +1,5 @@ # Development image with more bells and whistles -FROM ghcr.io/cosmostat/shapepipe:develop +FROM ghcr.io/cosmostat/shapepipe:im_sims RUN apt-get update -y --quiet --fix-missing && \ apt-get dist-upgrade -y --quiet --fix-missing && \ diff --git a/config/calibration/mask_v1.X.4_im_sim.yaml b/config/calibration/mask_v1.X.4_im_sim.yaml new file mode 100644 index 00000000..794566d4 --- /dev/null +++ b/config/calibration/mask_v1.X.4_im_sim.yaml @@ -0,0 +1,70 @@ +# Config file for masking and calibration. +# Standard cuts without coverage mask, type v1.X.4. + +# General parameters (can also given on command line) +params: + input_path: shape_catalog_comprehensive_ngmix.hdf5 + cmatrices: False + sky_regions: False + verbose: True + +# Masks +## Using columns in 'dat' group (ShapePipe flags) +dat: + # SExtractor flags + - col_name: FLAGS + label: SE FLAGS + kind: equal + value: 0 + + # Number of epochs + - col_name: N_EPOCH + label: r"$n_{\rm epoch}$" + kind: greater_equal + value: 2 + + # Magnitude range + - col_name: mag + label: mag range + kind: range + value: [15, 30] + + # ngmix flags + - col_name: NGMIX_MOM_FAIL + label: "ngmix moments failure" + kind: equal + value: 0 + + # invalid PSF ellipticities + - col_name: NGMIX_ELL_PSFo_NOSHEAR_0 + label: "bad PSF ellipticity comp 1" + kind: not_equal + value: -10 + - col_name: NGMIX_ELL_PSFo_NOSHEAR_1 + label: "bad PSF ellipticity comp 2" + kind: not_equal + value: -10 + +# Metacal parameters +metacal: + # Ellipticity dispersion + sigma_eps_prior: 0.34 + + # Signal-to-noise range + gal_snr_min: 10 + gal_snr_max: 500 + + # Relative-size (hlr / hlr_psf) range + gal_rel_size_min: 0.5 + gal_rel_size_max: 3 + + # Correct relative size for ellipticity? + gal_size_corr_ell: False + + # Weight for global response matrix, None for unweighted mean. + # Unweighted for image sims: no weights anywhere in sim m-bias (#227). + global_R_weight: null + + # Subtract additive bias (mean shear)? Use False for constant-shear + # image sims + additive_correction: False diff --git a/config/calibration/mask_v1.X.9_im_sim.overlay.yaml b/config/calibration/mask_v1.X.9_im_sim.overlay.yaml new file mode 100644 index 00000000..94c5b9ae --- /dev/null +++ b/config/calibration/mask_v1.X.9_im_sim.overlay.yaml @@ -0,0 +1,97 @@ +# Declared delta: image-sim mask/calibration config vs the data config. +# +# The image-sim calibration reuses the *data* mask config +# (mask_v1.X.9.yaml) and changes only what the sims genuinely differ on. +# Rather than maintain a second full copy that can silently drift from the +# base, this overlay states -- as a list of block operations on the base +# file -- exactly which pieces the sims drop or change, and one line of why +# for each. `im_compose_mask.py` applies these ops to mask_v1.X.9.yaml and +# reproduces mask_v1.X.9_im_sim.yaml byte-for-byte; a test locks that, so the +# runtime file and this declaration cannot diverge. +# +# Each op anchors to a block of base text (matched verbatim, and required to +# occur exactly once) and either drops it or replaces it. `why` is prose for +# the human reader; the compose ignores it. Ordering follows the base file. + +base: mask_v1.X.9.yaml + +ops: + # --- params ------------------------------------------------------------ + - why: >- + Sims read the ShapePipe FITS catalogue staged in the run dir, not the + survey-wide comprehensive HDF5 on /n17data. + replace: | + input_path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_comprehensive_struc_2024_v1.X.c.hdf5 + with: | + input_path: shape_catalog_comprehensive_ngmix.fits + + # --- dat cuts ---------------------------------------------------------- + - why: >- + No ShapePipe coverage/mask flags on sims: IMAFLAGS_ISO is a survey + artefact (external masks, bright-star haloes) the sims do not carry. + drop: |2 + + # ShapePipe flags + - col_name: IMAFLAGS_ISO + label: SP mask + kind: equal + value: 0 + + - why: >- + Same cuts, but flag the grammar: the sims run the ShapePipe-v2 PSF + columns (scalar G1/G2), so the comment is made explicit here. + replace: |2 + # invalid PSF ellipticities + with: |2 + # invalid PSF ellipticities (ShapePipe-v2 grammar: scalar G1/G2 components) + + # --- dat_ext (post-processing / coverage masks) ------------------------ + - why: >- + No coverage masks on sims: the whole dat_ext group (Stars, manual mask, + r-band footprint, Maximask) is survey post-processing with no analogue + in the simulated tiles. + drop: |2 + + ## Using columns in 'dat_ext' group (post-processing flags) + dat_ext: + + # Stars + - col_name: 4_Stars + label: "Stars" + kind: equal + value: False + + # Manual mask + - col_name: 8_Manual + label: "manual mask" + kind: equal + value: False + + # r-band footprint + - col_name: 64_r + label: "r-band imaging" + kind: equal + value: False + + # Maximask + - col_name: 1024_Maximask + label: "maximask" + kind: equal + value: False + + # --- metacal ----------------------------------------------------------- + - why: >- + Unweighted for image sims: no weights anywhere in the sim m-bias (#227), + and the w_des-weighted global R had only N_eff ~ 20-100 objects. + replace: |2 + # Weight for global response matrix, None for unweighted mean + global_R_weight: w + with: |2 + # Weight for global response matrix, None for unweighted mean. + # Unweighted for image sims: no weights anywhere in sim m-bias (#227), + # and the w_des-weighted R had N_eff ~ 20-100 objects. + global_R_weight: null + + # Subtract additive bias (mean shear)? Use False for constant-shear + # image sims + additive_correction: False diff --git a/config/calibration/mask_v1.X.9_im_sim.yaml b/config/calibration/mask_v1.X.9_im_sim.yaml new file mode 100644 index 00000000..1c6ff672 --- /dev/null +++ b/config/calibration/mask_v1.X.9_im_sim.yaml @@ -0,0 +1,77 @@ +# Config file for masking and calibration. +# Less conservative cuts, type v1.X.9 (e.g. for matching with spectroscopic sample). + +# General parameters (can also given on command line) +params: + input_path: shape_catalog_comprehensive_ngmix.fits + cmatrices: False + sky_regions: False + verbose: True + +# Masks +## Using columns in 'dat' group (ShapePipe flags) +dat: + # SExtractor flags + - col_name: FLAGS + label: SE FLAGS + kind: smaller_equal + value: 2 + + # Duplicate objects + - col_name: overlap + label: tile overlap + kind: equal + value: True + + # Number of epochs + - col_name: N_EPOCH + label: r"$n_{\rm epoch}$" + kind: greater_equal + value: 1 + + # Magnitude range + - col_name: mag + label: mag range + kind: range + value: [15, 30] + + # ngmix flags + - col_name: NGMIX_MCAL_TYPES_FAIL + label: "ngmix moments failure" + kind: equal + value: 0 + + # invalid PSF ellipticities (ShapePipe-v2 grammar: scalar G1/G2 components) + - col_name: NGMIX_G1_PSF_ORIG_NOSHEAR + label: "bad PSF ellipticity comp 1" + kind: not_equal + value: -10 + - col_name: NGMIX_G2_PSF_ORIG_NOSHEAR + label: "bad PSF ellipticity comp 2" + kind: not_equal + value: -10 + +# Metacal parameters +metacal: + # Ellipticity dispersion + sigma_eps_prior: 0.34 + + # Signal-to-noise range + gal_snr_min: 5 + gal_snr_max: 500 + + # Relative-size (hlr / hlr_psf) range + gal_rel_size_min: 0.25 + gal_rel_size_max: 10 + + # Correct relative size for ellipticity? + gal_size_corr_ell: False + + # Weight for global response matrix, None for unweighted mean. + # Unweighted for image sims: no weights anywhere in sim m-bias (#227), + # and the w_des-weighted R had N_eff ~ 20-100 objects. + global_R_weight: null + + # Subtract additive bias (mean shear)? Use False for constant-shear + # image sims + additive_correction: False diff --git a/scripts/calibration/calibrate_comprehensive_cat.py b/scripts/calibration/calibrate_comprehensive_cat.py index 4b15ed89..4928df9b 100644 --- a/scripts/calibration/calibrate_comprehensive_cat.py +++ b/scripts/calibration/calibrate_comprehensive_cat.py @@ -96,8 +96,13 @@ ) # %% +additive_correction = cm.get("additive_correction", True) +if not additive_correction: + print("Additive bias correction disabled (additive_correction: False)") + g_corr_mc, g_uncorr, w, mask_metacal, c, c_err = calibration.get_calibrated_m_c( - gal_metacal + gal_metacal, + additive_correction=additive_correction, ) num_ok = len(g_corr_mc[0]) diff --git a/scripts/calibration/extract_info.py b/scripts/calibration/extract_info.py index 1cbf836f..123c7437 100644 --- a/scripts/calibration/extract_info.py +++ b/scripts/calibration/extract_info.py @@ -32,6 +32,7 @@ import os import sys +import h5py import numpy as np from astropy.io import fits @@ -51,8 +52,11 @@ # ### Create and open output files and directories -make_out_dirs(output_dir, plot_dir, [], verbose=verbose) -stats_file = open_stats_file(plot_dir, stats_file_name) +os.makedirs(output_dir, exist_ok=True) +stats_file = open_stats_file(output_dir, stats_file_name) + +output_shape_cat_stem = output_shape_cat_base +output_ext = output_format # ## 2. Load data # @@ -107,6 +111,8 @@ tile_IDs, path_tile_ID, path_found_ID, path_missing_ID, verbose=verbose ) +print_stats(f"Tiles in input catalogue: {n_found}", stats_file, verbose=verbose) + # ### Load star catalogue if star_cat_path: @@ -146,42 +152,40 @@ verbose=verbose, ) -# #### Refine: Match to valid, unflagged galaxy sample - -# + -# Flags to indicate valid star sample - -m_star = ( - (dd["FLAGS"][ind_star] == 0) - & (dd["IMAFLAGS_ISO"][ind_star] == 0) - & (dd["NGMIX_MCAL_FLAGS"][ind_star] == 0) - & (dd["NGMIX_G1_PSF_ORIG_NOSHEAR"][ind_star] != -10) -) - -ra_star, dec_star, g_star_psf = spv_cat.match_subsample( - dd, - ind_star, - m_star, - [col_name_ra, col_name_dec], - key_PSF_g1, - key_PSF_g2, - n_star_tot, - stats_file, - verbose=verbose, -) -# - + # Flags to indicate valid star sample. + # Star matching, PSF-catalogue output and star metacalibration are all + # diagnostics that require an input star catalogue; the image-simulation + # pipeline has none (star_cat_path is None), so every star-dependent block + # below is guarded and simply skipped for the sims. + + m_star = ( + (dd["FLAGS"][ind_star] == 0) + & (dd["IMAFLAGS_ISO"][ind_star] == 0) + & (dd["NGMIX_MCAL_FLAGS"][ind_star] == 0) + & (dd["NGMIX_G1_PSF_ORIG_NOSHEAR"][ind_star] != -10) + ) -# MKDEBUG: Moved from end of this script + ra_star, dec_star, g_star_psf = spv_cat.match_subsample( + dd, + ind_star, + m_star, + [col_name_ra, col_name_dec], + key_PSF_g1, + key_PSF_g2, + n_star_tot, + stats_file, + verbose=verbose, + ) -# ### Write PSF catalogue with multi-epoch shapes from shape measurement methods + # ### Write PSF catalogue with multi-epoch shapes from shape measurement methods -spv_cat.write_PSF_cat( - f"{output_PSF_cat_base}_{shape}.fits", - ra_star, - dec_star, - g_star_psf[0], - g_star_psf[1], -) + spv_cat.write_PSF_cat( + f"{output_PSF_cat_base}_{shape}.fits", + ra_star, + dec_star, + g_star_psf[0], + g_star_psf[1], + ) # ## Check for objects with invalid PSF @@ -253,8 +257,9 @@ if verbose: print("Writing comprehensive catalogue...") +comprehensive_cat_path = f"{output_shape_cat_stem}_comprehensive_{shape}{output_ext}" spv_cat.write_shape_catalog( - f"{output_shape_cat_base}_comprehensive_{shape}.fits", + comprehensive_cat_path, ra_all, dec_all, iv_w, @@ -265,6 +270,17 @@ add_cols=ext_cols_pre_cal, add_cols_format=add_cols_pre_cal_format, ) + +# Write tile count to HDF5 attributes +if output_ext == ".hdf5": + try: + with h5py.File(comprehensive_cat_path, "a") as hf: + hf.attrs["n_tiles"] = n_found + if verbose: + print(f" Added n_tiles={n_found} to HDF5 attributes") + except Exception as e: + if verbose: + print(f" Warning: could not add n_tiles attribute: {e}") # - do_selection_calibration = False @@ -275,6 +291,7 @@ else: if verbose: print("Continuing with selection and calibration") + os.makedirs(os.path.join(output_dir, plot_dir), exist_ok=True) # ## 4. Select galaxies @@ -628,23 +645,24 @@ # ## Metacalibration for stars -star_metacal = metacal(dd[ind_star], m_star, masking_type="star", verbose=verbose) +if star_cat_path: + star_metacal = metacal(dd[ind_star], m_star, masking_type="star", verbose=verbose) -# #### Number density + # #### Number density -# + -# mask for 'no shear' images + # + + # mask for 'no shear' images -mask_ns_stars = star_metacal.mask_dict["ns"] -n_star = len(star_metacal.ns["g1"][mask_ns_stars]) + mask_ns_stars = star_metacal.mask_dict["ns"] + n_star = len(star_metacal.ns["g1"][mask_ns_stars]) -print_stats(f"Number of stars = {n_star}", stats_file, verbose=verbose) -print_stats( - "Star density = {:.2f} stars/deg2".format(n_star / area_deg2), - stats_file, - verbose=verbose, -) -# - + print_stats(f"Number of stars = {n_star}", stats_file, verbose=verbose) + print_stats( + "Star density = {:.2f} stars/deg2".format(n_star / area_deg2), + stats_file, + verbose=verbose, + ) + # - # ## Additive bias # Use raw, uncorrected ellipticities. @@ -730,20 +748,21 @@ print_stats(rs, stats_file, verbose=verbose) # + -print_stats("stars:", stats_file, verbose=verbose) +if star_cat_path: + print_stats("stars:", stats_file, verbose=verbose) -print_stats("total response matrix:", stats_file, verbose=verbose) -rs = np.array2string(star_metacal.R) -print_stats(rs, stats_file, verbose=verbose) + print_stats("total response matrix:", stats_file, verbose=verbose) + rs = np.array2string(star_metacal.R) + print_stats(rs, stats_file, verbose=verbose) -print_stats("shear response matrix:", stats_file, verbose=verbose) -R_shear_stars = np.mean(star_metacal.R_shear, 2) -rs = np.array2string(R_shear_stars) -print_stats(rs, stats_file, verbose=verbose) + print_stats("shear response matrix:", stats_file, verbose=verbose) + R_shear_stars = np.mean(star_metacal.R_shear, 2) + rs = np.array2string(R_shear_stars) + print_stats(rs, stats_file, verbose=verbose) -print_stats("selection response matrix:", stats_file, verbose=verbose) -rs = np.array2string(star_metacal.R_selection) -print_stats(rs, stats_file, verbose=verbose) + print_stats("selection response matrix:", stats_file, verbose=verbose) + rs = np.array2string(star_metacal.R_selection) + print_stats(rs, stats_file, verbose=verbose) # - # ### Plot distribution of response matrix elements @@ -757,14 +776,11 @@ linestyles = ["-", "-", ":", ":"] # + -labels = ["$R_{11}$ galaxies", "$R_{22}$ galaxies", "$R_{11}$ stars", "$R_{22}$ stars"] - -xs = [ - gal_metacal.R_shear[0, 0], - gal_metacal.R_shear[1, 1], - star_metacal.R_shear[0, 0], - star_metacal.R_shear[1, 1], -] +labels = ["$R_{11}$ galaxies", "$R_{22}$ galaxies"] +xs = [gal_metacal.R_shear[0, 0], gal_metacal.R_shear[1, 1]] +if star_cat_path: + labels += ["$R_{11}$ stars", "$R_{22}$ stars"] + xs += [star_metacal.R_shear[0, 0], star_metacal.R_shear[1, 1]] title = shape out_name = f"R_{shape}_diag.pdf" @@ -779,19 +795,16 @@ x_range, n_bin, out_path, - colors=colors, - linestyles=linestyles, + colors=colors[: len(xs)], + linestyles=linestyles[: len(xs)], ) # + -labels = ["$R_{12}$ galaxies", "$R_{21}$ galaxies", "$R_{12}$ stars", "$R_{21}$ stars"] - -xs = [ - gal_metacal.R_shear[0, 1], - gal_metacal.R_shear[1, 0], - star_metacal.R_shear[0, 1], - star_metacal.R_shear[1, 0], -] +labels = ["$R_{12}$ galaxies", "$R_{21}$ galaxies"] +xs = [gal_metacal.R_shear[0, 1], gal_metacal.R_shear[1, 0]] +if star_cat_path: + labels += ["$R_{12}$ stars", "$R_{21}$ stars"] + xs += [star_metacal.R_shear[0, 1], star_metacal.R_shear[1, 0]] title = shape out_name = f"R_{shape}_offdiag.pdf" out_path = os.path.join(plot_dir, out_name) @@ -805,8 +818,8 @@ x_range, n_bin, out_path, - colors=colors, - linestyles=linestyles, + colors=colors[: len(xs)], + linestyles=linestyles[: len(xs)], ) # - @@ -845,49 +858,50 @@ ) # + -xs = [star_metacal.ns["g1"][mask_ns_stars], star_metacal.ns["g2"][mask_ns_stars]] -weights = [star_metacal.ns["w"][mask_ns_stars]] * 2 - -title = "stars" -out_name = f"ell_stars_{shape}.pdf" -out_path = os.path.join(plot_dir, out_name) - -plot_histograms( - xs, - labels, - title, - x_label, - y_label, - x_range, - n_bin, - out_path, - weights=weights, - colors=colors, - linestyles=linestyles, -) -# - - -x_range = (-0.15, 0.15) -n_bin = 250 - -# + -xs = [dd[key_PSF_g1][mask_ns_stars], dd[key_PSF_g2][mask_ns_stars]] -title = "PSF" -out_name = f"ell_PSF_{shape}.pdf" -out_path = os.path.join(plot_dir, out_name) - -plot_histograms( - xs, - labels, - title, - x_label, - y_label, - x_range, - n_bin, - out_path, - colors=colors, - linestyles=linestyles, -) +if star_cat_path: + xs = [star_metacal.ns["g1"][mask_ns_stars], star_metacal.ns["g2"][mask_ns_stars]] + weights = [star_metacal.ns["w"][mask_ns_stars]] * 2 + + title = "stars" + out_name = f"ell_stars_{shape}.pdf" + out_path = os.path.join(plot_dir, out_name) + + plot_histograms( + xs, + labels, + title, + x_label, + y_label, + x_range, + n_bin, + out_path, + weights=weights, + colors=colors, + linestyles=linestyles, + ) + # - + + x_range = (-0.15, 0.15) + n_bin = 250 + + # + + xs = [dd[key_PSF_g1][mask_ns_stars], dd[key_PSF_g2][mask_ns_stars]] + title = "PSF" + out_name = f"ell_PSF_{shape}.pdf" + out_path = os.path.join(plot_dir, out_name) + + plot_histograms( + xs, + labels, + title, + x_label, + y_label, + x_range, + n_bin, + out_path, + colors=colors, + linestyles=linestyles, + ) # - # ## Magnitudes @@ -951,7 +965,7 @@ # ### Write basic shape catalogue spv_cat.write_shape_catalog( - f"{output_shape_cat_base}_{shape}.fits", + f"{output_shape_cat_stem}_{shape}{output_ext}", ra, dec, w, @@ -990,7 +1004,7 @@ # Extended catalogue with SNR, individual R matrices, ext_cols spv_cat.write_shape_catalog( - f"{output_shape_cat_base}_extended_{shape}.fits", + f"{output_shape_cat_stem}_extended_{shape}{output_ext}", ra, dec, w, @@ -1020,4 +1034,4 @@ ra = dd["RA"][cut_overlap] dec = dd["DEC"][cut_overlap] tile_id = dd["TILE_ID"][cut_overlap] - write_galaxy_cat(f"{output_shape_cat_base}.fits", ra, dec, tile_id) + write_galaxy_cat(f"{output_shape_cat_stem}{output_ext}", ra, dec, tile_id) diff --git a/scripts/calibration/params.py b/scripts/calibration/params.py index 2d0bd45f..5b2ef0c9 100644 --- a/scripts/calibration/params.py +++ b/scripts/calibration/params.py @@ -105,6 +105,9 @@ ## Output +### Output file format extension: '.fits' or '.hdf5' +output_format = ".hdf5" + ### Additional output columns add_cols = [ "FLUX_RADIUS", diff --git a/scripts/compute_m_bias_image_sims.py b/scripts/compute_m_bias_image_sims.py new file mode 100644 index 00000000..04ca49b7 --- /dev/null +++ b/scripts/compute_m_bias_image_sims.py @@ -0,0 +1,361 @@ +#!/usr/bin/env python +"""Compute multiplicative and additive shear bias from image simulations. + +Usage: + compute_m_bias_image_sims.py -c config.yaml [-v] [--cumulative] [--n_tiles N] +""" + +import argparse +import os +import sys + +# Configure matplotlib for non-interactive backend +import matplotlib +import numpy as np +import yaml + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt + +from sp_validation.image_sims import ImageSimMBias + + +def to_python(obj): + """Recursively cast numpy scalars/arrays to plain Python types. + + ``ImageSimMBias.run`` returns a nested dict of numpy floats; dumping those + to YAML with ``yaml.dump`` writes opaque ``!!python/object`` binary tags. A + recursive pass down the results tree (dicts, lists, arrays, scalars) leaves + a clean, human-readable, ``safe_load``-able document. + """ + if isinstance(obj, dict): + return {key: to_python(val) for key, val in obj.items()} + if isinstance(obj, (list, tuple)): + return [to_python(val) for val in obj] + if isinstance(obj, np.ndarray): + return float(obj.item()) if obj.size == 1 else obj.tolist() + if isinstance(obj, (np.integer, np.floating)): + return float(obj) + return obj + + +def parse_args(): + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("-c", "--config", required=True, help="config YAML file") + p.add_argument("-v", "--verbose", action="store_true", help="verbose output") + p.add_argument( + "--cumulative", + action="store_true", + default=True, + help="track convergence as tiles accumulate (default: True)", + ) + p.add_argument( + "--n_tiles", type=int, help="number of tiles (auto-detected if not given)" + ) + return p.parse_args() + + +def get_n_tiles(grids_dir, num): + """Detect number of tiles from final_cat HDF5 files.""" + try: + import h5py + + # Count tiles in first sim's final_cat + for sim in ["1z2z_grid", "1m2z_grid", "1p2z_grid", "1z2m_grid", "1z2p_grid"]: + sim_name = f"{sim}_{num}" + final_cat = os.path.join(grids_dir, sim_name, f"final_cat_{sim_name}.hdf5") + if os.path.isfile(final_cat): + with h5py.File(final_cat, "r") as hf: + if "patches" in hf: + n_tiles = sum( + 1 for patch in hf["patches"] for _ in hf[f"patches/{patch}"] + ) + return n_tiles + except Exception: + pass + return None + + +def update_cumulative_file(cumulative_path, n_tiles, results): + """Update the cumulative m/c bias tracking file. + + Writes ``results`` under the ``n_tiles`` key, *overwriting* an existing + entry for that count. The earlier behaviour silently skipped when the key + was already present, which meant a re-run against fixed catalogues left the + old (possibly wrong) number in place -- a stale value masquerading as + current. A fresh run is the authority for its tile count, so it overwrites. + + Returns ``True`` when a new key was added, ``False`` when an existing entry + was overwritten (the file is written either way). + """ + if os.path.isfile(cumulative_path): + with open(cumulative_path) as f: + try: + cumulative = yaml.safe_load(f) or {} + except yaml.YAMLError: + # Legacy file written before the to_python cleanup: it carries + # numpy python-object tags that safe_load rejects. Load it + # unsafely, then the to_python pass on write heals it in place. + f.seek(0) + cumulative = yaml.unsafe_load(f) or {} + else: + cumulative = {} + + is_new = str(n_tiles) not in cumulative + cumulative[str(n_tiles)] = results + with open(cumulative_path, "w") as f: + yaml.dump(to_python(cumulative), f, default_flow_style=False) + return is_new + + +def plot_convergence(cumulative_path, diagnostics_dir): + """Create convergence plots: m/c vs n_tiles and errors vs n_tiles.""" + os.makedirs(diagnostics_dir, exist_ok=True) + + try: + with open(cumulative_path) as f: + cumulative = yaml.safe_load(f) + except Exception as e: + print(f"Warning: could not read cumulative file {cumulative_path}: {e}") + return + + if not cumulative: + print("No cumulative data yet, skipping plots") + return + + # Sort by n_tiles + n_tiles_list = sorted([int(k) for k in cumulative.keys()]) + m1_vals = [] + m1_err_vals = [] + c1_vals = [] + c1_err_vals = [] + m2_vals = [] + m2_err_vals = [] + c2_vals = [] + c2_err_vals = [] + + for n in n_tiles_list: + res = cumulative[str(n)] + m1_vals.append(res["m1"]) + m1_err_vals.append(res["m1_err"]) + c1_vals.append(res["c1"]) + c1_err_vals.append(res["c1_err"]) + m2_vals.append(res["m2"]) + m2_err_vals.append(res["m2_err"]) + c2_vals.append(res["c2"]) + c2_err_vals.append(res["c2_err"]) + + n_tiles_str = ( + f"n_tiles = {n_tiles_list}" + if len(n_tiles_list) > 1 + else f"n_tiles = {n_tiles_list[0]}" + ) + + # Plot 1: m and c with error bars + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) + fig.suptitle(f"m and c convergence ({n_tiles_str})", fontsize=12) + + ax1.errorbar( + n_tiles_list, m1_vals, yerr=m1_err_vals, fmt="o-", label="m1", capsize=5 + ) + ax1.errorbar( + n_tiles_list, m2_vals, yerr=m2_err_vals, fmt="s-", label="m2", capsize=5 + ) + ax1.axhline(0, color="k", linestyle="--", alpha=0.3) + ax1.set_xlabel("Number of tiles") + ax1.set_ylabel("Multiplicative bias m") + ax1.legend() + ax1.grid(True, alpha=0.3) + + ax2.errorbar( + n_tiles_list, c1_vals, yerr=c1_err_vals, fmt="o-", label="c1", capsize=5 + ) + ax2.errorbar( + n_tiles_list, c2_vals, yerr=c2_err_vals, fmt="s-", label="c2", capsize=5 + ) + ax2.axhline(0, color="k", linestyle="--", alpha=0.3) + ax2.set_xlabel("Number of tiles") + ax2.set_ylabel("Additive bias c") + ax2.legend() + ax2.grid(True, alpha=0.3) + + plt.tight_layout() + plot1_path = os.path.join(diagnostics_dir, "mbias_convergence.png") + plt.savefig(plot1_path, dpi=150) + plt.close() + print(f"Saved convergence plot to {plot1_path}") + + # Plot 2: error bars only + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) + fig.suptitle(f"Error convergence ({n_tiles_str})", fontsize=12) + + ax1.errorbar( + n_tiles_list, + [0] * len(n_tiles_list), + yerr=m1_err_vals, + fmt="o-", + label="m1 error", + capsize=5, + alpha=0.7, + ) + ax1.errorbar( + n_tiles_list, + [0] * len(n_tiles_list), + yerr=m2_err_vals, + fmt="s-", + label="m2 error", + capsize=5, + alpha=0.7, + ) + ax1.set_xlabel("Number of tiles") + ax1.set_ylabel("Multiplicative bias error") + ax1.legend() + ax1.grid(True, alpha=0.3) + ax1.set_ylim(bottom=0) + + ax2.errorbar( + n_tiles_list, + [0] * len(n_tiles_list), + yerr=c1_err_vals, + fmt="o-", + label="c1 error", + capsize=5, + alpha=0.7, + ) + ax2.errorbar( + n_tiles_list, + [0] * len(n_tiles_list), + yerr=c2_err_vals, + fmt="s-", + label="c2 error", + capsize=5, + alpha=0.7, + ) + ax2.set_xlabel("Number of tiles") + ax2.set_ylabel("Additive bias error") + ax2.legend() + ax2.grid(True, alpha=0.3) + ax2.set_ylim(bottom=0) + + plt.tight_layout() + plot2_path = os.path.join(diagnostics_dir, "mbias_errors.png") + plt.savefig(plot2_path, dpi=150) + plt.close() + print(f"Saved errors plot to {plot2_path}") + + +def main(): + args = parse_args() + + with open(args.config) as f: + config = yaml.safe_load(f) + + print(f"Config: {args.config}") + print(f"Grids : {config['grids_dir']}") + print(f"Run : grid_{config['num']}") + print(f"g_in : ±{config['shear_amplitude']}") + print() + + # Auto-detect n_tiles if --cumulative + if args.cumulative and not args.n_tiles: + n_tiles = get_n_tiles(config["grids_dir"], config["num"]) + if n_tiles: + args.n_tiles = n_tiles + print(f"Auto-detected {n_tiles} tiles") + + mb = ImageSimMBias(config) + + print("Loading catalogues...") + mb.load_catalogs(verbose=args.verbose) + + # ``run`` returns a document with the primary scheme's m/c mirrored at the + # top level plus a per-scheme ``weights`` block. Cast the whole tree to + # plain Python floats so the YAML/text output is human-readable (raw numpy + # scalars serialise as !!python/object binary). + results = to_python(mb.run(verbose=True)) + + print() + print("=" * 40) + print(" Results") + print("=" * 40) + for scheme, res in results["weights"].items(): + print(f" weights: {scheme}") + print(f" m1 = {res['m1']:+.4f} +-{res['m1_err']:.4f}") + print(f" c1 = {res['c1']:+.4f} +-{res['c1_err']:.4f}") + print(f" m2 = {res['m2']:+.4f} +-{res['m2_err']:.4f}") + print(f" c2 = {res['c2']:+.4f} +-{res['c2_err']:.4f}") + print("=" * 40) + + # Cumulative tracking + if args.cumulative: + results_dir = config.get( + "diagnostics_dir", config.get("results_dir", "results") + ) + os.makedirs(results_dir, exist_ok=True) + else: + results_dir = None + + # Output path: in results dir if cumulative, else from config or current dir + if results_dir: + out_path = os.path.join(results_dir, "m_bias_results.yaml") + else: + out_path = config.get("output_path", "m_bias_results.yaml") + + # A result file describes itself: the provenance block the rule assembled + # (manifest hash, both repos' branch+commit, container sif + GHCR revision) + # rides verbatim from the config into the output yaml. It is appended as a + # separate top-level key, so the numeric m/c fields serialise byte-for-byte + # as before -- the reproduction gate sees only the added `provenance:` block. + output = dict(results) + if "provenance" in config: + output["provenance"] = config["provenance"] + + with open(out_path, "w") as f: + yaml.dump(output, f, default_flow_style=False) + print(f"Results written to {out_path}") + + # Also write to text file for readability + if results_dir: + txt_path = os.path.join(results_dir, "m_bias_results.txt") + else: + txt_path = config.get("output_path", "m_bias_results.yaml").replace( + ".yaml", ".txt" + ) + + with open(txt_path, "w") as f: + f.write("Multiplicative and additive shear bias from image simulations\n") + f.write("=" * 60 + "\n") + for scheme, res in results["weights"].items(): + f.write(f"\nweights: {scheme}\n") + f.write(f" m1 = {res['m1']:+.6f} ± {res['m1_err']:.6f}\n") + f.write(f" c1 = {res['c1']:+.6f} ± {res['c1_err']:.6f}\n") + f.write(f" m2 = {res['m2']:+.6f} ± {res['m2_err']:.6f}\n") + f.write(f" c2 = {res['c2']:+.6f} ± {res['c2_err']:.6f}\n") + f.write( + "\nErrors computed via bootstrap resampling " + f"(n={config['n_bootstrap']} resamples)\n" + ) + print(f"Results written to {txt_path}") + + if args.cumulative: + cumulative_path = os.path.join(results_dir, "mbias_cumulative.yaml") + if args.n_tiles: + added = update_cumulative_file(cumulative_path, args.n_tiles, results) + verb = "Added" if added else "Overwrote" + print(f"\n{verb} n_tiles={args.n_tiles} in {cumulative_path}") + # Regenerate plots after every update (an overwrite can shift the + # curve, so the plots must track it -- not just fresh additions). + try: + plot_convergence(cumulative_path, results_dir) + except Exception as e: + print( + f"Warning: could not generate convergence plots: {e}", + file=sys.stderr, + ) + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/diagnostics_image_sims.py b/scripts/diagnostics_image_sims.py new file mode 100644 index 00000000..699d5c49 --- /dev/null +++ b/scripts/diagnostics_image_sims.py @@ -0,0 +1,227 @@ +#!/usr/bin/env python +"""Per-sim diagnostics for image simulation catalogues. + +For each requested grid catalogue, produces: + - footprint (RA/Dec scatter) + - ellipticity histograms (e1, e2) + - weight histogram + - response matrix element histograms (R_g11, R_g22, R_g12, R_g21) + - PSF leakage scatter (e1 vs e1_PSF, e2 vs e2_PSF) + - additive bias (weighted mean e1, e2) + +Shares the estimator's config schema (``sp_validation.image_sims``): the same +``grids_dir`` / ``num`` / ``catalog_name`` keys, the ``branches`` list (the sim +map, not a hard-coded five), and the same ``w_col`` weight semantics -- a column +name, or ``null`` for unit weights. Reading ``w_col`` (rather than hard-coding +``w_des``) means the diagnostics never KeyError on a catalogue that lacks the +weight column, and they weight exactly as the m-bias run they accompany. + +Usage: + diagnostics_image_sims.py -c config.yaml [-v] +""" + +import argparse +import os +import sys + +import matplotlib +import numpy as np +import yaml + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from astropy.io import fits + +# Conventional campaign layout, used only when the config carries no branch map +# -- the same fallback the estimator uses. +_DEFAULT_BRANCHES = ["1z2z", "1m2z", "1p2z", "1z2m", "1z2p"] + + +def load(path): + with fits.open(path) as hdul: + return {col.name: hdul[1].data[col.name].copy() for col in hdul[1].columns} + + +def weights(cat, w_col): + """Per-object weights: the ``w_col`` column, or unit weights when null. + + Mirrors the estimator's ``w_col`` contract (image_sims._load_cat): ``None`` + -> every object unit weight (the no-weighting mode, #227). Reading it here + means the diagnostics never KeyError when the weight column is absent. + """ + return cat[w_col].copy() if w_col else np.ones(len(cat["RA"])) + + +def parse_args(): + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("-c", "--config", required=True) + p.add_argument("-v", "--verbose", action="store_true") + return p.parse_args() + + +def savefig(fig, out_dir, name): + path = f"{out_dir}/{name}.png" + fig.savefig(path, dpi=150, bbox_inches="tight") + plt.close(fig) + return path + + +def plot_footprints(cats, colors, out_dir): + fig, ax = plt.subplots(figsize=(8, 6)) + for name, d in cats.items(): + ax.scatter(d["RA"], d["Dec"], s=1, alpha=0.4, label=name, color=colors[name]) + ax.set_xlabel("RA [deg]") + ax.set_ylabel("Dec [deg]") + ax.legend(markerscale=5) + ax.set_title("Footprint") + return savefig(fig, out_dir, "footprint") + + +def plot_ellipticity(cats, colors, w_col, out_dir, nbins=100): + fig, axs = plt.subplots(1, 2, figsize=(14, 5)) + bins = np.linspace(-1.0, 1.0, nbins + 1) + for name, d in cats.items(): + w = weights(d, w_col) + for ax, col, label in zip(axs, ["e1", "e2"], [r"$e_1$", r"$e_2$"]): + ax.hist( + d[col], + bins=bins, + density=True, + weights=w, + histtype="step", + label=name, + color=colors[name], + ) + for ax, label in zip(axs, [r"$e_1$", r"$e_2$"]): + ax.set_xlabel(label) + ax.set_ylabel("normalised count") + ax.legend(fontsize=7) + wlabel = w_col if w_col else "unit" + fig.suptitle(f"Ellipticity histograms ({wlabel} weighted)") + return savefig(fig, out_dir, "ellipticity_hist") + + +def plot_weights(cats, colors, w_col, out_dir, nbins=50): + fig, ax = plt.subplots(figsize=(8, 5)) + for name, d in cats.items(): + ax.hist( + weights(d, w_col), + bins=nbins, + density=True, + histtype="step", + label=name, + color=colors[name], + ) + wlabel = w_col if w_col else "unit" + ax.set_xlabel(wlabel) + ax.set_ylabel("normalised count") + ax.legend() + ax.set_title("Weight distribution") + return savefig(fig, out_dir, "weight_hist") + + +def plot_response(cats, colors, out_dir, nbins=50): + cols = ["R_g11", "R_g22", "R_g12", "R_g21"] + fig, axs = plt.subplots(2, 2, figsize=(12, 10)) + for ax, col in zip(axs.flat, cols): + for name, d in cats.items(): + ax.hist( + d[col], + bins=nbins, + range=(-1, 2), + density=True, + histtype="step", + label=name, + color=colors[name], + ) + ax.set_xlim(-1, 2) + ax.set_xlabel(col) + ax.set_ylabel("normalised count") + ax.legend(fontsize=7) + fig.suptitle("Response matrix elements") + fig.tight_layout() + return savefig(fig, out_dir, "response_hist") + + +def plot_psf_leakage(cats, colors, out_dir): + fig, axs = plt.subplots(1, 2, figsize=(14, 5)) + for name, d in cats.items(): + for ax, eg, ep, label in zip( + axs, + ["e1", "e2"], + ["e1_PSF", "e2_PSF"], + [r"$e_1$", r"$e_2$"], + ): + ax.scatter(d[ep], d[eg], s=1, alpha=0.3, label=name, color=colors[name]) + for ax, xlab, ylab in zip( + axs, [r"$e_1^{\rm PSF}$", r"$e_2^{\rm PSF}$"], [r"$e_1$", r"$e_2$"] + ): + ax.set_xlabel(xlab) + ax.set_ylabel(ylab) + ax.legend(markerscale=5, fontsize=7) + fig.suptitle("Object-wise PSF leakage") + return savefig(fig, out_dir, "psf_leakage") + + +def calculate_additive_bias(cats, w_col, verbose=True): + print("\n--- Additive bias (weighted mean ellipticity) ---") + results = {} + for name, d in cats.items(): + w = weights(d, w_col) + c1 = np.average(d["e1"], weights=w) + c2 = np.average(d["e2"], weights=w) + results[name] = (c1, c2) + if verbose: + print(f" {name}: c1 = {c1:+.5f} c2 = {c2:+.5f}") + return results + + +def main(): + args = parse_args() + with open(args.config) as f: + config = yaml.safe_load(f) + + # Same config schema as the estimator: grids_dir (not base), num, + # catalog_name, the branch map, and the w_col weight contract. + grids_dir = config["grids_dir"] + num = config["num"] + cat_name = config.get("catalog_name", "shape_catalog_cut_ngmix.fits") + branches = list(config.get("branches", _DEFAULT_BRANCHES)) + w_col = config["w_col"] # required, like the estimator; null -> unit weights + out_dir = config.get("diagnostics_dir", f"{grids_dir}/diagnostics") + + # Colour per branch from a palette, so any branch list plots (no hard-coded + # five-branch colour map). + palette = plt.get_cmap("tab10") + colors = {name: palette(i % 10) for i, name in enumerate(branches)} + + os.makedirs(out_dir, exist_ok=True) + + print(f"Loading catalogues from {grids_dir}...") + cats = {} + for name in branches: + path = f"{grids_dir}/{name}_grid_{num}/{cat_name}" + if not os.path.exists(path): + print(f" WARNING: {path} not found, skipping") + continue + cats[name] = load(path) + if args.verbose: + print(f" {name}: {len(cats[name]['RA'])} objects") + + if not cats: + print("No catalogues found, exiting.") + return 1 + + print(f"\nSaving plots to {out_dir}/") + print(f" footprint -> {plot_footprints(cats, colors, out_dir)}") + print(f" ellipticity -> {plot_ellipticity(cats, colors, w_col, out_dir)}") + print(f" weights -> {plot_weights(cats, colors, w_col, out_dir)}") + print(f" response -> {plot_response(cats, colors, out_dir)}") + print(f" PSF leakage -> {plot_psf_leakage(cats, colors, out_dir)}") + calculate_additive_bias(cats, w_col, verbose=True) + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/sp_validation/calibration.py b/src/sp_validation/calibration.py index 9032bc63..e37fc4d0 100644 --- a/src/sp_validation/calibration.py +++ b/src/sp_validation/calibration.py @@ -56,7 +56,7 @@ def get_calibrated_quantities(gal_metacal): return g_corr, g_uncorr, w, mask -def get_calibrated_m_c(gal_metacal): +def get_calibrated_m_c(gal_metacal, additive_correction=True): """Get Calibrated C. Return catalogue quantities for objects calibrated for multiplicative and @@ -66,6 +66,11 @@ def get_calibrated_m_c(gal_metacal): ---------- gal_metacal : dict galaxy metacalibration catalogue + additive_correction : bool, optional, default=True + if False, do not subtract the additive bias c from the shear + estimates; use for constant-shear image sims, where the mean + shear is the signal (see issue #226). c and c_err are still + computed and returned Returns ------- @@ -99,10 +104,11 @@ def get_calibrated_m_c(gal_metacal): c_err[comp] = np.std(g_uncorr[comp]) # Shear estimate corrected for additive bias - g_corr_mc = np.zeros_like(g_corr) - c_corr = np.linalg.inv(gal_metacal.R).dot(c) - for comp in (0, 1): - g_corr_mc[comp] = g_corr[comp] - c_corr[comp] + g_corr_mc = np.copy(g_corr) + if additive_correction: + c_corr = np.linalg.inv(gal_metacal.R).dot(c) + for comp in (0, 1): + g_corr_mc[comp] = g_corr[comp] - c_corr[comp] return g_corr_mc, g_uncorr, w, mask_metacal, c, c_err diff --git a/src/sp_validation/catalog.py b/src/sp_validation/catalog.py index b1b175d9..e39d2076 100644 --- a/src/sp_validation/catalog.py +++ b/src/sp_validation/catalog.py @@ -12,6 +12,7 @@ """ import getpass +import os import h5py import numpy as np @@ -308,6 +309,35 @@ def match_subsample( return ra, dec, g +def match_catalogs_radec(ra1, dec1, ra2, dec2, thresh_deg=0.0002): + """Match two catalogues by RA/Dec. + + Match each object in catalogue 2 to the nearest in catalogue 1 + within a threshold. + + Parameters + ---------- + ra1, dec1 : array_like + coordinates of reference catalogue [deg] + ra2, dec2 : array_like + coordinates of catalogue to match [deg] + thresh_deg : float, optional + maximum separation [deg], default 0.0002 + + Returns + ------- + idx1 : ndarray of int + indices into catalogue 1 of matched objects + idx2 : ndarray of int + indices into catalogue 2 of matched objects + """ + coord1 = coords.SkyCoord(ra=ra1 * u.degree, dec=dec1 * u.degree) + coord2 = coords.SkyCoord(ra=ra2 * u.degree, dec=dec2 * u.degree) + idx1, sep, _ = coord2.match_to_catalog_sky(coord1) + mask = sep.deg < thresh_deg + return idx1[mask], np.where(mask)[0] + + def match_stars2(ra_gal, dec_gal, ra_star, dec_star, thresh=0.0002): """Add docstring. @@ -542,55 +572,91 @@ def write_shape_catalog( ) ) - # Write columns to FITS file - cols = [] - for col, _ in col_info_arr: - cols.append(col) - table_hdu = fits.BinTableHDU.from_columns(cols) - - # Add human-readable descriptions - for idx, col_info in enumerate(col_info_arr): - table_hdu.header[f"TTYPE{idx + 1}"] = ( - col_info[0].name, - col_info[1], - ) + ext = os.path.splitext(output_path)[1].lower() - # Primary HDU with information in header - primary_header = fits.Header() + if ext in (".hdf5", ".hdf", ".h5"): + # Build flat list of (name, 1d-array) pairs, splitting 2D columns + fields = [] + for col, _ in col_info_arr: + arr = np.asarray(col.array) + if arr.ndim == 2: + for idx in range(arr.shape[1]): + fields.append((f"{col.name}_{idx}", arr[:, idx])) + else: + fields.append((col.name, arr)) + + # Build structured numpy array and write as single "data" dataset + dtype = np.dtype([(name, arr.dtype) for name, arr in fields]) + structured = np.empty(len(fields[0][1]), dtype=dtype) + for name, arr in fields: + structured[name] = arr + + with h5py.File(output_path, "w") as f: + f.create_dataset("data", data=structured) + if add_header: + for key, val in add_header.items(): + f.attrs[key] = str(val) + if all(v is not None for v in (R, R_shear, R_select, c)): + f.attrs["R"] = R + f.attrs["R_shear"] = R_shear + f.attrs["R_select"] = R_select + f.attrs["c"] = c + if c_err is not None: + f.attrs["c1_err"] = c_err[0] + f.attrs["c2_err"] = c_err[1] + if sigma_epsilon is not None: + f.attrs["sig_eps"] = sigma_epsilon + if alpha_leakage is not None: + f.attrs["alpha"] = alpha_leakage - if add_header: - primary_header.update(add_header) + else: + # Write columns to FITS file + cols = [col for col, _ in col_info_arr] + table_hdu = fits.BinTableHDU.from_columns(cols) + + # Add human-readable descriptions + for idx, col_info in enumerate(col_info_arr): + table_hdu.header[f"TTYPE{idx + 1}"] = ( + col_info[0].name, + col_info[1], + ) - primary_header = cat.write_header_info_sp( - primary_header, - software_name="sp_validation", - software_version=__version__, - author=getpass.getuser(), - ) + # Primary HDU with information in header + primary_header = fits.Header() - if all(v is not None for v in (R, R_shear, R_select, c)): - cat.add_shear_bias_to_header(primary_header, R, R_shear, R_select, c) - if c_err is not None: - primary_header["c1_err"] = (c_err[0], "Standard deviation of c_1") - primary_header["c2_err"] = (c_err[1], "Standard deviation of c_2") + if add_header: + primary_header.update(add_header) - primary_header["w"] = "DES weight" + primary_header = cat.write_header_info_sp( + primary_header, + software_name="sp_validation", + software_version=__version__, + author=getpass.getuser(), + ) - if sigma_epsilon is not None: - primary_header["sig_eps"] = (sigma_epsilon, "Shape noise RMS") + if all(v is not None for v in (R, R_shear, R_select, c)): + cat.add_shear_bias_to_header(primary_header, R, R_shear, R_select, c) + if c_err is not None: + primary_header["c1_err"] = (c_err[0], "Standard deviation of c_1") + primary_header["c2_err"] = (c_err[1], "Standard deviation of c_2") - if alpha_leakage: - primary_header["alpha"] = ( - alpha_leakage, - "Mean scale-dependent PSF leakage", - ) + primary_header["w"] = "DES weight" + + if sigma_epsilon is not None: + primary_header["sig_eps"] = (sigma_epsilon, "Shape noise RMS") + + if alpha_leakage: + primary_header["alpha"] = ( + alpha_leakage, + "Mean scale-dependent PSF leakage", + ) - primary_hdu = fits.PrimaryHDU(header=primary_header) + primary_hdu = fits.PrimaryHDU(header=primary_header) - # Final file - hdu_list = fits.HDUList([primary_hdu, table_hdu]) + # Final file + hdu_list = fits.HDUList([primary_hdu, table_hdu]) - hdu_list.writeto(output_path, overwrite=True) + hdu_list.writeto(output_path, overwrite=True) def write_galaxy_cat(output_path, ra, dec, tile_id): diff --git a/src/sp_validation/catalog_builders.py b/src/sp_validation/catalog_builders.py index 1375ef8d..78dc7838 100644 --- a/src/sp_validation/catalog_builders.py +++ b/src/sp_validation/catalog_builders.py @@ -1124,6 +1124,22 @@ def read_cat(self, load_into_memory=False): fpath = self._params["input_path"] verbose = self._params["verbose"] + # Image-simulation path: a single per-run comprehensive catalogue in + # FITS, not the joined multi-patch HDF5 the data path builds. Read the + # FITS table directly into memory; there is no separate data_ext group. + extension = os.path.splitext(fpath)[1] + if extension == ".fits": + if verbose: + print(f"Reading FITS file {fpath}, HDU 1...") + dat = fits.getdata(fpath, 1) + dat_ext = None + if verbose: + print( + f"Found {len(dat)} (~{format.millify(len(dat))}) objects" + + " in catalogue" + ) + return dat, dat_ext + if verbose: print(f"Reading HDF5 file {fpath}...") diff --git a/src/sp_validation/image_sims.py b/src/sp_validation/image_sims.py new file mode 100644 index 00000000..e259109d --- /dev/null +++ b/src/sp_validation/image_sims.py @@ -0,0 +1,324 @@ +"""IMAGE_SIMS. + +:Description: Multiplicative and additive shear bias from image simulations. + +:Author: Martin Kilbinger + +""" + +import numpy as np +from astropy.io import fits + +from sp_validation.catalog import match_catalogs_radec + +# Conventional campaign layout, used only when the config carries no branch map +# (e.g. the synthetic-recovery tests). In a workflow run the branches and pairs +# come from manifest.yaml via the m_bias config; nothing about the injected +# shear is hard-coded on the estimator's side. +_DEFAULT_BRANCHES = ["1z2z", "1p2z", "1m2z", "1z2p", "1z2m"] +_DEFAULT_PAIRS = [ + ("1p2z", "1m2z", 0), # g1 component, index 0 → e1 + ("1z2p", "1z2m", 1), # g2 component, index 1 → e2 +] + + +# Weight-scheme name that means "no weighting": every object gets unit weight. +# ``None`` (from a YAML ``null``) is accepted as an alias, so the fiducial +# unweighted primary scheme can be written either ``none`` or ``null``. +_UNWEIGHTED = "none" + + +def _is_unweighted(scheme): + """True for the unit-weight scheme (``"none"`` or ``None``).""" + return scheme is None or scheme == _UNWEIGHTED + + +def _load_cat(path, w_cols): + """Load RA, Dec, ellipticities and per-scheme weights from a FITS catalogue. + + Reads the ``e1``/``e2`` columns, which the calibration stage writes as the + *calibrated* shear estimate ``g = R^-1 g_uncal - c`` (metacal response and + additive-bias corrected) -- not the raw ``e1_uncal``/``e2_uncal`` columns + that sit alongside them in the same catalogue. The bias this estimator + measures is therefore the *residual* m/c left after the chain's own metacal + calibration, not the raw pre-calibration bias. + + ``w_cols`` is the list of weight schemes to load. The scheme ``"none"`` + (equivalently a ``None``/``null`` entry) gives every object unit weight -- + the no-weighting mode for m-bias runs (#227: shape weights are excluded from + sim calibration); any other entry is read as a FITS column name. The weights + come back as a dict keyed by scheme so one catalogue load serves every + scheme in a multi-weight run. + """ + with fits.open(path) as hdul: + data = hdul[1].data + cat = { + "ra": data["RA"].copy(), + "dec": data["Dec"].copy(), + "e1": data["e1"].copy(), + "e2": data["e2"].copy(), + "w": {}, + } + for scheme in w_cols: + cat["w"][scheme] = ( + np.ones(len(cat["ra"])) + if _is_unweighted(scheme) + else data[scheme].copy() + ) + return cat + + +class ImageSimMBias: + """Compute multiplicative and additive shear bias from image simulations. + + The estimator consumes the *calibrated* ``e1``/``e2`` columns (the metacal + response- and additive-bias-corrected shear ``g = R^-1 g_uncal - c``), so + the headline m/c is the **residual** bias remaining after the chain's own + metacal calibration, not the raw pre-calibration bias. + + Parameters + ---------- + config : dict + Configuration dictionary with keys: + - grids_dir : str, path to the grids directory + - num : int, run number (e.g. 2 for *_grid_2) + - catalog_name : str, filename of the cut catalogue + (default 'shape_catalog_cut_ngmix.fits') + - shear_amplitude : float, input shear |g| (from manifest.yaml) + - branches : list of str, branch names in load order (incl. the + unsheared reference); defaults to the conventional 5-branch layout + - pairs : list of dicts {plus, minus, component}, the +/- sheared + branch pairing per component; defaults to the conventional pairs + - match_radius_deg : float, matching radius in degrees (required) + - pair_match : bool, match objects between the +g and -g sheared + catalogues (required); if False, use all objects of each + catalogue (the paired per-object cancellation is then unavailable) + - w_cols : list of str, weight schemes to compute in one run + (required); ``"none"`` (or a ``null`` entry) means unit weights, + any other entry is a FITS column name. The **first** entry is the + primary result surfaced at the top level of ``run()``'s output. Our + fiducial run leads with the unweighted scheme (``["none", ...]``), + per the #227 verdict that shape weights are excluded from sim + calibration; the unweighted m also avoids the ``cov(w, e)`` residual + weighted estimators carry on constant-shear sims. + - w_col : str or None, *deprecated* single weight scheme; accepted for + back-compat and used as ``[w_col]`` only when ``w_cols`` is absent. + - n_bootstrap : int, number of bootstrap resamples for errors (required) + - bootstrap_seed : int, seed for the per-pair bootstrap RNG (required); + makes the bootstrap errors bit-reproducible. The resample indices are + drawn once per pair and shared across every weight scheme, so the + schemes differ only in their weighting, never in their draws. + + The science knobs (``match_radius_deg``, ``pair_match``, ``w_cols``, + ``n_bootstrap``, ``bootstrap_seed``) are read with no in-code default: a + missing one is a config bug and raises ``KeyError`` at construction, per the + fail-fast contract (the workflow emits every one into the m_bias config). + The lone exception is the deprecated ``w_col``, which is honoured as a + fallback so pre-``w_cols`` configs still run. + """ + + def __init__(self, config): + self.cfg = config + self.g_in = config["shear_amplitude"] + self.thresh = config["match_radius_deg"] + self.pair_match = config["pair_match"] + # ``w_cols`` is the required science key. A pre-``w_cols`` config that + # still carries the deprecated scalar ``w_col`` is honoured as a + # single-scheme run; only a config with neither raises (fail-fast). + if "w_cols" in config: + w_cols = config["w_cols"] + else: + w_cols = [config["w_col"]] + # Normalise a ``None``/``null`` entry to the canonical "none" name so + # results key off a string; downstream still treats it as unit weights. + self.w_cols = [_UNWEIGHTED if _is_unweighted(w) else str(w) for w in w_cols] + self.n_boot = config["n_bootstrap"] + self.boot_seed = config["bootstrap_seed"] + # Branch list and pairing come from the manifest-derived config + # (``branches`` / ``pairs``); fall back to the conventional layout only + # when neither is given. ``branches`` fixes the catalogue load order; + # ``pairs`` fixes which sims difference into which component. + self.sim_names = list(config.get("branches", _DEFAULT_BRANCHES)) + if config.get("pairs"): + self.pairs = [ + (p["plus"], p["minus"], p["component"]) for p in config["pairs"] + ] + else: + self.pairs = list(_DEFAULT_PAIRS) + self.cats = {} + + def load_catalogs(self, verbose=True): + """Load the 5 sheared and reference catalogues.""" + grids_dir = self.cfg["grids_dir"] + num = self.cfg["num"] + cat_name = self.cfg.get("catalog_name", "shape_catalog_cut_ngmix.fits") + # ``sim_names`` (incl. the unsheared reference) comes from the config's + # branch map. The +g/-g pool estimator pairs the sheared sims directly; + # the reference is loaded for completeness and null-test diagnostics. + for name in self.sim_names: + path = f"{grids_dir}/{name}_grid_{num}/{cat_name}" + if verbose: + print(f" Loading {path}") + self.cats[name] = _load_cat(path, self.w_cols) + if verbose: + print(f" {len(self.cats[name]['ra'])} objects") + + def print_mean_ellipticities(self): + """Print the mean e1, e2 for each catalogue and weight scheme, as a check. + + The unweighted scheme (``"none"``) gives the plain unweighted means. + """ + for scheme in self.w_cols: + print(f"\nMean ellipticities (all objects, weights: {scheme}):") + for name, cat in self.cats.items(): + mean_e1 = np.average(cat["e1"], weights=cat["w"][scheme]) + mean_e2 = np.average(cat["e2"], weights=cat["w"][scheme]) + print(f" {name}: = {mean_e1:+.5f} = {mean_e2:+.5f}") + + def _m_c_pair(self, name_p, name_m, comp, verbose=True): + """Compute m and c for one shear pair and component (0=g1, 1=g2). + + Paired ("pool") estimator. The +g and -g simulations inject opposite + input shear on the *same* galaxies, so matching them directly by + RA/Dec yields a one-to-one correspondence. Differencing the two + ellipticities per object, + + m = <(e_+ - e_-) / (2 g_in) - 1> , c = <(e_+ + e_-) / 2> , + + cancels the intrinsic shape (sigma_e ~ 0.3) object-by-object in the + multiplicative term, leaving only measurement noise -- so sigma(m) + shrinks by ~sigma_e/sigma_meas relative to differencing two + independent means. (The additive term c is a *sum*, so intrinsic + shape does not cancel there and its error stays shape-noise limited.) + + With ``pair_match=False`` the +g and -g sims are *not* matched: every + object of each catalogue is used, so the per-object cancellation is + lost and m, c fall back to differencing/summing the two independent + weighted means. The paired bootstrap likewise cannot be applied (the + two arrays generally have different lengths), so each side is resampled + independently per replicate. + """ + e_key = f"e{comp + 1}" + + if self.pair_match: + # Match the +g and -g sims to each other: same galaxies, opposite + # shear. This is a nearest-neighbour match within `thresh`, not a + # strict bijection -- on grid sims galaxies are well separated so + # pairs are effectively 1:1 (verified ~99% co-located to <0.05" on + # SKiLLS grid_1); on denser fields a small fraction could share a + # +g partner and dilute the cancellation. + idx_p, idx_m = match_catalogs_radec( + self.cats[name_p]["ra"], + self.cats[name_p]["dec"], + self.cats[name_m]["ra"], + self.cats[name_m]["dec"], + thresh_deg=self.thresh, + ) + if verbose: + print(f" {name_p} <-> {name_m}: {len(idx_p)} paired objects") + else: + idx_p = slice(None) + idx_m = slice(None) + if verbose: + print( + f" no pair-matching: {name_p}: {len(self.cats[name_p][e_key])}" + f" | {name_m}: {len(self.cats[name_m][e_key])} objects" + ) + + e_p = self.cats[name_p][e_key][idx_p] + e_m = self.cats[name_m][e_key][idx_m] + + # Draw the bootstrap resample indices *once*, before the weight-scheme + # loop, and reuse them for every scheme -- the schemes then differ only + # in their weighting, never in their draws (so a scheme comparison is a + # clean weighting comparison). Pre-drawing the full ``(n_boot, n)`` block + # in one call is bit-identical to drawing ``rng.integers(0, n, n)`` once + # per replicate (numpy fills the block row-major), so the numbers match a + # single-scheme, per-iteration bootstrap to the last bit. + rng = np.random.default_rng(seed=self.boot_seed) + if self.pair_match: + n = len(e_p) + ib = rng.integers(0, n, (self.n_boot, n)) + else: + n_p, n_m = len(e_p), len(e_m) + ib_p = rng.integers(0, n_p, (self.n_boot, n_p)) + ib_m = rng.integers(0, n_m, (self.n_boot, n_m)) + + res = {} + for scheme in self.w_cols: + w_p = self.cats[name_p]["w"][scheme][idx_p] + w_m = self.cats[name_m]["w"][scheme][idx_m] + m_boot = np.empty(self.n_boot) + c_boot = np.empty(self.n_boot) + + if self.pair_match: + # Per-object shear-differenced (-> m) and summed (-> c) + # ellipticity, with a symmetric per-pair weight. + w = 0.5 * (w_p + w_m) + d = (e_p - e_m) / (2 * self.g_in) - 1 + s = (e_p + e_m) / 2 + + m = np.average(d, weights=w) + c = np.average(s, weights=w) + + # Paired bootstrap: the same object draw is applied to both + # sims, so the per-object cancellation in `d` is preserved in + # the error estimate. + for i in range(self.n_boot): + m_boot[i] = np.average(d[ib[i]], weights=w[ib[i]]) + c_boot[i] = np.average(s[ib[i]], weights=w[ib[i]]) + else: + # No matching: difference/sum the two independent weighted means. + mean_ep = np.average(e_p, weights=w_p) + mean_em = np.average(e_m, weights=w_m) + + m = (mean_ep - mean_em) / (2 * self.g_in) - 1 + c = (mean_ep + mean_em) / 2 + + # Unpaired bootstrap: the +g and -g arrays generally differ in + # length, so each side is resampled independently per replicate. + for i in range(self.n_boot): + ep_b = np.average(e_p[ib_p[i]], weights=w_p[ib_p[i]]) + em_b = np.average(e_m[ib_m[i]], weights=w_m[ib_m[i]]) + m_boot[i] = (ep_b - em_b) / (2 * self.g_in) - 1 + c_boot[i] = (ep_b + em_b) / 2 + + res[scheme] = (m, np.std(m_boot), c, np.std(c_boot)) + + return res + + def run(self, verbose=True): + """Compute m and c for both shear components and every weight scheme. + + Returns + ------- + dict + ``results["weights"][scheme]`` holds ``m1, m1_err, c1, c1_err, + m2, m2_err, c2, c2_err`` for each weight scheme. The primary + (first) scheme's keys are also mirrored at the top level, so a + reader that wants the headline m/c never has to know the scheme + name. + """ + results = {"weights": {scheme: {} for scheme in self.w_cols}} + for name_p, name_m, comp in self.pairs: + label = f"g{comp + 1}" + if verbose: + print(f"\n--- {label}: {name_p} / {name_m} ---") + res = self._m_c_pair(name_p, name_m, comp, verbose=verbose) + for scheme, (m, m_err, c, c_err) in res.items(): + w = results["weights"][scheme] + w[f"m{comp + 1}"] = m + w[f"m{comp + 1}_err"] = m_err + w[f"c{comp + 1}"] = c + w[f"c{comp + 1}_err"] = c_err + if verbose: + print( + f" [{scheme}] m{comp + 1} = {m:.4f} ± {m_err:.4f}" + f" c{comp + 1} = {c:.4f} ± {c_err:.4f}" + ) + + # Mirror the primary (first) scheme's m/c at the top level: the headline + # result reads out without knowing the scheme name, and a downstream + # gate keyed on the old flat keys still finds them. + results.update(results["weights"][self.w_cols[0]]) + return results diff --git a/src/sp_validation/tests/test_config_paths_exist.py b/src/sp_validation/tests/test_config_paths_exist.py index 9ae74ca0..2db03bc1 100644 --- a/src/sp_validation/tests/test_config_paths_exist.py +++ b/src/sp_validation/tests/test_config_paths_exist.py @@ -23,7 +23,14 @@ "catalog", "catalogue", ) -NON_PATH_KEYS = ("extra_output",) +# Keys whose values are never filesystem paths to check. ``extra_output`` is a +# flag, not a path. ``why``/``replace``/``with``/``drop`` are the declaration +# keys of a mask *overlay* (config/calibration/*.overlay.yaml): ``why`` is +# rationale prose and ``replace``/``with``/``drop`` are verbatim blocks of base +# config text -- content, not paths -- so the path walker must not treat them as +# files to stat. (The overlay's one real path, ``base:``, is deliberately not +# listed, so it is still validated.) +NON_PATH_KEYS = ("extra_output", "why", "replace", "with", "drop") PATH_PREFIX_KEYS = ("nz.dndz.path",) TEXT_SUFFIXES = ( ".fits", diff --git a/src/sp_validation/tests/test_image_sims.py b/src/sp_validation/tests/test_image_sims.py new file mode 100644 index 00000000..ef3e3b14 --- /dev/null +++ b/src/sp_validation/tests/test_image_sims.py @@ -0,0 +1,253 @@ +"""UNIT TESTS FOR THE IMAGE-SIMULATION m/c ESTIMATOR. + +Exercise ``sp_validation.image_sims.ImageSimMBias`` -- the multiplicative and +additive shear-bias estimator used by the image-simulation workflow -- and the +``sp_validation.catalog.match_catalogs_radec`` helper it relies on. + +The estimator recovers ``m`` and ``c`` from five calibrated catalogues named +``1z2z`` (reference, no input shear), ``1p2z``/``1m2z`` (input shear +``g1 = +-|g|``) and ``1z2p``/``1z2m`` (``g2 = +-|g|``). The +g and -g sims are +matched to *each other* by RA/Dec -- same galaxies, opposite input shear -- and +the bias is the object-paired ("pool") average + + m = <(e_+ - e_-) / (2 |g|) - 1> , c = <(e_+ + e_-) / 2> , + +so the intrinsic shape cancels object-by-object in ``m``. + +We build synthetic catalogues in which the measured ellipticity is exactly +``e = (1 + m_true) g_in + c_true`` at shared positions, so the recovered m/c +must equal the injected values to machine precision -- an analytic check of +the estimator maths that needs no pipeline run. + +:Author: cdaley + +""" + +import numpy as np +import numpy.testing as npt +from astropy.io import fits + +from sp_validation.catalog import match_catalogs_radec +from sp_validation.image_sims import ImageSimMBias + +# Injected truth, shared across the synthetic-recovery test. +A = 0.02 # input shear amplitude |g| +M_TRUE = 0.05 # multiplicative bias (same for both components) +C1_TRUE = 0.001 # additive bias, component 1 +C2_TRUE = -0.002 # additive bias, component 2 +N_GAL = 2000 + + +def _write_cat(path, ra, dec, e1, e2, w): + """Write a minimal calibrated shape catalogue (RA, Dec, e1, e2, w_des).""" + cols = [ + fits.Column(name=name, array=arr, format="D") + for name, arr in ( + ("RA", ra), + ("Dec", dec), + ("e1", e1), + ("e2", e2), + ("w_des", w), + ) + ] + fits.HDUList([fits.PrimaryHDU(), fits.BinTableHDU.from_columns(cols)]).writeto( + path, overwrite=True + ) + + +def _make_grid(grids_dir, num): + """Create the five sheared/reference catalogues with a known m/c.""" + rng = np.random.default_rng(0) + ra = 30.0 + rng.uniform(0, 0.1, N_GAL) + dec = rng.uniform(0, 0.1, N_GAL) + w = np.ones(N_GAL) + zero = np.zeros(N_GAL) + + def e(g_in): + return (1 + M_TRUE) * g_in + zero + + sims = { + "1z2z": (C1_TRUE + zero, C2_TRUE + zero), + "1p2z": (e(+A) + C1_TRUE, C2_TRUE + zero), + "1m2z": (e(-A) + C1_TRUE, C2_TRUE + zero), + "1z2p": (C1_TRUE + zero, e(+A) + C2_TRUE), + "1z2m": (C1_TRUE + zero, e(-A) + C2_TRUE), + } + for name, (e1, e2) in sims.items(): + sim_dir = grids_dir / f"{name}_grid_{num}" + sim_dir.mkdir(parents=True, exist_ok=True) + _write_cat(sim_dir / "cat.fits", ra, dec, e1, e2, w) + + +def test_match_catalogs_radec_identity(): + """Identical positions match one-to-one; a shifted object drops out.""" + ra = np.array([30.0, 30.01, 30.02]) + dec = np.array([10.0, 10.01, 10.02]) + # Second catalogue = first, but the last object nudged well past threshold. + ra2, dec2 = ra.copy(), dec.copy() + ra2[2] += 1.0 + idx1, idx2 = match_catalogs_radec(ra, dec, ra2, dec2, thresh_deg=0.0002) + npt.assert_array_equal(idx2, [0, 1]) + npt.assert_array_equal(idx1, [0, 1]) + + +def test_mbias_recovers_injected_values(tmp_path): + """ImageSimMBias recovers the injected m/c to machine precision.""" + num = 7 + _make_grid(tmp_path, num) + config = { + "grids_dir": str(tmp_path), + "num": num, + "catalog_name": "cat.fits", + "shear_amplitude": A, + "match_radius_deg": 0.0002, + "w_cols": ["w_des"], + "n_bootstrap": 50, + "pair_match": True, + "bootstrap_seed": 42, + } + mb = ImageSimMBias(config) + mb.load_catalogs(verbose=False) + res = mb.run(verbose=False) + + npt.assert_allclose(res["m1"], M_TRUE, atol=1e-9) + npt.assert_allclose(res["m2"], M_TRUE, atol=1e-9) + npt.assert_allclose(res["c1"], C1_TRUE, atol=1e-9) + npt.assert_allclose(res["c2"], C2_TRUE, atol=1e-9) + # Bootstrap errors are non-negative and finite. + for key in ("m1_err", "m2_err", "c1_err", "c2_err"): + assert np.isfinite(res[key]) and res[key] >= 0 + # Self-describing results: the primary scheme is mirrored at the top level + # *and* lives under ``weights[scheme]``, and the two agree exactly. + assert list(res["weights"]) == ["w_des"] + for key in ("m1", "m1_err", "c1", "c1_err", "m2", "m2_err", "c2", "c2_err"): + assert res[key] == res["weights"]["w_des"][key] + + +def test_mbias_multiple_weight_schemes_share_draws(tmp_path): + """Multi-scheme runs key results per scheme; the primary mirrors the first. + + ``none`` (unit weights) and a real weight column are computed in one run. + With uniform per-object weights in the synthetic grid the two schemes give + the *same* m/c (the weighting is a no-op), and the shared bootstrap indices + make even the errors identical -- the property that lets a scheme + comparison be a clean weighting comparison. The first entry (``none``) is + the primary result surfaced at the top level. + """ + num = 8 + _make_grid(tmp_path, num) + config = { + "grids_dir": str(tmp_path), + "num": num, + "catalog_name": "cat.fits", + "shear_amplitude": A, + "match_radius_deg": 0.0002, + "w_cols": ["none", "w_des"], + "n_bootstrap": 50, + "pair_match": True, + "bootstrap_seed": 42, + } + mb = ImageSimMBias(config) + mb.load_catalogs(verbose=False) + res = mb.run(verbose=False) + + assert list(res["weights"]) == ["none", "w_des"] + # Primary (first) scheme mirrored at the top level. + for key in ("m1", "m1_err", "c1", "c1_err", "m2", "m2_err", "c2", "c2_err"): + assert res[key] == res["weights"]["none"][key] + # Uniform grid weights make the schemes agree bit-for-bit, errors included + # (shared bootstrap draws). + assert res["weights"]["none"] == res["weights"]["w_des"] + + +def test_mbias_deprecated_w_col_still_runs(tmp_path): + """A pre-``w_cols`` config with the scalar ``w_col`` still runs. + + The deprecated single-scheme key is honoured as ``[w_col]`` when ``w_cols`` + is absent, so a legacy run config keeps working and produces the same + single-scheme result as the ``w_cols=[w_col]`` spelling. + """ + num = 9 + _make_grid(tmp_path, num) + base = { + "grids_dir": str(tmp_path), + "num": num, + "catalog_name": "cat.fits", + "shear_amplitude": A, + "match_radius_deg": 0.0002, + "n_bootstrap": 50, + "pair_match": True, + "bootstrap_seed": 42, + } + res_dep = ImageSimMBias({**base, "w_col": "w_des"}) + res_dep.load_catalogs(verbose=False) + out_dep = res_dep.run(verbose=False) + + res_new = ImageSimMBias({**base, "w_cols": ["w_des"]}) + res_new.load_catalogs(verbose=False) + out_new = res_new.run(verbose=False) + + assert list(out_dep["weights"]) == ["w_des"] + assert out_dep["weights"] == out_new["weights"] + + +def test_mbias_pool_cancels_shape_noise(tmp_path): + """The paired estimator cancels intrinsic shape noise in m. + + With realistic per-galaxy intrinsic ellipticity (sigma_e ~ 0.3) shared + between the +g and -g sims plus small independent measurement noise, the + object-paired difference cancels the intrinsic shape, so sigma(m) is set by + the measurement noise (~1e-2), not the shape noise. An *unpaired* estimator + (differencing two independently-drawn means) would instead return + sigma(m) ~ sigma_e / (2 |g| sqrt(N)) -- an order of magnitude larger. We + assert the recovered error sits well below that shape-noise floor, which is + the property the pooling exists to deliver. + """ + num = 3 + rng = np.random.default_rng(1) + ra = 30.0 + rng.uniform(0, 0.1, N_GAL) + dec = rng.uniform(0, 0.1, N_GAL) + w = np.ones(N_GAL) + sigma_e, sigma_meas = 0.3, 0.01 + e1_int = rng.normal(0, sigma_e, N_GAL) # intrinsic shape, shared across sims + e2_int = rng.normal(0, sigma_e, N_GAL) + + def measured(g1_in, g2_in): + """Measured ellipticity = intrinsic + (1 + m) * input shear + noise.""" + e1 = e1_int + (1 + M_TRUE) * g1_in + rng.normal(0, sigma_meas, N_GAL) + e2 = e2_int + (1 + M_TRUE) * g2_in + rng.normal(0, sigma_meas, N_GAL) + return e1, e2 + + sims = { + "1z2z": measured(0, 0), + "1p2z": measured(+A, 0), + "1m2z": measured(-A, 0), + "1z2p": measured(0, +A), + "1z2m": measured(0, -A), + } + for name, (e1, e2) in sims.items(): + sim_dir = tmp_path / f"{name}_grid_{num}" + sim_dir.mkdir(parents=True, exist_ok=True) + _write_cat(sim_dir / "cat.fits", ra, dec, e1, e2, w) + + config = { + "grids_dir": str(tmp_path), + "num": num, + "catalog_name": "cat.fits", + "shear_amplitude": A, + "match_radius_deg": 0.0002, + "w_cols": ["w_des"], + "n_bootstrap": 200, + "pair_match": True, + "bootstrap_seed": 42, + } + mb = ImageSimMBias(config) + mb.load_catalogs(verbose=False) + res = mb.run(verbose=False) + + shape_noise_floor = sigma_e / (2 * A * np.sqrt(N_GAL)) # the unpaired error + for comp in (1, 2): + # m recovered within a few sigma of truth... + assert abs(res[f"m{comp}"] - M_TRUE) < 5 * res[f"m{comp}_err"] + # ...and its error is far below what an unpaired estimator would give. + assert res[f"m{comp}_err"] < 0.1 * shape_noise_floor diff --git a/src/sp_validation/tests/test_mask_overlay.py b/src/sp_validation/tests/test_mask_overlay.py new file mode 100644 index 00000000..26244053 --- /dev/null +++ b/src/sp_validation/tests/test_mask_overlay.py @@ -0,0 +1,85 @@ +"""The image-sim mask config is a declared overlay on the data mask config. + +The image-sim calibration does not keep an independent copy of the mask / +calibration config: it keeps the *data* config (``mask_v1.X.9.yaml``) as the one +home for the shared cuts, and declares the sim-specific delta in an overlay +(``mask_v1.X.9_im_sim.overlay.yaml``). ``im_compose_mask.py`` applies the +overlay to the base and must reproduce the committed runtime file +(``mask_v1.X.9_im_sim.yaml``) **byte-for-byte**. + +This guard locks that equality, so the two artefacts cannot drift: + +* if someone edits the runtime file without updating the overlay (or vice + versa), :func:`test_compose_reproduces_runtime_byte_identical` goes red; +* if the base config changes such that an overlay anchor no longer matches, + the compose fails loudly rather than emitting a wrong file -- + :func:`test_compose_fails_loud_on_stale_anchor` locks that fail-fast. + +The runtime file is a tracked input to ``im_init``; keeping it byte-stable is +what keeps the reproduction gate bit-exact, so this test's unit is bytes, not +parsed YAML. +""" + +import importlib.util +from pathlib import Path + +import pytest + + +def _repo_root() -> Path: + """Locate the repo root by walking up to the ``pyproject.toml`` marker.""" + for parent in Path(__file__).resolve().parents: + if (parent / "pyproject.toml").exists(): + return parent + raise RuntimeError("could not locate repo root (no pyproject.toml above test)") + + +_CALIB_DIR = _repo_root() / "config" / "calibration" +_BASE = _CALIB_DIR / "mask_v1.X.9.yaml" +_OVERLAY = _CALIB_DIR / "mask_v1.X.9_im_sim.overlay.yaml" +_RUNTIME = _CALIB_DIR / "mask_v1.X.9_im_sim.yaml" + + +def _compose_module(): + """Import ``workflow/scripts/im_compose_mask.py`` (lives outside the package).""" + path = _repo_root() / "workflow" / "scripts" / "im_compose_mask.py" + spec = importlib.util.spec_from_file_location("im_compose_mask", path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def test_compose_reproduces_runtime_byte_identical(): + """compose(base, overlay) == the committed runtime file, byte-for-byte.""" + import yaml + + compose = _compose_module().compose + overlay = yaml.safe_load(_OVERLAY.read_text()) + base_text = _BASE.read_text() + + composed = compose(base_text, overlay) + + assert composed == _RUNTIME.read_text(), ( + "compose(mask_v1.X.9.yaml, overlay) no longer reproduces " + "mask_v1.X.9_im_sim.yaml byte-for-byte -- the runtime file and its " + "declared overlay have drifted; reconcile one against the other." + ) + + +def test_compose_fails_loud_on_stale_anchor(): + """A base whose text no longer carries an overlay anchor aborts, not composes. + + This is the drift-proofing: the overlay anchors to verbatim base text, so if + the base config is edited such that an anchor vanishes, the compose must die + with a clear message rather than silently emit a file missing that delta. + """ + import yaml + + module = _compose_module() + overlay = yaml.safe_load(_OVERLAY.read_text()) + # Drop the IMAFLAGS_ISO cut from the base text so its overlay anchor no + # longer matches; compose must abort (SystemExit from die()). + mangled = _BASE.read_text().replace("IMAFLAGS_ISO", "SOMETHING_ELSE") + + with pytest.raises(SystemExit, match="out of sync with the base"): + module.compose(mangled, overlay) diff --git a/workflow/README.md b/workflow/README.md index 47015158..f20ffe8f 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -24,3 +24,36 @@ What does *not* belong here: Runs stay modular, not monolithic: a paper or run composes these rules with Snakemake's `module` directive under its own config and an output `prefix`, so each namespaces cleanly under `results//`. + +## Running on the cluster — the candide profile + +`profiles/candide/config.yaml` is the committed SLURM profile: it hands +Snakemake the candide executor, account, partition, node excludes, and per-job +floor, so scheduling is repo state rather than an operator's shell. Drive any +target with one command: + +```bash +snakemake --profile workflow/profiles/candide \ + -s workflow/image_sims/Snakefile \ + --configfile +``` + +For example, the image-sim m-bias chain end to end (`im_mbias` fans out one +SLURM job per branch × tile, MPI-free): + +```bash +snakemake --profile workflow/profiles/candide \ + -s workflow/image_sims/Snakefile \ + im_mbias --configfile my_run.yaml +``` + +Give the target *before* `--configfile`: `--configfile` takes one-or-more +paths, so a target after it is read as a config file ("No such file: +im_mbias"). Always dry-run first with `-n`. + +The profile carries only cluster policy — no container settings (the image-sims +rules own their `apptainer exec` call) and no `OMP_NUM_THREADS` (pinned to 1 at +that same `apptainer exec` line, since the slurm executor's `--export=ALL` +propagates the driver's env, not a profile flag). Per-rule `mem_mb` / `runtime` +stay on the rules. Off-cluster, drop `--profile` and add `-j N`. See the +profile's own comments for the full rationale. diff --git a/workflow/Snakefile b/workflow/Snakefile index 6f59827b..7d91db81 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -47,3 +47,8 @@ include: "rules/glass_mock.smk" # guarded so paper configs without it (e.g. bmodes) don't trip on the lookups. if "cosmo_val" in config: include: "rules/cosmo_val.smk" + +# Image-simulation m/c-bias chain (image_sims.smk). Active only when the config +# carries an `image_sims` block; standalone runs use workflow/image_sims/Snakefile. +if "image_sims" in config: + include: "rules/image_sims.smk" diff --git a/workflow/image_sims/Snakefile b/workflow/image_sims/Snakefile new file mode 100644 index 00000000..d7a89065 --- /dev/null +++ b/workflow/image_sims/Snakefile @@ -0,0 +1,46 @@ +"""Standalone entry point for the image-simulation m-bias workflow. + +Run the sp_validation-side chain (merge -> extract -> calibrate -> m-bias), +optionally including the ShapePipe pipeline stage, without pulling in the +cosmology-validation config the top-level ``workflow/Snakefile`` requires. + +Layer a run config over the operational defaults -- the workflow config.yaml +carries operational defaults but *no* science keys, so it is incomplete on its +own (by design); the run config supplies the science knobs. The one drive +command on candide, with the committed SLURM profile owning all scheduling: + + snakemake --profile workflow/profiles/candide \\ + -s workflow/image_sims/Snakefile \\ + im_mbias --configfile my_run.yaml + +``configfile: "workflow/image_sims/config.yaml"`` below loads the operational +defaults automatically, so only ``my_run.yaml`` (the science knobs, and any +operational override that run wants) is passed on the command line; Snakemake +deep-merges the two. The profile supplies the executor, account, partition, +node excludes and job floor -- no ``-j`` needed (the slurm executor sets the +job cap). Off-cluster, drop ``--profile`` and add ``-j N`` to run locally. + +The target (``im_mbias``) is given *before* ``--configfile``: Snakemake's +``--configfile`` takes one-or-more paths, so a target placed after it is +swallowed as a config path ("No such file: im_mbias"). Put targets ahead of +``--configfile`` (or make ``--configfile`` the last flag on the line). Always +dry-run first with ``-n``. + +The same rules are also available inside the main workflow: they are included +there under ``if "image_sims" in config``. +""" + +configfile: "workflow/image_sims/config.yaml" + + +# The image-sims rules own their container invocation explicitly, so no +# top-level container is needed here. +container: None + + +include: "../rules/image_sims.smk" + + +rule all: + input: + f"{GRIDS_BASE}/results/m_bias_results.yaml", diff --git a/workflow/image_sims/config.yaml b/workflow/image_sims/config.yaml new file mode 100644 index 00000000..8f719c68 --- /dev/null +++ b/workflow/image_sims/config.yaml @@ -0,0 +1,89 @@ +# Image-simulation m-bias workflow configuration. +# +# Two kinds of keys live under `image_sims:`, and the split is the point: +# +# * OPERATIONAL keys default here (active lines below) and *nowhere else* -- +# the .smk reads them bare, so this file is their single home. Override in +# a run config only when a run genuinely differs from the shared setup. +# +# * SCIENCE keys have NO default -- not here, not in code. They fix the +# estimator's scientific behaviour and must be stated per run, so they +# appear below only as commented template lines. Supply them in a run +# config layered on top: +# +# snakemake -s workflow/image_sims/Snakefile \ +# --configfile workflow/image_sims/config.yaml \ +# --configfile my_run.yaml \ +# -j 4 im_mbias +# +# A run config that omits a science key fails at DAG parse, naming the key; an +# unknown key under `image_sims:` fails as a typo. The structural keys below +# (sif, repos, data roots, num, tile_ids) also have no default and must be set. + +image_sims: + + # --- containers ------------------------------------------------------- + # Two images, one per half of the chain (the split gate766 ran). One image + # is the eventual target -- the sp_validation image is FROM the ShapePipe + # image -- but until sp_validation is uv-locked with cosmo_numba declared, + # its published image can drift NumPy past numba's window (seen 2026-07-11: + # "Numba needs NumPy 2.4 or less. Got NumPy 2.5" at ngmix). PYTHONPATH + # shadows pure-Python code only, never binary deps. + sif: /n17data/cdaley/containers/sp_validation_im_sims.sif # extract/calibrate/m-bias + sif_pipeline: /n17data/cdaley/containers/shapepipe_im_sims-runtime.sif # pipeline/merge + # Apptainer bind mounts. /automnt is required when repos/data are + # automounted (candide gotcha); harmless otherwise. [operational] + binds: /n17data,/n09data,/home,/automnt + + # --- repositories ----------------------------------------------------- + # Bound into the image; both repos' src go on PYTHONPATH so this branch's + # code wins over the baked copies: ShapePipe's #766 build, and sp_validation's + # image_sims.py / catalog.match_catalogs_radec. + shapepipe_repo: /n17data/cdaley/unions/code/shapepipe + sp_validation_repo: /n17data/cdaley/unions/code/sp_validation + + # --- data and run directories ---------------------------------------- + # grids_base is the run/output root: one sub-directory per simulation. + grids_base: /n17data/cdaley/unions/scratch-wf/imsims-run/grids + input_sims_base: /n09data/hervas/skills_out + psf_dict: /home/hervas/fhervas/workdir_skills/input/psf_files/Full_psf_dict.pickle + + # --- simulation grid -------------------------------------------------- + sims_type: grid # 'grid' -> *_grid_{num}; anything else -> *_{num} [operational] + num: 1 + # Branches this run requests: the unsheared reference plus the four +/- + # sheared branches. Injected shear is NOT set here -- it is parsed from each + # branch's basic_info.txt by im_manifest. [operational] + branches: ["1z2z", "1p2z", "1m2z", "1z2p", "1z2m"] + # Tiles to process: an explicit list, the one tile-input mechanism. + tile_ids: ["233.293", "237.292", "238.292"] + + # --- calibration ------------------------------------------------------ + shape: ngmix # [operational] + # ShapePipe cfis configs (final_cat.param etc.); default is + # {shapepipe_repo}/example/cfis_image_sims once #766 lands. [operational] + config_dir: /n17data/cdaley/unions/scratch-wf/imsims-run/grids/_cfis_image_sims + psf_model: psfex # [operational] + n_smp: -1 # [operational] + # Extract/calibrate scripts run from the sp_validation repo checkout (branch + # code, not the baked copies). Point elsewhere for a different checkout. + # [operational] + extract_script: /n17data/cdaley/unions/code/sp_validation/scripts/calibration/extract_info.py + calibrate_script: /n17data/cdaley/unions/code/sp_validation/scripts/calibration/calibrate_comprehensive_cat.py + + # --- science knobs (REQUIRED in the run config; no default) ----------- + # Copy these into your run config and set them. There is deliberately no + # default: each fixes the estimator's scientific behaviour, so a run must + # state it. The injected |g| is separate -- parsed from basic_info.txt by + # im_manifest, never set here. + # + # mask_config: config/calibration/mask_v1.X.9_im_sim.yaml # relative to sp_validation_repo + # match_radius_deg: 0.0002 # RA/Dec pair-match radius, degrees + # w_cols: [none, w_iv] # weight schemes to compute in one run; "none" + # # (or null) = unit weights (#227); first entry is + # # the primary/headline result -- lead with "none" + # # for the fiducial unweighted estimator. The + # # deprecated scalar `w_col` is still accepted. + # pair_match: true # match +g/-g object-by-object (per-object cancellation) + # n_bootstrap: 500 # bootstrap resamples for the errors + # bootstrap_seed: 42 # seed for the bootstrap RNG (makes errors reproducible) diff --git a/workflow/image_sims/params_im_sim.py b/workflow/image_sims/params_im_sim.py new file mode 100644 index 00000000..e87124d0 --- /dev/null +++ b/workflow/image_sims/params_im_sim.py @@ -0,0 +1,228 @@ +""" + +:Name: params.py + +:Description: This script contains parameters to run the validation notebook. + +:Author: Martin Kilbinger + +:Date: 2021 + +:Package: sp_validation + +""" + +import os + +import numpy as np + +# Control + +## Verbose output +verbose = True + +## Math output +np.set_printoptions(precision=3, formatter={"float": "{: .3g}".format}) + + +# Survey parameters + +## Field or patch name -- derived from the run directory, which is named +## after the simulation (e.g. '1z2z_grid_1'), so one shared params file +## serves every sim. +name = os.path.basename(os.getcwd()) +print("Field name = {}".format(name)) + +## Area of a tile in deg^2 +area_tile = 0.25 + +## Pixel size in arcsec +pixel_size = 0.187 + +## Shape measurement method, implemented is +## 'ngix': multi-epoch model fitting +## 'galsim': stacked-image moments (experimental) +shape = "ngmix" + +# Paths + +## Input paths + +### Input data directory +data_dir = "." + +### Tile IDs +path_tile_ID = f"{data_dir}/tiles_{name}.txt" + +### Weak-lensing galaxy catalog name +galaxy_cat_path = f"{data_dir}/final_cat_{name}.hdf5" +print(f"Galaxy catalogue = {galaxy_cat_path}") + +## Parameter list; optional, set to `None` if not required +param_list_path = f"{data_dir}/cfis/final_cat.param" + +### Star and PSF catalog name; optional, set to `None` if not required +star_cat_path = None + +# HDU number of star and PSF catalogue +hdu_star_cat = 1 + +### External mask; optional, set to `None` if not required +mask_external_path = None + +## Output paths + +### Output base directory +output_dir = f"{data_dir}" + +### Galaxy shape catalogue base name +output_shape_cat_base = f"{output_dir}/shape_catalog" + +### PSF output catalogue base name. +output_PSF_cat_base = f"{output_dir}/psf_catalog" + +### File for found tile IDs +path_found_ID = f"{output_dir}/found_ID.txt" + +### File for missing tile IDs +path_missing_ID = f"{output_dir}/missing_ID.txt" + +### Plot directory and subdirs +plot_dir = f"{output_dir}/plots/" + +### Statistics text file +stats_file_name = "stats_file.txt" + +# Other IO options + +## Input + +### Coordinate column names +col_name_ra = "XWIN_WORLD" +col_name_dec = "YWIN_WORLD" + +### Memory mode, set to None unless very large file +mmap_mode = None + +## Output + +### Output file format extension: '.fits' or '.hdf5' +output_format = ".fits" + +### Additional output columns +add_cols = [ + "FLUX_RADIUS", + "FWHM_IMAGE", + "FWHM_WORLD", + "MAGERR_AUTO", + "MAG_WIN", + "MAGERR_WIN", + "FLUX_AUTO", + "FLUXERR_AUTO", + "FLUX_APER", + "FLUXERR_APER", + "NGMIX_T_NOSHEAR", + "NGMIX_T_PSF_RECONV_NOSHEAR", +] + +## Pre-calibration catalogue, including masked objects and mask flags. +## ShapePipe-v2 (post-#761) ngmix grammar: ellipticity in named scalar +## components NGMIX_G{1,2}_*, PSF size split into NGMIX_T_PSF_ORIG/RECONV. +## IMAFLAGS_ISO (present in the data-path params) is omitted: the simulation +## pipeline runs no imaging-flag masking stage, so the column does not exist. +## NGMIX_MCAL_TYPES_FAIL is kept -- it is the metacal moments-failure flag the +## calibration mask cuts on, identically to the data path. +add_cols_pre_cal = [ + "TILE_ID", + "NUMBER", + "FLAGS", + "NGMIX_MCAL_FLAGS", + "NGMIX_MCAL_TYPES_FAIL", + "N_EPOCH", + "NGMIX_N_EPOCH", + "NGMIX_G1_PSF_ORIG_NOSHEAR", + "NGMIX_G2_PSF_ORIG_NOSHEAR", + "NGMIX_G1_ERR_NOSHEAR", + "NGMIX_G2_ERR_NOSHEAR", +] + +### Set flag columns as integer format +add_cols_pre_cal_format = {} +for key in ( + "NUMBER", + "FLAGS", + "NGMIX_MCAL_FLAGS", + "NGMIX_MCAL_TYPES_FAIL", + "N_EPOCH", + "NGMIX_N_EPOCH", +): + add_cols_pre_cal_format[key] = "I" + +add_cols_pre_cal_format["TILE_ID"] = "A7" +add_cols_pre_cal_format["NUMBER"] = "J" + +# Create key names for metacal information +prefix = "NGMIX" +suffixes = ["1M", "1P", "2M", "2P", "NOSHEAR"] +centers = ["FLAGS", "G1", "G2", "FLUX", "FLUX_ERR", "T", "T_ERR", "T_PSF_RECONV"] +for center in centers: + for suffix in suffixes: + add_cols_pre_cal.append(f"{prefix}_{center}_{suffix}") + +for suffix in suffixes: + add_cols_pre_cal_format[f"FLAGS_{suffix}"] = "I" + + +# Catalog parameters + +## Star matching threshold [deg] +thresh = 0.0002 + +## Number of jackknife resamples for additive bias +## (0: no jackknife computation). +## If < 2000 the jackknife mean fluctuates a lot. +n_jack = 0 + + +## Galaxy selection + +# Flag to output selected and calibrated galaxy catalogue (<= SP v1.4.1). +# If False, only output comprehensive catalogue. +do_selection_calibration = False + +## Magnitude limits +gal_mag_bright = 15 +gal_mag_faint = 30 + +### Spread-model +do_spread_model = False + +### SExtractor flags to keep in addition to FLAGS=0 +### (bit-coded; list of powers of 2); +### Empty list if no flags +flags_keep = [] + +## Minimum number of epochs +n_epoch_min = 2 + +### Signal-to-noise (selection within metacal) +#### minimum to cut noisy objects +gal_snr_min = 10 +#### maximum to cut too bright objects, potentially too large for the postage stamp +gal_snr_max = 500 + +### Relative size, T_gal / T_psf (selection within metacal) +### to select objects that are not too small compared to the PSF, thus not likely to be point-like, +### or to big as they seem to bias the correlation functions +gal_rel_size_min = 0.5 +gal_rel_size_max = 3.0 + +### Correct galaxy size for ellipticity +gal_size_corr_ell = False + +### prior ellipticity dispersion (one component), *only* used for galaxy weight +sigma_eps_prior = 0.34 + + +## Wrap coordinates around this value [deg], set to != 0 if ra=0 is within coordinate range +wrap_ra = 0 diff --git a/workflow/profiles/candide/config.yaml b/workflow/profiles/candide/config.yaml new file mode 100644 index 00000000..0316b909 --- /dev/null +++ b/workflow/profiles/candide/config.yaml @@ -0,0 +1,85 @@ +# Committed SLURM profile for the candide cluster (IAP). +# +# This is the "one run command" half of the workflow: drive any target with +# +# snakemake --profile workflow/profiles/candide \ +# -s workflow/image_sims/Snakefile \ +# --configfile +# +# and Snakemake owns all scheduling -- it fans out one SLURM job per branch x +# tile and drives them against the cluster, MPI-free. Everything here is +# cluster policy (executor, account, partition, node excludes, per-job +# defaults); it carries no science and no workflow logic. +# +# What is deliberately NOT here: +# +# * Container / apptainer settings. The image-sims rules set +# ``container: None`` and own their ``apptainer exec`` call through the +# shared ``EXEC`` prefix (one image for every stage, with PYTHONPATH / +# PSF_DICT / OMP_NUM_THREADS injected there). So no +# ``software-deployment-method: apptainer`` / ``apptainer-args`` -- those +# would wrap a *second*, redundant container around jobs that already run +# inside one. +# +# * OMP_NUM_THREADS. It is pinned to 1 on the ``apptainer exec`` line in +# workflow/rules/image_sims.smk, not here. The slurm executor submits with +# ``--export=ALL``, which propagates the *driver's* ambient environment; a +# profile only sets CLI flags, never the driver's own env, so an +# ``OMP_NUM_THREADS`` set here would silently depend on the operator having +# exported it by hand. Injecting it at the container boundary puts it where +# the compute runs, committed and independent of the launching shell. +# +# * Per-rule resources (mem_mb, runtime). Those live on each rule in the +# .smk; the ``default-resources`` below are only the floor for rules that +# set none. + +executor: slurm + +# Cluster policy applied to every job unless a rule overrides it. The excludes +# are the flaky/no-internet candide nodes (n17 mount issues, n09 no internet, +# n36); ``slurm_extra`` is passed verbatim onto the sbatch line by the executor +# plugin, so the quoting is what sbatch must see. +# +# Three candide-specific SLURM lessons are baked into the values below (learned +# the hard way on the earlier hand-driven im-sims runs; see shapepipe's retired +# image_sims_pipeline/Snakefile docstring): +# +# * ``runtime`` MUST carry a unit (``60m``, ``6h``, ``2d``). Snakemake's +# resource parser reads a *bare* number as SECONDS, so ``runtime: 60`` would +# silently give every job a 60-second wall clock and kill it on start. The +# quoted-with-unit form here is deliberate; keep it that way, and prefer the +# same in any ``--default-resources`` passed on the command line. (A bare +# integer in a *rule's* ``resources: runtime=720`` is fine -- snakemake +# reads rule-level numeric runtime as minutes -- the seconds trap is only +# the CLI/default-resources parser.) +# +# * ``cpus_per_task`` is pinned to 12 to CAP JOBS PER NODE, not because a job +# needs 12 cores (the chain is MPI-free and pins ``OMP_NUM_THREADS=1`` at +# the container). candide's per-user process limit is ``ulimit -u 1200`` +# *per node*, and apptainer crashes ("can't start new thread") beyond ~4 +# concurrent jobs on a 48-core node. Requesting 12 CPUs/job holds SLURM to +# ~4 jobs per 48-core node, under the ceiling. Dropping this to 1 would let +# SLURM pack ~48 jobs onto a node and crash the compute-heavy im_pipeline +# stage (which inherits this default -- it sets mem/runtime but not cpus). +# +# * After launching a real fan-out, VERIFY the request actually landed: +# ``squeue -u $USER -o "%C %l"`` must show 12 (CPUs) and the wall clock you +# intended (e.g. 12:00:00 for im_pipeline). A silently-misparsed runtime or +# cpus shows up here before it wastes a queue slot. +default-resources: + slurm_account: "cusers" + slurm_partition: "comp,pscomp" + runtime: "60m" + cpus_per_task: 12 + slurm_extra: "'--exclude=n17,n09,n36'" + +# Give an appearing output file a moment on candide's automounted filesystems +# before Snakemake calls a job failed for a missing output, and retry a job +# once on transient node failure. +latency-wait: 5 +retries: 1 + +# Keep the SLURM logs of successful jobs (candide debugging), and rerun a job +# when its code / params / inputs change, not only on mtime. +slurm-keep-successful-logs: true +rerun-triggers: ["mtime", "params", "input", "code"] diff --git a/workflow/rules/image_sims.smk b/workflow/rules/image_sims.smk new file mode 100644 index 00000000..a4eb92d0 --- /dev/null +++ b/workflow/rules/image_sims.smk @@ -0,0 +1,568 @@ +"""Image-simulation orchestration: raw SKiLLS sim images -> shear m/c bias. + +This rule set drives the image-simulation validation chain end to end and is +the sp_validation-side half of the split described in +``UNIONS-WL/MultiBand_ImSim#1``: ShapePipe turns the simulated tiles into +per-tile shape catalogues, then sp_validation merges, extracts, calibrates and +finally measures the multiplicative/additive shear bias. + +Two images, one prefix shape. Architecturally one image could run every +stage -- the sp_validation image is built ``FROM`` the ShapePipe image, so it +carries both stacks -- but the *published* sp_validation image's environment +is not yet trustworthy for the ShapePipe half: sp_validation has no lockfile +and does not declare its numba-bearing dependency (``cosmo_numba``), so +unpinned install layers can drift NumPy past numba's window (a 2026-07-11 +gate run hit exactly this: ``Numba needs NumPy 2.4 or less. Got NumPy 2.5`` +at the ngmix stage). PYTHONPATH shadowing covers pure-Python *code*, never +binary deps, so until sp_validation is uv-locked with its deps declared +(spun off as its own task), each half runs in its own repo's image -- the +same split the gate766 baseline ran: + +* ShapePipe stages -> ``pipeline`` (raw images -> per-tile cats) and ``merge`` + (``create_final_cat`` -> ``final_cat_{sim}.hdf5``) run in ``sif_pipeline`` + (the ShapePipe image). +* sp_validation stages -> ``manifest``, ``extract`` (-> comprehensive cat), + ``calibrate`` (-> cut cat) and ``m_bias`` (-> ``m_bias_results.yaml``) run + in ``sif`` (the sp_validation image). + +Every rule sets ``container: None`` and calls ``apptainer exec`` explicitly +through a shared prefix template (``EXEC_PIPELINE`` / ``EXEC`` -- identical +env injections, different image), because the images are not the workflow's +top-level container. Everything is parameterised under +``config["image_sims"]`` -- the two ``sif`` keys, repository roots, data roots, +the PSF dictionary, the explicit ``tile_ids`` list and the sim/calibration +knobs -- so a fresh user drives it from config alone, with no hard-coded clone +layout. Configuration is fail-fast: a schema check at load rejects an unknown +key (typo) and a missing science key (see ``workflow/image_sims/config.yaml`` +for the operational/science split). The ``PYTHONPATH`` override injects both +repos' ``src`` so the *branch* source (ShapePipe's ``#766`` build; +sp_validation's ``image_sims.py``, ``catalog.match_catalogs_radec``) wins over +whatever is baked into the image. + +The five simulations per grid are the reference ``1z2z`` (no input shear) plus +the ``+/-`` shear pairs ``1p2z``/``1m2z`` (g1) and ``1z2p``/``1z2m`` (g2); the +m-bias estimator matches each to the reference by RA/Dec. +""" + +import os +from pathlib import Path + +IMSIM = config["image_sims"] + +# --- fail-fast schema check ---------------------------------------------- +# One home for every fact: the run config carries the science knobs, the +# workflow config.yaml carries the operational defaults, and *this* block is +# where a typo or a missing knob dies -- at DAG parse, before any compute. +# +# Every key must be declared below. An unknown key under ``image_sims:`` is a +# hard error (typo protection); a missing *science* key is a hard error naming +# the key (no silent code default anywhere). Operational keys default in the +# workflow config.yaml and nowhere else: the .smk reads them as bare +# ``IMSIM[key]`` (never ``.get`` with a second literal), so their value comes +# from config.yaml alone -- the single home for an operational default. +# +# Science keys: required from the *run* config; no default in config.yaml (only +# a commented template line) and no default in code. These fix the estimator's +# scientific behaviour, so they must be stated per run, never inherited. +_SCIENCE_KEYS = { + "w_cols", + "pair_match", + "match_radius_deg", + "n_bootstrap", + "bootstrap_seed", + "mask_config", +} +# Deprecated science keys: accepted (so a pre-``w_cols`` run config still parses +# and the estimator's back-compat path runs) but not *required* -- our configs +# state ``w_cols``. Listed here only to keep them out of the unknown-key error. +_DEPRECATED_KEYS = { + "w_col", +} +# Operational keys: default (visibly) in the workflow config.yaml; the .smk +# reads them bare, so config.yaml is their one home. +_OPERATIONAL_KEYS = { + "binds", + "sims_type", + "branches", + "shape", + "config_dir", + "psf_model", + "n_smp", + "extract_script", + "calibrate_script", +} +# Structural keys: paths/identifiers the run must supply (no sensible default). +_STRUCTURAL_KEYS = { + "sif", + "sif_pipeline", + "shapepipe_repo", + "sp_validation_repo", + "grids_base", + "input_sims_base", + "psf_dict", + "num", + "tile_ids", +} +_ALLOWED_KEYS = ( + _SCIENCE_KEYS | _DEPRECATED_KEYS | _OPERATIONAL_KEYS | _STRUCTURAL_KEYS +) + +_unknown = set(IMSIM) - _ALLOWED_KEYS +if _unknown: + raise ValueError( + "image_sims: unknown config key(s) " + f"{sorted(_unknown)} -- check for a typo (allowed keys: " + f"{sorted(_ALLOWED_KEYS)})" + ) +_missing_science = sorted(_SCIENCE_KEYS - set(IMSIM)) +if _missing_science: + raise ValueError( + "image_sims: missing required science key(s) " + f"{_missing_science} -- these have no default and must be set in the " + "run config (see the commented template in workflow/image_sims/config.yaml)" + ) +_missing_structural = sorted(_STRUCTURAL_KEYS - set(IMSIM)) +if _missing_structural: + raise ValueError( + "image_sims: missing required key(s) " + f"{_missing_structural} -- set them in the run config" + ) + +# --- containers ----------------------------------------------------------- +# Two images (see module docstring): the ShapePipe image for the pipeline and +# merge stages, the sp_validation image for everything downstream. Collapse +# back to one image once sp_validation's env is lock-managed. +SIF = IMSIM["sif"] # sp_validation stages +SIF_PIPELINE = IMSIM["sif_pipeline"] # ShapePipe stages +BINDS = IMSIM["binds"] + +# --- repositories (bound into the image; branch code overrides) ----------- +SHAPEPIPE_REPO = IMSIM["shapepipe_repo"] +SPV_REPO = IMSIM["sp_validation_repo"] + +# --- data and run directories -------------------------------------------- +GRIDS_BASE = IMSIM["grids_base"] # run/output root; one sub-dir per sim +INPUT_SIMS_BASE = IMSIM["input_sims_base"] # SKiLLS sim images +PSF_DICT = IMSIM["psf_dict"] # Herve's Full_psf_dict.pickle + +# --- simulation grid ------------------------------------------------------ +# SIM_BASES is the set of branches *this run requests* -- the reference plus the +# four +/- sheared branches. Their injected shear (amplitude, per-branch +# (g1,g2), pairing) is NOT a literal here: it lives only in manifest.yaml, built +# by im_manifest from each branch's basic_info.txt and read back by im_mbias. +NUM = IMSIM["num"] +SIMS_TYPE = IMSIM["sims_type"] +_SUFFIX = f"_{SIMS_TYPE}_{NUM}" if SIMS_TYPE == "grid" else f"_{NUM}" +SIM_BASES = list(IMSIM["branches"]) +SIMS = [f"{base}{_SUFFIX}" for base in SIM_BASES] +MANIFEST = f"{GRIDS_BASE}/manifest.yaml" +BUILD_MANIFEST = f"{SPV_REPO}/workflow/scripts/im_build_manifest.py" + +# --- tiles ---------------------------------------------------------------- +# tile_ids is the one tile-input mechanism: an explicit list in the run config. +TILE_IDS = list(IMSIM["tile_ids"]) + +# --- calibration / m-bias knobs ------------------------------------------ +SHAPE = IMSIM["shape"] +MASK_CONFIG = IMSIM["mask_config"] # e.g. config/calibration/mask_v1.X.9_im_sim.yaml +PARAMS_TEMPLATE = f"{SPV_REPO}/workflow/image_sims/params_im_sim.py" +# ShapePipe cfis_image_sims config dir (per-tile/exposure configs + final_cat.param). +CONFIG_DIR = IMSIM["config_dir"] + +# ShapePipe scripts live in the ShapePipe repo (also baked into its image). +CREATE_FINAL_CAT = f"{SHAPEPIPE_REPO}/scripts/python/create_final_cat.py" +RUN_JOB = f"{SHAPEPIPE_REPO}/scripts/sh/run_job_sp_canfar_v2.0.bash" +# Extract/calibrate run from the sp_validation *repo* checkout (bind-mounted), +# not the baked copies: the container tracks the branch but lags it, and the +# image-sims path needs branch-only fixes (star-catalogue-optional extract, +# FITS-aware CalibrateCat.read_cat). Overridable for a different checkout. +EXTRACT_INFO = IMSIM["extract_script"] +CALIBRATE = IMSIM["calibrate_script"] +# m-bias is *this branch's* extracted core, injected on PYTHONPATH. +COMPUTE_M_BIAS = f"{SPV_REPO}/scripts/compute_m_bias_image_sims.py" + +# --- container exec prefixes ---------------------------------------------- +# One prefix *shape* for every stage -- two instances, one per image. Three +# env injections make the on-disk branch +# code and the sim PSF win over the image's baked copies: +# +# * PYTHONPATH prepends BOTH repos' ``src`` (ShapePipe first, then +# sp_validation), so Python resolves the worktree build before +# ``/app``/``/sp_validation`` -- the local-testing counterpart of the +# git-ref deps, letting the branch code run without an image rebuild. This +# covers the Python *packages* only: the bash entry points (run_job) and +# the ShapePipe/sp_validation *scripts* are still invoked at the repo paths +# resolved from config (RUN_JOB, CREATE_FINAL_CAT, EXTRACT_INFO, ...), not +# shadowed by PYTHONPATH. +# * PSF_DICT points the fake_psf module (PSF_DICT_PATH = $PSF_DICT, expanded +# via getexpanded) at this run's PSF dictionary. +# +# The SLURM env vars are stripped (``env -u ...``) so that when the ShapePipe +# pipeline stage's OpenMPI initialises inside the image it does not try to +# attach to the host SLURM launcher (cf. apptainer_noslurm.sh). The strip is +# harmless for the pure-Python sp_validation stages, so one prefix serves all. +# +# ``OMP_NUM_THREADS=1`` is injected here, at the ``apptainer exec`` call, and +# not left to the SLURM profile. The chain is MPI-free: Snakemake fans out one +# job per branch x tile and each job's parallelism is ShapePipe's own internal +# multiprocessing (``-N n_smp``), so the OpenMP/BLAS thread pool inside the +# container must be pinned to 1 to avoid oversubscription. The SLURM profile +# cannot pin it reliably: the slurm executor submits with ``--export=ALL``, +# which propagates the *driver's* ambient environment -- but a Snakemake +# profile only sets CLI flags, never the driver's own env, so an +# ``OMP_NUM_THREADS`` there would depend on the operator having exported it by +# hand (the implicit, uncommitted state the "one run command" is meant to +# retire). Injecting it on the ``apptainer exec`` line puts it where the +# compute actually runs -- inside the container, independent of the driver's +# env -- the same lever this prefix already uses for PYTHONPATH/PSF_DICT. +_EXEC_PREFIX = ( + "env -u SLURM_JOBID -u SLURM_JOB_ID -u SLURM_PROCID " + f"apptainer exec --bind {BINDS} " + f"--env PYTHONPATH={SHAPEPIPE_REPO}/src:{SPV_REPO}/src " + f"--env PSF_DICT={PSF_DICT} --env OMP_NUM_THREADS=1 " +) +EXEC = _EXEC_PREFIX + SIF # sp_validation stages +EXEC_PIPELINE = _EXEC_PREFIX + SIF_PIPELINE # ShapePipe stages + +JOB_MASK = sum([1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048]) + + +wildcard_constraints: + sim="|".join(SIMS), + tile="|".join(t.replace(".", r"\.") for t in TILE_IDS), + + +# ========================================================================== +# Convenience targets (run in order) +# ========================================================================== +rule im_manifest_only: + input: + MANIFEST, + + +rule im_init_all: + input: + expand(f"{GRIDS_BASE}/{{sim}}/params.py", sim=SIMS), + + +rule im_pipeline_all: + input: + expand( + f"{GRIDS_BASE}/{{sim}}/logs/pipeline_{{tile}}.done", + sim=SIMS, + tile=TILE_IDS, + ), + + +rule im_merge_all: + input: + expand(f"{GRIDS_BASE}/{{sim}}/final_cat_{{sim}}.hdf5", sim=SIMS), + + +rule im_extract_all: + input: + expand( + f"{GRIDS_BASE}/{{sim}}/shape_catalog_comprehensive_{SHAPE}.fits", + sim=SIMS, + ), + + +rule im_calibrate_all: + input: + expand( + f"{GRIDS_BASE}/{{sim}}/shape_catalog_cut_{SHAPE}.fits", sim=SIMS + ), + + +# ========================================================================== +# Rules +# ========================================================================== +rule im_manifest: + """Build the campaign manifest at the head of the DAG. + + Parses ``g_cosmic`` from every requested branch's ``basic_info.txt``, + cross-checks each against its ``1{X}2{Y}`` name and the (0,0) reference, + derives the single injected amplitude, and writes ``manifest.yaml`` into the + run root. This is the one home for the injected-shear facts; im_mbias reads + the amplitude and branch map from here, nowhere else. Pure sp_validation + stage (stdlib parse of basic_info; PyYAML to write). + """ + input: + # basic_info.txt for each requested branch, so editing a sim's record + # rebuilds the manifest (and re-validates) rather than reusing a stale one. + basic_info=expand( + f"{INPUT_SIMS_BASE}/{{sim}}/basic_info.txt", sim=SIMS + ), + output: + manifest=MANIFEST, + params: + branch_args=lambda wc: " ".join(f"--branch {b}" for b in SIM_BASES), + input_sims_base=INPUT_SIMS_BASE, + sims_type=SIMS_TYPE, + num=NUM, + shell: + "{EXEC} python {BUILD_MANIFEST} " + "--input-sims-base {params.input_sims_base} " + "--sims-type {params.sims_type} --num {params.num} " + "{params.branch_args} -o {output.manifest}" + + +rule im_init: + """Stage per-sim run directory: params.py, mask config, ShapePipe configs, + and the raw SKiLLS image inputs. + + ``params_im_sim.py`` derives the field name from the directory basename, so + the same template serves every sim; ``config_mask.yaml`` and ``cfis`` are + symlinks the downstream calibration and merge steps read from cwd. + + ``input_tiles``/``input_exp`` are top-level symlinks to the raw SKiLLS tile + and exposure images; ShapePipe's ``get_images_runner`` resolves them via + ``$SP_DIR/input_{tiles,exp}`` (``$SP_DIR`` is the run dir). ``run_job`` does + not stage these, so ``im_init`` must -- this is what makes ``im_pipeline`` + runnable from raw images, not just from pre-staged intermediates. + """ + input: + # Tracked so that editing the params template or mask config re-stages + # them into every run dir (a plain params: value would not retrigger, + # silently leaving stale params.py behind after a grammar change). + template=PARAMS_TEMPLATE, + mask_src=os.path.join(SPV_REPO, MASK_CONFIG), + output: + params=f"{GRIDS_BASE}/{{sim}}/params.py", + mask=f"{GRIDS_BASE}/{{sim}}/config_mask.yaml", + params: + config_dir=CONFIG_DIR, + run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", + cfis=lambda wc: f"{GRIDS_BASE}/{wc.sim}/cfis", + sim_tiles=lambda wc: f"{INPUT_SIMS_BASE}/{wc.sim}/images/SP_tiles", + sim_exp=lambda wc: f"{INPUT_SIMS_BASE}/{wc.sim}/images/SP_exp", + shell: + # cfis / input_tiles / input_exp are stable read-only symlinks (used by + # get_images, merge, extract); created here but not tracked as outputs, + # which snakemake will not accept for a symlink/directory. + "mkdir -p $(dirname {output.params}) && " + "cp {input.template} {output.params} && " + "ln -sf {input.mask_src} {output.mask} && " + "ln -sfT {params.config_dir} {params.cfis} && " + "ln -sfT {params.sim_tiles} {params.run_dir}/input_tiles && " + "ln -sfT {params.sim_exp} {params.run_dir}/input_exp" + + +rule im_pipeline: + """Run ShapePipe on one simulated tile (ShapePipe stage). + + Delegates the module DAG to ShapePipe's own job runner; the sentinel log + marks tile completion for the merge step. This is the compute-heavy, + MPI-bearing stage. + """ + input: + # ``params.py`` is a *tracked* output of ``im_init``, so this one input + # supplies the im_init -> im_pipeline edge. The ``cfis`` symlink the + # shell reads (via {RUN_JOB}) is created by that same im_init shell block + # as an *untracked* side effect -- no rule declares it as an output + # (snakemake will not track a symlink/directory output). Declaring it an + # input here therefore asked the DAG for a file no rule produces: on a + # fresh grids_base it aborted the build with MissingInputException before + # any job ran. It is safe to drop -- cfis exists whenever params does, + # since im_init stages both together. + params=f"{GRIDS_BASE}/{{sim}}/params.py", + output: + done=touch(f"{GRIDS_BASE}/{{sim}}/logs/pipeline_{{tile}}.done"), + params: + run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", + psf=IMSIM["psf_model"], + n_smp=IMSIM["n_smp"], + resources: + mem_mb=16000, + runtime=720, + shell: + "cd {params.run_dir} && " + "{EXEC_PIPELINE} bash {RUN_JOB} " + "-e {wildcards.tile} -t image_sims -j {JOB_MASK} " + "-p {params.psf} -N {params.n_smp}" + + +rule im_merge: + """Merge per-tile ShapePipe catalogues into final_cat_{sim}.hdf5. + + ``create_final_cat.py`` lives in the ShapePipe repo/image; run in image_sims + mode (``-I``) it walks the per-tile output under the run directory. + """ + input: + tiles=expand( + f"{GRIDS_BASE}/{{{{sim}}}}/logs/pipeline_{{tile}}.done", + tile=TILE_IDS, + ), + output: + cat=f"{GRIDS_BASE}/{{sim}}/final_cat_{{sim}}.hdf5", + params: + run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", + shell: + "cd {params.run_dir} && " + "{EXEC_PIPELINE} python {CREATE_FINAL_CAT} " + "-I -m final_cat_{wildcards.sim}.hdf5 -i .. " + "-p cfis/final_cat.param -P {wildcards.sim} " + "-o n_tiles_final.txt -v" + + +rule im_extract: + """Extract the comprehensive ngmix catalogue (sp_validation stage). + + ``extract_info.py`` reads ``params.py`` from cwd and the merged catalogue, + writing ``shape_catalog_comprehensive_{shape}``. + """ + input: + cat=f"{GRIDS_BASE}/{{sim}}/final_cat_{{sim}}.hdf5", + params=f"{GRIDS_BASE}/{{sim}}/params.py", + output: + cat=f"{GRIDS_BASE}/{{sim}}/shape_catalog_comprehensive_{SHAPE}.fits", + params: + run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", + shell: + "cd {params.run_dir} && {EXEC} python {EXTRACT_INFO}" + + +rule im_calibrate: + """Calibrate and cut the comprehensive catalogue (sp_validation stage). + + ``calibrate_comprehensive_cat.py`` reads ``config_mask.yaml`` from cwd, + applies the metacal calibration and selection, and writes + ``shape_catalog_cut_{shape}.fits``. + """ + input: + cat=f"{GRIDS_BASE}/{{sim}}/shape_catalog_comprehensive_{SHAPE}.fits", + mask=f"{GRIDS_BASE}/{{sim}}/config_mask.yaml", + output: + cat=f"{GRIDS_BASE}/{{sim}}/shape_catalog_cut_{SHAPE}.fits", + params: + run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", + shell: + "cd {params.run_dir} && " + "{EXEC} python {CALIBRATE} -s calibrate" + + +rule im_mbias: + """Multiplicative/additive shear bias from the calibrated grids. + + Produces the workflow's headline artifact, ``m_bias_results.yaml``. The + injected shear (``shear_amplitude`` and the branch map) comes from + ``manifest.yaml`` alone -- no literal amplitude here or in config.yaml. The + generated ``m_bias_config.yaml`` carries the manifest's ``branches`` and + ``pairs``, so the estimator's sim list and pairing are the campaign's, not a + hard-coded default. + """ + input: + manifest=MANIFEST, + cats=expand( + f"{GRIDS_BASE}/{{sim}}/shape_catalog_cut_{SHAPE}.fits", sim=SIMS + ), + output: + results=f"{GRIDS_BASE}/results/m_bias_results.yaml", + params: + cfg=f"{GRIDS_BASE}/results/m_bias_config.yaml", + grids_base=GRIDS_BASE, + num=NUM, + cat_name=f"shape_catalog_cut_{SHAPE}.fits", + sif=SIF, + sif_pipeline=SIF_PIPELINE, + shapepipe_repo=SHAPEPIPE_REPO, + sp_validation_repo=SPV_REPO, + # Science knobs, read bare from the run config (no default here). + match_radius_deg=IMSIM["match_radius_deg"], + w_cols=IMSIM["w_cols"], + n_bootstrap=IMSIM["n_bootstrap"], + pair_match=IMSIM["pair_match"], + bootstrap_seed=IMSIM["bootstrap_seed"], + run: + import hashlib + import re + import subprocess + + import yaml + + with open(input.manifest) as fh: + manifest = yaml.safe_load(fh) + + def _git(repo, *args): + """Read a git fact from ``repo``; ``None`` if it is not a checkout.""" + try: + return subprocess.run( + ["git", "-C", repo, *args], + capture_output=True, + text=True, + check=True, + ).stdout.strip() + except (subprocess.CalledProcessError, FileNotFoundError): + return None + + def _sif_revision(sif_path): + """GHCR revision baked into the SIF's OCI labels. + + A plain-text scan of the image file (login-safe: no exec, no + container start), reading org.opencontainers.image.revision -- the + source commit GHCR built the image from. ``None`` if absent. + """ + try: + with open(sif_path, "rb") as fh: + blob = fh.read() + except OSError: + return None + m = re.search( + rb'org\.opencontainers\.image\.revision"?[:=]"?([0-9a-f]{7,40})', + blob, + ) + return m.group(1).decode() if m else None + + # Manifest hash: sha256 of the exact bytes im_manifest wrote, so the + # result records which injected-shear facts it was computed against. + with open(input.manifest, "rb") as fh: + manifest_sha256 = hashlib.sha256(fh.read()).hexdigest() + + provenance = { + "manifest_sha256": manifest_sha256, + "sp_validation": { + "branch": _git(params.sp_validation_repo, "rev-parse", "--abbrev-ref", "HEAD"), + "commit": _git(params.sp_validation_repo, "rev-parse", "HEAD"), + }, + "shapepipe": { + "branch": _git(params.shapepipe_repo, "rev-parse", "--abbrev-ref", "HEAD"), + "commit": _git(params.shapepipe_repo, "rev-parse", "HEAD"), + }, + "containers": { + "sif": params.sif, + "ghcr_revision": _sif_revision(params.sif), + "sif_pipeline": params.sif_pipeline, + "ghcr_revision_pipeline": _sif_revision(params.sif_pipeline), + }, + } + + os.makedirs(os.path.dirname(output.results), exist_ok=True) + # Emit *every* key the estimator requires -- pair_match and + # bootstrap_seed included. Requiring a key without emitting it would + # be a KeyError at run time, so the generated config is the complete + # contract between rule and estimator. ``provenance`` rides along as a + # top-level block: the compute script copies it verbatim into the output + # results yaml, so a result file is self-describing (which manifest, + # which repo commits, which container built the number). + mbias_cfg = { + "grids_dir": params.grids_base, + "num": params.num, + "catalog_name": params.cat_name, + # Injected shear: from the manifest, the single source of truth. + "shear_amplitude": manifest["shear_amplitude"], + "branches": list(manifest["branches"]), + "pairs": manifest["pairs"], + "match_radius_deg": params.match_radius_deg, + "w_cols": list(params.w_cols), + "pair_match": params.pair_match, + "n_bootstrap": params.n_bootstrap, + "bootstrap_seed": params.bootstrap_seed, + "results_dir": os.path.dirname(output.results), + "output_path": output.results, + "provenance": provenance, + } + with open(params.cfg, "w") as fh: + yaml.safe_dump(mbias_cfg, fh) + shell( + "{EXEC} python {COMPUTE_M_BIAS} -c {params.cfg} -v" + ) diff --git a/workflow/scripts/im_build_manifest.py b/workflow/scripts/im_build_manifest.py new file mode 100644 index 00000000..6f0a0faf --- /dev/null +++ b/workflow/scripts/im_build_manifest.py @@ -0,0 +1,244 @@ +#!/usr/bin/env python +"""Build the image-simulation campaign manifest from the sims' own records. + +This is the head of the image-simulation DAG: it reads the injected shear that +the sim campaign recorded for each requested branch, cross-checks it against the +branch's *name*, and writes a single ``manifest.yaml`` that every downstream +stage reads. The point is one home for the injected-shear facts -- amplitude, +per-branch ``(g1, g2)``, the reference branch, and the ``+/-`` pairing -- so no +literal amplitude or hard-coded branch list survives anywhere else. + +The source of truth is each branch's ``basic_info.txt``, written by the sim +campaign. The one line this parser needs looks like:: + + g_cosmic = 0.025 0.0 + +i.e. the literal key ``g_cosmic``, run-together whitespace, an ``=``, then the +two injected shear components ``g1 g2`` as space-separated floats (sign as a +leading ``-``; ``0.0`` for an un-sheared component). We parse *only* that line, +by stdlib string ops -- no YAML/regex dependency -- so the script runs inside +the container with nothing but the standard library. + +Branch names follow the ``1{X}2{Y}`` convention: the character after ``1`` is +the g1 sign, the character after ``2`` is the g2 sign, each one of ``p`` (+), +``m`` (-), ``z`` (0). So ``1p2z`` injects ``(+|g|, 0)``, ``1z2m`` injects +``(0, -|g|)``, ``1z2z`` is the un-sheared reference. The suffix +(``_grid_1`` etc.) is appended by the workflow and is not part of the sign code. + +Validation (all fail-loud, each message naming the offending file and field): + +* every branch's parsed ``(g1, g2)`` sign/axis matches its name's sign code; +* the reference branch parses to exactly ``(0, 0)``; +* the derived ``|g|`` (the single nonzero magnitude of a sheared branch) is + equal across all four sheared branches -- one injected amplitude for the + whole campaign. + +Only the branches this run requests are read and validated. +""" + +import argparse +import os +import sys + +import yaml + +# Branch-name sign code: the character after "1" (g1) and after "2" (g2). +_SIGN = {"p": +1, "m": -1, "z": 0} + + +def die(msg): + """Abort with a clear, prefixed message on stderr.""" + sys.exit(f"im_build_manifest: {msg}") + + +def basic_info_path(input_sims_base, branch): + """Path to a branch's basic_info.txt (the sim campaign's own record).""" + return os.path.join(input_sims_base, branch, "basic_info.txt") + + +def parse_g_cosmic(path): + """Parse the ``g_cosmic = g1 g2`` line from a basic_info.txt file. + + Returns ``(g1, g2)`` as floats. Fails loud, naming the file, if the line + is absent, malformed, or does not carry exactly two float components. + """ + if not os.path.isfile(path): + die(f"basic_info.txt not found: {path}") + + with open(path) as fh: + lines = fh.readlines() + + matches = [ln for ln in lines if ln.split("=", 1)[0].strip() == "g_cosmic"] + if not matches: + die(f"no 'g_cosmic' line in {path}") + if len(matches) > 1: + die(f"multiple 'g_cosmic' lines in {path}") + + rhs = matches[0].split("=", 1)[1].split() + if len(rhs) != 2: + die( + f"'g_cosmic' in {path}: expected two components 'g1 g2', " + f"got {len(rhs)}: {matches[0].strip()!r}" + ) + try: + return float(rhs[0]), float(rhs[1]) + except ValueError: + die(f"'g_cosmic' in {path}: components not floats: {matches[0].strip()!r}") + + +def sign_code(branch): + """Extract the ``(g1_sign, g2_sign)`` code from a ``1{X}2{Y}`` branch name. + + ``branch`` is the bare sign code (e.g. ``1p2z``), suffix already stripped. + Fails loud if the name does not match the convention. + """ + if ( + len(branch) != 4 + or branch[0] != "1" + or branch[2] != "2" + or branch[1] not in _SIGN + or branch[3] not in _SIGN + ): + die( + f"branch name {branch!r} does not match the 1{{X}}2{{Y}} convention " + f"(X, Y each one of p/m/z)" + ) + return _SIGN[branch[1]], _SIGN[branch[3]] + + +def build_manifest(input_sims_base, sims_type, num, branches): + """Parse + validate every requested branch, return the manifest dict. + + ``branches`` are bare sign codes (``1z2z``, ``1p2z``, ...); the on-disk + directory name is ``{branch}{suffix}`` with ``suffix`` derived from + ``sims_type``/``num`` exactly as the workflow builds it. + """ + suffix = f"_{sims_type}_{num}" if sims_type == "grid" else f"_{num}" + + parsed = {} # branch -> (g1, g2) from basic_info.txt + for branch in branches: + dirname = f"{branch}{suffix}" + g1, g2 = parse_g_cosmic(basic_info_path(input_sims_base, dirname)) + s1, s2 = sign_code(branch) + + # Sign/axis must agree with the name: a component is nonzero iff its + # sign code is nonzero, and its sign matches. + for comp, (g, s) in enumerate(((g1, s1), (g2, s2)), start=1): + path = basic_info_path(input_sims_base, dirname) + if s == 0 and g != 0.0: + die( + f"{path}: branch {branch!r} names g{comp} un-sheared (z) but " + f"g_cosmic gives g{comp} = {g}" + ) + if s != 0 and (g == 0.0 or (g > 0) != (s > 0)): + die( + f"{path}: branch {branch!r} names g{comp} sign {'+' if s > 0 else '-'} " + f"but g_cosmic gives g{comp} = {g}" + ) + parsed[branch] = (g1, g2) + + # Reference branch: the one whose name codes (0, 0). Must exist and be (0,0). + refs = [b for b in branches if sign_code(b) == (0, 0)] + if len(refs) != 1: + die( + f"expected exactly one reference branch (name code 1z2z) among " + f"{branches}, found {refs}" + ) + reference = refs[0] + if parsed[reference] != (0.0, 0.0): + die( + f"{basic_info_path(input_sims_base, f'{reference}{suffix}')}: reference " + f"branch {reference!r} must inject (0, 0), got {parsed[reference]}" + ) + + # Derived amplitude: the single nonzero magnitude of each sheared branch, + # cross-checked equal across all four. + amplitudes = {} # branch -> |g| + for branch in branches: + if branch == reference: + continue + g1, g2 = parsed[branch] + amplitudes[branch] = abs(g1) if g1 != 0.0 else abs(g2) + distinct = sorted(set(amplitudes.values())) + if len(distinct) != 1: + die( + "injected |g| differs across sheared branches (must be one campaign " + f"amplitude): {amplitudes} " + f"[files under {input_sims_base}/{suffix}/basic_info.txt]" + ) + shear_amplitude = distinct[0] + + # Pairs: (+component, -component) for each sheared axis, in branch order so + # the estimator's per-pair processing order is stable. + plus = {} # component (0/1) -> branch with +|g| on that component + minus = {} + for branch in branches: + if branch == reference: + continue + s1, s2 = sign_code(branch) + comp = 0 if s1 != 0 else 1 + (plus if (s1 or s2) > 0 else minus)[comp] = branch + pairs = [ + {"plus": plus[comp], "minus": minus[comp], "component": comp} + for comp in sorted(set(plus) & set(minus)) + ] + + return { + "input_sims_base": input_sims_base, + "sims_type": sims_type, + "num": num, + "shear_amplitude": shear_amplitude, + "reference": reference, + # Branch order preserved (dict insertion order round-trips through + # yaml.safe_dump with sort_keys=False) so downstream load order is fixed. + "branches": { + branch: {"g1": parsed[branch][0], "g2": parsed[branch][1]} + for branch in branches + }, + "pairs": pairs, + } + + +def parse_args(): + p = argparse.ArgumentParser(description=__doc__) + p.add_argument( + "--input-sims-base", + required=True, + help="root under which each branch dir holds basic_info.txt", + ) + p.add_argument( + "--sims-type", required=True, help="'grid' -> _grid_{num} suffix, else _{num}" + ) + p.add_argument("--num", required=True, type=int, help="run number") + p.add_argument( + "--branch", + required=True, + action="append", + dest="branches", + help="bare branch sign code (1z2z, 1p2z, ...); repeatable", + ) + p.add_argument("-o", "--output", required=True, help="manifest.yaml output path") + return p.parse_args() + + +def main(): + args = parse_args() + manifest = build_manifest( + args.input_sims_base, args.sims_type, args.num, args.branches + ) + os.makedirs(os.path.dirname(os.path.abspath(args.output)), exist_ok=True) + with open(args.output, "w") as fh: + yaml.safe_dump(manifest, fh, sort_keys=False) + print(f"im_build_manifest: wrote {args.output}") + print(f" shear_amplitude = {manifest['shear_amplitude']}") + print(f" reference = {manifest['reference']}") + print(f" branches = {list(manifest['branches'])}") + for pair in manifest["pairs"]: + print( + f" pair g{pair['component'] + 1}: {pair['plus']} (+) / {pair['minus']} (-)" + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/workflow/scripts/im_compose_mask.py b/workflow/scripts/im_compose_mask.py new file mode 100644 index 00000000..fb602a11 --- /dev/null +++ b/workflow/scripts/im_compose_mask.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python +"""Compose the image-sim mask/calibration config from base + declared overlay. + +The image-sim calibration differs from the data calibration in only a handful +of places (input path, dropped coverage cuts, unweighted global response, +additive-bias off). Instead of maintaining a second full copy of the config -- +which can silently drift from the base it was branched from -- we keep the +*data* config (``mask_v1.X.9.yaml``) as the single home for the shared cuts and +declare the sim-specific delta in an overlay +(``mask_v1.X.9_im_sim.overlay.yaml``). This script applies the overlay to the +base and emits the resolved config, which is byte-for-byte the committed runtime +file ``mask_v1.X.9_im_sim.yaml``. A test locks that equality, so the declaration +and the runtime file cannot diverge. + +The overlay is a list of block operations on the base *text* (not on parsed +YAML), so the resolved file preserves the base's exact formatting and comments +-- the property that makes byte-identity with a hand-maintained runtime file +achievable, and the delta legible as a plain diff. Each op: + +* ``drop:`` remove a verbatim block of base text; +* ``replace:`` / ``with:`` swap a verbatim block for new text. + +Every anchor (the ``drop`` block, or a ``replace`` block) must occur **exactly +once** in the base -- zero or multiple matches is a hard error, so an overlay +that has fallen out of sync with the base fails loudly instead of composing +something wrong. ``why`` is prose for the human reader and is ignored here. + +Stdlib + PyYAML only, so it runs inside the sp_validation container with nothing +extra. +""" + +import argparse +import os +import sys + +import yaml + + +def die(msg): + """Abort with a clear, prefixed message on stderr.""" + sys.exit(f"im_compose_mask: {msg}") + + +def _apply_once(text, anchor, replacement, *, kind, i): + """Replace the single occurrence of ``anchor`` in ``text`` with ``replacement``. + + ``anchor`` must occur exactly once; anything else (missing, or ambiguous) + means the overlay no longer matches the base and is a hard error naming the + offending op. + """ + n = text.count(anchor) + if n != 1: + die( + f"op {i} ({kind}): anchor block occurs {n} time(s) in the base, " + "expected exactly 1 -- the overlay is out of sync with the base.\n" + f"--- anchor ---\n{anchor}\n--------------" + ) + return text.replace(anchor, replacement) + + +def compose(base_text, overlay): + """Apply ``overlay['ops']`` to ``base_text`` and return the resolved text.""" + text = base_text + for i, op in enumerate(overlay["ops"]): + if "drop" in op: + text = _apply_once(text, op["drop"], "", kind="drop", i=i) + elif "replace" in op: + if "with" not in op: + die(f"op {i} (replace): missing 'with:' block") + text = _apply_once(text, op["replace"], op["with"], kind="replace", i=i) + else: + die(f"op {i}: needs a 'drop:' or 'replace:'/'with:' block") + return text + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + ap.add_argument( + "overlay", + help="overlay yaml declaring base + block ops " + "(e.g. mask_v1.X.9_im_sim.overlay.yaml)", + ) + ap.add_argument( + "-o", + "--output", + help="write resolved config here; default: stdout", + ) + args = ap.parse_args(argv) + + with open(args.overlay) as fh: + overlay = yaml.safe_load(fh) + + # The base path is stated in the overlay, relative to the overlay's own dir, + # so the pair travels together (both live in config/calibration/). + base_path = os.path.join(os.path.dirname(args.overlay), overlay["base"]) + with open(base_path) as fh: + base_text = fh.read() + + resolved = compose(base_text, overlay) + + if args.output: + with open(args.output, "w") as fh: + fh.write(resolved) + else: + sys.stdout.write(resolved) + + +if __name__ == "__main__": + main()