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/config/calibration/mask_v2.0.yaml b/config/calibration/mask_v2.0.yaml new file mode 100644 index 00000000..83aa72b1 --- /dev/null +++ b/config/calibration/mask_v2.0.yaml @@ -0,0 +1,164 @@ +# 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): +# +# 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_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 +# 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: + 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 + # n1/n2 are the two star-halo bits; the faint/bright assignment is + # unconfirmed for the Aug-2026 products. + - col_name: MASK_n1 + label: "star halos (n1)" + kind: equal + value: False + - col_name: MASK_n2 + label: "star halos (n2)" + 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 + + # 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: "reason bit n64" + kind: equal + value: False + + # 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: "u coverage" + # kind: equal + # value: False + # - col_name: MASK_n2048 + # 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 + 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 + + # 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" + 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/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/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/calibration/extract_info.py b/scripts/calibration/extract_info.py index 123c7437..d9de957e 100644 --- a/scripts/calibration/extract_info.py +++ b/scripts/calibration/extract_info.py @@ -34,10 +34,10 @@ import h5py import numpy as np -from astropy.io import fits # 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,9 +69,7 @@ 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) print_stats( @@ -116,7 +114,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) @@ -141,14 +139,13 @@ # #### 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"], [col_name_ra, col_name_dec], thresh, stats_file, - name=None, verbose=verbose, ) @@ -160,7 +157,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 +308,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( @@ -490,7 +488,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 ea871119..afb59560 100644 --- a/scripts/calibration/params.py +++ b/scripts/calibration/params.py @@ -25,9 +25,9 @@ # 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); required, it names +## the ShapePipe v2 products (final_cat_.hdf5, ...) +campaign = "W3" ## Area of a tile in deg^2 area_tile = 0.25 @@ -44,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 @@ -121,11 +122,32 @@ "NGMIX_T_PSF_RECONV_NOSHEAR", ] +## 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", +] + ## 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 +163,6 @@ add_cols_pre_cal_format = {} for key in ( "NUMBER", - "IMAFLAGS_ISO", "FLAGS", "NGMIX_MCAL_FLAGS", "NGMIX_MCAL_TYPES_FAIL", @@ -150,6 +171,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/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..5b679b98 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"): @@ -217,47 +217,55 @@ 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 patch 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): """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,16 +274,16 @@ 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: with open("area_deg2_tot.txt", "w") as f: @@ -285,10 +293,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 +308,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 +325,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 +342,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 +359,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 +376,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 +401,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 +432,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 +493,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 +554,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 +622,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 +695,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 +722,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 +732,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 +775,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 +796,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/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/examples/demo_comprehensive_to_minimal_cat.py b/scripts/examples/demo_comprehensive_to_minimal_cat.py index 11e5ed5a..5c99970e 100644 --- a/scripts/examples/demo_comprehensive_to_minimal_cat.py +++ b/scripts/examples/demo_comprehensive_to_minimal_cat.py @@ -50,13 +50,19 @@ # + # 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_n64", + "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/scripts/glass_mock/compute_leakage_harmony.py b/scripts/glass_mock/compute_leakage_harmony.py index a63c6bd4..302f822e 100644 --- a/scripts/glass_mock/compute_leakage_harmony.py +++ b/scripts/glass_mock/compute_leakage_harmony.py @@ -8,6 +8,7 @@ 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 +54,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/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/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/scripts/merge_psf_cat.py b/scripts/merge_psf_cat.py index 65c1e891..5bf089f6 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,21 +138,24 @@ 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"], [campaign] * 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" + # 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) @@ -177,12 +169,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_leakage.py b/scripts/plot_leakage.py index 629d81dc..ea6194aa 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,10 @@ 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 +274,7 @@ def main(argv=None): config.theta_min_amin, config.theta_max_amin, config.leakage_alpha_ylim, + labels=labels, ) return 0 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 e39d2076..6ea273bb 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 @@ -23,7 +24,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 +154,6 @@ def check_matching( keys_2, thresh, stats_file, - name=None, verbose=False, ): """Check matching. @@ -178,28 +177,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 """ - 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)) - # 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%}" @@ -215,7 +206,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): @@ -758,71 +749,448 @@ 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. + + 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 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. + + 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 + 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 ())] + 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}{where}." + + f" Available columns: {sorted(dtype.names or ())}" + ) + + +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, + 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 + 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`` + 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: + 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_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( + f"Reading {len(keys)} datasets," + + 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_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: + 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_units(hdf5_file, group, file_path, attr="n_tiles", unit="tile"): + """Check Number of Units. + + 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-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 declared count + + """ + n_declared = hdf5_file.attrs.get(attr) + if n_declared is None: + return + n_found = len(group) + if int(n_declared) != n_found: + raise ValueError( + f"Catalogue {file_path} declares {attr} = {int(n_declared)} but" + + f" holds {n_found} {unit} dataset(s); the file is incomplete." + ) + - Read hdf5 file and return contained data. +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 - 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_list : list of str, optional + columns to keep; default is ``None`` (keep all) Returns ------- - dict - data + tuple + number of rows (int) and dtype (numpy.dtype) """ - param_list = read_param_file(param_path, verbose=True) if param_path else None + with h5py.File(file_path, "r") as hdf5_file: + group = find_dataset_group(hdf5_file) + check_n_units(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) + if param_list is not None: + for key in keys: + _check_columns(group[key].dtype, param_list, file_path, key) + 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. + + 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: - # 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}"] + group = find_dataset_group(hdf5_file) + check_n_units(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( + file_path, + param_path=None, + param_list=None, + verbose=True, +): + """Read Campaign Catalogue. - # 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) + Read a campaign galaxy catalogue (``final_cat_.hdf5``) and + return its per-tile datasets concatenated into one structured array. - print( - f"Estimating {num_cols * num_rows * 8 / 1024**3:.1f}" - + f" Gb memory for the ({num_cols} x {num_rows}) data array ..." + Parameters + ---------- + file_path : str + input file path + 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 + ------- + numpy.ndarray + catalogue data + + """ + 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: + group = find_dataset_group(hdf5_file) + check_n_units(hdf5_file, group, file_path) + 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. + + 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``. - if not check_only: - # Add new to existing data - data_list.append(data) + 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). - print("Combine tile catalogues") - data_comb = np.concatenate(data_list, axis=0) - print("Done") + 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. - # Print problematic tile IDs - for ID in ID_pbl: - print("Tile IDs with missing keys:", file=stats_file) - print(ID, file=stats_file) + 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`` - return data_comb + Returns + ------- + numpy.ndarray + star catalogue data, with an added ``EXPID`` column on the HDF5 path + + """ + if str(file_path).endswith(".fits"): + return fits.getdata(file_path, hdu) + + 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") + 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/catalog_builders.py b/src/sp_validation/catalog_builders.py index 78dc7838..c0e58106 100644 --- a/src/sp_validation/catalog_builders.py +++ b/src/sp_validation/catalog_builders.py @@ -240,6 +240,28 @@ def close_hd5(self): self._hd5file.close() +# 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): + """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. @@ -279,147 +301,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 - - """ - 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 + list of str + input file paths """ - 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,77 +401,85 @@ def dtype_out(self, name, dtype_in): cols_keep_dtype = [ "RA", "Dec", + "DEC", "FLAGS", - "IMAFLAGS_ISO", "NUMBER", ] if dtype_in.kind == "U": # 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 + 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 init_data(self, n_col, n_obj, ndim, dat): - """Init Data. + def output_dtype(self, dtypes_in, n_char_campaign): + """Output Dtype. - Initialize empty structured data. + Return the merged-catalogue dtype: the input columns (possibly + reduced in precision, and promoted to a common type across all input + campaigns) 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 + 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 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) + Raises + ------ + ValueError + if the inputs have different column sets, or a column is + multi-dimensional (campaign catalogues are scalar-column only) - 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="", - ) + """ + if isinstance(dtypes_in, np.dtype): + dtypes_in = [dtypes_in] - dat_all = np.empty((n_obj,), dtype=dtype_tmp_struct) + 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)}" + ) - if self._params["verbose"]: - print("done") + 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 dat_all + return np.dtype(fields) - 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. @@ -521,8 +488,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 = ( @@ -532,12 +499,11 @@ 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 + f.create_dataset("data", data=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. @@ -546,92 +512,84 @@ 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 + # 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): + # 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 + n_campaign += len(dat) + del dat - 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: - raise ValueError( - "Inconsistent number of columns, {i_col + 1}" + f" != {n_col}" - ) if self._params["verbose"]: print( - f"{patch}: Added {len(dat)} (~{format.millify(len(dat))})" - + f" objects (from {start} to {end - 1})." + f"{campaign}: added {n_campaign}" + + f" (~{format.millify(n_campaign)}) objects." ) - start = end - del dat + if self._params["verbose"]: + print( + f"Merged {len(dat_all)} (~{format.millify(len(dat_all))})" + + f" objects from {len(campaigns)} campaign(s)." + ) - self.write_hdf5_file(dat_all, patches) + return dat_all def run(self): """Run. @@ -639,11 +597,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): @@ -1048,8 +1007,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) @@ -1125,7 +1084,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": diff --git a/src/sp_validation/galaxy.py b/src/sp_validation/galaxy.py index 3c49a96b..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 @@ -28,6 +29,129 @@ # required square root: FWHM = 2.35482 sqrt(T / 2) from sp_validation import io +#: All mask columns written by ShapePipe v2 (bool, ``True`` = masked). +#: 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", + "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. 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", +) + + +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) + 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] + 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) + n_undefined = 0 + for col in columns: + 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 + def classification_galaxy_overlap_ra_dec(dd, ra_key="XWIN_WORLD", dec_key="YWIN_WORLD"): """Classification Galaxy Overlap Ra Dec. @@ -141,11 +265,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 +303,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) ) @@ -189,8 +320,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/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 diff --git a/src/sp_validation/rho_tau.py b/src/sp_validation/rho_tau.py index 6874a08c..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} (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_campaign_readers.py b/src/sp_validation/tests/test_campaign_readers.py new file mode 100644 index 00000000..d30431d1 --- /dev/null +++ b/src/sp_validation/tests/test_campaign_readers.py @@ -0,0 +1,512 @@ +"""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 +import numpy.testing as npt + +from sp_validation import catalog, galaxy +from sp_validation.catalog_builders import JointCat + +GAL_DTYPE = np.dtype( + [ + ("RA", "f8"), + ("Dec", "f8"), + ("MAG_AUTO", "f4"), + ("MASK_n4", "?"), + ("MASK_n1", "?"), + ("MASK_n2", "?"), + ("MASK_n8", "?"), + ("MASK_n64", "?"), + ("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): + """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" + write_campaign(path, "legacy", self._tiles) + + dat = catalog.read_campaign_catalogue(path, verbose=False) + + self.assertEqual(len(dat), 5) + npt.assert_array_equal(dat["RA"], 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(dat["RA"], 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_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_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: + 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: + 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) + + 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) + # 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 + + 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) + + 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.""" + + 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_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")]) + 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() 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]) diff --git a/src/sp_validation/tests/test_survey.py b/src/sp_validation/tests/test_survey.py index f42f884c..e786c8f8 100644 --- a/src/sp_validation/tests/test_survey.py +++ b/src/sp_validation/tests/test_survey.py @@ -28,10 +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): self.number_tile = None @@ -57,9 +53,3 @@ def test_get_area(self): sorted(tile_IDs) == sorted(self._tile_IDs), 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..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 or patch 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 @@ -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",