From c90ae4ff5bab50936a63c29443d9176bfb684705 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Wed, 25 Feb 2026 16:23:46 +0100 Subject: [PATCH 1/9] 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 2/9] 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 3/9] 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 4/9] 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 5/9] 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 9798047c676ce0cadfd48774fc301c5fe141c34d Mon Sep 17 00:00:00 2001 From: "martin.kilbinger@cea.fr" Date: Fri, 26 Jun 2026 10:53:20 +0200 Subject: [PATCH 6/9] 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 7/9] 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 8/9] 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 9/9] 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