diff --git a/cosmo_val/cat_config.yaml b/cosmo_val/cat_config.yaml index cc1326df..cae6411b 100644 --- a/cosmo_val/cat_config.yaml +++ b/cosmo_val/cat_config.yaml @@ -980,7 +980,98 @@ SP_v1.4.11.3: e2_star_col: HSM_G2_STAR shear: R: 1.0 - path: /n17data/UNIONS/WL/v1.4.x/v1.4.11.3/unions_shapepipe_cut_struc_2024_v1.4.11.3.fits + 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 @@ -1038,11 +1129,12 @@ SP_v1.4.11.3_ecut07: 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_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 e1_col: e1_uncal @@ -1053,11 +1145,12 @@ SP_v1.4.6_uncal: path: unions_shapepipe_psf_2024_v1.4.a.fits 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 e1_col: e1_uncal @@ -1068,11 +1161,12 @@ SP_v1.4.6_uncal_w_iv: path: unions_shapepipe_psf_2024_v1.4.a.fits 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 e1_col: e1_uncal diff --git a/scripts/check_filled_fields.py b/scripts/check_filled_fields.py new file mode 100644 index 00000000..675ad116 --- /dev/null +++ b/scripts/check_filled_fields.py @@ -0,0 +1,74 @@ +#!/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 new file mode 100755 index 00000000..36a0a1b2 --- /dev/null +++ b/scripts/fill_photoz_bands.py @@ -0,0 +1,874 @@ +#!/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. + +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. + +: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 +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", + "Z_B_MAX", + "T_B", + "Z_ML", + "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", + "M_0", + "BPZ_FILT", + "BPZ_NONDETFILT", + "BPZ_FLAGFILT", +] + + +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": "UNIONS_DR6", + "checkpoint": "fill_photoz_bands_checkpoint.json", + "n_check_rows": 10, + "check_tol_arcsec": 0.5, + "verbose": False, + } + + short_options = { + "input": "-i", + "output": "-o", + "fits_dir": "-d", + "checkpoint": "-c", + } + + 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 + + +# --------------------------------------------------------------------------- +# 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. + ``('= 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). ``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 False, 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]] + + 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( + ra_hdf5[valid], dec_hdf5[valid], ra_fits[valid], dec_fits[valid] + ) + max_sep = float(np.max(sep)) + + return bool(max_sep <= tol_arcsec), n_checked, 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 + Compound HDF5 dataset opened in r+ mode. + hdf5_indices : numpy.ndarray + Row indices in dset corresponding to this tile (will be sorted). + fits_data : numpy.recarray + Data from the PhotoPipe FITS HDU. + valid_keys : list of str + Field names to copy from fits_data into dset. + + """ + sorted_idx = np.sort(hdf5_indices) + idx_min = int(sorted_idx[0]) + idx_max = int(sorted_idx[-1]) + n_range = idx_max - idx_min + 1 + + if n_range == len(sorted_idx): + # Contiguous block: single read-modify-write + chunk = dset[idx_min : idx_max + 1] + for key in valid_keys: + chunk[key] = fits_data[key] + dset[idx_min : idx_max + 1] = chunk + else: + # Non-contiguous: split into contiguous sub-blocks + gaps = np.where(np.diff(sorted_idx) > 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) + + # 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") + + +# --------------------------------------------------------------------------- +# 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 + 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 (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) + + 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 + + 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: + 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") + + 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) + + 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 + n_consec_fails += 1 + pbar.set_postfix( + {"done": n_processed, "size_err": n_skipped_size} + ) + continue + + # 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: + 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( + { + "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.") + 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: + 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}") + 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 or n_skipped_order > 0): + print( + "WARNING: no tiles were filled; all available tiles failed.", + file=sys.stderr, + ) + + return 0 + + +if __name__ == "__main__": + sys.exit(main(sys.argv)) 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..ae45a637 --- /dev/null +++ b/src/sp_validation/tests/test_fill_photoz_bands.py @@ -0,0 +1,304 @@ +"""``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 +import warnings +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) + + +@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