From c90ae4ff5bab50936a63c29443d9176bfb684705 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Wed, 25 Feb 2026 16:23:46 +0100 Subject: [PATCH 01/13] additive bias calculation for paper updated to v1.4.6.3 --- notebooks/cosmo_val/cat_config.yaml | 276 ++++++++++++++++++++++++++++ pyproject.toml | 4 + 2 files changed, 280 insertions(+) diff --git a/notebooks/cosmo_val/cat_config.yaml b/notebooks/cosmo_val/cat_config.yaml index a1a6b756..6596de22 100644 --- a/notebooks/cosmo_val/cat_config.yaml +++ b/notebooks/cosmo_val/cat_config.yaml @@ -726,6 +726,7 @@ SP_v1.4.8: e1_col: e1 e2_col: e2 path: unions_shapepipe_star_2024_v1.4.a.fits +<<<<<<< Updated upstream # For additive bias only SP_v1.4.6_uncal: @@ -830,6 +831,281 @@ SP_v1.4.8_uncal: hdu: 1 patch_number: 100 +======= +SP_v1.4.11.2: + subdir: /n17data/UNIONS/WL/v1.4.x + pipeline: SP + colour: orange + getdist_colour: 0.0, 0.5, 1.0 + ls: dashed + marker: s + cov_th: + A: 2894.0303815287743 + n_e: 6.25812596800506 + n_psf: 0.752316232272063 + sigma_e: 0.380499698062307 + mask: /home/guerrini/sp_validation/cosmo_inference/data/mask/mask_map_footprint_nside_4096.fits + psf: + PSF_flag: FLAG_PSF_HSM + PSF_size: SIGMA_PSF_HSM + square_size: true + star_flag: FLAG_STAR_HSM + star_size: SIGMA_STAR_HSM + hdu: 1 + path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + ra_col: RA + dec_col: Dec + e1_PSF_col: E1_PSF_HSM + e1_star_col: E1_STAR_HSM + e2_PSF_col: E2_PSF_HSM + e2_star_col: E2_STAR_HSM + shear: + R: 1.0 + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + path: /n17data/UNIONS/WL/v1.4.x/v1.4.11.2/unions_shapepipe_cut_struc_2024_v1.4.11.2.fits + redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + w_col: w_des + e1_col: e1 + e1_col_corrected: e1_leak_corrected + e1_PSF_col: e1_PSF + e2_col: e2 + e2_col_corrected: e2_leak_corrected + e2_PSF_col: e2_PSF + star: + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + path: unions_shapepipe_star_2024_v1.4.a.fits +SP_v1.4.11.3: + subdir: /n17data/UNIONS/WL/v1.4.x + pipeline: SP + colour: red + getdist_colour: 0.0, 0.5, 1.0 + ls: dashdot + marker: v + cov_th: + A: 2405.3892055695346 + n_e: 6.128201234871523 + n_psf: 0.752316232272063 + sigma_e: 0.379587601488189 + mask: /home/guerrini/sp_validation/cosmo_inference/data/mask/mask_map_v1.4.6_nside_8192.fits + psf: + PSF_flag: FLAG_PSF_HSM + PSF_size: SIGMA_PSF_HSM + square_size: true + star_flag: FLAG_STAR_HSM + star_size: SIGMA_STAR_HSM + hdu: 1 + path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + ra_col: RA + dec_col: Dec + e1_PSF_col: E1_PSF_HSM + e1_star_col: E1_STAR_HSM + e2_PSF_col: E2_PSF_HSM + e2_star_col: E2_STAR_HSM + shear: + R: 1.0 + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + path: v1.4.11.3/unions_shapepipe_cut_struc_2024_v1.4.11.3.fits + redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + w_col: w_des + e1_col: e1 + e1_col_corrected: e1_leak_corrected + e1_PSF_col: e1_PSF + e2_col: e2 + e2_col_corrected: e2_leak_corrected + e2_PSF_col: e2_PSF + star: + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits +SP_v1.4.12.3: + subdir: /n17data/UNIONS/WL/v1.4.x + pipeline: SP + colour: lightblue + getdist_colour: 0.0, 0.5, 1.0 + ls: dashdot + marker: ^ + cov_th: + A: 2405.3892055695346 + n_e: 6.128201234871523 + n_psf: 0.752316232272063 + sigma_e: 0.379587601488189 + mask: /home/guerrini/sp_validation/cosmo_inference/data/mask/mask_map_v1.4.6_nside_8192.fits + psf: + PSF_flag: FLAG_PSF_HSM + PSF_size: SIGMA_PSF_HSM + square_size: true + star_flag: FLAG_STAR_HSM + star_size: SIGMA_STAR_HSM + hdu: 1 + path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + ra_col: RA + dec_col: Dec + e1_PSF_col: E1_PSF_HSM + e1_star_col: E1_STAR_HSM + e2_PSF_col: E2_PSF_HSM + e2_star_col: E2_STAR_HSM + shear: + R: 1.0 + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + path: v1.4.12.3/unions_shapepipe_cut_struc_2024_v1.4.12.3.fits + redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + w_col: w_des + e1_col: e1 + e1_col_corrected: e1_leak_corrected + e1_PSF_col: e1_PSF + e2_col: e2 + e2_col_corrected: e2_leak_corrected + e2_PSF_col: e2_PSF + star: + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits +SP_v1.4.13.3: + subdir: /n17data/UNIONS/WL/v1.4.x + pipeline: SP + colour: cyan + getdist_colour: 0.0, 0.5, 1.0 + ls: dashdot + marker: h + cov_th: + A: 2405.3892055695346 + n_e: 6.128201234871523 + n_psf: 0.752316232272063 + sigma_e: 0.379587601488189 + mask: /home/guerrini/sp_validation/cosmo_inference/data/mask/mask_map_v1.4.6_nside_8192.fits + psf: + PSF_flag: FLAG_PSF_HSM + PSF_size: SIGMA_PSF_HSM + square_size: true + star_flag: FLAG_STAR_HSM + star_size: SIGMA_STAR_HSM + hdu: 1 + path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + ra_col: RA + dec_col: Dec + e1_PSF_col: E1_PSF_HSM + e1_star_col: E1_STAR_HSM + e2_PSF_col: E2_PSF_HSM + e2_star_col: E2_STAR_HSM + shear: + R: 1.0 + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + path: v1.4.13.3/unions_shapepipe_cut_struc_2024_v1.4.13.3.fits + redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + w_col: w_des + e1_col: e1 + e1_col_corrected: e1_leak_corrected + e1_PSF_col: e1_PSF + e2_col: e2 + e2_col_corrected: e2_leak_corrected + e2_PSF_col: e2_PSF + star: + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits +SP_v1.4.11.3_ecut07: + subdir: /n17data/UNIONS/WL/v1.4.x + pipeline: SP + colour: darkblue + getdist_colour: 0.0, 0.3, 0.7 + ls: dotted + marker: s + cov_th: + A: 2894.0303815287743 + n_e: 5.830114732534133 + n_psf: 0.752316232272063 + sigma_e: 0.34058513153783426 + mask: /home/guerrini/sp_validation/cosmo_inference/data/mask/mask_map_footprint_nside_4096.fits + psf: + PSF_flag: FLAG_PSF_HSM + PSF_size: SIGMA_PSF_HSM + square_size: true + star_flag: FLAG_STAR_HSM + star_size: SIGMA_STAR_HSM + hdu: 1 + path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + ra_col: RA + dec_col: Dec + e1_PSF_col: E1_PSF_HSM + e1_star_col: E1_STAR_HSM + e2_PSF_col: E2_PSF_HSM + e2_star_col: E2_STAR_HSM + shear: + R: 1.0 + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + path: /n17data/cdaley/unions/pure_eb/results/ecut/SP_v1.4.11.3_ecut07.fits + redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + w_col: w_des + e1_col: e1 + e1_col_corrected: e1_leak_corrected + e1_PSF_col: e1_PSF + e2_col: e2 + e2_col_corrected: e2_leak_corrected + e2_PSF_col: e2_PSF + star: + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits +SP_v1.4.6.3_uncal: + pipeline: SP + subdir: /n17data/UNIONS/WL/v1.4.x + shear: + path: v1.4.6.3/unions_shapepipe_cut_struc_2024_v1.4.6.3.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1_uncal + e2_col: e2_uncal + w_col: w_des + R: 1.0 + psf: + path: unions_shapepipe_psf_2024_v1.4.a.fits + hdu: 1 + patch_number: 100 +SP_v1.4.6.3_uncal_w_iv: + pipeline: SP + subdir: /n17data/UNIONS/WL/v1.4.x + shear: + path: v1.4.6.3/unions_shapepipe_cut_struc_2024_v1.4.6.3.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1_uncal + e2_col: e2_uncal + w_col: w_iv + R: 1.0 + psf: + path: unions_shapepipe_psf_2024_v1.4.a.fits + hdu: 1 + patch_number: 100 +SP_v1.4.6.3_uncal_w_1: + pipeline: SP + subdir: /n17data/UNIONS/WL/v1.4.x + shear: + path: v1.4.6.3/unions_shapepipe_cut_struc_2024_v1.4.6.3.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1_uncal + e2_col: e2_uncal + w_col: one + R: 1.0 + psf: + path: unions_shapepipe_psf_2024_v1.4.a.fits + hdu: 1 + patch_number: 100 +>>>>>>> Stashed changes SP_v1.4_LFmask_8k: subdir: /n17data/mkilbing/astro/data/CFIS/v1.0/SP_LFmask pipeline: SP diff --git a/pyproject.toml b/pyproject.toml index d0ee1d58..f706548b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,7 +35,11 @@ dependencies = [ "lenspack", "lmfit", "numexpr", +<<<<<<< Updated upstream "numpy>=1.14", +======= + "numpy>=1.26", +>>>>>>> Stashed changes "opencv-python-headless", "pyccl", "pyarrow", From 392461e7eb3d79f3c93849aad74aed12bb6f3b40 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 20 Apr 2026 11:23:30 +0200 Subject: [PATCH 02/13] added config --- notebooks/cosmo_val/cat_config.yaml | 69 +++-------------------------- 1 file changed, 6 insertions(+), 63 deletions(-) diff --git a/notebooks/cosmo_val/cat_config.yaml b/notebooks/cosmo_val/cat_config.yaml index 6596de22..2c95f9f8 100644 --- a/notebooks/cosmo_val/cat_config.yaml +++ b/notebooks/cosmo_val/cat_config.yaml @@ -726,14 +726,11 @@ SP_v1.4.8: e1_col: e1 e2_col: e2 path: unions_shapepipe_star_2024_v1.4.a.fits -<<<<<<< Updated upstream - -# For additive bias only -SP_v1.4.6_uncal: +SP_v1.4.6.3_uncal: pipeline: SP subdir: /n17data/UNIONS/WL/v1.4.x shear: - path: v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + path: v1.4.6.3/unions_shapepipe_cut_struc_2024_v1.4.6.3.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt ra_col: RA dec_col: Dec @@ -746,11 +743,11 @@ SP_v1.4.6_uncal: hdu: 1 patch_number: 100 -SP_v1.4.6_uncal_w_iv: +SP_v1.4.6.3_uncal_w_iv: pipeline: SP subdir: /n17data/UNIONS/WL/v1.4.x shear: - path: v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + path: v1.4.6.3/unions_shapepipe_cut_struc_2024_v1.4.6.3.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt ra_col: RA dec_col: Dec @@ -763,11 +760,11 @@ SP_v1.4.6_uncal_w_iv: hdu: 1 patch_number: 100 -SP_v1.4.6_uncal_w_1: +SP_v1.4.6.3_uncal_w_1: pipeline: SP subdir: /n17data/UNIONS/WL/v1.4.x shear: - path: v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + path: v1.4.6.3/unions_shapepipe_cut_struc_2024_v1.4.6.3.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt ra_col: RA dec_col: Dec @@ -779,59 +776,6 @@ SP_v1.4.6_uncal_w_1: path: unions_shapepipe_psf_2024_v1.4.a.fits hdu: 1 patch_number: 100 - -SP_v1.4.5_uncal: - pipeline: SP - subdir: /n17data/UNIONS/WL/v1.4.x - shear: - path: v1.4.5/unions_shapepipe_cut_struc_2024_v1.4.5.fits - covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt - ra_col: RA - dec_col: Dec - e1_col: e1_uncal - e2_col: e2_uncal - w_col: w_des - R: 1.0 - psf: - path: unions_shapepipe_psf_2024_v1.4.a.fits - hdu: 1 - patch_number: 100 - -SP_v1.4.7_uncal: - pipeline: SP - subdir: /n17data/UNIONS/WL/v1.4.x - shear: - path: v1.4.7/unions_shapepipe_cut_struc_2024_v1.4.7.fits - covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt - ra_col: RA - dec_col: Dec - e1_col: e1_uncal - e2_col: e2_uncal - w_col: w_des - R: 1.0 - psf: - path: unions_shapepipe_psf_2024_v1.4.a.fits - hdu: 1 - patch_number: 100 - -SP_v1.4.8_uncal: - pipeline: SP - subdir: /n17data/UNIONS/WL/v1.4.x - shear: - path: v1.4.8/unions_shapepipe_cut_struc_2024_v1.4.8.fits - covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt - ra_col: RA - dec_col: Dec - e1_col: e1_uncal - e2_col: e2_uncal - w_col: w_des - R: 1.0 - psf: - path: unions_shapepipe_psf_2024_v1.4.a.fits - hdu: 1 - patch_number: 100 - -======= SP_v1.4.11.2: subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP @@ -1105,7 +1049,6 @@ SP_v1.4.6.3_uncal_w_1: path: unions_shapepipe_psf_2024_v1.4.a.fits hdu: 1 patch_number: 100 ->>>>>>> Stashed changes SP_v1.4_LFmask_8k: subdir: /n17data/mkilbing/astro/data/CFIS/v1.0/SP_LFmask pipeline: SP From 0ec3e3f21577708c474e5ff69450db9549dfe6f2 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 20 Apr 2026 11:30:53 +0200 Subject: [PATCH 03/13] added fill_photoz script --- scripts/fill_photoz_bands.py | 509 +++++++++++++++++++++++++++++++++++ 1 file changed, 509 insertions(+) create mode 100755 scripts/fill_photoz_bands.py diff --git a/scripts/fill_photoz_bands.py b/scripts/fill_photoz_bands.py new file mode 100755 index 00000000..83865b15 --- /dev/null +++ b/scripts/fill_photoz_bands.py @@ -0,0 +1,509 @@ +#!/usr/bin/env python + +"""fill_photoz_bands.py + +Add PhotoPipe photo-z and multi-band magnitude fields to a ShapePipe +comprehensive HDF5 catalogue and fill them from PhotoPipe FITS tiles. + +Two phases, both resumable via a checkpoint file: + + Phase 1 — Create output: copy the input HDF5 and append the new + columns initialised to EMPTY_VALUE (-199). + Phase 2 — Fill tiles: for each PhotoPipe FITS tile found in fits_dir, + write the field values into the output file. + +Skips missing FITS files. Warns (does not error) on missing PhotoPipe +keys. Supports interrupt + restart at any point. + +:Authors: Martin Kilbinger + +""" + +import json +import os +import sys +import warnings +from collections import defaultdict +from timeit import default_timer as timer + +import h5py +import numpy as np +import tqdm +from astropy.io import fits + +from cs_util import args as cs_args +from cs_util import logging + + +FITS_HDU = 1 +EMPTY_VALUE = -199 +COPY_CHUNK = 2_000_000 # rows per chunk when copying input → output +SCAN_CHUNK = 5_000_000 # rows per chunk when scanning TILE_ID + +REQUESTED_KEYS = [ + "Z_B", + "Z_B_MIN", + "Z_B_MAX", + "T_B", + "MAG_GAAP_0p7_u", + "MAG_GAAP_0p7_g", + "MAG_GAAP_0p7_r", + "MAG_GAAP_0p7_i", + "MAG_GAAP_0p7_z", + "MAG_GAAP_0p7_z2", +] + + +def params_default(): + """Params Default. + + Return default parameter values and additional information + about type and command line options. + + Returns + ------- + tuple + parameter dict, short_options dict, types dict, help_strings dict + + """ + params = { + "input": "unions_shapepipe_comprehensive_struc_2024_v1.5.c.hdf5", + "output": "unions_shapepipe_comprehensive_struc_ugriz_2024_v1.5.c.hdf5", + "fits_dir": "UNIONS5000", + "checkpoint": "fill_photoz_bands_checkpoint.json", + "verbose": False, + } + + short_options = { + "input": "-i", + "output": "-o", + "fits_dir": "-d", + "checkpoint": "-c", + } + + types = {} + + help_strings = { + "input": "input HDF5 catalogue (no PhotoPipe fields), default={}", + "output": "output HDF5 catalogue (PhotoPipe fields added and filled), default={}", + "fits_dir": "directory with PhotoPipe FITS tiles, default={}", + "checkpoint": "checkpoint JSON file for resume support, default={}", + } + + return params, short_options, types, help_strings + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def detect_dataset_name(hf): + """Return the first dataset name in an HDF5 file. + + Tries common names first: dat, dat_comb, data. + + Parameters + ---------- + hf : h5py.File + + Returns + ------- + str + Dataset name. + + Raises + ------ + KeyError + If no dataset is found. + + """ + for name in ("dat", "dat_comb", "data"): + if name in hf: + return name + # Fall back to first key + keys = list(hf.keys()) + if not keys: + raise KeyError("HDF5 file contains no datasets.") + return keys[0] + + +def strip_dtype(dtype): + """Strip h5py metadata from a structured numpy dtype descriptor. + + h5py sometimes adds encoding metadata to dtype descriptions, e.g. + ``(' 1)[0] + 1 + blocks = np.split(sorted_idx, gaps) + fits_offset = 0 + for block in blocks: + b_min, b_max = int(block[0]), int(block[-1]) + n_block = b_max - b_min + 1 + chunk = dset[b_min : b_max + 1] + for key in valid_keys: + chunk[key] = fits_data[key][fits_offset : fits_offset + n_block] + dset[b_min : b_max + 1] = chunk + fits_offset += n_block + + +# --------------------------------------------------------------------------- +# Phase 1: create output file +# --------------------------------------------------------------------------- + +def create_output_file(input_path, output_path, dataset_name, verbose=False): + """Create output HDF5 by copying input and appending empty PhotoPipe fields. + + Parameters + ---------- + input_path : str + output_path : str + dataset_name : str + Dataset name in the input file (used for output too). + verbose : bool + + """ + print(f"Phase 1: creating output file '{output_path}'") + t0 = timer() + + with h5py.File(input_path, "r") as hf_in: + dset_in = hf_in[dataset_name] + n_total = dset_in.shape[0] + dtype_out = build_output_dtype(dset_in.dtype, REQUESTED_KEYS) + new_keys = [k for k in REQUESTED_KEYS if k not in set(dset_in.dtype.names)] + + print(f" Input rows : {n_total:,}") + print(f" Input fields : {len(dset_in.dtype.names)}") + print(f" New fields : {new_keys}") + size_mb = n_total * dtype_out.itemsize / 1_048_576 + print(f" Output size : ~{size_mb:,.0f} MB") + + with h5py.File(output_path, "w") as hf_out: + dset_out = hf_out.create_dataset( + dataset_name, + shape=(n_total,), + dtype=dtype_out, + ) + + # Copy input fields chunk by chunk + input_fields = dset_in.dtype.names + with tqdm.tqdm( + total=n_total, unit="rows", unit_scale=True, desc=" Copying" + ) as pbar: + for start in range(0, n_total, COPY_CHUNK): + end = min(start + COPY_CHUNK, n_total) + chunk_in = dset_in[start:end] + chunk_out = np.empty(end - start, dtype=dtype_out) + for field in input_fields: + chunk_out[field] = chunk_in[field] + for key in new_keys: + chunk_out[key] = EMPTY_VALUE + dset_out[start:end] = chunk_out + pbar.update(end - start) + + elapsed = timer() - t0 + print(f" Done in {elapsed:.1f}s\n") + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(argv=None): + """Main. + + Main program. + + """ + params, short_options, types, help_strings = params_default() + + options = cs_args.parse_options(params, short_options, types, help_strings) + params.update(options) + + logging.log_command(argv) + + verbose = params["verbose"] + + if not os.path.exists(params["input"]): + print(f"ERROR: input file not found: {params['input']}", file=sys.stderr) + return 1 + if not os.path.isdir(params["fits_dir"]): + print(f"ERROR: FITS directory not found: {params['fits_dir']}", file=sys.stderr) + return 1 + + # ------------------------------------------------------------------ + # Load checkpoint + # ------------------------------------------------------------------ + if os.path.exists(params["checkpoint"]): + with open(params["checkpoint"]) as f: + checkpoint = json.load(f) + done_tiles = set(checkpoint.get("done_tiles", [])) + print(f"Resuming: {len(done_tiles)} tiles already completed.") + else: + done_tiles = set() + checkpoint = {} + + # ------------------------------------------------------------------ + # Detect input dataset name (cache in checkpoint) + # ------------------------------------------------------------------ + if "dataset_name" in checkpoint: + dataset_name = checkpoint["dataset_name"] + else: + with h5py.File(params["input"], "r") as hf: + dataset_name = detect_dataset_name(hf) + checkpoint["dataset_name"] = dataset_name + with open(params["checkpoint"], "w") as cf: + json.dump(checkpoint, cf) + + if verbose: + print(f"Input dataset: '{dataset_name}'") + + # ------------------------------------------------------------------ + # Phase 1: create output file if needed + # ------------------------------------------------------------------ + if not checkpoint.get("output_created", False): + if os.path.exists(params["output"]): + print( + f"WARNING: output file '{params['output']}' exists but checkpoint " + "does not mark it as complete. Overwriting." + ) + create_output_file( + params["input"], params["output"], dataset_name, verbose=verbose + ) + checkpoint["output_created"] = True + # Invalidate any stale tile_index_map from a previous partial run + checkpoint.pop("tile_index_map", None) + with open(params["checkpoint"], "w") as cf: + json.dump(checkpoint, cf) + else: + if verbose: + print("Phase 1 already done (output file exists in checkpoint).") + + # ------------------------------------------------------------------ + # Phase 2: fill tiles + # ------------------------------------------------------------------ + t0 = timer() + with h5py.File(params["output"], "r+") as hf: + dset = hf[dataset_name] + n_total = dset.shape[0] + print(f"Phase 2: filling tiles in '{params['output']}'") + print(f" {n_total:,} rows, dataset '{dataset_name}'") + + hdf5_fields = set(dset.dtype.names) + + # Check all requested keys are present + for key in REQUESTED_KEYS: + if key not in hdf5_fields: + warnings.warn( + f"Key '{key}' missing from output dataset — was Phase 1 complete?" + ) + + # Build TILE_ID → output row indices (scan once; cached in checkpoint) + if "tile_index_map" in checkpoint: + if verbose: + print(" Loading tile→index map from checkpoint...") + tile_index_map = { + k.encode(): np.array(v, dtype=np.int64) + for k, v in checkpoint["tile_index_map"].items() + } + else: + print(" Building tile→index map (scans all rows)...") + tile_index_map_lists = defaultdict(list) + with tqdm.tqdm( + total=n_total, unit="rows", unit_scale=True, desc=" Scanning TILE_ID" + ) as pbar: + for start in range(0, n_total, SCAN_CHUNK): + end = min(start + SCAN_CHUNK, n_total) + tile_chunk = dset["TILE_ID"][start:end] + for local_i, tid in enumerate(tile_chunk): + tile_index_map_lists[tid].append(start + local_i) + pbar.update(end - start) + + tile_index_map = { + tid: np.array(idxs, dtype=np.int64) + for tid, idxs in tile_index_map_lists.items() + } + checkpoint["tile_index_map"] = { + tid.decode(): idxs.tolist() + for tid, idxs in tile_index_map.items() + } + with open(params["checkpoint"], "w") as cf: + json.dump(checkpoint, cf) + print(f" Map saved: {len(tile_index_map)} unique tiles.") + + unique_tiles = sorted(tile_index_map.keys()) + n_tiles = len(unique_tiles) + + valid_keys = None # determined from first available FITS tile + n_skipped_missing = 0 + n_skipped_done = 0 + n_skipped_size = 0 + n_processed = 0 + + print(f"\n Processing {n_tiles} tiles ({len(done_tiles)} already done)...\n") + pbar = tqdm.tqdm(unique_tiles, total=n_tiles, unit="tile") + + for tile_id in pbar: + tile_str = tile_id.decode() if isinstance(tile_id, bytes) else tile_id + + if tile_str in done_tiles: + n_skipped_done += 1 + continue + + fits_path = tile_id_to_fits_path(tile_id, params["fits_dir"]) + if not os.path.exists(fits_path): + n_skipped_missing += 1 + done_tiles.add(tile_str) + pbar.set_postfix({"done": n_processed, "missing": n_skipped_missing}) + continue + + try: + with fits.open(fits_path, memmap=True) as hdu_list: + fits_data = hdu_list[FITS_HDU].data + + # Validate keys on first successfully opened tile + if valid_keys is None: + valid_keys, missing_keys = check_fits_keys( + fits_data.dtype.names, REQUESTED_KEYS + ) + valid_keys = [k for k in valid_keys if k in hdf5_fields] + if missing_keys: + warnings.warn( + f"Keys absent from PhotoPipe FITS (skipped): {missing_keys}" + ) + tqdm.tqdm.write(f"\n Keys to fill: {valid_keys}\n") + + hdf5_indices = tile_index_map[tile_id] + n_hdf5 = len(hdf5_indices) + n_fits = len(fits_data) + + if n_hdf5 != n_fits: + warnings.warn( + f"Tile {tile_str}: HDF5 has {n_hdf5} rows, FITS has " + f"{n_fits} — size mismatch, skipping." + ) + n_skipped_size += 1 + pbar.set_postfix( + {"done": n_processed, "size_err": n_skipped_size} + ) + continue + + write_tile_to_hdf5(dset, hdf5_indices, fits_data, valid_keys) + + except Exception as e: + warnings.warn(f"Tile {tile_str}: error ({e}), skipping.") + continue + + n_processed += 1 + done_tiles.add(tile_str) + + if n_processed % 50 == 0: + checkpoint["done_tiles"] = list(done_tiles) + with open(params["checkpoint"], "w") as cf: + json.dump(checkpoint, cf) + + pbar.set_postfix({"done": n_processed, "missing": n_skipped_missing}) + + # Final checkpoint flush + checkpoint["done_tiles"] = list(done_tiles) + with open(params["checkpoint"], "w") as cf: + json.dump(checkpoint, cf) + + elapsed = timer() - t0 + print(f"\nDone in {elapsed:.1f}s") + print(f" Tiles processed : {n_processed}") + print(f" Tiles skipped (done) : {n_skipped_done}") + print(f" FITS files missing : {n_skipped_missing}") + print(f" Size mismatches : {n_skipped_size}") + + return 0 + + +if __name__ == "__main__": + sys.exit(main(sys.argv)) From 7f29a8b58dcbdfc00d601a8cd489ca3b99357225 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Tue, 16 Jun 2026 13:41:30 +0200 Subject: [PATCH 04/13] fill photoz bands fixes --- notebooks/create_shear_mb_empty.py | 172 --------------------------- notebooks/demo_add_bands_to_empty.py | 17 +-- scripts/fill_photoz_bands.py | 59 ++++----- 3 files changed, 34 insertions(+), 214 deletions(-) delete mode 100644 notebooks/create_shear_mb_empty.py diff --git a/notebooks/create_shear_mb_empty.py b/notebooks/create_shear_mb_empty.py deleted file mode 100644 index 6d29bfc2..00000000 --- a/notebooks/create_shear_mb_empty.py +++ /dev/null @@ -1,172 +0,0 @@ -# --- -# jupyter: -# jupytext: -# text_representation: -# extension: .py -# format_name: light -# format_version: '1.5' -# jupytext_version: 1.15.1 -# kernelspec: -# display_name: sp_validation -# language: python -# name: sp_validation -# --- - -# # Demo notebook to add (u,g,i,z,z2) bands to an r-band catalogue - -# %reload_ext autoreload -# %autoreload 2 - -# + -import os -import numpy as np -import numpy.lib.recfunctions as rfn -import h5py - -from timeit import default_timer as timer -import tqdm -import healsparse as hsp -from astropy.io import fits - -from sp_validation import run_joint_cat as sp_joint -# - - -# Create instance of object -obj = sp_joint.BaseCat() - - -# + -# Set parameters -base = "unions_shapepipe_comprehensive_struc" -year = 2024 -ver = "v1.6.c" - -obj._params = {} - -obj._params["input_path"] = f"{base}_{year}_{ver}.hdf5" -obj._params["verbose"] = True - -# + -path_bands = "./UNIONS5000" -subdir_base = "UNIONS." - -path_base = subdir_base -path_suff = "_SP_ugriz_photoz_ext.cat" - -# NUMBER key in photo-z catalogue -key_num = "SeqNr" - -keys_mag = [f"MAG_GAAP_0p7_{band}" for band in ("u", "g", "r", "i", "z", "z2")] - -keys = ["Z_B", "Z_B_MIN", "Z_B_MAX", "T_B"] + keys_mag - -hdu_no = 1 -# - - -# ## Run - -# + -# Check parameter validity -#obj.check_params() - -# Update parameters (here: strings to list) -#obj.update_params() -# - - -# Read catalogue -dat = obj.read_cat(load_into_memory=False, mode="r") -n_rows = len(dat) - - -def get_dtype_keys(keys,path=None, hdu_no=1): - - if path is None: - - dtype = np.dtype([(key, np.float32) for key in keys]) - - else: - - print(" Read data from file:", path, end=" ") - start = timer() - hdu_list = fits.open(path) - dat_mb = hdu_list[hdu_no].data - dtype = np.dtype([dt for dt in dat_mb.dtype.descr if dt[0] in keys]) - end = timer() - print(f" {end - start:.1f}s") - - return dtype - - -# + -# Get dtype of new keys - -path = None -#path = os.path.join(path_bands, f"{path_base}{tile_ID}", f"{path_base}{tile_IDs[0]}{path_suff}") - -dtype_keys = get_dtype_keys(keys, path=path, hdu_no=hdu_no) - - -# - - -def strip_h5py_metadata_dtype(dat_dtype, dat_ext_dtype): - cleaned_fields = [] - for name, dt in dat_dtype.descr + dat_ext_dtype.descr: - # If dt is a tuple (e.g., ('S7', {'h5py_encoding': 'ascii'})) - if isinstance(dt, tuple): - cleaned_fields.append((name, dt[0])) # keep only the base dtype string - else: - cleaned_fields.append((name, dt)) # use as-is - return cleaned_fields - - -# + -# Create empty array with new keys -# Initialise with -199 to later be able to check for unfilled values - -total_bytes = n_rows * np.dtype(dtype_keys).itemsize -print(" Create new combined array.", end=" ") -print(f"Expected size = {total_bytes / 1_048_576:.2f} MB", end=" ") -start = timer() - -obj._params["output_path"] = f"{base}_empty_ugriz_{year}_{ver}.hdf5" -dtype_sp = dat.dtype -dtype_comb = strip_h5py_metadata_dtype(dtype_sp, dtype_keys) -with h5py.File(obj._params["output_path"], "w") as f: - - # Create new dataset - dset_comb = f.create_dataset( - "dat_comb", - shape=(n_rows,), - dtype=dtype_comb, - ) - - # Copy old data field-by-field - for name in dtype_sp.names: - dset_comb[name] = dat[name] - - # Fill new fields with default value (-199) - for name in dtype_keys.names: - dset_comb[name] = -199 - -#new_empty = np.full(n_rows, -199, dtype=dtype_keys) - -end = timer() -print(f" {end - start:.1f}s") - - - -# + -# Merge with original data - -#print(" Merge empty to original", end=" ") -#start = timer() -#combined = rfn.merge_arrays([dat, new_empty], flatten=True) -#end = timer() -#print(f" {end - start:.1f}s") -# - - -# obj.write_hdf5_file(combined) - - -# Close input HDF5 catalogue file -# obj.close_hd5() diff --git a/notebooks/demo_add_bands_to_empty.py b/notebooks/demo_add_bands_to_empty.py index 14788826..3ea24ed6 100644 --- a/notebooks/demo_add_bands_to_empty.py +++ b/notebooks/demo_add_bands_to_empty.py @@ -43,7 +43,7 @@ # Set parameters base = "unions_shapepipe_comprehensive_struc_empty_ugriz" year = 2024 -ver = "v1.5.c" +ver = "v1.6.c" obj._params = {} @@ -59,19 +59,20 @@ bands = ("u", "g", "r", "i", "z", "z2") -base_keys = ["MAGERR_GAAP", "FLUX_GAAP", "FLUXERR_GAAP", "FLAG_GAAP", "MAG_LIM"] keys_mag = [f"MAG_GAAP_0p7_{band}" for band in bands] -for base_key in base_keys: - keys_mag.extend([f"_{base_key}_{band}" for band in bands]) +#base_keys = ["MAGERR_GAAP", "FLUX_GAAP", "FLUXERR_GAAP", "FLAG_GAAP", "MAG_LIM"] +#for base_key in base_keys: + #keys_mag.extend([f"{base_key}_{band}" for band in bands]) + +# "EXTINCTION", +# "MP_NAME", +# "ODDS", keys = [ - "EXTINCTION", - "MP_NAME", "Z_B", "Z_B_MIN", - "Z_B_MAX" + "Z_B_MAX", "T_B", - "ODDS", ] + keys_mag hdu_no = 1 diff --git a/scripts/fill_photoz_bands.py b/scripts/fill_photoz_bands.py index 83865b15..31ab290e 100755 --- a/scripts/fill_photoz_bands.py +++ b/scripts/fill_photoz_bands.py @@ -287,6 +287,13 @@ def create_output_file(input_path, output_path, dataset_name, verbose=False): dset_out[start:end] = chunk_out pbar.update(end - start) + # Copy all other top-level datasets/groups unchanged + for key in hf_in.keys(): + if key != dataset_name: + hf_in.copy(key, hf_out) + if verbose: + print(f" Copied group/dataset '{key}' unchanged.") + elapsed = timer() - t0 print(f" Done in {elapsed:.1f}s\n") @@ -357,8 +364,6 @@ def main(argv=None): params["input"], params["output"], dataset_name, verbose=verbose ) checkpoint["output_created"] = True - # Invalidate any stale tile_index_map from a previous partial run - checkpoint.pop("tile_index_map", None) with open(params["checkpoint"], "w") as cf: json.dump(checkpoint, cf) else: @@ -384,38 +389,24 @@ def main(argv=None): f"Key '{key}' missing from output dataset — was Phase 1 complete?" ) - # Build TILE_ID → output row indices (scan once; cached in checkpoint) - if "tile_index_map" in checkpoint: - if verbose: - print(" Loading tile→index map from checkpoint...") - tile_index_map = { - k.encode(): np.array(v, dtype=np.int64) - for k, v in checkpoint["tile_index_map"].items() - } - else: - print(" Building tile→index map (scans all rows)...") - tile_index_map_lists = defaultdict(list) - with tqdm.tqdm( - total=n_total, unit="rows", unit_scale=True, desc=" Scanning TILE_ID" - ) as pbar: - for start in range(0, n_total, SCAN_CHUNK): - end = min(start + SCAN_CHUNK, n_total) - tile_chunk = dset["TILE_ID"][start:end] - for local_i, tid in enumerate(tile_chunk): - tile_index_map_lists[tid].append(start + local_i) - pbar.update(end - start) - - tile_index_map = { - tid: np.array(idxs, dtype=np.int64) - for tid, idxs in tile_index_map_lists.items() - } - checkpoint["tile_index_map"] = { - tid.decode(): idxs.tolist() - for tid, idxs in tile_index_map.items() - } - with open(params["checkpoint"], "w") as cf: - json.dump(checkpoint, cf) - print(f" Map saved: {len(tile_index_map)} unique tiles.") + # Build TILE_ID → output row indices (always scanned; too large for checkpoint) + print(" Building tile→index map (scans all rows)...") + tile_index_map_lists = defaultdict(list) + with tqdm.tqdm( + total=n_total, unit="rows", unit_scale=True, desc=" Scanning TILE_ID" + ) as pbar: + for start in range(0, n_total, SCAN_CHUNK): + end = min(start + SCAN_CHUNK, n_total) + tile_chunk = dset[start:end]["TILE_ID"] + for local_i, tid in enumerate(tile_chunk): + tile_index_map_lists[tid].append(start + local_i) + pbar.update(end - start) + + tile_index_map = { + tid: np.array(idxs, dtype=np.int64) + for tid, idxs in tile_index_map_lists.items() + } + print(f" Map built: {len(tile_index_map)} unique tiles.") unique_tiles = sorted(tile_index_map.keys()) n_tiles = len(unique_tiles) From 7e88b056432f1a4a88d1785cd2ffe8e3bb3e8317 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Thu, 25 Jun 2026 15:30:28 +0200 Subject: [PATCH 05/13] adding mag errors to fill_photoz --- scripts/fill_photoz_bands.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/scripts/fill_photoz_bands.py b/scripts/fill_photoz_bands.py index 31ab290e..627ac3cf 100755 --- a/scripts/fill_photoz_bands.py +++ b/scripts/fill_photoz_bands.py @@ -45,12 +45,18 @@ "Z_B_MIN", "Z_B_MAX", "T_B", - "MAG_GAAP_0p7_u", - "MAG_GAAP_0p7_g", - "MAG_GAAP_0p7_r", - "MAG_GAAP_0p7_i", - "MAG_GAAP_0p7_z", - "MAG_GAAP_0p7_z2", + "MAG_GAAP_u", + "MAGERR_GAAP_u", + "MAG_GAAP_g", + "MAGERR_GAAP_g", + "MAG_GAAP_r", + "MAGERR_GAAP_r", + "MAG_GAAP_i", + "MAGERR_GAAP_i", + "MAG_GAAP_z", + "MAGERR_GAAP_z", + "MAG_GAAP_z2", + "MAGERR_GAAP_z2", ] From 310598d17b7b17756e28b227fd29a10bb6a51d33 Mon Sep 17 00:00:00 2001 From: "martin.kilbinger" Date: Thu, 25 Jun 2026 15:55:14 +0200 Subject: [PATCH 06/13] added Z_ML to fill_photoz --- scripts/fill_photoz_bands.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/fill_photoz_bands.py b/scripts/fill_photoz_bands.py index 627ac3cf..e69abd17 100755 --- a/scripts/fill_photoz_bands.py +++ b/scripts/fill_photoz_bands.py @@ -45,6 +45,7 @@ "Z_B_MIN", "Z_B_MAX", "T_B", + "Z_ML", "MAG_GAAP_u", "MAGERR_GAAP_u", "MAG_GAAP_g", From 9798047c676ce0cadfd48774fc301c5fe141c34d Mon Sep 17 00:00:00 2001 From: "martin.kilbinger@cea.fr" Date: Fri, 26 Jun 2026 10:53:20 +0200 Subject: [PATCH 07/13] added more flags to fill_photoz --- scripts/fill_photoz_bands.py | 38 ++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/scripts/fill_photoz_bands.py b/scripts/fill_photoz_bands.py index 627ac3cf..d11af42d 100755 --- a/scripts/fill_photoz_bands.py +++ b/scripts/fill_photoz_bands.py @@ -47,16 +47,54 @@ "T_B", "MAG_GAAP_u", "MAGERR_GAAP_u", + "FLAG_GAAP_u", + "MAG_LIM_u", + "FLUX_GAAP_u", + "FLUXERR_GAAP_u", + "EXTINCTION_u", "MAG_GAAP_g", "MAGERR_GAAP_g", + "FLAG_GAAP_g", + "MAG_LIM_g", + "FLUX_GAAP_g", + "FLUXERR_GAAP_g", + "EXTINCTION_g", "MAG_GAAP_r", "MAGERR_GAAP_r", + "FLAG_GAAP_r", + "MAG_LIM_r", + "FLUX_GAAP_r", + "FLUXERR_GAAP_r", + "EXTINCTION_r", "MAG_GAAP_i", "MAGERR_GAAP_i", + "FLAG_GAAP_i", + "MAG_LIM_i", + "FLUX_GAAP_i", + "FLUXERR_GAAP_i", + "EXTINCTION_i", "MAG_GAAP_z", "MAGERR_GAAP_z", + "FLAG_GAAP_z", + "MAG_LIM_z", + "FLUX_GAAP_z", + "FLUXERR_GAAP_z", + "EXTINCTION_z", "MAG_GAAP_z2", "MAGERR_GAAP_z2", + "FLAG_GAAP_z2", + "MAG_LIM_z2", + "FLUX_GAAP_z2", + "FLUXERR_GAAP_z2", + "EXTINCTION_z2", + "EXTINCTION", + "ODDS", + "CHI_SQUARED_BPZ", + "M_0", + "BPZ_FILT", + "BPZ_NONDETFILT", + "BPZ_FLAGFILT", + "MP_NAME", ] From 7cdcfdfa27555bb5d51b5c36c66b2165f22fd9a0 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 29 Jun 2026 08:51:13 +0200 Subject: [PATCH 08/13] added 0p7 and 1p0 aperture magnitudes to PhotoPipe + SP output --- scripts/fill_photoz_bands.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/scripts/fill_photoz_bands.py b/scripts/fill_photoz_bands.py index d11af42d..761b0308 100755 --- a/scripts/fill_photoz_bands.py +++ b/scripts/fill_photoz_bands.py @@ -45,48 +45,79 @@ "Z_B_MIN", "Z_B_MAX", "T_B", + "MAG_GAAP_u", "MAGERR_GAAP_u", + "MAG_GAAP_0p7_u", + "MAGERR_GAAP_0p7_u", + "MAG_GAAP_1p0_u", + "MAGERR_GAAP_1p0_u", "FLAG_GAAP_u", "MAG_LIM_u", "FLUX_GAAP_u", "FLUXERR_GAAP_u", "EXTINCTION_u", + "MAG_GAAP_g", "MAGERR_GAAP_g", + "MAG_GAAP_0p7_g", + "MAGERR_GAAP_0p7_g", + "MAG_GAAP_1p0_g", + "MAGERR_GAAP_1p0_g", "FLAG_GAAP_g", "MAG_LIM_g", "FLUX_GAAP_g", "FLUXERR_GAAP_g", "EXTINCTION_g", + "MAG_GAAP_r", "MAGERR_GAAP_r", + "MAG_GAAP_0p7_r", + "MAGERR_GAAP_0p7_r", + "MAG_GAAP_1p0_r", + "MAGERR_GAAP_1p0_r", "FLAG_GAAP_r", "MAG_LIM_r", "FLUX_GAAP_r", "FLUXERR_GAAP_r", "EXTINCTION_r", + "MAG_GAAP_i", "MAGERR_GAAP_i", + "MAG_GAAP_0p7_i", + "MAGERR_GAAP_0p7_i", + "MAG_GAAP_1p0_i", + "MAGERR_GAAP_1p0_i", "FLAG_GAAP_i", "MAG_LIM_i", "FLUX_GAAP_i", "FLUXERR_GAAP_i", "EXTINCTION_i", + "MAG_GAAP_z", "MAGERR_GAAP_z", + "MAG_GAAP_0p7_z", + "MAGERR_GAAP_0p7_z", + "MAG_GAAP_1p0_z", + "MAGERR_GAAP_1p0_z", "FLAG_GAAP_z", "MAG_LIM_z", "FLUX_GAAP_z", "FLUXERR_GAAP_z", "EXTINCTION_z", + "MAG_GAAP_z2", "MAGERR_GAAP_z2", + "MAG_GAAP_0p7_z2", + "MAGERR_GAAP_0p7_z2", + "MAG_GAAP_1p0_z2", + "MAGERR_GAAP_1p0_z2", "FLAG_GAAP_z2", "MAG_LIM_z2", "FLUX_GAAP_z2", "FLUXERR_GAAP_z2", "EXTINCTION_z2", + "EXTINCTION", "ODDS", "CHI_SQUARED_BPZ", From bfae66522c6d62ebe937b02ee3fc68c554c60627 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 6 Jul 2026 14:41:06 +0200 Subject: [PATCH 09/13] Fixed fill photoz script with new MP_NAME type --- notebooks/calibrate_comprehensive_cat.py | 4 +- scripts/check_filled_fields.py | 72 ++++++++++++++++++++++++ scripts/fill_photoz_bands.py | 34 ++++++++--- 3 files changed, 99 insertions(+), 11 deletions(-) create mode 100644 scripts/check_filled_fields.py diff --git a/notebooks/calibrate_comprehensive_cat.py b/notebooks/calibrate_comprehensive_cat.py index 404a7bec..5a0de543 100644 --- a/notebooks/calibrate_comprehensive_cat.py +++ b/notebooks/calibrate_comprehensive_cat.py @@ -47,8 +47,8 @@ dat, dat_ext = obj.read_cat(load_into_memory=False) # %% -n_test = -1 -#n_test = 100000 +#n_test = -1 +n_test = 100000 if n_test > 0: print(f"MKDEBUG testing only first {n_test} objects") dat = dat[:n_test] diff --git a/scripts/check_filled_fields.py b/scripts/check_filled_fields.py new file mode 100644 index 00000000..c998c26f --- /dev/null +++ b/scripts/check_filled_fields.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 +"""Check how many entries in selected HDF5 fields are filled (value != -199). + +Reads in chunks with a progress bar and periodic fill-fraction updates. +""" + +import h5py +import numpy as np +from tqdm import tqdm + +HDF5_FILE = "unions_shapepipe_comprehensive_struc_ugriz_2024_v1.6.c.DR6.hdf5" +EMPTY_VALUE = -199 +CHUNK_SIZE = 5_000_000 # rows per chunk +REPORT_EVERY = 10 # print running fractions every N chunks + +FIELDS = [ + "Z_B", + "Z_B_MIN", + "Z_B_MAX", + "T_B", + "MAG_GAAP_0p7_u", + "MAG_GAAP_1p0_u", + "MAG_GAAP_0p7_g", + "MAG_GAAP_1p0_g", + "MAG_GAAP_0p7_r", + "MAG_GAAP_1p0_r", + "MAG_GAAP_0p7_i", + "MAG_GAAP_1p0_i", + "MAG_GAAP_0p7_z", + "MAG_GAAP_1p0_z", + "MAG_GAAP_0p7_z2", + "MAG_GAAP_1p0_z2", +] + +with h5py.File(HDF5_FILE, "r") as f: + data = f["data"] + n_total = data.shape[0] + n_chunks = (n_total + CHUNK_SIZE - 1) // CHUNK_SIZE + print(f"Total entries : {n_total:,}") + print(f"Chunk size : {CHUNK_SIZE:,} ({n_chunks} chunks)\n") + + counts = {field: 0 for field in FIELDS} + + with tqdm(total=n_total, unit="rows", unit_scale=True, desc="Reading") as pbar: + for chunk_idx in range(n_chunks): + start = chunk_idx * CHUNK_SIZE + end = min(start + CHUNK_SIZE, n_total) + + for field in FIELDS: + counts[field] += int(np.sum(data[field, start:end] != EMPTY_VALUE)) + + pbar.update(end - start) + + # Periodic running-fraction report + if (chunk_idx + 1) % REPORT_EVERY == 0 or (chunk_idx + 1) == n_chunks: + rows_done = end + tqdm.write(f"\n --- after {rows_done:,} rows ({100*rows_done/n_total:.1f}%) ---") + tqdm.write(f" {'Field':<22} {'Filled %':>9}") + for field in FIELDS: + pct = 100.0 * counts[field] / rows_done + tqdm.write(f" {field:<22} {pct:>8.2f}%") + +# Final summary +print(f"\n{'='*56}") +print(f"FINAL SUMMARY (total rows: {n_total:,})") +print(f"{'Field':<22} {'Filled':>12} {'Empty':>12} {'Filled %':>10}") +print("-" * 60) +for field in FIELDS: + n_filled = counts[field] + n_empty = n_total - n_filled + pct = 100.0 * n_filled / n_total + print(f"{field:<22} {n_filled:>12,} {n_empty:>12,} {pct:>9.2f}%") diff --git a/scripts/fill_photoz_bands.py b/scripts/fill_photoz_bands.py index 761b0308..f29c9e4a 100755 --- a/scripts/fill_photoz_bands.py +++ b/scripts/fill_photoz_bands.py @@ -39,13 +39,13 @@ EMPTY_VALUE = -199 COPY_CHUNK = 2_000_000 # rows per chunk when copying input → output SCAN_CHUNK = 5_000_000 # rows per chunk when scanning TILE_ID +MAX_CONSEC_FAILS = 10 # abort if this many tiles in a row fail REQUESTED_KEYS = [ "Z_B", "Z_B_MIN", "Z_B_MAX", "T_B", - "MAG_GAAP_u", "MAGERR_GAAP_u", "MAG_GAAP_0p7_u", @@ -57,7 +57,6 @@ "FLUX_GAAP_u", "FLUXERR_GAAP_u", "EXTINCTION_u", - "MAG_GAAP_g", "MAGERR_GAAP_g", "MAG_GAAP_0p7_g", @@ -69,7 +68,6 @@ "FLUX_GAAP_g", "FLUXERR_GAAP_g", "EXTINCTION_g", - "MAG_GAAP_r", "MAGERR_GAAP_r", "MAG_GAAP_0p7_r", @@ -81,7 +79,6 @@ "FLUX_GAAP_r", "FLUXERR_GAAP_r", "EXTINCTION_r", - "MAG_GAAP_i", "MAGERR_GAAP_i", "MAG_GAAP_0p7_i", @@ -93,7 +90,6 @@ "FLUX_GAAP_i", "FLUXERR_GAAP_i", "EXTINCTION_i", - "MAG_GAAP_z", "MAGERR_GAAP_z", "MAG_GAAP_0p7_z", @@ -105,7 +101,6 @@ "FLUX_GAAP_z", "FLUXERR_GAAP_z", "EXTINCTION_z", - "MAG_GAAP_z2", "MAGERR_GAAP_z2", "MAG_GAAP_0p7_z2", @@ -117,7 +112,6 @@ "FLUX_GAAP_z2", "FLUXERR_GAAP_z2", "EXTINCTION_z2", - "EXTINCTION", "ODDS", "CHI_SQUARED_BPZ", @@ -125,7 +119,6 @@ "BPZ_FILT", "BPZ_NONDETFILT", "BPZ_FLAGFILT", - "MP_NAME", ] @@ -144,7 +137,7 @@ def params_default(): params = { "input": "unions_shapepipe_comprehensive_struc_2024_v1.5.c.hdf5", "output": "unions_shapepipe_comprehensive_struc_ugriz_2024_v1.5.c.hdf5", - "fits_dir": "UNIONS5000", + "fits_dir": "UNIONS_DR6", "checkpoint": "fill_photoz_bands_checkpoint.json", "verbose": False, } @@ -490,12 +483,25 @@ def main(argv=None): n_skipped_missing = 0 n_skipped_done = 0 n_skipped_size = 0 + n_errors = 0 + n_consec_fails = 0 n_processed = 0 print(f"\n Processing {n_tiles} tiles ({len(done_tiles)} already done)...\n") pbar = tqdm.tqdm(unique_tiles, total=n_tiles, unit="tile") for tile_id in pbar: + if n_consec_fails >= MAX_CONSEC_FAILS: + checkpoint["done_tiles"] = list(done_tiles) + with open(params["checkpoint"], "w") as cf: + json.dump(checkpoint, cf) + print( + f"\nERROR: {n_consec_fails} tiles failed in a row, " + "likely a systematic problem; aborting.", + file=sys.stderr, + ) + sys.exit(1) + tile_str = tile_id.decode() if isinstance(tile_id, bytes) else tile_id if tile_str in done_tiles: @@ -535,6 +541,7 @@ def main(argv=None): f"{n_fits} — size mismatch, skipping." ) n_skipped_size += 1 + n_consec_fails += 1 pbar.set_postfix( {"done": n_processed, "size_err": n_skipped_size} ) @@ -544,9 +551,12 @@ def main(argv=None): except Exception as e: warnings.warn(f"Tile {tile_str}: error ({e}), skipping.") + n_errors += 1 + n_consec_fails += 1 continue n_processed += 1 + n_consec_fails = 0 done_tiles.add(tile_str) if n_processed % 50 == 0: @@ -567,6 +577,12 @@ def main(argv=None): print(f" Tiles skipped (done) : {n_skipped_done}") print(f" FITS files missing : {n_skipped_missing}") print(f" Size mismatches : {n_skipped_size}") + print(f" Tiles failed (error) : {n_errors}") + if n_processed == 0 and (n_errors > 0 or n_skipped_size > 0): + print( + "WARNING: no tiles were filled; all available tiles failed.", + file=sys.stderr, + ) return 0 From 2d5ef0745f37b09130fe83946b40c5c5f65f4170 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Fri, 28 Aug 2026 05:39:27 +0000 Subject: [PATCH 10/13] ruff autofix (format + safe lint fixes) Pushed by the lint gate. --- check.txt | 4 +++- format.txt | 5 +++-- report.md | 11 ++++++++--- scripts/check_filled_fields.py | 10 ++++++---- scripts/fill_photoz_bands.py | 13 +++++++------ 5 files changed, 27 insertions(+), 16 deletions(-) diff --git a/check.txt b/check.txt index 1f5f344d..6589f995 100644 --- a/check.txt +++ b/check.txt @@ -1 +1,3 @@ -All checks passed! +scripts/fill_photoz_bands.py:22:1: I001 [*] Import block is un-sorted or un-formatted +Found 1 error. +[*] 1 fixable with the `--fix` option. diff --git a/format.txt b/format.txt index a0b12d74..9e4589eb 100644 --- a/format.txt +++ b/format.txt @@ -1,2 +1,3 @@ -Would reformat: src/sp_validation/tests/test_masks.py -1 file would be reformatted, 225 files already formatted +Would reformat: scripts/check_filled_fields.py +Would reformat: scripts/fill_photoz_bands.py +2 files would be reformatted, 231 files already formatted diff --git a/report.md b/report.md index e681a67c..de9975ce 100644 --- a/report.md +++ b/report.md @@ -1,10 +1,15 @@ ### `ruff check .` -✅ clean +``` +scripts/fill_photoz_bands.py:22:1: I001 [*] Import block is un-sorted or un-formatted +Found 1 error. +[*] 1 fixable with the `--fix` option. +``` ### `ruff format --check .` ``` -Would reformat: src/sp_validation/tests/test_masks.py -1 file would be reformatted, 225 files already formatted +Would reformat: scripts/check_filled_fields.py +Would reformat: scripts/fill_photoz_bands.py +2 files would be reformatted, 231 files already formatted ``` diff --git a/scripts/check_filled_fields.py b/scripts/check_filled_fields.py index c998c26f..675ad116 100644 --- a/scripts/check_filled_fields.py +++ b/scripts/check_filled_fields.py @@ -10,8 +10,8 @@ HDF5_FILE = "unions_shapepipe_comprehensive_struc_ugriz_2024_v1.6.c.DR6.hdf5" EMPTY_VALUE = -199 -CHUNK_SIZE = 5_000_000 # rows per chunk -REPORT_EVERY = 10 # print running fractions every N chunks +CHUNK_SIZE = 5_000_000 # rows per chunk +REPORT_EVERY = 10 # print running fractions every N chunks FIELDS = [ "Z_B", @@ -54,14 +54,16 @@ # Periodic running-fraction report if (chunk_idx + 1) % REPORT_EVERY == 0 or (chunk_idx + 1) == n_chunks: rows_done = end - tqdm.write(f"\n --- after {rows_done:,} rows ({100*rows_done/n_total:.1f}%) ---") + tqdm.write( + f"\n --- after {rows_done:,} rows ({100 * rows_done / n_total:.1f}%) ---" + ) tqdm.write(f" {'Field':<22} {'Filled %':>9}") for field in FIELDS: pct = 100.0 * counts[field] / rows_done tqdm.write(f" {field:<22} {pct:>8.2f}%") # Final summary -print(f"\n{'='*56}") +print(f"\n{'=' * 56}") print(f"FINAL SUMMARY (total rows: {n_total:,})") print(f"{'Field':<22} {'Filled':>12} {'Empty':>12} {'Filled %':>10}") print("-" * 60) diff --git a/scripts/fill_photoz_bands.py b/scripts/fill_photoz_bands.py index f29c9e4a..442012a9 100755 --- a/scripts/fill_photoz_bands.py +++ b/scripts/fill_photoz_bands.py @@ -30,16 +30,14 @@ import numpy as np import tqdm from astropy.io import fits - from cs_util import args as cs_args from cs_util import logging - FITS_HDU = 1 EMPTY_VALUE = -199 -COPY_CHUNK = 2_000_000 # rows per chunk when copying input → output -SCAN_CHUNK = 5_000_000 # rows per chunk when scanning TILE_ID -MAX_CONSEC_FAILS = 10 # abort if this many tiles in a row fail +COPY_CHUNK = 2_000_000 # rows per chunk when copying input → output +SCAN_CHUNK = 5_000_000 # rows per chunk when scanning TILE_ID +MAX_CONSEC_FAILS = 10 # abort if this many tiles in a row fail REQUESTED_KEYS = [ "Z_B", @@ -165,6 +163,7 @@ def params_default(): # Helpers # --------------------------------------------------------------------------- + def detect_dataset_name(hf): """Return the first dataset name in an HDF5 file. @@ -305,6 +304,7 @@ def write_tile_to_hdf5(dset, hdf5_indices, fits_data, valid_keys): # Phase 1: create output file # --------------------------------------------------------------------------- + def create_output_file(input_path, output_path, dataset_name, verbose=False): """Create output HDF5 by copying input and appending empty PhotoPipe fields. @@ -370,6 +370,7 @@ def create_output_file(input_path, output_path, dataset_name, verbose=False): # Main # --------------------------------------------------------------------------- + def main(argv=None): """Main. @@ -479,7 +480,7 @@ def main(argv=None): unique_tiles = sorted(tile_index_map.keys()) n_tiles = len(unique_tiles) - valid_keys = None # determined from first available FITS tile + valid_keys = None # determined from first available FITS tile n_skipped_missing = 0 n_skipped_done = 0 n_skipped_size = 0 From c43688d01a1c57aa3eaba9a73bf1ddb7d0c9e4c6 Mon Sep 17 00:00:00 2001 From: "martin.kilbinger" Date: Mon, 31 Aug 2026 10:32:00 +0200 Subject: [PATCH 11/13] removed leftover git marker --- cosmo_val/cat_config.yaml | 1 - 1 file changed, 1 deletion(-) diff --git a/cosmo_val/cat_config.yaml b/cosmo_val/cat_config.yaml index b8d1a50a..5837b20f 100644 --- a/cosmo_val/cat_config.yaml +++ b/cosmo_val/cat_config.yaml @@ -1271,7 +1271,6 @@ SP_v1.4.8_uncal: path: unions_shapepipe_psf_2024_v1.4.a.fits hdu: 1 patch_number: 100 ->>>>>>> upstream/develop SP_v1.4_LFmask_8k: subdir: /n17data/mkilbing/astro/data/CFIS/v1.0/SP_LFmask pipeline: SP From 4577a37f76dca30ba26a4d504af33aa3f1db61ef Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 4 Sep 2026 23:47:04 +0200 Subject: [PATCH 12/13] fill_photoz_bands: spot-check FITS/HDF5 row order before filling a tile The fill pairs FITS row k with the k-th HDF5 row of the tile (sorted-index order) and only ever verified the row *count*, so a PhotoPipe tile ordered differently from the comprehensive catalogue would be filled with silently mismatched photo-z. After the size check, compare RA/Dec for a small sample of rows -- up to five at each end plus evenly spaced interior rows, --n_check_rows (default 10), --check_tol_arcsec (default 0.5). Only the sampled HDF5 rows are read (one fancy-index into the already-sorted index array), so the cost is a handful of point reads per tile rather than a full per-row match. A failing tile is warned about, counted as a row-order mismatch in the end-of-run summary, counted towards the consecutive-failure abort, and skipped without being added to done_tiles -- exactly as a size mismatch is, so a resume retries it. --n_check_rows 0 disables the check. HDF5 columns are RA/Dec (cat_config.yaml ra_col/dec_col for SP_v1.4.x). The FITS names are resolved at runtime from a candidate list; ALPHA_J2000 / DELTA_J2000 is first, verified against a real DR6 tile (/n17data/UNIONS/WL/photometry/UNIONS_DR6/UNIONS.001.227_SP_ugriz_photoz_ext.cat). If no candidate pair matches, the check disables itself with a warning rather than skipping tiles. Test: synthetic 3-tile HDF5 + FITS pair, tiles interleaved so the non-contiguous write path runs too; the two aligned tiles fill, the tile with reversed FITS rows is skipped and left empty, and --n_check_rows 0 reproduces the old unchecked behaviour. --- scripts/fill_photoz_bands.py | 260 +++++++++++++++++- .../tests/test_fill_photoz_bands.py | 237 ++++++++++++++++ 2 files changed, 494 insertions(+), 3 deletions(-) create mode 100644 src/sp_validation/tests/test_fill_photoz_bands.py diff --git a/scripts/fill_photoz_bands.py b/scripts/fill_photoz_bands.py index 88c9e16a..3eaf32dd 100755 --- a/scripts/fill_photoz_bands.py +++ b/scripts/fill_photoz_bands.py @@ -12,6 +12,11 @@ Phase 2 — Fill tiles: for each PhotoPipe FITS tile found in fits_dir, write the field values into the output file. +The fill assumes that, within a tile, the PhotoPipe FITS rows are in the +same order as the HDF5 rows of that tile. Before writing, this is checked +on a small sample of rows by comparing their sky positions (see +``check_tile_row_order``); tiles that fail are skipped, not written. + Skips missing FITS files. Warns (does not error) on missing PhotoPipe keys. Supports interrupt + restart at any point. @@ -39,6 +44,21 @@ SCAN_CHUNK = 5_000_000 # rows per chunk when scanning TILE_ID MAX_CONSEC_FAILS = 10 # abort if this many tiles in a row fail +# Sky-position columns of the ShapePipe comprehensive HDF5 catalogue +# (cosmo_val/cat_config.yaml: ra_col: RA, dec_col: Dec for all SP_v1.4.x). +RA_COL_HDF5 = "RA" +DEC_COL_HDF5 = "Dec" + +# Candidate (RA, Dec) column names in the PhotoPipe FITS tiles, most likely +# first. DR6 tiles (UNIONS.*_SP_ugriz_photoz_ext.cat) carry ALPHA_J2000 / +# DELTA_J2000; the remaining pairs are fallbacks for other PhotoPipe outputs. +FITS_RADEC_CANDIDATES = [ + ("ALPHA_J2000", "DELTA_J2000"), + ("RA", "DEC"), + ("RA", "Dec"), + ("X_WORLD", "Y_WORLD"), +] + REQUESTED_KEYS = [ "Z_B", "Z_B_MIN", @@ -138,6 +158,8 @@ def params_default(): "output": "unions_shapepipe_comprehensive_struc_ugriz_2024_v1.5.c.hdf5", "fits_dir": "UNIONS_DR6", "checkpoint": "fill_photoz_bands_checkpoint.json", + "n_check_rows": 10, + "check_tol_arcsec": 0.5, "verbose": False, } @@ -148,13 +170,24 @@ def params_default(): "checkpoint": "-c", } - types = {} + types = { + "n_check_rows": "int", + "check_tol_arcsec": "float", + } help_strings = { "input": "input HDF5 catalogue (no PhotoPipe fields), default={}", "output": "output HDF5 catalogue (PhotoPipe fields added and filled), default={}", "fits_dir": "directory with PhotoPipe FITS tiles, default={}", "checkpoint": "checkpoint JSON file for resume support, default={}", + "n_check_rows": ( + "number of rows per tile whose RA/Dec are compared between HDF5 and" + " FITS to verify row order, 0 to disable, default={}" + ), + "check_tol_arcsec": ( + "maximum angular separation [arcsec] for a row-order check to pass," + " default={}" + ), } return params, short_options, types, help_strings @@ -257,12 +290,173 @@ def check_fits_keys(fits_columns, requested_keys): return valid, missing +def find_fits_radec_columns(fits_columns): + """Find RA/Dec Columns. + + Return the (RA, Dec) column names of a PhotoPipe FITS table. + + Parameters + ---------- + fits_columns : iterable of str + Column names of the FITS HDU. + + Returns + ------- + tuple + (ra_name, dec_name), or (None, None) if no candidate pair matches. + + """ + columns = set(fits_columns) + for ra_name, dec_name in FITS_RADEC_CANDIDATES: + if ra_name in columns and dec_name in columns: + return ra_name, dec_name + return None, None + + +def sample_row_positions(n_rows, n_sample): + """Sample Row Positions. + + Pick a small, cheap-to-read set of row positions spanning a tile: up to + five rows at each end (where a shifted or reversed ordering shows up + first) plus evenly spaced rows in between. + + Parameters + ---------- + n_rows : int + Number of rows in the tile. + n_sample : int + Requested number of sampled rows. + + Returns + ------- + numpy.ndarray + Sorted, unique row positions, of length ``min(n_sample, n_rows)`` + or less. + + """ + if n_rows <= 0 or n_sample <= 0: + return np.empty(0, dtype=np.int64) + if n_sample >= n_rows: + return np.arange(n_rows, dtype=np.int64) + + n_edge = min(5, n_sample // 3) + head = np.arange(n_edge, dtype=np.int64) + tail = np.arange(n_rows - n_edge, n_rows, dtype=np.int64) + + n_mid = n_sample - 2 * n_edge + if n_mid > 0: + mid = np.linspace(n_edge, n_rows - n_edge - 1, n_mid).astype(np.int64) + else: + mid = np.empty(0, dtype=np.int64) + + return np.unique(np.concatenate([head, mid, tail])) + + +def angular_separation_arcsec(ra_1, dec_1, ra_2, dec_2): + """Angular Separation Arcsec. + + Great-circle separation between two sets of sky coordinates, computed + with the haversine formula (numerically stable at small separations and + correct across the RA=0 wrap). + + Parameters + ---------- + ra_1 : numpy.ndarray + Right ascensions of the first set [deg]. + dec_1 : numpy.ndarray + Declinations of the first set [deg]. + ra_2 : numpy.ndarray + Right ascensions of the second set [deg]. + dec_2 : numpy.ndarray + Declinations of the second set [deg]. + + Returns + ------- + numpy.ndarray + Angular separations [arcsec]. + + """ + ra_1 = np.radians(np.asarray(ra_1, dtype=np.float64)) + dec_1 = np.radians(np.asarray(dec_1, dtype=np.float64)) + ra_2 = np.radians(np.asarray(ra_2, dtype=np.float64)) + dec_2 = np.radians(np.asarray(dec_2, dtype=np.float64)) + + d_ra = ra_2 - ra_1 + d_dec = dec_2 - dec_1 + hav = np.sin(d_dec / 2) ** 2 + np.cos(dec_1) * np.cos(dec_2) * np.sin(d_ra / 2) ** 2 + sep_rad = 2 * np.arcsin(np.sqrt(np.clip(hav, 0, 1))) + + return np.degrees(sep_rad) * 3600 + + +def check_tile_row_order( + dset, + sorted_idx, + fits_data, + fits_radec_cols, + n_sample, + tol_arcsec, +): + """Check Tile Row Order. + + Spot-check that the FITS rows of a tile line up with the HDF5 rows they + are about to be written into. ``write_tile_to_hdf5`` assigns FITS row + ``k`` to HDF5 row ``sorted_idx[k]``; this compares the sky positions of + a sample of those pairs. Only the sampled HDF5 rows are read. + + Parameters + ---------- + dset : h5py.Dataset + Compound HDF5 dataset. + sorted_idx : numpy.ndarray + Sorted HDF5 row indices of this tile. + fits_data : numpy.recarray + Data from the PhotoPipe FITS HDU, same length as ``sorted_idx``. + fits_radec_cols : tuple + (RA, Dec) column names in ``fits_data``. + n_sample : int + Number of rows to compare. + tol_arcsec : float + Maximum tolerated angular separation [arcsec]. + + Returns + ------- + tuple + (ok, n_checked, max_sep_arcsec). ``max_sep_arcsec`` is ``numpy.nan`` + if no row could be checked. + + """ + positions = sample_row_positions(len(sorted_idx), n_sample) + if len(positions) == 0: + return True, 0, np.nan + + ra_col_fits, dec_col_fits = fits_radec_cols + + # Fancy-index a handful of rows only; positions is sorted and unique, and + # sorted_idx is sorted, so the h5py selection is strictly increasing. + rows_hdf5 = dset[sorted_idx[positions]] + + sep = angular_separation_arcsec( + rows_hdf5[RA_COL_HDF5], + rows_hdf5[DEC_COL_HDF5], + fits_data[ra_col_fits][positions], + fits_data[dec_col_fits][positions], + ) + max_sep = float(np.nanmax(sep)) + + return bool(np.all(sep <= tol_arcsec)), len(positions), max_sep + + def write_tile_to_hdf5(dset, hdf5_indices, fits_data, valid_keys): """Write valid_keys from fits_data into dset at hdf5_indices. Reads the HDF5 range in one chunk, fills fields in memory, writes back. Handles both contiguous and non-contiguous index ranges. + Assumes row ``k`` of ``fits_data`` corresponds to HDF5 row + ``numpy.sort(hdf5_indices)[k]``; ``check_tile_row_order`` spot-checks + this before the write. + Parameters ---------- dset : h5py.Dataset @@ -482,9 +676,20 @@ def main(argv=None): n_tiles = len(unique_tiles) valid_keys = None # determined from first available FITS tile + fits_radec_cols = None # idem, (RA, Dec) column names in the FITS tiles + + # The row-order check needs sky positions on both sides. + check_rows = params["n_check_rows"] + if check_rows > 0 and not {RA_COL_HDF5, DEC_COL_HDF5} <= hdf5_fields: + warnings.warn( + f"Columns '{RA_COL_HDF5}'/'{DEC_COL_HDF5}' absent from the HDF5" + " dataset — row-order check disabled." + ) + check_rows = 0 n_skipped_missing = 0 n_skipped_done = 0 n_skipped_size = 0 + n_skipped_order = 0 n_errors = 0 n_consec_fails = 0 n_processed = 0 @@ -533,6 +738,24 @@ def main(argv=None): ) tqdm.tqdm.write(f"\n Keys to fill: {valid_keys}\n") + fits_radec_cols = find_fits_radec_columns(fits_data.dtype.names) + if check_rows > 0 and fits_radec_cols[0] is None: + warnings.warn( + "No known RA/Dec column pair in the PhotoPipe" + f" FITS (tried {FITS_RADEC_CANDIDATES}) —" + " row-order check disabled." + ) + check_rows = 0 + elif check_rows > 0: + tqdm.tqdm.write( + " Row-order check: comparing" + f" {check_rows} rows/tile," + f" {RA_COL_HDF5}/{DEC_COL_HDF5} (HDF5) vs" + f" {fits_radec_cols[0]}/{fits_radec_cols[1]}" + f" (FITS), tolerance" + f" {params['check_tol_arcsec']} arcsec\n" + ) + hdf5_indices = tile_index_map[tile_id] n_hdf5 = len(hdf5_indices) n_fits = len(fits_data) @@ -549,7 +772,37 @@ def main(argv=None): ) continue - write_tile_to_hdf5(dset, hdf5_indices, fits_data, valid_keys) + # write_tile_to_hdf5 pairs FITS row k with HDF5 row + # sorted_idx[k]; spot-check that pairing before writing. + sorted_idx = np.sort(hdf5_indices) + if check_rows > 0: + ok, n_checked, max_sep = check_tile_row_order( + dset, + sorted_idx, + fits_data, + fits_radec_cols, + check_rows, + params["check_tol_arcsec"], + ) + if not ok: + warnings.warn( + f"Tile {tile_str}: RA/Dec disagree for" + f" {n_checked} checked rows (max separation" + f" {max_sep:.3g} arcsec >" + f" {params['check_tol_arcsec']} arcsec) —" + " row order mismatch, skipping." + ) + n_skipped_order += 1 + n_consec_fails += 1 + pbar.set_postfix( + { + "done": n_processed, + "order_err": n_skipped_order, + } + ) + continue + + write_tile_to_hdf5(dset, sorted_idx, fits_data, valid_keys) except Exception as e: warnings.warn(f"Tile {tile_str}: error ({e}), skipping.") @@ -579,8 +832,9 @@ def main(argv=None): print(f" Tiles skipped (done) : {n_skipped_done}") print(f" FITS files missing : {n_skipped_missing}") print(f" Size mismatches : {n_skipped_size}") + print(f" Row-order mismatches : {n_skipped_order}") print(f" Tiles failed (error) : {n_errors}") - if n_processed == 0 and (n_errors > 0 or n_skipped_size > 0): + if n_processed == 0 and (n_errors > 0 or n_skipped_size > 0 or n_skipped_order > 0): print( "WARNING: no tiles were filled; all available tiles failed.", file=sys.stderr, diff --git a/src/sp_validation/tests/test_fill_photoz_bands.py b/src/sp_validation/tests/test_fill_photoz_bands.py new file mode 100644 index 00000000..f45c1161 --- /dev/null +++ b/src/sp_validation/tests/test_fill_photoz_bands.py @@ -0,0 +1,237 @@ +"""``scripts/fill_photoz_bands.py`` must not trust FITS/HDF5 row order blindly. + +The script writes PhotoPipe columns into the ShapePipe comprehensive HDF5 +catalogue tile by tile, pairing FITS row ``k`` with the ``k``-th HDF5 row of +that tile (in sorted-index order). Nothing in the file formats guarantees +that pairing -- only the row *count* used to be checked -- so a tile whose +PhotoPipe catalogue happens to be ordered differently would be filled with +silently mismatched photo-z. + +These tests build a tiny synthetic pair of catalogues (three tiles, +deliberately interleaved so the non-contiguous write path is exercised) and +run the script end to end: + +* the two tiles whose FITS rows line up are filled; +* the tile whose FITS rows are reversed is *skipped*, counted under + ``Row-order mismatches``, and left at ``EMPTY_VALUE``, so a resumed run + retries it rather than treating it as done. +""" + +import importlib.util +import json +import subprocess +import sys +from pathlib import Path + +import numpy as np +import pytest + +h5py = pytest.importorskip("h5py") +fits = pytest.importorskip("astropy.io.fits") +pytest.importorskip("cs_util") +pytest.importorskip("tqdm") + + +def _repo_root() -> Path: + """Locate the repo root by walking up to the ``pyproject.toml`` marker.""" + for parent in Path(__file__).resolve().parents: + if (parent / "pyproject.toml").exists(): + return parent + raise RuntimeError("could not locate repo root (no pyproject.toml above test)") + + +_SCRIPT = _repo_root() / "scripts" / "fill_photoz_bands.py" + +# Tile layout: three tiles, rows interleaved round-robin over the catalogue so +# each tile's HDF5 indices are non-contiguous. +TILES = ("100.100", "200.200", "300.300") +N_PER_TILE = 20 +N_ROWS = len(TILES) * N_PER_TILE +FILL_KEYS = ("Z_B", "MAG_GAAP_r") + + +def _module(): + """Import ``scripts/fill_photoz_bands.py`` (lives outside the package).""" + spec = importlib.util.spec_from_file_location("fill_photoz_bands", _SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def _synthetic_catalogues(tmp_path, reversed_tile): + """Write a synthetic HDF5 catalogue and matching PhotoPipe FITS tiles. + + Parameters + ---------- + tmp_path : pathlib.Path + Directory to write into. + reversed_tile : str + Tile whose FITS rows are written in reverse order (row-order breakage). + + Returns + ------- + tuple + (input HDF5 path, FITS directory, truth dict tile -> (Z_B, MAG_GAAP_r)) + + """ + row = np.arange(N_ROWS) + tile_of_row = np.array([TILES[i % len(TILES)] for i in row]) + + dtype = np.dtype([("TILE_ID", "S7"), ("RA", " 0) + # Not all at the edges: something is sampled from the bulk. + assert np.any((positions > 10) & (positions < 989)) + + +@pytest.mark.fast +def test_find_fits_radec_columns(): + """The DR6 PhotoPipe names are recognised; unknown tables yield None.""" + module = _module() + + assert module.find_fits_radec_columns( + ["SeqNr", "ALPHA_J2000", "DELTA_J2000", "Z_B"] + ) == ("ALPHA_J2000", "DELTA_J2000") + assert module.find_fits_radec_columns(["RA", "DEC"]) == ("RA", "DEC") + assert module.find_fits_radec_columns(["X", "Y"]) == (None, None) + + +@pytest.mark.fast +def test_angular_separation_arcsec(): + """Separation is correct at the pole-free small-angle limit and wraps RA.""" + module = _module() + + sep = module.angular_separation_arcsec([0.0], [0.0], [1.0 / 3600], [0.0]) + np.testing.assert_allclose(sep, [1.0], rtol=1e-6) + + # Across the RA=0 wrap: 359.999 deg vs 0.001 deg is 0.002 deg, not 360. + sep = module.angular_separation_arcsec([359.999], [0.0], [0.001], [0.0]) + np.testing.assert_allclose(sep, [0.002 * 3600], rtol=1e-6) From 5e05634a95d3bb4d9a21192c9441880203997f16 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 5 Sep 2026 11:25:04 +0200 Subject: [PATCH 13/13] fill_photoz_bands: row-order check ignores invalid positions, fails if none A sampled row with a non-finite or sentinel (|Dec| > 90) coordinate on either side carries no information about the FITS/HDF5 pairing. Before, such a row made `sep <= tol` False (a false row-order failure), and an all-NaN sample raised a RuntimeWarning from np.nanmax. Now those rows are excluded from the comparison, n_checked reports only the rows actually compared, and a tile whose sample has no comparable row is treated as unverifiable and skipped (distinct warning), so it is retried on resume rather than written blind. Adds a direct unit test covering partial/all-invalid samples, a reversed remainder, and a single-row tile (one-element h5py fancy index). Co-Authored-By: Claude Fable 5.1 --- scripts/fill_photoz_bands.py | 59 +++++++++++----- .../tests/test_fill_photoz_bands.py | 67 +++++++++++++++++++ 2 files changed, 110 insertions(+), 16 deletions(-) diff --git a/scripts/fill_photoz_bands.py b/scripts/fill_photoz_bands.py index 3eaf32dd..36a0a1b2 100755 --- a/scripts/fill_photoz_bands.py +++ b/scripts/fill_photoz_bands.py @@ -422,13 +422,17 @@ def check_tile_row_order( Returns ------- tuple - (ok, n_checked, max_sep_arcsec). ``max_sep_arcsec`` is ``numpy.nan`` - if no row could be checked. + (ok, n_checked, max_sep_arcsec). ``n_checked`` counts the sampled + rows with a valid sky position on both sides; only those are + compared. ``ok`` is ``False`` if any compared pair is further apart + than ``tol_arcsec``, or if the sample holds no comparable row at all + (``n_checked == 0``, ``max_sep_arcsec`` is ``numpy.nan``): a tile + whose order cannot be verified is not written. """ positions = sample_row_positions(len(sorted_idx), n_sample) if len(positions) == 0: - return True, 0, np.nan + return False, 0, np.nan ra_col_fits, dec_col_fits = fits_radec_cols @@ -436,15 +440,31 @@ def check_tile_row_order( # sorted_idx is sorted, so the h5py selection is strictly increasing. rows_hdf5 = dset[sorted_idx[positions]] + ra_hdf5 = np.asarray(rows_hdf5[RA_COL_HDF5], dtype=np.float64) + dec_hdf5 = np.asarray(rows_hdf5[DEC_COL_HDF5], dtype=np.float64) + ra_fits = np.asarray(fits_data[ra_col_fits][positions], dtype=np.float64) + dec_fits = np.asarray(fits_data[dec_col_fits][positions], dtype=np.float64) + + # A row carries no information about the pairing if either side has no + # valid position (NaN, or a sentinel such as -199 outside the sky). + valid = ( + np.isfinite(ra_hdf5) + & np.isfinite(dec_hdf5) + & np.isfinite(ra_fits) + & np.isfinite(dec_fits) + & (np.abs(dec_hdf5) <= 90) + & (np.abs(dec_fits) <= 90) + ) + n_checked = int(valid.sum()) + if n_checked == 0: + return False, 0, np.nan + sep = angular_separation_arcsec( - rows_hdf5[RA_COL_HDF5], - rows_hdf5[DEC_COL_HDF5], - fits_data[ra_col_fits][positions], - fits_data[dec_col_fits][positions], + ra_hdf5[valid], dec_hdf5[valid], ra_fits[valid], dec_fits[valid] ) - max_sep = float(np.nanmax(sep)) + max_sep = float(np.max(sep)) - return bool(np.all(sep <= tol_arcsec)), len(positions), max_sep + return bool(max_sep <= tol_arcsec), n_checked, max_sep def write_tile_to_hdf5(dset, hdf5_indices, fits_data, valid_keys): @@ -785,13 +805,20 @@ def main(argv=None): params["check_tol_arcsec"], ) if not ok: - warnings.warn( - f"Tile {tile_str}: RA/Dec disagree for" - f" {n_checked} checked rows (max separation" - f" {max_sep:.3g} arcsec >" - f" {params['check_tol_arcsec']} arcsec) —" - " row order mismatch, skipping." - ) + if n_checked == 0: + warnings.warn( + f"Tile {tile_str}: no sampled row has a" + " valid RA/Dec on both sides —" + " row order unverifiable, skipping." + ) + else: + warnings.warn( + f"Tile {tile_str}: RA/Dec disagree for" + f" {n_checked} checked rows (max" + f" separation {max_sep:.3g} arcsec >" + f" {params['check_tol_arcsec']} arcsec)" + " — row order mismatch, skipping." + ) n_skipped_order += 1 n_consec_fails += 1 pbar.set_postfix( diff --git a/src/sp_validation/tests/test_fill_photoz_bands.py b/src/sp_validation/tests/test_fill_photoz_bands.py index f45c1161..ae45a637 100644 --- a/src/sp_validation/tests/test_fill_photoz_bands.py +++ b/src/sp_validation/tests/test_fill_photoz_bands.py @@ -21,6 +21,7 @@ import json import subprocess import sys +import warnings from pathlib import Path import numpy as np @@ -235,3 +236,69 @@ def test_angular_separation_arcsec(): # Across the RA=0 wrap: 359.999 deg vs 0.001 deg is 0.002 deg, not 360. sep = module.angular_separation_arcsec([359.999], [0.0], [0.001], [0.0]) np.testing.assert_allclose(sep, [0.002 * 3600], rtol=1e-6) + + +@pytest.mark.fast +def test_check_tile_row_order_ignores_invalid_positions(tmp_path): + """Rows without a valid position on both sides are not compared. + + A NaN or sentinel coordinate carries no information about the pairing, + so it must neither fail the check (NaN <= tol is False) nor pass it; a + sample with no comparable row at all is an unverifiable tile and fails. + """ + module = _module() + + n = 8 + ra = 10.0 + 0.01 * np.arange(n) + dec = 20.0 + 0.01 * np.arange(n) + dtype = np.dtype([("RA", ">f8"), ("Dec", ">f8")]) + data = np.empty(n, dtype=dtype) + data["RA"], data["Dec"] = ra, dec + with h5py.File(tmp_path / "cat.hdf5", "w") as hf: + dset = hf.create_dataset("dat", data=data) + sorted_idx = np.arange(n) + cols = ("ALPHA_J2000", "DELTA_J2000") + + def fits_like(ra_f, dec_f): + rec = np.empty(n, dtype=[(cols[0], " 0.5 + + # Nothing comparable at all: fail, no RuntimeWarning from an all-NaN max. + with warnings.catch_warnings(): + warnings.simplefilter("error") + ok, n_checked, max_sep = module.check_tile_row_order( + dset, sorted_idx, fits_like(np.full(n, np.nan), dec), cols, n, 0.5 + ) + assert not ok and n_checked == 0 and np.isnan(max_sep) + + # A single-row tile is a valid (one-element) fancy index. + ok, n_checked, _ = module.check_tile_row_order( + dset, sorted_idx[3:4], fits_like(ra, dec)[3:4], cols, 10, 0.5 + ) + assert ok and n_checked == 1