diff --git a/Dockerfile b/Dockerfile index c791877c..ff0ea389 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,5 +1,5 @@ # Development image with more bells and whistles -FROM ghcr.io/cosmostat/shapepipe:develop +FROM ghcr.io/cosmostat/shapepipe:im_sims RUN apt-get update -y --quiet --fix-missing && \ apt-get dist-upgrade -y --quiet --fix-missing && \ diff --git a/config/calibration/mask_v1.X.4_im_sim.yaml b/config/calibration/mask_v1.X.4_im_sim.yaml new file mode 100644 index 00000000..794566d4 --- /dev/null +++ b/config/calibration/mask_v1.X.4_im_sim.yaml @@ -0,0 +1,70 @@ +# Config file for masking and calibration. +# Standard cuts without coverage mask, type v1.X.4. + +# General parameters (can also given on command line) +params: + input_path: shape_catalog_comprehensive_ngmix.hdf5 + cmatrices: False + sky_regions: False + verbose: True + +# Masks +## Using columns in 'dat' group (ShapePipe flags) +dat: + # SExtractor flags + - col_name: FLAGS + label: SE FLAGS + kind: equal + value: 0 + + # Number of epochs + - col_name: N_EPOCH + label: r"$n_{\rm epoch}$" + kind: greater_equal + value: 2 + + # Magnitude range + - col_name: mag + label: mag range + kind: range + value: [15, 30] + + # ngmix flags + - col_name: NGMIX_MOM_FAIL + label: "ngmix moments failure" + kind: equal + value: 0 + + # invalid PSF ellipticities + - col_name: NGMIX_ELL_PSFo_NOSHEAR_0 + label: "bad PSF ellipticity comp 1" + kind: not_equal + value: -10 + - col_name: NGMIX_ELL_PSFo_NOSHEAR_1 + label: "bad PSF ellipticity comp 2" + kind: not_equal + value: -10 + +# Metacal parameters +metacal: + # Ellipticity dispersion + sigma_eps_prior: 0.34 + + # Signal-to-noise range + gal_snr_min: 10 + gal_snr_max: 500 + + # Relative-size (hlr / hlr_psf) range + gal_rel_size_min: 0.5 + gal_rel_size_max: 3 + + # Correct relative size for ellipticity? + gal_size_corr_ell: False + + # Weight for global response matrix, None for unweighted mean. + # Unweighted for image sims: no weights anywhere in sim m-bias (#227). + global_R_weight: null + + # Subtract additive bias (mean shear)? Use False for constant-shear + # image sims + additive_correction: False diff --git a/config/calibration/mask_v1.X.9_im_sim.yaml b/config/calibration/mask_v1.X.9_im_sim.yaml new file mode 100644 index 00000000..1c6ff672 --- /dev/null +++ b/config/calibration/mask_v1.X.9_im_sim.yaml @@ -0,0 +1,77 @@ +# Config file for masking and calibration. +# Less conservative cuts, type v1.X.9 (e.g. for matching with spectroscopic sample). + +# General parameters (can also given on command line) +params: + input_path: shape_catalog_comprehensive_ngmix.fits + cmatrices: False + sky_regions: False + verbose: True + +# Masks +## Using columns in 'dat' group (ShapePipe flags) +dat: + # SExtractor flags + - col_name: FLAGS + label: SE FLAGS + kind: smaller_equal + value: 2 + + # Duplicate objects + - col_name: overlap + label: tile overlap + kind: equal + value: True + + # Number of epochs + - col_name: N_EPOCH + label: r"$n_{\rm epoch}$" + kind: greater_equal + value: 1 + + # Magnitude range + - col_name: mag + label: mag range + kind: range + value: [15, 30] + + # ngmix flags + - col_name: NGMIX_MCAL_TYPES_FAIL + label: "ngmix moments failure" + kind: equal + value: 0 + + # invalid PSF ellipticities (ShapePipe-v2 grammar: scalar G1/G2 components) + - col_name: NGMIX_G1_PSF_ORIG_NOSHEAR + label: "bad PSF ellipticity comp 1" + kind: not_equal + value: -10 + - col_name: NGMIX_G2_PSF_ORIG_NOSHEAR + label: "bad PSF ellipticity comp 2" + kind: not_equal + value: -10 + +# Metacal parameters +metacal: + # Ellipticity dispersion + sigma_eps_prior: 0.34 + + # Signal-to-noise range + gal_snr_min: 5 + gal_snr_max: 500 + + # Relative-size (hlr / hlr_psf) range + gal_rel_size_min: 0.25 + gal_rel_size_max: 10 + + # Correct relative size for ellipticity? + gal_size_corr_ell: False + + # Weight for global response matrix, None for unweighted mean. + # Unweighted for image sims: no weights anywhere in sim m-bias (#227), + # and the w_des-weighted R had N_eff ~ 20-100 objects. + global_R_weight: null + + # Subtract additive bias (mean shear)? Use False for constant-shear + # image sims + additive_correction: False diff --git a/scripts/calibration/calibrate_comprehensive_cat.py b/scripts/calibration/calibrate_comprehensive_cat.py index 4b15ed89..4928df9b 100644 --- a/scripts/calibration/calibrate_comprehensive_cat.py +++ b/scripts/calibration/calibrate_comprehensive_cat.py @@ -96,8 +96,13 @@ ) # %% +additive_correction = cm.get("additive_correction", True) +if not additive_correction: + print("Additive bias correction disabled (additive_correction: False)") + g_corr_mc, g_uncorr, w, mask_metacal, c, c_err = calibration.get_calibrated_m_c( - gal_metacal + gal_metacal, + additive_correction=additive_correction, ) num_ok = len(g_corr_mc[0]) diff --git a/scripts/calibration/extract_info.py b/scripts/calibration/extract_info.py index 1cbf836f..123c7437 100644 --- a/scripts/calibration/extract_info.py +++ b/scripts/calibration/extract_info.py @@ -32,6 +32,7 @@ import os import sys +import h5py import numpy as np from astropy.io import fits @@ -51,8 +52,11 @@ # ### Create and open output files and directories -make_out_dirs(output_dir, plot_dir, [], verbose=verbose) -stats_file = open_stats_file(plot_dir, stats_file_name) +os.makedirs(output_dir, exist_ok=True) +stats_file = open_stats_file(output_dir, stats_file_name) + +output_shape_cat_stem = output_shape_cat_base +output_ext = output_format # ## 2. Load data # @@ -107,6 +111,8 @@ tile_IDs, path_tile_ID, path_found_ID, path_missing_ID, verbose=verbose ) +print_stats(f"Tiles in input catalogue: {n_found}", stats_file, verbose=verbose) + # ### Load star catalogue if star_cat_path: @@ -146,42 +152,40 @@ verbose=verbose, ) -# #### Refine: Match to valid, unflagged galaxy sample - -# + -# Flags to indicate valid star sample - -m_star = ( - (dd["FLAGS"][ind_star] == 0) - & (dd["IMAFLAGS_ISO"][ind_star] == 0) - & (dd["NGMIX_MCAL_FLAGS"][ind_star] == 0) - & (dd["NGMIX_G1_PSF_ORIG_NOSHEAR"][ind_star] != -10) -) - -ra_star, dec_star, g_star_psf = spv_cat.match_subsample( - dd, - ind_star, - m_star, - [col_name_ra, col_name_dec], - key_PSF_g1, - key_PSF_g2, - n_star_tot, - stats_file, - verbose=verbose, -) -# - + # Flags to indicate valid star sample. + # Star matching, PSF-catalogue output and star metacalibration are all + # diagnostics that require an input star catalogue; the image-simulation + # pipeline has none (star_cat_path is None), so every star-dependent block + # below is guarded and simply skipped for the sims. + + m_star = ( + (dd["FLAGS"][ind_star] == 0) + & (dd["IMAFLAGS_ISO"][ind_star] == 0) + & (dd["NGMIX_MCAL_FLAGS"][ind_star] == 0) + & (dd["NGMIX_G1_PSF_ORIG_NOSHEAR"][ind_star] != -10) + ) -# MKDEBUG: Moved from end of this script + ra_star, dec_star, g_star_psf = spv_cat.match_subsample( + dd, + ind_star, + m_star, + [col_name_ra, col_name_dec], + key_PSF_g1, + key_PSF_g2, + n_star_tot, + stats_file, + verbose=verbose, + ) -# ### Write PSF catalogue with multi-epoch shapes from shape measurement methods + # ### Write PSF catalogue with multi-epoch shapes from shape measurement methods -spv_cat.write_PSF_cat( - f"{output_PSF_cat_base}_{shape}.fits", - ra_star, - dec_star, - g_star_psf[0], - g_star_psf[1], -) + spv_cat.write_PSF_cat( + f"{output_PSF_cat_base}_{shape}.fits", + ra_star, + dec_star, + g_star_psf[0], + g_star_psf[1], + ) # ## Check for objects with invalid PSF @@ -253,8 +257,9 @@ if verbose: print("Writing comprehensive catalogue...") +comprehensive_cat_path = f"{output_shape_cat_stem}_comprehensive_{shape}{output_ext}" spv_cat.write_shape_catalog( - f"{output_shape_cat_base}_comprehensive_{shape}.fits", + comprehensive_cat_path, ra_all, dec_all, iv_w, @@ -265,6 +270,17 @@ add_cols=ext_cols_pre_cal, add_cols_format=add_cols_pre_cal_format, ) + +# Write tile count to HDF5 attributes +if output_ext == ".hdf5": + try: + with h5py.File(comprehensive_cat_path, "a") as hf: + hf.attrs["n_tiles"] = n_found + if verbose: + print(f" Added n_tiles={n_found} to HDF5 attributes") + except Exception as e: + if verbose: + print(f" Warning: could not add n_tiles attribute: {e}") # - do_selection_calibration = False @@ -275,6 +291,7 @@ else: if verbose: print("Continuing with selection and calibration") + os.makedirs(os.path.join(output_dir, plot_dir), exist_ok=True) # ## 4. Select galaxies @@ -628,23 +645,24 @@ # ## Metacalibration for stars -star_metacal = metacal(dd[ind_star], m_star, masking_type="star", verbose=verbose) +if star_cat_path: + star_metacal = metacal(dd[ind_star], m_star, masking_type="star", verbose=verbose) -# #### Number density + # #### Number density -# + -# mask for 'no shear' images + # + + # mask for 'no shear' images -mask_ns_stars = star_metacal.mask_dict["ns"] -n_star = len(star_metacal.ns["g1"][mask_ns_stars]) + mask_ns_stars = star_metacal.mask_dict["ns"] + n_star = len(star_metacal.ns["g1"][mask_ns_stars]) -print_stats(f"Number of stars = {n_star}", stats_file, verbose=verbose) -print_stats( - "Star density = {:.2f} stars/deg2".format(n_star / area_deg2), - stats_file, - verbose=verbose, -) -# - + print_stats(f"Number of stars = {n_star}", stats_file, verbose=verbose) + print_stats( + "Star density = {:.2f} stars/deg2".format(n_star / area_deg2), + stats_file, + verbose=verbose, + ) + # - # ## Additive bias # Use raw, uncorrected ellipticities. @@ -730,20 +748,21 @@ print_stats(rs, stats_file, verbose=verbose) # + -print_stats("stars:", stats_file, verbose=verbose) +if star_cat_path: + print_stats("stars:", stats_file, verbose=verbose) -print_stats("total response matrix:", stats_file, verbose=verbose) -rs = np.array2string(star_metacal.R) -print_stats(rs, stats_file, verbose=verbose) + print_stats("total response matrix:", stats_file, verbose=verbose) + rs = np.array2string(star_metacal.R) + print_stats(rs, stats_file, verbose=verbose) -print_stats("shear response matrix:", stats_file, verbose=verbose) -R_shear_stars = np.mean(star_metacal.R_shear, 2) -rs = np.array2string(R_shear_stars) -print_stats(rs, stats_file, verbose=verbose) + print_stats("shear response matrix:", stats_file, verbose=verbose) + R_shear_stars = np.mean(star_metacal.R_shear, 2) + rs = np.array2string(R_shear_stars) + print_stats(rs, stats_file, verbose=verbose) -print_stats("selection response matrix:", stats_file, verbose=verbose) -rs = np.array2string(star_metacal.R_selection) -print_stats(rs, stats_file, verbose=verbose) + print_stats("selection response matrix:", stats_file, verbose=verbose) + rs = np.array2string(star_metacal.R_selection) + print_stats(rs, stats_file, verbose=verbose) # - # ### Plot distribution of response matrix elements @@ -757,14 +776,11 @@ linestyles = ["-", "-", ":", ":"] # + -labels = ["$R_{11}$ galaxies", "$R_{22}$ galaxies", "$R_{11}$ stars", "$R_{22}$ stars"] - -xs = [ - gal_metacal.R_shear[0, 0], - gal_metacal.R_shear[1, 1], - star_metacal.R_shear[0, 0], - star_metacal.R_shear[1, 1], -] +labels = ["$R_{11}$ galaxies", "$R_{22}$ galaxies"] +xs = [gal_metacal.R_shear[0, 0], gal_metacal.R_shear[1, 1]] +if star_cat_path: + labels += ["$R_{11}$ stars", "$R_{22}$ stars"] + xs += [star_metacal.R_shear[0, 0], star_metacal.R_shear[1, 1]] title = shape out_name = f"R_{shape}_diag.pdf" @@ -779,19 +795,16 @@ x_range, n_bin, out_path, - colors=colors, - linestyles=linestyles, + colors=colors[: len(xs)], + linestyles=linestyles[: len(xs)], ) # + -labels = ["$R_{12}$ galaxies", "$R_{21}$ galaxies", "$R_{12}$ stars", "$R_{21}$ stars"] - -xs = [ - gal_metacal.R_shear[0, 1], - gal_metacal.R_shear[1, 0], - star_metacal.R_shear[0, 1], - star_metacal.R_shear[1, 0], -] +labels = ["$R_{12}$ galaxies", "$R_{21}$ galaxies"] +xs = [gal_metacal.R_shear[0, 1], gal_metacal.R_shear[1, 0]] +if star_cat_path: + labels += ["$R_{12}$ stars", "$R_{21}$ stars"] + xs += [star_metacal.R_shear[0, 1], star_metacal.R_shear[1, 0]] title = shape out_name = f"R_{shape}_offdiag.pdf" out_path = os.path.join(plot_dir, out_name) @@ -805,8 +818,8 @@ x_range, n_bin, out_path, - colors=colors, - linestyles=linestyles, + colors=colors[: len(xs)], + linestyles=linestyles[: len(xs)], ) # - @@ -845,49 +858,50 @@ ) # + -xs = [star_metacal.ns["g1"][mask_ns_stars], star_metacal.ns["g2"][mask_ns_stars]] -weights = [star_metacal.ns["w"][mask_ns_stars]] * 2 - -title = "stars" -out_name = f"ell_stars_{shape}.pdf" -out_path = os.path.join(plot_dir, out_name) - -plot_histograms( - xs, - labels, - title, - x_label, - y_label, - x_range, - n_bin, - out_path, - weights=weights, - colors=colors, - linestyles=linestyles, -) -# - - -x_range = (-0.15, 0.15) -n_bin = 250 - -# + -xs = [dd[key_PSF_g1][mask_ns_stars], dd[key_PSF_g2][mask_ns_stars]] -title = "PSF" -out_name = f"ell_PSF_{shape}.pdf" -out_path = os.path.join(plot_dir, out_name) - -plot_histograms( - xs, - labels, - title, - x_label, - y_label, - x_range, - n_bin, - out_path, - colors=colors, - linestyles=linestyles, -) +if star_cat_path: + xs = [star_metacal.ns["g1"][mask_ns_stars], star_metacal.ns["g2"][mask_ns_stars]] + weights = [star_metacal.ns["w"][mask_ns_stars]] * 2 + + title = "stars" + out_name = f"ell_stars_{shape}.pdf" + out_path = os.path.join(plot_dir, out_name) + + plot_histograms( + xs, + labels, + title, + x_label, + y_label, + x_range, + n_bin, + out_path, + weights=weights, + colors=colors, + linestyles=linestyles, + ) + # - + + x_range = (-0.15, 0.15) + n_bin = 250 + + # + + xs = [dd[key_PSF_g1][mask_ns_stars], dd[key_PSF_g2][mask_ns_stars]] + title = "PSF" + out_name = f"ell_PSF_{shape}.pdf" + out_path = os.path.join(plot_dir, out_name) + + plot_histograms( + xs, + labels, + title, + x_label, + y_label, + x_range, + n_bin, + out_path, + colors=colors, + linestyles=linestyles, + ) # - # ## Magnitudes @@ -951,7 +965,7 @@ # ### Write basic shape catalogue spv_cat.write_shape_catalog( - f"{output_shape_cat_base}_{shape}.fits", + f"{output_shape_cat_stem}_{shape}{output_ext}", ra, dec, w, @@ -990,7 +1004,7 @@ # Extended catalogue with SNR, individual R matrices, ext_cols spv_cat.write_shape_catalog( - f"{output_shape_cat_base}_extended_{shape}.fits", + f"{output_shape_cat_stem}_extended_{shape}{output_ext}", ra, dec, w, @@ -1020,4 +1034,4 @@ ra = dd["RA"][cut_overlap] dec = dd["DEC"][cut_overlap] tile_id = dd["TILE_ID"][cut_overlap] - write_galaxy_cat(f"{output_shape_cat_base}.fits", ra, dec, tile_id) + write_galaxy_cat(f"{output_shape_cat_stem}{output_ext}", ra, dec, tile_id) diff --git a/scripts/calibration/params.py b/scripts/calibration/params.py index 2d0bd45f..5b2ef0c9 100644 --- a/scripts/calibration/params.py +++ b/scripts/calibration/params.py @@ -105,6 +105,9 @@ ## Output +### Output file format extension: '.fits' or '.hdf5' +output_format = ".hdf5" + ### Additional output columns add_cols = [ "FLUX_RADIUS", diff --git a/scripts/compute_m_bias_image_sims.py b/scripts/compute_m_bias_image_sims.py new file mode 100644 index 00000000..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/scripts/diagnostics_image_sims.py b/scripts/diagnostics_image_sims.py new file mode 100644 index 00000000..4773a62f --- /dev/null +++ b/scripts/diagnostics_image_sims.py @@ -0,0 +1,176 @@ +#!/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, range=(-1, 2), density=True, + histtype="step", label=name, color=COLORS[name]) + ax.set_xlim(-1, 2) + ax.set_xlabel(col) + ax.set_ylabel("normalised count") + ax.legend(fontsize=7) + fig.suptitle("Response matrix elements") + fig.tight_layout() + return savefig(fig, out_dir, "response_hist") + + +def plot_psf_leakage(cats, 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["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") + + 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/calibration.py b/src/sp_validation/calibration.py index 9032bc63..e37fc4d0 100644 --- a/src/sp_validation/calibration.py +++ b/src/sp_validation/calibration.py @@ -56,7 +56,7 @@ def get_calibrated_quantities(gal_metacal): return g_corr, g_uncorr, w, mask -def get_calibrated_m_c(gal_metacal): +def get_calibrated_m_c(gal_metacal, additive_correction=True): """Get Calibrated C. Return catalogue quantities for objects calibrated for multiplicative and @@ -66,6 +66,11 @@ def get_calibrated_m_c(gal_metacal): ---------- gal_metacal : dict galaxy metacalibration catalogue + additive_correction : bool, optional, default=True + if False, do not subtract the additive bias c from the shear + estimates; use for constant-shear image sims, where the mean + shear is the signal (see issue #226). c and c_err are still + computed and returned Returns ------- @@ -99,10 +104,11 @@ def get_calibrated_m_c(gal_metacal): c_err[comp] = np.std(g_uncorr[comp]) # Shear estimate corrected for additive bias - g_corr_mc = np.zeros_like(g_corr) - c_corr = np.linalg.inv(gal_metacal.R).dot(c) - for comp in (0, 1): - g_corr_mc[comp] = g_corr[comp] - c_corr[comp] + g_corr_mc = np.copy(g_corr) + if additive_correction: + c_corr = np.linalg.inv(gal_metacal.R).dot(c) + for comp in (0, 1): + g_corr_mc[comp] = g_corr[comp] - c_corr[comp] return g_corr_mc, g_uncorr, w, mask_metacal, c, c_err diff --git a/src/sp_validation/catalog.py b/src/sp_validation/catalog.py index b1b175d9..e39d2076 100644 --- a/src/sp_validation/catalog.py +++ b/src/sp_validation/catalog.py @@ -12,6 +12,7 @@ """ import getpass +import os import h5py import numpy as np @@ -308,6 +309,35 @@ def match_subsample( return ra, dec, g +def match_catalogs_radec(ra1, dec1, ra2, dec2, thresh_deg=0.0002): + """Match two catalogues by RA/Dec. + + Match each object in catalogue 2 to the nearest in catalogue 1 + within a threshold. + + Parameters + ---------- + ra1, dec1 : array_like + coordinates of reference catalogue [deg] + ra2, dec2 : array_like + coordinates of catalogue to match [deg] + thresh_deg : float, optional + maximum separation [deg], default 0.0002 + + Returns + ------- + idx1 : ndarray of int + indices into catalogue 1 of matched objects + idx2 : ndarray of int + indices into catalogue 2 of matched objects + """ + coord1 = coords.SkyCoord(ra=ra1 * u.degree, dec=dec1 * u.degree) + coord2 = coords.SkyCoord(ra=ra2 * u.degree, dec=dec2 * u.degree) + idx1, sep, _ = coord2.match_to_catalog_sky(coord1) + mask = sep.deg < thresh_deg + return idx1[mask], np.where(mask)[0] + + def match_stars2(ra_gal, dec_gal, ra_star, dec_star, thresh=0.0002): """Add docstring. @@ -542,55 +572,91 @@ def write_shape_catalog( ) ) - # Write columns to FITS file - cols = [] - for col, _ in col_info_arr: - cols.append(col) - table_hdu = fits.BinTableHDU.from_columns(cols) - - # Add human-readable descriptions - for idx, col_info in enumerate(col_info_arr): - table_hdu.header[f"TTYPE{idx + 1}"] = ( - col_info[0].name, - col_info[1], - ) + ext = os.path.splitext(output_path)[1].lower() - # Primary HDU with information in header - primary_header = fits.Header() + if ext in (".hdf5", ".hdf", ".h5"): + # Build flat list of (name, 1d-array) pairs, splitting 2D columns + fields = [] + for col, _ in col_info_arr: + arr = np.asarray(col.array) + if arr.ndim == 2: + for idx in range(arr.shape[1]): + fields.append((f"{col.name}_{idx}", arr[:, idx])) + else: + fields.append((col.name, arr)) + + # Build structured numpy array and write as single "data" dataset + dtype = np.dtype([(name, arr.dtype) for name, arr in fields]) + structured = np.empty(len(fields[0][1]), dtype=dtype) + for name, arr in fields: + structured[name] = arr + + with h5py.File(output_path, "w") as f: + f.create_dataset("data", data=structured) + if add_header: + for key, val in add_header.items(): + f.attrs[key] = str(val) + if all(v is not None for v in (R, R_shear, R_select, c)): + f.attrs["R"] = R + f.attrs["R_shear"] = R_shear + f.attrs["R_select"] = R_select + f.attrs["c"] = c + if c_err is not None: + f.attrs["c1_err"] = c_err[0] + f.attrs["c2_err"] = c_err[1] + if sigma_epsilon is not None: + f.attrs["sig_eps"] = sigma_epsilon + if alpha_leakage is not None: + f.attrs["alpha"] = alpha_leakage - if add_header: - primary_header.update(add_header) + else: + # Write columns to FITS file + cols = [col for col, _ in col_info_arr] + table_hdu = fits.BinTableHDU.from_columns(cols) + + # Add human-readable descriptions + for idx, col_info in enumerate(col_info_arr): + table_hdu.header[f"TTYPE{idx + 1}"] = ( + col_info[0].name, + col_info[1], + ) - primary_header = cat.write_header_info_sp( - primary_header, - software_name="sp_validation", - software_version=__version__, - author=getpass.getuser(), - ) + # Primary HDU with information in header + primary_header = fits.Header() - if all(v is not None for v in (R, R_shear, R_select, c)): - cat.add_shear_bias_to_header(primary_header, R, R_shear, R_select, c) - if c_err is not None: - primary_header["c1_err"] = (c_err[0], "Standard deviation of c_1") - primary_header["c2_err"] = (c_err[1], "Standard deviation of c_2") + if add_header: + primary_header.update(add_header) - primary_header["w"] = "DES weight" + primary_header = cat.write_header_info_sp( + primary_header, + software_name="sp_validation", + software_version=__version__, + author=getpass.getuser(), + ) - if sigma_epsilon is not None: - primary_header["sig_eps"] = (sigma_epsilon, "Shape noise RMS") + if all(v is not None for v in (R, R_shear, R_select, c)): + cat.add_shear_bias_to_header(primary_header, R, R_shear, R_select, c) + if c_err is not None: + primary_header["c1_err"] = (c_err[0], "Standard deviation of c_1") + primary_header["c2_err"] = (c_err[1], "Standard deviation of c_2") - if alpha_leakage: - primary_header["alpha"] = ( - alpha_leakage, - "Mean scale-dependent PSF leakage", - ) + primary_header["w"] = "DES weight" + + if sigma_epsilon is not None: + primary_header["sig_eps"] = (sigma_epsilon, "Shape noise RMS") + + if alpha_leakage: + primary_header["alpha"] = ( + alpha_leakage, + "Mean scale-dependent PSF leakage", + ) - primary_hdu = fits.PrimaryHDU(header=primary_header) + primary_hdu = fits.PrimaryHDU(header=primary_header) - # Final file - hdu_list = fits.HDUList([primary_hdu, table_hdu]) + # Final file + hdu_list = fits.HDUList([primary_hdu, table_hdu]) - hdu_list.writeto(output_path, overwrite=True) + hdu_list.writeto(output_path, overwrite=True) def write_galaxy_cat(output_path, ra, dec, tile_id): diff --git a/src/sp_validation/catalog_builders.py b/src/sp_validation/catalog_builders.py index 1375ef8d..78dc7838 100644 --- a/src/sp_validation/catalog_builders.py +++ b/src/sp_validation/catalog_builders.py @@ -1124,6 +1124,22 @@ def read_cat(self, load_into_memory=False): fpath = self._params["input_path"] verbose = self._params["verbose"] + # Image-simulation path: a single per-run comprehensive catalogue in + # FITS, not the joined multi-patch HDF5 the data path builds. Read the + # FITS table directly into memory; there is no separate data_ext group. + extension = os.path.splitext(fpath)[1] + if extension == ".fits": + if verbose: + print(f"Reading FITS file {fpath}, HDU 1...") + dat = fits.getdata(fpath, 1) + dat_ext = None + if verbose: + print( + f"Found {len(dat)} (~{format.millify(len(dat))}) objects" + + " in catalogue" + ) + return dat, dat_ext + if verbose: print(f"Reading HDF5 file {fpath}...") diff --git a/src/sp_validation/image_sims.py b/src/sp_validation/image_sims.py new file mode 100644 index 00000000..1980d308 --- /dev/null +++ b/src/sp_validation/image_sims.py @@ -0,0 +1,216 @@ +"""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. + + ``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 { + "ra": data["RA"].copy(), + "dec": data["Dec"].copy(), + "e1": data["e1"].copy(), + "e2": data["e2"].copy(), + "w": data[w_col].copy() if w_col else np.ones(len(data["RA"])), + } + + +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 + (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 the +g and -g sheared + catalogues (default True); if False, use all objects of each + catalogue (the paired per-object cancellation is then unavailable) + - 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 + """ + + 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 5 sheared and reference catalogues.""" + grids_dir = self.cfg["grids_dir"] + num = self.cfg["num"] + cat_name = self.cfg.get("catalog_name", "shape_catalog_cut_ngmix.fits") + # 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: + 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 print_mean_ellipticities(self): + """Print weighted mean e1, e2 for each catalogue, as a check. + + Works with unit weights too (``w_col=None`` → w all ones), in which + case these are the plain unweighted means. + """ + 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). + + Paired ("pool") estimator. The +g and -g simulations inject opposite + input shear on the *same* galaxies, so matching them directly by + RA/Dec yields a one-to-one correspondence. Differencing the two + ellipticities per object, + + m = <(e_+ - e_-) / (2 g_in) - 1> , c = <(e_+ + e_-) / 2> , + + cancels the intrinsic shape (sigma_e ~ 0.3) object-by-object in the + multiplicative term, leaving only measurement noise -- so sigma(m) + shrinks by ~sigma_e/sigma_meas relative to differencing two + independent means. (The additive term c is a *sum*, so intrinsic + shape does not cancel there and its error stays shape-noise limited.) + + With ``pair_match=False`` the +g and -g sims are *not* matched: every + object of each catalogue is used, so the per-object cancellation is + lost and m, c fall back to differencing/summing the two independent + weighted means. The paired bootstrap likewise cannot be applied (the + two arrays generally have different lengths), so each side is resampled + independently per replicate. + """ + e_key = f"e{comp + 1}" + + if self.pair_match: + # Match the +g and -g sims to each other: same galaxies, opposite + # shear. This is a nearest-neighbour match within `thresh`, not a + # strict bijection -- on grid sims galaxies are well separated so + # pairs are effectively 1:1 (verified ~99% co-located to <0.05" on + # SKiLLS grid_1); on denser fields a small fraction could share a + # +g partner and dilute the cancellation. + idx_p, idx_m = match_catalogs_radec( + self.cats[name_p]["ra"], + self.cats[name_p]["dec"], + self.cats[name_m]["ra"], + self.cats[name_m]["dec"], + thresh_deg=self.thresh, + ) + if verbose: + print(f" {name_p} <-> {name_m}: {len(idx_p)} paired objects") + else: + idx_p = slice(None) + idx_m = slice(None) + if verbose: + print( + f" no pair-matching: {name_p}: {len(self.cats[name_p][e_key])}" + f" | {name_m}: {len(self.cats[name_m][e_key])} objects" + ) + + e_p = self.cats[name_p][e_key][idx_p] + 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] + + rng = np.random.default_rng(seed=42) + m_boot = np.empty(self.n_boot) + c_boot = np.empty(self.n_boot) + + if self.pair_match: + # Per-object shear-differenced (-> m) and summed (-> c) + # ellipticity, with a symmetric per-pair weight. + w = 0.5 * (w_p + w_m) + d = (e_p - e_m) / (2 * self.g_in) - 1 + s = (e_p + e_m) / 2 + + m = np.average(d, weights=w) + c = np.average(s, weights=w) + + # Paired bootstrap: resample objects once and apply the same draw + # to both sims, so the per-object cancellation in `d` is preserved + # in the error estimate. + n = len(d) + for i in range(self.n_boot): + 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]) + else: + # No matching: difference/sum the two independent weighted means. + mean_ep = np.average(e_p, weights=w_p) + mean_em = np.average(e_m, weights=w_m) + + m = (mean_ep - mean_em) / (2 * self.g_in) - 1 + c = (mean_ep + mean_em) / 2 + + # Unpaired bootstrap: the +g and -g arrays generally differ in + # length, so resample each side independently per replicate. + 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..d1bc6e9e --- /dev/null +++ b/src/sp_validation/tests/test_image_sims.py @@ -0,0 +1,177 @@ +"""UNIT TESTS FOR THE IMAGE-SIMULATION m/c ESTIMATOR. + +Exercise ``sp_validation.image_sims.ImageSimMBias`` -- the multiplicative and +additive shear-bias estimator used by the image-simulation workflow -- and the +``sp_validation.catalog.match_catalogs_radec`` helper it relies on. + +The estimator recovers ``m`` and ``c`` from five calibrated catalogues named +``1z2z`` (reference, no input shear), ``1p2z``/``1m2z`` (input shear +``g1 = +-|g|``) and ``1z2p``/``1z2m`` (``g2 = +-|g|``). The +g and -g sims are +matched to *each other* by RA/Dec -- same galaxies, opposite input shear -- and +the bias is the object-paired ("pool") average + + m = <(e_+ - e_-) / (2 |g|) - 1> , c = <(e_+ + e_-) / 2> , + +so the intrinsic shape cancels object-by-object in ``m``. + +We build synthetic catalogues in which the measured ellipticity is exactly +``e = (1 + m_true) g_in + c_true`` at shared positions, so the recovered m/c +must equal the injected values to machine precision -- an analytic check of +the estimator maths that needs no pipeline run. + +:Author: cdaley + +""" + +import numpy as np +import numpy.testing as npt +from astropy.io import fits + +from sp_validation.catalog import match_catalogs_radec +from sp_validation.image_sims import ImageSimMBias + +# Injected truth, shared across the synthetic-recovery test. +A = 0.02 # input shear amplitude |g| +M_TRUE = 0.05 # multiplicative bias (same for both components) +C1_TRUE = 0.001 # additive bias, component 1 +C2_TRUE = -0.002 # additive bias, component 2 +N_GAL = 2000 + + +def _write_cat(path, ra, dec, e1, e2, w): + """Write a minimal calibrated shape catalogue (RA, Dec, e1, e2, w_des).""" + cols = [ + fits.Column(name=name, array=arr, format="D") + for name, arr in ( + ("RA", ra), + ("Dec", dec), + ("e1", e1), + ("e2", e2), + ("w_des", w), + ) + ] + fits.HDUList([fits.PrimaryHDU(), fits.BinTableHDU.from_columns(cols)]).writeto( + path, overwrite=True + ) + + +def _make_grid(grids_dir, num): + """Create the five sheared/reference catalogues with a known m/c.""" + rng = np.random.default_rng(0) + ra = 30.0 + rng.uniform(0, 0.1, N_GAL) + dec = rng.uniform(0, 0.1, N_GAL) + w = np.ones(N_GAL) + zero = np.zeros(N_GAL) + + def e(g_in): + return (1 + M_TRUE) * g_in + zero + + sims = { + "1z2z": (C1_TRUE + zero, C2_TRUE + zero), + "1p2z": (e(+A) + C1_TRUE, C2_TRUE + zero), + "1m2z": (e(-A) + C1_TRUE, C2_TRUE + zero), + "1z2p": (C1_TRUE + zero, e(+A) + C2_TRUE), + "1z2m": (C1_TRUE + zero, e(-A) + C2_TRUE), + } + for name, (e1, e2) in sims.items(): + sim_dir = grids_dir / f"{name}_grid_{num}" + sim_dir.mkdir(parents=True, exist_ok=True) + _write_cat(sim_dir / "cat.fits", ra, dec, e1, e2, w) + + +def test_match_catalogs_radec_identity(): + """Identical positions match one-to-one; a shifted object drops out.""" + ra = np.array([30.0, 30.01, 30.02]) + dec = np.array([10.0, 10.01, 10.02]) + # Second catalogue = first, but the last object nudged well past threshold. + ra2, dec2 = ra.copy(), dec.copy() + ra2[2] += 1.0 + idx1, idx2 = match_catalogs_radec(ra, dec, ra2, dec2, thresh_deg=0.0002) + npt.assert_array_equal(idx2, [0, 1]) + npt.assert_array_equal(idx1, [0, 1]) + + +def test_mbias_recovers_injected_values(tmp_path): + """ImageSimMBias recovers the injected m/c to machine precision.""" + num = 7 + _make_grid(tmp_path, num) + config = { + "grids_dir": str(tmp_path), + "num": num, + "catalog_name": "cat.fits", + "shear_amplitude": A, + "match_radius_deg": 0.0002, + "w_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 + + +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 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..1ba0b5ca --- /dev/null +++ b/workflow/image_sims/params_im_sim.py @@ -0,0 +1,228 @@ +""" + +:Name: params.py + +:Description: This script contains parameters to run the validation notebook. + +:Author: Martin Kilbinger + +:Date: 2021 + +:Package: sp_validation + +""" + +import os + +import numpy as np + +# Control + +## Verbose output +verbose = True + +## Math output +np.set_printoptions(precision=3, formatter={"float": "{: .3g}".format}) + + +# Survey parameters + +## Field or patch name -- derived from the run directory, which is named +## after the simulation (e.g. '1z2z_grid_1'), so one shared params file +## serves every sim. +name = os.path.basename(os.getcwd()) +print("Field name = {}".format(name)) + +## Area of a tile in deg^2 +area_tile = 0.25 + +## Pixel size in arcsec +pixel_size = 0.187 + +## Shape measurement method, implemented is +## 'ngix': multi-epoch model fitting +## 'galsim': stacked-image moments (experimental) +shape = "ngmix" + +# Paths + +## Input paths + +### Input data directory +data_dir = "." + +### Tile IDs +path_tile_ID = f"{data_dir}/tiles_{name}.txt" + +### Weak-lensing galaxy catalog name +galaxy_cat_path = f"{data_dir}/final_cat_{name}.hdf5" +print(f"Galaxy catalogue = {galaxy_cat_path}") + +## Parameter list; optional, set to `None` if not required +param_list_path = f"{data_dir}/cfis/final_cat.param" + +### Star and PSF catalog name; optional, set to `None` if not required +star_cat_path = None + +# HDU number of star and PSF catalogue +hdu_star_cat = 1 + +### External mask; optional, set to `None` if not required +mask_external_path = None + +## Output paths + +### Output base directory +output_dir = f"{data_dir}" + +### Galaxy shape catalogue base name +output_shape_cat_base = f"{output_dir}/shape_catalog" + +### PSF output catalogue base name. +output_PSF_cat_base = f"{output_dir}/psf_catalog" + +### File for found tile IDs +path_found_ID = f"{output_dir}/found_ID.txt" + +### File for missing tile IDs +path_missing_ID = f"{output_dir}/missing_ID.txt" + +### Plot directory and subdirs +plot_dir = f"{output_dir}/plots/" + +### Statistics text file +stats_file_name = "stats_file.txt" + +# Other IO options + +## Input + +### Coordinate column names +col_name_ra = "XWIN_WORLD" +col_name_dec = "YWIN_WORLD" + +### Memory mode, set to None unless very large file +mmap_mode = None + +## Output + +### Output file format extension: '.fits' or '.hdf5' +output_format = ".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_T_PSF_RECONV_NOSHEAR", +] + +## Pre-calibration catalogue, including masked objects and mask flags. +## ShapePipe-v2 (post-#761) ngmix grammar: ellipticity in named scalar +## components NGMIX_G{1,2}_*, PSF size split into NGMIX_T_PSF_ORIG/RECONV. +## IMAFLAGS_ISO (present in the data-path params) is omitted: the simulation +## pipeline runs no imaging-flag masking stage, so the column does not exist. +## NGMIX_MCAL_TYPES_FAIL is kept -- it is the metacal moments-failure flag the +## calibration mask cuts on, identically to the data path. +add_cols_pre_cal = [ + "TILE_ID", + "NUMBER", + "FLAGS", + "NGMIX_MCAL_FLAGS", + "NGMIX_MCAL_TYPES_FAIL", + "N_EPOCH", + "NGMIX_N_EPOCH", + "NGMIX_G1_PSF_ORIG_NOSHEAR", + "NGMIX_G2_PSF_ORIG_NOSHEAR", + "NGMIX_G1_ERR_NOSHEAR", + "NGMIX_G2_ERR_NOSHEAR", +] + +### Set flag columns as integer format +add_cols_pre_cal_format = {} +for key in ( + "NUMBER", + "FLAGS", + "NGMIX_MCAL_FLAGS", + "NGMIX_MCAL_TYPES_FAIL", + "N_EPOCH", + "NGMIX_N_EPOCH", +): + add_cols_pre_cal_format[key] = "I" + +add_cols_pre_cal_format["TILE_ID"] = "A7" +add_cols_pre_cal_format["NUMBER"] = "J" + +# Create key names for metacal information +prefix = "NGMIX" +suffixes = ["1M", "1P", "2M", "2P", "NOSHEAR"] +centers = ["FLAGS", "G1", "G2", "FLUX", "FLUX_ERR", "T", "T_ERR", "T_PSF_RECONV"] +for center in centers: + for suffix in suffixes: + add_cols_pre_cal.append(f"{prefix}_{center}_{suffix}") + +for suffix in suffixes: + add_cols_pre_cal_format[f"FLAGS_{suffix}"] = "I" + + +# Catalog parameters + +## Star matching threshold [deg] +thresh = 0.0002 + +## Number of jackknife resamples for additive bias +## (0: no jackknife computation). +## If < 2000 the jackknife mean fluctuates a lot. +n_jack = 0 + + +## Galaxy selection + +# Flag to output selected and calibrated galaxy catalogue (<= SP v1.4.1). +# If False, only output comprehensive catalogue. +do_selection_calibration = False + +## Magnitude limits +gal_mag_bright = 15 +gal_mag_faint = 30 + +### Spread-model +do_spread_model = False + +### SExtractor flags to keep in addition to FLAGS=0 +### (bit-coded; list of powers of 2); +### Empty list if no flags +flags_keep = [] + +## Minimum number of epochs +n_epoch_min = 2 + +### Signal-to-noise (selection within metacal) +#### minimum to cut noisy objects +gal_snr_min = 10 +#### maximum to cut too bright objects, potentially too large for the postage stamp +gal_snr_max = 500 + +### Relative size, T_gal / T_psf (selection within metacal) +### to select objects that are not too small compared to the PSF, thus not likely to be point-like, +### or to big as they seem to bias the correlation functions +gal_rel_size_min = 0.5 +gal_rel_size_max = 3.0 + +### Correct galaxy size for ellipticity +gal_size_corr_ell = False + +### prior ellipticity dispersion (one component), *only* used for galaxy weight +sigma_eps_prior = 0.34 + + +## Wrap coordinates around this value [deg], set to != 0 if ra=0 is within coordinate range +wrap_ra = 0 diff --git a/workflow/rules/image_sims.smk b/workflow/rules/image_sims.smk new file mode 100644 index 00000000..7f464e19 --- /dev/null +++ b/workflow/rules/image_sims.smk @@ -0,0 +1,314 @@ +"""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 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", 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" + +# --- 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: 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]) + + +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}.fits", + 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, + and the raw SKiLLS image inputs. + + ``params_im_sim.py`` derives the field name from the directory basename, so + the same template serves every sim; ``config_mask.yaml`` and ``cfis`` are + symlinks the downstream calibration and merge steps read from cwd. + + ``input_tiles``/``input_exp`` are top-level symlinks to the raw SKiLLS tile + and exposure images; ShapePipe's ``get_images_runner`` resolves them via + ``$SP_DIR/input_{tiles,exp}`` (``$SP_DIR`` is the run dir). ``run_job`` does + not stage these, so ``im_init`` must -- this is what makes ``im_pipeline`` + runnable from raw images, not just from pre-staged intermediates. + """ + input: + # Tracked so that editing the params template or mask config re-stages + # them into every run dir (a plain params: value would not retrigger, + # silently leaving stale params.py behind after a grammar change). + template=PARAMS_TEMPLATE, + mask_src=os.path.join(SPV_REPO, MASK_CONFIG), + output: + params=f"{GRIDS_BASE}/{{sim}}/params.py", + mask=f"{GRIDS_BASE}/{{sim}}/config_mask.yaml", + params: + config_dir=CONFIG_DIR, + run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", + cfis=lambda wc: f"{GRIDS_BASE}/{wc.sim}/cfis", + sim_tiles=lambda wc: f"{INPUT_SIMS_BASE}/{wc.sim}/images/SP_tiles", + sim_exp=lambda wc: f"{INPUT_SIMS_BASE}/{wc.sim}/images/SP_exp", + shell: + # cfis / input_tiles / input_exp are stable read-only symlinks (used by + # get_images, merge, extract); created here but not tracked as outputs, + # which snakemake will not accept for a symlink/directory. + "mkdir -p $(dirname {output.params}) && " + "cp {input.template} {output.params} && " + "ln -sf {input.mask_src} {output.mask} && " + "ln -sfT {params.config_dir} {params.cfis} && " + "ln -sfT {params.sim_tiles} {params.run_dir}/input_tiles && " + "ln -sfT {params.sim_exp} {params.run_dir}/input_exp" + + +rule im_pipeline: + """Run ShapePipe on one simulated tile (ShapePipe 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}.fits", + 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}.fits", + mask=f"{GRIDS_BASE}/{{sim}}/config_mask.yaml", + output: + cat=f"{GRIDS_BASE}/{{sim}}/shape_catalog_cut_{SHAPE}.fits", + params: + run_dir=lambda wc: f"{GRIDS_BASE}/{wc.sim}", + shell: + "cd {params.run_dir} && " + "{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" + )