From 72b4159276e096b967d870b94455d625d0b66697 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Tue, 16 Jun 2026 13:46:54 +0200 Subject: [PATCH 01/15] Using shapepipe image_sims docker as base --- Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Dockerfile b/Dockerfile index b0c0bad1..b0886a4c 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:image_sims-dev RUN apt-get update -y --quiet --fix-missing && \ apt-get dist-upgrade -y --quiet --fix-missing && \ From 75b51b1dabd055a62bab8cb4f44173752349a807 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Tue, 16 Jun 2026 13:50:02 +0200 Subject: [PATCH 02/15] Using shapepipe image_sims docker as base --- Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Dockerfile b/Dockerfile index b0886a4c..0fc67d6f 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,5 +1,5 @@ # Development image with more bells and whistles -FROM ghcr.io/cosmostat/shapepipe:image_sims-dev +FROM ghcr.io/cosmostat/shapepipe:im_sims RUN apt-get update -y --quiet --fix-missing && \ apt-get dist-upgrade -y --quiet --fix-missing && \ From 003058c689db0abaf461325c417b3bf4549bce48 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 29 Jun 2026 08:39:22 +0200 Subject: [PATCH 03/15] sp_val on image sims; more hdf5 support, n_tiles in header --- config/calibration/mask_v1.X.4_im_sim.yaml | 62 +++++ config/calibration/mask_v1.X.9_im_sim.yaml | 71 ++++++ notebooks/calibrate_comprehensive_cat.py | 2 +- notebooks/extract_info.py | 113 +++++---- notebooks/params.py | 13 +- scripts/compute_m_bias_image_sims.py | 263 +++++++++++++++++++++ src/sp_validation/cat.py | 142 +++++++---- src/sp_validation/image_sims.py | 141 +++++++++++ src/sp_validation/io.py | 30 +-- 9 files changed, 701 insertions(+), 136 deletions(-) create mode 100644 config/calibration/mask_v1.X.4_im_sim.yaml create mode 100644 config/calibration/mask_v1.X.9_im_sim.yaml create mode 100644 scripts/compute_m_bias_image_sims.py create mode 100644 src/sp_validation/image_sims.py 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..f8045551 --- /dev/null +++ b/config/calibration/mask_v1.X.4_im_sim.yaml @@ -0,0 +1,62 @@ +# 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 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..dd8734fd --- /dev/null +++ b/config/calibration/mask_v1.X.9_im_sim.yaml @@ -0,0 +1,71 @@ +# 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.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: 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_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: 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 + global_R_weight: w diff --git a/notebooks/calibrate_comprehensive_cat.py b/notebooks/calibrate_comprehensive_cat.py index 404a7bec..bb529c78 100644 --- a/notebooks/calibrate_comprehensive_cat.py +++ b/notebooks/calibrate_comprehensive_cat.py @@ -100,7 +100,7 @@ rel_size_max=cm["gal_rel_size_max"], size_corr_ell=cm["gal_size_corr_ell"], sigma_eps=cm["sigma_eps_prior"], - global_R_weight=cm["global_R_weight"], + global_R_weight=cm.get("global_R_weight", None), col_2d=False, verbose=True, ) diff --git a/notebooks/extract_info.py b/notebooks/extract_info.py index a4e33451..61cd6280 100644 --- a/notebooks/extract_info.py +++ b/notebooks/extract_info.py @@ -33,6 +33,7 @@ import os import numpy as np from astropy.io import fits +import h5py from cs_util import canfar from sp_validation.io import * @@ -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 # @@ -100,6 +104,8 @@ n_found, n_missing = missing_tiles(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: @@ -139,47 +145,35 @@ 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_ELL_PSFo_NOSHEAR'][:,0][ind_star] != -10) -) + m_star = ( + (dd['FLAGS'][ind_star] == 0) + & (dd['IMAFLAGS_ISO'][ind_star] == 0) + & (dd['NGMIX_MCAL_FLAGS'][ind_star] == 0) + & (dd['NGMIX_ELL_PSFo_NOSHEAR'][:,0][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_ell, - n_star_tot, - stats_file, - verbose=verbose -) -# - + ra_star, dec_star, g_star_psf = spv_cat.match_subsample( + dd, + ind_star, + m_star, + [col_name_ra, col_name_dec], + key_PSF_ell, + n_star_tot, + stats_file, + verbose=verbose + ) -# MKDEBUG: Moved from end of this script # ### 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], -) - -#### Refine: Match to SPREAD_CLASS samples -if "SPREAD_CLASS" in dd.dtype.names: - spv_cat.match_spread_class(dd, ind_star, m_star, stats_file, len(ra_star), verbose=verbose) -else: - print_stats("No SPREAD_CLASS in input, skipping star-gal matching", stats_file, verbose=verbose) +if star_cat_path: + 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 @@ -248,8 +242,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, @@ -260,6 +255,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 @@ -270,18 +276,10 @@ else: if verbose: print("Continuing with selection and calibration") + os.makedirs(os.path.join(output_dir, plot_dir), exist_ok=True) # ## 4. Select galaxies -# #### Common flags and cuts -# First, set cuts common to ngmix and galsim: -# - spread model: select objects well larger than the PSF -# - magnitude: cut galaxies that are too faint (= too noisy, likely to be -# artefacts), and too bright (might be too large for postage stamp) -# - flags: cut objects that were flagged as invalid or masked -# - n_epoch: select objects observed on at leatst one epoch (for safety, -# to avoid potential errors with empty data) - # + cut_common = classification_galaxy_base( dd, @@ -620,12 +618,13 @@ # ## 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 @@ -950,7 +949,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, @@ -989,7 +988,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, @@ -1019,5 +1018,5 @@ 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/notebooks/params.py b/notebooks/params.py index 45478543..cac7e191 100644 --- a/notebooks/params.py +++ b/notebooks/params.py @@ -70,16 +70,12 @@ ## Output paths ### Output base directory -output_dir = f'{data_dir}/sp_output' +output_dir = f'{data_dir}' -### Galaxy shape catalogue base name. -### Will be appended by -### - '_{sh}.fits' for the basic catalogue -### - 'extended_{sh}.fits' for the extended catalogue -output_shape_cat_base= f'{output_dir}/shape_catalog' +### Galaxy shape catalogue base name +output_shape_cat_base = f'{output_dir}/shape_catalog' ### PSF output catalogue base name. -### Will be appended by '_{sh}.fits' output_PSF_cat_base = f'{output_dir}/psf_catalog' ### File for found tile IDs @@ -107,6 +103,9 @@ ## Output +### Output file format extension: '.fits' or '.hdf5' +output_format = '.hdf5' + ### 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_Tpsf_NOSHEAR"] diff --git a/scripts/compute_m_bias_image_sims.py b/scripts/compute_m_bias_image_sims.py new file mode 100644 index 00000000..d58ca0a5 --- /dev/null +++ b/scripts/compute_m_bias_image_sims.py @@ -0,0 +1,263 @@ +#!/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 sys +import os +import argparse +import yaml +import numpy as np + +# Configure matplotlib for non-interactive backend +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from collections import OrderedDict + +from sp_validation.image_sims import ImageSimMBias + + +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 cumulative m/c bias tracking file.""" + if os.path.isfile(cumulative_path): + with open(cumulative_path) as f: + cumulative = yaml.safe_load(f) or {} + else: + cumulative = {} + + # Check if this n_tiles already exists + if str(n_tiles) in cumulative: + return False + + # Convert numpy types to Python floats for clean YAML + clean_results = {} + for key, val in results.items(): + if isinstance(val, np.ndarray): + clean_results[key] = float(val.item()) if val.size == 1 else val.tolist() + elif isinstance(val, (np.integer, np.floating)): + clean_results[key] = float(val) + else: + clean_results[key] = val + + cumulative[str(n_tiles)] = clean_results + with open(cumulative_path, 'w') as f: + yaml.dump(cumulative, f, default_flow_style=False) + return True + + +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) + + results = mb.run(verbose=True) + + print() + print("=" * 40) + print(" Results") + print("=" * 40) + print(f" m1 = {results['m1']:+.4f} +-{results['m1_err']:.4f}") + print(f" c1 = {results['c1']:+.4f} +-{results['c1_err']:.4f}") + print(f" m2 = {results['m2']:+.4f} +-{results['m2_err']:.4f}") + print(f" c2 = {results['c2']:+.4f} +-{results['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") + + with open(out_path, "w") as f: + yaml.dump(results, 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\n") + f.write(f"m1 = {results['m1']:+.6f} ± {results['m1_err']:.6f}\n") + f.write(f"c1 = {results['c1']:+.6f} ± {results['c1_err']:.6f}\n\n") + f.write(f"m2 = {results['m2']:+.6f} ± {results['m2_err']:.6f}\n") + f.write(f"c2 = {results['c2']:+.6f} ± {results['c2_err']:.6f}\n\n") + f.write("Errors computed via bootstrap resampling (n=500 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) + if added: + print(f"\nAdded n_tiles={args.n_tiles} to {cumulative_path}") + # Generate plots after update + try: + plot_convergence(cumulative_path, results_dir) + except Exception as e: + print(f"Warning: could not generate convergence plots: {e}", file=sys.stderr) + else: + print(f"\nn_tiles={args.n_tiles} already in {cumulative_path}") + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/sp_validation/cat.py b/src/sp_validation/cat.py index 20342828..75125559 100644 --- a/src/sp_validation/cat.py +++ b/src/sp_validation/cat.py @@ -344,6 +344,35 @@ def match_spread_class(dd, ind, mask, stats_file, n_ref, verbose=False): io.print_stats(msg, stats_file, verbose=verbose) +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. @@ -583,55 +612,76 @@ 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 add_header: - primary_header.update(add_header) - - primary_header = cat.write_header_info_sp( - primary_header, - software_name="sp_validation", - software_version=__version__, - author=getpass.getuser(), - ) - - 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") - - primary_header["w"] = "DES weight" - - if sigma_epsilon is not None: - primary_header["sig_eps"] = (sigma_epsilon, "Shape noise RMS") + if ext in (".hdf5", ".hdf", ".h5"): + # Build flat list of (name, 1d-array) pairs, splitting 2D columns + fields = [] + for col, descr 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 alpha_leakage: - primary_header["alpha"] = ( - alpha_leakage, - "Mean scale-dependent PSF leakage", + 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 HDU with information in header + primary_header = fits.Header() + if add_header: + primary_header.update(add_header) + primary_header = cat.write_header_info_sp( + primary_header, + software_name="sp_validation", + software_version=__version__, + author=getpass.getuser(), ) - - primary_hdu = fits.PrimaryHDU(header=primary_header) - - # Final file - hdu_list = fits.HDUList([primary_hdu, table_hdu]) - - hdu_list.writeto(output_path, overwrite=True) + 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") + 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) + hdu_list = fits.HDUList([primary_hdu, table_hdu]) + hdu_list.writeto(output_path, overwrite=True) def write_galaxy_cat(output_path, ra, dec, tile_id): diff --git a/src/sp_validation/image_sims.py b/src/sp_validation/image_sims.py new file mode 100644 index 00000000..4d09aa59 --- /dev/null +++ b/src/sp_validation/image_sims.py @@ -0,0 +1,141 @@ +"""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.cat import match_catalogs_radec + + +# Shear component for each simulation pair +_PAIRS = [ + ("1p2z", "1m2z", 0), # g1 component, index 0 → e1 + ("1z2p", "1z2m", 1), # g2 component, index 1 → e2 +] + + +def _load_cat(path, e_col, w_col): + """Load RA, Dec, ellipticity component and weight from a FITS catalogue.""" + with fits.open(path) as hdul: + data = hdul[1].data + return { + "ra": data["RA"].copy(), + "dec": data["Dec"].copy(), + "e1": data["e1"].copy(), + "e2": data["e2"].copy(), + "w": data[w_col].copy(), + } + + +class ImageSimMBias: + """Compute multiplicative and additive shear bias from image simulations. + + 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 + - shear_amplitude : float, input shear |g| (e.g. 0.02) + - match_radius_deg : float, matching radius in degrees + - w_col : str, weight column name (default 'w_des') + - n_bootstrap : int, number of bootstrap resamples for errors + """ + + def __init__(self, config): + self.cfg = config + self.g_in = config["shear_amplitude"] + self.thresh = config.get("match_radius_deg", 0.0002) + self.w_col = config.get("w_col", "w_des") + self.n_boot = config.get("n_bootstrap", 500) + 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["catalog_name"] + sim_names = ["1z2z", "1p2z", "1m2z", "1z2p", "1z2m"] + + for name in sim_names: + path = f"{grids_dir}/{name}_grid_{num}/{cat_name}" + if verbose: + print(f" Loading {path}") + self.cats[name] = _load_cat(path, "e1", self.w_col) + if verbose: + print(f" {len(self.cats[name]['ra'])} objects") + + def _match_to_ref(self, name): + """Return (idx_ref, idx_sim) matched indices between name and 1z2z.""" + ref = self.cats["1z2z"] + sim = self.cats[name] + idx_ref, idx_sim = match_catalogs_radec( + ref["ra"], ref["dec"], + sim["ra"], sim["dec"], + thresh_deg=self.thresh, + ) + return idx_ref, idx_sim + + 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).""" + e_key = f"e{comp + 1}" + + idx_ref_p, idx_p = self._match_to_ref(name_p) + idx_ref_m, idx_m = self._match_to_ref(name_m) + + if verbose: + print(f" {name_p}: {len(idx_p)} matched | {name_m}: {len(idx_m)} matched") + + e_p = self.cats[name_p][e_key][idx_p] + w_p = self.cats[name_p]["w"][idx_p] + e_m = self.cats[name_m][e_key][idx_m] + w_m = self.cats[name_m]["w"][idx_m] + + 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 + + # Bootstrap errors + rng = np.random.default_rng(seed=42) + m_boot = np.empty(self.n_boot) + c_boot = np.empty(self.n_boot) + n_p, n_m = len(e_p), len(e_m) + for i in range(self.n_boot): + ib_p = rng.integers(0, n_p, n_p) + ib_m = rng.integers(0, n_m, n_m) + ep_b = np.average(e_p[ib_p], weights=w_p[ib_p]) + em_b = np.average(e_m[ib_m], weights=w_m[ib_m]) + m_boot[i] = (ep_b - em_b) / (2 * self.g_in) - 1 + c_boot[i] = (ep_b + em_b) / 2 + + return m, np.std(m_boot), c, np.std(c_boot) + + def run(self, verbose=True): + """Compute m and c for both shear components. + + Returns + ------- + dict with keys m1, m1_err, c1, c1_err, m2, m2_err, c2, c2_err + """ + results = {} + for name_p, name_m, comp in _PAIRS: + label = f"g{comp + 1}" + if verbose: + print(f"\n--- {label}: {name_p} / {name_m} ---") + m, m_err, c, c_err = self._m_c_pair(name_p, name_m, comp, verbose=verbose) + results[f"m{comp + 1}"] = m + results[f"m{comp + 1}_err"] = m_err + results[f"c{comp + 1}"] = c + results[f"c{comp + 1}_err"] = c_err + if verbose: + print(f" m{comp+1} = {m:.4f} ± {m_err:.4f}") + print(f" c{comp+1} = {c:.4f} ± {c_err:.4f}") + return results diff --git a/src/sp_validation/io.py b/src/sp_validation/io.py index bc038b99..9c97ef64 100644 --- a/src/sp_validation/io.py +++ b/src/sp_validation/io.py @@ -14,31 +14,11 @@ import numpy as np -def make_out_dirs(output_dir, plot_dir, plot_subdirs, verbose=False): - """Make output directories. - - Create output directories and subdirs - - Parameters - ---------- - plot_dir: string - directory name - plot_subdirs: array of string - subdirectory names - verbose: bool, optional, default=False - verbose output if True - """ - for d in (output_dir, plot_dir): - if not os.path.isdir(d): - if verbose: - print('Creating dir {}'.format(d)) - os.mkdir(d) - for sd in plot_subdirs: - dsd = '{}/{}'.format(plot_dir, sd) - if not os.path.isdir(dsd): - if verbose: - print('Creating dir {}'.format(dsd)) - os.mkdir(dsd) +def make_out_dirs(output_dir, verbose=False): + """Make output directory.""" + os.makedirs(output_dir, exist_ok=True) + if verbose: + print(f'Created dir {output_dir}') def open_stats_file(directory, file_name): From 10baa5267522754ae60b9e0795b0b3f383aa516c Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 3 Jul 2026 23:20:32 +0200 Subject: [PATCH 04/15] feat(image-sims): extract self-contained m/c bias estimator Bring the image-simulation shear-bias core over from the divergent image_sims_val branch as a clean, self-contained unit: - src/sp_validation/image_sims.py: ImageSimMBias reads the five calibrated grid catalogues (reference + g1/g2 +-shear pairs), matches each to the reference by RA/Dec, and returns m1/m2/c1/c2 with bootstrap errors. Its only sp_validation dependency is match_catalogs_radec. - src/sp_validation/catalog.py: add match_catalogs_radec (nearest- neighbour RA/Dec match within a threshold) alongside match_stars2 -- develop had no such helper, and image_sims_val's cat.py is not on develop. - scripts/compute_m_bias_image_sims.py: CLI driver; now casts numpy scalars to plain floats before dumping so m_bias_results.yaml is human-readable rather than !!python/object binary. - config/calibration/mask_v1.X.{4,9}_im_sim.yaml: the image-sim calibration/cut configs the calibrate step consumes. - tests/test_image_sims.py: analytic synthetic-recovery test (injects a known m/c, asserts exact recovery) plus a match_catalogs_radec check. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01QcYFjEWXxyCtUGzQJ7AhLH --- config/calibration/mask_v1.X.4_im_sim.yaml | 62 ++++ config/calibration/mask_v1.X.9_im_sim.yaml | 71 +++++ scripts/compute_m_bias_image_sims.py | 326 +++++++++++++++++++++ src/sp_validation/catalog.py | 29 ++ src/sp_validation/image_sims.py | 144 +++++++++ src/sp_validation/tests/test_image_sims.py | 114 +++++++ 6 files changed, 746 insertions(+) create mode 100644 config/calibration/mask_v1.X.4_im_sim.yaml create mode 100644 config/calibration/mask_v1.X.9_im_sim.yaml create mode 100644 scripts/compute_m_bias_image_sims.py create mode 100644 src/sp_validation/image_sims.py create mode 100644 src/sp_validation/tests/test_image_sims.py 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..f8045551 --- /dev/null +++ b/config/calibration/mask_v1.X.4_im_sim.yaml @@ -0,0 +1,62 @@ +# 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 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..dd8734fd --- /dev/null +++ b/config/calibration/mask_v1.X.9_im_sim.yaml @@ -0,0 +1,71 @@ +# 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.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: 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_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: 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 + global_R_weight: w diff --git a/scripts/compute_m_bias_image_sims.py b/scripts/compute_m_bias_image_sims.py new file mode 100644 index 00000000..23ba4821 --- /dev/null +++ b/scripts/compute_m_bias_image_sims.py @@ -0,0 +1,326 @@ +#!/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 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 cumulative m/c bias tracking file.""" + if os.path.isfile(cumulative_path): + with open(cumulative_path) as f: + cumulative = yaml.safe_load(f) or {} + else: + cumulative = {} + + # Check if this n_tiles already exists + if str(n_tiles) in cumulative: + return False + + # Convert numpy types to Python floats for clean YAML + clean_results = {} + for key, val in results.items(): + if isinstance(val, np.ndarray): + clean_results[key] = float(val.item()) if val.size == 1 else val.tolist() + elif isinstance(val, (np.integer, np.floating)): + clean_results[key] = float(val) + else: + clean_results[key] = val + + cumulative[str(n_tiles)] = clean_results + with open(cumulative_path, "w") as f: + yaml.dump(cumulative, f, default_flow_style=False) + return True + + +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) + + results = mb.run(verbose=True) + + # Cast numpy scalars to plain Python floats so the YAML/text output is + # human-readable (raw numpy scalars serialise as !!python/object binary). + results = { + key: (float(val) if isinstance(val, (np.integer, np.floating)) else val) + for key, val in results.items() + } + + print() + print("=" * 40) + print(" Results") + print("=" * 40) + print(f" m1 = {results['m1']:+.4f} +-{results['m1_err']:.4f}") + print(f" c1 = {results['c1']:+.4f} +-{results['c1_err']:.4f}") + print(f" m2 = {results['m2']:+.4f} +-{results['m2_err']:.4f}") + print(f" c2 = {results['c2']:+.4f} +-{results['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") + + with open(out_path, "w") as f: + yaml.dump(results, 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\n") + f.write(f"m1 = {results['m1']:+.6f} ± {results['m1_err']:.6f}\n") + f.write(f"c1 = {results['c1']:+.6f} ± {results['c1_err']:.6f}\n\n") + f.write(f"m2 = {results['m2']:+.6f} ± {results['m2_err']:.6f}\n") + f.write(f"c2 = {results['c2']:+.6f} ± {results['c2_err']:.6f}\n\n") + f.write("Errors computed via bootstrap resampling (n=500 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) + if added: + print(f"\nAdded n_tiles={args.n_tiles} to {cumulative_path}") + # Generate plots after update + try: + plot_convergence(cumulative_path, results_dir) + except Exception as e: + print( + f"Warning: could not generate convergence plots: {e}", + file=sys.stderr, + ) + else: + print(f"\nn_tiles={args.n_tiles} already in {cumulative_path}") + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/sp_validation/catalog.py b/src/sp_validation/catalog.py index b1b175d9..ca7f1d83 100644 --- a/src/sp_validation/catalog.py +++ b/src/sp_validation/catalog.py @@ -308,6 +308,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. diff --git a/src/sp_validation/image_sims.py b/src/sp_validation/image_sims.py new file mode 100644 index 00000000..703ba128 --- /dev/null +++ b/src/sp_validation/image_sims.py @@ -0,0 +1,144 @@ +"""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 + +# Shear component for each simulation pair +_PAIRS = [ + ("1p2z", "1m2z", 0), # g1 component, index 0 → e1 + ("1z2p", "1z2m", 1), # g2 component, index 1 → e2 +] + + +def _load_cat(path, e_col, w_col): + """Load RA, Dec, ellipticity component and weight from a FITS catalogue.""" + with fits.open(path) as hdul: + data = hdul[1].data + return { + "ra": data["RA"].copy(), + "dec": data["Dec"].copy(), + "e1": data["e1"].copy(), + "e2": data["e2"].copy(), + "w": data[w_col].copy(), + } + + +class ImageSimMBias: + """Compute multiplicative and additive shear bias from image simulations. + + 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 + - shear_amplitude : float, input shear |g| (e.g. 0.02) + - match_radius_deg : float, matching radius in degrees + - w_col : str, weight column name (default 'w_des') + - n_bootstrap : int, number of bootstrap resamples for errors + """ + + def __init__(self, config): + self.cfg = config + self.g_in = config["shear_amplitude"] + self.thresh = config.get("match_radius_deg", 0.0002) + self.w_col = config.get("w_col", "w_des") + self.n_boot = config.get("n_bootstrap", 500) + 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["catalog_name"] + sim_names = ["1z2z", "1p2z", "1m2z", "1z2p", "1z2m"] + + for name in sim_names: + path = f"{grids_dir}/{name}_grid_{num}/{cat_name}" + if verbose: + print(f" Loading {path}") + self.cats[name] = _load_cat(path, "e1", self.w_col) + if verbose: + print(f" {len(self.cats[name]['ra'])} objects") + + def _match_to_ref(self, name): + """Return (idx_ref, idx_sim) matched indices between name and 1z2z.""" + ref = self.cats["1z2z"] + sim = self.cats[name] + idx_ref, idx_sim = match_catalogs_radec( + ref["ra"], + ref["dec"], + sim["ra"], + sim["dec"], + thresh_deg=self.thresh, + ) + return idx_ref, idx_sim + + 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).""" + e_key = f"e{comp + 1}" + + idx_ref_p, idx_p = self._match_to_ref(name_p) + idx_ref_m, idx_m = self._match_to_ref(name_m) + + if verbose: + print( + f" {name_p}: {len(idx_p)} matched | {name_m}: {len(idx_m)} matched" + ) + + e_p = self.cats[name_p][e_key][idx_p] + w_p = self.cats[name_p]["w"][idx_p] + e_m = self.cats[name_m][e_key][idx_m] + w_m = self.cats[name_m]["w"][idx_m] + + 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 + + # Bootstrap errors + rng = np.random.default_rng(seed=42) + m_boot = np.empty(self.n_boot) + c_boot = np.empty(self.n_boot) + n_p, n_m = len(e_p), len(e_m) + for i in range(self.n_boot): + ib_p = rng.integers(0, n_p, n_p) + ib_m = rng.integers(0, n_m, n_m) + ep_b = np.average(e_p[ib_p], weights=w_p[ib_p]) + em_b = np.average(e_m[ib_m], weights=w_m[ib_m]) + m_boot[i] = (ep_b - em_b) / (2 * self.g_in) - 1 + c_boot[i] = (ep_b + em_b) / 2 + + return m, np.std(m_boot), c, np.std(c_boot) + + def run(self, verbose=True): + """Compute m and c for both shear components. + + Returns + ------- + dict with keys m1, m1_err, c1, c1_err, m2, m2_err, c2, c2_err + """ + results = {} + for name_p, name_m, comp in _PAIRS: + label = f"g{comp + 1}" + if verbose: + print(f"\n--- {label}: {name_p} / {name_m} ---") + m, m_err, c, c_err = self._m_c_pair(name_p, name_m, comp, verbose=verbose) + results[f"m{comp + 1}"] = m + results[f"m{comp + 1}_err"] = m_err + results[f"c{comp + 1}"] = c + results[f"c{comp + 1}_err"] = c_err + if verbose: + print(f" m{comp + 1} = {m:.4f} ± {m_err:.4f}") + print(f" c{comp + 1} = {c:.4f} ± {c_err:.4f}") + return results 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..8988347c --- /dev/null +++ b/src/sp_validation/tests/test_image_sims.py @@ -0,0 +1,114 @@ +"""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|``). Objects are matched to +the reference by RA/Dec, then, per component, + + m = ( - ) / (2 |g|) - 1 , c = ( + ) / 2 . + +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_col": "w_des", + "n_bootstrap": 50, + } + 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 From f586466d30cf599458f0b3ac50cd07f6e4272da9 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 3 Jul 2026 23:55:25 +0200 Subject: [PATCH 05/15] feat(image-sims): snakemake orchestration for the m-bias chain Add workflow/rules/image_sims.smk and a standalone entry point that drive the image-simulation validation chain end to end across two containers, the sp_validation-side half of UNIONS-WL/MultiBand_ImSim#1: ShapePipe container: im_pipeline (raw sim tiles -> per-tile cats), im_merge (create_final_cat -> final_cat_{sim}) sp_validation container: im_extract -> comprehensive, im_calibrate -> cut, im_mbias -> m_bias_results.yaml Design: - Native container pattern: every rule sets `container: None` and calls `apptainer exec` explicitly, since the two halves live in different images. Extract/calibrate use the sp_validation image's own self-consistent calibration environment; m-bias injects this branch's extracted core on PYTHONPATH. - Fully parameterised under config["image_sims"] -- container paths, repo roots, data roots, PSF dict, tile list, sim/calibration knobs -- so a fresh user drives it from config alone, no hard-coded clone layout (the "runnable by someone other than Fabian" gap). - Shell blocks avoid f-strings, so snakemake's {input}/{wildcards} placeholders interpolate correctly -- fixing the {input.tiles} proxy- repr bug from the original Snakefile. - params_im_sim.py derives the field name from the run-dir basename, so a single shared template serves every sim. - Included in the main workflow under `if "image_sims" in config`; standalone via workflow/image_sims/Snakefile. Verified on candide: DAG resolves (36 jobs, correct chaining); im_init and im_calibrate run in-container and produce a real cut catalogue; the extracted m-bias core produces a real m2/c2 from grid_2's calibrated ShapePipe catalogues. Full end-to-end from Martin's on-disk intermediates is blocked by data-vintage inconsistencies (see fiber), not the workflow. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01QcYFjEWXxyCtUGzQJ7AhLH --- workflow/Snakefile | 5 + workflow/image_sims/Snakefile | 28 +++ workflow/image_sims/config.yaml | 52 +++++ workflow/image_sims/params_im_sim.py | 220 ++++++++++++++++++++ workflow/rules/image_sims.smk | 297 +++++++++++++++++++++++++++ 5 files changed, 602 insertions(+) create mode 100644 workflow/image_sims/Snakefile create mode 100644 workflow/image_sims/config.yaml create mode 100644 workflow/image_sims/params_im_sim.py create mode 100644 workflow/rules/image_sims.smk 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..2c44a406 --- /dev/null +++ b/workflow/image_sims/Snakefile @@ -0,0 +1,28 @@ +"""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. + + snakemake -s workflow/image_sims/Snakefile \\ + --configfile workflow/image_sims/config.yaml \\ + -j 4 im_mbias + +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..fee24c63 --- /dev/null +++ b/workflow/image_sims/config.yaml @@ -0,0 +1,52 @@ +# Image-simulation m-bias workflow configuration. +# +# Every path the workflow touches is set here: a fresh user edits this file +# and nothing else. Defaults point at the shared candide assets (world- +# readable) used to develop the workflow. + +image_sims: + + # --- containers ------------------------------------------------------- + # ShapePipe image (pipeline + merge) and sp_validation image (extract, + # calibrate, m-bias). Rebuild from the retargeted #766 / this branch when + # newer images are available. + shapepipe_sif: /n17data/mkilbing/shapepipe_im_sims.sif + sp_validation_sif: /n17data/mkilbing/sp_validation_im_sims.sif + # Apptainer bind mounts. /automnt is required when repos/data are + # automounted (candide gotcha); harmless otherwise. + binds: /n17data,/n09data,/home,/automnt + + # --- repositories ----------------------------------------------------- + # Bound into the images; the sp_validation src is put on PYTHONPATH so this + # branch's code (image_sims.py, catalog.match_catalogs_radec) is used. + 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} + num: 1 + # Tiles to process. Give an explicit list (example below) or point + # tile_ids_file at a one-tile-per-line file for a production run. + tile_ids: ["233.293", "237.292", "238.292"] + # tile_ids_file: /path/to/tile_numbers.txt + + # --- calibration ------------------------------------------------------ + shape: ngmix + mask_config: config/calibration/mask_v1.X.9_im_sim.yaml # relative to sp_validation_repo + # ShapePipe cfis configs (final_cat.param etc.); default is + # {shapepipe_repo}/example/cfis_image_sims once #766 lands. + config_dir: /n17data/cdaley/unions/scratch-wf/imsims-run/grids/_cfis_image_sims + psf_model: psfex + n_smp: -1 + + # --- m-bias ----------------------------------------------------------- + shear_amplitude: 0.02 # input |g| of the +/- shear grids + match_radius_deg: 0.0002 + w_col: w_des + n_bootstrap: 500 diff --git a/workflow/image_sims/params_im_sim.py b/workflow/image_sims/params_im_sim.py new file mode 100644 index 00000000..d929751b --- /dev/null +++ b/workflow/image_sims/params_im_sim.py @@ -0,0 +1,220 @@ +""" + +: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 = ".hdf5" + +### 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_Tpsf_NOSHEAR", +] + +## Pre-calibration catalogue, including masked objects and mask flags +add_cols_pre_cal = [ + "TILE_ID", + "NUMBER", + "FLAGS", + "NGMIX_MCAL_FLAGS", + "NGMIX_MOM_FAIL", + "N_EPOCH", + "NGMIX_N_EPOCH", + "NGMIX_ELL_PSFo_NOSHEAR", + "NGMIX_ELL_ERR_NOSHEAR", +] + +### Set flag columns as integer format +add_cols_pre_cal_format = {} +for key in ( + "NUMBER", + "FLAGS", + "NGMIX_MCAL_FLAGS", + "NGMIX_MOM_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", "ELL", "FLUX", "FLUX_ERR", "T", "T_ERR", "Tpsf"] +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/rules/image_sims.smk b/workflow/rules/image_sims.smk new file mode 100644 index 00000000..af6001cc --- /dev/null +++ b/workflow/rules/image_sims.smk @@ -0,0 +1,297 @@ +"""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 (in *its* container) turns the +simulated tiles into per-tile shape catalogues, then sp_validation (in *its* +container) merges, extracts, calibrates and finally measures the +multiplicative/additive shear bias. + +Two containers are wired through one workflow. Every rule sets +``container: None`` and calls ``apptainer exec`` explicitly, because neither +image is the workflow's top-level container and the two halves run in +different images: + +* ShapePipe container -> ``pipeline`` (raw images -> per-tile cats) and + ``merge`` (``create_final_cat`` -> ``final_cat_{sim}.hdf5``). +* sp_validation container -> ``extract`` (-> comprehensive cat), + ``calibrate`` (-> cut cat) and ``m_bias`` (-> ``m_bias_results.yaml``). + +Everything is parameterised under ``config["image_sims"]`` -- container paths, +repository roots, data roots, the PSF dictionary, tile list and sim/calibration +knobs -- so a fresh user drives it from config alone, with no hard-coded clone +layout. The ``PYTHONPATH`` override on the sp_validation exec makes the +*branch* source (``image_sims.py``, ``catalog.match_catalogs_radec``) win 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"] + +# --- containers ----------------------------------------------------------- +SHAPEPIPE_SIF = IMSIM["shapepipe_sif"] +SPV_SIF = IMSIM["sp_validation_sif"] +BINDS = IMSIM.get("binds", "/n17data,/n09data,/home,/automnt") + +# --- repositories (bound into the images; 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 ------------------------------------------------------ +NUM = IMSIM["num"] +SIMS_TYPE = IMSIM.get("sims_type", "grid") +_SUFFIX = f"_{SIMS_TYPE}_{NUM}" if SIMS_TYPE == "grid" else f"_{NUM}" +SIM_BASES = ["1z2z", "1p2z", "1m2z", "1z2p", "1z2m"] +SIMS = [f"{base}{_SUFFIX}" for base in SIM_BASES] + +# --- tiles ---------------------------------------------------------------- +if IMSIM.get("tile_ids"): + TILE_IDS = list(IMSIM["tile_ids"]) +else: + with open(IMSIM["tile_ids_file"]) as fh: + TILE_IDS = [line.strip() for line in fh if line.strip()] + +# --- calibration / m-bias knobs ------------------------------------------ +SHAPE = IMSIM.get("shape", "ngmix") +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.get( + "config_dir", f"{SHAPEPIPE_REPO}/example/cfis_image_sims" +) + +# 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 inside the sp_validation image's own, self-consistent +# calibration environment (its baked scripts + baked sp_validation), so the +# ngmix column grammar matches whatever produced the catalogues. Overridable +# for a host checkout. +EXTRACT_INFO = IMSIM.get("extract_script", "/sp_validation/notebooks/extract_info.py") +CALIBRATE = IMSIM.get( + "calibrate_script", "/sp_validation/notebooks/calibrate_comprehensive_cat.py" +) +# 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 --------------------------------------------- +# ShapePipe stages. MPI/SLURM env vars are stripped so OpenMPI inside the +# image does not try to attach to the host launcher (cf. apptainer_noslurm.sh). +SP_EXEC = f"env -u SLURM_JOBID -u SLURM_JOB_ID -u SLURM_PROCID apptainer exec --bind {BINDS} {SHAPEPIPE_SIF}" +# sp_validation calibration stages: the image's own environment, no override. +SPV_EXEC = f"apptainer exec --bind {BINDS} {SPV_SIF}" +# m-bias stage: inject this branch's src so the extracted ImageSimMBias + +# catalog.match_catalogs_radec win over anything baked into the image. +SPV_EXEC_MBIAS = ( + f"apptainer exec --bind {BINDS} --env PYTHONPATH={SPV_REPO}/src {SPV_SIF}" +) + +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_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}.hdf5", + sim=SIMS, + ), + + +rule im_calibrate_all: + input: + expand( + f"{GRIDS_BASE}/{{sim}}/shape_catalog_cut_{SHAPE}.fits", sim=SIMS + ), + + +# ========================================================================== +# Rules +# ========================================================================== +rule im_init: + """Stage per-sim run directory: params.py, mask config, ShapePipe configs. + + ``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. + """ + output: + params=f"{GRIDS_BASE}/{{sim}}/params.py", + mask=f"{GRIDS_BASE}/{{sim}}/config_mask.yaml", + params: + template=PARAMS_TEMPLATE, + mask_src=lambda wc: os.path.join(SPV_REPO, MASK_CONFIG), + config_dir=CONFIG_DIR, + cfis=lambda wc: f"{GRIDS_BASE}/{wc.sim}/cfis", + shell: + # cfis is a stable read-only symlink to the ShapePipe configs (used by + # merge/extract); created here but not tracked as a directory output, + # which snakemake will not accept for a symlink. + "mkdir -p $(dirname {output.params}) && " + "cp {params.template} {output.params} && " + "ln -sf {params.mask_src} {output.mask} && " + "ln -sfT {params.config_dir} {params.cfis}" + + +rule im_pipeline: + """Run ShapePipe on one simulated tile (ShapePipe container). + + 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=f"{GRIDS_BASE}/{{sim}}/params.py", + cfis=f"{GRIDS_BASE}/{{sim}}/cfis", + output: + done=touch(f"{GRIDS_BASE}/{{sim}}/logs/pipeline_{{tile}}.done"), + params: + run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", + psf=IMSIM.get("psf_model", "psfex"), + n_smp=IMSIM.get("n_smp", -1), + resources: + mem_mb=16000, + runtime=720, + shell: + "cd {params.run_dir} && " + "{SP_EXEC} 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} && " + "{SP_EXEC} 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 container). + + ``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}.hdf5", + params: + run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", + shell: + "cd {params.run_dir} && {SPV_EXEC} python {EXTRACT_INFO}" + + +rule im_calibrate: + """Calibrate and cut the comprehensive catalogue (sp_validation container). + + ``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}.hdf5", + 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} && " + "{SPV_EXEC} python {CALIBRATE} -s calibrate" + + +rule im_mbias: + """Multiplicative/additive shear bias from the five calibrated grids. + + Produces the workflow's headline artifact, ``m_bias_results.yaml``. + """ + input: + 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", + shear_amplitude=IMSIM.get("shear_amplitude", 0.02), + match_radius_deg=IMSIM.get("match_radius_deg", 0.0002), + w_col=IMSIM.get("w_col", "w_des"), + n_bootstrap=IMSIM.get("n_bootstrap", 500), + run: + import yaml + + os.makedirs(os.path.dirname(output.results), exist_ok=True) + mbias_cfg = { + "grids_dir": params.grids_base, + "num": params.num, + "catalog_name": params.cat_name, + "shear_amplitude": params.shear_amplitude, + "match_radius_deg": params.match_radius_deg, + "w_col": params.w_col, + "n_bootstrap": params.n_bootstrap, + "results_dir": os.path.dirname(output.results), + "output_path": output.results, + } + with open(params.cfg, "w") as fh: + yaml.safe_dump(mbias_cfg, fh) + shell( + "{SPV_EXEC_MBIAS} python {COMPUTE_M_BIAS} -c {params.cfg} -v" + ) From 58756eb7a7e809cb97df2827d02a124e691cd85f Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 4 Jul 2026 03:17:15 +0200 Subject: [PATCH 06/15] fix(image-sims): im_init stages raw SKiLLS inputs so im_pipeline runs from scratch im_init only staged params.py / mask / cfis, but ShapePipe's get_images_runner resolves the raw tile+exposure images via $SP_DIR/input_tiles and $SP_DIR/input_exp -- top-level symlinks that neither run_job nor any init script creates. Without them im_pipeline can only run off pre-staged intermediates, not from raw images, which defeats the "runnable from raw SKiLLS" purpose of the workflow. Stage input_tiles -> {input_sims_base}/{sim}/images/SP_tiles and input_exp -> .../SP_exp in im_init, matching the ground-truth run-dir layout of Martin's working image-sims runs. Verified: init produces resolving symlinks to the real CFIS_simu_image-*.fits for all five grid_1 sims. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_013uTmBAtdP5Wvav2M76MgCA --- workflow/rules/image_sims.smk | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/workflow/rules/image_sims.smk b/workflow/rules/image_sims.smk index af6001cc..6d4cb3dc 100644 --- a/workflow/rules/image_sims.smk +++ b/workflow/rules/image_sims.smk @@ -146,11 +146,18 @@ rule im_calibrate_all: # Rules # ========================================================================== rule im_init: - """Stage per-sim run directory: params.py, mask config, ShapePipe configs. + """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. """ output: params=f"{GRIDS_BASE}/{{sim}}/params.py", @@ -159,15 +166,20 @@ rule im_init: template=PARAMS_TEMPLATE, mask_src=lambda wc: os.path.join(SPV_REPO, MASK_CONFIG), 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 is a stable read-only symlink to the ShapePipe configs (used by - # merge/extract); created here but not tracked as a directory output, - # which snakemake will not accept for a symlink. + # 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 {params.template} {output.params} && " "ln -sf {params.mask_src} {output.mask} && " - "ln -sfT {params.config_dir} {params.cfis}" + "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: From 00ed2c422fbc413e6c6e2b50e471451e228dd6ac Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 4 Jul 2026 14:55:46 +0200 Subject: [PATCH 07/15] fix(image-sims): run the m-bias downstream end-to-end on the ngmix-v2 grammar The image-sim merge->extract->calibrate chain had never been run end-to-end on the post-#761 ngmix column grammar; doing so surfaced a cascade of old-grammar and interface defects, all fixed here (validated on 1z2z_grid_1: 2345 -> 1684 post-cut objects, , ~ 0 for the no-shear reference): - params_im_sim.py: migrate add_cols / add_cols_pre_cal / centers to the v2 grammar (NGMIX_ELL_* -> NGMIX_G{1,2}_*, NGMIX_Tpsf_* -> NGMIX_T_PSF_RECONV_*, drop NGMIX_MOM_FAIL, keep NGMIX_MCAL_TYPES_FAIL; IMAFLAGS_ISO omitted -- the sims run no imaging-flag masking). - extract_info.py: guard every star-catalogue-dependent block under `if star_cat_path`. The sims have no star catalogue (star_cat_path=None), so the star match/metacal/PSF-leakage diagnostics -- all downstream of the comprehensive-cat write -- must be skipped rather than crash on undefined `ind_star`. Galaxy diagnostics preserved; mixed gal+star R histograms fall back to galaxy-only. - catalog_builders.py: make CalibrateCat.read_cat FITS-aware. The sims produce a single per-run comprehensive catalogue (FITS, from write_shape_catalog), not the joined multi-patch HDF5 the data path builds; read it directly with dat_ext=None instead of forcing h5py. - mask_v1.X.9_im_sim.yaml: v2 grammar (NGMIX_MOM_FAIL -> NGMIX_MCAL_TYPES_FAIL, NGMIX_ELL_PSFo_NOSHEAR_{0,1} -> NGMIX_G{1,2}_PSF_ORIG_NOSHEAR) and .fits input. - image_sims.smk: comprehensive cat is .fits (write_shape_catalog always writes FITS); default extract/calibrate to the repo checkout and inject PYTHONPATH={sp_validation_repo}/src on every SPV stage, so the branch code (these fixes) wins over the lagging baked container. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01D9atA9VPrg7VMC7oBArny9 --- config/calibration/mask_v1.X.9_im_sim.yaml | 10 +- scripts/calibration/extract_info.py | 232 +++++++++++---------- src/sp_validation/catalog_builders.py | 16 ++ workflow/image_sims/params_im_sim.py | 22 +- workflow/rules/image_sims.smk | 33 +-- 5 files changed, 170 insertions(+), 143 deletions(-) diff --git a/config/calibration/mask_v1.X.9_im_sim.yaml b/config/calibration/mask_v1.X.9_im_sim.yaml index dd8734fd..4e0377f4 100644 --- a/config/calibration/mask_v1.X.9_im_sim.yaml +++ b/config/calibration/mask_v1.X.9_im_sim.yaml @@ -3,7 +3,7 @@ # General parameters (can also given on command line) params: - input_path: shape_catalog_comprehensive_ngmix.hdf5 + input_path: shape_catalog_comprehensive_ngmix.fits cmatrices: False sky_regions: False verbose: True @@ -36,17 +36,17 @@ dat: value: [15, 30] # ngmix flags - - col_name: NGMIX_MOM_FAIL + - col_name: NGMIX_MCAL_TYPES_FAIL label: "ngmix moments failure" kind: equal value: 0 - # invalid PSF ellipticities - - col_name: NGMIX_ELL_PSFo_NOSHEAR_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_ELL_PSFo_NOSHEAR_1 + - col_name: NGMIX_G2_PSF_ORIG_NOSHEAR label: "bad PSF ellipticity comp 2" kind: not_equal value: -10 diff --git a/scripts/calibration/extract_info.py b/scripts/calibration/extract_info.py index 1cbf836f..865505a7 100644 --- a/scripts/calibration/extract_info.py +++ b/scripts/calibration/extract_info.py @@ -149,39 +149,44 @@ # #### Refine: Match to valid, unflagged galaxy sample # + -# Flags to indicate valid star sample +# 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) -) +if star_cat_path: + 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, -) -# - + 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, + ) + # - -# MKDEBUG: Moved from end of this script + # MKDEBUG: Moved from end of this script -# ### 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 @@ -628,23 +633,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 +736,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 +764,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 +783,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 +806,8 @@ x_range, n_bin, out_path, - colors=colors, - linestyles=linestyles, + colors=colors[: len(xs)], + linestyles=linestyles[: len(xs)], ) # - @@ -845,49 +846,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 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/workflow/image_sims/params_im_sim.py b/workflow/image_sims/params_im_sim.py index d929751b..1ba0b5ca 100644 --- a/workflow/image_sims/params_im_sim.py +++ b/workflow/image_sims/params_im_sim.py @@ -122,20 +122,28 @@ "FLUX_APER", "FLUXERR_APER", "NGMIX_T_NOSHEAR", - "NGMIX_Tpsf_NOSHEAR", + "NGMIX_T_PSF_RECONV_NOSHEAR", ] -## Pre-calibration catalogue, including masked objects and mask flags +## 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_MOM_FAIL", + "NGMIX_MCAL_TYPES_FAIL", "N_EPOCH", "NGMIX_N_EPOCH", - "NGMIX_ELL_PSFo_NOSHEAR", - "NGMIX_ELL_ERR_NOSHEAR", + "NGMIX_G1_PSF_ORIG_NOSHEAR", + "NGMIX_G2_PSF_ORIG_NOSHEAR", + "NGMIX_G1_ERR_NOSHEAR", + "NGMIX_G2_ERR_NOSHEAR", ] ### Set flag columns as integer format @@ -144,7 +152,7 @@ "NUMBER", "FLAGS", "NGMIX_MCAL_FLAGS", - "NGMIX_MOM_FAIL", + "NGMIX_MCAL_TYPES_FAIL", "N_EPOCH", "NGMIX_N_EPOCH", ): @@ -156,7 +164,7 @@ # Create key names for metacal information prefix = "NGMIX" suffixes = ["1M", "1P", "2M", "2P", "NOSHEAR"] -centers = ["FLAGS", "ELL", "FLUX", "FLUX_ERR", "T", "T_ERR", "Tpsf"] +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}") diff --git a/workflow/rules/image_sims.smk b/workflow/rules/image_sims.smk index 6d4cb3dc..16c362db 100644 --- a/workflow/rules/image_sims.smk +++ b/workflow/rules/image_sims.smk @@ -74,13 +74,15 @@ CONFIG_DIR = IMSIM.get( # 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 inside the sp_validation image's own, self-consistent -# calibration environment (its baked scripts + baked sp_validation), so the -# ngmix column grammar matches whatever produced the catalogues. Overridable -# for a host checkout. -EXTRACT_INFO = IMSIM.get("extract_script", "/sp_validation/notebooks/extract_info.py") +# 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.get( + "extract_script", f"{SPV_REPO}/scripts/calibration/extract_info.py" +) CALIBRATE = IMSIM.get( - "calibrate_script", "/sp_validation/notebooks/calibrate_comprehensive_cat.py" + "calibrate_script", f"{SPV_REPO}/scripts/calibration/calibrate_comprehensive_cat.py" ) # m-bias is *this branch's* extracted core, injected on PYTHONPATH. COMPUTE_M_BIAS = f"{SPV_REPO}/scripts/compute_m_bias_image_sims.py" @@ -89,13 +91,12 @@ COMPUTE_M_BIAS = f"{SPV_REPO}/scripts/compute_m_bias_image_sims.py" # ShapePipe stages. MPI/SLURM env vars are stripped so OpenMPI inside the # image does not try to attach to the host launcher (cf. apptainer_noslurm.sh). SP_EXEC = f"env -u SLURM_JOBID -u SLURM_JOB_ID -u SLURM_PROCID apptainer exec --bind {BINDS} {SHAPEPIPE_SIF}" -# sp_validation calibration stages: the image's own environment, no override. -SPV_EXEC = f"apptainer exec --bind {BINDS} {SPV_SIF}" -# m-bias stage: inject this branch's src so the extracted ImageSimMBias + -# catalog.match_catalogs_radec win over anything baked into the image. -SPV_EXEC_MBIAS = ( - f"apptainer exec --bind {BINDS} --env PYTHONPATH={SPV_REPO}/src {SPV_SIF}" -) +# sp_validation calibration stages: inject the branch source on PYTHONPATH so +# the repo's sp_validation package (newer than the baked one) wins -- the +# image-sims path depends on branch-only fixes to catalog_builders/extract. +SPV_EXEC = f"apptainer exec --bind {BINDS} --env PYTHONPATH={SPV_REPO}/src {SPV_SIF}" +# m-bias stage uses the same injected environment. +SPV_EXEC_MBIAS = SPV_EXEC JOB_MASK = sum([1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048]) @@ -130,7 +131,7 @@ rule im_merge_all: rule im_extract_all: input: expand( - f"{GRIDS_BASE}/{{sim}}/shape_catalog_comprehensive_{SHAPE}.hdf5", + f"{GRIDS_BASE}/{{sim}}/shape_catalog_comprehensive_{SHAPE}.fits", sim=SIMS, ), @@ -241,7 +242,7 @@ rule im_extract: 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}.hdf5", + cat=f"{GRIDS_BASE}/{{sim}}/shape_catalog_comprehensive_{SHAPE}.fits", params: run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", shell: @@ -256,7 +257,7 @@ rule im_calibrate: ``shape_catalog_cut_{shape}.fits``. """ input: - cat=f"{GRIDS_BASE}/{{sim}}/shape_catalog_comprehensive_{SHAPE}.hdf5", + 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", From ddb95f02cc8037ccef5db4b1b34ad917beb74ae8 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Sat, 4 Jul 2026 18:33:41 +0200 Subject: [PATCH 08/15] im sims: m-bias computation --- scripts/compute_m_bias_image_sims.py | 2 +- scripts/diagnostics_image_sims.py | 175 +++++++++++++++++++++++++++ src/sp_validation/image_sims.py | 22 ++-- 3 files changed, 183 insertions(+), 16 deletions(-) create mode 100644 scripts/diagnostics_image_sims.py diff --git a/scripts/compute_m_bias_image_sims.py b/scripts/compute_m_bias_image_sims.py index d58ca0a5..18421cde 100644 --- a/scripts/compute_m_bias_image_sims.py +++ b/scripts/compute_m_bias_image_sims.py @@ -35,7 +35,7 @@ def get_n_tiles(grids_dir, num): try: import h5py # Count tiles in first sim's final_cat - for sim in ['1z2z_grid', '1m2z_grid', '1p2z_grid', '1z2m_grid', '1z2p_grid']: + for sim in ['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): diff --git a/scripts/diagnostics_image_sims.py b/scripts/diagnostics_image_sims.py new file mode 100644 index 00000000..866220f6 --- /dev/null +++ b/scripts/diagnostics_image_sims.py @@ -0,0 +1,175 @@ +#!/usr/bin/env python +"""Per-sim diagnostics for image simulation catalogues. + +For each of the 5 sheared grid catalogues, produces: + - footprint (RA/Dec scatter) + - ellipticity histograms (e1, e2) + - weight histogram (w_des) + - 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) + +Usage: + diagnostics_image_sims.py -c config.yaml [-v] +""" + +import sys +import argparse +import yaml +import numpy as np +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from astropy.io import fits + + +SIM_NAMES = ["1z2z", "1m2z", "1p2z", "1z2m", "1z2p"] +COLORS = {"1z2z": "black", "1m2z": "C0", "1p2z": "C1", "1z2m": "C2", "1z2p": "C3"} + + +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 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, 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, 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 = d["w_des"] + 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) + fig.suptitle("Ellipticity histograms (w_des weighted)") + return savefig(fig, out_dir, "ellipticity_hist") + + +def plot_weights(cats, out_dir, nbins=50): + fig, ax = plt.subplots(figsize=(8, 5)) + for name, d in cats.items(): + ax.hist(d["w_des"], bins=nbins, density=True, + histtype="step", label=name, color=COLORS[name]) + ax.set_xlabel("w_des") + ax.set_ylabel("normalised count") + ax.legend() + ax.set_title("Weight distribution") + return savefig(fig, out_dir, "weight_hist") + + +def plot_response(cats, 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, density=True, + histtype="step", label=name, color=COLORS[name]) + 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, 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, verbose=True): + print("\n--- Additive bias (weighted mean ellipticity) ---") + results = {} + for name, d in cats.items(): + w = d["w_des"] + 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) + + grids_dir = config["grids_dir"] + num = config["num"] + cat_name = config.get("catalog_name", "shape_catalog_cut_ngmix.fits") + out_dir = config.get("diagnostics_dir", f"{grids_dir}/diagnostics") + + import os + os.makedirs(out_dir, exist_ok=True) + + print(f"Loading catalogues from {grids_dir}...") + cats = {} + for name in SIM_NAMES: + 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, out_dir)}") + print(f" ellipticity -> {plot_ellipticity(cats, out_dir)}") + print(f" weights -> {plot_weights(cats, out_dir)}") + print(f" response -> {plot_response(cats, out_dir)}") + print(f" PSF leakage -> {plot_psf_leakage(cats, out_dir)}") + calculate_additive_bias(cats, verbose=True) + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/sp_validation/image_sims.py b/src/sp_validation/image_sims.py index 4d09aa59..92a8cd7f 100644 --- a/src/sp_validation/image_sims.py +++ b/src/sp_validation/image_sims.py @@ -57,11 +57,11 @@ def __init__(self, config): self.cats = {} def load_catalogs(self, verbose=True): - """Load the 5 sheared and reference catalogues.""" + """Load the 4 sheared catalogues (1p2z, 1m2z, 1z2p, 1z2m).""" grids_dir = self.cfg["grids_dir"] num = self.cfg["num"] cat_name = self.cfg["catalog_name"] - sim_names = ["1z2z", "1p2z", "1m2z", "1z2p", "1z2m"] + sim_names = ["1p2z", "1m2z", "1z2p", "1z2m"] for name in sim_names: path = f"{grids_dir}/{name}_grid_{num}/{cat_name}" @@ -71,23 +71,15 @@ def load_catalogs(self, verbose=True): if verbose: print(f" {len(self.cats[name]['ra'])} objects") - def _match_to_ref(self, name): - """Return (idx_ref, idx_sim) matched indices between name and 1z2z.""" - ref = self.cats["1z2z"] - sim = self.cats[name] - idx_ref, idx_sim = match_catalogs_radec( - ref["ra"], ref["dec"], - sim["ra"], sim["dec"], - thresh_deg=self.thresh, - ) - return idx_ref, idx_sim - 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).""" e_key = f"e{comp + 1}" - idx_ref_p, idx_p = self._match_to_ref(name_p) - idx_ref_m, idx_m = self._match_to_ref(name_m) + 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}: {len(idx_p)} matched | {name_m}: {len(idx_m)} matched") From c67d94e1f36fc8eb68e622113dc729c0ac70e63c Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 4 Jul 2026 18:54:54 +0200 Subject: [PATCH 09/15] fix(image-sims): track params template + mask as im_init inputs im_init staged params.py by `cp {params.template}` where the template was a snakemake `params:` value, not a tracked `input:`. Editing the template (e.g. the ngmix-v2 grammar migration) therefore did NOT retrigger im_init, so run dirs kept stale old-grammar params.py and extract died on `no field of name NGMIX_Tpsf_NOSHEAR` -- only for sims whose run dir predated the edit. Move template + mask_src to `input:` so a template change re-stages into every run dir. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01D9atA9VPrg7VMC7oBArny9 --- workflow/rules/image_sims.smk | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/workflow/rules/image_sims.smk b/workflow/rules/image_sims.smk index 16c362db..7f464e19 100644 --- a/workflow/rules/image_sims.smk +++ b/workflow/rules/image_sims.smk @@ -160,12 +160,16 @@ rule im_init: 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: - template=PARAMS_TEMPLATE, - mask_src=lambda wc: os.path.join(SPV_REPO, MASK_CONFIG), config_dir=CONFIG_DIR, run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", cfis=lambda wc: f"{GRIDS_BASE}/{wc.sim}/cfis", @@ -176,8 +180,8 @@ rule im_init: # 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 {params.template} {output.params} && " - "ln -sf {params.mask_src} {output.mask} && " + "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" From 4e38326731df149f76878491971126eaf42fbae7 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 4 Jul 2026 19:24:54 +0200 Subject: [PATCH 10/15] fix(image-sims): paired "pool" m-bias estimator (shape-noise cancellation) The extracted ImageSimMBias._m_c_pair was UNPAIRED: it matched the +g and -g sims each to the reference separately, meaned them independently, and subtracted. Intrinsic shape (sigma_e ~ 0.3) never cancelled, so a single tile gave sigma(m) ~ 1.5 -- unconstrained. Restore the paired estimator the SKiLLS sims are built for: match the +g and -g sims directly to each other (same galaxies, opposite input shear; verified 98.8% co-located to <0.05" with e-correlation 0.81), difference per object, m = <(e_+ - e_-)/(2 g_in) - 1> , c = <(e_+ + e_-)/2> , and bootstrap the paired objects (one resample applied to both sims) so the per-object cancellation carries into the error. Intrinsic shape cancels in m; c is a sum so it stays shape-noise limited (as it must). This also restores Martin's direct +g<->-g matching, which the extraction (10baa52) regressed, and goes beyond it (his bootstrap still resampled the two means independently). New test test_mbias_pool_cancels_shape_noise builds catalogues with realistic intrinsic scatter and asserts sigma(m) sits far below the unpaired shape-noise floor -- the property the pooling exists to deliver. Existing exact-recovery test still passes to machine precision. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01HmxPiKJEechqueT3YU1yrx --- src/sp_validation/image_sims.py | 69 ++++++++++++---------- src/sp_validation/tests/test_image_sims.py | 69 +++++++++++++++++++++- 2 files changed, 104 insertions(+), 34 deletions(-) diff --git a/src/sp_validation/image_sims.py b/src/sp_validation/image_sims.py index 703ba128..9b79d07e 100644 --- a/src/sp_validation/image_sims.py +++ b/src/sp_validation/image_sims.py @@ -70,54 +70,61 @@ def load_catalogs(self, verbose=True): if verbose: print(f" {len(self.cats[name]['ra'])} objects") - def _match_to_ref(self, name): - """Return (idx_ref, idx_sim) matched indices between name and 1z2z.""" - ref = self.cats["1z2z"] - sim = self.cats[name] - idx_ref, idx_sim = match_catalogs_radec( - ref["ra"], - ref["dec"], - sim["ra"], - sim["dec"], - thresh_deg=self.thresh, - ) - return idx_ref, idx_sim - 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).""" + """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.) + """ e_key = f"e{comp + 1}" - idx_ref_p, idx_p = self._match_to_ref(name_p) - idx_ref_m, idx_m = self._match_to_ref(name_m) + # Match the +g and -g sims to each other: same galaxies, opposite shear. + 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}: {len(idx_p)} matched | {name_m}: {len(idx_m)} matched" - ) + print(f" {name_p} <-> {name_m}: {len(idx_p)} paired objects") e_p = self.cats[name_p][e_key][idx_p] w_p = self.cats[name_p]["w"][idx_p] e_m = self.cats[name_m][e_key][idx_m] w_m = self.cats[name_m]["w"][idx_m] - mean_ep = np.average(e_p, weights=w_p) - mean_em = np.average(e_m, weights=w_m) + # 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 = (mean_ep - mean_em) / (2 * self.g_in) - 1 - c = (mean_ep + mean_em) / 2 + m = np.average(d, weights=w) + c = np.average(s, weights=w) - # Bootstrap errors + # Paired bootstrap: resample objects once and apply the same draw to + # both sims, so the per-object cancellation in `d` is preserved in + # the error estimate. rng = np.random.default_rng(seed=42) + n = len(d) m_boot = np.empty(self.n_boot) c_boot = np.empty(self.n_boot) - n_p, n_m = len(e_p), len(e_m) for i in range(self.n_boot): - ib_p = rng.integers(0, n_p, n_p) - ib_m = rng.integers(0, n_m, n_m) - ep_b = np.average(e_p[ib_p], weights=w_p[ib_p]) - em_b = np.average(e_m[ib_m], weights=w_m[ib_m]) - m_boot[i] = (ep_b - em_b) / (2 * self.g_in) - 1 - c_boot[i] = (ep_b + em_b) / 2 + ib = rng.integers(0, n, n) + m_boot[i] = np.average(d[ib], weights=w[ib]) + c_boot[i] = np.average(s[ib], weights=w[ib]) return m, np.std(m_boot), c, np.std(c_boot) diff --git a/src/sp_validation/tests/test_image_sims.py b/src/sp_validation/tests/test_image_sims.py index 8988347c..d1bc6e9e 100644 --- a/src/sp_validation/tests/test_image_sims.py +++ b/src/sp_validation/tests/test_image_sims.py @@ -6,10 +6,13 @@ 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|``). Objects are matched to -the reference by RA/Dec, then, per component, +``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 = ( - ) / (2 |g|) - 1 , c = ( + ) / 2 . + 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 @@ -112,3 +115,63 @@ def test_mbias_recovers_injected_values(tmp_path): # 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 + + +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_col": "w_des", + "n_bootstrap": 200, + } + 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 From 7b7d15d42b6e2ae0a822a2c9d38409ee91ff5be9 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 4 Jul 2026 19:43:04 +0200 Subject: [PATCH 11/15] docs(image-sims): note reference-load + non-bijective matching caveats Fresh-eyes review notes on 4e38326: document that 1z2z is the unsheared reference (loaded for null-test diagnostics, not consumed by the +g/-g pool estimator) and that the RA/Dec match is nearest-neighbour, not a strict bijection -- effectively 1:1 on well-separated grid sims (~99% co-located to <0.05" on SKiLLS grid_1) but could dilute cancellation on denser fields. No behavioural change. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01HmxPiKJEechqueT3YU1yrx --- src/sp_validation/image_sims.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/sp_validation/image_sims.py b/src/sp_validation/image_sims.py index 9b79d07e..8ade3b91 100644 --- a/src/sp_validation/image_sims.py +++ b/src/sp_validation/image_sims.py @@ -60,6 +60,9 @@ def load_catalogs(self, verbose=True): grids_dir = self.cfg["grids_dir"] num = self.cfg["num"] cat_name = self.cfg["catalog_name"] + # 1z2z is the unsheared reference. The +g/-g pool estimator does not use + # it (it pairs the sheared sims directly); it is loaded for completeness + # and for null-test diagnostics on the zero-shear catalogue. sim_names = ["1z2z", "1p2z", "1m2z", "1z2p", "1z2m"] for name in sim_names: @@ -89,6 +92,11 @@ def _m_c_pair(self, name_p, name_m, comp, verbose=True): e_key = f"e{comp + 1}" # 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"], From e38ace9cfd3911ff4add90deb1bae123195a2189 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Tue, 7 Jul 2026 16:21:45 +0200 Subject: [PATCH 12/15] implementing additive bias correction flag, set to False for image simulations with constant shear --- config/calibration/mask_v1.X.4_im_sim.yaml | 4 +++ config/calibration/mask_v1.X.9_im_sim.yaml | 4 +++ notebooks/calibrate_comprehensive_cat.py | 9 ++++- scripts/compute_m_bias_image_sims.py | 12 +++++-- scripts/diagnostics_image_sims.py | 5 +-- src/sp_validation/calibration.py | 18 ++++++---- src/sp_validation/image_sims.py | 42 ++++++++++++++++------ 7 files changed, 71 insertions(+), 23 deletions(-) diff --git a/config/calibration/mask_v1.X.4_im_sim.yaml b/config/calibration/mask_v1.X.4_im_sim.yaml index f8045551..c1e1fd49 100644 --- a/config/calibration/mask_v1.X.4_im_sim.yaml +++ b/config/calibration/mask_v1.X.4_im_sim.yaml @@ -60,3 +60,7 @@ metacal: # Correct relative size for ellipticity? gal_size_corr_ell: False + + # Subtract additive bias (mean shear)? Keep False for constant-shear + # image sims, where the mean shear is the signal (issue #226) + 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 index dd8734fd..a88acd82 100644 --- a/config/calibration/mask_v1.X.9_im_sim.yaml +++ b/config/calibration/mask_v1.X.9_im_sim.yaml @@ -69,3 +69,7 @@ metacal: # Weight for global response matrix, None for unweighted mean global_R_weight: w + + # Subtract additive bias (mean shear)? Use False for constant-shear + # image sims + additive_correction: False diff --git a/notebooks/calibrate_comprehensive_cat.py b/notebooks/calibrate_comprehensive_cat.py index bb529c78..bd0cd184 100644 --- a/notebooks/calibrate_comprehensive_cat.py +++ b/notebooks/calibrate_comprehensive_cat.py @@ -106,8 +106,15 @@ ) # %% +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) + calibration.get_calibrated_m_c( + gal_metacal, + additive_correction=additive_correction, + ) ) num_ok = len(g_corr_mc[0]) diff --git a/scripts/compute_m_bias_image_sims.py b/scripts/compute_m_bias_image_sims.py index 18421cde..b699c10c 100644 --- a/scripts/compute_m_bias_image_sims.py +++ b/scripts/compute_m_bias_image_sims.py @@ -180,14 +180,14 @@ def main(): config = yaml.safe_load(f) print(f"Config: {args.config}") - print(f"Grids : {config['grids_dir']}") + print(f"Grids : {config['base']}") 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']) + n_tiles = get_n_tiles(config['base'], config['num']) if n_tiles: args.n_tiles = n_tiles print(f"Auto-detected {n_tiles} tiles") @@ -197,6 +197,9 @@ def main(): print("Loading catalogues...") mb.load_catalogs(verbose=args.verbose) + if args.verbose: + mb.print_mean_ellipticities() + results = mb.run(verbose=True) print() @@ -211,7 +214,10 @@ def main(): # Cumulative tracking if args.cumulative: - results_dir = config.get("diagnostics_dir", config.get("results_dir", "results")) + results_dir = config.get( + "diagnostics_dir", + config.get("results_dir", os.path.join(config["base"], "results")), + ) os.makedirs(results_dir, exist_ok=True) else: results_dir = None diff --git a/scripts/diagnostics_image_sims.py b/scripts/diagnostics_image_sims.py index 866220f6..4773a62f 100644 --- a/scripts/diagnostics_image_sims.py +++ b/scripts/diagnostics_image_sims.py @@ -90,8 +90,9 @@ def plot_response(cats, out_dir, nbins=50): 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, density=True, + 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) @@ -137,7 +138,7 @@ def main(): with open(args.config) as f: config = yaml.safe_load(f) - grids_dir = config["grids_dir"] + grids_dir = config["base"] num = config["num"] cat_name = config.get("catalog_name", "shape_catalog_cut_ngmix.fits") out_dir = config.get("diagnostics_dir", f"{grids_dir}/diagnostics") diff --git a/src/sp_validation/calibration.py b/src/sp_validation/calibration.py index 0272822a..f87b1722 100644 --- a/src/sp_validation/calibration.py +++ b/src/sp_validation/calibration.py @@ -67,7 +67,7 @@ def get_calibrated_quantities(gal_metacal, shape_method='ngmix'): return g_corr, g_uncorr, w, mask -def get_calibrated_m_c(gal_metacal, shape_method='ngmix'): +def get_calibrated_m_c(gal_metacal, shape_method='ngmix', additive_correction=True): """Get Calibrated C. Return catalogue quantities for objects calibrated for multiplicative and @@ -79,6 +79,11 @@ def get_calibrated_m_c(gal_metacal, shape_method='ngmix'): galaxy metacalibration catalogue shape_method : string, optional, default='ngmix' shape measurement method, one in 'ngmix', 'galsim' + 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 ------- @@ -114,11 +119,12 @@ def get_calibrated_m_c(gal_metacal, shape_method='ngmix'): 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/image_sims.py b/src/sp_validation/image_sims.py index 92a8cd7f..00738a26 100644 --- a/src/sp_validation/image_sims.py +++ b/src/sp_validation/image_sims.py @@ -39,11 +39,14 @@ class ImageSimMBias: ---------- config : dict Configuration dictionary with keys: - - grids_dir : str, path to the grids directory + - base : str, path to the grids base 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| (e.g. 0.02) - match_radius_deg : float, matching radius in degrees + - pair_match : bool, match objects between sheared catalogues + (default True); if False, use all objects of each catalogue - w_col : str, weight column name (default 'w_des') - n_bootstrap : int, number of bootstrap resamples for errors """ @@ -52,15 +55,16 @@ def __init__(self, config): self.cfg = config self.g_in = config["shear_amplitude"] self.thresh = config.get("match_radius_deg", 0.0002) + self.pair_match = config.get("pair_match", True) self.w_col = config.get("w_col", "w_des") self.n_boot = config.get("n_bootstrap", 500) self.cats = {} def load_catalogs(self, verbose=True): """Load the 4 sheared catalogues (1p2z, 1m2z, 1z2p, 1z2m).""" - grids_dir = self.cfg["grids_dir"] + grids_dir = self.cfg["base"] num = self.cfg["num"] - cat_name = self.cfg["catalog_name"] + cat_name = self.cfg.get("catalog_name", "shape_catalog_cut_ngmix.fits") sim_names = ["1p2z", "1m2z", "1z2p", "1z2m"] for name in sim_names: @@ -71,18 +75,34 @@ def load_catalogs(self, verbose=True): if verbose: print(f" {len(self.cats[name]['ra'])} objects") + def print_mean_ellipticities(self): + """Print weighted mean e1, e2 for each sheared catalogue, as check.""" + print("\nMean weighted ellipticities (all objects):") + for name, cat in self.cats.items(): + mean_e1 = np.average(cat["e1"], weights=cat["w"]) + mean_e2 = np.average(cat["e2"], weights=cat["w"]) + 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).""" e_key = f"e{comp + 1}" - 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}: {len(idx_p)} matched | {name_m}: {len(idx_m)} matched") + if self.pair_match: + 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}: {len(idx_p)} matched | {name_m}: {len(idx_m)} matched") + 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] w_p = self.cats[name_p]["w"][idx_p] From d0984f209a56c76a604706363fbcd500b2c5b2c4 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 07:57:26 +0200 Subject: [PATCH 13/15] feat(image-sims): w_col=None gives unit weights (no-weighting mode, #227) Fabian's #227 verdict excludes shape weights from sim m-bias calibration (w_des is pathological on sims, N_eff ~1-5%; w_iv behaves but is also excluded). w_col: null in the workflow config now flows through the .smk untouched and lands here as unit weights. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01UvQvyJTEwYcLCwDEgpWMJ6 --- src/sp_validation/image_sims.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/sp_validation/image_sims.py b/src/sp_validation/image_sims.py index 8ade3b91..be3bbde1 100644 --- a/src/sp_validation/image_sims.py +++ b/src/sp_validation/image_sims.py @@ -19,7 +19,11 @@ def _load_cat(path, e_col, w_col): - """Load RA, Dec, ellipticity component and weight from a FITS catalogue.""" + """Load RA, Dec, ellipticity component and weight from a FITS catalogue. + + ``w_col=None`` gives every object unit weight — the no-weighting mode for + m-bias runs (#227: shape weights are excluded from sim calibration). + """ with fits.open(path) as hdul: data = hdul[1].data return { @@ -27,7 +31,7 @@ def _load_cat(path, e_col, w_col): "dec": data["Dec"].copy(), "e1": data["e1"].copy(), "e2": data["e2"].copy(), - "w": data[w_col].copy(), + "w": data[w_col].copy() if w_col else np.ones(len(data["RA"])), } @@ -43,7 +47,8 @@ class ImageSimMBias: - catalog_name : str, filename of the cut catalogue - shear_amplitude : float, input shear |g| (e.g. 0.02) - match_radius_deg : float, matching radius in degrees - - w_col : str, weight column name (default 'w_des') + - w_col : str or None, weight column name (default 'w_des'); + None → unit weights (no weighting, per the #227 verdict) - n_bootstrap : int, number of bootstrap resamples for errors """ From 64ec1d7877f6d9a6aad2a42fa48f469bd3a75969 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 08:03:47 +0200 Subject: [PATCH 14/15] feat(calibration): mirror Martin's #226 additive_correction flag into the branch layout Martin implemented the #226 fix (skip additive-bias subtraction for constant-shear sims) on his image_sims_val branch (e38ace9c) against the old notebooks/ layout; #225 moved calibration to scripts/calibration/. This ports the identical semantics: get_calibrated_m_c gains additive_correction (default True, data path unchanged); the calibrate script reads it from the mask config's metacal: section; both im_sim mask configs set False. c/c_err still computed as diagnostics. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01UvQvyJTEwYcLCwDEgpWMJ6 --- config/calibration/mask_v1.X.4_im_sim.yaml | 4 ++++ config/calibration/mask_v1.X.9_im_sim.yaml | 4 ++++ .../calibration/calibrate_comprehensive_cat.py | 7 ++++++- src/sp_validation/calibration.py | 16 +++++++++++----- 4 files changed, 25 insertions(+), 6 deletions(-) diff --git a/config/calibration/mask_v1.X.4_im_sim.yaml b/config/calibration/mask_v1.X.4_im_sim.yaml index f8045551..c3e23491 100644 --- a/config/calibration/mask_v1.X.4_im_sim.yaml +++ b/config/calibration/mask_v1.X.4_im_sim.yaml @@ -60,3 +60,7 @@ metacal: # Correct relative size for ellipticity? gal_size_corr_ell: False + + # 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 index 4e0377f4..aeb36d40 100644 --- a/config/calibration/mask_v1.X.9_im_sim.yaml +++ b/config/calibration/mask_v1.X.9_im_sim.yaml @@ -69,3 +69,7 @@ metacal: # Weight for global response matrix, None for unweighted mean global_R_weight: w + + # 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/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 From 04464e172896e4b6721c6edff90ae611d930c3c3 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 15:08:08 +0200 Subject: [PATCH 15/15] config(image-sims): unweighted global response for sim calibration (#227) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The w_des-weighted global R was measured on N_eff ~ 20-100 objects and differs from the unweighted mean by up to 17% (0.780 vs 0.669 on tile 233.293) — no weights anywhere in sim m-bias, per Fabian's #227 verdict. calibrate_comprehensive_cat reads the key fail-fast, so both im_sim mask configs now carry it explicitly. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01UvQvyJTEwYcLCwDEgpWMJ6 --- config/calibration/mask_v1.X.4_im_sim.yaml | 4 ++++ config/calibration/mask_v1.X.9_im_sim.yaml | 6 ++++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/config/calibration/mask_v1.X.4_im_sim.yaml b/config/calibration/mask_v1.X.4_im_sim.yaml index c3e23491..794566d4 100644 --- a/config/calibration/mask_v1.X.4_im_sim.yaml +++ b/config/calibration/mask_v1.X.4_im_sim.yaml @@ -61,6 +61,10 @@ metacal: # 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.yaml b/config/calibration/mask_v1.X.9_im_sim.yaml index aeb36d40..1c6ff672 100644 --- a/config/calibration/mask_v1.X.9_im_sim.yaml +++ b/config/calibration/mask_v1.X.9_im_sim.yaml @@ -67,8 +67,10 @@ metacal: # Correct relative size for ellipticity? gal_size_corr_ell: False - # Weight for global response matrix, None for unweighted mean - global_R_weight: w + # 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