From 4fd191f00a1afd0aa5d8f884259b74c14f106134 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 20:13:03 -0400 Subject: [PATCH 01/17] catalog: layout-agnostic campaign and star catalogue readers Replace read_hdf5_file's hardcoded patches// lookup with find_dataset_group(), which descends from the file root through single container groups until it reaches the per-unit datasets. This reads the legacy patches// layout that ShapePipe still writes as a compatibility shim, a future flat tiles/ layout, and the exposures/ layout of full_starcat_.hdf5 with the same code. read_star_catalogue() keeps the FITS path for files ending in .fits. Requested columns missing from the data now raise a clear KeyError. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- src/sp_validation/catalog.py | 213 +++++++++++++++++++++++++++-------- 1 file changed, 168 insertions(+), 45 deletions(-) diff --git a/src/sp_validation/catalog.py b/src/sp_validation/catalog.py index e39d2076..2c4ebadf 100644 --- a/src/sp_validation/catalog.py +++ b/src/sp_validation/catalog.py @@ -758,71 +758,194 @@ def read_param_file(path, verbose=False): return param_list_unique -def read_hdf5_file(file_path, name, stats_file, check_only=False, param_path=None): - """Read HDF5 File. +#: Columns written by ShapePipe v2's ``MergeStarCatPSFEX`` into +#: ``full_starcat_.hdf5`` (one dataset per exposure). +STAR_CAT_COLUMNS = ( + "X", + "Y", + "RA", + "DEC", + "HSM_G1_PSF", + "HSM_G2_PSF", + "HSM_T_PSF", + "HSM_G1_STAR", + "HSM_G2_STAR", + "HSM_T_STAR", + "HSM_FLAG_PSF", + "HSM_FLAG_STAR", + "MAG", + "SNR", + "ACCEPTED", + "CCD_NB", +) + + +def find_dataset_group(hdf5_file): + """Find Dataset Group. + + Descend from the root of an open HDF5 file to the single group whose + members are the per-unit datasets (one per tile, or one per exposure). + + This makes the reader independent of how deeply the products nest that + group: it walks down as long as the current node holds exactly one + sub-group, and stops as soon as the members are datasets. It therefore + reads both the legacy ``patches//`` layout (the + "patches" key is a ShapePipe-side compatibility shim, not a concept) and + a flat ``tiles/`` or ``exposures/`` layout. - Read hdf5 file and return contained data. + Parameters + ---------- + hdf5_file : h5py.File or h5py.Group + open input file + + Returns + ------- + h5py.Group + group whose members are the per-unit datasets + + Raises + ------ + ValueError + if the file is empty, or a level holds more than one sub-group + + """ + node = hdf5_file + while True: + keys = list(node) + if not keys: + raise ValueError(f"No data found under {node.name!r} in {hdf5_file.file.filename}") + if all(isinstance(node[key], h5py.Dataset) for key in keys): + return node + if len(keys) != 1: + raise ValueError( + f"Expected a single container group under {node.name!r} in" + + f" {hdf5_file.file.filename}, found {len(keys)}: {keys[:5]}" + ) + node = node[keys[0]] + + +def _check_columns(dtype, param_list, file_path): + """Raise a clear error if requested columns are absent from the data.""" + missing = [col for col in param_list if col not in (dtype.names or ())] + if missing: + raise KeyError( + f"Column(s) {missing} not found in catalogue {file_path}." + + f" Available columns: {sorted(dtype.names or ())}" + ) + + +def concatenate_datasets(group, param_list=None, file_path="", verbose=True): + """Concatenate Datasets. + + Concatenate every dataset of an HDF5 group into one structured array, + optionally restricted to a list of columns. + + Parameters + ---------- + group : h5py.Group + group whose members are structured datasets + param_list : list of str, optional + columns to keep; default is ``None`` (keep all) + file_path : str, optional + input file path, for error messages + verbose : bool, optional + verbose output if ``True`` + + Returns + ------- + numpy.ndarray + concatenated structured array + + """ + keys = list(group) + if param_list is not None: + _check_columns(group[keys[0]].dtype, param_list, file_path) + + n_rows = sum(group[key].shape[0] for key in keys) + if verbose: + n_cols = len(param_list) if param_list is not None else len(group[keys[0]].dtype) + print( + f"Reading {len(keys)} datasets," + + f" estimating {n_cols * n_rows * 8 / 1024**3:.1f}" + + f" Gb memory for the ({n_cols} x {n_rows}) data array ..." + ) + + data_list = [] + for key in tqdm.tqdm(keys, disable=not verbose): + data = group[key][()] + if param_list is not None: + data = data[param_list] + data_list.append(data) + + return np.concatenate(data_list, axis=0) + + +def read_campaign_catalogue( + file_path, + param_path=None, + param_list=None, + verbose=True, +): + """Read Campaign Catalogue. + + Read a campaign galaxy catalogue (``final_cat_.hdf5``) and + return its per-tile datasets concatenated into one structured array. Parameters ---------- file_path : str input file path - name : str - patch name - stats_file : file handler - summary statistics output file handler - check_only : bool, optional - If True only check, not return data + param_path : str, optional + path to a parameter file listing the columns to keep + param_list : list of str, optional + columns to keep; takes precedence over ``param_path`` + verbose : bool, optional + verbose output if ``True`` Returns ------- - dict - data + numpy.ndarray + catalogue data """ - param_list = read_param_file(param_path, verbose=True) if param_path else None + if param_list is None and param_path: + param_list = read_param_file(param_path, verbose=verbose) with h5py.File(file_path, "r") as hdf5_file: - # Find patch group in hierarchical structure - if f"patches/{name}" not in hdf5_file: - raise KeyError(f"Entry patches/{name} not found in file {file_path}") - patch_group = hdf5_file[f"patches/{name}"] - - # Get size of data array - num_rows = sum(patch_group[ID].shape[0] for ID in patch_group) - # num_cols = patch_group[next(iter(patch_group))].shape[1] - num_cols = len(param_list) - - print( - f"Estimating {num_cols * num_rows * 8 / 1024**3:.1f}" - + f" Gb memory for the ({num_cols} x {num_rows}) data array ..." + group = find_dataset_group(hdf5_file) + return concatenate_datasets( + group, param_list=param_list, file_path=file_path, verbose=verbose ) - # data_comb = np.memmap(output_file, dtype=patch_group[next(iter(patch_group))].dtype, - # mode="w+", shape=(num_rows, num_cols)) - data_list = [] - ID_pbl = set() - for ID in tqdm.tqdm(patch_group): - # Get data for this ID from file - data = patch_group[ID][()] - # Restrict to parameter list if given - data = data[param_list] if param_list is not None else data +def read_star_catalogue(file_path, hdu=1, verbose=True): + """Read Star Catalogue. - if not check_only: - # Add new to existing data - data_list.append(data) + Read a campaign star/PSF catalogue. Reads the ShapePipe v2 + ``full_starcat_.hdf5`` (one dataset per exposure), or a legacy + FITS star catalogue when ``file_path`` ends in ``.fits``. + + Parameters + ---------- + file_path : str + input file path + hdu : int, optional + HDU number for the FITS path; default is 1 + verbose : bool, optional + verbose output if ``True`` - print("Combine tile catalogues") - data_comb = np.concatenate(data_list, axis=0) - print("Done") + Returns + ------- + numpy.ndarray + star catalogue data - # Print problematic tile IDs - for ID in ID_pbl: - print("Tile IDs with missing keys:", file=stats_file) - print(ID, file=stats_file) + """ + if str(file_path).endswith(".fits"): + return fits.getdata(file_path, hdu) - return data_comb + with h5py.File(file_path, "r") as hdf5_file: + group = find_dataset_group(hdf5_file) + return concatenate_datasets(group, file_path=file_path, verbose=verbose) def get_maked_col(dat, col, mask): From f5c3b2e079feac9d43bd3ebd8e861f73e45dff2a Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 20:14:27 -0400 Subject: [PATCH 02/17] galaxy: replace IMAFLAGS_ISO cut with config-driven mask columns MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ShapePipe v2 drops IMAFLAGS_ISO for eleven boolean MASK_n* columns. galaxy.mask_cut() ORs a configurable list of them (default MASK_n4, MASK_n1, MASK_n2, MASK_n8, MASK_n1024 — stars, star halos, manual galaxy mask, MaxiMask) and returns the keep mask; a catalogue missing any of the requested columns raises a KeyError naming them. classification_galaxy_base takes mask_columns; extract_info.py passes the params.py mask_columns list and uses it for the star-sample cut too, and now reads both catalogues through the new readers. Column lists in params.py and masking.py's SPATIAL_CUTS updated to the MASK_n* names. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- scripts/calibration/extract_info.py | 10 ++- scripts/calibration/params.py | 25 +++++- scripts/glass_mock/compute_leakage_harmony.py | 4 +- scripts/masking.py | 12 ++- src/sp_validation/catalog_builders.py | 1 - src/sp_validation/galaxy.py | 90 ++++++++++++++++++- 6 files changed, 132 insertions(+), 10 deletions(-) diff --git a/scripts/calibration/extract_info.py b/scripts/calibration/extract_info.py index 123c7437..4d449416 100644 --- a/scripts/calibration/extract_info.py +++ b/scripts/calibration/extract_info.py @@ -38,6 +38,7 @@ # from sp_validation.catalog import * from sp_validation import catalog as spv_cat +from sp_validation import galaxy from sp_validation.calibration import * from sp_validation.calibration import metacal from sp_validation.galaxy import * @@ -69,8 +70,8 @@ dd = np.load(galaxy_cat_path, mmap_mode=mmap_mode) else: print("Loading galaxy .hdf5 file...") - dd = spv_cat.read_hdf5_file( - galaxy_cat_path, name, stats_file, param_path=param_list_path + dd = spv_cat.read_campaign_catalogue( + galaxy_cat_path, param_path=param_list_path ) n_obj = len(dd) @@ -116,7 +117,7 @@ # ### Load star catalogue if star_cat_path: - d_star = fits.getdata(star_cat_path, hdu_star_cat) + d_star = spv_cat.read_star_catalogue(star_cat_path, hdu=hdu_star_cat) if star_cat_path: print_stats("Stars:", stats_file, verbose=verbose) @@ -160,7 +161,7 @@ m_star = ( (dd["FLAGS"][ind_star] == 0) - & (dd["IMAFLAGS_ISO"][ind_star] == 0) + & galaxy.mask_cut(dd, mask_columns)[ind_star] & (dd["NGMIX_MCAL_FLAGS"][ind_star] == 0) & (dd["NGMIX_G1_PSF_ORIG_NOSHEAR"][ind_star] != -10) ) @@ -311,6 +312,7 @@ gal_mag_faint=gal_mag_faint, flags_keep=flags_keep, n_epoch_min=n_epoch_min, + mask_columns=mask_columns, ) if shape == "ngmix": m_gal = classification_galaxy_ngmix( diff --git a/scripts/calibration/params.py b/scripts/calibration/params.py index ea871119..db45e55b 100644 --- a/scripts/calibration/params.py +++ b/scripts/calibration/params.py @@ -121,11 +121,30 @@ "NGMIX_T_PSF_RECONV_NOSHEAR", ] +## ShapePipe v2 mask columns OR'd together for the galaxy selection cut +mask_columns = [ + "MASK_n4", + "MASK_n1", + "MASK_n2", + "MASK_n8", + "MASK_n1024", +] + ## Pre-calibration catalogue, including masked objects and mask flags add_cols_pre_cal = [ "TILE_ID", "NUMBER", - "IMAFLAGS_ISO", + "MASK_n1", + "MASK_n2", + "MASK_n4", + "MASK_n8", + "MASK_n16", + "MASK_n32", + "MASK_n64", + "MASK_n128", + "MASK_n256", + "MASK_n1024", + "MASK_n2048", "FLAGS", "NGMIX_MCAL_FLAGS", "NGMIX_MCAL_TYPES_FAIL", @@ -141,7 +160,6 @@ add_cols_pre_cal_format = {} for key in ( "NUMBER", - "IMAFLAGS_ISO", "FLAGS", "NGMIX_MCAL_FLAGS", "NGMIX_MCAL_TYPES_FAIL", @@ -150,6 +168,9 @@ ): add_cols_pre_cal_format[key] = "I" +for key in mask_columns: + add_cols_pre_cal_format[key] = "L" + add_cols_pre_cal_format["TILE_ID"] = "A7" add_cols_pre_cal_format["NUMBER"] = "J" diff --git a/scripts/glass_mock/compute_leakage_harmony.py b/scripts/glass_mock/compute_leakage_harmony.py index a63c6bd4..22c2439d 100644 --- a/scripts/glass_mock/compute_leakage_harmony.py +++ b/scripts/glass_mock/compute_leakage_harmony.py @@ -8,6 +8,8 @@ import numpy as np from astropy.io import fits +from sp_validation import catalog as spv_cat + from sp_validation.glass_mock import compute_leakage_harmony @@ -53,7 +55,7 @@ def get_parser(): print("Catalog data loaded successfully.") print("Loading the star catalog data...") - cat_star = fits.getdata(f"{args.star_cat_path}") + cat_star = spv_cat.read_star_catalogue(args.star_cat_path) print("Star catalog data loaded successfully.") print("Computing leakage in harmonic space...") diff --git a/scripts/masking.py b/scripts/masking.py index 32a970ed..0ef0491f 100644 --- a/scripts/masking.py +++ b/scripts/masking.py @@ -16,7 +16,17 @@ # the footprint definition. SPATIAL_CUTS = { "overlap", - "IMAFLAGS_ISO", + "MASK_n1", + "MASK_n2", + "MASK_n4", + "MASK_n8", + "MASK_n16", + "MASK_n32", + "MASK_n64", + "MASK_n128", + "MASK_n256", + "MASK_n1024", + "MASK_n2048", "N_EPOCH", "4_Stars", "8_Manual", diff --git a/src/sp_validation/catalog_builders.py b/src/sp_validation/catalog_builders.py index 78dc7838..eaea4fb5 100644 --- a/src/sp_validation/catalog_builders.py +++ b/src/sp_validation/catalog_builders.py @@ -443,7 +443,6 @@ def dtype_out(self, name, dtype_in): "RA", "Dec", "FLAGS", - "IMAFLAGS_ISO", "NUMBER", ] if dtype_in.kind == "U": diff --git a/src/sp_validation/galaxy.py b/src/sp_validation/galaxy.py index 3c49a96b..3eda60af 100644 --- a/src/sp_validation/galaxy.py +++ b/src/sp_validation/galaxy.py @@ -28,6 +28,87 @@ # required square root: FWHM = 2.35482 sqrt(T / 2) from sp_validation import io +#: All mask columns written by ShapePipe v2 (bool, ``True`` = masked). +#: n4 stars; n1/n2 faint/bright star halos; n8 manual galaxy mask; +#: n1024 MaxiMask; n16..n256 per-band coverage; n2048 no PS-z2 coverage. +#: These replace the single IMAFLAGS_ISO bitmask of ShapePipe v1. +MASK_COLUMNS = ( + "MASK_n1", + "MASK_n2", + "MASK_n4", + "MASK_n8", + "MASK_n16", + "MASK_n32", + "MASK_n64", + "MASK_n128", + "MASK_n256", + "MASK_n1024", + "MASK_n2048", +) + +#: Mask columns OR'd together for the default galaxy selection. Deliberately +#: not a blanket OR over MASK_COLUMNS: the per-band coverage columns +#: (n16..n256) and n2048 would mask essentially the whole catalogue. +DEFAULT_MASK_COLUMNS = ( + "MASK_n4", + "MASK_n1", + "MASK_n2", + "MASK_n8", + "MASK_n1024", +) + + +def _column_names(dd): + """Return the column names of a structured array or mapping.""" + dtype = getattr(dd, "dtype", None) + if dtype is not None and dtype.names is not None: + return tuple(dtype.names) + return tuple(dd.keys()) + + +def mask_cut(dd, mask_columns=None): + """Mask Cut. + + Return a boolean mask that is ``True`` for objects *not* flagged by any + of the requested ShapePipe mask columns. + + Parameters + ---------- + dd : numpy.ndarray or dict + input catalogue + mask_columns : list of str, optional + mask columns to OR together; default is ``DEFAULT_MASK_COLUMNS`` + + Returns + ------- + numpy.ndarray + boolean mask, ``True`` = keep + + Raises + ------ + KeyError + if any requested mask column is absent from the catalogue + + """ + columns = list(DEFAULT_MASK_COLUMNS if mask_columns is None else mask_columns) + + available = _column_names(dd) + missing = [col for col in columns if col not in available] + if missing: + raise KeyError( + f"Mask column(s) {missing} not found in catalogue." + + " ShapePipe v2 catalogues carry the boolean columns" + + f" {list(MASK_COLUMNS)}; ShapePipe v1 catalogues carry" + + " IMAFLAGS_ISO instead and are not supported." + + f" Available columns: {sorted(available)}" + ) + + masked = np.zeros(len(dd[columns[0]]), dtype=bool) + for col in columns: + masked |= np.asarray(dd[col], dtype=bool) + + return ~masked + def classification_galaxy_overlap_ra_dec(dd, ra_key="XWIN_WORLD", dec_key="YWIN_WORLD"): """Classification Galaxy Overlap Ra Dec. @@ -141,11 +222,18 @@ def classification_galaxy_base( gal_mag_faint=26, flags_keep=None, n_epoch_min=1, + mask_columns=None, ): """Classification Galaxy Base. Return mask corresponding to basic classification for galaxies. + Parameters + ---------- + mask_columns : list of str, optional + ShapePipe mask columns OR'd together to reject masked objects; + default is ``DEFAULT_MASK_COLUMNS`` + """ # SExtractor flags # Keep some flags if specified @@ -172,7 +260,7 @@ def classification_galaxy_base( & cut_flags & (dd["MAG_AUTO"] <= gal_mag_faint) & (dd["MAG_AUTO"] >= gal_mag_bright) - & (dd["IMAFLAGS_ISO"] == 0) + & mask_cut(dd, mask_columns) & (dd["N_EPOCH"] >= n_epoch_min) ) From af9949127b7091a27be77e8ca0b84c29f9d354d7 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 20:22:06 -0400 Subject: [PATCH 03/17] retire patch logic (#340) ShapePipe v2 processes a campaign (a tile list); P1-P7 no longer exist. - catalog_builders.JointCat: get_patches()/get_n_obj() and the per-patch FITS merge are replaced by a merge over a list of campaign hdf5 files (-i final_cat_A.hdf5+final_cat_B.hdf5), read through read_campaign_catalogue. The 'patch' int8 column becomes a 'campaign' string column; the hdf5 root attr becomes 'campaigns'. - survey.get_footprint() deleted: it was a lookup table of P1-P7 (plus W3) RA/Dec boundaries, meaningless for a campaign. Its only caller, catalog.check_matching, used it as an optional pre-filter via a 'name' argument that every caller passed as None; the argument goes too, along with the test_survey test that exercised P5. - merge_psf_cat.py, combine_results.py, stats_tile_id_gal_counts.py, compute_area.py: patch vocabulary and v1/v1.5/v1.6 P-name shortcuts generalised to an explicit list of campaigns. - params.py: 'name = "P7"' becomes 'campaign = None'. Deleted (only ever meaningful for the P1-P7 era): - scripts/prepare_patch_for_spval.sh: symlinks a v1 per-patch run tree (~/psfex/${patch}/output/run_sp_Ms/.../full_starcat-0000000.fits, tiles_${patch}.txt) into a working dir; neither the layout nor the file names exist in v2. - scripts/plot_rho_stats_patches.py: globs P* directories and reads P*/output/run_sp_Pl/mccd_plots_runner/output/rho_stats_id.fits, one curve per patch. No campaign analogue. - scripts/survey_stats_all.sh: hardcoded 'for patch in P1 ... P7' over v1 run-directory bookkeeping. - scripts/star_match_stats.py: sums stats_file.txt over the seven patches. - scripts/check_tile_IDs_SP_LF.py: compares per-patch ShapePipe tile IDs against CFIS3500_THELI_P.list Lensfit files; both sides P-named. The word 'patch' survives only in catalog.py's comment naming the legacy hdf5 group, and in the treecorr jackknife sense (npatch/patch_number). Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- docs/source/post_processing.md | 18 +- docs/source/run_validation.md | 7 +- docs/source/using_the_catalogues.md | 4 +- scripts/calibration/extract_info.py | 3 +- scripts/calibration/params.py | 5 +- scripts/check_tile_IDs_SP_LF.py | 83 ------- scripts/combine_results.py | 272 ++++++++++---------- scripts/compute_area.py | 10 +- scripts/merge_psf_cat.py | 71 +++--- scripts/plot_rho_stats_patches.py | 98 -------- scripts/prepare_patch_for_spval.sh | 31 --- scripts/star_match_stats.py | 44 ---- scripts/stats_tile_id_gal_counts.py | 43 ++-- scripts/survey_stats_all.sh | 94 ------- src/sp_validation/catalog.py | 10 +- src/sp_validation/catalog_builders.py | 327 +++++++++---------------- src/sp_validation/galaxy.py | 4 + src/sp_validation/rho_tau.py | 2 +- src/sp_validation/survey.py | 71 ------ src/sp_validation/tests/test_survey.py | 8 - workflow/image_sims/params_im_sim.py | 10 +- 21 files changed, 338 insertions(+), 877 deletions(-) delete mode 100644 scripts/check_tile_IDs_SP_LF.py delete mode 100755 scripts/plot_rho_stats_patches.py delete mode 100755 scripts/prepare_patch_for_spval.sh delete mode 100644 scripts/star_match_stats.py delete mode 100644 scripts/survey_stats_all.sh diff --git a/docs/source/post_processing.md b/docs/source/post_processing.md index 6d2f2423..d16e2e33 100644 --- a/docs/source/post_processing.md +++ b/docs/source/post_processing.md @@ -3,8 +3,8 @@ ## Science-ready catalogue production Processing steps of `ShapePipe` output catalogues carried out by the `sp_validation` package to produce science-ready catalogues are: -1. Extract relevant information from a final `ShapePipe` output catalogue per patch; run basic diagnostic tests, create pre-calibration shear catalogues. -2. Merge pre-calibration catalogues created in the previous step, e.g. processed by individual patches, into one or more joint catalogues; +1. Extract relevant information from a final `ShapePipe` output catalogue per campaign; run basic diagnostic tests, create pre-calibration shear catalogues. +2. Merge pre-calibration catalogues created in the previous step, e.g. processed as individual campaigns, into one or more joint catalogues; 3. Apply external area and footprint masks. These are the "structural" and the coverage masks. 4. Create calibrated galaxy shear catalogue. This step includes the tasks: a. Mask objects using flags and criteria in `ShapePipe` output catalogues and external (e.g. mask) files; @@ -19,7 +19,7 @@ This is performed (version > v1.4.1, < v2.0) with the python script `scripts/cal This script creates three shear catalogues in FITS format: - _Basic_ catalogue containing - positions, shapes (calibrated + PSF-leakage corrected), weights (DES), magnitude, patch ID. Masking and galaxy selection are applied. + positions, shapes (calibrated + PSF-leakage corrected), weights (DES), magnitude, campaign ID. Masking and galaxy selection are applied. - _Extended_ catalogue containing **in addition** uncalibrated shapes inverse-variance weights, shear response matrices, SNR, flux, size, PSF quantities. Masking and galaxy selection are applied. - _Comprehensive_ catalogue containing **in addition** @@ -27,11 +27,11 @@ This script creates three shear catalogues in FITS format: This catalogue does not contain calibrated shear estimates, since the calibration is carried out after applying masking and selection. This is the main output catalogue that will be processed further. -This step is carried out per patch. Parameters have to be set via the python configuration file `params.py` (template at `scripts/calibration/params.py`). +This step is carried out per campaign. Parameters have to be set via the python configuration file `params.py` (template at `scripts/calibration/params.py`). ### 2. Merge catalogues -The patch-wise comprehensive catalogues extracted in the previous step are merged using the script `scripts/calibration/create_joint_comprehensive_cat.py`, which is a front-end +The per-campaign comprehensive catalogues extracted in the previous step are merged using the script `scripts/calibration/create_joint_comprehensive_cat.py`, which is a front-end of the `sp_validation` library class `catalog_builders:JointCat`. ### 3. Apply external masks @@ -57,8 +57,8 @@ The following describes the pre-v1.4.2 method to create a joint, calibrated shea Summary statistics created by shear validation runs of sub-areas of a survey can be combined to create joint summary statistics. This is useful in cases where the galaxy catalogue of an entire survey is too large to process, and -needs to be broken down in smaller patches. This step provides global summary -statistics from those patches. +needs to be broken down into smaller campaigns. This step provides global +summary statistics from those campaigns. Depending on the type of summary, their combination can be the sum (e.g. for number of objects), average, weighted average (e.g. for the additive bias), the @@ -66,7 +66,7 @@ weighted average of the square (e.g. the ellipticity dispersion), the weighted variance (to combine variance estimates), or the weighted variance of the mean (to combine mean variance estimates). -In a directory containing the subpatches as subdirectories, and within each +In a directory containing the campaigns as subdirectories, and within each their own output directory (`sp_output`by default in `params.py`) with results of the validation runs, type ```bash @@ -89,7 +89,7 @@ calibration outputs can be used to create a combined, globally calibrated shear catalogue. The calibration is obtained from the files `R.txt` and `c.txt` created above. -In the same directory containing the subpatches as above, type +In the same directory containing the campaigns as above, type ```bash create_joint_shape_cat.py ``` diff --git a/docs/source/run_validation.md b/docs/source/run_validation.md index 7ee06767..8d3fded1 100644 --- a/docs/source/run_validation.md +++ b/docs/source/run_validation.md @@ -12,7 +12,7 @@ including the sheared values for metacalibration. All inputs and settings are contained in the python configuration script `scripts/calibration/params.py`, that needs to be edited accordingly. The main parameters are: -- `name`: field or patch name, can be any string. E.g. `P3` for patch 3. +- `campaign`: campaign name (the ShapePipe tile list), can be any string. - `data_dir`: input directory for data. Set to `.` for validation run in current directory. - `galaxy_cat_path`: path to galaxy catalogue, format `.fits`. or `.hdf5`. @@ -27,8 +27,9 @@ Optional parameters are: - `mask_external_path`: path to external mask file, format `.reg`. Set to `None` if not required. -See the script `prepare_patch_for_spval.sh` for an example of copying -the required input files to where the validation is to be run. +Link or copy the campaign's merged products -- `final_cat_.hdf5` +and `full_starcat_.hdf5` -- into the directory where the validation +is to be run. ### Run diff --git a/docs/source/using_the_catalogues.md b/docs/source/using_the_catalogues.md index 72781c12..db408c65 100644 --- a/docs/source/using_the_catalogues.md +++ b/docs/source/using_the_catalogues.md @@ -16,7 +16,7 @@ how to apply the metacalibration corrections yourself. ```{note} The examples below target catalogue **v1.0** (April 2022), which is distributed as FITS. From ShapePipe catalogue **v1.4.1** onward the merged catalogues ship -as HDF5 instead; open those with {func}`sp_validation.io.read_hdf5_file` (or +as HDF5 instead; open those with {func}`sp_validation.catalog.read_campaign_catalogue` (or `h5py` / `astropy`) in place of `astropy.io.fits` below — the column names and the calibration recipe are unchanged. ``` @@ -167,7 +167,7 @@ mask = np.full(len(data_ext), True) # Other examples: # mask = data_ext['mask_extern'] == 0 # LensFit-unmasked regions -# mask = data_ext['patch'] == 3 # patch P3 +# mask = data_ext['campaign'] == b'W3' # objects from campaign W3 # mask = data_ext['mag'] < 23.5 # r-band magnitude cut n_kept, n_all = np.count_nonzero(mask), len(data_ext) diff --git a/scripts/calibration/extract_info.py b/scripts/calibration/extract_info.py index 4d449416..b65c1975 100644 --- a/scripts/calibration/extract_info.py +++ b/scripts/calibration/extract_info.py @@ -149,7 +149,6 @@ [col_name_ra, col_name_dec], thresh, stats_file, - name=None, verbose=verbose, ) @@ -492,7 +491,7 @@ write_tile_id_gal_counts(detection_IDs, galaxy_IDs, shape_IDs, fname) # + -# Add all weights (for combining weighted averages of subpatches) +# Add all weights (for combining weighted averages of sub-samples) w_tot = np.sum(w) diff --git a/scripts/calibration/params.py b/scripts/calibration/params.py index db45e55b..e4866942 100644 --- a/scripts/calibration/params.py +++ b/scripts/calibration/params.py @@ -25,9 +25,8 @@ # Survey parameters -## Field or patch name. Put None if n/a -name = "P7" -print("Field name = {}".format(name)) +## Campaign name (the tile list processed by ShapePipe). Put None if n/a +campaign = None ## Area of a tile in deg^2 area_tile = 0.25 diff --git a/scripts/check_tile_IDs_SP_LF.py b/scripts/check_tile_IDs_SP_LF.py deleted file mode 100644 index 54f23fda..00000000 --- a/scripts/check_tile_IDs_SP_LF.py +++ /dev/null @@ -1,83 +0,0 @@ -import re -import sys - -import numpy as np - - -def main(argv=None): - - survey = "v1" - - if survey == "v1": - n_patch = 7 - patches = [f"P{x}" for x in np.arange(n_patch) + 1] - - IDs_sp_base = "found_ID_wshapes.txt" - tile_ID_gal_counts_sp_base = "tile_id_gal_counts_ngmix.txt" - - n_SP_not_in_LF_all = 0 - n_LF_not_in_SP_all = 0 - - for patch in patches: - print(patch) - - # ShapePipe - path = f"{patch}/sp_output/{IDs_sp_base}" - - with open(path) as f: - dat = f.readlines() - ID_SP = [] - for line in dat: - ID_SP.append(line.rstrip()) - print(f" #SP = {len(ID_SP)}") - - # LensFit - path = f"CFIS3500_THELI_{patch}.list" - with open(path) as f: - dat = f.readlines() - ID_LF = [] - for line in dat: - m = re.match(r".*CFIS\.(\d{3}\.\d{3})\.r", line) - if m: - ID_LF.append(m[1].rstrip()) - print(f" #LF = {len(ID_LF)}") - - # tile stats for SP - dat = np.loadtxt(f"{patch}/sp_output/{tile_ID_gal_counts_sp_base}") - tile_ID = [] - for my_ID in dat[:, 0]: - tile_ID.append(f"{my_ID:07.3f}") - - # ShapePipe tiles not contained in LensFit - n_SP_not_in_LF = 0 - for ID in ID_SP: - if ID not in ID_LF: - # Print number of galaxies on those tiles not contained in LF - n_SP_not_in_LF += 1 - if n_SP_not_in_LF == -1: - print(f" SP {ID} not in LF") - n_SP_not_in_LF_all += n_SP_not_in_LF - - # LensFit tiles not contained in ShapePipe - n_LF_not_in_SP = 0 - for ID in ID_LF: - if ID not in ID_SP: - n_LF_not_in_SP += 1 - if n_LF_not_in_SP == -1: - print(f" LF {ID} not in SP") - print(f" [{ID}]") - n_LF_not_in_SP_all += n_LF_not_in_SP - - print(f" # SP not in LF = {n_SP_not_in_LF}") - print(f" # LF not in SP = {n_LF_not_in_SP}") - - print() - print("All patches") - print(f" # SP not in LF = {n_SP_not_in_LF_all}") - print(f" # LF not in SP = {n_LF_not_in_SP_all}") - - return 0 - - -if __name__ == "__main__": - sys.exit(main(sys.argv)) diff --git a/scripts/combine_results.py b/scripts/combine_results.py index 636d5c07..50cc91c5 100755 --- a/scripts/combine_results.py +++ b/scripts/combine_results.py @@ -10,11 +10,11 @@ from astropy.io import ascii -def get_match(stats_files, patch, pattern, previous=None, n_previous=[1], typ=str): +def get_match(stats_files, campaign, pattern, previous=None, n_previous=[1], typ=str): prev_ok = False - for idx, line in enumerate(stats_files[patch]): + for idx, line in enumerate(stats_files[campaign]): m = re.search(pattern, line) if m: if (previous and prev_ok) or not previous: @@ -29,19 +29,19 @@ def get_match(stats_files, patch, pattern, previous=None, n_previous=[1], typ=st for prev, n_prev in zip(previous, n_previous): # Line index n_previous earlier, +1 since next line will be read in # next loop; look for pattern in previous line - m_prev = re.search(prev, stats_files[patch][idx - n_prev + 1]) + m_prev = re.search(prev, stats_files[campaign][idx - n_prev + 1]) if m_prev: prev_ok = True raise ValueError( - f"No match of '{pattern}' in patch {patch} (prev='{previous}'), n_prev={n_previous}" + f"No match of '{pattern}' in campaign {campaign} (prev='{previous}'), n_prev={n_previous}" ) -def read_stats_files(patches, path, verbose=False): +def read_stats_files(campaigns, path, verbose=False): stats_files = {} - for p in patches: + for p in campaigns: fname = f"{p}/{path}" if os.path.exists(fname): if verbose: @@ -115,7 +115,7 @@ def combine(results): # Weight values w = np.array(list(results["value"][key_w].values())) - # Patch mean values + # Campaign mean values m = np.array(list(results["value"][key_m].values())) # Overall mean @@ -142,7 +142,7 @@ def combine(results): # Weight values w = np.array(list(results["value"][key_w].values())) - # Patch mean values + # Campaign mean values m = np.array(list(results["value"][key_m].values())) # Overall mean @@ -177,18 +177,18 @@ def print_all( # Header if header: - print("# patch", " " * 3, end=" ", file=fout) + print("# campaign", " " * 3, end=" ", file=fout) for key in keys: print(f"{key:>11s}", end=" ", file=fout) print(file=fout) - # Loop over patches - for patch in stats_files: - print(f"{patch:11s}", end=" ", file=fout) + # Loop over campaigns + for campaign in stats_files: + print(f"{campaign:11s}", end=" ", file=fout) # Write value for each key for key in keys: - val = results["value"][key][patch] + val = results["value"][key][campaign] if key == "N_gal": print(f"{val:>11.0f}", end=" ", file=fout) elif key in ("w_tot", "n_gal_am2"): @@ -222,7 +222,7 @@ def get_area(fname): with open(fname) as f: lines = f.readlines() for line in lines: - m = re.search("nmasked patch area without overlap = (.*) deg", line) + m = re.search("nmasked campaign area without overlap = (.*) deg", line) if m: return float(m[1]) @@ -234,30 +234,30 @@ def get_area(fname): def get_values(results, stats_files, shape, use_keys, area_deg2=-1): """Get Values - Get values from stats files for all patches. + Get values from stats files for all campaigns. Parameters ---------- results : dict results dictionary stats_files : dict of array of str - stats files content for each patch + stats files content for each campaign shape : str shape measurement method use_keys : dict keys to include area_deg2 : float, optional area in square degree, optional is -1 (to be retrieved - for each patch) + for each campaign) """ # Number of galaxies key = "N_gal" if use_keys[key]: init_key(results, key, "sum") - for patch in stats_files: - results["value"][key][patch] = get_match( + for campaign in stats_files: + results["value"][key][campaign] = get_match( stats_files, - patch, + campaign, r"Number of galaxies after metacal = (\d+)/", previous=[f"^{shape}$"], typ=int, @@ -266,15 +266,15 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): area_deg2_tot = 0 key_der = "n_gal_am2" init_key(results, key_der, "w_avg", extra="N_gal") - for patch in stats_files: + for campaign in stats_files: if area_deg2 < 0: - area_deg2_patch = get_area(f"{patch}/area.txt") - area_deg2_tot += area_deg2_patch - print(f"area({patch}) = {area_deg2_patch} deg^2") + area_deg2_campaign = get_area(f"{campaign}/area.txt") + area_deg2_tot += area_deg2_campaign + print(f"area({campaign}) = {area_deg2_campaign} deg^2") else: - area_deg2_patch = area_deg2 - results["value"][key_der][patch] = results["value"]["N_gal"][patch] / ( - area_deg2_patch * 3600 + area_deg2_campaign = area_deg2 + results["value"][key_der][campaign] = results["value"]["N_gal"][campaign] / ( + area_deg2_campaign * 3600 ) if area_deg2 < 0: @@ -285,10 +285,10 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): key = "w_tot" if use_keys[key]: init_key(results, key, "sum") - for patch in stats_files: - results["value"][key][patch] = get_match( + for campaign in stats_files: + results["value"][key][campaign] = get_match( stats_files, - patch, + campaign, r"Sum of weights = (\S+)", typ=float, previous=[f"^{shape}$"], @@ -300,16 +300,16 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): for comp in (1, 2): key = f"{key_base}{comp}" init_key(results, key, "w_avg", extra="N_gal") - for patch in stats_files: + for campaign in stats_files: c = get_match( stats_files, - patch, + campaign, rf"{key_base}{comp} = (\S+)", previous=[f"^{shape}:$"], n_previous=[2 * comp - 1], typ=float, ) - results["value"][key][patch] = c + results["value"][key][campaign] = c # Additive bias (unweighted error) key_base = "dc_" @@ -317,16 +317,16 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): for comp in (1, 2): key = f"{key_base}{comp}" init_key(results, key, "var", extra=["N_gal", f"c_{comp}"]) - for patch in stats_files: + for campaign in stats_files: dc = get_match( stats_files, - patch, + campaign, rf"{key_base}{comp} = (\S+)", typ=float, previous=[f"^{shape}:$"], n_previous=[2 * comp + 7], ) - results["value"][key][patch] = dc + results["value"][key][campaign] = dc # Additive bias (unweighted error of mean) key_base = "dmc_" @@ -334,16 +334,16 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): for comp in (1, 2): key = f"{key_base}{comp}" init_key(results, key, "var_m", extra=["N_gal", f"c_{comp}"]) - for patch in stats_files: + for campaign in stats_files: dmc = get_match( stats_files, - patch, + campaign, rf"{key_base}{comp} = (\S+)", typ=float, previous=[f"^{shape}:$"], n_previous=[2 * comp + 7], ) - results["value"][key][patch] = dmc + results["value"][key][campaign] = dmc # Additive bias (weighted mean) key_base = "cw_" @@ -351,16 +351,16 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): for comp in (1, 2): key = f"{key_base}{comp}" init_key(results, key, "w_avg", extra="w_tot") - for patch in stats_files: + for campaign in stats_files: c = get_match( stats_files, - patch, + campaign, rf"{key_base}{comp} = (\S+)", previous=[f"^{shape}:$"], n_previous=[comp * 2], typ=float, ) - results["value"][key][patch] = c + results["value"][key][campaign] = c # Additive bias (weighted error of mean) key_base = "dmcw_" @@ -368,16 +368,16 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): for comp in (1, 2): key = f"{key_base}{comp}" init_key(results, key, "var_m", extra=["N_gal", f"cw_{comp}"]) - for patch in stats_files: + for campaign in stats_files: dmc = get_match( stats_files, - patch, + campaign, rf"{key_base}{comp} = (\S+)", typ=float, previous=[f"^{shape}:$"], n_previous=[2 * comp + 8], ) - results["value"][key][patch] = dmc + results["value"][key][campaign] = dmc # Additive bias (jackknife) key_base = "cjk_" @@ -393,26 +393,26 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): key_s = key init_key(results, key, "var_m", extra=["w_tot", f"cjk_{comp}"]) - for patch in stats_files: + for campaign in stats_files: c, dc = get_match( stats_files, - patch, + campaign, rf"{key_base}{comp} = (\S+)", previous=[f"^{shape}:$"], n_previous=[comp], typ="ufloat", ) - results["value"][key_m][patch] = c - results["value"][key_s][patch] = dc + results["value"][key_m][campaign] = c + results["value"][key_s][campaign] = dc # Ellipticity dispersion key = "sigma2_epsilon" if use_keys[key]: init_key(results, key, "w_avg", extra="N_gal") - for patch in stats_files: - results["value"][key][patch] = get_match( + for campaign in stats_files: + results["value"][key][campaign] = get_match( stats_files, - patch, + campaign, r"Dispersion of complex ellipticity = (\S+)", previous=[f"^{shape}$"], typ=float, @@ -424,60 +424,60 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): keys = [f"{key_base}11", f"{key_base}12", f"{key_base}21", f"{key_base}22"] for key in keys: init_key(results, key, "w_avg", extra="N_gal") - for patch in stats_files: + for campaign in stats_files: tmp = get_match( stats_files, - patch, + campaign, r"\[\[(\s?\S+)\s+\S+]", previous=["ngmix galaxies:", "total response matrix:"], n_previous=[2, 1], typ=float, ) - results["value"][keys[0]][patch] = tmp + results["value"][keys[0]][campaign] = tmp tmp = get_match( stats_files, - patch, + campaign, r"\[\[\s?\S+\s+(\S+)]", previous=["ngmix galaxies:", "total response matrix:"], n_previous=[2, 1], typ=float, ) - results["value"][keys[1]][patch] = tmp + results["value"][keys[1]][campaign] = tmp tmp = get_match( stats_files, - patch, + campaign, r"\[(\s?\S+)\s+\S+\]\]", previous=["ngmix galaxies", "total response matrix:"], n_previous=[3, 2], typ=float, ) - results["value"][keys[2]][patch] = tmp + results["value"][keys[2]][campaign] = tmp tmp = get_match( stats_files, - patch, + campaign, r" \[\s?\S+\s+(\S+)\]\]", previous=["ngmix galaxies:", "total response matrix:"], n_previous=[3, 2], typ=float, ) - results["value"][keys[3]][patch] = tmp + results["value"][keys[3]][campaign] = tmp # Normalised trace = mean diagonal key_der = "trN_R_tot" init_key(results, key_der, "w_avg", extra="N_gal") - for patch in stats_files: - results["value"][key_der][patch] = ( - results["value"]["R_tot_11"][patch] - + results["value"]["R_tot_22"][patch] + for campaign in stats_files: + results["value"][key_der][campaign] = ( + results["value"]["R_tot_11"][campaign] + + results["value"]["R_tot_22"][campaign] ) / 2 # Sum of absolute off-diagonal key_der = "abs_off_R_tot" init_key(results, key_der, "w_avg", extra="N_gal") - for patch in stats_files: - results["value"][key_der][patch] = np.abs( - results["value"]["R_tot_12"][patch] - ) + np.abs(results["value"]["R_tot_21"][patch]) + for campaign in stats_files: + results["value"][key_der][campaign] = np.abs( + results["value"]["R_tot_12"][campaign] + ) + np.abs(results["value"]["R_tot_21"][campaign]) # Galaxy shear response matrix key_base = "R_shear_" @@ -485,60 +485,60 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): keys = [f"{key_base}11", f"{key_base}12", f"{key_base}21", f"{key_base}22"] for key in keys: init_key(results, key, "w_avg", extra="N_gal") - for patch in stats_files: + for campaign in stats_files: tmp = get_match( stats_files, - patch, + campaign, r"\[\[(\s?\S+)\s+\S+]", previous=["ngmix galaxies:", "shear response matrix:"], n_previous=[5, 1], typ=float, ) - results["value"][keys[0]][patch] = tmp + results["value"][keys[0]][campaign] = tmp tmp = get_match( stats_files, - patch, + campaign, r"\[\[\s?\S+\s+(\S+)]", previous=["ngmix galaxies:", "shear response matrix:"], n_previous=[5, 1], typ=float, ) - results["value"][keys[1]][patch] = tmp + results["value"][keys[1]][campaign] = tmp tmp = get_match( stats_files, - patch, + campaign, r"\[(\s?\S+)\s+\S+\]\]", previous=["ngmix galaxies", "shear response matrix:"], n_previous=[6, 2], typ=float, ) - results["value"][keys[2]][patch] = tmp + results["value"][keys[2]][campaign] = tmp tmp = get_match( stats_files, - patch, + campaign, r" \[\s?\S+\s+(\S+)\]\]", previous=["ngmix galaxies:", "shear response matrix:"], n_previous=[6, 2], typ=float, ) - results["value"][keys[3]][patch] = tmp + results["value"][keys[3]][campaign] = tmp # Normalised trace = mean diagonal key_der = "trN_R_shear" init_key(results, key_der, "w_avg", extra="N_gal") - for patch in stats_files: - results["value"][key_der][patch] = ( - results["value"]["R_shear_11"][patch] - + results["value"]["R_shear_22"][patch] + for campaign in stats_files: + results["value"][key_der][campaign] = ( + results["value"]["R_shear_11"][campaign] + + results["value"]["R_shear_22"][campaign] ) / 2 # Sum of absolute off-diagonal key_der = "abs_off_R_shear" init_key(results, key_der, "w_avg", extra="N_gal") - for patch in stats_files: - results["value"][key_der][patch] = np.abs( - results["value"]["R_shear_12"][patch] - ) + np.abs(results["value"]["R_shear_21"][patch]) + for campaign in stats_files: + results["value"][key_der][campaign] = np.abs( + results["value"]["R_shear_12"][campaign] + ) + np.abs(results["value"]["R_shear_21"][campaign]) # Galaxy selection response matrix key_base = "R_select_" @@ -546,60 +546,60 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): keys = [f"{key_base}11", f"{key_base}12", f"{key_base}21", f"{key_base}22"] for key in keys: init_key(results, key, "w_avg", extra="N_gal") - for patch in stats_files: + for campaign in stats_files: tmp = get_match( stats_files, - patch, + campaign, r"\[\[(\s?\S+)\s+\S+]", previous=["ngmix galaxies:", "selection response matrix:"], n_previous=[8, 1], typ=float, ) - results["value"][keys[0]][patch] = tmp + results["value"][keys[0]][campaign] = tmp tmp = get_match( stats_files, - patch, + campaign, r"\[\[\s?\S+\s+(\S+)]", previous=["ngmix galaxies:", "selection response matrix:"], n_previous=[8, 1], typ=float, ) - results["value"][keys[1]][patch] = tmp + results["value"][keys[1]][campaign] = tmp tmp = get_match( stats_files, - patch, + campaign, r"\[(\s?\S+)\s+\S+\]\]", previous=["ngmix galaxies", "selection response matrix:"], n_previous=[9, 2], typ=float, ) - results["value"][keys[2]][patch] = tmp + results["value"][keys[2]][campaign] = tmp tmp = get_match( stats_files, - patch, + campaign, r" \[\s?\S+\s+(\S+)\]\]", previous=["ngmix galaxies:", "selection response matrix:"], n_previous=[9, 2], typ=float, ) - results["value"][keys[3]][patch] = tmp + results["value"][keys[3]][campaign] = tmp # Normalised trace = mean diagonal key_der = "trN_R_select" init_key(results, key_der, "w_avg", extra="N_gal") - for patch in stats_files: - results["value"][key_der][patch] = ( - results["value"]["R_select_11"][patch] - + results["value"]["R_select_22"][patch] + for campaign in stats_files: + results["value"][key_der][campaign] = ( + results["value"]["R_select_11"][campaign] + + results["value"]["R_select_22"][campaign] ) / 2 # Sum of absolute off-diagonal key_der = "abs_off_R_select" init_key(results, key_der, "w_avg", extra="N_gal") - for patch in stats_files: - results["value"][key_der][patch] = np.abs( - results["value"]["R_select_12"][patch] - ) + np.abs(results["value"]["R_select_21"][patch]) + for campaign in stats_files: + results["value"][key_der][campaign] = np.abs( + results["value"]["R_select_12"][campaign] + ) + np.abs(results["value"]["R_select_21"][campaign]) # Object-wise PSF leakage key_base = "m_" @@ -614,70 +614,70 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): ] for key in keys: init_key(results, key, "w_avg", extra="N_gal") - for patch in stats_files: + for campaign in stats_files: m, dm = get_match( stats_files, - patch, + campaign, "\\$e_\\{1\\}\\^\\{\\\\rm PSF\\}\\$: m_1=(\\S*)", previous=["ngmix"], n_previous=[1], typ="ufloat", ) - results["value"]["m_11"][patch] = m + results["value"]["m_11"][campaign] = m m, dm = get_match( stats_files, - patch, + campaign, "\\$e_\\{1\\}\\^\\{\\\\rm PSF\\}\\$: m_2=(\\S*)", previous=["ngmix"], n_previous=[2], typ="ufloat", ) - results["value"]["m_12"][patch] = m + results["value"]["m_12"][campaign] = m m, dm = get_match( stats_files, - patch, + campaign, "\\$e_\\{2\\}\\^\\{\\\\rm PSF\\}\\$: m_1=(\\S*)", previous=["ngmix"], n_previous=[3], typ="ufloat", ) - results["value"]["m_21"][patch] = m + results["value"]["m_21"][campaign] = m m, dm = get_match( stats_files, - patch, + campaign, "\\$e_\\{2\\}\\^\\{\\\\rm PSF\\}\\$: m_2=(\\S*)", previous=["ngmix"], n_previous=[4], typ="ufloat", ) - results["value"]["m_22"][patch] = m + results["value"]["m_22"][campaign] = m m, dm = get_match( stats_files, - patch, + campaign, "\\$\\\\mathrm\\{FWHM\\}\\^\\{\\\\rm PSF\\}\\$ \\[arcsec]: m_1=(\\S+)", previous=["ngmix"], n_previous=[5], typ="ufloat", ) - results["value"]["m_s1"][patch] = m + results["value"]["m_s1"][campaign] = m m, dm = get_match( stats_files, - patch, + campaign, "\\$\\\\mathrm\\{FWHM\\}\\^\\{\\\\rm PSF\\}\\$ \\[arcsec]: m_2=(\\S+)", previous=["ngmix"], n_previous=[6], typ="ufloat", ) - results["value"]["m_s2"][patch] = m + results["value"]["m_s2"][campaign] = m # Scale-dependent PSF leakage key = "alpha" if use_keys[key]: init_key(results, key, "w_avg", extra="N_gal") - for patch in stats_files: - results["value"][key][patch] = get_match( + for campaign in stats_files: + results["value"][key][campaign] = get_match( stats_files, - patch, + campaign, r"ngmix: Weighted average alpha =(\s?\S+)", typ=float, ) @@ -687,15 +687,15 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): if use_keys[key_base]: init_key(results, "xi_sys_p", "w_avg", extra="N_gal") init_key(results, "xi_sys_m", "w_avg", extra="N_gal") - for patch in stats_files: + for campaign in stats_files: tmp = get_match( - stats_files, patch, r"ngmix: <\|xi_sys_\+\|> = (\S*)", typ=float + stats_files, campaign, r"ngmix: <\|xi_sys_\+\|> = (\S*)", typ=float ) - results["value"]["xi_sys_p"][patch] = tmp + results["value"]["xi_sys_p"][campaign] = tmp tmp = get_match( - stats_files, patch, r"ngmix: <\|xi_sys_\-\|> = (\S*)", typ=float + stats_files, campaign, r"ngmix: <\|xi_sys_\-\|> = (\S*)", typ=float ) - results["value"]["xi_sys_m"][patch] = tmp + results["value"]["xi_sys_m"][campaign] = tmp def latex_table(file_base, cols=None, col_names=None): @@ -714,7 +714,7 @@ def latex_table(file_base, cols=None, col_names=None): print(r"}\hline\hline", file=fout) # Table header - str_line = "patch\t&" + str_line = "campaign\t&" for name in col_names: str_line = f"{str_line} ${name}$\t&" # str_line = f'{str_line} \\multicolumn{{2}}{{c}}{{${name}$}}\t&' @@ -724,7 +724,7 @@ def latex_table(file_base, cols=None, col_names=None): for nl in range(n_lines): str_line = "" - str_line = f"{str_line}{dat['patch'][nl]}\t&" + str_line = f"{str_line}{dat['campaign'][nl]}\t&" for col in cols: if len(col) == 2: @@ -767,25 +767,17 @@ def main(argv=None): if argv[1] == "snr": # All directories - patches = [f.path for f in os.scandir(".") if f.is_dir()] + campaigns = [f.path for f in os.scandir(".") if f.is_dir()] all = False - elif argv[1] == "v1": - n_patch = 7 - patches = [f"P{x}" for x in np.arange(n_patch) + 1] - elif argv[1] == "v1.5": - n_patch = 8 - patches = [f"P{x}" for x in np.arange(n_patch) + 1] - elif argv[1] == "test": - patches = ["P7", "W3", "S4"] elif argv[1] == "comb": # Validate with combined catalogue - patches = ["comb"] + campaigns = ["comb"] else: - patches = argv[1].split("+") + campaigns = argv[1].split("+") - n_patch = len(patches) + n_campaign = len(campaigns) - print("combine_results.py:", patches) + print("combine_results.py:", campaigns) directory = "sp_output/plots" fbase = "stats_file" @@ -796,7 +788,7 @@ def main(argv=None): verbose = False - stats_files = read_stats_files(patches, path, verbose=verbose) + stats_files = read_stats_files(campaigns, path, verbose=verbose) results = {"value": {}, "type": {}, "extra": {}, "all": {}} diff --git a/scripts/compute_area.py b/scripts/compute_area.py index afe797ad..1e9bfa50 100755 --- a/scripts/compute_area.py +++ b/scripts/compute_area.py @@ -21,7 +21,7 @@ def main(argv=None): random_log_path = "output/run_sp_Rc/random_cat_runner/logs" log_file_base = "process" - # Get expected number of tiles in patch + # Get expected number of tiles in campaign num_lines = sum(1 for _ in open(tile_ID_path)) print(f"Found {num_lines} tiles in ID file {tile_ID_path}") @@ -108,9 +108,9 @@ def main(argv=None): area_deg2_non_overl_tile = ufloat( np.mean(area_deg2_non_overl), np.std(area_deg2_non_overl) ) - print(f"Patch area without overlap = {area_deg2_non_overl_total:.3f} deg^2") + print(f"Campaign area without overlap = {area_deg2_non_overl_total:.3f} deg^2") print( - f"Patch area without overlap and no 0 gal = {area_deg2_non_overl_total_wgal:.3f} deg^2" + f"Campaign area without overlap and no 0 gal = {area_deg2_non_overl_total_wgal:.3f} deg^2" ) print(f"Tile area without overlap = {area_deg2_non_overl_tile:.3fP} deg^2") @@ -120,10 +120,10 @@ def main(argv=None): np.mean(area_deg2_eff_non_overl), np.std(area_deg2_eff_non_overl) ) print( - f"Unmasked patch area without overlap = {area_deg2_eff_non_overl_total:.3f} deg^2" + f"Unmasked campaign area without overlap = {area_deg2_eff_non_overl_total:.3f} deg^2" ) print( - f"Unmasked patch area without overlap and no 0 gal = {area_deg2_eff_non_overl_total_wgal:.3f} deg^2" + f"Unmasked campaign area without overlap and no 0 gal = {area_deg2_eff_non_overl_total_wgal:.3f} deg^2" ) print( f"Unmasked tile area without overlap = {area_deg2_eff_non_overl_tile:.3fP} deg^2" diff --git a/scripts/merge_psf_cat.py b/scripts/merge_psf_cat.py index 65c1e891..adda9d22 100644 --- a/scripts/merge_psf_cat.py +++ b/scripts/merge_psf_cat.py @@ -1,7 +1,7 @@ #!/usr/bin/env python3 """MERGE PSF CAT. -Merge PSF catalogues (psf_catalog_ngmix.fits) from different patches +Merge PSF catalogues (psf_catalog_ngmix.fits) from different campaigns into a single FITS file. :Author: Martin Kilbinger @@ -18,7 +18,7 @@ class MergePsfCat: """Merge Psf Cat. - Class to merge PSF catalogues from multiple patches. + Class to merge PSF catalogues from multiple campaigns. """ @@ -32,7 +32,7 @@ def params_default(self): """ self._params = { - "patches": "v1", + "campaigns": None, "sh": "ngmix", "survey": "unions", "year": "2024", @@ -43,7 +43,7 @@ def params_default(self): "verbose": False, } self._short_options = { - "patches": "-p", + "campaigns": "-p", "sh": "-g", "survey": "-s", "year": "-y", @@ -55,15 +55,12 @@ def params_default(self): "hdu": "int", } self._help_strings = { - "patches": ( - "list of patches separated by '+', or shortcut " - "(allowed are 'v1', 'v1.5'), default={}" - ), + "campaigns": "list of campaigns separated by '+'", "sh": "shape measurement method, default={}", "survey": "survey name, default={}", "year": "year of processing, default={}", "version": "catalogue version, default={}", - "base_path": "base path containing patch directories, default={}", + "base_path": "base path containing campaign directories, default={}", "hdu": "HDU number to read from input FITS files, default={}", } @@ -85,40 +82,32 @@ def set_params_from_command_line(self, args): logging.log_command(args) - def get_patches(self): - """Get Patches. + def get_campaigns(self): + """Get Campaigns. - Return list of patches according to option parameter value. + Return list of campaigns according to option parameter value. Returns ------- list - patches, list of str + campaigns, list of str """ - if self._params["patches"] == "v1": - n_patch = 7 - patches = [f"P{x}" for x in np.arange(n_patch) + 1] - elif self._params["patches"] == "v1.5": - n_patch = 8 - patches = [f"P{x}" for x in np.arange(n_patch) + 1] - elif self._params["patches"] == "v1.6": - n_patch = 9 - patches = [f"P{x}" for x in np.arange(n_patch) + 1] - else: - patches = self._params["patches"].split("+") - - return patches - - def merge_catalogues(self, patches): + campaigns = self._params["campaigns"] + if not campaigns: + raise ValueError("No campaigns given; set 'campaigns'") + + return campaigns.split("+") + + def merge_catalogues(self, campaigns): """Merge Catalogues. - Merge PSF catalogues from sub-patches into one FITS file. + Merge PSF catalogues from campaigns into one FITS file. Parameters ---------- - patches : list of str - list of patches/sub-directories + campaigns : list of str + list of campaigns / sub-directories """ base_path = self._params["base_path"] @@ -133,11 +122,11 @@ def merge_catalogues(self, patches): ) dat_all = {} - for idx, patch in enumerate(patches): + for idx, campaign in enumerate(campaigns): if verbose: - print(f" {patch}") + print(f" {campaign}") - input_path = f"{base_path}/{patch}/{input_sub_path}" + input_path = f"{base_path}/{campaign}/{input_sub_path}" try: dat = fits.getdata(input_path, hdu_in) except Exception: @@ -149,18 +138,18 @@ def merge_catalogues(self, patches): col_names = dat.dtype.names for name in col_names: dat_all[name] = [] - dat_all["patch"] = [] + dat_all["campaign"] = [] for name in col_names: dat_all[name] = np.append(dat_all[name], dat[name]) - dat_all["patch"] = np.append(dat_all["patch"], [idx + 1] * len(dat)) + dat_all["campaign"] = np.append(dat_all["campaign"], [idx + 1] * len(dat)) - col_names = col_names + ("patch",) + col_names = col_names + ("campaign",) column_all = [] for name in col_names: - if name != "patch": + if name != "campaign": my_format = "D" else: my_format = "I" @@ -177,12 +166,12 @@ def run(self): Main processing function. """ - patches = self.get_patches() + campaigns = self.get_campaigns() if self._params["verbose"]: - print("Merging PSF catalogues from patches:", patches) + print("Merging PSF catalogues from campaigns:", campaigns) - self.merge_catalogues(patches) + self.merge_catalogues(campaigns) def main(argv=None): diff --git a/scripts/plot_rho_stats_patches.py b/scripts/plot_rho_stats_patches.py deleted file mode 100755 index 915c07d1..00000000 --- a/scripts/plot_rho_stats_patches.py +++ /dev/null @@ -1,98 +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: python3 -# --- - -import glob -import itertools -import os - -import matplotlib.pylab as plt -from cs_util import args -from shear_psf_leakage.rho_tau_stat import RhoStat - - -# + -# Set parameters from file or user input -class dummy(object): - def __init__(self): - - self._params = { - "in_dir_base": ".", - "title": None, - } - - -obj = dummy() -params_upd = args.read_param_script("params_rho.py", obj._params, verbose=True) -for key in params_upd: - obj._params[key] = params_upd[key] - -# patches = [f'P{x}' for x in np.arange(n_patch) + 1] -patches = glob.glob("P*") -print(patches) - -default_colors = plt.rcParams["axes.prop_cycle"].by_key()["color"] -color_cycle = itertools.cycle(default_colors) -col = {} -for patch in patches: - col[patch] = next(color_cycle) - -coord_units = "deg" -theta_min = 0.1 -theta_max = 250 -sep_units = "arcmin" -nbins = 20 - -# ## Set up -TreeCorrConfig = { - "ra_units": coord_units, - "dec_units": coord_units, - "min_sep": theta_min, - "max_sep": theta_max, - "sep_units": sep_units, - "nbins": nbins, - "var_method": "bootstrap", -} -# - - -rho_stat_handler = RhoStat( - output=obj._params["in_dir_base"], treecorr_config=TreeCorrConfig, verbose=True -) - -# + -filenames = [] -colors = [] - -for patch in patches: - path = f"{patch}/output/run_sp_Pl/mccd_plots_runner/output/rho_stats_id.fits" - if os.path.exists(f"{obj._params['in_dir_base']}/{path}"): - print(f"Reading rho stats {obj._params['in_dir_base']}/{path}...") - rho_stat_handler.load_rho_stats(path) - filenames.append(path) - colors.append(col[patch]) - else: - print( - f"File rho stats {obj._params['in_dir_base']}/{path} not found, skipping.." - ) -# - - -# Create plot -rho_stat_handler.plot_rho_stats( - filenames, - colors, - patches, - abs=False, - savefig="rho_stats.png", - legend="outside", - title=obj._params["title"], -) diff --git a/scripts/prepare_patch_for_spval.sh b/scripts/prepare_patch_for_spval.sh deleted file mode 100755 index 634d901f..00000000 --- a/scripts/prepare_patch_for_spval.sh +++ /dev/null @@ -1,31 +0,0 @@ -#!/usr/bin/env bash - -patch=$1 - -spdir=$HOME/astro/repositories/github/sp_validation - -# Galaxy catalogue -#cp ~/psfex/final_cat_${patch}.hdf5 . -ln -s ~/psfex/final_cat_${patch}.hdf5 - -# Parameter file, to avoid read errors for hdf5 file -ln -sf ~/shapepipe/example/cfis/final_cat.param - -# Star catalogue -## Ellipticities in pixel coordinates, MCCD output -# ln -sf $HOME/psfex/${patch}/output/run_sp_Ms/merge_starcat_runner/output/full_starcat-0000000.fits - -## Projected back to world coordinates -ln -sf $HOME/psfex/star_cat/${patch}/output/run_sp_Ms/merge_starcat_runner/output/full_starcat-0000000.fits - -# Tile number list -ln -sf ~/shapepipe/auxdir/CFIS/tiles_202106/tiles_${patch}.txt - -# Parameter file -#cp $spdir/scripts/calibration/params.py . -echo "Diff:" -diff $spdir/scripts/calibration/params.py params.py -echo "Run?" -echo "cp $spdir/scripts/calibration/params.py params.py" -echo "Run?" -echo "ipython $spdir/scripts/calibration/extract_info.py" diff --git a/scripts/star_match_stats.py b/scripts/star_match_stats.py deleted file mode 100644 index 5c74ed5c..00000000 --- a/scripts/star_match_stats.py +++ /dev/null @@ -1,44 +0,0 @@ -#!/usr/bin/env python - -import re -import sys - -import numpy as np - - -def main(argv=None): - - types = ["star", "gal", "other"] - text = "Number of stars selected as" - - n_patch = 7 - patches = [f"P{x}" for x in np.arange(n_patch) + 1] - - ntyp = {} - ntot = {} - for typ in types: - ntyp[typ] = 0 - ntot[typ] = 0 - - for patch in patches: - # print(patch) - path = f"{patch}/sp_output/plots/stats_file.txt" - with open(path, "r") as fin: - lines = fin.readlines() - for typ in types: - for line in lines: - pattern = rf"{text} {typ}.*= (\d+)/(\d+)" - m = re.search(pattern, line) - if m: - ntyp_patch = int(m.group(1)) - ntot_patch = int(m.group(2)) - # print(typ, m.group(1), m.group(2)) - ntyp[typ] += ntyp_patch - ntot[typ] += ntot_patch - - for typ in types: - print(f"{text} {typ} = {ntyp[typ]}/{ntot[typ]} = {ntyp[typ] / ntot[typ]:.2%}") - - -if __name__ == "__main__": - sys.exit(main(sys.argv)) diff --git a/scripts/stats_tile_id_gal_counts.py b/scripts/stats_tile_id_gal_counts.py index caed59ae..26adca7f 100755 --- a/scripts/stats_tile_id_gal_counts.py +++ b/scripts/stats_tile_id_gal_counts.py @@ -1,6 +1,7 @@ #!/usr/bin/env python3 import copy +import os import sys from optparse import OptionParser @@ -40,7 +41,7 @@ def params_default(): parameter values """ - p_def = param(survey="v1") + p_def = param(campaigns=None) return p_def @@ -66,7 +67,13 @@ def parse_options(p_def): # I/O parser.add_option("-i", "--input", dest="input", type="string", help="input file") - parser.add_option("-s", "--survey", dest="survey", type="string", help="survey") + parser.add_option( + "-c", + "--campaigns", + dest="campaigns", + type="string", + help="campaigns separated by '+'", + ) options, args = parser.parse_args() @@ -87,8 +94,8 @@ def check_options(options): Result of option check. False if invalid option value. """ - if not options.input and not options.survey: - print("Either input or survey need to be specified") + if not options.input and not options.campaigns: + print("Either input or campaigns need to be specified") return False return True @@ -148,13 +155,13 @@ def main(argv=None): # Get input file paths print("Retrieving input file paths") input_files = [] - if param.survey == "v1": - n_patch = 7 - patches = [f"P{x}" for x in np.arange(n_patch) + 1] - for patch in patches: - path = f"{patch}/sp_output/tile_id_gal_counts_{sh}.txt" + if param.campaigns: + campaigns = param.campaigns.split("+") + for campaign in campaigns: + path = f"{campaign}/sp_output/tile_id_gal_counts_{sh}.txt" input_files.append(path) else: + campaigns = [os.path.dirname(param.input) or "."] input_files = [param.input] print(f"Found {len(input_files)} input files") @@ -166,11 +173,11 @@ def main(argv=None): n_gal_arr = [] n_shape_arr = [] - for patch, input_path in zip(patches, input_files): - dat[patch] = np.loadtxt(input_path) - n_det = dat[patch][:, 1] - n_gal = dat[patch][:, 2] - n_shape = dat[patch][:, 3] + for campaign, input_path in zip(campaigns, input_files): + dat[campaign] = np.loadtxt(input_path) + n_det = dat[campaign][:, 1] + n_gal = dat[campaign][:, 2] + n_shape = dat[campaign][:, 3] n_det_arr.extend(n_det) n_gal_arr.extend(n_gal) @@ -178,11 +185,11 @@ def main(argv=None): # Write tile IDs with number of shapes > 0 print("Writing tile IDs with n_shapes>0") - for patch in patches: - out_path = f"{patch}/sp_output/found_ID_wshapes.txt" + for campaign in campaigns: + out_path = f"{campaign}/sp_output/found_ID_wshapes.txt" with open(out_path, "w") as f_out: - mask_n_shape = dat[patch][:, 3] > 0 - tile_ID_masked = dat[patch][mask_n_shape, 0] + mask_n_shape = dat[campaign][:, 3] > 0 + tile_ID_masked = dat[campaign][mask_n_shape, 0] for ID in tile_ID_masked: print(f"{ID:07.3f}", file=f_out) diff --git a/scripts/survey_stats_all.sh b/scripts/survey_stats_all.sh deleted file mode 100644 index f6c4cd5f..00000000 --- a/scripts/survey_stats_all.sh +++ /dev/null @@ -1,94 +0,0 @@ -#!/usr/bin/bash - -fbase_found='found_ID' -fbase_found_wsh='found_ID_wshapes' -fbase_lf='CFIS3500_THELI' -rm -f ${fbase_found}_all.txt -rm -f ${fbase_found_wsh}_all.txt -rm -f ${fbase_found_random}_all.txt -rm -f ${fbase_lf}_all.txt - -for patch in P1 P2 P3 P4 P5 P6 P7; do - - # Patch - - ## Total number of tiles - wc -l $patch/tiles_P?.txt - - - # Final catalogue - #echo "Final catalogue" - - ## Number of final catalogues (.tgz) - ntgz=`ls -rtl $patch/final*.tgz | wc -l` - echo "$ntgz final .tgz cats" - - ## Number of final catalogues (.fits) - nfits=`ls -rtl $patch/output/run_sp_combined/make_catalog_runner/output/final* | wc -l` - echo "$nfits final .fits cats" - - ## Number of merged final catalogues - if [ -e $patch/log_merge_final_gal_cat ]; then - tail -n 1 $patch/log_merge_final_gal_cat - fi - - ## Number of tile IDs found in merged catalogue - if [ -e $patch/sp_output/${fbase_found}.txt ]; then - wc -l $patch/sp_output/${fbase_found}.txt - wc -l $patch/sp_output/${fbase_found}.txt >> ${fbase_found}_all.txt - fi - - ## Number of tile IDs found in merged catalogue with shapes - if [ -e $patch/sp_output/${fbase_found_wsh}.txt ]; then - wc -l $patch/sp_output/${fbase_found_wsh}.txt - wc -l $patch/sp_output/${fbase_found_wsh}.txt >> ${fbase_found_wsh}_all.txt - fi - - - # Random catalogue - - ## Number of final catalogues (.tgz) - ntgz=`ls -rtl $patch/pipeline_flag*.tgz | wc -l` - echo "$ntgz random .tgz cats" - - - ## Number of random catalogues (.fits) - if [ -d $patch/output/run_sp_combined_flag ]; then - nfits=`ls -rtl $patch/output/run_sp_combined_flag/mask_runner/output/pip* | wc -l` - echo "$nfits random .fits cats" - fi - - ## Number of merged random catalogues - if [ -e $patch/log_merge_final_rand_cat ]; then - tail -n 1 $patch/log_merge_final_rand_cat - fi - - ### Number of tiles for random catalogue validation - if [ -e $patch/sp_output_random/${fbase_found}.txt ]; then - wc -l $patch/sp_output_random/${fbase_found}.txt - wc -l $patch/sp_output_random/${fbase_found}.txt >> ${fbase_found_random}_all.txt - fi - - # LensFit tile IDs - if [ -e ${fbase_lf}_$patch.list ]; then - wc -l ${fbase_lf}_$patch.list - wc -l ${fbase_lf}_$patch.list >> ${fbase_lf}_all.txt - fi - - echo - -done - -echo -n "number of tiles in ${fbase_found}_all.txt = " -summe.pl ${fbase_found}_all.txt 0 - -echo -n "number of tiles in ${fbase_found_wsh}_all.txt = " -summe.pl ${fbase_found_wsh}_all.txt 0 - -echo -n "number of tiles in ${fbase_found_random}_all.txt = " -summe.pl ${fbase_found_random}_all.txt 0 - -if [ -e ${fbase_lf} ]; then - echo -n "number of tiles in ${fbase_lf}_all.txt = " - summe.pl ${fbase_lf}_all.txt 0 -fi diff --git a/src/sp_validation/catalog.py b/src/sp_validation/catalog.py index 2c4ebadf..52e850e4 100644 --- a/src/sp_validation/catalog.py +++ b/src/sp_validation/catalog.py @@ -23,7 +23,6 @@ from cs_util import cat from sp_validation import format, io -from sp_validation.survey import get_footprint from sp_validation.version import __version__ @@ -154,7 +153,6 @@ def check_matching( keys_2, thresh, stats_file, - name=None, verbose=False, ): """Check matching. @@ -182,13 +180,7 @@ def check_matching( index list of tiles in footprint """ - if name is not None: - # Filter stars outside footprint for efficiency - mask_area_tiles = get_footprint(name, d1[keys_1[0]], d1[keys_1[1]]) - if len(np.where(mask_area_tiles)[0]) == 0: - raise ValueError(f"Error: no object found in field '{name}'") - else: - mask_area_tiles = np.arange(len(d1)) + mask_area_tiles = np.arange(len(d1)) # Match stars from exposure (PSF) catalogue to total catalogue ind = match_stars2( diff --git a/src/sp_validation/catalog_builders.py b/src/sp_validation/catalog_builders.py index eaea4fb5..00b00694 100644 --- a/src/sp_validation/catalog_builders.py +++ b/src/sp_validation/catalog_builders.py @@ -279,147 +279,84 @@ def params_default(self): """ self._params = { - "patches": "v1", + "input_paths": None, "sh": "ngmix", "survey": "unions", "year": "2024", "version": "1.4.2", "pipeline": "shapepipe", - "hdu": 1, + "param_path": None, "reduce_mem": False, "verbose": False, } self._short_options = { - "patches": "-p", + "input_paths": "-i", "sh": "-g", "survey": "-s", "year": "-y", "version": "-V", + "param_path": "-p", "reduce_mem": "-r", } self._types = { - "hdu": "int", "reduce_mem": "bool", } self._help_strings = { - "patches": "list of patches separated by '+', or shortcut (allowed are 'v1'), default={}", + "input_paths": ( + "campaign catalogue files (final_cat_.hdf5) to merge," + + " separated by '+'" + ), "sh": "shape measurement method, default={}", "survey": "survey name, default={}", "year": "year of processing, default={}", "version": "catalogue version, default={}", + "param_path": "path to parameter file listing columns to keep", "reduce_mem": "output some columns in lower precision to reduce memory", } - def get_patches(self): - """Get Patches. + def get_input_paths(self): + """Get Input Paths. - Return list of patches according to option parameter value. + Return the list of campaign catalogue files to merge. Returns ------- - list - patches, list of str + list of str + input file paths """ - if self._params["patches"] == "v1": - n_patch = 7 - patches = [f"P{x}" for x in np.arange(n_patch) + 1] - elif self._params["patches"] == "v1.5": - n_patch = 8 - patches = [f"P{x}" for x in np.arange(n_patch) + 1] - - else: - patches = self._params["patches"].split("+") - - return patches - - def get_n_obj(self, patches, base_path, input_sub_path): - """Get N Obj. - - Get number of objects from FITS file headers. - - Parameters - ---------- - patches : list - input patches, type is str - base_path : str - input base directory, root dir of patches - input_sub_path : str - input file name; input path is base_path/patch/input_sub_path - - Raises: - ValueError: if input file canont be read - - Returns: - list - HDUs - list - number of objects per file - int - total number of objects - - """ - if self._params["verbose"]: - print("Getting number of objects") - n_obj_list = [] - n_obj = 0 - hdu_lists = [] - for patch in patches: - input_path = f"{base_path}/{patch}/{input_sub_path}" - try: - hdu_list = fits.open(input_path) - except Exception as err: - raise ValueError( - f"Could not open file {input_path} at HDU" - + f" #{self._params['hdu']}" - ) from err - hdu_lists.append(hdu_list) - - this_n = int(hdu_list[self._params["hdu"]].header["NAXIS2"]) - n_obj_list.append(this_n) - n_obj += this_n - - if self._params["verbose"]: - print(f"Found a total of {n_obj} (~{format.millify(n_obj)}) objects.") + input_paths = self._params["input_paths"] + if not input_paths: + raise ValueError( + "No input campaign catalogues given; set 'input_paths' to one" + + " or more final_cat_.hdf5 files separated by '+'" + ) + if isinstance(input_paths, str): + input_paths = input_paths.split("+") - return hdu_lists, n_obj_list, n_obj + return [path.strip() for path in input_paths if path.strip()] - def get_col_info(self, dat): - """Get Col Info. + @staticmethod + def campaign_name(input_path): + """Campaign Name. - Return information of input columns. + Return the campaign name encoded in a catalogue file name, + ``final_cat_.hdf5`` -> ````. Parameters ---------- - dat : numpy.ndarray - input data + input_path : str + input file path Returns ------- - list - column names - list - column formats - int - number of columns + str + campaign name """ - col_names = dat.dtype.names - - n_col = 0 - formats = {} - ndim = {} - for name in col_names: - formats[name] = dat.dtype.fields[name][0] - ndim[name] = dat[name].ndim - n_col += ndim[name] - # Add one for patch - n_col += 1 - - if self._params["verbose"]: - print(f"Number of input (output) columns = {len(col_names)} ({n_col})") - - return col_names, formats, ndim, n_col + stem = os.path.splitext(os.path.basename(input_path))[0] + prefix = "final_cat_" + return stem[len(prefix) :] if stem.startswith(prefix) else stem def dtype_out(self, name, dtype_in): """Set output dtype. @@ -442,6 +379,7 @@ def dtype_out(self, name, dtype_in): cols_keep_dtype = [ "RA", "Dec", + "DEC", "FLAGS", "NUMBER", ] @@ -449,9 +387,9 @@ def dtype_out(self, name, dtype_in): # Transform unicode to string of equal length return np.dtype(f"S{dtype_in.itemsize // 4}") - if self._params["reduce_mem"] == False: + if not self._params["reduce_mem"]: return dtype_in - elif name not in cols_keep_dtype: + if name not in cols_keep_dtype: if dtype_in.kind == "f" and dtype_in.itemsize == 8: return np.float32 if dtype_in.kind == "i" and dtype_in.itemsize == 4: @@ -459,59 +397,33 @@ def dtype_out(self, name, dtype_in): return dtype_in - def init_data(self, n_col, n_obj, ndim, dat): - """Init Data. + def output_dtype(self, dtype_in, n_char_campaign): + """Output Dtype. - Initialize empty structured data. + Return the merged-catalogue dtype: the input columns (possibly + reduced in precision) plus a ``campaign`` column. Parameters ---------- - n_col : int - number of columns - n_obj : int - number of objects (rows) - ndim : dict - dimension of input columns - dat : numpy.ndarray - example data + dtype_in : numpy.dtype + structured dtype of an input campaign catalogue + n_char_campaign : int + width of the campaign name column Returns ------- - numpy.ndarray - combined structure data, (n_col x n_obj) array + numpy.dtype + output structured dtype """ - # Create dtypes from input column names and types. - # Reduce memory if flag set. - # Transform multi-D columns into 1D columns - dtype_tmp_list = [] - for name in ndim: - if ndim[name] == 1: - dtype_tmp_list.append((name, self.dtype_out(name, dat[name].dtype))) - else: - for jdx in range(ndim[name]): - dtype_tmp_list.append( - (f"{name}_{jdx}", self.dtype_out(name, dat[name].dtype)) - ) - dtype_tmp_list.append(("patch", np.int8)) - dtype_tmp_struct = np.dtype(dtype_tmp_list) - - if self._params["verbose"]: - memory = n_obj * dtype_tmp_struct.itemsize - print( - f"Allocating <= {memory / 1024**3:.1f}" - + f" Gb memory for the ({n_col} x {n_obj}) input data array ...", - end="", - ) - - dat_all = np.empty((n_obj,), dtype=dtype_tmp_struct) + fields = [ + (name, self.dtype_out(name, dtype_in[name])) for name in dtype_in.names + ] + fields.append(("campaign", np.dtype(f"S{n_char_campaign}"))) - if self._params["verbose"]: - print("done") + return np.dtype(fields) - return dat_all - - def write_hdf5_file(self, dat_all, patches): + def write_hdf5_file(self, dat_all, campaigns=None): """Write HDF5 File. Write data to HDF5 file. @@ -520,8 +432,8 @@ def write_hdf5_file(self, dat_all, patches): ---------- dat_all : numpy.ndarray input data - patches : list - input patches, list of str + campaigns : list, optional + input campaign names, list of str """ output_path = ( @@ -531,12 +443,12 @@ def write_hdf5_file(self, dat_all, patches): ) with h5py.File(output_path, "w") as f: - self.write_hdf5_header(f) + self.write_hdf5_header(f, campaigns=campaigns) dset = f.create_dataset("data", data=dat_all) dset[:] = dat_all - def write_hdf5_header(self, hd5file, patches=None): + def write_hdf5_header(self, hd5file, campaigns=None): """Write HDF5 Header. Write header information to HDF5 file. @@ -545,92 +457,82 @@ def write_hdf5_header(self, hd5file, patches=None): ---------- hd5file : h5py.File input HDF5 file - patches : list, optional - input patches, list of str, default is ``None`` + campaigns : list, optional + input campaign names, list of str, default is ``None`` """ super().write_hdf5_header(hd5file) - if patches is not None: - patches_str = " ".join(patches) - hd5file.attrs["patches"] = patches_str + if campaigns is not None: + hd5file.attrs["campaigns"] = " ".join(campaigns) - def merge_catalogues(self, patches, base_path="."): + def merge_catalogues(self, input_paths): """Merge Catalogues. - Merge individual patch-based catalogues. + Merge a list of campaign catalogues into one joint catalogue, adding + a ``campaign`` column that records each object's origin. Parameters ---------- - patches : list - input patches; list of `str` - base_path : str, optional - input base directory path; default is "." + input_paths : list of str + campaign catalogue files (final_cat_.hdf5) + + Returns + ------- + numpy.ndarray + merged catalogue """ - input_sub_path = ( - f"sp_output/shape_catalog_comprehensive_{self._params['sh']}.fits" + param_list = ( + sp_cat.read_param_file( + self._params["param_path"], verbose=self._params["verbose"] + ) + if self._params["param_path"] + else None ) - # Get input FITS files - hdu_lists, n_obj_list, n_obj = self.get_n_obj( - patches, - base_path, - input_sub_path, - ) + campaigns = [self.campaign_name(path) for path in input_paths] + n_char_campaign = max(len(name) for name in campaigns) - # Read data - start = end = 0 - for idx, patch in enumerate(patches): - input_path = f"{base_path}/{patch}/{input_sub_path}" - try: - dat = fits.getdata(input_path, self._params["hdu"]) - # dat = hdu_lists[idx][self._params["hdu"]].data + data_list = [] + dtype_out = None + for input_path, campaign in zip(input_paths, campaigns): + dat = sp_cat.read_campaign_catalogue( + input_path, + param_list=param_list, + verbose=self._params["verbose"], + ) - hdu_lists[idx].close() - except Exception as err: - raise ValueError( - f"Could not read data of file {input_path} at HDU" - + f" #{self._params['hdu']}" - ) from err - - # Create empty lists if first patch - if idx == 0: - col_names, formats, ndim, n_col = self.get_col_info(dat) - dat_all = self.init_data(n_col, n_obj, ndim, dat) - - # Append new data for that patch (between start and end) - end += n_obj_list[idx] - - # Copy data - i_col = 0 - names_out = dat_all.dtype.names - for name in col_names: - if ndim[name] == 1: - # Copy 1D column - dat_all[names_out[i_col]][start:end] = dat[name] - else: - # Copy all components of multi-D column - for jdx in range(ndim[name]): - dat_all[names_out[i_col + jdx]][start:end] = dat[name][:, jdx] - i_col += ndim[name] - # Add patch number - dat_all["patch"][start:end] = patch[1:] - - if i_col + 1 != n_col: + if dtype_out is None: + dtype_out = self.output_dtype(dat.dtype, n_char_campaign) + elif set(dat.dtype.names) != set(dtype_out.names) - {"campaign"}: raise ValueError( - "Inconsistent number of columns, {i_col + 1}" + f" != {n_col}" + f"Campaign catalogue {input_path} has columns" + + f" {sorted(dat.dtype.names)}, incompatible with" + + f" {sorted(set(dtype_out.names) - {'campaign'})}" ) + + dat_out = np.empty(len(dat), dtype=dtype_out) + for name in dat.dtype.names: + dat_out[name] = dat[name] + dat_out["campaign"] = campaign.encode() + data_list.append(dat_out) + if self._params["verbose"]: print( - f"{patch}: Added {len(dat)} (~{format.millify(len(dat))})" - + f" objects (from {start} to {end - 1})." + f"{campaign}: added {len(dat)}" + + f" (~{format.millify(len(dat))}) objects." ) - start = end - del dat + dat_all = np.concatenate(data_list, axis=0) - self.write_hdf5_file(dat_all, patches) + if self._params["verbose"]: + print( + f"Merged {len(dat_all)} (~{format.millify(len(dat_all))})" + + f" objects from {len(campaigns)} campaign(s)." + ) + + return dat_all def run(self): """Run. @@ -638,11 +540,12 @@ def run(self): Main processing function. """ - patches = self.get_patches() + input_paths = self.get_input_paths() if self._params["verbose"]: - print("Merging patches", patches) + print("Merging campaigns", input_paths) - self.merge_catalogues(patches) + dat_all = self.merge_catalogues(input_paths) + self.write_hdf5_file(dat_all, [self.campaign_name(p) for p in input_paths]) class ApplyHspMasks(BaseCat): diff --git a/src/sp_validation/galaxy.py b/src/sp_validation/galaxy.py index 3eda60af..7a8b8b1c 100644 --- a/src/sp_validation/galaxy.py +++ b/src/sp_validation/galaxy.py @@ -91,6 +91,10 @@ def mask_cut(dd, mask_columns=None): """ columns = list(DEFAULT_MASK_COLUMNS if mask_columns is None else mask_columns) + if not columns: + # No masking requested (e.g. the image simulations, which run no + # imaging-flag masking stage and carry no mask columns). + return np.ones(len(dd[_column_names(dd)[0]]), dtype=bool) available = _column_names(dd) missing = [col for col in columns if col not in available] diff --git a/src/sp_validation/rho_tau.py b/src/sp_validation/rho_tau.py index 6874a08c..13be5849 100644 --- a/src/sp_validation/rho_tau.py +++ b/src/sp_validation/rho_tau.py @@ -351,7 +351,7 @@ def get_jackknife_cov( tau_chunk = outdir + f"/cov_tau_{version}{i}.npy" rho_chunk = outdir + f"/cov_rho_{version}{i}.npy" if not (os.path.exists(tau_chunk) and os.path.exists(rho_chunk)): - print(f"Computing rho-statistics for {version} (patch {i + 1}/{ncov})") + print(f"Computing rho-statistics for {version} (jackknife realisation {i + 1}/{ncov})") if f"psf_{version}{i}" not in rho_stat_handler.catalogs.catalogs_dict: # Build catalogues diff --git a/src/sp_validation/survey.py b/src/sp_validation/survey.py index fe35cde9..e9d61025 100644 --- a/src/sp_validation/survey.py +++ b/src/sp_validation/survey.py @@ -183,77 +183,6 @@ def write_tile_id_gal_counts(detection_IDs, galaxy_IDs, shape_IDs, fname): print(file=f) -def get_footprint(patch, ra, dec): - """Get Footprint. - - Return coordinates within footprint of patch. - - Parameters - ---------- - patch : str - patch name - ra : array of float - R,A, coordintates - dec : array of float - DEC coordinates - - Returns - ------- - list of float - list of coordinates withint footprint - - """ - # Set boundary coordinates between some of the patches - ra_14 = 157.5 - ra_45 = 207 - ra2_45 = 220 - ra_36 = 230 - dec_3456 = 48 - - ra2_34 = 190 - - dec_min = 29 - dec_max = 60 - - # Check whether input matches one of the seven CFIS patch name. - # Return coordinates within the patch - if patch == "P1": - return (ra > 100) & (ra < ra_14) & (dec > dec_min) & (dec < dec_max) - - elif patch == "P2": - # -30 < ra < 60 - return ((ra > 0) & (ra < 60)) | ((ra > 330) & (ra < 360)) & (dec > dec_min) & ( - dec < dec_max - ) - - elif patch == "P3": - return (ra > ra2_34) & (ra < ra_36) & (dec > dec_3456) & (dec < 70) - - elif patch == "P4": - return ( - ((ra > ra_14) & (ra < ra_45) & (dec > dec_min) & (dec < dec_3456)) - | ((ra > ra_14) & (ra < ra2_34) & (dec > dec_min) & (dec < 70)) - | ((ra > ra_45) & (ra < ra2_45) * (dec > dec_min) & (dec < 36)) - ) - - elif patch == "P5": - return ((ra > ra2_45) & (ra < 330) & (dec > dec_min) & (dec < dec_3456)) | ( - (ra > ra_45) & (ra < ra2_45) & (dec > 36) & (dec < dec_3456) - ) - - elif patch == "P6": - return (ra > ra_36) & (ra < 330) & (dec > dec_3456) & (dec > 70) - - elif patch == "P7": - return (ra > 60) & (ra < 180) & (dec > 60) & (dec < 90) - - elif patch == "W3": - return (ra > 208) & (ra < 221) & (dec > 51) & (dec < 58) - - else: - return dec > dec_min - - def area_from_coords(ra, dec, nside): """Survey area from galaxy coordinates via HEALPix pixel counting. diff --git a/src/sp_validation/tests/test_survey.py b/src/sp_validation/tests/test_survey.py index f42f884c..c31f75cc 100644 --- a/src/sp_validation/tests/test_survey.py +++ b/src/sp_validation/tests/test_survey.py @@ -28,9 +28,6 @@ def setUp(self): self._area_amin2 = 3600 self._tile_IDs = (270.283, 188.308) - self._ra = np.array([240.0]) - self._dec = np.array([32.0]) - self._patch = ["P5"] def tearDown(self): @@ -58,8 +55,3 @@ def test_get_area(self): msg=f"{tile_IDs}!={self._tile_IDs}", ) - def test_get_footprint(self): - """Test ``sp_validation.survey_get_footprint`` method.""" - for patch in self._patch: - coords = survey.get_footprint(patch, self._ra, self._dec) - self.assertTrue(coords[0]) diff --git a/workflow/image_sims/params_im_sim.py b/workflow/image_sims/params_im_sim.py index 1cd395c6..9c614e33 100644 --- a/workflow/image_sims/params_im_sim.py +++ b/workflow/image_sims/params_im_sim.py @@ -27,7 +27,7 @@ # Survey parameters -## Field or patch name -- derived from the run directory, which is named +## Field name -- derived from the run directory, which is named ## after the simulation (e.g. '1z2z_grid_1'), so one shared params file ## serves every sim. name = os.path.basename(os.getcwd()) @@ -125,10 +125,14 @@ ## Pre-calibration catalogue, including masked objects and mask flags. ## ShapePipe-v2 (post-#761) ngmix grammar: ellipticity in named scalar ## components NGMIX_G{1,2}_*, PSF size split into NGMIX_T_PSF_ORIG/RECONV. -## IMAFLAGS_ISO (present in the data-path params) is omitted: the simulation -## pipeline runs no imaging-flag masking stage, so the column does not exist. +## The MASK_n* columns (present in the data-path params) are omitted: the +## simulation pipeline runs no imaging-flag masking stage, so they do not +## exist -- hence the empty mask_columns below. ## NGMIX_MCAL_TYPES_FAIL is kept -- it is the metacal moments-failure flag the ## calibration mask cuts on, identically to the data path. +## No mask columns to OR: no masking stage in the simulation pipeline +mask_columns = [] + add_cols_pre_cal = [ "TILE_ID", "NUMBER", From d053ce8856d9555920080c2dda199a5bc38452ae Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 20:25:50 -0400 Subject: [PATCH 04/17] tests: campaign readers, mask cut, and campaign merge Seventeen unit tests on tiny synthetic hdf5 fixtures built in a temp dir: both campaign layouts (legacy patches// and flat tiles/) read identically, param-list restriction, missing-column and ambiguous-layout errors; the star reader on exposures/ hdf5 and on FITS; galaxy.mask_cut defaults, explicit list, empty list, missing column, and a v1 IMAFLAGS_ISO-only catalogue; JointCat.merge_catalogues across two campaigns of different layouts, plus its error paths. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- .../tests/test_campaign_readers.py | 259 ++++++++++++++++++ 1 file changed, 259 insertions(+) create mode 100644 src/sp_validation/tests/test_campaign_readers.py diff --git a/src/sp_validation/tests/test_campaign_readers.py b/src/sp_validation/tests/test_campaign_readers.py new file mode 100644 index 00000000..ccba4fe5 --- /dev/null +++ b/src/sp_validation/tests/test_campaign_readers.py @@ -0,0 +1,259 @@ +"""Tests for the ShapePipe v2 campaign catalogue readers and mask cut.""" + +import unittest + +import h5py +import numpy as np +import numpy.testing as npt + +from sp_validation import catalog, galaxy +from sp_validation.catalog_builders import JointCat + +import tempfile +from pathlib import Path + +GAL_DTYPE = np.dtype( + [ + ("RA", "f8"), + ("Dec", "f8"), + ("MAG_AUTO", "f4"), + ("MASK_n4", "?"), + ("MASK_n1", "?"), + ("MASK_n2", "?"), + ("MASK_n8", "?"), + ("MASK_n1024", "?"), + ] +) + + +def make_galaxy_data(n_obj, offset=0): + dat = np.zeros(n_obj, dtype=GAL_DTYPE) + dat["RA"] = np.arange(n_obj) + offset + dat["Dec"] = np.arange(n_obj) + offset + 0.5 + dat["MAG_AUTO"] = 22.0 + return dat + + +def write_campaign(path, layout, tiles): + """Write a campaign hdf5 file in the legacy or flat layout.""" + with h5py.File(path, "w") as f: + if layout == "legacy": + group = f.create_group("patches").create_group("CAMPAIGN") + else: + group = f.create_group("tiles") + for tile_id, dat in tiles.items(): + group.create_dataset(tile_id, data=dat) + f.attrs["n_tiles"] = len(tiles) + + +class TestCampaignReader(unittest.TestCase): + """Campaign galaxy catalogue reader.""" + + def setUp(self): + self._tmp = tempfile.TemporaryDirectory() + self._dir = Path(self._tmp.name) + self._tiles = { + "000.000": make_galaxy_data(3), + "001.000": make_galaxy_data(2, offset=100), + } + + def tearDown(self): + self._tmp.cleanup() + + def _expected(self): + return np.concatenate(list(self._tiles.values())) + + def test_legacy_layout(self): + path = self._dir / "final_cat_CAMPAIGN.hdf5" + write_campaign(path, "legacy", self._tiles) + + dat = catalog.read_campaign_catalogue(path, verbose=False) + + self.assertEqual(len(dat), 5) + npt.assert_array_equal(np.sort(dat["RA"]), np.sort(self._expected()["RA"])) + + def test_flat_layout(self): + path = self._dir / "final_cat_CAMPAIGN.hdf5" + write_campaign(path, "flat", self._tiles) + + dat = catalog.read_campaign_catalogue(path, verbose=False) + + self.assertEqual(len(dat), 5) + npt.assert_array_equal(np.sort(dat["RA"]), np.sort(self._expected()["RA"])) + + def test_layouts_agree(self): + legacy = self._dir / "legacy.hdf5" + flat = self._dir / "flat.hdf5" + write_campaign(legacy, "legacy", self._tiles) + write_campaign(flat, "flat", self._tiles) + + npt.assert_array_equal( + catalog.read_campaign_catalogue(legacy, verbose=False), + catalog.read_campaign_catalogue(flat, verbose=False), + ) + + def test_param_list_restriction(self): + path = self._dir / "final_cat_CAMPAIGN.hdf5" + write_campaign(path, "flat", self._tiles) + + dat = catalog.read_campaign_catalogue( + path, param_list=["RA", "Dec"], verbose=False + ) + + self.assertEqual(tuple(dat.dtype.names), ("RA", "Dec")) + + def test_missing_column_raises(self): + path = self._dir / "final_cat_CAMPAIGN.hdf5" + write_campaign(path, "flat", self._tiles) + + with self.assertRaises(KeyError) as ctx: + catalog.read_campaign_catalogue( + path, param_list=["RA", "NOT_A_COLUMN"], verbose=False + ) + self.assertIn("NOT_A_COLUMN", str(ctx.exception)) + + def test_ambiguous_layout_raises(self): + path = self._dir / "ambiguous.hdf5" + with h5py.File(path, "w") as f: + f.create_group("tiles") + f.create_group("other") + + with self.assertRaises(ValueError): + catalog.read_campaign_catalogue(path, verbose=False) + + +class TestStarCatalogueReader(unittest.TestCase): + """Campaign star catalogue reader.""" + + def setUp(self): + self._tmp = tempfile.TemporaryDirectory() + self._dir = Path(self._tmp.name) + + dtype = np.dtype([(name, "f8") for name in catalog.STAR_CAT_COLUMNS]) + self._exposures = { + "2110000p": np.zeros(4, dtype=dtype), + "2110001p": np.ones(6, dtype=dtype), + } + + self._path = self._dir / "full_starcat_CAMPAIGN.hdf5" + with h5py.File(self._path, "w") as f: + group = f.create_group("exposures") + for exp, dat in self._exposures.items(): + group.create_dataset(exp, data=dat) + + def tearDown(self): + self._tmp.cleanup() + + def test_read_hdf5(self): + dat = catalog.read_star_catalogue(self._path, verbose=False) + + self.assertEqual(len(dat), 10) + self.assertEqual(tuple(dat.dtype.names), catalog.STAR_CAT_COLUMNS) + npt.assert_array_equal(dat["MAG"][:4], np.zeros(4)) + npt.assert_array_equal(dat["MAG"][4:], np.ones(6)) + + def test_read_fits(self): + from astropy.io import fits + + fits_path = self._dir / "full_starcat-0000000.fits" + dat = np.concatenate(list(self._exposures.values())) + fits.BinTableHDU(data=dat).writeto(fits_path) + + out = catalog.read_star_catalogue(str(fits_path), verbose=False) + + self.assertEqual(len(out), 10) + npt.assert_array_equal(np.asarray(out["MAG"]), dat["MAG"]) + + +class TestMaskCut(unittest.TestCase): + """Mask-column galaxy selection cut.""" + + def setUp(self): + self._dat = make_galaxy_data(5) + + def test_default_columns(self): + self._dat["MASK_n4"][0] = True + self._dat["MASK_n1024"][3] = True + + npt.assert_array_equal( + galaxy.mask_cut(self._dat), + np.array([False, True, True, False, True]), + ) + + def test_explicit_column_list(self): + self._dat["MASK_n4"][0] = True + self._dat["MASK_n8"][1] = True + + npt.assert_array_equal( + galaxy.mask_cut(self._dat, ["MASK_n8"]), + np.array([True, False, True, True, True]), + ) + + def test_empty_column_list_keeps_everything(self): + self._dat["MASK_n4"][:] = True + + npt.assert_array_equal( + galaxy.mask_cut(self._dat, []), np.ones(len(self._dat), dtype=bool) + ) + + def test_missing_column_raises(self): + with self.assertRaises(KeyError) as ctx: + galaxy.mask_cut(self._dat, ["MASK_n16"]) + self.assertIn("MASK_n16", str(ctx.exception)) + + def test_v1_catalogue_raises(self): + dat = np.zeros(3, dtype=[("IMAFLAGS_ISO", "i2")]) + + with self.assertRaises(KeyError): + galaxy.mask_cut(dat) + + +class TestCampaignMerge(unittest.TestCase): + """Merge of several campaign catalogues into a joint catalogue.""" + + def setUp(self): + self._tmp = tempfile.TemporaryDirectory() + self._dir = Path(self._tmp.name) + + self._paths = [] + for name, layout, n_obj in (("W3", "legacy", 3), ("SGC", "flat", 2)): + path = self._dir / f"final_cat_{name}.hdf5" + write_campaign(path, layout, {"000.000": make_galaxy_data(n_obj)}) + self._paths.append(str(path)) + + self._obj = JointCat() + self._obj._params["verbose"] = False + + def tearDown(self): + self._tmp.cleanup() + + def test_campaign_name(self): + self.assertEqual( + JointCat.campaign_name("/some/dir/final_cat_W3.hdf5"), "W3" + ) + + def test_merge(self): + dat = self._obj.merge_catalogues(self._paths) + + self.assertEqual(len(dat), 5) + self.assertIn("campaign", dat.dtype.names) + npt.assert_array_equal( + dat["campaign"], np.array([b"W3"] * 3 + [b"SGC"] * 2) + ) + npt.assert_array_equal(dat["RA"][:3], np.arange(3)) + + def test_merge_incompatible_columns_raises(self): + other = self._dir / "final_cat_X.hdf5" + dat = np.zeros(2, dtype=[("RA", "f8")]) + write_campaign(other, "flat", {"000.000": dat}) + + with self.assertRaises(ValueError): + self._obj.merge_catalogues(self._paths + [str(other)]) + + def test_no_input_raises(self): + with self.assertRaises(ValueError): + self._obj.get_input_paths() + + +if __name__ == "__main__": + unittest.main() From 7ec470b772df1c6c62eafda7fd57972c5ae82fc3 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 20:41:51 -0400 Subject: [PATCH 05/17] catalog: safe, low-memory campaign reading and merging - concatenate_datasets preallocates the output and fills it column by column, so peak memory is the packed output plus one tile instead of the full-width uncut catalogue; the verbose estimate uses the real itemsize. - validate the requested columns against every tile dataset, not only the first, and name the offending dataset in the error. - read_campaign_catalogue checks the root n_tiles attribute and refuses a truncated file; campaign_shape reports row count and dtype from metadata. - JointCat.merge_catalogues preallocates the merged array from that first pass (no more accumulate-then-concatenate, which doubled peak memory), promotes each column's dtype across all campaigns so a wider string or integer column in a later file is no longer silently truncated, and refuses multi-dimensional columns explicitly. - reduce_mem reduces int32 to int16 (int8 wrapped N_EPOCH/CCD_NB values) and every assignment is range-checked, raising instead of wrapping. - check_matching drops the identity index over d1 that the retired footprint prefilter left behind, and extract_info applies mask_cut to the matched subset instead of the whole catalogue. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- scripts/calibration/extract_info.py | 4 +- src/sp_validation/catalog.py | 125 ++++++++++++++++++++++---- src/sp_validation/catalog_builders.py | 123 +++++++++++++++++++------ 3 files changed, 202 insertions(+), 50 deletions(-) diff --git a/scripts/calibration/extract_info.py b/scripts/calibration/extract_info.py index b65c1975..ed4fb145 100644 --- a/scripts/calibration/extract_info.py +++ b/scripts/calibration/extract_info.py @@ -142,7 +142,7 @@ # #### Match to all objects if star_cat_path: - ind_star, mask_area_tiles, n_star_tot = spv_cat.check_matching( + ind_star, n_star_tot = spv_cat.check_matching( d_star, dd, ["RA", "DEC"], @@ -160,7 +160,7 @@ m_star = ( (dd["FLAGS"][ind_star] == 0) - & galaxy.mask_cut(dd, mask_columns)[ind_star] + & galaxy.mask_cut(dd[ind_star], mask_columns) & (dd["NGMIX_MCAL_FLAGS"][ind_star] == 0) & (dd["NGMIX_G1_PSF_ORIG_NOSHEAR"][ind_star] != -10) ) diff --git a/src/sp_validation/catalog.py b/src/sp_validation/catalog.py index 52e850e4..b353f8ef 100644 --- a/src/sp_validation/catalog.py +++ b/src/sp_validation/catalog.py @@ -176,22 +176,20 @@ def check_matching( ------- ind : array of int index list of d2 of objects that were matched to d1 - mask_area_tiles : array of int - index list of tiles in footprint + n_tot : int + number of objects in d1 """ - mask_area_tiles = np.arange(len(d1)) - # Match stars from exposure (PSF) catalogue to total catalogue ind = match_stars2( d2[keys_2[0]], d2[keys_2[1]], - d1[keys_1[0]][mask_area_tiles], - d1[keys_1[1]][mask_area_tiles], + d1[keys_1[0]], + d1[keys_1[1]], thresh=thresh, ) - n_tot = len(d1[keys_1[0]][mask_area_tiles]) + n_tot = len(d1[keys_1[0]]) msg = ( "Number of matched stars from exposures to total catalogue = " + f"{len(ind)}/{n_tot} = {len(ind) / n_tot:.1%}" @@ -207,7 +205,7 @@ def check_matching( ) io.print_stats(msg, stats_file, verbose=verbose) - return ind, mask_area_tiles, n_tot + return ind, n_tot def check_invalid(dd, key, val, stats_file, name=None, verbose=False): @@ -816,12 +814,13 @@ def find_dataset_group(hdf5_file): node = node[keys[0]] -def _check_columns(dtype, param_list, file_path): +def _check_columns(dtype, param_list, file_path, dataset_key=None): """Raise a clear error if requested columns are absent from the data.""" missing = [col for col in param_list if col not in (dtype.names or ())] if missing: + where = f" (dataset {dataset_key!r})" if dataset_key is not None else "" raise KeyError( - f"Column(s) {missing} not found in catalogue {file_path}." + f"Column(s) {missing} not found in catalogue {file_path}{where}." + f" Available columns: {sorted(dtype.names or ())}" ) @@ -832,6 +831,10 @@ def concatenate_datasets(group, param_list=None, file_path="", verbose=True): Concatenate every dataset of an HDF5 group into one structured array, optionally restricted to a list of columns. + The output array is preallocated and filled slice by slice, so peak + memory is the output catalogue plus one tile, not the full-width + uncut catalogue. + Parameters ---------- group : h5py.Group @@ -849,27 +852,110 @@ def concatenate_datasets(group, param_list=None, file_path="", verbose=True): concatenated structured array """ - keys = list(group) + keys = sorted(group) + if not keys: + raise ValueError(f"No datasets found in catalogue {file_path}") + + # Validate every dataset up front: a column may be missing from any tile, + # not only the first one. + for key in keys: + if param_list is not None: + _check_columns(group[key].dtype, param_list, file_path, key) + + dtype_in = group[keys[0]].dtype if param_list is not None: - _check_columns(group[keys[0]].dtype, param_list, file_path) + dtype_out = np.dtype([(name, dtype_in[name]) for name in param_list]) + else: + dtype_out = dtype_in n_rows = sum(group[key].shape[0] for key in keys) if verbose: - n_cols = len(param_list) if param_list is not None else len(group[keys[0]].dtype) print( f"Reading {len(keys)} datasets," - + f" estimating {n_cols * n_rows * 8 / 1024**3:.1f}" - + f" Gb memory for the ({n_cols} x {n_rows}) data array ..." + + f" estimating {dtype_out.itemsize * n_rows / 1024**3:.1f}" + + f" Gb memory for the ({len(dtype_out.names)} x {n_rows}) data array ..." ) - data_list = [] + data_out = np.empty(n_rows, dtype=dtype_out) + start = 0 for key in tqdm.tqdm(keys, disable=not verbose): data = group[key][()] + end = start + len(data) + for name in dtype_out.names: + data_out[name][start:end] = data[name] + start = end + del data + + return data_out + + +def check_n_tiles(hdf5_file, group, file_path): + """Check Number of Tiles. + + Compare the number of datasets found against the ``n_tiles`` root + attribute the ShapePipe v2 product carries, to catch a catalogue that + was truncated by an interrupted merge job or file transfer. + + Parameters + ---------- + hdf5_file : h5py.File + open input file + group : h5py.Group + group holding the per-tile datasets + file_path : str + input file path, for the error message + + Raises + ------ + ValueError + if the number of datasets differs from the ``n_tiles`` attribute + + """ + n_tiles = hdf5_file.attrs.get("n_tiles") + if n_tiles is None: + return + n_found = len(group) + if int(n_tiles) != n_found: + raise ValueError( + f"Catalogue {file_path} declares n_tiles = {int(n_tiles)} but holds" + + f" {n_found} tile dataset(s); the file is incomplete." + ) + + +def campaign_shape(file_path, param_list=None): + """Campaign Shape. + + Return the number of rows and the structured dtype of a campaign + catalogue without reading its data, so that a merged output array can + be preallocated. + + Parameters + ---------- + file_path : str + input file path + param_list : list of str, optional + columns to keep; default is ``None`` (keep all) + + Returns + ------- + tuple + number of rows (int) and dtype (numpy.dtype) + + """ + with h5py.File(file_path, "r") as hdf5_file: + group = find_dataset_group(hdf5_file) + check_n_tiles(hdf5_file, group, file_path) + keys = sorted(group) + if not keys: + raise ValueError(f"No datasets found in catalogue {file_path}") + n_rows = sum(group[key].shape[0] for key in keys) + dtype_in = group[keys[0]].dtype if param_list is not None: - data = data[param_list] - data_list.append(data) + for key in keys: + _check_columns(group[key].dtype, param_list, file_path, key) + dtype_in = np.dtype([(name, dtype_in[name]) for name in param_list]) - return np.concatenate(data_list, axis=0) + return n_rows, dtype_in def read_campaign_catalogue( @@ -905,6 +991,7 @@ def read_campaign_catalogue( with h5py.File(file_path, "r") as hdf5_file: group = find_dataset_group(hdf5_file) + check_n_tiles(hdf5_file, group, file_path) return concatenate_datasets( group, param_list=param_list, file_path=file_path, verbose=verbose ) diff --git a/src/sp_validation/catalog_builders.py b/src/sp_validation/catalog_builders.py index 00b00694..899d4db5 100644 --- a/src/sp_validation/catalog_builders.py +++ b/src/sp_validation/catalog_builders.py @@ -240,6 +240,37 @@ def close_hd5(self): self._hd5file.close() +def _promote(dtype_a, dtype_b): + """Return a dtype that holds both input dtypes without truncation.""" + if dtype_a == dtype_b: + return dtype_a + if dtype_a.kind in "SU" and dtype_b.kind in "SU": + kind = "U" if "U" in (dtype_a.kind, dtype_b.kind) else "S" + size_a = dtype_a.itemsize // (4 if dtype_a.kind == "U" else 1) + size_b = dtype_b.itemsize // (4 if dtype_b.kind == "U" else 1) + return np.dtype(f"{kind}{max(size_a, size_b)}") + + return np.promote_types(dtype_a, dtype_b) + + +def _checked_assign(target, start, end, values, name): + """Assign `values` into `target[start:end]`, refusing a lossy cast.""" + target[start:end] = values + written = target[start:end] + if written.dtype == values.dtype: + return + if written.dtype.kind in "fc": + bad = np.isfinite(values) & ~np.isfinite(written) + else: + bad = written != values + if np.any(bad): + raise ValueError( + f"Column {name!r} cannot be stored as {written.dtype}:" + + f" {int(np.sum(bad))} value(s) overflow or are truncated." + + " Disable reduce_mem or widen the output dtype." + ) + + class JointCat(BaseCat): """Joint Cat. @@ -391,22 +422,26 @@ def dtype_out(self, name, dtype_in): return dtype_in if name not in cols_keep_dtype: if dtype_in.kind == "f" and dtype_in.itemsize == 8: - return np.float32 + return np.dtype(np.float32) if dtype_in.kind == "i" and dtype_in.itemsize == 4: - return np.int8 + # int32 -> int16, not int8: int8 cannot hold e.g. N_EPOCH or + # CCD_NB values and wrapped them silently. Values that do not + # fit are caught at assignment time by ``_checked_assign``. + return np.dtype(np.int16) return dtype_in - def output_dtype(self, dtype_in, n_char_campaign): + def output_dtype(self, dtypes_in, n_char_campaign): """Output Dtype. Return the merged-catalogue dtype: the input columns (possibly - reduced in precision) plus a ``campaign`` column. + reduced in precision, and promoted to a common type across all input + campaigns) plus a ``campaign`` column. Parameters ---------- - dtype_in : numpy.dtype - structured dtype of an input campaign catalogue + dtypes_in : numpy.dtype or list of numpy.dtype + structured dtype(s) of the input campaign catalogues n_char_campaign : int width of the campaign name column @@ -415,10 +450,40 @@ def output_dtype(self, dtype_in, n_char_campaign): numpy.dtype output structured dtype + Raises + ------ + ValueError + if the inputs have different column sets, or a column is + multi-dimensional (campaign catalogues are scalar-column only) + """ - fields = [ - (name, self.dtype_out(name, dtype_in[name])) for name in dtype_in.names - ] + if isinstance(dtypes_in, np.dtype): + dtypes_in = [dtypes_in] + + names = dtypes_in[0].names + for dtype_in in dtypes_in[1:]: + if set(dtype_in.names) != set(names): + raise ValueError( + "Campaign catalogues have incompatible column sets:" + + f" {sorted(names)} vs {sorted(dtype_in.names)}" + ) + + fields = [] + for name in names: + subs = [dtype_in[name] for dtype_in in dtypes_in] + for sub in subs: + if sub.subdtype is not None: + raise ValueError( + f"Column {name!r} is multi-dimensional (shape" + + f" {sub.subdtype[1]}); campaign catalogues are" + + " expected to hold scalar columns only." + ) + # Promote across campaigns so a wider string or integer column in + # a later file is not silently truncated or overflowed. + promoted = subs[0] + for sub in subs[1:]: + promoted = _promote(promoted, sub) + fields.append((name, self.dtype_out(name, promoted))) fields.append(("campaign", np.dtype(f"S{n_char_campaign}"))) return np.dtype(fields) @@ -494,8 +559,18 @@ def merge_catalogues(self, input_paths): campaigns = [self.campaign_name(path) for path in input_paths] n_char_campaign = max(len(name) for name in campaigns) - data_list = [] - dtype_out = None + # First pass over file metadata only (row counts and dtypes), so the + # merged array is allocated once and filled in place, instead of + # concatenating per-campaign copies (peak memory 2x the output). + shapes = [ + sp_cat.campaign_shape(path, param_list=param_list) + for path in input_paths + ] + n_total = sum(n_rows for n_rows, _ in shapes) + dtype_out = self.output_dtype([dtype for _, dtype in shapes], n_char_campaign) + + dat_all = np.empty(n_total, dtype=dtype_out) + start = 0 for input_path, campaign in zip(input_paths, campaigns): dat = sp_cat.read_campaign_catalogue( input_path, @@ -503,28 +578,18 @@ def merge_catalogues(self, input_paths): verbose=self._params["verbose"], ) - if dtype_out is None: - dtype_out = self.output_dtype(dat.dtype, n_char_campaign) - elif set(dat.dtype.names) != set(dtype_out.names) - {"campaign"}: - raise ValueError( - f"Campaign catalogue {input_path} has columns" - + f" {sorted(dat.dtype.names)}, incompatible with" - + f" {sorted(set(dtype_out.names) - {'campaign'})}" - ) - - dat_out = np.empty(len(dat), dtype=dtype_out) + end = start + len(dat) for name in dat.dtype.names: - dat_out[name] = dat[name] - dat_out["campaign"] = campaign.encode() - data_list.append(dat_out) + _checked_assign(dat_all[name], start, end, dat[name], name) + dat_all["campaign"][start:end] = campaign.encode() + start = end if self._params["verbose"]: print( f"{campaign}: added {len(dat)}" + f" (~{format.millify(len(dat))}) objects." ) - - dat_all = np.concatenate(data_list, axis=0) + del dat if self._params["verbose"]: print( @@ -950,8 +1015,8 @@ def write_hdf5_header(self, hd5file): ---------- hd5file : h5py.File input HDF5 file - patches : list, optional - input patches, list of str, default is ``None`` + campaigns : list, optional + input campaign names, list of str, default is ``None`` """ super().write_hdf5_header(hd5file) @@ -1027,7 +1092,7 @@ def read_cat(self, load_into_memory=False): verbose = self._params["verbose"] # Image-simulation path: a single per-run comprehensive catalogue in - # FITS, not the joined multi-patch HDF5 the data path builds. Read the + # FITS, not the joined multi-campaign HDF5 the data path builds. Read the # FITS table directly into memory; there is no separate data_ext group. extension = os.path.splitext(fpath)[1] if extension == ".fits": From bf2ea8614a76bf4758a1eb17c682d89091d17dfc Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 20:42:07 -0400 Subject: [PATCH 06/17] config: v2 mask configuration with the MASK_n* columns Add config/calibration/mask_v2.0.yaml: the v1.X.11 cut set with IMAFLAGS_ISO and the v1 post-processing masks replaced by the boolean MASK_n columns (True = masked, hence kind: equal, value: False), the coverage bits and MASK_n2048 listed but commented out. The v1.X configs are left untouched: each describes a legacy catalogue that really has IMAFLAGS_ISO, and rewriting them would break reproducing published versions. plots.sky_plots no longer hardcodes the v1 label set; it combines whichever of IMAFLAGS_ISO / MASK_n* / npoint3 / 1024_Maximask the config declared. The comprehensive-to-minimal demo asks for the v2 mask labels. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- config/calibration/mask_v2.0.yaml | 122 ++++++++++++++++++ .../demo_comprehensive_to_minimal_cat.py | 9 +- src/sp_validation/plots.py | 20 ++- 3 files changed, 145 insertions(+), 6 deletions(-) create mode 100644 config/calibration/mask_v2.0.yaml diff --git a/config/calibration/mask_v2.0.yaml b/config/calibration/mask_v2.0.yaml new file mode 100644 index 00000000..db824f59 --- /dev/null +++ b/config/calibration/mask_v2.0.yaml @@ -0,0 +1,122 @@ +# Config file for masking and calibration, ShapePipe v2 catalogues. +# +# ShapePipe v2 replaces the single IMAFLAGS_ISO column with per-reason boolean +# mask columns MASK_n, where True means MASKED. They are therefore cut +# with `kind: equal, value: False` (keep the un-masked objects): +# +# MASK_n1, MASK_n2 bright/faint star halos +# MASK_n4 stars +# MASK_n8 manual galaxy mask +# MASK_n16 .. MASK_n256 per-band coverage +# MASK_n1024 MaxiMask +# MASK_n2048 no Pan-STARRS z2 coverage +# +# The default selection (see sp_validation.galaxy.DEFAULT_MASK_COLUMNS) is +# n4 + n1 + n2 + n8 + n1024; the coverage bits and n2048 are listed here but +# commented out, since ORing all of them masks essentially everything. + +# General parameters (can also given on command line) +params: + input_path: unions_shapepipe_comprehensive_2025_v2.0.hdf5 + cmatrices: False + sky_regions: False + verbose: True + +# Masks +## Using columns in 'dat' group (ShapePipe flags) +dat: + # SExtractor flags + - col_name: FLAGS + label: SE FLAGS + kind: smaller_equal + value: 3 + + # Duplicate objects + - col_name: overlap + label: tile overlap + kind: equal + value: True + + # ShapePipe masks (boolean, True = masked) + - col_name: MASK_n4 + label: "stars" + kind: equal + value: False + - col_name: MASK_n1 + label: "faint star halos" + kind: equal + value: False + - col_name: MASK_n2 + label: "bright star halos" + kind: equal + value: False + - col_name: MASK_n8 + label: "manual mask" + kind: equal + value: False + - col_name: MASK_n1024 + label: "maximask" + kind: equal + value: False + + # Per-band coverage and Pan-STARRS z2; enable as required + # - col_name: MASK_n16 + # label: "coverage (n16)" + # kind: equal + # value: False + # - col_name: MASK_n2048 + # label: "no PS-z2" + # kind: equal + # value: False + + # Number of epochs + - col_name: N_EPOCH + label: r"$n_{\rm epoch}$" + kind: greater_equal + value: 2 + + # Magnitude range + - col_name: mag + label: mag range + kind: range + value: [15, 30] + + # ngmix flags + - col_name: NGMIX_MCAL_TYPES_FAIL + label: "ngmix moments failure" + kind: equal + value: 0 + + # invalid PSF ellipticities + - col_name: NGMIX_G1_PSF_ORIG_NOSHEAR + label: "bad PSF ellipticity comp 1" + kind: not_equal + value: -10 + - col_name: NGMIX_G2_PSF_ORIG_NOSHEAR + label: "bad PSF ellipticity comp 2" + kind: not_equal + value: -10 + +## Using columns in 'dat_ext' group (post-processing flags) +## ShapePipe v2 carries the imaging masks in the 'dat' group above, so this +## group is empty unless external masks are added. +dat_ext: [] + +# Metacal parameters +metacal: + # Ellipticity dispersion + sigma_eps_prior: 0.34 + + # Signal-to-noise range + gal_snr_min: 10 + gal_snr_max: 500 + + # Relative-size (hlr / hlr_psf) range + gal_rel_size_min: 0.707 + gal_rel_size_max: 3 + + # Correct relative size for ellipticity? + gal_size_corr_ell: False + + # Weight for global response matrix, None for unweighted mean + global_R_weight: w diff --git a/scripts/examples/demo_comprehensive_to_minimal_cat.py b/scripts/examples/demo_comprehensive_to_minimal_cat.py index 11e5ed5a..30b6b1df 100644 --- a/scripts/examples/demo_comprehensive_to_minimal_cat.py +++ b/scripts/examples/demo_comprehensive_to_minimal_cat.py @@ -50,13 +50,18 @@ # + # List of masks to apply +# (labels as declared in config/calibration/mask_v2.0.yaml; ShapePipe v2 +# replaces IMAFLAGS_ISO and the v1 post-processing masks by MASK_n) masks_to_apply = [ "overlap", - "IMAFLAGS_ISO", + "MASK_n4", + "MASK_n1", + "MASK_n2", + "MASK_n8", + "MASK_n1024", "NGMIX_MCAL_TYPES_FAIL", "NGMIX_G1_PSF_ORIG_NOSHEAR", "NGMIX_G2_PSF_ORIG_NOSHEAR", - "8_Manual", ] # List of masks not to apply and not to copy to minimal catalogue diff --git a/src/sp_validation/plots.py b/src/sp_validation/plots.py index b21599c5..60b5b2b4 100644 --- a/src/sp_validation/plots.py +++ b/src/sp_validation/plots.py @@ -20,6 +20,7 @@ # Imported for its import-time side effect: sets matplotlib rcParams (plot style). import sp_validation.plot_style # noqa: F401 +from sp_validation.galaxy import MASK_COLUMNS from sp_validation.masks import Mask @@ -477,8 +478,13 @@ def sky_plots(dat, masks, labels, zoom_ra, zoom_dec): # No mask plot_area_mask(ra, dec, zoom) - # SExtractor and SP flags - m_flags = masks[labels["FLAGS"]]._mask & masks[labels["IMAFLAGS_ISO"]]._mask + # SExtractor and SP flags. ShapePipe v2 splits the single IMAFLAGS_ISO + # mask into per-reason MASK_n columns, so combine whichever of them + # the mask config declared. + m_flags = masks[labels["FLAGS"]]._mask + for col in ("IMAFLAGS_ISO",) + tuple(MASK_COLUMNS): + if col in labels: + m_flags = m_flags & masks[labels[col]]._mask plot_area_mask(ra, dec, zoom, mask=m_flags) # Overlap regions @@ -486,11 +492,17 @@ def sky_plots(dat, masks, labels, zoom_ra, zoom_dec): plot_area_mask(ra, dec, zoom, mask=m_over) # Coverage mask - m_point = masks[labels["npoint3"]]._mask & m_over + # Rough pointing coverage; v2 encodes coverage in the MASK_n* columns + m_point = m_over + if "npoint3" in labels: + m_point = masks[labels["npoint3"]]._mask & m_over plot_area_mask(ra, dec, zoom, mask=m_point) # Maximask - m_maxi = masks[labels["1024_Maximask"]]._mask & m_point + m_maxi = m_point + for maxi_key in ("1024_Maximask", "MASK_n1024"): + if maxi_key in labels: + m_maxi = masks[labels[maxi_key]]._mask & m_point plot_area_mask(ra, dec, zoom, mask=m_maxi) # Combined mask over all supplied masks (was passed in by the caller before From ce2dfd511d313065fca256d8b8c71afcb883d973 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 20:42:08 -0400 Subject: [PATCH 07/17] scripts: finish the campaign rename in configs and entry points - params.py interpolated the removed 'name' into two paths (NameError on import); use 'campaign', give it a real default, and point star_cat_path at full_starcat_.hdf5 (hdu_star_cat now documented as legacy FITS only). params_im_sim.py renames 'name' to 'campaign' likewise. - combine_results.get_area matches both the campaign and the legacy patch wording, and raises on a missing file or unmatched pattern instead of returning None / a 1 deg^2 placeholder that silently rescales densities. - merge_psf_cat writes the campaign *name* as a string column (FITS 'A'), matching JointCat; the old 1-based ordinal depended on -p argument order. - compute_m_bias_image_sims counts tiles via the n_tiles attribute or find_dataset_group, not a hardcoded legacy group inside a bare except. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- scripts/calibration/params.py | 14 ++++++----- scripts/combine_results.py | 28 ++++++++++++++-------- scripts/compute_m_bias_image_sims.py | 36 +++++++++++++++------------- scripts/merge_psf_cat.py | 9 +++++-- workflow/image_sims/params_im_sim.py | 10 ++++---- 5 files changed, 57 insertions(+), 40 deletions(-) diff --git a/scripts/calibration/params.py b/scripts/calibration/params.py index e4866942..03504946 100644 --- a/scripts/calibration/params.py +++ b/scripts/calibration/params.py @@ -25,8 +25,9 @@ # Survey parameters -## Campaign name (the tile list processed by ShapePipe). Put None if n/a -campaign = None +## Campaign name (the tile list processed by ShapePipe); required, it names +## the ShapePipe v2 products (final_cat_.hdf5, ...) +campaign = "W3" ## Area of a tile in deg^2 area_tile = 0.25 @@ -43,19 +44,20 @@ data_dir = "." ### Tile IDs -path_tile_ID = f"{data_dir}/tiles_{name}.txt" +path_tile_ID = f"{data_dir}/tiles_{campaign}.txt" ### Weak-lensing galaxy catalog name -galaxy_cat_path = f"{data_dir}/final_cat_{name}.hdf5" +galaxy_cat_path = f"{data_dir}/final_cat_{campaign}.hdf5" print(f"Galaxy catalogue = {galaxy_cat_path}") ## Parameter list; optional, set to `None` if not required param_list_path = f"{data_dir}/cfis/final_cat.param" ### Star and PSF catalog name; optional, set to `None` if not required -star_cat_path = f"{data_dir}/full_starcat-0000000.fits" +star_cat_path = f"{data_dir}/full_starcat_{campaign}.hdf5" -# HDU number of star and PSF catalogue +# HDU number of star and PSF catalogue; only used for the legacy FITS star +# catalogue (a path ending in .fits), ignored for the v2 HDF5 product hdu_star_cat = 1 ### External mask; optional, set to `None` if not required diff --git a/scripts/combine_results.py b/scripts/combine_results.py index 50cc91c5..1b3b4f85 100755 --- a/scripts/combine_results.py +++ b/scripts/combine_results.py @@ -217,18 +217,26 @@ def print_all( def get_area(fname): + """Return the unmasked area in deg^2 read from an area.txt file. - if os.path.exists(fname): - with open(fname) as f: - lines = f.readlines() - for line in lines: - m = re.search("nmasked campaign area without overlap = (.*) deg", line) - if m: - return float(m[1]) + Accepts both the v2 wording ("campaign") and the legacy one ("patch"), + so results computed before the campaign rename can still be combined. + Raises rather than returning a placeholder: a wrong area silently + rescales every density. + """ + if not os.path.exists(fname): + raise FileNotFoundError(f"No file {fname} found to obtain area") + + with open(fname) as f: + lines = f.readlines() + for line in lines: + m = re.search( + r"nmasked (?:campaign|patch) area without overlap = (.*) deg", line + ) + if m: + return float(m[1]) - else: - print(f"Warning: No file {fname} found to obtain area") - return 1 + raise ValueError(f"No unmasked area found in file {fname}") def get_values(results, stats_files, shape, use_keys, area_deg2=-1): diff --git a/scripts/compute_m_bias_image_sims.py b/scripts/compute_m_bias_image_sims.py index 04ca49b7..eb0a76a8 100644 --- a/scripts/compute_m_bias_image_sims.py +++ b/scripts/compute_m_bias_image_sims.py @@ -57,23 +57,25 @@ def parse_args(): def get_n_tiles(grids_dir, num): - """Detect number of tiles from final_cat HDF5 files.""" - try: - import h5py - - # Count tiles in first sim's final_cat - for sim in ["1z2z_grid", "1m2z_grid", "1p2z_grid", "1z2m_grid", "1z2p_grid"]: - sim_name = f"{sim}_{num}" - final_cat = os.path.join(grids_dir, sim_name, f"final_cat_{sim_name}.hdf5") - if os.path.isfile(final_cat): - with h5py.File(final_cat, "r") as hf: - if "patches" in hf: - n_tiles = sum( - 1 for patch in hf["patches"] for _ in hf[f"patches/{patch}"] - ) - return n_tiles - except Exception: - pass + """Detect number of tiles from final_cat HDF5 files. + + Layout-agnostic: uses ``sp_validation.catalog.find_dataset_group``, so it + works on both the legacy nested and the flat per-tile HDF5 layouts. + """ + import h5py + + from sp_validation.catalog import find_dataset_group + + for sim in ["1z2z_grid", "1m2z_grid", "1p2z_grid", "1z2m_grid", "1z2p_grid"]: + sim_name = f"{sim}_{num}" + final_cat = os.path.join(grids_dir, sim_name, f"final_cat_{sim_name}.hdf5") + if os.path.isfile(final_cat): + with h5py.File(final_cat, "r") as hf: + n_tiles = hf.attrs.get("n_tiles") + if n_tiles is not None: + return int(n_tiles) + return len(find_dataset_group(hf)) + return None diff --git a/scripts/merge_psf_cat.py b/scripts/merge_psf_cat.py index adda9d22..8ca34cf4 100644 --- a/scripts/merge_psf_cat.py +++ b/scripts/merge_psf_cat.py @@ -143,7 +143,9 @@ def merge_catalogues(self, campaigns): for name in col_names: dat_all[name] = np.append(dat_all[name], dat[name]) - dat_all["campaign"] = np.append(dat_all["campaign"], [idx + 1] * len(dat)) + dat_all["campaign"] = np.append( + dat_all["campaign"], [campaign] * len(dat) + ) col_names = col_names + ("campaign",) @@ -152,7 +154,10 @@ def merge_catalogues(self, campaigns): if name != "campaign": my_format = "D" else: - my_format = "I" + # Store the campaign name, matching JointCat's string + # `campaign` column; an ordinal would depend on the order of + # the -p argument and could not be joined back. + my_format = f"A{max(len(c) for c in campaigns)}" column = fits.Column(name=name, array=dat_all[name], format=my_format) column_all.append(column) diff --git a/workflow/image_sims/params_im_sim.py b/workflow/image_sims/params_im_sim.py index 9c614e33..7ba7b7ba 100644 --- a/workflow/image_sims/params_im_sim.py +++ b/workflow/image_sims/params_im_sim.py @@ -27,11 +27,11 @@ # Survey parameters -## Field name -- derived from the run directory, which is named +## Campaign name -- derived from the run directory, which is named ## after the simulation (e.g. '1z2z_grid_1'), so one shared params file ## serves every sim. -name = os.path.basename(os.getcwd()) -print("Field name = {}".format(name)) +campaign = os.path.basename(os.getcwd()) +print("Campaign name = {}".format(campaign)) ## Area of a tile in deg^2 area_tile = 0.25 @@ -49,10 +49,10 @@ data_dir = "." ### Tile IDs -path_tile_ID = f"{data_dir}/tiles_{name}.txt" +path_tile_ID = f"{data_dir}/tiles_{campaign}.txt" ### Weak-lensing galaxy catalog name -galaxy_cat_path = f"{data_dir}/final_cat_{name}.hdf5" +galaxy_cat_path = f"{data_dir}/final_cat_{campaign}.hdf5" print(f"Galaxy catalogue = {galaxy_cat_path}") ## Parameter list; optional, set to `None` if not required From 2dbe5010e4562c3ebcf50ac02e9af9c6e61ce53e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 20:42:08 -0400 Subject: [PATCH 08/17] tests, docs: pin row order, cover the new failure modes Tests: assert the exact concatenation instead of sorted values, add a fixture whose keys are inserted out of order, and cover the truncated n_tiles file, a column missing from a later tile, dtype promotion across campaigns in both argument orders, reduce_mem overflow, and a multi-dimensional column. Docs: CLAUDE.md, scripts/calibration/README.md, homogenize_cat_extended.py and the catalog_builders docstrings now say campaign. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- CLAUDE.md | 4 +- scripts/calibration/README.md | 9 +- scripts/homogenize_cat_extended.py | 2 +- .../tests/test_campaign_readers.py | 115 +++++++++++++++++- 4 files changed, 119 insertions(+), 11 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index c35c20be..d1aa1c3d 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -65,10 +65,10 @@ Run via `./pipeline.sh` with flags: ### Configuration Main configuration in `scripts/calibration/params.py` with parameters: -- `name`: Field/patch identifier +- `campaign`: Campaign name (the ShapePipe tile list); names the input products - `data_dir`: Input data directory - `galaxy_cat_path`: Galaxy catalogue path (.fits/.hdf5) -- `star_cat_path`: Star catalogue path (.fits) +- `star_cat_path`: Star catalogue path (.hdf5, or legacy .fits) ### Key Dependencies - astropy, numpy, scipy for core calculations diff --git a/scripts/calibration/README.md b/scripts/calibration/README.md index 870c87ce..cd0a49e3 100644 --- a/scripts/calibration/README.md +++ b/scripts/calibration/README.md @@ -6,14 +6,13 @@ in order. See `docs/source/post_processing.md` for the full prose. | Step | Script | Does | |------|--------|------| -| 1 | `extract_info.py` | Extract metacal + diagnostic info per patch; create pre-calibration shear catalogues. Configured via `params.py`. | -| 2 | `create_joint_comprehensive_cat.py` | Merge the patch-wise comprehensive catalogues into one joint catalogue (front-end of `catalog_builders.JointCat`). | +| 1 | `extract_info.py` | Extract metacal + diagnostic info for one campaign; create pre-calibration shear catalogues. Configured via `params.py`. | +| 2 | `create_joint_comprehensive_cat.py` | Merge the per-campaign comprehensive catalogues into one joint catalogue (front-end of `catalog_builders.JointCat`). | | 3 | `demo_apply_hsp_masks.py` | Add the structural and coverage (HealSparse) masks. | | 4 | `calibrate_comprehensive_cat.py` | Galaxy selection + metacalibration. Uses the mask configs in `config/calibration/`. | `params.py` is the shared parameter template (paths, column names, survey constants) imported by `extract_info.py`; copy and edit it per run. -> **v2.0 note (Martin):** this clustering reflects the current (v1.4.x) reduction -> flow. v2.0 no longer has patches, so step 2 (and the per-patch structure of -> steps 1/3) will change. +> **v2.0 note:** ShapePipe v2 has no patches. Step 1 runs per campaign and +> step 2 merges a list of campaign catalogues (`final_cat_.hdf5`). diff --git a/scripts/homogenize_cat_extended.py b/scripts/homogenize_cat_extended.py index dd54922e..f197a941 100644 --- a/scripts/homogenize_cat_extended.py +++ b/scripts/homogenize_cat_extended.py @@ -3,7 +3,7 @@ Overwrite the extended catalog to replace the columns e1 and e2 from the extended catalog with the columns e1 and e2 from the non-extended one. Currently, the calibration of the columns e1 and e2 of the extended catalog -are calibrated per patch, while the non-extended catalog is calibrated on the whole footprint. +are calibrated per campaign, while the non-extended catalog is calibrated on the whole footprint. :Authors: Sacha Guerrini, Martin Kilbinger """ diff --git a/src/sp_validation/tests/test_campaign_readers.py b/src/sp_validation/tests/test_campaign_readers.py index ccba4fe5..7272b4a7 100644 --- a/src/sp_validation/tests/test_campaign_readers.py +++ b/src/sp_validation/tests/test_campaign_readers.py @@ -61,7 +61,8 @@ def tearDown(self): self._tmp.cleanup() def _expected(self): - return np.concatenate(list(self._tiles.values())) + """Expected concatenation: datasets in sorted-key order.""" + return np.concatenate([self._tiles[key] for key in sorted(self._tiles)]) def test_legacy_layout(self): path = self._dir / "final_cat_CAMPAIGN.hdf5" @@ -70,7 +71,7 @@ def test_legacy_layout(self): dat = catalog.read_campaign_catalogue(path, verbose=False) self.assertEqual(len(dat), 5) - npt.assert_array_equal(np.sort(dat["RA"]), np.sort(self._expected()["RA"])) + npt.assert_array_equal(dat["RA"], self._expected()["RA"]) def test_flat_layout(self): path = self._dir / "final_cat_CAMPAIGN.hdf5" @@ -79,7 +80,7 @@ def test_flat_layout(self): dat = catalog.read_campaign_catalogue(path, verbose=False) self.assertEqual(len(dat), 5) - npt.assert_array_equal(np.sort(dat["RA"]), np.sort(self._expected()["RA"])) + npt.assert_array_equal(dat["RA"], self._expected()["RA"]) def test_layouts_agree(self): legacy = self._dir / "legacy.hdf5" @@ -112,6 +113,54 @@ def test_missing_column_raises(self): ) self.assertIn("NOT_A_COLUMN", str(ctx.exception)) + + def test_row_order_is_tile_name_order(self): + """Keys inserted out of order still concatenate in name order.""" + path = self._dir / "unordered.hdf5" + tiles = { + "222.000": make_galaxy_data(2, offset=200), + "000.000": make_galaxy_data(2, offset=0), + "111.000": make_galaxy_data(2, offset=100), + } + with h5py.File(path, "w") as f: + group = f.create_group("tiles") + for tile_id, dat in tiles.items(): + group.create_dataset(tile_id, data=dat) + f.attrs["n_tiles"] = len(tiles) + + dat = catalog.read_campaign_catalogue(path, verbose=False) + + npt.assert_array_equal(dat["RA"], [0, 1, 100, 101, 200, 201]) + + def test_truncated_file_raises(self): + """n_tiles attribute larger than the number of datasets is fatal.""" + path = self._dir / "truncated.hdf5" + write_campaign(path, "flat", self._tiles) + with h5py.File(path, "a") as f: + f.attrs["n_tiles"] = 10 + + with self.assertRaises(ValueError) as ctx: + catalog.read_campaign_catalogue(path, verbose=False) + self.assertIn("incomplete", str(ctx.exception)) + + def test_column_missing_from_later_tile_raises_clear_error(self): + """A column absent from a non-first tile is named, with its dataset.""" + path = self._dir / "ragged.hdf5" + with h5py.File(path, "w") as f: + group = f.create_group("tiles") + group.create_dataset("000.000", data=make_galaxy_data(2)) + group.create_dataset( + "001.000", data=np.zeros(2, dtype=[("RA", "f8"), ("Dec", "f8")]) + ) + + with self.assertRaises(KeyError) as ctx: + catalog.read_campaign_catalogue( + path, param_list=["RA", "MAG_AUTO"], verbose=False + ) + message = str(ctx.exception) + self.assertIn("MAG_AUTO", message) + self.assertIn("001.000", message) + def test_ambiguous_layout_raises(self): path = self._dir / "ambiguous.hdf5" with h5py.File(path, "w") as f: @@ -242,6 +291,66 @@ def test_merge(self): ) npt.assert_array_equal(dat["RA"][:3], np.arange(3)) + + def test_merge_promotes_column_widths(self): + """A wider string/int column in a later file is not truncated.""" + dtype_narrow = np.dtype([("RA", "f8"), ("TILE_ID", "S7"), ("N", "i4")]) + dtype_wide = np.dtype([("RA", "f8"), ("TILE_ID", "S12"), ("N", "i8")]) + + narrow = np.zeros(1, dtype=dtype_narrow) + narrow["TILE_ID"] = b"123.456" + narrow["N"] = 7 + wide = np.zeros(1, dtype=dtype_wide) + wide["TILE_ID"] = b"999888.7776" + wide["N"] = 2**40 + + path_a = self._dir / "final_cat_AA.hdf5" + path_b = self._dir / "final_cat_BBBBBBBB.hdf5" + write_campaign(path_a, "flat", {"000.000": narrow}) + write_campaign(path_b, "flat", {"000.000": wide}) + + for paths in ([path_a, path_b], [path_b, path_a]): + dat = self._obj.merge_catalogues([str(path) for path in paths]) + by_campaign = { + name: row for name, row in zip(dat["campaign"], dat) + } + self.assertEqual(by_campaign[b"BBBBBBBB"]["TILE_ID"], b"999888.7776") + self.assertEqual(by_campaign[b"BBBBBBBB"]["N"], 2**40) + self.assertEqual(by_campaign[b"AA"]["TILE_ID"], b"123.456") + + def test_merge_reduce_mem_overflow_raises(self): + """reduce_mem never silently wraps out-of-range values.""" + dat = np.zeros(3, dtype=[("RA", "f8"), ("N_EPOCH", "i4")]) + dat["N_EPOCH"] = [3, 200, 300000] + path = self._dir / "final_cat_Z.hdf5" + write_campaign(path, "flat", {"000.000": dat}) + + self._obj._params["reduce_mem"] = True + with self.assertRaises(ValueError) as ctx: + self._obj.merge_catalogues([str(path)]) + self.assertIn("N_EPOCH", str(ctx.exception)) + + def test_merge_reduce_mem_in_range_ok(self): + dat = np.zeros(2, dtype=[("RA", "f8"), ("N_EPOCH", "i4")]) + dat["N_EPOCH"] = [3, 200] + path = self._dir / "final_cat_Y.hdf5" + write_campaign(path, "flat", {"000.000": dat}) + + self._obj._params["reduce_mem"] = True + out = self._obj.merge_catalogues([str(path)]) + + npt.assert_array_equal(out["N_EPOCH"], [3, 200]) + self.assertEqual(out.dtype["RA"], np.dtype("f8")) # RA keeps precision + + def test_merge_rejects_multidimensional_column(self): + dat = np.zeros(2, dtype=[("RA", "f8"), ("XY", "f8", (2,))]) + path = self._dir / "final_cat_M.hdf5" + write_campaign(path, "flat", {"000.000": dat}) + + with self.assertRaises(ValueError) as ctx: + self._obj.merge_catalogues([str(path)]) + self.assertIn("multi-dimensional", str(ctx.exception)) + def test_merge_incompatible_columns_raises(self): other = self._dir / "final_cat_X.hdf5" dat = np.zeros(2, dtype=[("RA", "f8")]) From d772b3fd1e07f43f994747eb6a95a06b7130e36b Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 21:31:49 -0400 Subject: [PATCH 09/17] catalog: promote dtypes across tiles, stream the campaign merge Three defects in the campaign reading path, found in review: - concatenate_datasets and campaign_shape both took the output dtype from the first dataset alone, so a campaign whose tiles differ (S7 next to S12 tile IDs, i2 next to i4, f4 next to f8 after a partial reprocessing) had the later tiles silently truncated, downcast or wrapped. np.concatenate, which this code replaced, promoted. Both now build the dtype with group_dtype(), promoting every column across every dataset; catalog_builders._promote becomes an alias of the shared catalog.promote_dtypes rather than a second copy of it. - merge_catalogues still held a whole campaign in memory next to the preallocated output, the very thing its comment claimed the rewrite avoided -- and with one campaign per merge in v2, that is the normal case, ~2x the merged catalogue at DR6 scale. It now fills the output tile by tile through the new catalog.iter_campaign_tiles(), so peak memory is the output plus a single tile. - write_hdf5_file wrote the merged array twice, create_dataset(data=...) followed by an immediate dset[:] = dat_all. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- src/sp_validation/catalog.py | 113 ++++++++++++++++-- src/sp_validation/catalog_builders.py | 41 +++---- .../tests/test_campaign_readers.py | 61 ++++++++++ 3 files changed, 183 insertions(+), 32 deletions(-) diff --git a/src/sp_validation/catalog.py b/src/sp_validation/catalog.py index b353f8ef..5f846ef4 100644 --- a/src/sp_validation/catalog.py +++ b/src/sp_validation/catalog.py @@ -814,6 +814,72 @@ def find_dataset_group(hdf5_file): node = node[keys[0]] +def promote_dtypes(dtype_a, dtype_b): + """Promote Dtypes. + + Return a scalar dtype that holds both input dtypes without truncation + or overflow. + + Parameters + ---------- + dtype_a : numpy.dtype + first input dtype + dtype_b : numpy.dtype + second input dtype + + Returns + ------- + numpy.dtype + promoted dtype + + """ + if dtype_a == dtype_b: + return dtype_a + if dtype_a.kind in "SU" and dtype_b.kind in "SU": + kind = "U" if "U" in (dtype_a.kind, dtype_b.kind) else "S" + size_a = dtype_a.itemsize // (4 if dtype_a.kind == "U" else 1) + size_b = dtype_b.itemsize // (4 if dtype_b.kind == "U" else 1) + return np.dtype(f"{kind}{max(size_a, size_b)}") + + return np.promote_types(dtype_a, dtype_b) + + +def group_dtype(group, keys, param_list=None): + """Group Dtype. + + Build the structured output dtype of a group of per-tile datasets, + promoting each column across *every* dataset. A campaign that was + partially reprocessed can carry e.g. ``S7`` tile IDs in one tile and + ``S12`` in another, or ``f4`` next to ``f8``; taking the dtype of the + first dataset alone would silently truncate the others. + + Parameters + ---------- + group : h5py.Group + group whose members are structured datasets + keys : list of str + dataset names to consider + param_list : list of str, optional + columns to keep; default is ``None`` (keep all) + + Returns + ------- + numpy.dtype + structured output dtype + + """ + names = param_list if param_list is not None else list(group[keys[0]].dtype.names) + + fields = [] + for name in names: + promoted = group[keys[0]].dtype[name] + for key in keys[1:]: + promoted = promote_dtypes(promoted, group[key].dtype[name]) + fields.append((name, promoted)) + + return np.dtype(fields) + + def _check_columns(dtype, param_list, file_path, dataset_key=None): """Raise a clear error if requested columns are absent from the data.""" missing = [col for col in param_list if col not in (dtype.names or ())] @@ -862,11 +928,7 @@ def concatenate_datasets(group, param_list=None, file_path="", verbose=True): if param_list is not None: _check_columns(group[key].dtype, param_list, file_path, key) - dtype_in = group[keys[0]].dtype - if param_list is not None: - dtype_out = np.dtype([(name, dtype_in[name]) for name in param_list]) - else: - dtype_out = dtype_in + dtype_out = group_dtype(group, keys, param_list=param_list) n_rows = sum(group[key].shape[0] for key in keys) if verbose: @@ -949,13 +1011,48 @@ def campaign_shape(file_path, param_list=None): if not keys: raise ValueError(f"No datasets found in catalogue {file_path}") n_rows = sum(group[key].shape[0] for key in keys) - dtype_in = group[keys[0]].dtype if param_list is not None: for key in keys: _check_columns(group[key].dtype, param_list, file_path, key) - dtype_in = np.dtype([(name, dtype_in[name]) for name in param_list]) + dtype_out = group_dtype(group, keys, param_list=param_list) + + return n_rows, dtype_out + + +def iter_campaign_tiles(file_path, param_list=None, verbose=True): + """Iter Campaign Tiles. + + Yield the per-tile datasets of a campaign catalogue one at a time, so a + caller filling a preallocated output array never holds more than one + tile in memory in addition to that output. - return n_rows, dtype_in + Parameters + ---------- + file_path : str + input file path + param_list : list of str, optional + columns to keep; default is ``None`` (keep all) + verbose : bool, optional + verbose output if ``True`` + + Yields + ------ + numpy.ndarray + one tile as a structured array + + """ + with h5py.File(file_path, "r") as hdf5_file: + group = find_dataset_group(hdf5_file) + check_n_tiles(hdf5_file, group, file_path) + keys = sorted(group) + if not keys: + raise ValueError(f"No datasets found in catalogue {file_path}") + for key in keys: + if param_list is not None: + _check_columns(group[key].dtype, param_list, file_path, key) + for key in tqdm.tqdm(keys, disable=not verbose): + data = group[key][()] + yield data if param_list is None else data[param_list] def read_campaign_catalogue( diff --git a/src/sp_validation/catalog_builders.py b/src/sp_validation/catalog_builders.py index 899d4db5..54759439 100644 --- a/src/sp_validation/catalog_builders.py +++ b/src/sp_validation/catalog_builders.py @@ -240,17 +240,8 @@ def close_hd5(self): self._hd5file.close() -def _promote(dtype_a, dtype_b): - """Return a dtype that holds both input dtypes without truncation.""" - if dtype_a == dtype_b: - return dtype_a - if dtype_a.kind in "SU" and dtype_b.kind in "SU": - kind = "U" if "U" in (dtype_a.kind, dtype_b.kind) else "S" - size_a = dtype_a.itemsize // (4 if dtype_a.kind == "U" else 1) - size_b = dtype_b.itemsize // (4 if dtype_b.kind == "U" else 1) - return np.dtype(f"{kind}{max(size_a, size_b)}") - - return np.promote_types(dtype_a, dtype_b) +# Column-dtype promotion is shared with the per-campaign reader in ``catalog``. +_promote = sp_cat.promote_dtypes def _checked_assign(target, start, end, values, name): @@ -510,8 +501,7 @@ def write_hdf5_file(self, dat_all, campaigns=None): with h5py.File(output_path, "w") as f: self.write_hdf5_header(f, campaigns=campaigns) - dset = f.create_dataset("data", data=dat_all) - dset[:] = dat_all + f.create_dataset("data", data=dat_all) def write_hdf5_header(self, hd5file, campaigns=None): """Write HDF5 Header. @@ -572,24 +562,27 @@ def merge_catalogues(self, input_paths): dat_all = np.empty(n_total, dtype=dtype_out) start = 0 for input_path, campaign in zip(input_paths, campaigns): - dat = sp_cat.read_campaign_catalogue( + # Fill tile by tile: peak memory is the merged output plus a + # single tile, never a whole campaign copy on top of it. + n_campaign = 0 + for dat in sp_cat.iter_campaign_tiles( input_path, param_list=param_list, verbose=self._params["verbose"], - ) - - end = start + len(dat) - for name in dat.dtype.names: - _checked_assign(dat_all[name], start, end, dat[name], name) - dat_all["campaign"][start:end] = campaign.encode() - start = end + ): + end = start + len(dat) + for name in dat.dtype.names: + _checked_assign(dat_all[name], start, end, dat[name], name) + dat_all["campaign"][start:end] = campaign.encode() + start = end + n_campaign += len(dat) + del dat if self._params["verbose"]: print( - f"{campaign}: added {len(dat)}" - + f" (~{format.millify(len(dat))}) objects." + f"{campaign}: added {n_campaign}" + + f" (~{format.millify(n_campaign)}) objects." ) - del dat if self._params["verbose"]: print( diff --git a/src/sp_validation/tests/test_campaign_readers.py b/src/sp_validation/tests/test_campaign_readers.py index 7272b4a7..be998876 100644 --- a/src/sp_validation/tests/test_campaign_readers.py +++ b/src/sp_validation/tests/test_campaign_readers.py @@ -161,6 +161,67 @@ def test_column_missing_from_later_tile_raises_clear_error(self): self.assertIn("MAG_AUTO", message) self.assertIn("001.000", message) + def test_dtype_promoted_across_tiles(self): + """A per-tile dtype difference within a campaign is not truncated.""" + narrow = np.zeros( + 2, dtype=[("N_EPOCH", "i2"), ("TILE_ID", "S7"), ("RA", "f4")] + ) + narrow["N_EPOCH"] = [1, 2] + narrow["TILE_ID"] = [b"123.456", b"123.457"] + narrow["RA"] = [1.5, 2.5] + + wide = np.zeros( + 2, dtype=[("N_EPOCH", "i4"), ("TILE_ID", "S12"), ("RA", "f8")] + ) + wide["N_EPOCH"] = [70000, 3] + wide["TILE_ID"] = [b"999888.7776", b"123.458"] + wide["RA"] = [3.123456789, 4.0] + + path = self._dir / "final_cat_MIX.hdf5" + write_campaign(path, "flat", {"000.000": narrow, "000.001": wide}) + param_list = ["N_EPOCH", "TILE_ID", "RA"] + + dat = catalog.read_campaign_catalogue( + str(path), param_list=param_list, verbose=False + ) + + self.assertEqual(dat.dtype["N_EPOCH"], np.dtype("i4")) + self.assertEqual(dat.dtype["TILE_ID"], np.dtype("S12")) + self.assertEqual(dat.dtype["RA"], np.dtype("f8")) + npt.assert_array_equal(dat["N_EPOCH"], [1, 2, 70000, 3]) + npt.assert_array_equal( + dat["TILE_ID"], + [b"123.456", b"123.457", b"999888.7776", b"123.458"], + ) + self.assertEqual(dat["RA"][2], 3.123456789) + + # campaign_shape must report the same promoted dtype, since the merge + # preallocates from it. + n_rows, dtype_out = catalog.campaign_shape( + str(path), param_list=param_list + ) + self.assertEqual(n_rows, 4) + self.assertEqual(dtype_out, dat.dtype) + + def test_iter_campaign_tiles(self): + """The streaming reader yields one restricted tile at a time.""" + path = self._dir / "final_cat_CAMPAIGN.hdf5" + write_campaign(path, "legacy", self._tiles) + + tiles = list( + catalog.iter_campaign_tiles( + str(path), param_list=["RA"], verbose=False + ) + ) + + self.assertEqual([len(tile) for tile in tiles], [3, 2]) + for tile in tiles: + self.assertEqual(tile.dtype.names, ("RA",)) + npt.assert_array_equal( + np.concatenate([tile["RA"] for tile in tiles]), + self._expected()["RA"], + ) + def test_ambiguous_layout_raises(self): path = self._dir / "ambiguous.hdf5" with h5py.File(path, "w") as f: From b2077ba01bb25af7ed8c2746ef2a31719d4b99b6 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 21:31:49 -0400 Subject: [PATCH 10/17] config: restore the r-band coverage cut in the v2 mask config mask_v2.0.yaml dropped v1.X's '64_r' r-band imaging cut without a replacement, so a v2 calibration run admitted objects outside the r-band footprint that v1 rejected -- a silent change of effective area, n(z) and galaxy-density normalisation. Enable MASK_n64, v1's '64_r' equivalent. v1's other coverage cut, npoint3 >= 3, came from an external post-processing catalogue and has no v2 counterpart; say so in the config rather than leaving its absence unexplained. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- config/calibration/mask_v2.0.yaml | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/config/calibration/mask_v2.0.yaml b/config/calibration/mask_v2.0.yaml index db824f59..201388eb 100644 --- a/config/calibration/mask_v2.0.yaml +++ b/config/calibration/mask_v2.0.yaml @@ -59,7 +59,15 @@ dat: kind: equal value: False - # Per-band coverage and Pan-STARRS z2; enable as required + # r-band coverage: the v2 equivalent of v1's '64_r' cut, kept so that the + # v2 selection covers the same imaging footprint as v1.X. + - col_name: MASK_n64 + label: "r-band imaging" + kind: equal + value: False + + # Remaining per-band coverage bits and Pan-STARRS z2; enable as required. + # Cutting on all of them at once leaves essentially no objects. # - col_name: MASK_n16 # label: "coverage (n16)" # kind: equal @@ -68,6 +76,12 @@ dat: # label: "no PS-z2" # kind: equal # value: False + # + # v1.X also applied a rough pointing-coverage cut, 'npoint3' >= 3, from an + # external post-processing catalogue. ShapePipe v2 emits no such column and + # no MASK_n* bit encodes it, so that cut has no v2 counterpart; add it back + # here if an external pointing-coverage map is ever joined onto the + # catalogue (it would belong in the 'dat_ext' group below). # Number of epochs - col_name: N_EPOCH From 204286edfe414946ca1bf3e89ac2cbee4c70d5f3 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 21:31:49 -0400 Subject: [PATCH 11/17] scripts: cheaper star mask cut, campaign-derived leakage labels extract_info fancy-indexed the full-width catalogue, dd[ind_star], only to read ~5 boolean mask columns from it -- a copy of every column for every matched star (~GB at DR6 scale). Mask first, index the resulting bool array. plot_leakage still labelled its curves "all", "P1" ... "P7": a P-named survivor of the #340 retirement, and a fixed length that silently mismatched the number of input files. Labels and colours now follow the input files. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- scripts/calibration/extract_info.py | 2 +- scripts/plot_leakage.py | 14 +++++++++++++- 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/scripts/calibration/extract_info.py b/scripts/calibration/extract_info.py index ed4fb145..d0994f07 100644 --- a/scripts/calibration/extract_info.py +++ b/scripts/calibration/extract_info.py @@ -160,7 +160,7 @@ m_star = ( (dd["FLAGS"][ind_star] == 0) - & galaxy.mask_cut(dd[ind_star], mask_columns) + & galaxy.mask_cut(dd, mask_columns)[ind_star] & (dd["NGMIX_MCAL_FLAGS"][ind_star] == 0) & (dd["NGMIX_G1_PSF_ORIG_NOSHEAR"][ind_star] != -10) ) diff --git a/scripts/plot_leakage.py b/scripts/plot_leakage.py index 629d81dc..c5b8f956 100755 --- a/scripts/plot_leakage.py +++ b/scripts/plot_leakage.py @@ -1,6 +1,7 @@ #!/usr/bin/env python3 import copy +import os import sys from optparse import OptionParser @@ -152,6 +153,7 @@ def plot_alpha_leakage( xmin, xmax, ylim=None, + labels=None, ): """Plot Alpha Leakage. @@ -173,6 +175,9 @@ def plot_alpha_leakage( largest angular scale, interpreted in arcmin ylim : list, optional y-axis plot limits, default is `Ǹone` + labels : list, optional + curve labels, one per input curve; default is ``None``, for + unlabelled curves """ theta = meanr @@ -187,7 +192,7 @@ def plot_alpha_leakage( linewidths[0] = 3 colors = ["grey", "k", "b", "r", "c", "m", "g", "orange"] - labels = ["all", "P1", "P2", "P3", "P4", "P5", "P6", "P7"] + colors = [colors[idx % len(colors)] for idx in range(len(meanr))] plot_data_1d( theta, @@ -239,6 +244,12 @@ def main(argv=None): if param.verbose: print("Input files: ", fnames) + # Curve labels come from the input file names: the first file is the + # reference (all objects), the others whatever selection they hold. + labels = ["all"] + [ + os.path.splitext(os.path.basename(fn))[0] for fn in fnames[1:] + ] + # read input files, append data theta = [] alpha_leak = [] @@ -265,6 +276,7 @@ def main(argv=None): config.theta_min_amin, config.theta_max_amin, config.leakage_alpha_ylim, + labels=labels, ) return 0 From d3741e6f21712a7ef3ab3c7e0c0c532176650877 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 21:55:17 -0400 Subject: [PATCH 12/17] masks: treat n64 as a reason bit, not r-band coverage The v2 branch inferred from v1's '64_r' naming that MASK_n64 flags r-band imaging coverage. It does not: bit 64 is an undocumented *reason* bit of the r-band default bitmask, and OR{n1,n2,n4,n8,n64,n1024} reproduces mask_r, the v1 r-band mask, exactly on the P3 region. Add MASK_n64 to DEFAULT_MASK_COLUMNS so the default galaxy cut is exactly that set, documented as reproducing mask_r, and mirror it in the calibration params, the v2 mask config and the minimal-catalogue demo. Document n16/n32/n128/n256 as the u/g/i/z coverage flags (no r flag: the catalogue is r-selected) and n2048 as absent Pan-STARRS z2. The faint vs bright assignment of n1/n2 is unconfirmed for the Aug-2026 products, so the labels no longer claim one. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- config/calibration/mask_v2.0.yaml | 34 +++++++++++++------ scripts/calibration/params.py | 4 ++- .../demo_comprehensive_to_minimal_cat.py | 1 + src/sp_validation/galaxy.py | 23 +++++++++---- .../tests/test_campaign_readers.py | 1 + 5 files changed, 45 insertions(+), 18 deletions(-) diff --git a/config/calibration/mask_v2.0.yaml b/config/calibration/mask_v2.0.yaml index 201388eb..37b0f5c7 100644 --- a/config/calibration/mask_v2.0.yaml +++ b/config/calibration/mask_v2.0.yaml @@ -4,16 +4,26 @@ # mask columns MASK_n, where True means MASKED. They are therefore cut # with `kind: equal, value: False` (keep the un-masked objects): # -# MASK_n1, MASK_n2 bright/faint star halos +# Reason bits of the ShapePipe r-band default bitmask: +# +# MASK_n1, MASK_n2 star halos (which is faint and which is bright is +# unconfirmed for the Aug-2026 products) # MASK_n4 stars # MASK_n8 manual galaxy mask -# MASK_n16 .. MASK_n256 per-band coverage +# MASK_n64 undocumented reason bit # MASK_n1024 MaxiMask +# +# Per-band coverage flags and Pan-STARRS: +# +# MASK_n16, MASK_n32, MASK_n128, MASK_n256 u, g, i, z coverage (no r +# flag: the catalogue is r-selected) # MASK_n2048 no Pan-STARRS z2 coverage # # The default selection (see sp_validation.galaxy.DEFAULT_MASK_COLUMNS) is -# n4 + n1 + n2 + n8 + n1024; the coverage bits and n2048 are listed here but -# commented out, since ORing all of them masks essentially everything. +# n1 + n2 + n4 + n8 + n64 + n1024, whose OR reproduces mask_r, the v1 r-band +# mask, exactly on the P3 region. The coverage flags and n2048 are listed +# here but commented out, since ORing all of them masks essentially +# everything. # General parameters (can also given on command line) params: @@ -42,12 +52,14 @@ dat: label: "stars" kind: equal value: False + # n1/n2 are the two star-halo bits; the faint/bright assignment is + # unconfirmed for the Aug-2026 products. - col_name: MASK_n1 - label: "faint star halos" + label: "star halos (n1)" kind: equal value: False - col_name: MASK_n2 - label: "bright star halos" + label: "star halos (n2)" kind: equal value: False - col_name: MASK_n8 @@ -59,17 +71,17 @@ dat: kind: equal value: False - # r-band coverage: the v2 equivalent of v1's '64_r' cut, kept so that the - # v2 selection covers the same imaging footprint as v1.X. + # Undocumented reason bit of the r-band default bitmask. Required for the + # OR above to reproduce mask_r; it is not a coverage flag. - col_name: MASK_n64 - label: "r-band imaging" + label: "reason bit n64" kind: equal value: False - # Remaining per-band coverage bits and Pan-STARRS z2; enable as required. + # Per-band coverage flags and Pan-STARRS z2; enable as required. # Cutting on all of them at once leaves essentially no objects. # - col_name: MASK_n16 - # label: "coverage (n16)" + # label: "u coverage" # kind: equal # value: False # - col_name: MASK_n2048 diff --git a/scripts/calibration/params.py b/scripts/calibration/params.py index 03504946..afb59560 100644 --- a/scripts/calibration/params.py +++ b/scripts/calibration/params.py @@ -122,12 +122,14 @@ "NGMIX_T_PSF_RECONV_NOSHEAR", ] -## ShapePipe v2 mask columns OR'd together for the galaxy selection cut +## ShapePipe v2 mask columns OR'd together for the galaxy selection cut: +## the reason bits of the r-band default bitmask, whose OR reproduces mask_r mask_columns = [ "MASK_n4", "MASK_n1", "MASK_n2", "MASK_n8", + "MASK_n64", "MASK_n1024", ] diff --git a/scripts/examples/demo_comprehensive_to_minimal_cat.py b/scripts/examples/demo_comprehensive_to_minimal_cat.py index 30b6b1df..5c99970e 100644 --- a/scripts/examples/demo_comprehensive_to_minimal_cat.py +++ b/scripts/examples/demo_comprehensive_to_minimal_cat.py @@ -58,6 +58,7 @@ "MASK_n1", "MASK_n2", "MASK_n8", + "MASK_n64", "MASK_n1024", "NGMIX_MCAL_TYPES_FAIL", "NGMIX_G1_PSF_ORIG_NOSHEAR", diff --git a/src/sp_validation/galaxy.py b/src/sp_validation/galaxy.py index 7a8b8b1c..df901b3e 100644 --- a/src/sp_validation/galaxy.py +++ b/src/sp_validation/galaxy.py @@ -29,9 +29,16 @@ from sp_validation import io #: All mask columns written by ShapePipe v2 (bool, ``True`` = masked). -#: n4 stars; n1/n2 faint/bright star halos; n8 manual galaxy mask; -#: n1024 MaxiMask; n16..n256 per-band coverage; n2048 no PS-z2 coverage. -#: These replace the single IMAFLAGS_ISO bitmask of ShapePipe v1. +#: They replace the single IMAFLAGS_ISO bitmask of ShapePipe v1. +#: +#: Reason bits making up the r-band default bitmask: n1/n2 star halos +#: (which of the two is faint and which bright is unconfirmed for the +#: Aug-2026 products), n4 stars, n8 manual galaxy mask, n64 (an +#: undocumented reason bit), n1024 MaxiMask. +#: +#: Per-band coverage flags: n16 (u), n32 (g), n128 (i), n256 (z). There is +#: no r coverage flag because the catalogue is r-selected. n2048 is ``True`` +#: where Pan-STARRS z2 coverage is absent. MASK_COLUMNS = ( "MASK_n1", "MASK_n2", @@ -46,14 +53,18 @@ "MASK_n2048", ) -#: Mask columns OR'd together for the default galaxy selection. Deliberately -#: not a blanket OR over MASK_COLUMNS: the per-band coverage columns -#: (n16..n256) and n2048 would mask essentially the whole catalogue. +#: Mask columns OR'd together for the default galaxy selection. This set is +#: exactly the reason bits of the ShapePipe r-band default bitmask: their OR +#: reproduces ``mask_r``, the v1 r-band mask, on the P3 region. Deliberately +#: not a blanket OR over MASK_COLUMNS: the per-band coverage flags +#: (n16, n32, n128, n256) and n2048 would mask essentially the whole +#: catalogue. DEFAULT_MASK_COLUMNS = ( "MASK_n4", "MASK_n1", "MASK_n2", "MASK_n8", + "MASK_n64", "MASK_n1024", ) diff --git a/src/sp_validation/tests/test_campaign_readers.py b/src/sp_validation/tests/test_campaign_readers.py index be998876..69374fdf 100644 --- a/src/sp_validation/tests/test_campaign_readers.py +++ b/src/sp_validation/tests/test_campaign_readers.py @@ -21,6 +21,7 @@ ("MASK_n1", "?"), ("MASK_n2", "?"), ("MASK_n8", "?"), + ("MASK_n64", "?"), ("MASK_n1024", "?"), ] ) From 76e2d7a43cfc5e95263baffc350e608960c467c2 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 10 Sep 2026 01:57:08 +0000 Subject: [PATCH 13/17] ruff autofix (format + safe lint fixes) Pushed by the lint gate. --- scripts/calibration/extract_info.py | 5 +-- scripts/combine_results.py | 6 ++-- scripts/glass_mock/compute_leakage_harmony.py | 1 - scripts/merge_psf_cat.py | 4 +-- scripts/plot_leakage.py | 4 +-- src/sp_validation/catalog.py | 4 ++- src/sp_validation/catalog_builders.py | 3 +- src/sp_validation/rho_tau.py | 4 ++- .../tests/test_campaign_readers.py | 35 +++++-------------- src/sp_validation/tests/test_survey.py | 2 -- 10 files changed, 22 insertions(+), 46 deletions(-) diff --git a/scripts/calibration/extract_info.py b/scripts/calibration/extract_info.py index d0994f07..d9de957e 100644 --- a/scripts/calibration/extract_info.py +++ b/scripts/calibration/extract_info.py @@ -34,7 +34,6 @@ import h5py import numpy as np -from astropy.io import fits # from sp_validation.catalog import * from sp_validation import catalog as spv_cat @@ -70,9 +69,7 @@ dd = np.load(galaxy_cat_path, mmap_mode=mmap_mode) else: print("Loading galaxy .hdf5 file...") - dd = spv_cat.read_campaign_catalogue( - galaxy_cat_path, param_path=param_list_path - ) + dd = spv_cat.read_campaign_catalogue(galaxy_cat_path, param_path=param_list_path) n_obj = len(dd) print_stats( diff --git a/scripts/combine_results.py b/scripts/combine_results.py index 1b3b4f85..5b679b98 100755 --- a/scripts/combine_results.py +++ b/scripts/combine_results.py @@ -281,9 +281,9 @@ def get_values(results, stats_files, shape, use_keys, area_deg2=-1): print(f"area({campaign}) = {area_deg2_campaign} deg^2") else: area_deg2_campaign = area_deg2 - results["value"][key_der][campaign] = results["value"]["N_gal"][campaign] / ( - area_deg2_campaign * 3600 - ) + results["value"][key_der][campaign] = results["value"]["N_gal"][ + campaign + ] / (area_deg2_campaign * 3600) if area_deg2 < 0: with open("area_deg2_tot.txt", "w") as f: diff --git a/scripts/glass_mock/compute_leakage_harmony.py b/scripts/glass_mock/compute_leakage_harmony.py index 22c2439d..302f822e 100644 --- a/scripts/glass_mock/compute_leakage_harmony.py +++ b/scripts/glass_mock/compute_leakage_harmony.py @@ -9,7 +9,6 @@ from astropy.io import fits from sp_validation import catalog as spv_cat - from sp_validation.glass_mock import compute_leakage_harmony diff --git a/scripts/merge_psf_cat.py b/scripts/merge_psf_cat.py index 8ca34cf4..5bf089f6 100644 --- a/scripts/merge_psf_cat.py +++ b/scripts/merge_psf_cat.py @@ -143,9 +143,7 @@ def merge_catalogues(self, campaigns): for name in col_names: dat_all[name] = np.append(dat_all[name], dat[name]) - dat_all["campaign"] = np.append( - dat_all["campaign"], [campaign] * len(dat) - ) + dat_all["campaign"] = np.append(dat_all["campaign"], [campaign] * len(dat)) col_names = col_names + ("campaign",) diff --git a/scripts/plot_leakage.py b/scripts/plot_leakage.py index c5b8f956..ea6194aa 100755 --- a/scripts/plot_leakage.py +++ b/scripts/plot_leakage.py @@ -246,9 +246,7 @@ def main(argv=None): # Curve labels come from the input file names: the first file is the # reference (all objects), the others whatever selection they hold. - labels = ["all"] + [ - os.path.splitext(os.path.basename(fn))[0] for fn in fnames[1:] - ] + labels = ["all"] + [os.path.splitext(os.path.basename(fn))[0] for fn in fnames[1:]] # read input files, append data theta = [] diff --git a/src/sp_validation/catalog.py b/src/sp_validation/catalog.py index 5f846ef4..037a2850 100644 --- a/src/sp_validation/catalog.py +++ b/src/sp_validation/catalog.py @@ -803,7 +803,9 @@ def find_dataset_group(hdf5_file): while True: keys = list(node) if not keys: - raise ValueError(f"No data found under {node.name!r} in {hdf5_file.file.filename}") + raise ValueError( + f"No data found under {node.name!r} in {hdf5_file.file.filename}" + ) if all(isinstance(node[key], h5py.Dataset) for key in keys): return node if len(keys) != 1: diff --git a/src/sp_validation/catalog_builders.py b/src/sp_validation/catalog_builders.py index 54759439..c0e58106 100644 --- a/src/sp_validation/catalog_builders.py +++ b/src/sp_validation/catalog_builders.py @@ -553,8 +553,7 @@ def merge_catalogues(self, input_paths): # merged array is allocated once and filled in place, instead of # concatenating per-campaign copies (peak memory 2x the output). shapes = [ - sp_cat.campaign_shape(path, param_list=param_list) - for path in input_paths + sp_cat.campaign_shape(path, param_list=param_list) for path in input_paths ] n_total = sum(n_rows for n_rows, _ in shapes) dtype_out = self.output_dtype([dtype for _, dtype in shapes], n_char_campaign) diff --git a/src/sp_validation/rho_tau.py b/src/sp_validation/rho_tau.py index 13be5849..3895a9e2 100644 --- a/src/sp_validation/rho_tau.py +++ b/src/sp_validation/rho_tau.py @@ -351,7 +351,9 @@ def get_jackknife_cov( tau_chunk = outdir + f"/cov_tau_{version}{i}.npy" rho_chunk = outdir + f"/cov_rho_{version}{i}.npy" if not (os.path.exists(tau_chunk) and os.path.exists(rho_chunk)): - print(f"Computing rho-statistics for {version} (jackknife realisation {i + 1}/{ncov})") + print( + f"Computing rho-statistics for {version} (jackknife realisation {i + 1}/{ncov})" + ) if f"psf_{version}{i}" not in rho_stat_handler.catalogs.catalogs_dict: # Build catalogues diff --git a/src/sp_validation/tests/test_campaign_readers.py b/src/sp_validation/tests/test_campaign_readers.py index 69374fdf..bf794e9d 100644 --- a/src/sp_validation/tests/test_campaign_readers.py +++ b/src/sp_validation/tests/test_campaign_readers.py @@ -1,6 +1,8 @@ """Tests for the ShapePipe v2 campaign catalogue readers and mask cut.""" +import tempfile import unittest +from pathlib import Path import h5py import numpy as np @@ -9,9 +11,6 @@ from sp_validation import catalog, galaxy from sp_validation.catalog_builders import JointCat -import tempfile -from pathlib import Path - GAL_DTYPE = np.dtype( [ ("RA", "f8"), @@ -114,7 +113,6 @@ def test_missing_column_raises(self): ) self.assertIn("NOT_A_COLUMN", str(ctx.exception)) - def test_row_order_is_tile_name_order(self): """Keys inserted out of order still concatenate in name order.""" path = self._dir / "unordered.hdf5" @@ -164,16 +162,12 @@ def test_column_missing_from_later_tile_raises_clear_error(self): def test_dtype_promoted_across_tiles(self): """A per-tile dtype difference within a campaign is not truncated.""" - narrow = np.zeros( - 2, dtype=[("N_EPOCH", "i2"), ("TILE_ID", "S7"), ("RA", "f4")] - ) + narrow = np.zeros(2, dtype=[("N_EPOCH", "i2"), ("TILE_ID", "S7"), ("RA", "f4")]) narrow["N_EPOCH"] = [1, 2] narrow["TILE_ID"] = [b"123.456", b"123.457"] narrow["RA"] = [1.5, 2.5] - wide = np.zeros( - 2, dtype=[("N_EPOCH", "i4"), ("TILE_ID", "S12"), ("RA", "f8")] - ) + wide = np.zeros(2, dtype=[("N_EPOCH", "i4"), ("TILE_ID", "S12"), ("RA", "f8")]) wide["N_EPOCH"] = [70000, 3] wide["TILE_ID"] = [b"999888.7776", b"123.458"] wide["RA"] = [3.123456789, 4.0] @@ -198,9 +192,7 @@ def test_dtype_promoted_across_tiles(self): # campaign_shape must report the same promoted dtype, since the merge # preallocates from it. - n_rows, dtype_out = catalog.campaign_shape( - str(path), param_list=param_list - ) + n_rows, dtype_out = catalog.campaign_shape(str(path), param_list=param_list) self.assertEqual(n_rows, 4) self.assertEqual(dtype_out, dat.dtype) @@ -210,9 +202,7 @@ def test_iter_campaign_tiles(self): write_campaign(path, "legacy", self._tiles) tiles = list( - catalog.iter_campaign_tiles( - str(path), param_list=["RA"], verbose=False - ) + catalog.iter_campaign_tiles(str(path), param_list=["RA"], verbose=False) ) self.assertEqual([len(tile) for tile in tiles], [3, 2]) @@ -339,21 +329,16 @@ def tearDown(self): self._tmp.cleanup() def test_campaign_name(self): - self.assertEqual( - JointCat.campaign_name("/some/dir/final_cat_W3.hdf5"), "W3" - ) + self.assertEqual(JointCat.campaign_name("/some/dir/final_cat_W3.hdf5"), "W3") def test_merge(self): dat = self._obj.merge_catalogues(self._paths) self.assertEqual(len(dat), 5) self.assertIn("campaign", dat.dtype.names) - npt.assert_array_equal( - dat["campaign"], np.array([b"W3"] * 3 + [b"SGC"] * 2) - ) + npt.assert_array_equal(dat["campaign"], np.array([b"W3"] * 3 + [b"SGC"] * 2)) npt.assert_array_equal(dat["RA"][:3], np.arange(3)) - def test_merge_promotes_column_widths(self): """A wider string/int column in a later file is not truncated.""" dtype_narrow = np.dtype([("RA", "f8"), ("TILE_ID", "S7"), ("N", "i4")]) @@ -373,9 +358,7 @@ def test_merge_promotes_column_widths(self): for paths in ([path_a, path_b], [path_b, path_a]): dat = self._obj.merge_catalogues([str(path) for path in paths]) - by_campaign = { - name: row for name, row in zip(dat["campaign"], dat) - } + by_campaign = {name: row for name, row in zip(dat["campaign"], dat)} self.assertEqual(by_campaign[b"BBBBBBBB"]["TILE_ID"], b"999888.7776") self.assertEqual(by_campaign[b"BBBBBBBB"]["N"], 2**40) self.assertEqual(by_campaign[b"AA"]["TILE_ID"], b"123.456") diff --git a/src/sp_validation/tests/test_survey.py b/src/sp_validation/tests/test_survey.py index c31f75cc..e786c8f8 100644 --- a/src/sp_validation/tests/test_survey.py +++ b/src/sp_validation/tests/test_survey.py @@ -28,7 +28,6 @@ def setUp(self): self._area_amin2 = 3600 self._tile_IDs = (270.283, 188.308) - def tearDown(self): self.number_tile = None @@ -54,4 +53,3 @@ def test_get_area(self): sorted(tile_IDs) == sorted(self._tile_IDs), msg=f"{tile_IDs}!={self._tile_IDs}", ) - From 1c199aee55eeaf3dfd15a887a473f013c16e9199 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 10 Sep 2026 00:49:47 -0400 Subject: [PATCH 14/17] galaxy: cut never-fit objects on NGMIX_N_EPOCH explicitly ShapePipe's make_cat pre-fills the NGMIX_* columns with sentinels (G1/G2 = -10, T/FLUX = 0) and overwrites them only for objects present in the ngmix output, so an object ngmix never fit keeps NGMIX_MCAL_FLAGS == 0 and passes a flag-only cut. In final_cat_smk-g7.hdf5 that is 18,983 of 1,851,100 objects (1.03%); admitting them drags mean e1 to -0.096 (std 0.98) from +0.0001. classification_galaxy_ngmix already rejected all 18,983 via the NGMIX_G1_PSF_ORIG_NOSHEAR != -10 guard, so the production selection was never affected -- verified on the real file, which gives 1,105,851 rows out with and without the new cut. But that protection was incidental: it is an exact float equality against a sentinel ShapePipe may change, and the coadd N_EPOCH >= 2 cut in classification_galaxy_base does not substitute for it (18,750 of the 18,983 have N_EPOCH >= 1). Cut on NGMIX_N_EPOCH > 0 explicitly so the guarantee is stated, not inferred. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- src/sp_validation/galaxy.py | 11 +++++++++ src/sp_validation/tests/test_galaxy.py | 34 ++++++++++++++++++++++++++ 2 files changed, 45 insertions(+) diff --git a/src/sp_validation/galaxy.py b/src/sp_validation/galaxy.py index df901b3e..5baf6c1a 100644 --- a/src/sp_validation/galaxy.py +++ b/src/sp_validation/galaxy.py @@ -292,8 +292,19 @@ def classification_galaxy_ngmix( Return mask corresponding to ngmix classification of galaxies """ + # NGMIX_N_EPOCH == 0 marks objects ngmix never fit: ShapePipe's make_cat + # pre-fills every NGMIX_* column with sentinels (G1/G2 = -10, T/FLUX = 0) + # and only overwrites them for objects present in the ngmix output, so a + # never-fit object keeps NGMIX_MCAL_FLAGS == 0 and passes a flag-only cut. + # In final_cat_smk-g7.hdf5 this is 18,983 / 1,851,100 objects (1.03%); + # admitting them drags mean e1 to -0.096 (std 0.98) from +0.0001 (std + # 0.24). The coadd N_EPOCH cut in classification_galaxy_base does not + # catch them (18,750 of the 18,983 have N_EPOCH >= 1). Cut on N_EPOCH + # explicitly rather than relying on the -10 sentinel comparison below, + # which is an exact float equality against a value ShapePipe may change. m_gal_ngmix = ( cut_common + & (dd["NGMIX_N_EPOCH"] > 0) & (dd["NGMIX_MCAL_FLAGS"] == 0) & (dd["NGMIX_G1_PSF_ORIG_NOSHEAR"] != -10) & (dd["NGMIX_MCAL_TYPES_FAIL"] == 0) diff --git a/src/sp_validation/tests/test_galaxy.py b/src/sp_validation/tests/test_galaxy.py index 529fd0ad..5ea6a2b7 100644 --- a/src/sp_validation/tests/test_galaxy.py +++ b/src/sp_validation/tests/test_galaxy.py @@ -37,3 +37,37 @@ def test_T_to_fwhm_is_dimensionally_correct(self): npt.assert_allclose(sigma_to_fwhm(1.0), 2.3548200450, rtol=1e-6) # the old linear form would give 4.0 here npt.assert_allclose(T_to_fwhm(2.0), sigma_to_fwhm(1.0)) + + def test_never_fit_objects_are_rejected(self): + """Test that NGMIX_N_EPOCH == 0 objects are cut. + + ShapePipe's make_cat pre-fills the NGMIX_* columns with + sentinels (G1/G2 = -10, T/FLUX = 0) and overwrites them only + for objects present in the ngmix output. An object ngmix never + fit therefore keeps NGMIX_MCAL_FLAGS == 0 and passes a + flag-only selection, carrying e1 = -10 into the shear + statistics. In final_cat_smk-g7.hdf5 that is 18,983 of + 1,851,100 objects (1.03%). + """ + import numpy as np + + from sp_validation.galaxy import classification_galaxy_ngmix + + # row 0: never fit (sentinels, flags clean); row 1: a good fit + dd = { + "NGMIX_N_EPOCH": np.array([0.0, 3.0]), + "NGMIX_MCAL_FLAGS": np.array([0.0, 0.0]), + "NGMIX_G1_PSF_ORIG_NOSHEAR": np.array([-10.0, 0.02]), + "NGMIX_MCAL_TYPES_FAIL": np.array([0.0, 0.0]), + "NGMIX_G1_NOSHEAR": np.array([-10.0, 0.1]), + } + cut_common = np.array([True, True]) + + keep = classification_galaxy_ngmix(dd, cut_common) + npt.assert_array_equal(keep, [False, True]) + + # the N_EPOCH cut must stand on its own: even if the -10 + # sentinel were to change, the never-fit row stays rejected + dd["NGMIX_G1_PSF_ORIG_NOSHEAR"] = np.array([-99.0, 0.02]) + keep = classification_galaxy_ngmix(dd, cut_common) + npt.assert_array_equal(keep, [False, True]) From de7a837b5a900015e5a70f35b61cbd5ad4e470d2 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 10 Sep 2026 00:55:21 -0400 Subject: [PATCH 15/17] readers: NaN-safe mask cut, star-catalogue provenance and count check mask_cut decided on truthiness via astype(bool). ShapePipe writes the MASK_n* columns as float64 {0, 1} rather than bool (being fixed upstream), and astype(bool) reads NaN as True, so an incomplete mask column would have silently deleted sky. Decide on the value instead (masked iff > 0.5), accept bool, int and float alike, and treat NaN as "no verdict recorded" -- keep the object, but count and warn, since a nonzero count means the product is defective. final_cat_smk-g7.hdf5 carries no NaNs and only exact 0.0/1.0, so this is defensive: the real file gives 1,105,851 rows out before and after. read_star_catalogue silently accepted a truncated file and threw away exposure provenance. It now validates the n_exposures root attribute against the datasets found, as the galaxy reader validates n_tiles (check_n_tiles generalised to check_n_units), and adds an EXPID column carrying the exposure number each star came from. The datasets are named by that number and concatenating them discarded it, leaving no way to group stars by exposure downstream. Names may be bare ("2086324", as smk-g7 writes them) or carry the CFIS suffix ("2110000p"), so EXPID takes the leading digits. On the real star catalogue: 53,264 stars over 127 exposures, matching n_exposures. Also note in group_dtype that ShapePipe writes TILE_ID as f8, so the string-promotion branch is for a future string-valued TILE_ID, with a TODO recording that as an open schema decision. No behaviour change. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- src/sp_validation/catalog.py | 107 ++++++++++++++---- src/sp_validation/galaxy.py | 30 ++++- .../tests/test_campaign_readers.py | 103 ++++++++++++++++- 3 files changed, 219 insertions(+), 21 deletions(-) diff --git a/src/sp_validation/catalog.py b/src/sp_validation/catalog.py index 037a2850..40d516d4 100644 --- a/src/sp_validation/catalog.py +++ b/src/sp_validation/catalog.py @@ -13,6 +13,7 @@ import getpass import os +import re import h5py import numpy as np @@ -855,6 +856,16 @@ def group_dtype(group, keys, param_list=None): ``S12`` in another, or ``f4`` next to ``f8``; taking the dtype of the first dataset alone would silently truncate the others. + Note that ShapePipe currently writes ``TILE_ID`` as ``f8`` (the tile + ``183.307`` arrives as the float ``183.307``), not as a string, so the + string-promotion branch above is for a future string-valued ``TILE_ID`` + and is not exercised by today's products. + + TODO: whether ``TILE_ID`` should be a string is an open schema decision. + A float cannot represent the ID exactly and cannot be compared for + equality safely; changing it is a breaking product change, so it is left + as-is here and this reader deliberately handles both. + Parameters ---------- group : h5py.Group @@ -893,7 +904,9 @@ def _check_columns(dtype, param_list, file_path, dataset_key=None): ) -def concatenate_datasets(group, param_list=None, file_path="", verbose=True): +def concatenate_datasets( + group, param_list=None, file_path="", verbose=True, key_column=None +): """Concatenate Datasets. Concatenate every dataset of an HDF5 group into one structured array, @@ -913,12 +926,23 @@ def concatenate_datasets(group, param_list=None, file_path="", verbose=True): input file path, for error messages verbose : bool, optional verbose output if ``True`` + key_column : str, optional + name of an extra integer column to add, filled per row with the + integer-valued name of the dataset the row came from. Concatenating + the group otherwise throws that name away; the star catalogue needs + it to keep exposure provenance. Default is ``None`` (add no column) Returns ------- numpy.ndarray concatenated structured array + Raises + ------ + ValueError + if ``key_column`` is given but a dataset name is not an integer, or + collides with an existing column + """ keys = sorted(group) if not keys: @@ -932,6 +956,25 @@ def concatenate_datasets(group, param_list=None, file_path="", verbose=True): dtype_out = group_dtype(group, keys, param_list=param_list) + if key_column is not None: + if key_column in (dtype_out.names or ()): + raise ValueError( + f"Cannot add column {key_column!r} to catalogue {file_path}:" + + " a column of that name is already present." + ) + # Dataset names are exposure numbers, either bare ("2086324", as the + # smk-g7 products write them) or carrying the CFIS processed-exposure + # suffix ("2110000p"), so take the leading run of digits. + matches = {key: re.match(r"\d+", key) for key in keys} + bad = [key for key, match in matches.items() if match is None] + if bad: + raise ValueError( + f"Cannot derive {key_column!r} for catalogue {file_path}:" + + f" dataset name(s) {bad} do not begin with an integer." + ) + key_values = {key: int(match.group()) for key, match in matches.items()} + dtype_out = np.dtype(dtype_out.descr + [(key_column, "i8")]) + n_rows = sum(group[key].shape[0] for key in keys) if verbose: print( @@ -946,43 +989,54 @@ def concatenate_datasets(group, param_list=None, file_path="", verbose=True): data = group[key][()] end = start + len(data) for name in dtype_out.names: - data_out[name][start:end] = data[name] + if key_column is not None and name == key_column: + data_out[name][start:end] = key_values[key] + else: + data_out[name][start:end] = data[name] start = end del data return data_out -def check_n_tiles(hdf5_file, group, file_path): - """Check Number of Tiles. +def check_n_units(hdf5_file, group, file_path, attr="n_tiles", unit="tile"): + """Check Number of Units. - Compare the number of datasets found against the ``n_tiles`` root - attribute the ShapePipe v2 product carries, to catch a catalogue that - was truncated by an interrupted merge job or file transfer. + Compare the number of datasets found against the count the ShapePipe v2 + product declares in a root attribute, to catch a catalogue that was + truncated by an interrupted merge job or file transfer. Galaxy + catalogues declare ``n_tiles`` and hold one dataset per tile; star + catalogues declare ``n_exposures`` and hold one per exposure. + + A missing attribute is not an error: older products carry none. Parameters ---------- hdf5_file : h5py.File open input file group : h5py.Group - group holding the per-tile datasets + group holding the per-unit datasets file_path : str input file path, for the error message + attr : str, optional + root attribute holding the declared count; default is ``n_tiles`` + unit : str, optional + name of one unit, for the error message; default is ``tile`` Raises ------ ValueError - if the number of datasets differs from the ``n_tiles`` attribute + if the number of datasets differs from the declared count """ - n_tiles = hdf5_file.attrs.get("n_tiles") - if n_tiles is None: + n_declared = hdf5_file.attrs.get(attr) + if n_declared is None: return n_found = len(group) - if int(n_tiles) != n_found: + if int(n_declared) != n_found: raise ValueError( - f"Catalogue {file_path} declares n_tiles = {int(n_tiles)} but holds" - + f" {n_found} tile dataset(s); the file is incomplete." + f"Catalogue {file_path} declares {attr} = {int(n_declared)} but" + + f" holds {n_found} {unit} dataset(s); the file is incomplete." ) @@ -1008,7 +1062,7 @@ def campaign_shape(file_path, param_list=None): """ with h5py.File(file_path, "r") as hdf5_file: group = find_dataset_group(hdf5_file) - check_n_tiles(hdf5_file, group, file_path) + check_n_units(hdf5_file, group, file_path) keys = sorted(group) if not keys: raise ValueError(f"No datasets found in catalogue {file_path}") @@ -1045,7 +1099,7 @@ def iter_campaign_tiles(file_path, param_list=None, verbose=True): """ with h5py.File(file_path, "r") as hdf5_file: group = find_dataset_group(hdf5_file) - check_n_tiles(hdf5_file, group, file_path) + check_n_units(hdf5_file, group, file_path) keys = sorted(group) if not keys: raise ValueError(f"No datasets found in catalogue {file_path}") @@ -1090,7 +1144,7 @@ def read_campaign_catalogue( with h5py.File(file_path, "r") as hdf5_file: group = find_dataset_group(hdf5_file) - check_n_tiles(hdf5_file, group, file_path) + check_n_units(hdf5_file, group, file_path) return concatenate_datasets( group, param_list=param_list, file_path=file_path, verbose=verbose ) @@ -1103,6 +1157,16 @@ def read_star_catalogue(file_path, hdu=1, verbose=True): ``full_starcat_.hdf5`` (one dataset per exposure), or a legacy FITS star catalogue when ``file_path`` ends in ``.fits``. + The HDF5 path adds an ``EXPID`` column carrying the exposure number each + star came from. The datasets are named by that number and concatenating + them would otherwise discard it, leaving no way to group stars by + exposure downstream (e.g. for per-exposure PSF residuals). + + The dataset count is validated against the ``n_exposures`` root + attribute, as the galaxy reader validates ``n_tiles``, so a catalogue + truncated by an interrupted merge is caught on read rather than showing + up as a quietly short star sample. + Parameters ---------- file_path : str @@ -1115,7 +1179,7 @@ def read_star_catalogue(file_path, hdu=1, verbose=True): Returns ------- numpy.ndarray - star catalogue data + star catalogue data, with an added ``EXPID`` column on the HDF5 path """ if str(file_path).endswith(".fits"): @@ -1123,7 +1187,12 @@ def read_star_catalogue(file_path, hdu=1, verbose=True): with h5py.File(file_path, "r") as hdf5_file: group = find_dataset_group(hdf5_file) - return concatenate_datasets(group, file_path=file_path, verbose=verbose) + check_n_units( + hdf5_file, group, file_path, attr="n_exposures", unit="exposure" + ) + return concatenate_datasets( + group, file_path=file_path, verbose=verbose, key_column="EXPID" + ) def get_maked_col(dat, col, mask): diff --git a/src/sp_validation/galaxy.py b/src/sp_validation/galaxy.py index 5baf6c1a..60e170fa 100644 --- a/src/sp_validation/galaxy.py +++ b/src/sp_validation/galaxy.py @@ -11,6 +11,7 @@ """ import re +import warnings import numpy as np import regions @@ -119,8 +120,35 @@ def mask_cut(dd, mask_columns=None): ) masked = np.zeros(len(dd[columns[0]]), dtype=bool) + n_undefined = 0 for col in columns: - masked |= np.asarray(dd[col], dtype=bool) + values = np.asarray(dd[col]) + if values.dtype == bool: + flagged = values + else: + # ShapePipe's writer emits the MASK_n* columns as float64 {0, 1} + # rather than bool (being fixed upstream), so decide on the value + # rather than on truthiness: a bare astype(bool) would silently + # read NaN as True, i.e. masked. Threshold at 0.5 so an integer, + # a float and a bool column all behave identically. + values = values.astype(float) + undefined = np.isnan(values) + n_undefined += int(undefined.sum()) + flagged = np.where(undefined, False, values > 0.5) + masked |= flagged + + if n_undefined: + # NaN means the masking stage recorded no verdict for this object. + # Treat it as un-masked (keep the object) so an incomplete mask + # column cannot silently delete sky, but say so loudly: a nonzero + # count here means the input product is defective. + warnings.warn( + f"{n_undefined} NaN value(s) in mask column(s) {columns};" + + " treated as not masked. The mask columns of a complete" + + " ShapePipe product hold only 0 and 1.", + RuntimeWarning, + stacklevel=2, + ) return ~masked diff --git a/src/sp_validation/tests/test_campaign_readers.py b/src/sp_validation/tests/test_campaign_readers.py index bf794e9d..6b22de99 100644 --- a/src/sp_validation/tests/test_campaign_readers.py +++ b/src/sp_validation/tests/test_campaign_readers.py @@ -238,6 +238,7 @@ def setUp(self): self._path = self._dir / "full_starcat_CAMPAIGN.hdf5" with h5py.File(self._path, "w") as f: + f.attrs["n_exposures"] = len(self._exposures) group = f.create_group("exposures") for exp, dat in self._exposures.items(): group.create_dataset(exp, data=dat) @@ -249,10 +250,61 @@ def test_read_hdf5(self): dat = catalog.read_star_catalogue(self._path, verbose=False) self.assertEqual(len(dat), 10) - self.assertEqual(tuple(dat.dtype.names), catalog.STAR_CAT_COLUMNS) + # the reader appends EXPID to the ShapePipe star columns + self.assertEqual( + tuple(dat.dtype.names), catalog.STAR_CAT_COLUMNS + ("EXPID",) + ) npt.assert_array_equal(dat["MAG"][:4], np.zeros(4)) npt.assert_array_equal(dat["MAG"][4:], np.ones(6)) + def test_expid_records_exposure_provenance(self): + """Test that EXPID keeps the exposure each star came from. + + The star catalogue holds one dataset per exposure, named by the + exposure number; concatenating them otherwise throws that number + away, leaving no way to group stars by exposure downstream (for + per-exposure PSF residuals, say). The name may be bare + ("2086324", as the smk-g7 products write it) or carry the CFIS + processed-exposure suffix ("2110000p"). + """ + dat = catalog.read_star_catalogue(self._path, verbose=False) + + self.assertEqual(dat.dtype["EXPID"].kind, "i") + npt.assert_array_equal(dat["EXPID"][:4], np.full(4, 2110000)) + npt.assert_array_equal(dat["EXPID"][4:], np.full(6, 2110001)) + + def test_truncated_star_catalogue_is_rejected(self): + """Test that n_exposures is validated against the datasets found. + + A merge job killed part-way through leaves a readable file with + fewer exposures than it declares. Without this check that shows + up only as a quietly short star sample, never as an error. The + galaxy reader already validates n_tiles this way. + """ + path = self._dir / "truncated.hdf5" + with h5py.File(path, "w") as f: + f.attrs["n_exposures"] = 7 # but only two datasets written + group = f.create_group("exposures") + for exp, dat in self._exposures.items(): + group.create_dataset(exp, data=dat) + + with self.assertRaises(ValueError) as ctx: + catalog.read_star_catalogue(path, verbose=False) + self.assertIn("n_exposures", str(ctx.exception)) + + def test_missing_n_exposures_attr_is_allowed(self): + """Test that a product carrying no n_exposures still reads. + + Older products declare no count; that is not a defect. + """ + path = self._dir / "no_attr.hdf5" + with h5py.File(path, "w") as f: + group = f.create_group("exposures") + for exp, dat in self._exposures.items(): + group.create_dataset(exp, data=dat) + + self.assertEqual(len(catalog.read_star_catalogue(path, verbose=False)), 10) + def test_read_fits(self): from astropy.io import fits @@ -308,6 +360,55 @@ def test_v1_catalogue_raises(self): with self.assertRaises(KeyError): galaxy.mask_cut(dat) + def test_float_and_int_mask_columns(self): + """Test that float and int mask columns cut like bool ones. + + ShapePipe's writer emits the real MASK_n* columns as float64 + {0.0, 1.0} rather than bool (being fixed upstream), so the cut + must decide on the value, not on the dtype. + """ + for dtype in ("f8", "i4"): + dat = np.zeros( + 4, dtype=[(col, dtype) for col in galaxy.DEFAULT_MASK_COLUMNS] + ) + dat["MASK_n4"][0] = 1 + dat["MASK_n1024"][2] = 1 + + npt.assert_array_equal( + galaxy.mask_cut(dat), + np.array([False, True, False, True]), + err_msg=f"mask column dtype {dtype}", + ) + + def test_nan_mask_value_is_not_masked_and_warns(self): + """Test that a NaN mask value keeps the object, loudly. + + A bare astype(bool) reads NaN as True, i.e. masked, so an + incomplete mask column would silently delete sky. NaN means the + masking stage recorded no verdict, so keep the object and warn: + a nonzero count means the input product is defective. The real + smk-g7 catalogue carries no NaNs, so this is defensive. + """ + dat = np.zeros(3, dtype=[(col, "f8") for col in galaxy.DEFAULT_MASK_COLUMNS]) + dat["MASK_n4"][0] = np.nan + dat["MASK_n4"][1] = 1.0 + + with self.assertWarns(RuntimeWarning) as ctx: + keep = galaxy.mask_cut(dat) + + # row 0 NaN -> kept, row 1 masked, row 2 clean -> kept + npt.assert_array_equal(keep, np.array([True, False, True])) + self.assertIn("NaN", str(ctx.warning)) + + def test_no_warning_when_no_nan(self): + dat = np.zeros(2, dtype=[(col, "f8") for col in galaxy.DEFAULT_MASK_COLUMNS]) + + import warnings + + with warnings.catch_warnings(): + warnings.simplefilter("error", RuntimeWarning) + npt.assert_array_equal(galaxy.mask_cut(dat), np.array([True, True])) + class TestCampaignMerge(unittest.TestCase): """Merge of several campaign catalogues into a joint catalogue.""" From 5313efac17f0f6dffc7f2e4a698375b653514ff8 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 10 Sep 2026 04:56:41 +0000 Subject: [PATCH 16/17] ruff autofix (format + safe lint fixes) Pushed by the lint gate. --- src/sp_validation/catalog.py | 4 +--- src/sp_validation/tests/test_campaign_readers.py | 4 +--- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/sp_validation/catalog.py b/src/sp_validation/catalog.py index 40d516d4..6ea273bb 100644 --- a/src/sp_validation/catalog.py +++ b/src/sp_validation/catalog.py @@ -1187,9 +1187,7 @@ def read_star_catalogue(file_path, hdu=1, verbose=True): with h5py.File(file_path, "r") as hdf5_file: group = find_dataset_group(hdf5_file) - check_n_units( - hdf5_file, group, file_path, attr="n_exposures", unit="exposure" - ) + check_n_units(hdf5_file, group, file_path, attr="n_exposures", unit="exposure") return concatenate_datasets( group, file_path=file_path, verbose=verbose, key_column="EXPID" ) diff --git a/src/sp_validation/tests/test_campaign_readers.py b/src/sp_validation/tests/test_campaign_readers.py index 6b22de99..d30431d1 100644 --- a/src/sp_validation/tests/test_campaign_readers.py +++ b/src/sp_validation/tests/test_campaign_readers.py @@ -251,9 +251,7 @@ def test_read_hdf5(self): self.assertEqual(len(dat), 10) # the reader appends EXPID to the ShapePipe star columns - self.assertEqual( - tuple(dat.dtype.names), catalog.STAR_CAT_COLUMNS + ("EXPID",) - ) + self.assertEqual(tuple(dat.dtype.names), catalog.STAR_CAT_COLUMNS + ("EXPID",)) npt.assert_array_equal(dat["MAG"][:4], np.zeros(4)) npt.assert_array_equal(dat["MAG"][4:], np.ones(6)) From b2d8e32934bec9b92c786009cc0c28d9db59e3eb Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 10 Sep 2026 00:57:50 -0400 Subject: [PATCH 17/17] config: cut never-fit objects in the v2 mask config too mask_v2.0.yaml, applied downstream to the comprehensive catalogue, cut on neither NGMIX_MCAL_FLAGS nor any epoch column. It rejected the never-fit objects only through its NGMIX_G1/G2_PSF_ORIG_NOSHEAR != -10 cuts, and only because make_cat happens to fill the PSF columns from the same -10 literal it uses for the galaxy ellipticities. That is the same accidental immunity just removed from classification_galaxy_ngmix, one stage further downstream. Add NGMIX_N_EPOCH >= 1. The column is already carried into the comprehensive catalogue via add_cols_pre_cal in params.py. On final_cat_smk-g7.hdf5 the cut keeps 1,832,117 of 1,851,100 objects, removing exactly the 18,983 (1.03%) never-fit rows. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- config/calibration/mask_v2.0.yaml | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/config/calibration/mask_v2.0.yaml b/config/calibration/mask_v2.0.yaml index 37b0f5c7..83aa72b1 100644 --- a/config/calibration/mask_v2.0.yaml +++ b/config/calibration/mask_v2.0.yaml @@ -113,6 +113,22 @@ dat: kind: equal value: 0 + # Objects ngmix never fit. ShapePipe's make_cat pre-fills the NGMIX_* + # columns with sentinels (G1/G2 = -10, T/FLUX = 0) and overwrites them only + # for objects present in the ngmix output, so a never-fit object also keeps + # NGMIX_MCAL_FLAGS = 0 and is not caught by any flag cut. In + # final_cat_smk-g7.hdf5 that is 18,983 of 1,851,100 objects (1.03%); + # admitting them gives mean e1 = -0.107 (std 1.03) against -0.004 (std + # 0.22). The -10 cuts below do reject them today, but only because make_cat + # happens to use the same literal for the PSF columns: that is an exact + # float equality against a sentinel ShapePipe may change, so state the + # condition directly. Mirrors the guard in + # sp_validation.galaxy.classification_galaxy_ngmix. + - col_name: NGMIX_N_EPOCH + label: "ngmix never fit" + kind: greater_equal + value: 1 + # invalid PSF ellipticities - col_name: NGMIX_G1_PSF_ORIG_NOSHEAR label: "bad PSF ellipticity comp 1"