diff --git a/astra.yaml b/astra.yaml index f339eb4f4..24dd350d3 100644 --- a/astra.yaml +++ b/astra.yaml @@ -329,10 +329,11 @@ analyses: # ═════════════════════════════════════════════════════════════════════════ detection: description: >- - Object detection with SExtractor on r-band tiles (the galaxy sample) and - on single-exposure CCDs (PSF-star candidates). Tiles follow the MegaPipe - parameters of Gwyn's UNIONS tile catalogue; exposures keep ShapePipe's - stock values. + Object detection on r-band tiles (the galaxy sample) and on + single-exposure CCDs (PSF-star candidates, always SExtractor). The tile + sample is Gwyn's UNIONS tile catalogue or ShapePipe's own SExtractor run + (tile_detection); the latter follows the MegaPipe parameters of that + catalogue, while exposures keep ShapePipe's stock values. inputs: - id: tile_stack type: data @@ -344,12 +345,34 @@ analyses: description: Per-tile SExtractor LDAC catalogue with per-epoch CCD membership. inputs: [tile_stack] decisions: - [detection_threshold_policy, deblending_policy, background_model, + [tile_detection, detection_threshold_policy, deblending_policy, + background_model, weight_map_usage, zero_weight_interpolation, detection_source_mode, - epoch_membership_ccd_bounds, photometry_parameters, + epoch_membership_ccd_bounds, catalogue_neighbour_marking, + photometry_parameters, spurious_detection_cleaning, blend_photometry_mask_type, saturation_level] decisions: + tile_detection: + label: Source of the tile galaxy sample + rationale: >- + Real data takes its tile sample from Gwyn's UNIONS per-tile + catalogue, converted to the sexcat the chain reads with its own + NUMBER kept, so shape, photometry and photo-z catalogues share one + object list and one ID (TILE_UNIQUE_ID). PSF stars still come from + ShapePipe's exposure-level SExtractor run, because external star + catalogues are too shallow. detection_source_mode and the tile + SExtractor settings of the other decisions here govern only the + sextractor option; epoch_membership_ccd_bounds applies to both. + Image simulations, which have no UNIONS catalogue, use sextractor. + Values: + input_types.data.tile_detection = unions_catalogue. + default: unions_catalogue + options: + unions_catalogue: + label: Gwyn's UNIONS per-tile catalogue, converted in place + sextractor: + label: SExtractor on the tile image detection_threshold_policy: label: Detection significance, minimum area, matched filter rationale: >- @@ -599,15 +622,48 @@ analyses: admits its upper endpoint, column 2080. A WCS inversion failure skips that CCD, lowering N_EPOCH; this sets how many exposures enter each galaxy's multi-epoch fit. + Both tile_detection options run this post-processing. Values: SEXTRACTOR_RUNNER.CCD_SIZE = 33,2080,1,4612; - SEXTRACTOR_RUNNER.MAKE_POST_PROCESS = True. + SEXTRACTOR_RUNNER.MAKE_POST_PROCESS = True; + READ_EXT_SEXCAT_RUNNER.CCD_SIZE = 33,2080,1,4612; + READ_EXT_SEXCAT_RUNNER.MAKE_POST_PROCESS = True. default: trimmed_bounds_33_2080 options: trimmed_bounds_33_2080: label: "x in (33,2080), y in (1,4612), strict" inclusive_bounds: label: Same bounds, inclusive + catalogue_neighbour_marking: + label: Neighbour pixels in the catalogue path's VIGNET + rationale: >- + ngmix masks a neighbour only where the tile VIGNET is -1e30 (flag + 2**10: zero weight, noise-filled under noisefill), and SExtractor + writes -1e30 on neighbours' footprints and off the image. The + unions_catalogue converter reproduces that from the catalogue's own + r-band segmentation map (CFIS..r.seg.fits.fz, on the tile's + pixel grid), so both tile_detection options mask the same way. + Segmentation labels are not the catalogue NUMBER; each object claims + the footprint under its centre pixel (the VIGNET centre pixel), which + on 202.301 matches 36,064 of 36,065 objects with none shared. The map + is relabelled to NUMBER (unclaimed footprints -1; an object on sky or + on a claimed footprint gets a 3-pixel-radius disc of its own), and + VIGNET pixels whose relabelled value is neither 0 nor the object's + NUMBER become -1e30. Without it, catalogue-mode ngmix fits neighbour + light: in pilot smk-g10, NGMIX_MAG_NOSHEAR was over 1 mag brighter + than MAG_AUTO for 6.1% of objects (2.7% with SExtractor, g9). + Values: + READ_EXT_SEXCAT_RUNNER.SEGMENTATION = True. + default: segmentation_map + options: + segmentation_map: + label: -1e30 on other objects' segmentation footprints, as SExtractor + none: + label: Image pixels only; ngmix masks no neighbours + excluded: true + excluded_reason: >- + Leaves neighbour light in the fit, unlike the sextractor option + it replaces; the cause of the smk-g10 magnitude outliers. prior_insights: guinot22_sextractor_params: claim: >- @@ -913,9 +969,8 @@ analyses: psf_modelling_software: label: PSF model, PSFEx per CCD or MCCD over the focal plane rationale: >- - psf_model in workflow/config.yaml's machines: table (per machine and - input type; image sims use fake) selects the exposure and tile - config pair; the committed value is psfex, which fits each CCD + psf_model in workflow/config.yaml's input_types: table (image sims + use fake) selects the exposure and tile config pair; the committed value is psfex, which fits each CCD independently. MCCD (Liaudat+2021) fits one hybrid local+global model over the focal plane; the completeness table treats its counts as warnings because no campaign has run it. Up to the model, diff --git a/scripts/python/create_final_cat.py b/scripts/python/create_final_cat.py index d68a20eda..8e989b627 100755 --- a/scripts/python/create_final_cat.py +++ b/scripts/python/create_final_cat.py @@ -33,8 +33,8 @@ def params_from_run_config(params, defaults): The workflow already knows where a campaign writes, so a manual merge should not have to restate it. Resolution goes through the workflow's own resolver (workflow/scripts/run_config.py), layering the run config on - workflow/config.yaml and then the machines: table, so what lands here is - what the rules would have used. + workflow/config.yaml and then the input_types: and machines: tables, so + what lands here is what the rules would have used. Only values still at their default are filled -- an explicit flag always wins. Nothing is derived for the data path: its patch naming differs and diff --git a/src/shapepipe/modules/fake_psf_package/fake_psf.py b/src/shapepipe/modules/fake_psf_package/fake_psf.py index 1313b0917..14e5de71f 100644 --- a/src/shapepipe/modules/fake_psf_package/fake_psf.py +++ b/src/shapepipe/modules/fake_psf_package/fake_psf.py @@ -61,7 +61,8 @@ def process(self): raise try: - n_gal = len(sex[3].data.field("NUMBER")) + numbers = np.asarray(sex[3].data.field("NUMBER")) + n_gal = len(numbers) except Exception as e: self._w_log.error( f"Error reading catalogue data from HDU 3 in {self._sexcat_path}: {e}" @@ -91,7 +92,7 @@ def process(self): output_file = SqliteDict(self._output_path) missing = 0 for idx, gal_row in enumerate(masked): - galaxy_number = idx + 1 # 1-based, matches NUMBER field + galaxy_number = int(numbers[idx]) gal_dict = {} for exp_ccd in gal_row.compressed(): if exp_ccd not in psf_dict: diff --git a/src/shapepipe/modules/make_cat_package/__init__.py b/src/shapepipe/modules/make_cat_package/__init__.py index 3525abdc8..cf1241168 100644 --- a/src/shapepipe/modules/make_cat_package/__init__.py +++ b/src/shapepipe/modules/make_cat_package/__init__.py @@ -23,7 +23,10 @@ basic measurement parameters, the PSF model at galaxy positions, and the shape measurement. Every detected object is kept: the catalogue carries no star/galaxy classification, which is done downstream. Each object is -tagged with its source tile via the ``TILE_ID`` column; objects duplicated +tagged with its source tile via the ``TILE_ID`` column and carries the +survey-wide object ID ``TILE_UNIQUE_ID = tile_id * 10**6 + NUMBER`` +(``tile_id = RRR * 1000 + DDD``, see +:func:`shapepipe.utilities.cfis.get_tile_unique_id`); objects duplicated across overlapping tiles are not deduplicated here, and are left to downstream selection. diff --git a/src/shapepipe/modules/make_cat_package/make_cat.py b/src/shapepipe/modules/make_cat_package/make_cat.py index 5c8cca748..b417f2842 100644 --- a/src/shapepipe/modules/make_cat_package/make_cat.py +++ b/src/shapepipe/modules/make_cat_package/make_cat.py @@ -21,7 +21,7 @@ get_type_flags, ) from shapepipe.pipeline import file_io -from shapepipe.utilities import mask_query +from shapepipe.utilities import cfis, mask_query def get_output_name(output_dir, file_number_string): @@ -100,7 +100,11 @@ def remove_field_name(arr, name): def save_sextractor_data(final_cat_file, sexcat_path, remove_vignet=True): """Save SExtractor Data. - Save the SExtractor catalogue into the final one. + Save the SExtractor catalogue into the final one, adding the tile as + ``TILE_ID`` (float ``RRR.DDD``) and the survey-wide object ID + ``TILE_UNIQUE_ID`` (:func:`shapepipe.utilities.cfis.get_tile_unique_id` + of the tile and ``NUMBER``). The tile is read from the catalogue's file + name, e.g. ``sexcat-301-279.fits``. Parameters ---------- @@ -123,22 +127,19 @@ def save_sextractor_data(final_cat_file, sexcat_path, remove_vignet=True): data = np.copy(sexcat_file.get_data()) if remove_vignet: data = remove_field_name(data, "VIGNET") - - final_cat_file.save_as_fits(data, ext_name="RESULTS") - cat_size = len(data) - tile_id = float( - ".".join( - re.split("-", os.path.splitext(os.path.split(sexcat_path)[1])[0])[ - 1: - ] - ) + tile_name = os.path.basename(sexcat_path) + nix, niy = cfis.get_tile_number(tile_name) + tile_id_array = np.full(cat_size, float(f"{nix}.{niy}")) + unique_id = cfis.get_tile_unique_id( + cfis.get_tile_id(tile_name), data["NUMBER"] ) - tile_id_array = np.ones(cat_size) * tile_id + final_cat_file.save_as_fits(data, ext_name="RESULTS") final_cat_file.open() final_cat_file.add_col("TILE_ID", tile_id_array) + final_cat_file.add_col("TILE_UNIQUE_ID", unique_id) sexcat_file.close() diff --git a/src/shapepipe/modules/ngmix_package/__init__.py b/src/shapepipe/modules/ngmix_package/__init__.py index 66323822c..e30d36bda 100644 --- a/src/shapepipe/modules/ngmix_package/__init__.py +++ b/src/shapepipe/modules/ngmix_package/__init__.py @@ -42,13 +42,14 @@ Save the output catalogue in batches of this size; default is ``-1`` (no batch saving) ID_OBJ_MIN : int - ID of first galaxy object to be processed; not used if set to ``-1`` - (default). Environment variables are expanded, so an orchestrator can - set the object range per chunk, for example - ``ID_OBJ_MIN = $SP_NGMIX_ID_OBJ_MIN``. + First tile-catalogue row to process, as a 1-based row position (not a + ``NUMBER`` value); not used if set to ``-1`` (default). Environment + variables are expanded, so an orchestrator can set the row range per + chunk, for example ``ID_OBJ_MIN = $NGMIX_ROW_MIN``. ID_OBJ_MAX : int - ID of last galaxy object to be processed; not used if set to ``-1`` - (default). Environment variables are expanded, as for ``ID_OBJ_MIN``. + Last tile-catalogue row to process (1-based, inclusive); not used if + set to ``-1`` (default). Environment variables are expanded, as for + ``ID_OBJ_MIN``. BKG_RMS_VIGNET_PATH : str, optional Path to a ``background_rms_vignet*.sqlite`` file produced by ``vignetmaker_runner``. The string may contain diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 707985b80..6c05ef62e 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -420,13 +420,40 @@ def get_prior(pixel_scale, rng, T_range=None, F_range=None): ) +def chunk_rows(n_obj, row_min, row_max): + """Catalogue rows of one ngmix chunk. + + A chunk is a closed range of 1-based row positions in the tile + catalogue, independent of the ``NUMBER`` values those rows carry, so a + partition of ``1..n_obj`` covers every object once however ``NUMBER`` + is ordered or spaced. A bound ``<= 0`` is unbounded on that side. + + Parameters + ---------- + n_obj : int + Number of rows in the tile catalogue + row_min, row_max : int + First and last row of the chunk (1-based, inclusive) + + Returns + ------- + range + 0-based row indices of the chunk + + """ + start = row_min - 1 if row_min > 0 else 0 + stop = min(row_max, n_obj) if row_max > 0 else n_obj + return range(start, max(start, stop)) + + def position_seed(ra, dec, ccd): """Deterministic RNG seed from an object's sky position (ngmix#796). Position seeding gives the same object the same RNG stream in each image branch, provided its sky position falls in the same seed box. It also makes the result independent of how the tile is split into - ``ID_OBJ_MIN``/``ID_OBJ_MAX`` chunks, which is why it is now the only mode. + ``ID_OBJ_MIN``/``ID_OBJ_MAX`` row chunks, which is why it is now the only + mode. Box math (kept exactly as Fabian's issue #796):: @@ -711,11 +738,11 @@ class Ngmix(object): Save output catalogue in batches of this size; detaul is ``-1`` (no batch save) id_obj_min : int, optional - First galaxy ID to process, not used if the value is set to ``-1``; - the default is ``-1`` + First catalogue row to process (1-based, see :func:`chunk_rows`), + not used if the value is set to ``-1``; the default is ``-1`` id_obj_max : int, optional - Last galaxy ID to process, not used if the value is set to ``-1``; - the default is ``-1`` + Last catalogue row to process (1-based, inclusive), not used if the + value is set to ``-1``; the default is ``-1`` centroid_source : {"wcs", "hsm"}, optional How to place the galaxy Jacobian origin for the centroid prior. The default ``"wcs"`` places it at the coadd centroid: the sub-pixel @@ -1256,11 +1283,11 @@ def process(self): count_batch = 0 saved_batch_cumul = 0 - for i_tile, obj_id in enumerate(tile_cat.obj_id): - if self._id_obj_min > 0 and obj_id < self._id_obj_min: - continue - if self._id_obj_max > 0 and obj_id > self._id_obj_max: - continue + rows = chunk_rows( + len(tile_cat.obj_id), self._id_obj_min, self._id_obj_max + ) + for i_tile in rows: + obj_id = tile_cat.obj_id[i_tile] if id_first == -1: id_first = obj_id id_last = obj_id diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 1509803ef..75d26f60a 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -122,10 +122,10 @@ def ngmix_runner( # No batch saving save_batch = -1 - # First and last galaxy ID to process. Read via ``getexpanded`` so an - # orchestrator can drive the chunk bounds from environment variables - # (``$SP_NGMIX_ID_OBJ_MIN`` and friends); ``getexpanded`` is the only - # accessor in ShapePipe's config that expands ``$VAR``. + # First and last catalogue row (1-based) to process. Read via + # ``getexpanded`` so an orchestrator can drive the chunk bounds from + # environment variables (``$NGMIX_ROW_MIN`` and friends); ``getexpanded`` + # is the only accessor in ShapePipe's config that expands ``$VAR``. id_obj_min = int(config.getexpanded(module_config_sec, "ID_OBJ_MIN")) id_obj_max = int(config.getexpanded(module_config_sec, "ID_OBJ_MAX")) diff --git a/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py b/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py index f1224d791..ff4f58cc8 100644 --- a/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py +++ b/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py @@ -47,12 +47,114 @@ def _build_ldac_imhead(img_header): return hdu -def _extract_vignets(image_data, x_pos, y_pos, stamp_size): +# SExtractor's VIGNET value for pixels that are not the object's: off the +# image and on a neighbour's segmentation footprint. ngmix flags these pixels +# (tile VIGNET == -1e30) and noise-fills them at zero weight. +BIG = -1e30 + +# The label of a footprint no catalogue object claims in the relabelled +# segmentation map. Negative, so it never collides with a NUMBER; VIGNET +# marking only asks "self or not self", so neighbours need no identity. +NEIGHBOUR_LABEL = -1 + +# Radius, in pixels, of the disc of its own NUMBER painted for an object with +# no footprint of its own (on sky, or on a footprint claimed by another). +FALLBACK_RADIUS = 3 + + +def _centre_pixels(x_pos, y_pos): + """0-based (column, row) of the pixel holding each 1-based position.""" + col = np.rint(np.asarray(x_pos, dtype=float)).astype(np.int64) - 1 + row = np.rint(np.asarray(y_pos, dtype=float)).astype(np.int64) - 1 + return col, row + + +def relabel_seg(seg, number, x_image, y_image, fallback_radius=FALLBACK_RADIUS): + """Relabel a segmentation map into a catalogue's NUMBER. + + Each object claims, in NUMBER order and first come first served, the + footprint its centre pixel falls in; that footprint becomes its NUMBER. + Footprints nobody claims become ``NEIGHBOUR_LABEL``; sky stays 0. An + object on sky or on an already-claimed footprint gets a disc of radius + ``fallback_radius`` of its NUMBER, and every object on the image ends + holding its own centre pixel (the lower NUMBER wins a shared pixel). + + Parameters + ---------- + seg : numpy.ndarray + Segmentation map, 0 for sky and one positive label per footprint + number : array_like + Catalogue ``NUMBER`` + x_image, y_image : array_like + 1-based pixel positions on the grid of ``seg`` + fallback_radius : int, optional + Radius of the disc painted for an object without a footprint + + Returns + ------- + numpy.ndarray + Relabelled map, ``int32`` + dict + Counts of the claim outcomes ``matched``, ``unclaimed`` (on sky), + ``shared`` and ``off_image``, which partition the catalogue, plus + ``shared_pixel`` (objects rounding to a lower NUMBER's pixel) + + @sc [decision:detection.catalogue_neighbour_marking] + """ + number = np.asarray(number) + col, row = _centre_pixels(x_image, y_image) + n_row, n_col = seg.shape + inside = (col >= 0) & (col < n_col) & (row >= 0) & (row < n_row) + order = np.argsort(number, kind="stable") + + counts = dict(matched=0, unclaimed=0, shared=0, off_image=0, + shared_pixel=0) + table = np.full(max(int(seg.max()), 0) + 1, NEIGHBOUR_LABEL, np.int32) + table[0] = 0 + claimed, fallback = set(), [] + for i in order: + if not inside[i]: + counts["off_image"] += 1 + continue + label = int(seg[row[i], col[i]]) + if label == 0 or label in claimed: + counts["unclaimed" if label == 0 else "shared"] += 1 + fallback.append(i) + else: + claimed.add(label) + table[label] = number[i] + counts["matched"] += 1 + out = table[np.clip(seg, 0, None)] + + r = np.arange(-fallback_radius, fallback_radius + 1) + disc = np.argwhere(np.add.outer(r**2, r**2) <= fallback_radius**2) + disc -= fallback_radius + for i in fallback: + out[np.clip(row[i] + disc[:, 0], 0, n_row - 1), + np.clip(col[i] + disc[:, 1], 0, n_col - 1)] = number[i] + + # Centres last, so no disc can take an object's centre away. + taken = set() + for i in order[inside[order]]: + pixel = (row[i], col[i]) + if pixel in taken: + counts["shared_pixel"] += 1 + continue + taken.add(pixel) + out[pixel] = number[i] + return out, counts + + +def _extract_vignets(image_data, x_pos, y_pos, stamp_size, seg=None, + number=None): """Extract postage stamps from a tile image array. For each object position, a ``stamp_size x stamp_size`` cutout is - extracted. Objects whose stamp falls partially outside the image are - padded with zeros, matching SExtractor behaviour. + extracted, centred on the pixel holding the position. Pixels off the + image are set to ``BIG``, as SExtractor does. With a segmentation map in + the catalogue's numbering (:func:`relabel_seg`), pixels on any footprint + other than the object's own are set to ``BIG`` too, which is how + SExtractor's VIGNET marks neighbours. Parameters ---------- @@ -64,45 +166,47 @@ def _extract_vignets(image_data, x_pos, y_pos, stamp_size): Y pixel positions, 1-based (SExtractor convention) stamp_size : int Side length of the square postage stamp (should be odd) + seg : numpy.ndarray, optional + Segmentation map on the image grid, labelled with ``number`` + number : array_like, optional + Catalogue ``NUMBER``, required with ``seg`` Returns ------- numpy.ndarray Array of shape ``(n_obj, stamp_size, stamp_size)``, dtype float32 + @sc [decision:detection.catalogue_neighbour_marking] """ ny, nx = image_data.shape half = stamp_size // 2 - n_obj = len(x_pos) - vignets = np.zeros((n_obj, stamp_size, stamp_size), dtype=np.float32) - - for i, (x, y) in enumerate(zip(x_pos, y_pos)): - xi = int(round(float(x))) - 1 - yi = int(round(float(y))) - 1 - - x0, x1 = xi - half, xi + half + 1 - y0, y1 = yi - half, yi + half + 1 - - xc0, xc1 = max(0, x0), min(nx, x1) - yc0, yc1 = max(0, y0), min(ny, y1) - - dx0, dx1 = xc0 - x0, xc0 - x0 + (xc1 - xc0) - dy0, dy1 = yc0 - y0, yc0 - y0 + (yc1 - yc0) - - vignets[i, dy0:dy1, dx0:dx1] = image_data[yc0:yc1, xc0:xc1] + col, row = _centre_pixels(x_pos, y_pos) + vignets = np.full((len(col), stamp_size, stamp_size), BIG, np.float32) + seg_stamp = np.zeros((stamp_size, stamp_size), np.int32) + + for i, (xi, yi) in enumerate(zip(col, row)): + x0, y0 = xi - half, yi - half + xc0, xc1 = max(0, x0), min(nx, x0 + stamp_size) + yc0, yc1 = max(0, y0), min(ny, y0 + stamp_size) + if xc0 >= xc1 or yc0 >= yc1: + continue + inner = np.s_[yc0 - y0:yc1 - y0, xc0 - x0:xc1 - x0] + vignets[i][inner] = image_data[yc0:yc1, xc0:xc1] + if seg is not None: + seg_stamp[:] = 0 + seg_stamp[inner] = seg[yc0:yc1, xc0:xc1] + vignets[i][(seg_stamp != 0) & (seg_stamp != number[i])] = BIG return vignets -def _build_ldac_objects(cat_data, unique_id, vignets): +def _build_ldac_objects(cat_data, vignets): """Build LDAC_OBJECTS extension from an astropy table. Parameters ---------- cat_data : astropy.table.Table Catalogue data read from the ASCII SExtractor file - unique_id : numpy.ndarray - 1-D int64 array of per-object unique IDs across all tiles vignets : numpy.ndarray Array of shape ``(n_obj, stamp_size, stamp_size)`` @@ -127,10 +231,6 @@ def _build_ldac_objects(cat_data, unique_id, vignets): fmt = f"{arr.dtype.itemsize}A" fits_cols.append(fits.Column(name=colname, format=fmt, array=arr)) - fits_cols.append( - fits.Column(name="TILE_UNIQUE_ID", format="K", array=unique_id) - ) - _aliases = { "XWIN_IMAGE": "X_IMAGE", "YWIN_IMAGE": "Y_IMAGE", @@ -162,47 +262,12 @@ def _build_ldac_objects(cat_data, unique_id, vignets): return hdu -def _tile_id_from_file_number_string(file_number_string): - """Decode file-number string into a collision-free tile ID. - - CFIS tile file-number strings are of the form ``'-RRR-DDD'`` where - ``RRR`` is the RA grid index and ``DDD`` the dec grid index (both - 3-digit). The encoded tile_id is ``RRR * 1000 + DDD``; the 3-digit - assumption is an invariant of the CFIS grid, not a data assumption, - so we assert it rather than silently overflowing. - - Parameters - ---------- - file_number_string : str - Pipeline file-number string, e.g. ``'-301-279'`` - - Returns - ------- - int - Collision-free tile_id (e.g. ``301279``) - - Raises - ------ - ValueError - If either grid component exceeds three digits. - - """ - parts = file_number_string.lstrip("-").split("-") - ra_idx, dec_idx = int(parts[0]), int(parts[1]) - if not (0 <= ra_idx < 1000 and 0 <= dec_idx < 1000): - raise ValueError( - f"file_number_string {file_number_string!r}: both grid " - f"components must be in [0, 1000); got ra={ra_idx}, dec={dec_idx}" - ) - return ra_idx * 1000 + dec_idx - - def make_ldac_from_ascii( input_cat_path, image_path, output_cat_path, - file_number_string, stamp_size=51, + seg_path=None, w_log=None, ): """Convert an external ASCII catalogue to FITS-LDAC format. @@ -214,12 +279,12 @@ def make_ldac_from_ascii( writes a standard FITS-LDAC file compatible with all downstream ShapePipe modules. - A ``TILE_UNIQUE_ID`` column and a ``VIGNET`` column (postage stamps - extracted from the tile image) are added to ``LDAC_OBJECTS``. - - The unique ID is computed as ``tile_id * 10**6 + NUMBER`` where - ``tile_id`` is derived from the tile RA/Dec grid coordinates encoded in - ``file_number_string`` (e.g. ``'-301-279'`` → ``tile_id = 301279``). + The input columns, ``NUMBER`` included, are copied unchanged; a + ``VIGNET`` column (postage stamps extracted from the tile image) is + added to ``LDAC_OBJECTS``. Given the catalogue's segmentation map, which + shares the tile's pixel grid, the map is relabelled to the catalogue's + ``NUMBER`` (:func:`relabel_seg`) and used to set neighbours' pixels in each ``VIGNET`` to ``BIG``, as + SExtractor does. Parameters ---------- @@ -229,10 +294,10 @@ def make_ldac_from_ascii( Path to tile image FITS file output_cat_path : str Path to the output FITS-LDAC catalogue - file_number_string : str - Pipeline file-number string, e.g. ``'-301-279'`` stamp_size : int, optional Side length of the square postage stamp in pixels, default 51 + seg_path : str, optional + Path to the catalogue's segmentation map (FITS, compressed or not) w_log : logging.Logger, optional Pipeline logger @@ -242,15 +307,31 @@ def make_ldac_from_ascii( if w_log: w_log.info(f"Read {n_obj} objects from {input_cat_path}") - tile_id = _tile_id_from_file_number_string(file_number_string) - unique_id = ( - tile_id * 10**6 + np.array(cat_data["NUMBER"], dtype=np.int64) - ) - with fits.open(image_path) as hdul: img_header = hdul[0].header image_data = hdul[0].data.astype(np.float32) + seg = None + if seg_path is not None: + with fits.open(seg_path) as hdul: + hdu = next(h for h in hdul if h.data is not None) + seg_raw = hdu.data + if seg_raw.shape != image_data.shape: + raise ValueError( + f"Segmentation map {seg_path} has shape {seg_raw.shape}, the" + + f" image {image_data.shape}; they must share one grid." + ) + seg, counts = relabel_seg( + seg_raw, cat_data["NUMBER"], cat_data["X_IMAGE"], + cat_data["Y_IMAGE"], + ) + del seg_raw + if w_log: + w_log.info( + f"Relabelled {seg_path} to NUMBER: " + + ", ".join(f"{k}={v}" for k, v in counts.items()) + ) + if w_log: w_log.info( f"Extracting {stamp_size}x{stamp_size} vignets from {image_path}" @@ -260,10 +341,12 @@ def make_ldac_from_ascii( cat_data["X_IMAGE"], cat_data["Y_IMAGE"], stamp_size, + seg=seg, + number=np.asarray(cat_data["NUMBER"]), ) ldac_imhead = _build_ldac_imhead(img_header) - ldac_objects = _build_ldac_objects(cat_data, unique_id, vignets) + ldac_objects = _build_ldac_objects(cat_data, vignets) hdul_out = fits.HDUList([fits.PrimaryHDU(), ldac_imhead, ldac_objects]) hdul_out.writeto(output_cat_path, overwrite=True) diff --git a/src/shapepipe/modules/read_ext_sexcat_runner.py b/src/shapepipe/modules/read_ext_sexcat_runner.py index 0e79c16bf..17e23163d 100644 --- a/src/shapepipe/modules/read_ext_sexcat_runner.py +++ b/src/shapepipe/modules/read_ext_sexcat_runner.py @@ -30,11 +30,18 @@ def read_ext_sexcat_runner( Reads an external ASCII catalogue (SExtractor format), converts it to a FITS-LDAC catalogue compatible with downstream ShapePipe modules. - If MAKE_POST_PROCESS = True, runs multi-epoch post-processing to add - per-exposure HDUs. + The inputs are the catalogue and the tile image, then the catalogue's + segmentation map if SEGMENTATION = True, then the WCS log if + MAKE_POST_PROCESS = True. With the segmentation map, neighbours' pixels + in VIGNET are set to -1e30; the map itself is not written out; the + only output is the FITS-LDAC ``sexcat.fits``. MAKE_POST_PROCESS runs the + multi-epoch post-processing that adds per-exposure HDUs. """ - cat_path = input_file_list[0] - image_path = input_file_list[1] + cat_path, image_path, *extra_inputs = input_file_list + use_seg = config.has_option( + module_config_sec, "SEGMENTATION" + ) and config.getboolean(module_config_sec, "SEGMENTATION") + seg_path = extra_inputs.pop(0) if use_seg else None if config.has_option(module_config_sec, "SUFFIX"): suffix = config.get(module_config_sec, "SUFFIX") @@ -55,21 +62,21 @@ def read_ext_sexcat_runner( cat_path, image_path, output_path, - file_number_string, stamp_size=stamp_size, + seg_path=seg_path, w_log=w_log, ) if config.getboolean(module_config_sec, "MAKE_POST_PROCESS"): # The WCS log is supplied as a positional input (via FILE_PATTERN) # when post-processing is enabled, not by the decorator default. - if len(input_file_list) < 3: + if not extra_inputs: raise ValueError( - "MAKE_POST_PROCESS requires the WCS log file as a third" + "MAKE_POST_PROCESS requires the WCS log file as the last" + " input; add 'log_exp_headers' to FILE_PATTERN and" + f" FILE_EXT in the [{module_config_sec}] config section." ) - f_wcs_path = input_file_list[2] + f_wcs_path = extra_inputs[0] pos_params = config.getlist(module_config_sec, "WORLD_POSITION") ccd_size = config.getlist(module_config_sec, "CCD_SIZE") w_log.info("Running post-processing") diff --git a/src/shapepipe/utilities/cfis.py b/src/shapepipe/utilities/cfis.py index 0d94b3be3..205d44a6f 100644 --- a/src/shapepipe/utilities/cfis.py +++ b/src/shapepipe/utilities/cfis.py @@ -563,8 +563,8 @@ def get_tile_number(tile_name): tile number for x and tile number for y """ - m = re.search(r"(\d{3})[\.-](\d{3})", tile_name) - if m is None or len(m.groups()) != 2: + m = re.search(r"(?= TILE_UNIQUE_ID_BASE + ): + raise CfisError( + f"Object number outside [0, {TILE_UNIQUE_ID_BASE}) in tile " + f"{tile_id}: range [{number.min()}, {number.max()}]" + ) + unique_id = np.int64(tile_id) * TILE_UNIQUE_ID_BASE + number + return unique_id[()] if unique_id.ndim == 0 else unique_id + + +def split_tile_unique_id(unique_id): + """Split Tile Unique ID. + + Inverse of :func:`get_tile_unique_id`. + + Parameters + ---------- + unique_id : int or array_like + Unique ID(s) + + Returns + ------- + tuple + tile ID(s) and object number(s) + + """ + unique_id = np.asarray(unique_id, dtype=np.int64) + return np.divmod(unique_id, TILE_UNIQUE_ID_BASE) + + def get_log_file(path, verbose=False): """Get Log File. diff --git a/tests/module/data/dr6_202.301_seg_patch.fits b/tests/module/data/dr6_202.301_seg_patch.fits new file mode 100644 index 000000000..c465e6446 Binary files /dev/null and b/tests/module/data/dr6_202.301_seg_patch.fits differ diff --git a/tests/module/test_final_cat_columns.py b/tests/module/test_final_cat_columns.py new file mode 100644 index 000000000..51ea46772 --- /dev/null +++ b/tests/module/test_final_cat_columns.py @@ -0,0 +1,141 @@ +"""The final-catalogue columns each tile_detection mode requests. + +``final_cat_merge`` asks every tile catalogue for the columns of +``workflow/config/cfis/final_cat.param`` and stops on any one missing. That +list is written against SExtractor-mode tiles; under ``tile_detection: +unions_catalogue`` the detection columns are whatever the UNIONS per-tile +catalogue carries, copied through ``read_ext_sexcat``. So the catalogue-mode +request is the param list less ``merge_final_cat.SEXTRACTOR_ONLY_COLUMNS``. + +Here the catalogue-mode detection columns are not written down but derived: +the converter runs on a catalogue with the UNIONS DR6 header (copied verbatim +from ``vos:cfis/tiles_DR6/CFIS.202.301.r.cat``) and its output is what a tile +can carry. The requested SExtractor columns are those of the tile SExtractor +parameter file. Every other requested column comes from stages that run the +same way in both modes (post-processing, ngmix, make_cat), so the detection +columns are the whole of the difference between them. +""" + +import importlib.util +import re +import sys +from pathlib import Path + +import numpy as np +import pytest +from astropy.io import fits + +from shapepipe.modules.read_ext_sexcat_package import read_ext_sexcat as rs + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPTS = REPO_ROOT / "workflow" / "scripts" +CFIS = REPO_ROOT / "workflow" / "config" / "cfis" + +# The header of a UNIONS DR6 per-tile catalogue and its first object, moved to +# pixel (4, 4) so that its stamp falls on the small test image. +DR6_CATALOGUE = """\ +# 1 NUMBER Running object number +# 2 X_IMAGE Object position along x [pixel] +# 3 Y_IMAGE Object position along y [pixel] +# 4 ALPHA_J2000 Right ascension of barycenter (J2000) [deg] +# 5 DELTA_J2000 Declination of barycenter (J2000) [deg] +# 6 MAG_AUTO Kron-like elliptical aperture magnitude [mag] +# 7 MAGERR_AUTO RMS error for AUTO magnitude [mag] +# 8 MAG_BEST Best of MAG_AUTO and MAG_ISOCOR [mag] +# 9 MAGERR_BEST RMS error for MAG_BEST [mag] +# 10 MAG_APER Fixed aperture magnitude vector [mag] +# 11 MAGERR_APER RMS error vector for fixed aperture mag. [mag] +# 12 A_WORLD Profile RMS along major axis (world units) [deg] +# 13 ERRA_WORLD World RMS position error along major axis [deg] +# 14 B_WORLD Profile RMS along minor axis (world units) [deg] +# 15 ERRB_WORLD World RMS position error along minor axis [deg] +# 16 THETA_J2000 Position angle (east of north) (J2000) [deg] +# 17 ERRTHETA_J2000 J2000 error ellipse pos. angle (east of north) [deg] +# 18 ISOAREA_IMAGE Isophotal area above Analysis threshold [pixel**2] +# 19 MU_MAX Peak surface brightness above background [mag * arcsec**(-2)] +# 20 FLUX_RADIUS Fraction-of-light radii [pixel] +# 21 FLAGS Extraction flags + 1 4.0000 4.0000 205.3679556 +60.2576198 14.8792 0.0002 14.8792 0.0002 14.9100 0.0002 0.000389709 2.80088e-07 0.0003095117 1.974931e-07 -1.48 -0.72 5707 16.3802 4.474 0 +""" + + +def _load_script(name, path): + sys.path.insert(0, str(SCRIPTS)) + try: + spec = importlib.util.spec_from_file_location(f"_{name}", path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + finally: + sys.path.remove(str(SCRIPTS)) + return module + + +merge_final_cat = _load_script("merge_final_cat", + SCRIPTS / "merge_final_cat.py") +create_final_cat = _load_script( + "create_final_cat", + REPO_ROOT / "scripts" / "python" / "create_final_cat.py") + + +@pytest.fixture(scope="module") +def param_list(): + return create_final_cat.read_param_file(str(CFIS / "final_cat.param")) + + +@pytest.fixture(scope="module") +def sextractor_columns(): + """Column names the tile SExtractor run writes (vector sizes stripped).""" + names = set() + for line in (CFIS / "default_noimaflags.param").read_text().splitlines(): + entry = line.split("#")[0].strip() + if entry: + names.add(re.sub(r"\(.*\)$", "", entry)) + return names + + +@pytest.fixture(scope="module") +def catalogue_columns(tmp_path_factory): + """Detection columns of a tile under tile_detection: unions_catalogue.""" + tmp = tmp_path_factory.mktemp("dr6") + cat, img, out = (tmp / "CFIS.202.301.r.cat", tmp / "CFIS.202.301.r.fits", + tmp / "sexcat-202-301.fits") + cat.write_text(DR6_CATALOGUE) + fits.PrimaryHDU(np.zeros((8, 8), dtype=np.float32)).writeto(img) + rs.make_ldac_from_ascii(str(cat), str(img), str(out), stamp_size=3) + with fits.open(out) as hdul: + return set(hdul["LDAC_OBJECTS"].columns.names) - {"VIGNET"} + + +def test_sextractor_only_is_exactly_what_the_catalogue_lacks( + param_list, sextractor_columns, catalogue_columns): + """Every requested SExtractor column the catalogue lacks, and no other.""" + lacking = {c for c in param_list + if c in sextractor_columns and c not in catalogue_columns} + assert set(merge_final_cat.SEXTRACTOR_ONLY_COLUMNS) == lacking + + +def test_catalogue_mode_request_is_carried( + param_list, sextractor_columns, catalogue_columns): + """Each requested detection column is one the converter writes.""" + requested = merge_final_cat.requested_columns(param_list, + "unions_catalogue") + detection = [c for c in requested if c in sextractor_columns] + assert detection and set(detection) <= catalogue_columns + # Only detection columns are dropped, and the order is the param file's. + assert requested == [c for c in param_list + if c not in merge_final_cat.SEXTRACTOR_ONLY_COLUMNS] + assert {"MAG_AUTO", "MAGERR_AUTO", "FLUX_RADIUS"} <= set(requested) + + +def test_sextractor_mode_requests_the_param_file(): + for input_type in ("cfis", "cfis_image_sims"): + params = create_final_cat.read_param_file( + str(REPO_ROOT / "workflow" / "config" / input_type + / "final_cat.param")) + assert merge_final_cat.requested_columns( + params, "sextractor") == params + + +def test_unknown_mode_is_refused(param_list): + with pytest.raises(ValueError): + merge_final_cat.requested_columns(param_list, "sextractr") diff --git a/tests/module/test_make_cat.py b/tests/module/test_make_cat.py index bdcd13f84..d1a7f6a11 100644 --- a/tests/module/test_make_cat.py +++ b/tests/module/test_make_cat.py @@ -30,6 +30,7 @@ from shapepipe.modules.ngmix_package.ngmix import Ngmix from shapepipe.pipeline import file_io from shapepipe.pipeline.config import CustomParser +from shapepipe.utilities import cfis class _NullLogger: @@ -832,3 +833,62 @@ def test_save_psf_data_carries_fourth_moments_per_epoch(tmp_path): npt.assert_allclose(out[f"HSM_M4_2_PSF_{n}"], [-10.0]) npt.assert_allclose(out[f"HSM_RHO4_PSF_{n}"], [-1.0]) npt.assert_allclose(out["HSM_G1_PSF_2"], [0.03]) + + +def _write_sexcat(path, number): + """Write a minimal FITS-LDAC tile catalogue (``LDAC_OBJECTS`` at HDU 2).""" + n_obj = len(number) + cols = [ + fits.Column(name="NUMBER", format="J", array=np.asarray(number)), + fits.Column(name="XWIN_WORLD", format="D", array=np.zeros(n_obj)), + fits.Column( + name="VIGNET", format="4E", dim="(2,2)", + array=np.zeros((n_obj, 2, 2), dtype=np.float32), + ), + ] + fits.HDUList([ + fits.PrimaryHDU(), + fits.BinTableHDU.from_columns( + [fits.Column(name="Field Header Card", format="80A", + array=np.array([""]))], + name="LDAC_IMHEAD", + ), + fits.BinTableHDU.from_columns(cols, name="LDAC_OBJECTS"), + ]).writeto(path, overwrite=True) + + +def test_save_sextractor_data_writes_tile_unique_id(tmp_path): + """SExtractor-mode final catalogue carries tile_id * 10**6 + NUMBER. + + ``NUMBER`` is deliberately gapped and unsorted: the ID is built from the + value each row carries, not from its position. + """ + number = np.array([7, 3, 12, 999999]) + sexcat = tmp_path / "sexcat-301-279.fits" + _write_sexcat(sexcat, number) + + final_cat = make_cat.prepare_final_cat_file(str(tmp_path), "-301-279") + n_obj = make_cat.save_sextractor_data(final_cat, str(sexcat)) + final_cat.close() + + assert n_obj == len(number) + with fits.open(tmp_path / "final_cat-301-279.fits") as hdul: + data = hdul["RESULTS"].data + assert "VIGNET" not in data.names + npt.assert_array_equal(data["NUMBER"], number) + assert data["TILE_UNIQUE_ID"].dtype.newbyteorder("=") == np.int64 + npt.assert_array_equal( + data["TILE_UNIQUE_ID"], 301279 * 10**6 + number + ) + npt.assert_allclose(data["TILE_ID"], 301.279) + + +def test_save_sextractor_data_refuses_number_beyond_id_range(tmp_path): + """A NUMBER that would overflow into the tile digits raises, writes nothing.""" + sexcat = tmp_path / "sexcat-301-279.fits" + _write_sexcat(sexcat, np.array([1, 10**6])) + + final_cat = make_cat.prepare_final_cat_file(str(tmp_path), "-301-279") + with pytest.raises(cfis.CfisError): + make_cat.save_sextractor_data(final_cat, str(sexcat)) + assert not (tmp_path / "final_cat-301-279.fits").exists() diff --git a/tests/module/test_ngmix_chunk_rows.py b/tests/module/test_ngmix_chunk_rows.py new file mode 100644 index 000000000..e2f0c8581 --- /dev/null +++ b/tests/module/test_ngmix_chunk_rows.py @@ -0,0 +1,110 @@ +"""ngmix chunks select catalogue rows, so they cover every object once. + +The workflow partitions a tile into chunks with ``ngmix_range.row_ranges`` +(1-based closed row ranges) and each ngmix chunk keeps the rows +``ngmix.chunk_rows`` returns. The pair must visit every row exactly once +whatever the ``NUMBER`` column holds: an external detection catalogue keeps +its own, possibly gapped and unsorted, ``NUMBER``. +""" + +import importlib.util +from pathlib import Path + +import numpy as np +import pytest +from astropy.io import fits +from hypothesis import given, settings +from hypothesis import strategies as st + +from shapepipe.modules.ngmix_package.ngmix import chunk_rows + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPT = REPO_ROOT / "workflow" / "scripts" / "ngmix_range.py" + + +def _load_ngmix_range(): + spec = importlib.util.spec_from_file_location("_ngmix_range", SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +ngmix_range = _load_ngmix_range() + + +@settings(max_examples=200, deadline=None) +@given( + st.integers(min_value=1, max_value=300), + st.integers(min_value=1, max_value=12), + st.data(), +) +def test_chunks_cover_every_row_once(n_obj, n_chunks, data): + epochs = data.draw( + st.lists(st.integers(0, 20), min_size=n_obj, max_size=n_obj) + ) + # Gapped, unsorted NUMBER: distinct values drawn from a wide range. + number = np.array( + data.draw( + st.lists( + st.integers(1, 10**6 - 1), + min_size=n_obj, + max_size=n_obj, + unique=True, + ) + ) + ) + + selected = [ + number[i] + for lo, hi in ngmix_range.row_ranges(epochs, n_chunks) + for i in chunk_rows(n_obj, lo, hi) + ] + assert selected == list(number) + + +def test_unbounded_and_empty_chunks(): + assert chunk_rows(5, -1, -1) == range(0, 5) + assert chunk_rows(5, 3, -1) == range(2, 5) + assert chunk_rows(5, -1, 2) == range(0, 2) + # The splitter's canonical empty range (n_obj + 1, n_obj). + assert len(chunk_rows(5, 6, 5)) == 0 + + +def _write_sexcat(run_dir, number, epoch_numbers, ccd_n): + out = run_dir / "output/run_sp_tile_Sx/sextractor_runner/output" + out.mkdir(parents=True) + hdus = [ + fits.PrimaryHDU(), + fits.BinTableHDU.from_columns( + [fits.Column(name="NUMBER", format="K", array=number)], + name="LDAC_OBJECTS", + ), + ] + for k, (enum, ccd) in enumerate(zip(epoch_numbers, ccd_n)): + hdus.append( + fits.BinTableHDU.from_columns( + [ + fits.Column(name="NUMBER", format="K", array=enum), + fits.Column(name="CCD_N", format="K", array=ccd), + ], + name=f"EPOCH_{k}", + ) + ) + fits.HDUList(hdus).writeto(out / "sexcat-301-279.fits") + + +def test_object_epochs_accepts_gapped_unsorted_number(tmp_path): + number = np.array([40, 7, 1000, 3]) + ccd_n = [np.array([0, -1, 5, 2]), np.array([1, 1, -1, -1])] + _write_sexcat(tmp_path, number, [number, number], ccd_n) + np.testing.assert_array_equal( + ngmix_range.object_epochs(tmp_path), [2, 1, 1, 1] + ) + + +def test_object_epochs_refuses_misaligned_epoch_rows(tmp_path): + number = np.array([40, 7, 1000, 3]) + ccd_n = [np.zeros(4, dtype=int), np.zeros(4, dtype=int)] + _write_sexcat(tmp_path, number, [number, number[::-1]], ccd_n) + with pytest.raises(SystemExit, match="aligned"): + ngmix_range.object_epochs(tmp_path) diff --git a/tests/module/test_read_ext_sexcat.py b/tests/module/test_read_ext_sexcat.py new file mode 100644 index 000000000..fa406b9a6 --- /dev/null +++ b/tests/module/test_read_ext_sexcat.py @@ -0,0 +1,424 @@ +"""UNIT TESTS FOR MODULE PACKAGE: READ_EXT_SEXCAT. + +Drives ``make_ldac_from_ascii`` on a synthetic ASCII SExtractor-format +catalogue and a synthetic tile image, and checks the FITS-LDAC it writes is +what the tile chain downstream of ``tile_detect`` reads: the LDAC_IMHEAD +extension carrying the tile header, the SExtractor column aliases, one +``VIGNET`` stamp per object cut from the image, and the input ``NUMBER`` +kept as is. It follows the catalogue through +``make_cat.save_sextractor_data``, which builds ``TILE_UNIQUE_ID``. The rest +covers the segmentation map: relabelling it to the catalogue's ``NUMBER`` and +setting neighbours' ``VIGNET`` pixels to -1e30, as SExtractor does, which is +all ngmix reads to mask neighbours. +""" + +from pathlib import Path + +import numpy as np +import numpy.testing as npt +import pytest +from astropy.io import fits + +from shapepipe.modules.make_cat_package import make_cat +from shapepipe.modules.read_ext_sexcat_package import read_ext_sexcat as rs + +NX, NY = 40, 30 +STAMP = 5 +# (NUMBER, X_IMAGE, Y_IMAGE): an interior object, one on the left edge, one +# in the top-right corner. +OBJECTS = [(1, 10.0, 12.0), (2, 1.0, 20.0), (7, 40.0, 30.0)] + + +def _write_ascii_cat(path): + lines = [ + "# 1 NUMBER Running object number", + "# 2 X_IMAGE Object position along x [pixel]", + "# 3 Y_IMAGE Object position along y [pixel]", + "# 4 ALPHA_J2000 Right ascension of barycenter [deg]", + "# 5 DELTA_J2000 Declination of barycenter [deg]", + "# 6 MAG_AUTO Kron-like elliptical aperture magnitude [mag]", + ] + for num, x, y in OBJECTS: + lines.append(f"{num} {x} {y} {150.0 + num} {30.0 + num} {20.0 + num}") + path.write_text("\n".join(lines) + "\n") + + +def _write_image(path): + # pixel value = 1000*row + column (0-based), so every stamp pixel names + # where it came from. + data = (np.arange(NY)[:, None] * 1000 + np.arange(NX)[None, :]).astype( + np.float32 + ) + hdu = fits.PrimaryHDU(data) + hdu.header["HISTORY"] = "input image 2605805p.fits" + hdu.header["TILEKEY"] = "kept" + hdu.writeto(path, overwrite=True) + + +@pytest.fixture +def ldac(tmp_path): + cat = tmp_path / "CFIS_cat-301-279.cat" + img = tmp_path / "CFIS_image-301-279.fits" + out = tmp_path / "sexcat-301-279.fits" + _write_ascii_cat(cat) + _write_image(img) + rs.make_ldac_from_ascii( + str(cat), str(img), str(out), stamp_size=STAMP + ) + return out + + +def test_ldac_layout_and_header(ldac): + with fits.open(ldac) as hdul: + assert [h.name for h in hdul] == ["PRIMARY", "LDAC_IMHEAD", "LDAC_OBJECTS"] + cards = hdul["LDAC_IMHEAD"].data[0][0] + assert isinstance(cards, str) or cards.ndim == 1 + text = "".join(cards) if not isinstance(cards, str) else cards + assert "TILEKEY" in text and "2605805p" in text + + +def test_number_is_kept_and_aliases_added(ldac): + with fits.open(ldac) as hdul: + data = hdul["LDAC_OBJECTS"].data + # The input NUMBER (gapped here: 1, 2, 7) is the object's identity and is + # copied unchanged; the converter builds no ID of its own. + npt.assert_array_equal(data["NUMBER"], [o[0] for o in OBJECTS]) + assert "TILE_UNIQUE_ID" not in data.names + npt.assert_array_equal(data["XWIN_IMAGE"], data["X_IMAGE"]) + npt.assert_array_equal(data["YWIN_IMAGE"], data["Y_IMAGE"]) + npt.assert_array_equal(data["XWIN_WORLD"], data["ALPHA_J2000"]) + npt.assert_array_equal(data["YWIN_WORLD"], data["DELTA_J2000"]) + + +def test_vignets_are_cut_from_the_image_and_padded_as_sextractor(ldac): + with fits.open(ldac) as hdul: + vignets = hdul["LDAC_OBJECTS"].data["VIGNET"] + assert vignets.shape == (len(OBJECTS), STAMP, STAMP) + + # Interior object at (10, 12), 1-based: centre pixel is (row 11, col 9). + assert vignets[0, STAMP // 2, STAMP // 2] == 11 * 1000 + 9 + assert vignets[0, 0, 0] == 9 * 1000 + 7 + + # Left edge, x = 1: the two columns left of the image are -1e30, the value + # SExtractor writes off the image. + assert (vignets[1, :, :2] == rs.BIG).all() + assert vignets[1, STAMP // 2, STAMP // 2] == 19 * 1000 + 0 + + # Top-right corner: only the lower-left quadrant of the stamp is in the + # image. + assert (vignets[2, STAMP // 2 + 1:, :] == rs.BIG).all() + assert (vignets[2, :, STAMP // 2 + 1:] == rs.BIG).all() + assert vignets[2, STAMP // 2, STAMP // 2] == 29 * 1000 + 39 + + +def test_tile_unique_id_reaches_the_final_catalogue(ldac, tmp_path): + """make_cat builds the ID from the tile and the catalogue's own NUMBER.""" + final = make_cat.prepare_final_cat_file(str(tmp_path), "-301-279") + n_obj = make_cat.save_sextractor_data(final, str(ldac)) + assert n_obj == len(OBJECTS) + with fits.open(tmp_path / "final_cat-301-279.fits") as hdul: + data = hdul["RESULTS"].data + assert "VIGNET" not in data.names + npt.assert_array_equal( + data["TILE_UNIQUE_ID"], 301279 * 10**6 + np.array([1, 2, 7]) + ) + npt.assert_allclose(data["TILE_ID"], 301.279) + + +# --- the segmentation map ------------------------------------------------- + + +def _seg_map(): + """Two footprints, labelled 7 and 9, on a 20x20 sky.""" + seg = np.zeros((20, 20), dtype=np.int32) + seg[2:6, 2:6] = 7 + seg[12:18, 12:18] = 9 + return seg + + +@pytest.mark.decision("detection.catalogue_neighbour_marking") +class TestRelabel: + """The relabelled map carries each object's NUMBER on its own footprint. + + Segmentation labels are not the catalogue's NUMBER; the claim is by the + pixel under the object's position, which is all both maps share. + """ + + def test_each_object_owns_the_footprint_it_sits_in(self): + seg = _seg_map() + # Positions in an order that is NOT the label order, so a + # relabelling that merely renumbered would fail. + out, counts = rs.relabel_seg( + seg, number=np.array([1, 2]), + x_image=np.array([15.0, 4.0]), y_image=np.array([15.0, 4.0])) + assert counts["matched"] == 2 + assert set(np.unique(out[seg == 9])) == {1} + assert set(np.unique(out[seg == 7])) == {2} + assert np.all(out[seg == 0] == 0) + + def test_an_unclaimed_footprint_becomes_a_neighbour(self): + seg = _seg_map() + out, _ = rs.relabel_seg( + seg, number=np.array([1]), + x_image=np.array([4.0]), y_image=np.array([4.0])) + assert set(np.unique(out[seg == 9])) == {rs.NEIGHBOUR_LABEL} + assert set(np.unique(out[seg == 7])) == {1} + + def test_an_object_on_sky_gets_a_disc(self): + seg = _seg_map() + out, counts = rs.relabel_seg( + seg, number=np.array([1, 5]), + x_image=np.array([4.0, 10.0]), y_image=np.array([4.0, 10.0]), + fallback_radius=2) + assert counts["unclaimed"] == 1 + assert out[9, 9] == 5 + assert np.count_nonzero(out == 5) == np.count_nonzero( + np.add.outer(np.arange(-2, 3) ** 2, np.arange(-2, 3) ** 2) <= 4) + + def test_two_objects_in_one_footprint_both_keep_a_centre(self): + seg = _seg_map() + out, counts = rs.relabel_seg( + seg, number=np.array([1, 2]), + x_image=np.array([14.0, 16.0]), y_image=np.array([14.0, 16.0]), + fallback_radius=1) + assert counts == dict(matched=1, unclaimed=0, shared=1, off_image=0, + shared_pixel=0) + assert out[13, 13] == 1 and out[15, 15] == 2 + assert np.count_nonzero(out == 1) > np.count_nonzero(out == 2) + + def test_every_object_is_self_somewhere(self): + rng = np.random.default_rng(0) + seg = np.zeros((60, 60), dtype=np.int32) + for label in range(1, 12): + row, col = rng.integers(0, 55, size=2) + seg[row:row + 5, col:col + 5] = label + number = np.arange(1, 31) + x_image = rng.uniform(1, 60, size=30) + y_image = rng.uniform(1, 60, size=30) + out, counts = rs.relabel_seg(seg, number, x_image, y_image) + assert sum(counts[k] for k in + ("matched", "unclaimed", "shared", "off_image")) == 30 + for num in number: + assert np.any(out == num), f"object {num} has no self pixels" + assert set(np.unique(out)) <= set(number) | {0, rs.NEIGHBOUR_LABEL} + + def test_two_objects_on_one_pixel(self): + seg = _seg_map() + out, counts = rs.relabel_seg( + seg, number=np.array([4, 6]), x_image=np.array([10.0, 10.2]), + y_image=np.array([10.0, 10.1]), fallback_radius=1) + assert counts["shared_pixel"] == 1 + assert out[9, 9] == 4 + assert np.any(out == 6) + + @pytest.mark.parametrize("x, y", [(0.4, 5.0), (5.0, 21.0)]) + def test_a_position_off_the_image_is_counted(self, x, y): + out, counts = rs.relabel_seg( + _seg_map(), number=np.array([1]), x_image=np.array([x]), + y_image=np.array([y])) + assert counts["off_image"] == 1 + assert not np.any(out == 1) + + +# Objects on _seg_map(): 1 in footprint 7, 2 on sky between the footprints; +# footprint 9 is claimed by nobody. +SEG_OBJECTS = [(1, 4.0, 4.0), (2, 10.0, 9.0)] +SEG_STAMP = 21 + + +def _marked(number, x, y, seg=None): + image = np.ones(_seg_map().shape, np.float32) + seg = _seg_map() if seg is None else seg + relabelled, _ = rs.relabel_seg(seg, number, x, y) + return rs._extract_vignets(image, x, y, SEG_STAMP, seg=relabelled, + number=number) + + +def _stamp_of(array, x, y, fill): + """The SEG_STAMP stamp of ``array`` centred on 1-based (x, y).""" + half = SEG_STAMP // 2 + padded = np.pad(array, half, constant_values=fill) + row, col = int(round(y)) - 1, int(round(x)) - 1 + return padded[row:row + SEG_STAMP, col:col + SEG_STAMP] + + +@pytest.mark.decision("detection.catalogue_neighbour_marking") +def test_neighbour_footprints_become_big_and_nothing_else(): + number = np.array([o[0] for o in SEG_OBJECTS]) + x = np.array([o[1] for o in SEG_OBJECTS]) + y = np.array([o[2] for o in SEG_OBJECTS]) + vignets = _marked(number, x, y) + seg = _seg_map() + + # Object 1 in footprint 7: footprint 9 (unclaimed) is a neighbour; its + # own footprint and the sky keep their image values. + s = _stamp_of(seg, x[0], y[0], fill=-99) + v = vignets[0] + assert (v[s == 9] == rs.BIG).all() + assert (v[s == 7] == 1).all() + assert (v[s == -99] == rs.BIG).all() # off the image + # Object 2's disc is a catalogue object's pixels, so a neighbour too; the + # rest of the sky is untouched. + relabelled, _ = rs.relabel_seg(seg, number, x, y) + disc = _stamp_of(relabelled, x[0], y[0], fill=-99) == 2 + assert disc.any() and (v[disc] == rs.BIG).all() + assert (v[(s == 0) & ~disc] == 1).all() + + # Object 2 on sky: every footprint in its stamp is a neighbour; the sky, + # its own centre included, keeps its image values. + s = _stamp_of(seg, x[1], y[1], fill=-99) + v = vignets[1] + assert (v[(s == 7) | (s == 9)] == rs.BIG).all() + assert (v[s == 0] == 1).all() + assert v[SEG_STAMP // 2, SEG_STAMP // 2] == 1 + + +@pytest.mark.decision("detection.catalogue_neighbour_marking") +def test_a_claimed_neighbour_is_masked_by_its_number(): + """With both footprints claimed, each object masks the other's.""" + number, x, y = np.array([3, 8]), np.array([4.0, 15.0]), np.array([4.0, 15.0]) + vignets = _marked(number, x, y) + seg = _seg_map() + for i, (own, other) in enumerate([(7, 9), (9, 7)]): + s = _stamp_of(seg, x[i], y[i], fill=-99) + assert (vignets[i][s == other] == rs.BIG).all() + assert (vignets[i][s == own] == 1).all() + + +@pytest.mark.decision("detection.catalogue_neighbour_marking") +def test_converter_relabels_and_marks_from_compressed_map(tmp_path): + """End to end from a compressed map, as fetched from vos.""" + cat = tmp_path / "CFIS_cat-301-279.cat" + img = tmp_path / "CFIS_image-301-279.fits" + seg_in = tmp_path / "CFIS_seg-301-279.fitsfz" + out = tmp_path / "sexcat-301-279.fits" + lines = ["# 1 NUMBER", "# 2 X_IMAGE", "# 3 Y_IMAGE", + "# 4 ALPHA_J2000", "# 5 DELTA_J2000"] + lines += [f"{n} {x} {y} 150.0 30.0" for n, x, y in SEG_OBJECTS] + cat.write_text("\n".join(lines) + "\n") + fits.PrimaryHDU(np.ones((20, 20), np.float32)).writeto(img) + fits.HDUList([fits.PrimaryHDU(), + fits.CompImageHDU(_seg_map())]).writeto(seg_in) + + rs.make_ldac_from_ascii(str(cat), str(img), str(out), stamp_size=SEG_STAMP, + seg_path=str(seg_in)) + + seg = _seg_map() + number, x, y = (np.array(c) for c in zip(*SEG_OBJECTS)) + relabelled, _ = rs.relabel_seg(seg, number, x, y) + assert set(np.unique(relabelled[seg == 7])) == {1} + assert set(np.unique(relabelled[seg == 9])) == {rs.NEIGHBOUR_LABEL} + with fits.open(out) as hdul: + v = hdul["LDAC_OBJECTS"].data["VIGNET"][0] + s = _stamp_of(seg, 4.0, 4.0, fill=-99) + assert (v[s == 9] == rs.BIG).all() and (v[s == 7] == 1).all() + + fits.HDUList([fits.PrimaryHDU(), + fits.CompImageHDU(_seg_map()[:10])]).writeto( + seg_in, overwrite=True) + with pytest.raises(ValueError, match="one grid"): + rs.make_ldac_from_ascii(str(cat), str(img), str(out), + stamp_size=SEG_STAMP, seg_path=str(seg_in)) + + +DR6_PATCH = Path(__file__).parent / "data" / "dr6_202.301_seg_patch.fits" + + +@pytest.mark.decision("detection.catalogue_neighbour_marking") +def test_dr6_marks_every_neighbour_pixel_and_no_own_pixel(): + """On a crowded 200x200 patch of the real 202.301 map and catalogue. + + The raw labels are not NUMBER, so this checks the whole chain on real + data: 1-based positions, the centre-pixel claim, and the stamp geometry. + Every stamp fully on the patch is checked against the RAW map: pixels of + footprints other than the one under the object are all -1e30, its own + footprint and the sky are untouched. (SExtractor's own VIGNET, on a + SExtractor seg map of an image-sim tile, marks 93% of neighbour-label + pixels and 0.14% of own pixels: check_sex_vignet.py.) + """ + with fits.open(DR6_PATCH) as hdul: + seg = hdul["SEG"].data + objects = hdul["OBJECTS"].data + number = np.array(objects["NUMBER"]) + x, y = np.array(objects["X_IMAGE"]), np.array(objects["Y_IMAGE"]) + relabelled, counts = rs.relabel_seg(seg, number, x, y) + assert counts["matched"] == len(number) + + stamp, half = 51, 25 + image = np.ones(seg.shape, np.float32) + vignets = rs._extract_vignets(image, x, y, stamp, seg=relabelled, + number=number) + col, row = np.rint(x).astype(int) - 1, np.rint(y).astype(int) - 1 + full = ((col >= half) & (col < seg.shape[1] - half) + & (row >= half) & (row < seg.shape[0] - half)) + assert full.sum() >= 10 + + marked = {"neighbour": [0, 0], "own": [0, 0], "sky": [0, 0]} + for i in np.flatnonzero(full): + s = seg[row[i] - half:row[i] + half + 1, col[i] - half:col[i] + half + 1] + big = vignets[i] == rs.BIG + own = s[half, half] + for kind, where in (("neighbour", (s != 0) & (s != own)), + ("own", s == own), ("sky", s == 0)): + marked[kind][0] += (big & where).sum() + marked[kind][1] += where.sum() + assert marked["neighbour"][1] > 1000 + assert marked["neighbour"][0] == marked["neighbour"][1] + assert marked["own"][0] == 0 + assert marked["sky"][0] == 0 + + +# --- the runner's output against the completeness table ------------------- + + +def test_runner_output_matches_tile_detect_completeness(tmp_path, monkeypatch): + """The runner writes exactly the files ``tile_detect`` expects. + + Runs the real runner, segmentation map on, into a run dir and checks it + with ``completeness.check_counts`` under ``SP_TILE_DETECTION=unions_catalogue``, + so the table and the converter's outputs cannot drift apart. + """ + import configparser + import importlib.util + import logging + + from shapepipe.modules.read_ext_sexcat_runner import read_ext_sexcat_runner + + scripts = Path(__file__).resolve().parents[2] / "workflow" / "scripts" + spec = importlib.util.spec_from_file_location( + "_completeness", scripts / "completeness.py") + completeness = importlib.util.module_from_spec(spec) + spec.loader.exec_module(completeness) + + cat = tmp_path / "CFIS_cat-301-279.cat" + img = tmp_path / "CFIS_image-301-279.fits" + seg_in = tmp_path / "CFIS_seg-301-279.fitsfz" + lines = ["# 1 NUMBER", "# 2 X_IMAGE", "# 3 Y_IMAGE", + "# 4 ALPHA_J2000", "# 5 DELTA_J2000"] + lines += [f"{n} {x} {y} 150.0 30.0" for n, x, y in SEG_OBJECTS] + cat.write_text("\n".join(lines) + "\n") + fits.PrimaryHDU(np.ones((20, 20), np.float32)).writeto(img) + fits.HDUList([fits.PrimaryHDU(), + fits.CompImageHDU(_seg_map())]).writeto(seg_in) + + run_dir = tmp_path / "run_sp_tile_Rx" + out_dir = run_dir / "read_ext_sexcat_runner" / "output" + out_dir.mkdir(parents=True) + config = configparser.ConfigParser() + config["READ_EXT_SEXCAT_RUNNER"] = { + "SEGMENTATION": "True", "MAKE_POST_PROCESS": "False", + "VIGNET_SIZE": str(SEG_STAMP), + } + read_ext_sexcat_runner( + [str(cat), str(img), str(seg_in)], {"output": str(out_dir)}, + "-301-279", config, "READ_EXT_SEXCAT_RUNNER", + logging.getLogger("test"), + ) + + monkeypatch.setenv("SP_TILE_DETECTION", "unions_catalogue") + table = completeness.COMPLETENESS["tile_detect"]["unions_catalogue"] + expect = table["read_ext_sexcat_runner"]["expect"] + written = sorted(p.name for p in out_dir.iterdir()) + assert len(written) == expect, written + ok, details = completeness.check_counts("tile_detect", run_dir) + assert ok, details diff --git a/tests/unit/test_final_cat_merge_invariants.py b/tests/unit/test_final_cat_merge_invariants.py index 7dcfc6183..21eb6adad 100644 --- a/tests/unit/test_final_cat_merge_invariants.py +++ b/tests/unit/test_final_cat_merge_invariants.py @@ -130,7 +130,7 @@ def _campaign(root: Path, drop=None): argv = [sys.executable, str(SCRIPT), "--products-dir", str(products), "--tile-list", str(tile_list), "--index-db", str(index), "--output", str(output), "--campaign", CAMPAIGN, - "--param-file", str(param)] + "--param-file", str(param), "--tile-detection", "sextractor"] return argv, output, sources diff --git a/tests/unit/test_ngmix_range.py b/tests/unit/test_ngmix_range.py index 30027cd29..fe2bea719 100644 --- a/tests/unit/test_ngmix_range.py +++ b/tests/unit/test_ngmix_range.py @@ -77,7 +77,7 @@ def test_ranges_tile_the_catalogue(n_obj, n_chunks, data): max_size=n_obj, ) ) - assert_tiles(ngmix_range.id_ranges(epochs, n_chunks), n_obj, n_chunks) + assert_tiles(ngmix_range.row_ranges(epochs, n_chunks), n_obj, n_chunks) @settings(deadline=None) @@ -99,8 +99,8 @@ def test_ranges_are_deterministic(n_obj, n_chunks, data): max_size=n_obj, ) ) - first = ngmix_range.id_ranges(epochs, n_chunks) - assert first == ngmix_range.id_ranges(list(epochs), n_chunks) + first = ngmix_range.row_ranges(epochs, n_chunks) + assert first == ngmix_range.row_ranges(list(epochs), n_chunks) assert all(isinstance(b, int) for lo, hi in first for b in (lo, hi)) @@ -124,7 +124,7 @@ def test_written_rows_reproduce_the_per_chunk_computation(n_obj, n_chunks, data) """THE INVARIANCE TEST for materialising the split. Writing the whole partition and reading chunk ``k``'s row back must give - byte-for-byte what the old per-chunk ``id_ranges(...)[k - 1]`` returned. If + byte-for-byte what the old per-chunk ``row_ranges(...)[k - 1]`` returned. If this ever fails, the two mechanisms have drifted and a tile's coverage is what pays. """ @@ -137,7 +137,7 @@ def test_written_rows_reproduce_the_per_chunk_computation(n_obj, n_chunks, data) ) doc = _roundtrip(epochs, n_chunks) assert doc["n_obj"] == n_obj and doc["n_chunks"] == n_chunks - expected = ngmix_range.id_ranges(epochs, n_chunks) + expected = ngmix_range.row_ranges(epochs, n_chunks) got = [ngmix_range.chunk_range(doc, k) for k in range(1, n_chunks + 1)] assert got == expected assert_tiles(got, n_obj, n_chunks) @@ -167,12 +167,12 @@ def test_reading_an_absent_ranges_file_is_fatal(tmp_path): def test_single_chunk_takes_everything(): """n_chunks == 1: one range over the whole catalogue.""" - assert ngmix_range.id_ranges([3] * 17, 1) == [(1, 17)] + assert ngmix_range.row_ranges([3] * 17, 1) == [(1, 17)] def test_one_object_per_chunk(): """n_obj == n_chunks: one object each, however lopsided the weights.""" - assert ngmix_range.id_ranges([0, 9, 1, 40], 4) == [ + assert ngmix_range.row_ranges([0, 9, 1, 40], 4) == [ (1, 1), (2, 2), (3, 3), (4, 4) ] @@ -184,27 +184,27 @@ def test_fewer_objects_than_chunks_pads_with_empty_ranges(): reads ``ID_OBJ_MAX = 0`` as unbounded — so each of them would have measured the entire tile rather than nothing. """ - assert ngmix_range.id_ranges([2, 5, 1], 6) == [ + assert ngmix_range.row_ranges([2, 5, 1], 6) == [ (1, 1), (2, 2), (3, 3), (4, 3), (4, 3), (4, 3) ] - assert_tiles(ngmix_range.id_ranges([2, 5, 1], 6), 3, 6) + assert_tiles(ngmix_range.row_ranges([2, 5, 1], 6), 3, 6) def test_zero_objects_is_fatal(): """No split is meaningful, and (1, 0) is ngmix's unbounded sentinel.""" with pytest.raises(ValueError, match="zero objects"): - ngmix_range.id_ranges([], 8) + ngmix_range.row_ranges([], 8) def test_zero_chunks_is_fatal(): """There is no zeroth chunk to hand a range to.""" with pytest.raises(ValueError, match="n_chunks"): - ngmix_range.id_ranges([1, 2, 3], 0) + ngmix_range.row_ranges([1, 2, 3], 0) def test_equal_weights_reproduce_an_equal_count_split(): """Uniform epochs: chunk sizes differ by at most one object.""" - ranges = ngmix_range.id_ranges([3] * 1000, 8) + ranges = ngmix_range.row_ranges([3] * 1000, 8) assert_tiles(ranges, 1000, 8) sizes = [hi - lo + 1 for lo, hi in ranges] assert max(sizes) - min(sizes) <= 1 @@ -217,7 +217,7 @@ def test_zero_epoch_objects_still_weigh_something(): Ninety-six zero-epoch objects and four 1-epoch ones. Weighing only epochs would let a single chunk swallow all ninety-six. """ - ranges = ngmix_range.id_ranges([0] * 96 + [1] * 4, 4) + ranges = ngmix_range.row_ranges([0] * 96 + [1] * 4, 4) assert_tiles(ranges, 100, 4) assert max(hi - lo + 1 for lo, hi in ranges) < 96 @@ -225,7 +225,7 @@ def test_zero_epoch_objects_still_weigh_something(): def test_one_enormously_heavy_object_is_isolated(): """The heavy object gets a chunk to itself; the tail still tiles.""" epochs = [1] * 20 + [10_000] + [1] * 20 - ranges = ngmix_range.id_ranges(epochs, 4) + ranges = ngmix_range.row_ranges(epochs, 4) assert_tiles(ranges, 41, 4) assert (21, 21) in ranges @@ -268,6 +268,6 @@ def test_slowest_chunk_is_minimal(epochs, n_chunks): ngmix_range.MILLI_EPOCH * e + ngmix_range.ALPHA_MILLI_EPOCHS for e in epochs ] - ranges = ngmix_range.id_ranges(epochs, n_chunks) + ranges = ngmix_range.row_ranges(epochs, n_chunks) loads = [sum(weights[lo - 1:hi]) for lo, hi in ranges] assert max(loads) == _brute_force_min_max(weights, n_chunks) diff --git a/tests/unit/test_run_config.py b/tests/unit/test_run_config.py index 3290f9fc3..2fea284ba 100644 --- a/tests/unit/test_run_config.py +++ b/tests/unit/test_run_config.py @@ -1,27 +1,119 @@ -"""``workflow/scripts/run_config.py``: `run:` is a required key, and every -machine-key path expands fully or is reported. +"""``workflow/scripts/run_config.py``: the resolver's layering (config.yaml < +input_types < machines < run config, with $variables expanded across every +layer), `run:` as a required key, and every path expanding fully or being +reported. `run:` names the campaign's merged catalogues (``final_cat_.hdf5``, -``full_starcat_.hdf5``). ``unresolved()`` already reports a ``$run`` left -unexpanded in a path; these tests pin that a run config whose paths never -mention ``$run`` is refused too, and that the shipped ``config.yaml`` leaves the -name to the run config. A ``$base_dir`` whose value holds ``$run`` expands -through both, and an optional path left holding a ``$`` is reported. +``full_starcat_.hdf5``), so a run config whose paths never mention +``$run`` is refused too, and the shipped ``config.yaml`` leaves the name to +the run config. + +Container-free: run_config.py needs only PyYAML. """ import importlib.util from pathlib import Path -import yaml +import pytest + +yaml = pytest.importorskip("yaml") REPO_ROOT = Path(__file__).resolve().parents[2] -SCRIPT = REPO_ROOT / "workflow" / "scripts" / "run_config.py" CONFIG_YAML = REPO_ROOT / "workflow" / "config.yaml" -_spec = importlib.util.spec_from_file_location("_run_config", SCRIPT) +_spec = importlib.util.spec_from_file_location( + "_run_config", REPO_ROOT / "workflow" / "scripts" / "run_config.py") run_config = importlib.util.module_from_spec(_spec) _spec.loader.exec_module(run_config) +BASE = { + "input_type": "data", + "run": "r1", + "input_types": { + "data": {"psf_model": "psfex", "tile_detection": "unions_catalogue"}, + "image_sims": {"psf_model": "fake", "tile_detection": "sextractor"}, + }, + "machines": { + "m": { + "base_dir": "/base", + "data": { + "tile_list": "$base_dir/$run/tiles.txt", + "inputs": {"tiles": "$base_dir/tiles", + "exposures": "$base_dir/exp"}, + "outputs": {"run_dir": "/scratch/$run", + "index_db": "/idx/$run.sqlite"}, + }, + "image_sims": {"psf_model": "mccd", + "inputs": {"tiles": "/sims/$run"}}, + }, + }, +} + + +def _resolve(tmp_path, monkeypatch, over, base=BASE): + monkeypatch.setenv("SP_PROFILE", "m") + base_path = tmp_path / "config.yaml" + base_path.write_text(yaml.safe_dump(base)) + run_path = tmp_path / "run.yaml" + run_path.write_text(yaml.safe_dump(over)) + return run_config.load(str(base_path), str(run_path)) + + +def test_input_type_table_supplies_its_defaults(tmp_path, monkeypatch): + config = _resolve(tmp_path, monkeypatch, {}) + assert config["psf_model"] == "psfex" + assert config["tile_detection"] == "unions_catalogue" + + +def test_machine_entry_beats_input_type_table(tmp_path, monkeypatch): + config = _resolve(tmp_path, monkeypatch, {"input_type": "image_sims"}) + assert config["psf_model"] == "mccd" + assert config["tile_detection"] == "sextractor" + + +def test_run_config_beats_both_tables(tmp_path, monkeypatch): + config = _resolve(tmp_path, monkeypatch, { + "input_type": "image_sims", "psf_model": "psfex", + "tile_detection": "unions_catalogue"}) + assert config["psf_model"] == "psfex" + assert config["tile_detection"] == "unions_catalogue" + + +def test_dicts_merge_per_subkey_and_expand(tmp_path, monkeypatch): + config = _resolve(tmp_path, monkeypatch, { + "run": "r2", "inputs": {"exposures": "$base_dir/own_exp"}}) + assert config["inputs"] == {"tiles": "/base/tiles", + "exposures": "/base/own_exp"} + assert config["tile_list"] == "/base/r2/tiles.txt" + assert run_config.unresolved(config) == [] + + +def test_input_type_table_values_expand(tmp_path, monkeypatch): + base = dict(BASE, input_types={"data": {"psf_dict": "$base_dir/$run.pkl"}}) + config = _resolve(tmp_path, monkeypatch, {}, base=base) + assert config["psf_dict"] == "/base/r1.pkl" + + +def test_unexpanded_variable_is_unresolved(tmp_path, monkeypatch): + config = _resolve(tmp_path, monkeypatch, + {"outputs": {"run_dir": "$nowhere/x"}}) + assert run_config.unresolved(config) == ["outputs.run_dir"] + + +def test_committed_top_level_sets_no_table_key(): + """config.yaml's top level sits beneath both tables only because it sets + none of their keys: one that did would shadow them, since the tables fill + only what is unset.""" + config = yaml.safe_load(CONFIG_YAML.read_text()) + table_keys = set() + for entry in config["input_types"].values(): + table_keys |= set(entry) + for machine in config["machines"].values(): + for input_type in config["input_types"]: + table_keys |= set(machine.get(input_type) or {}) + assert not table_keys & set(config) + + # Every other REQUIRED key, as literal paths with no `$run` in them. LITERAL = { "tile_list": "/data/tiles.txt", @@ -70,14 +162,14 @@ def _machine_config(base_dir, **over): def test_base_dir_holding_run_expands_fully(): - cfg = run_config.apply_machine_defaults(_machine_config("/b/$run")) + cfg = run_config.apply_defaults(_machine_config("/b/$run")) assert cfg["outputs"]["run_dir"] == "/b/smk-g6/run" assert run_config.unresolved(cfg) == [] def test_run_config_shorthands_nest_and_run_is_written_back(): """Top-level scalars are variables, ${name} delimits, nesting resolves.""" - cfg = run_config.apply_machine_defaults(_machine_config( + cfg = run_config.apply_defaults(_machine_config( "/b", shear="1p2z", grid="grid_2", run="${shear}_${grid}", outputs={"run_dir": "/o/$run/scratch"})) assert cfg["run"] == "1p2z_grid_2" @@ -86,12 +178,12 @@ def test_run_config_shorthands_nest_and_run_is_written_back(): def test_run_holding_an_unknown_variable_is_reported(): - cfg = run_config.apply_machine_defaults(_machine_config("/b", run="$nope")) + cfg = run_config.apply_defaults(_machine_config("/b", run="$nope")) assert "run" in run_config.unresolved(cfg) def test_optional_path_with_an_unknown_variable_is_reported(): - cfg = run_config.apply_machine_defaults(_machine_config( + cfg = run_config.apply_defaults(_machine_config( "/b", outputs={"products_dir": "/p/$nope/products"}, inputs={"masks": "/m/$nope"}, container="/c/$nope.sif")) assert set(run_config.unresolved(cfg)) == { diff --git a/tests/unit/test_tile_unique_id.py b/tests/unit/test_tile_unique_id.py new file mode 100644 index 000000000..15cb9125d --- /dev/null +++ b/tests/unit/test_tile_unique_id.py @@ -0,0 +1,49 @@ +"""Survey-wide object ID: ``tile_id * 10**6 + NUMBER``.""" + +import numpy as np +import pytest + +from shapepipe.utilities import cfis + + +@pytest.mark.parametrize( + "name", ["-301-279", "CFIS.301.279.r", "CFIS.301.279.r.weight.fits.fz"] +) +def test_tile_id_from_name(name): + assert cfis.get_tile_id(name) == 301279 + + +def test_tile_id_keeps_leading_zeros(): + assert cfis.get_tile_id("-004-012") == 4012 + + +@pytest.mark.parametrize("name", ["-1-2", "-1301-279", "-301-2790"]) +def test_tile_id_rejects_non_three_digit_components(name): + with pytest.raises(cfis.CfisError): + cfis.get_tile_id(name) + + +def test_unique_id_example(): + assert cfis.get_tile_unique_id(301279, 42) == 301279000042 + + +def test_unique_id_round_trip(): + tile_id = 999999 + number = np.array([0, 1, 12345, 999999]) + unique_id = cfis.get_tile_unique_id(tile_id, number) + assert unique_id.dtype == np.int64 + assert len(np.unique(unique_id)) == len(number) + tile_back, number_back = cfis.split_tile_unique_id(unique_id) + np.testing.assert_array_equal(tile_back, tile_id) + np.testing.assert_array_equal(number_back, number) + + +@pytest.mark.parametrize("number", [10**6, -1, [5, 10**6]]) +def test_unique_id_rejects_number_out_of_range(number): + with pytest.raises(cfis.CfisError): + cfis.get_tile_unique_id(301279, number) + + +def test_unique_id_rejects_tile_id_out_of_range(): + with pytest.raises(cfis.CfisError): + cfis.get_tile_unique_id(10**6, 1) diff --git a/tests/unit/test_workflow_tile_detection.py b/tests/unit/test_workflow_tile_detection.py new file mode 100644 index 000000000..24775304e --- /dev/null +++ b/tests/unit/test_workflow_tile_detection.py @@ -0,0 +1,178 @@ +"""The two tile_detect modes agree on where the sexcat is. + +``tile_detection: unions_catalogue`` swaps the SExtractor rule for a fetch plus +a conversion, and everything downstream of ``tile_detect`` is unchanged only +because both modes write ``run_sp_tile_Sx`` and both are checked as the +``tile_detect`` stage. The agreements that make that true are between files +that never see each other at run time: the two inis' RUN_NAMEs, +``completeness.STAGE_DIR``, the flavoured ``COMPLETENESS['tile_detect']`` +table, and ``run_report``'s stage list. They are asserted here, statically. +Container-free: the scripts are stdlib-only. +""" + +import configparser +import importlib.util +import sys +from pathlib import Path + +import pytest + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPTS = REPO_ROOT / "workflow" / "scripts" +CONFIG_DIR = REPO_ROOT / "workflow" / "config" / "cfis" + + +def _load(name): + sys.path.insert(0, str(SCRIPTS)) + try: + spec = importlib.util.spec_from_file_location(f"_{name}", SCRIPTS / f"{name}.py") + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + finally: + sys.path.remove(str(SCRIPTS)) + return module + + +def _ini(name): + parser = configparser.ConfigParser() + assert parser.read(CONFIG_DIR / name) == [str(CONFIG_DIR / name)] + return parser + + +completeness = _load("completeness") + + +def test_modes_are_the_two_the_rules_know(): + assert completeness.TILE_DETECTIONS == ("sextractor", "unions_catalogue") + assert set(completeness.COMPLETENESS["tile_detect"]) == set( + completeness.TILE_DETECTIONS) + + +def test_both_detection_inis_write_the_stage_dir(): + """config_tile_Sx and config_tile_Uc share the run dir STAGE_DIR names.""" + level, subdir = completeness.STAGE_DIR["tile_detect"] + assert level == "tile" + for name in ("config_tile_Sx.ini", "config_tile_Uc.ini"): + assert _ini(name)["DEFAULT"]["RUN_NAME"].strip() == subdir, ( + f"{name} writes a run dir other than {subdir}: unit_pre would clear " + "the wrong directory and the chain downstream would read nothing.") + assert _ini("config_tile_Uc.ini")["DEFAULT"]["RUN_DATETIME"].strip() == "False" + + +def test_fetch_ini_matches_its_stage_and_feeds_the_converter(): + level, gic = completeness.STAGE_DIR["tile_get_catalogue"] + assert level == "tile" + assert _ini("config_tile_Gic.ini")["DEFAULT"]["RUN_NAME"].strip() == gic + uc = _ini("config_tile_Uc.ini")["READ_EXT_SEXCAT_RUNNER"] + assert f"$SP_RUN/output/{gic}/get_images_runner/output" in uc["INPUT_DIR"] + assert uc["FILE_PATTERN"].split(",")[0].strip() == "CFIS_cat" + # The segmentation map rides the same fetch and is the converter's third + # input, the position SEGMENTATION = True reads it from. + gi = _ini("config_tile_Gic.ini")["GET_IMAGES_RUNNER"] + assert [p.strip() for p in gi["OUTPUT_FILE_PATTERN"].split(",")] == [ + "CFIS_cat-", "CFIS_seg-"] + assert [e.strip() for e in gi["INPUT_FILE_EXT"].split(",")] == [ + ".cat", ".fits.fz"] + third = [[v.strip() for v in uc[k].split(",")][2] + for k in ("INPUT_DIR", "FILE_PATTERN", "FILE_EXT")] + assert third == [f"$SP_RUN/output/{gic}/get_images_runner/output", + "CFIS_seg", ".fitsfz"] + assert uc["SEGMENTATION"].strip() == "True" + # The multi-epoch post-processing is what gives the sexcat its EPOCH_k + # extensions; ngmix_range.py refuses a sexcat without them. + assert uc["MAKE_POST_PROCESS"].strip() == "True" + + +def _stage_dir(tmp_path, runner, n): + out = tmp_path / runner / "output" + out.mkdir(parents=True) + for k in range(n): + (out / f"f{k}").write_text("") + return tmp_path + + +@pytest.mark.parametrize("mode, runner, expect", [ + ("sextractor", "sextractor_runner", 2), + ("unions_catalogue", "read_ext_sexcat_runner", 1), +]) +def test_tile_detect_is_checked_per_mode(tmp_path, monkeypatch, mode, runner, expect): + monkeypatch.setenv("SP_TILE_DETECTION", mode) + ok, details = completeness.check_counts( + "tile_detect", _stage_dir(tmp_path, runner, expect)) + assert ok and details == [(runner, expect, expect, False)] + + +def test_unset_mode_is_sextractor(tmp_path, monkeypatch): + """A prologue without the export checks as before: data runs are unchanged.""" + monkeypatch.delenv("SP_TILE_DETECTION", raising=False) + ok, details = completeness.check_counts( + "tile_detect", _stage_dir(tmp_path, "sextractor_runner", 2)) + assert ok and details[0][0] == "sextractor_runner" + + +def test_invalid_mode_is_fatal(tmp_path, monkeypatch): + monkeypatch.setenv("SP_TILE_DETECTION", "steven") + with pytest.raises(ValueError, match="SP_TILE_DETECTION"): + completeness.check_counts("tile_detect", tmp_path) + + +@pytest.mark.parametrize("mode, present", [ + ("sextractor", False), ("unions_catalogue", True), (None, False), +]) +def test_report_lists_the_fetch_stage_only_for_catalogue_runs(monkeypatch, mode, present): + if mode is None: + monkeypatch.delenv("SP_TILE_DETECTION", raising=False) + else: + monkeypatch.setenv("SP_TILE_DETECTION", mode) + stages = _load("run_report").TILE_STAGES + assert ("tile_get_catalogue" in stages) is present + if present: + assert stages.index("tile_get_catalogue") == stages.index("tile_detect") - 1 + + +@pytest.mark.parametrize("machine", ["nibi", "candide"]) +@pytest.mark.parametrize("input_type", ["data", "image_sims"]) +def test_defaults_pair_the_catalogue_with_its_source( + machine, input_type, monkeypatch, tmp_path): + """Real data defaults to the catalogue, and every machine says where it is. + + tile_detection follows from the input type (the input_types: table), but + `tile_detection: unions_catalogue` is refused by the Snakefile without + `inputs.catalogues`, which is a per-machine path: every machine with a + data entry has to declare one. Image sims have no UNIONS catalogue and + use SExtractor. + """ + pytest.importorskip("yaml") + run_config = _load("run_config") + monkeypatch.setenv("SP_PROFILE", machine) + over = tmp_path / "run.yaml" + over.write_text(f"input_type: {input_type}\n") + config = run_config.load( + str(REPO_ROOT / "workflow" / "config.yaml"), str(over)) + if input_type == "data": + assert config["tile_detection"] == "unions_catalogue" + assert run_config.catalogue_source(config)[0] + else: + assert config["tile_detection"] == "sextractor" + + +@pytest.mark.parametrize("catalogues", [None, "", "TBD"]) +def test_catalogue_run_without_a_source_fails_at_parse(catalogues): + """An unset or placeholder `inputs.catalogues` is refused before any job runs.""" + run_config = _load("run_config") + config = {"tile_detection": "unions_catalogue", + "inputs": {} if catalogues is None else {"catalogues": catalogues}} + with pytest.raises(ValueError, match="inputs.catalogues"): + run_config.catalogue_source(config) + config["tile_detection"] = "sextractor" + assert run_config.catalogue_source(config) == ("", "symlink") + + +@pytest.mark.parametrize("source, retrieve", [ + ("vos:cfis/tiles_DR6", "vos"), ("/data/tiles_DR6", "symlink"), +]) +def test_catalogue_retrieve_follows_the_prefix(source, retrieve): + run_config = _load("run_config") + config = {"tile_detection": "unions_catalogue", + "inputs": {"catalogues": source}} + assert run_config.catalogue_source(config) == (source, retrieve) diff --git a/tests/workflow/README.md b/tests/workflow/README.md index 062b4ce49..c3576f08c 100644 --- a/tests/workflow/README.md +++ b/tests/workflow/README.md @@ -28,6 +28,7 @@ Scratch and persistent roots differ, and the products directory's basename delib The three modes are `data+psfex`, `data+mccd`, and `image_sims+fake`. `final_cat_merge` is present in all three: it merges galaxy catalogues regardless of the PSF source. Only `exp_persist` and `star_cat_merge` disappear for fake PSFs. +The run config leaves `tile_detection` to config.yaml's `input_types:` table, so `data` plans the UNIONS-catalogue detection (adding `tile_get_catalogue`, fed from a fixture-local `inputs.catalogues`) and `image_sims` plans SExtractor. `resolve()` uses the same `SP_PHASE`, `SP_PROFILE`, `SP_RUN_CONFIG`, image selection, and state-directory conventions as `workflow/bin/sp`. It isolates the environment, source cache, bare script imports, and Snakemake's shared global namespace. diff --git a/tests/workflow/harness.py b/tests/workflow/harness.py index e3e487b28..1c8420f26 100644 --- a/tests/workflow/harness.py +++ b/tests/workflow/harness.py @@ -88,6 +88,7 @@ def __post_init__(self): "inputs": { "tiles": "$base_dir/inputs/$run/tiles", "exposures": "$base_dir/inputs/$run/exposures", + "catalogues": "$base_dir/inputs/$run/catalogues", }, "outputs": { "run_dir": "$base_dir/scratch/$run", @@ -115,6 +116,16 @@ def __post_init__(self): } self.write_config() + @property + def tile_detection(self): + """Return the committed tile_detection default for this input_type. + + The run config leaves it unset, so the DAG follows config.yaml's + input_types: table and a change to that default moves the pin. + """ + committed = yaml.safe_load((REPO / "workflow" / "config.yaml").read_text()) + return committed["input_types"][self.input_type]["tile_detection"] + def _manifest(self, tile, stage): path = self.tile_manifest(tile, stage) path.parent.mkdir(parents=True, exist_ok=True) diff --git a/tests/workflow/params_pin.json b/tests/workflow/params_pin.json index 1924265dc..eeb61d950 100644 --- a/tests/workflow/params_pin.json +++ b/tests/workflow/params_pin.json @@ -8,12 +8,13 @@ "exp_persist": "302e2837542bc1102430c27c81c600b7cda32e8bddcb5fd60d33950987609fff", "exp_psf": "2c4f6d00f1939ccbf05b4982a0202a0ff92a4373a727f4e00f3aaab4eba03352", "exp_split": "6e954f8f3d06f44d3f164675912ce27d6216648855d04968f9168bd7d0f2c4fa", - "final_cat_merge": "e7f46859c4503a2220713d7bb2507555515d0a9632d780b20f14c59e32210023", + "final_cat_merge": "f6222c98be9777cbab4131d8857bd3817719df9302d0a321c6cd99f67d5bdc76", "prepare_all_tiles": "b8f872a22adf014e25a7fa5198f49b71a6fe9e56042ed82b682bc8763970a844", "star_cat_merge": "6277450958474af5270982fa35360f2f237a29f7533c526ee9265dfd5acc07a0", - "tile_detect": "457456f4f71a500e4166a63933b88d14eed34a019a9814f44546e137b0e71d7a", + "tile_detect": "1b4b0871bf28296f8e5d6855ae47ff2cb84542c93593a5ae5c0ffc3c84c45479", "tile_exp_forest": "7447ab4a1049de5f0b5c81e5f9ed2a8644c7bdab85cb0a060bfde81261a89b28", "tile_find_exposures": "8704317871744996c44351c2836fcb222d7986a602e9e046a90d684ba7b3c184", + "tile_get_catalogue": "7e3f889a955a14b0b917a015c2a8bc90e433a89b1cc5e00b8ca94695c4c93d3c", "tile_get_images": "331a67e747f211ebf4c14b946a7af7f9fc9f55243d69d2791d74aecc3ca228c3", "tile_make_cat": "a57518b04c11f70bf41b320532fddffdd1115fdac604d7dd3928eb61cca28f23", "tile_merge_cats": "ff21216ea804dccc2d2c290d2b2499d5d05f0c34c0a56993c233f43fe3c06bdb", @@ -23,13 +24,14 @@ "tile_vignets": "9d4ae0d99c18217f2185f245281312454c8a219ec1628176e08c71a5efc4dc91" }, "schema": 1, - "sha256": "74adbcd14bc304f82c3f2d1a50374bc607e5313fa5f3f8b3db6302653bad67ea", + "sha256": "f9772987e517e9827c8d501d310ec323be251a62b7b89aad0725a03f8f4f120e", "unit_pre": { "exp_get_images": "8dec850af212879f225fcf27a5f1281e1a075264c7b97d38c2214395d360168c", "exp_psf": "f2358ddf7385918dc5033d10b37f6dc97a15d02b071a3ea0a4619a5f7e6f5bec", "exp_split": "358fa8bbe59680d4f9839e007b343dd25fed733cf3156dccd30d305d33ac9480", "tile_detect": "adad5d671fa65dd04433e1c82a725b845635e9d2e368a1e70833ed990df2d55a", "tile_find_exposures": "c5922fb507f6fd040a179b53fba0818697661016c9688984b6ce249186dc6986", + "tile_get_catalogue": "b124d12252617abb8df8a98d6234ddfd6cc503e19fa72c9a0c48c21559d6f48c", "tile_get_images": "45b47c44bfeb34ac973b89d4e752c8c82e78028f9f0da379c44627005b45e279", "tile_make_cat": "7579a52e32e76523c0b77d75c468a47d7f2e2ca0c55865f40e9628a7481cd79d", "tile_merge_cats": "69cfe94ba941d2ba7b1ce24883961b47b191b80aff381da1865c6bc38de0c354", diff --git a/tests/workflow/test_dag.py b/tests/workflow/test_dag.py index e582138bc..46ae93fc0 100644 --- a/tests/workflow/test_dag.py +++ b/tests/workflow/test_dag.py @@ -16,13 +16,17 @@ "final_cat_merge", } PSF_RULES = {"exp_persist", "star_cat_merge"} +CATALOGUE_RULES = {"tile_get_catalogue"} def test_rule_set_matches_input_mode(campaign, dag): - """A PSF gate cannot remove real-PSF products or add them to fake PSFs.""" + """A PSF gate cannot remove real-PSF products or add them to fake PSFs; + the UNIONS-catalogue detection adds its fetch rule and nothing else.""" expected = BASE_RULES.copy() if campaign.psf_model != "fake": expected |= PSF_RULES + if campaign.tile_detection == "unions_catalogue": + expected |= CATALOGUE_RULES assert dag.rule_names == expected assert "merge_final_cats" not in dag.declared_rule_names diff --git a/universes/committed.yaml b/universes/committed.yaml index 63b6e795b..bc34d7601 100644 --- a/universes/committed.yaml +++ b/universes/committed.yaml @@ -12,6 +12,7 @@ analyses: sky_mask_application: deferred_downstream detection: decisions: + tile_detection: unions_catalogue detection_threshold_policy: megapipe_tiles deblending_policy: megapipe_tiles background_model: auto_megapipe_tiles @@ -23,6 +24,7 @@ analyses: photometry_parameters: kron_25_35 detection_source_mode: sx_nomask_single_image epoch_membership_ccd_bounds: trimmed_bounds_33_2080 + catalogue_neighbour_marking: segmentation_map preparation: decisions: astrometric_solution_source: delivered_headers diff --git a/workflow/README.md b/workflow/README.md index 02eb6c418..6e91b344f 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -23,10 +23,17 @@ source /project/def-mjhudson/cdaley/snakemake-env/bin/activate uv pip install 'snakemake>=9,<10' 'snakemake-executor-plugin-slurm>=2.7,<3' # Write a run config (see Run configuration below) that sets at least `run:`, -# the campaign's name; workflow/config.yaml's machines: table supplies the rest. +# the campaign's name; workflow/config.yaml's input_types: and machines: tables supply the rest. # `psf_model` is `psfex` or `mccd`. psfex is exercised by smk-g4 through smk-g6; mccd has run the full chain on # an image-sim star tile (one focal-plane model per exposure, ~1.5 CPU-hours each). +# `tile_detection` is `unions_catalogue` (the input_types default for data: +# the UNIONS per-tile catalogue at `inputs.catalogues` is fetched and +# converted in place, keeping its NUMBER, and its segmentation map sets +# neighbours' VIGNET pixels to -1e30 as SExtractor does) or `sextractor` (the +# tile is detected with SExtractor; the default for image sims). Either way +# make_cat writes TILE_UNIQUE_ID = tile_id * 10**6 + NUMBER and ngmix masks +# the same neighbours. # The committed launcher loads apptainer/1.4.5 + the /project venv, so a # fresh shell always has the right state. @@ -77,9 +84,11 @@ A run config passed with `-c/--config-file` is merged on top of `workflow/config.yaml` and snapshotted with the code. (`-c` is `sp`'s own flag; pass snakemake's cores as `--cores`/`-j`. `SP_RUN_CONFIG` still works and is what the jobs read.) `SP_PROFILE` (default `nibi`, or `machine:` in the run config, which must -agree with it) and `input_type:` then select an entry of the `machines:` table, which supplies -`tile_list`, `retrieve` (`symlink` or `vos`), `inputs`, `outputs` and -`container` for any of these the run config leaves unset (`$base_dir` expands +agree with it) and `input_type:` then select defaults for whatever the run +config leaves unset: first the `input_types:` entry, which supplies what follows +from the kind of input (`tile_detection`, `psf_model`), then, overriding it, +the `machines:` entry, which supplies `tile_list`, `retrieve` (`symlink` or +`vos`), `inputs`, `outputs`, `container` and `psf_dict` (`$base_dir` expands to that machine's `base_dir`, and `$name` or `${name}` to any top-level scalar of the run config, e.g. `$run` to `run:`; write `${name}` when word characters follow). `run:` is @@ -92,8 +101,6 @@ SKiLLS shear branch on candide: machine: candide input_type: image_sims run: 1z2z_grid_3 -psf_model: fake -psf_dict: /home/hervas/fhervas/workdir_skills/input/psf_files/Full_psf_dict.pickle tile_list: /path/to/tiles.txt inputs: tiles: /n09data/hervas/skills_out/1z2z_grid_3/images/SP_tiles @@ -231,7 +238,7 @@ workflow/ rules/ prepare.smk tile get_images/uncompress/find_exposures exposure.smk per-exposure: get_images, split, psf, persist (no temp()); campaign star_cat_merge - tile.smk per-tile: exp forest, merge_headers, detect, vignets, ngmix, merge, make_cat; campaign final_cat_merge + tile.smk per-tile: exp forest, merge_headers, detect (SExtractor, or fetch + convert the UNIONS catalogue), vignets, ngmix, merge, make_cat; campaign final_cat_merge scripts/ build_index.py prepare-phase run_index.sqlite builder (plain script) build_forest.py per-tile exposure symlink forest (group-compatible shell) diff --git a/workflow/Snakefile b/workflow/Snakefile index 913b84543..c9124fcd2 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -50,7 +50,8 @@ from snakemake.logging import logger # --directory on /scratch (bin/sp) so .snakemake/ state never lands on /project. # # config.yaml is always loaded; SP_RUN_CONFIG (bin/sp) is merged on top, then -# the `machines:` entry for (machine, input_type) fills whatever is still unset +# the `input_types:` entry for input_type and the `machines:` entry for +# (machine, input_type) fill whatever is still unset, the machine entry winning # (scripts/run_config.py, shared with bin/sp and container.py). sys.path.insert(0, str(Path(workflow.snakefile).parent / "scripts")) import run_config # noqa: E402 @@ -81,7 +82,7 @@ if config.get("machine") is not None and config["machine"] != SP_PROFILE: f"SP_PROFILE={SP_PROFILE!r}. Set SP_PROFILE={MACHINE}, or fix " f"`machine:`.") -run_config.apply_machine_defaults(config) +run_config.apply_defaults(config) _unresolved = run_config.unresolved(config) if _unresolved: raise WorkflowError( @@ -172,7 +173,20 @@ SNAPSHOT_JSON = Path(workflow.basedir).parent / "snapshot.json" sys.path.insert(0, str(SCRIPTS)) import build_index # noqa: E402 -from completeness import STAGE_DIR # noqa: E402 +from completeness import STAGE_DIR, TILE_DETECTIONS # noqa: E402 + +# Where the tile's galaxy sample comes from: SExtractor on the tile image, or +# the UNIONS per-tile catalogue fetched and converted in place (tile.smk). +# Both write run_sp_tile_Sx, so nothing downstream of tile_detect branches. +TILE_DETECTION = config.get("tile_detection", "sextractor") +if TILE_DETECTION not in TILE_DETECTIONS: + raise WorkflowError( + f"Invalid tile_detection={TILE_DETECTION!r}; expected one of " + f"{sorted(TILE_DETECTIONS)}.") +try: + CATALOGUES, CATALOGUE_RETRIEVE = run_config.catalogue_source(config) +except ValueError as err: + raise WorkflowError(str(err)) from None # SP_PHASE is set by bin/sp and NOWHERE else: `prepare`/`compute` on the two # invocations of `sp run`, `passthrough` on direct commands (`sp --unlock`, `sp @@ -452,7 +466,7 @@ MERGE_FINAL_HASH = ":".join(( path_hash(CONFIG_DIR / "final_cat.param"))) # ngmix_range.py earns a hash for a stronger reason than the others. What it # emits is not a stale RESULT but a stale BOUNDARY, and a tile's eight chunks are -# a PARTITION of its object IDs: resume a tile across an edit to the split and +# a PARTITION of its catalogue rows: resume a tile across an edit to the split and # the chunks that already succeeded keep the old ranges while the reruns take the # new ones, so within one tile some objects are measured twice and others by # nobody. merge_sep_cats concatenates whatever it is handed, so the tile @@ -1265,7 +1279,8 @@ rule prepare_all_tiles: # emitted automatically at the end of the COMPUTE invocation, runnable any time # via `sp report`. def _report(status): - shell(f"python {SCRIPTS}/run_report.py --run-dir {RUN_DIR} " + shell(f"SP_TILE_DETECTION='{TILE_DETECTION}' " + f"python {SCRIPTS}/run_report.py --run-dir {RUN_DIR} " f"--index {INDEX_DB} --status {status} || true") if PHASE == "compute": diff --git a/workflow/bin/sp b/workflow/bin/sp index 944595172..cb6754108 100755 --- a/workflow/bin/sp +++ b/workflow/bin/sp @@ -37,7 +37,8 @@ # # Run config: -c FILE, also --config-file. Read after workflow/config.yaml and # merged on top; anything still unset comes from the machines: entry for -# SP_PROFILE and input_type (workflow/scripts/run_config.py). +# SP_PROFILE and input_type, then the input_types: entry +# (workflow/scripts/run_config.py). # To check a resolved value without running: # # python workflow/scripts/run_config.py workflow/config.yaml ~/my_run.yaml outputs.run_dir @@ -348,6 +349,7 @@ case "$cmd" in # call it through the Snakefile's SCRIPTS, which is the snapshot). Read-only # either way -- this is consistency, not safety. shift + SP_TILE_DETECTION="$(cfg tile_detection)" \ python "$(code_root)/workflow/scripts/run_report.py" \ --run-dir "$RUN_DIR" --index "$INDEX_DB" \ --status manual "$@" diff --git a/workflow/config.yaml b/workflow/config.yaml index 95a0f3d88..b6f22a01e 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -16,15 +16,27 @@ # input_tyle: allowed are data or image_sims input_type: data -# psf_model (psfex, mccd, or fake) and psf_dict are set per (machine, -# input_type) in the machines: table below, NOT here: a top-level key -# shadows that table. The Snakefile falls back to psfex if no entry -# supplies one. - # Run name, available as `$run` in the paths below; it also names the # campaign's merged catalogues. Required, and set by the run config. # run: smk-g6 +# Default settings per input_type, can be overridden by a user-defined run config file. +# - tile_detection: allowed are `unions_catalogue` (fetch official UNIONS catalogue +# from vos; recommended for data); `sextractor` (detect objects via a SExtractor +# module run; required for image_sims) +# - psf_model: allowed are `psfex` (default), `mccd`; `fake` (image_sims only, +# PSF read from psf_dict) +input_types: + data: + # @sc [decision:detection.tile_detection] + tile_detection: unions_catalogue + # @sc [decision:star_selection_psf.psf_modelling_software] + psf_model: psfex + image_sims: + tile_detection: sextractor + # @sc [decision:star_selection_psf.psf_modelling_software] + psf_model: fake + # Entries per machine (SP_PROFILE, default nibi; implemented: nibi, candide) # and input_type. A run config may instead state `machine:`, which must # then agree with SP_PROFILE. @@ -32,6 +44,9 @@ input_type: data # - tile_list: text file with tile IDs # - retrieve: method to get input data, allowed are symlink, vos # - inputs.tiles,, .exposures: path for input tile and exposure +# - inputs.catalogues: UNIONS per-tile catalogues, a local directory +# (symlinked in) or a vos: URL (tile_detection: unions_catalogue) +# - psf_dict: simulation PSF stamps (psf_model: fake) # - outputs.run_dir: (scratch) run directory, where tmp files will be stored # - outputs.products_dir: path to final products (default run_dir; one root # needs clean: false and clean_tiles: false) @@ -42,13 +57,15 @@ machines: nibi: base_dir: /project/def-mjhudson data: - # @sc [decision:star_selection_psf.psf_modelling_software] - psf_model: psfex tile_list: $base_dir/cdaley/sp-products/$run/tiles.txt retrieve: symlink inputs: tiles: $base_dir/unions-wl/tiles exposures: $base_dir/unions-wl/exposures + # UNIONS per-tile catalogues and segmentation maps (CFIS..r.cat, + # CFIS..r.seg.fits.fz), fetched by vcp in + # the job: needs network from compute nodes and ~/.ssl/cadcproxy.pem. + catalogues: vos:cfis/tiles_DR6 outputs: run_dir: /scratch/cdaley/shapepipe-output/$run products_dir: $base_dir/cdaley/sp-products/$run @@ -58,22 +75,19 @@ machines: candide: base_dir: /n17data/UNIONS/WL data: - # @sc [decision:star_selection_psf.psf_modelling_software] - psf_model: psfex retrieve: vos + inputs: + catalogues: vos:cfis/tiles_DR6 container: /n17data/cdaley/containers/shapepipe_develop-runtime-20260718.sif image_sims: retrieve: symlink - # True simulation PSF: no exposure PSF fit; fake_interp_runner - # reads psf_dict below. - # @sc [decision:star_selection_psf.psf_modelling_software] - psf_model: fake inputs: # $run is the branch dir (e.g. 1p2z_grid_1), so a run config # picks the grid by setting `run:` alone. tiles: /n09data/hervas/skills_out/$run/images/SP_tiles exposures: /n09data/hervas/skills_out/$run/images/SP_exp - # Dictionary of PSF model stamps + # Dictionary of PSF model stamps, which psf_model: fake reads in place + # of an exposure PSF fit. psf_dict: /home/hervas/fhervas/workdir_skills/input/psf_files/Full_psf_dict.pickle container: /n17data/cdaley/containers/shapepipe_im_sims-runtime.sif diff --git a/example/cfis/config_tile_Git_cat_vos.ini b/workflow/config/cfis/config_tile_Gic.ini similarity index 68% rename from example/cfis/config_tile_Git_cat_vos.ini rename to workflow/config/cfis/config_tile_Gic.ini index d2a18d816..1e6bdb7f5 100644 --- a/example/cfis/config_tile_Git_cat_vos.ini +++ b/workflow/config/cfis/config_tile_Gic.ini @@ -1,4 +1,6 @@ -# ShapePipe configuration file for: get external tile catalogue +# ShapePipe configuration file for: get the UNIONS per-tile object catalogue +# and its r-band segmentation map +# (tile_detection: unions_catalogue in the workflow run config) ## Default ShapePipe options @@ -11,7 +13,7 @@ VERBOSE = False RUN_NAME = run_sp_tile_Gic # Add date and time to RUN_NAME, optional, default: False -RUN_DATETIME = True +RUN_DATETIME = False ## ShapePipe execution options @@ -52,7 +54,7 @@ TIMEOUT = 96:00:00 ## Module options -# Get external tile catalogue +# Get the external tile catalogue and its segmentation map [GET_IMAGES_RUNNER] FILE_PATTERN = tile_numbers @@ -64,23 +66,25 @@ NUMBERING_SCHEME = # Paths -# Input path where catalogues are stored -INPUT_PATH = vos:cfis/tiles_DR6 +# Where the catalogues are: the run config's inputs.catalogues, a local +# directory or a vos: URL (vos:cfis/tiles_DR6) +INPUT_PATH = $SP_INPUT_CATALOGUES, $SP_INPUT_CATALOGUES # Input file pattern including tile number as dummy template -INPUT_FILE_PATTERN = CFIS.000.000.r +INPUT_FILE_PATTERN = CFIS.000.000.r, CFIS.000.000.r.seg # Input file extensions -INPUT_FILE_EXT = .cat +INPUT_FILE_EXT = .cat, .fits.fz # Input numbering scheme, python regexp INPUT_NUMBERING = \d{3}\.\d{3} # Output file pattern without number -OUTPUT_FILE_PATTERN = CFIS_cat- +OUTPUT_FILE_PATTERN = CFIS_cat-, CFIS_seg- -# Copy/download method, one in 'vos', 'symlink' -RETRIEVE = vos +# Copy/download method, one in 'vos', 'symlink'; the workflow derives it from +# inputs.catalogues (vos: URL -> vos, local directory -> symlink) +RETRIEVE = $SP_RETRIEVE_CATALOGUES # If RETRIEVE=vos, number of attempts to download # Optional, default=3 diff --git a/workflow/config/cfis/config_tile_Ng_template.ini b/workflow/config/cfis/config_tile_Ng_template.ini index 8ec39f5f1..61addb4f2 100644 --- a/workflow/config/cfis/config_tile_Ng_template.ini +++ b/workflow/config/cfis/config_tile_Ng_template.ini @@ -120,9 +120,8 @@ MAG_ZP = 30.0 # @sc [decision:shape_measurement.fit_priors] PIXEL_SCALE = 0.186 -# ID_OBJ_MIN/MAX: this chunk's closed SExtractor NUMBER-column range, -# computed at execution time from the tile's own object count and expanded -# via ShapePipe's getexpanded (ngmix_runner.py verified: env-expanded, not -# plain getint). -ID_OBJ_MIN = $NGMIX_ID_MIN -ID_OBJ_MAX = $NGMIX_ID_MAX +# ID_OBJ_MIN/MAX: this chunk's closed range of 1-based catalogue ROWS (not +# NUMBER values), computed at execution time by ngmix_range.py and expanded +# via ShapePipe's getexpanded (env-expanded, not plain getint). +ID_OBJ_MIN = $NGMIX_ROW_MIN +ID_OBJ_MAX = $NGMIX_ROW_MAX diff --git a/example/cfis/config_tile_Uc.ini b/workflow/config/cfis/config_tile_Uc.ini similarity index 59% rename from example/cfis/config_tile_Uc.ini rename to workflow/config/cfis/config_tile_Uc.ini index 628e81a0c..ab9d04742 100644 --- a/example/cfis/config_tile_Uc.ini +++ b/workflow/config/cfis/config_tile_Uc.ini @@ -1,5 +1,7 @@ -# ShapePipe configuration file for tile object selection using -# an (external) catalogue +# ShapePipe configuration file for tile object selection from the UNIONS +# per-tile catalogue (tile_detection: unions_catalogue in the workflow run +# config). Writes the same run dir as config_tile_Sx.ini, so that the tile +# chain downstream reads one path whichever detection mode produced it. ## Default ShapePipe options @@ -9,10 +11,10 @@ VERBOSE = True # Name of run (optional) default: shapepipe_run -RUN_NAME = run_sp_tile_Uc +RUN_NAME = run_sp_tile_Sx # Add date and time to RUN_NAME, optional, default: True -; RUN_DATETIME = False +RUN_DATETIME = False ## ShapePipe execution options @@ -20,7 +22,6 @@ RUN_NAME = run_sp_tile_Uc # Module name, single string or comma-separated list of valid module runner names MODULE = read_ext_sexcat_runner - # Run mode, SMP or MPI MODE = SMP @@ -35,6 +36,10 @@ LOG_NAME = log_sp # Runner log file name, optional, default: shapepipe_runs RUN_LOG_NAME = log_run_sp +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed tile ID (e.g. -210-282). +NUMBER_LIST = $SP_UNIT_NUM + # Input directory, containing input files, single string or list of names with length matching FILE_PATTERN INPUT_DIR = $SP_RUN/output @@ -56,11 +61,11 @@ TIMEOUT = 96:00:00 [READ_EXT_SEXCAT_RUNNER] -INPUT_DIR = run_sp_tile_Gic:get_images_runner, run_sp_tile_Git:get_images_runner, run_sp_tile_Mh_exp:merge_headers_runner +INPUT_DIR = $SP_RUN/output/run_sp_tile_Gic/get_images_runner/output, $SP_RUN/output/run_sp_tile_Git/get_images_runner/output, $SP_RUN/output/run_sp_tile_Gic/get_images_runner/output, $SP_RUN/output/run_sp_tile_Mh_exp/merge_headers_runner/output -FILE_PATTERN = CFIS_cat, CFIS_image, log_exp_headers +FILE_PATTERN = CFIS_cat, CFIS_image, CFIS_seg, log_exp_headers -FILE_EXT = .cat, .fits, .sqlite +FILE_EXT = .cat, .fits, .fitsfz, .sqlite # NUMBERING_SCHEME (optional) string with numbering pattern for input files NUMBERING_SCHEME = -000-000 @@ -72,13 +77,20 @@ SUFFIX = sexcat # image, in pixels (must be odd). Default: 51 VIGNET_SIZE = 51 +# The third input is the catalogue's r-band segmentation map: set neighbours' +# VIGNET pixels to -1e30, as SExtractor does, so ngmix masks them +# @sc [decision:detection.catalogue_neighbour_marking] +SEGMENTATION = True + ## Post-processing # Necessary for tiles, to enable multi-exposure processing +# @sc [decision:detection.epoch_membership_ccd_bounds] MAKE_POST_PROCESS = True # World coordinate keywords, SExtractor output. Format: KEY_X,KEY_Y WORLD_POSITION = ALPHA_J2000,DELTA_J2000 # Number of pixels in x,y of a CCD. Format: Nx,Ny +# @sc [decision:detection.epoch_membership_ccd_bounds] CCD_SIZE = 33,2080,1,4612 diff --git a/workflow/config/cfis/final_cat.param b/workflow/config/cfis/final_cat.param index c3719fc54..cc5feaabf 100644 --- a/workflow/config/cfis/final_cat.param +++ b/workflow/config/cfis/final_cat.param @@ -8,10 +8,14 @@ YWIN_WORLD # @sc [decision:catalogue_assembly.tile_overlap_handling] TILE_ID -# SExtractor object number within the tile: sp_validation's extraction joins -# on it, and the per-tile catalogue always carries it. +# Object number within the tile: sp_validation's extraction joins on it, and +# the per-tile catalogue always carries it. NUMBER +# Survey-wide object ID, tile_id * 10**6 + NUMBER (make_cat, both +# tile_detection modes): the join key to other catalogues. +TILE_UNIQUE_ID + # flags FLAGS # NO IMAFLAGS_ISO, AND NO MASK COLUMN AT ALL — READ THIS BEFORE ADDING ONE. @@ -222,6 +226,10 @@ NGMIX_FLUX_ERR_2P NGMIX_FLUX_ERR_NOSHEAR # magnitudes +# SExtractor detection columns. Under tile_detection: unions_catalogue the +# UNIONS catalogue supplies them, and the merge drops the ones it lacks +# (merge_final_cat.SEXTRACTOR_ONLY_COLUMNS: the WIN, AUTO-flux, APER and FWHM +# columns and SNR_WIN); MAG_AUTO, MAGERR_AUTO and FLUX_RADIUS are in both. MAG_AUTO MAGERR_AUTO MAG_WIN diff --git a/workflow/config/cfis_image_sims/final_cat.param b/workflow/config/cfis_image_sims/final_cat.param index 40c7c9faa..9db095b66 100644 --- a/workflow/config/cfis_image_sims/final_cat.param +++ b/workflow/config/cfis_image_sims/final_cat.param @@ -19,6 +19,9 @@ YWIN_WORLD # tile ID, for plot of tile-dependent additive bias. TILE_ID +# Survey-wide object ID, tile_id * 10**6 + NUMBER (make_cat). +TILE_UNIQUE_ID + # SExtractor number NUMBER diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index 8f2e21f24..f105d23b4 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -41,6 +41,19 @@ The heavy middle (tile_detect) stays out: it is a 16 GB / 8 thread SExtractor run that the shape chain does not need co-scheduled, and folding it in would add its runtime to a sum that has no room. +``tile_detection: unions_catalogue`` (config.yaml) replaces that SExtractor run +with two rules: tile_get_catalogue fetches the UNIONS per-tile catalogue and +its r-band segmentation map (get_images_runner, config_tile_Gic.ini) and +tile_detect converts it to the FITS-LDAC sexcat SExtractor would have written +(read_ext_sexcat_runner, config_tile_Uc.ini), with the tile image's header, +VIGNET stamps cut from the tile image with neighbours' footprints set to -1e30 +as SExtractor sets them (so ngmix masks neighbours in both modes), and the +multi-epoch post-processing. It keeps the catalogue's own +NUMBER, from which make_cat builds ``TILE_UNIQUE_ID`` exactly as in SExtractor +mode. The converter writes run_sp_tile_Sx/read_ext_sexcat_runner, and the rule links it as +sextractor_runner, the one path every downstream config reads; the manifest is +tile_detect.json in both modes, so tile_vignets onwards is the same DAG. + There is no `tile_mask` rule, and there will not be one (PR #847). ShapePipe generates no masks: tiles have no instrument flag image of their own, so tile_detect runs SExtractor with FLAG_IMAGE = False against @@ -459,24 +472,87 @@ rule tile_merge_headers: shell: sp_shell("tile_merge_headers", "config_tile_Mh_exp.ini") -# SExtractor object detection on the tile. -rule tile_detect: - input: - uz = f"{TILE_DIR}/manifests/tile_uncompress.json", - mh = rules.tile_merge_headers.output.manifest, - output: - manifest = f"{TILE_DIR}/manifests/tile_detect.json" - log: - f"{TILE_DIR}/logs/tile_detect.json" - params: - pre = lambda wc: unit_pre("tile_detect", wc.tile), - script_hash = SCRIPT_HASH - threads: 8 - resources: - mem_mb = lambda wc, attempt: 16000 * attempt, - runtime = 180 - shell: - sp_shell("tile_detect", "config_tile_Sx.ini") +# The tile's galaxy sample. One of two definitions of tile_detect, chosen at +# parse time by the run config; both produce manifests/tile_detect.json and +# run_sp_tile_Sx/sextractor_runner/output/sexcat.fits. +if TILE_DETECTION == "sextractor": + + # SExtractor object detection on the tile. + rule tile_detect: + input: + uz = f"{TILE_DIR}/manifests/tile_uncompress.json", + mh = rules.tile_merge_headers.output.manifest, + output: + manifest = f"{TILE_DIR}/manifests/tile_detect.json" + log: + f"{TILE_DIR}/logs/tile_detect.json" + params: + pre = lambda wc: unit_pre("tile_detect", wc.tile), + script_hash = SCRIPT_HASH + threads: 8 + resources: + mem_mb = lambda wc, attempt: 16000 * attempt, + runtime = 180 + shell: + sp_shell("tile_detect", "config_tile_Sx.ini") + +else: + + # Fetch the UNIONS per-tile catalogue (CFIS..r.cat) and segmentation + # map (CFIS..r.seg.fits.fz), from a local mirror or from vos, the way tile_get_images fetches the image. Reads + # only tile_numbers.txt, which unit_pre writes; the edge on the image + # manifest is what puts it after the prepare phase. + rule tile_get_catalogue: + input: + git = f"{TILE_DIR}/manifests/tile_get_images.json", + output: + manifest = f"{TILE_DIR}/manifests/tile_get_catalogue.json" + log: + f"{TILE_DIR}/logs/tile_get_catalogue.json" + params: + pre = lambda wc: unit_pre( + "tile_get_catalogue", wc.tile, + env={"SP_INPUT_CATALOGUES": CATALOGUES, + "SP_RETRIEVE_CATALOGUES": CATALOGUE_RETRIEVE}), + script_hash = SCRIPT_HASH + threads: 1 + retries: 2 + resources: + mem_mb = lambda wc, attempt: 4000 * attempt, + runtime = 60 + shell: + sp_shell("tile_get_catalogue", "config_tile_Gic.ini") + + # Convert the catalogue to the FITS-LDAC sexcat the chain expects. The + # completeness check counts read_ext_sexcat_runner's own output; the link + # after it is what makes that output readable at the sextractor_runner path + # of every downstream config, and it is made only on success so a failed + # run leaves nothing at that path. + rule tile_detect: + input: + cat = rules.tile_get_catalogue.output.manifest, + git = f"{TILE_DIR}/manifests/tile_get_images.json", + mh = rules.tile_merge_headers.output.manifest, + output: + manifest = f"{TILE_DIR}/manifests/tile_detect.json" + log: + f"{TILE_DIR}/logs/tile_detect.json" + params: + pre = lambda wc: unit_pre( + "tile_detect", wc.tile, + env={"SP_TILE_DETECTION": TILE_DETECTION}), + script_hash = SCRIPT_HASH + threads: 1 + resources: + # The whole tile image plus one 51x51 float32 stamp per object. + mem_mb = lambda wc, attempt: 8000 * attempt, + runtime = 60 + shell: + sp_shell("tile_detect", "config_tile_Uc.ini", + post="if [ $rc -eq 0 ]; then\n" + ' ln -s read_ext_sexcat_runner ' + '"$SP_RUN/output/run_sp_tile_Sx/sextractor_runner"\n' + "fi\n") # Configured PSF interpolation to galaxies + vignet postage stamps: the last # stage that reads exposure products, and the bulk intra-tile intermediate. The store it @@ -535,7 +611,7 @@ rule tile_vignets: check_args=' --run-dir "$SP_LOCAL" --unit {wildcards.tile}') # ngmix shape measurement — N chunks per tile (D4). Each chunk LOOKS UP its own -# CLOSED object-ID range in the file tile_vignets materialised at the top of this +# CLOSED catalogue-row range in the file tile_vignets materialised at the top of this # group job (TILE_NGMIX_RANGES); the ranges are knowable only at EXECUTION time, # from this tile's own sexcat, which is why a params function cannot supply them. # Closed, not open-ended: `ID_OBJ_MAX = -1` on the last chunk was the 13-hour @@ -964,6 +1040,7 @@ rule final_cat_merge: tile_list = str(config["tile_list"]), index_db = str(INDEX_DB), param_file = str(CONFIG_DIR / "final_cat.param"), + tile_detection = TILE_DETECTION, campaign = CAMPAIGN, snapshot = str(SNAPSHOT_JSON), inputs = unit_fingerprint(TILES_READY), @@ -988,4 +1065,5 @@ rule final_cat_merge: " --output {output.merged}" " --campaign '{params.campaign}'" " --param-file '{params.param_file}'" + " --tile-detection {params.tile_detection}" " --snapshot-json '{params.snapshot}'" diff --git a/workflow/scripts/completeness.py b/workflow/scripts/completeness.py index 89cd63601..52720b724 100644 --- a/workflow/scripts/completeness.py +++ b/workflow/scripts/completeness.py @@ -72,8 +72,15 @@ import sys from pathlib import Path +# How the tile's galaxy sample is produced: SExtractor on the tile image, or +# the UNIONS per-tile catalogue converted in place. The run config's +# `tile_detection:` picks one; the rules export it as $SP_TILE_DETECTION only +# when it is not the default, so a SExtractor run's prologue is unchanged. +TILE_DETECTIONS = ("sextractor", "unions_catalogue") + # stage -> {runner_subdir: {expect, [warn], [subpath]}} -# exp_psf and tile_vignets are selected by $SP_PSF at check time. +# exp_psf and tile_vignets are selected by $SP_PSF at check time, tile_detect +# by $SP_TILE_DETECTION. # @sc [decision:per_unit_completeness] COMPLETENESS = { # --- tile prepare (phase A) --- @@ -135,7 +142,13 @@ # --- tile post --- "tile_merge_headers": {"merge_headers_runner": dict(expect=1)}, - "tile_detect": {"sextractor_runner": dict(expect=2)}, + # The fetched UNIONS catalogue and its r-band segmentation map. + "tile_get_catalogue": {"get_images_runner": dict(expect=2)}, + "tile_detect": { + "sextractor": {"sextractor_runner": dict(expect=2)}, + # The FITS-LDAC sexcat converted from the fetched catalogue. + "unions_catalogue": {"read_ext_sexcat_runner": dict(expect=1)}, + }, "tile_vignets": { "psfex": { "psfex_interp_runner": dict(expect=1), @@ -210,6 +223,13 @@ def check_counts(stage, run_dir): raise ValueError( f"Invalid SP_PSF={psf_model!r}; expected one of {sorted(table)}." ) from exc + elif stage == "tile_detect": + detection = os.environ.get("SP_TILE_DETECTION", TILE_DETECTIONS[0]) + if detection not in TILE_DETECTIONS: + raise ValueError( + f"Invalid SP_TILE_DETECTION={detection!r}; expected one of " + f"{', '.join(TILE_DETECTIONS)}.") + table = table[detection] details, ok = [], True for runner, spec in table.items(): n = count_products(run_dir, runner, spec) @@ -241,6 +261,9 @@ def check_counts(stage, run_dir): "exp_split": ("exp", "run_sp_exp_Sp"), "exp_psf": ("exp", "run_sp_exp_SxSePsf"), "tile_merge_headers": ("tile", "run_sp_tile_Mh_exp"), + "tile_get_catalogue": ("tile", "run_sp_tile_Gic"), + # Both detection modes write here (config_tile_Sx.ini / config_tile_Uc.ini + # share the RUN_NAME): the chain downstream reads one path. "tile_detect": ("tile", "run_sp_tile_Sx"), "tile_vignets": ("tile", "run_sp_tile_PiViVi"), "tile_ngmix": ("tile", "run_sp_tile_ngmix_Ng${SP_NGMIX_CHUNK}u"), diff --git a/workflow/scripts/merge_final_cat.py b/workflow/scripts/merge_final_cat.py index 1ad6fbdff..f4be98e06 100644 --- a/workflow/scripts/merge_final_cat.py +++ b/workflow/scripts/merge_final_cat.py @@ -52,6 +52,18 @@ argument validator and then falls through to the ordinary walk, so it is not a way to add one tile by hand.) +WHICH COLUMNS DEPEND ON WHERE THE TILE'S DETECTIONS CAME FROM. The param file +lists the columns of a SExtractor-mode tile catalogue. Under ``tile_detection: +unions_catalogue`` the detection columns are instead the UNIONS catalogue's own +(read_ext_sexcat copies them as they are), and that catalogue lacks some of the +SExtractor quantities the param file names; ``SEXTRACTOR_ONLY_COLUMNS`` lists +them, and ``requested_columns`` subtracts exactly that list in that mode. The +merge stays strict on everything else: a column still requested and absent +from a tile stops the merge. tests/module/test_final_cat_columns.py derives the +catalogue-mode columns by running the converter on the UNIONS catalogue's +header and asserts the list is exactly the requested SExtractor columns the +catalogue does not carry, so it cannot silently go stale in either direction. + WHICH TILES — AND WHY THE JOB DERIVES THE SET RATHER THAN BEING TOLD IT. The set is the CAMPAIGN's: every tile both declared in ``tile_list`` and present in the index, which is exactly the Snakefile's TILES_READY, rebuilt here from the same @@ -88,12 +100,41 @@ # Same directory; the rule invokes this file by path, so it is sys.path[0]. import build_index import hdf5_reconcile +from completeness import TILE_DETECTIONS # /scripts/python/create_final_cat.py, from /workflow/scripts/this. CFC_PATH = (Path(__file__).resolve().parents[2] / "scripts" / "python" / "create_final_cat.py") +# Requested SExtractor columns the UNIONS per-tile catalogue (DR6) does not +# carry, so a unions_catalogue-mode tile catalogue cannot have them. None has a +# DR6 column measuring the same quantity: DR6's MAG_APER is a magnitude through +# its own aperture, not FLUX_APER; MAG_AUTO, MAGERR_AUTO and FLUX_RADIUS are +# the same SExtractor measurements in both modes and stay requested. +SEXTRACTOR_ONLY_COLUMNS = ( + "MAG_WIN", "MAGERR_WIN", # Gaussian-windowed magnitude + "FLUX_AUTO", "FLUXERR_AUTO", # DR6 carries the AUTO magnitude only + "FLUX_APER", "FLUXERR_APER", # ShapePipe's aperture, in flux + "SNR_WIN", # Gaussian-windowed SNR + "FWHM_IMAGE", "FWHM_WORLD", # Gaussian-core FWHM +) + + +def requested_columns(param_list: list, tile_detection: str) -> list: + """The columns a tile catalogue of this detection mode must carry. + + The param file's list, in its order, less ``SEXTRACTOR_ONLY_COLUMNS`` when + the detections are the UNIONS catalogue's (see the module docstring). + """ + if tile_detection not in TILE_DETECTIONS: + raise ValueError(f"tile_detection {tile_detection!r} is not one of " + f"{TILE_DETECTIONS}") + if tile_detection == "sextractor": + return list(param_list) + return [c for c in param_list if c not in SEXTRACTOR_ONLY_COLUMNS] + + def spval_group(campaign: str) -> str: """The hdf5 group the campaign's per-tile datasets live under. @@ -152,6 +193,9 @@ def main() -> None: help="names the campaign's group in the output file") p.add_argument("--param-file", required=True, type=Path, help="the input type's final_cat.param — the column list") + p.add_argument("--tile-detection", required=True, choices=TILE_DETECTIONS, + help="where the tile detections came from (config " + "tile_detection); selects the columns requested") p.add_argument("--hdu", type=int, default=1) p.add_argument("--snapshot-json", type=Path, default=None, help="sp run's code snapshot (bin/sp's " @@ -159,7 +203,9 @@ def main() -> None: args = p.parse_args() cfc = load_create_final_cat() - param_list = cfc.read_param_file(str(args.param_file), verbose=False) + param_list = requested_columns( + cfc.read_param_file(str(args.param_file), verbose=False), + args.tile_detection) if not param_list: sys.exit(f"merge_final_cat: no columns read from {args.param_file}") # read_data/copy_data read their knobs out of this dict, exactly as diff --git a/workflow/scripts/ngmix_range.py b/workflow/scripts/ngmix_range.py index aa745c2ce..6f2cbbf03 100644 --- a/workflow/scripts/ngmix_range.py +++ b/workflow/scripts/ngmix_range.py @@ -11,7 +11,7 @@ # each chunk, at its own start eval "$(ngmix_range.py --read $SP_LOCAL/ngmix_ranges.json --chunk 3)" - # -> export NGMIX_ID_MIN=751; export NGMIX_ID_MAX=1125 + # -> export NGMIX_ROW_MIN=751; export NGMIX_ROW_MAX=1125 The file is group-internal plumbing, NOT DAG currency: it lives on $SP_LOCAL and dies with the group job, exactly like the vignette store. It is deliberately not @@ -24,10 +24,12 @@ The range is still only knowable at EXECUTION time (PRD D4) — a params function cannot compute it, because params evaluate before the sexcat exists. -SExtractor's NUMBER column (ngmix's obj_id) runs 1..N contiguous, so covering -[1, N] processes every object exactly once. The bounds are CLOSED, never -ID_OBJ_MAX = -1 — ngmix treats ``id_obj_max <= 0`` as unbounded -(``ngmix_package/ngmix.py:804-806``), so an open-ended last chunk silently +The ranges are 1-based ROW positions in the sexcat, not NUMBER values: ngmix +selects its chunk by row (``ngmix_package.ngmix.chunk_rows``), so covering +[1, N] processes every object exactly once however NUMBER is ordered or +spaced (an external detection catalogue keeps its own NUMBER). The bounds are +CLOSED, never ID_OBJ_MAX = -1 — ngmix treats ``id_obj_max <= 0`` as +unbounded, so an open-ended last chunk silently re-measures the whole tile instead of its share; rule tile_ngmix carries the straggler that taught this. @@ -51,7 +53,7 @@ measurement (see ``ngmix_package.ngmix.position_seed``). DETERMINISM HERE IS CORRECTNESS, NOT TIDINESS. A tile's chunk ranges are a -PARTITION of its object IDs: if two chunks disagree about the boundaries, +PARTITION of its catalogue rows: if two chunks disagree about the boundaries, objects are silently measured twice or silently dropped and nothing downstream notices — merge_sep_cats concatenates whatever it is given. @@ -103,11 +105,11 @@ ALPHA_MILLI_EPOCHS = 184 -def id_ranges(epochs, n_chunks: int) -> list[tuple[int, int]]: - """Split object IDs ``1..len(epochs)`` into ``n_chunks`` closed ranges. +def row_ranges(epochs, n_chunks: int) -> list[tuple[int, int]]: + """Split catalogue rows ``1..len(epochs)`` into ``n_chunks`` closed ranges. - ``epochs[i]`` is object ``i + 1``'s geometric epoch count. The ranges are - CONTIGUOUS — ``NGMIX_ID_MIN``/``NGMIX_ID_MAX`` is an interval, not a set — + ``epochs[i]`` is row ``i + 1``'s geometric epoch count. The ranges are + CONTIGUOUS — ``NGMIX_ROW_MIN``/``NGMIX_ROW_MAX`` is an interval, not a set — and tile ``[1, n_obj]`` exactly, so every object is measured once. The objective is the slowest chunk, not the average one, because the @@ -188,8 +190,9 @@ def object_epochs(run_dir: Path): tile_detect's SExtractor post-process (``MAKE_POST_PROCESS`` in ``config_tile_Sx.ini``) writes one ``EPOCH_`` extension per exposure overlapping the tile — so the extension COUNT is tile-specific and is - discovered by name, never assumed — each with ``n_obj`` rows in NUMBER - order and ``CCD_N < 0`` where the object misses that exposure. Summing + discovered by name, never assumed — each with ``n_obj`` rows in + ``LDAC_OBJECTS`` row order and ``CCD_N < 0`` where the object misses that + exposure. Summing ``CCD_N >= 0`` across them reproduces the final catalogue's ``N_EPOCH`` column exactly — checked row by row against 186.307's ``run_sp_tile_Mc/.../final_cat-186-307.fits``, all 35,298 of them, 7 extensions, @@ -234,19 +237,20 @@ def object_epochs(run_dir: Path): f"[ngmix_range] FATAL: no LDAC_OBJECTS in {cats[0]}" ) n_obj = int(hdul["LDAC_OBJECTS"].header["NAXIS2"]) - expected = np.arange(1, n_obj + 1, dtype=np.int64) counts = np.zeros(n_obj, dtype=np.int64) + number = None for hdu in epoch_hdus: data = hdu.data - # NUMBER is asserted, not assumed: the whole scheme is an ID - # INTERVAL, so a permuted or gappy NUMBER column would make the - # weights describe different objects than the bounds select. - if len(data) != n_obj or not np.array_equal( - np.asarray(data["NUMBER"], dtype=np.int64), expected - ): + # Row alignment is asserted, not assumed: the weights are indexed + # by row, so every EPOCH extension must carry the same NUMBER + # sequence, row for row, as the first one. + this_number = np.asarray(data["NUMBER"], dtype=np.int64) + if number is None: + number = this_number + if len(data) != n_obj or not np.array_equal(this_number, number): raise SystemExit( f"[ngmix_range] FATAL: {cats[0]}[{hdu.name}] is not " - f"{n_obj} rows of NUMBER = 1..{n_obj} in order" + f"{n_obj} rows aligned with {epoch_hdus[0].name}" ) counts += np.asarray(data["CCD_N"]) >= 0 return counts @@ -261,12 +265,12 @@ def partition(epochs, n_chunks: int) -> dict: against what was actually written, not against what the caller believes). ``chunk`` is 1-based, matching SP_NGMIX_CHUNK and the run-directory suffix. """ - ranges = id_ranges(epochs, n_chunks) + ranges = row_ranges(epochs, n_chunks) return { "n_obj": len(epochs), "n_chunks": n_chunks, "chunks": [ - {"chunk": k, "id_min": lo, "id_max": hi} + {"chunk": k, "row_min": lo, "row_max": hi} for k, (lo, hi) in enumerate(ranges, start=1) ], } @@ -289,7 +293,7 @@ def chunk_range(doc: dict, chunk: int) -> tuple[int, int]: f"[ngmix_range] FATAL: ranges file row {chunk - 1} is chunk " f"{row['chunk']}, not {chunk}" ) - return int(row["id_min"]), int(row["id_max"]) + return int(row["row_min"]), int(row["row_max"]) def read_ranges(path: Path) -> dict: @@ -342,7 +346,7 @@ def main() -> None: if a.chunk is None: raise SystemExit("[ngmix_range] FATAL: --read needs --chunk") lo, hi = chunk_range(read_ranges(a.read), a.chunk) - print(f"export NGMIX_ID_MIN={lo}; export NGMIX_ID_MAX={hi}") + print(f"export NGMIX_ROW_MIN={lo}; export NGMIX_ROW_MAX={hi}") if __name__ == "__main__": diff --git a/workflow/scripts/run_config.py b/workflow/scripts/run_config.py index 831147986..669a3ec60 100644 --- a/workflow/scripts/run_config.py +++ b/workflow/scripts/run_config.py @@ -1,7 +1,14 @@ #!/usr/bin/env python3 -"""Resolve a run config: config.yaml, then SP_RUN_CONFIG merged on top, then -the `machines:` entry for (machine, input_type) filling whatever is still -unset. One definition shared by the Snakefile, bin/sp and container.py. +"""Resolve a run config: config.yaml with SP_RUN_CONFIG merged on top, then +two tables of defaults filling whatever is still unset -- `input_types:` +entry for input_type, overridden by the `machines:` entry for +(machine, input_type). One definition shared by the Snakefile, bin/sp and +container.py. + +Precedence, lowest first: config.yaml's top level < input_types[input_type] +< machines[machine][input_type] < the run config. config.yaml's top level +sits beneath both tables only because it sets none of the keys they carry; +tests/unit/test_run_config.py holds it to that. CLI (used by bin/sp): run_config.py CONFIG_YAML RUN_CONFIG KEY[.SUBKEY] prints the resolved value, or an empty line if unset. RUN_CONFIG may be "". @@ -14,15 +21,8 @@ import yaml PLACEHOLDER = "TBD" -# Keys the machines: table may default, per (machine, input_type). -# psf_model/psf_dict belong here because they are per-input_type facts, -# not per-run ones: psf_model=fake is only legal with -# input_type=image_sims, and psf_dict is the sim PSF it reads. A key -# also present at the TOP level of config.yaml shadows the table (the -# setdefault below only fires when the key is absent), so a key listed -# here must not carry a top-level default as well. -MACHINE_KEYS = ("tile_list", "retrieve", "container", "inputs", "outputs", - "psf_model", "psf_dict") +# The tables themselves: read here, never defaults or $-expanded. +TABLES = ("input_types", "machines") # `run` is required in its own right, not only through `$run` in the paths: # it names the campaign's merged catalogues, so a run config whose paths # never mention `$run` must still set it. @@ -67,19 +67,25 @@ def _expand(value, variables): return value -def apply_machine_defaults(config): - """Fill unset MACHINE_KEYS from machines[machine][input_type], in place. +def apply_defaults(config): + """Fill unset keys from the input_types: and machines: tables, in place. The machine is `machine:` when the run config states one, else SP_PROFILE (default nibi) -- the same value bin/sp picks the SLURM profile with. - - A key already in `config` wins; for `inputs`/`outputs` the merge is per - sub-key. In all of these, `$base_dir` expands to machines[machine].base_dir - and `$name` / `${name}` to any top-level scalar (`$run` to `run:`). + input_types[input_type] holds what follows from the kind of input alone + (the PSF model, the tile detection); machines[machine][input_type] holds + where things are on that machine, and wins where the two overlap. + + A key already in `config` wins; for dict values (`inputs`, `outputs`) the + merge is per sub-key. Every resolved key then has `$name` / `${name}` + expanded: `$base_dir` is machines[machine].base_dir, and any top-level + scalar is a variable too (`$run` is the top-level `run:`). """ + input_type = config.get("input_type", "data") machine = config.get("machine") or os.environ.get("SP_PROFILE", "nibi") entry = (config.get("machines") or {}).get(machine) or {} - defaults = entry.get(config.get("input_type", "data")) or {} + defaults = merge((config.get("input_types") or {}).get(input_type) or {}, + entry.get(input_type) or {}) # Every top-level scalar is a variable, so a run config can define its own # shorthands. They are resolved AGAINST EACH OTHER first, to a fixpoint, so # one shorthand may be written in terms of another @@ -97,17 +103,13 @@ def apply_machine_defaults(config): if resolved == variables: break variables = resolved - # `run` is consumed downstream (paths, the hdf5 group name), so the - # resolved value has to go back into the config, not just the table. - if isinstance(config.get("run"), str): - config["run"] = variables.get("run", config["run"]) - for key in MACHINE_KEYS: - default = defaults.get(key) + for key, default in defaults.items(): if isinstance(default, dict): config[key] = merge(default, config.get(key) or {}) - elif default is not None: + else: config.setdefault(key, default) - if key in config: + for key in config: + if key not in TABLES: config[key] = _expand(config[key], variables) return config @@ -129,23 +131,41 @@ def _dollar_keys(value, prefix): def unresolved(config): """REQUIRED keys that are unset, the placeholder, or hold an unexpanded - $variable, then every MACHINE_KEYS value (recursively through - inputs/outputs) that still holds one (e.g. `$run` with no `run:` set, or - a misspelt name).""" + $variable, then every other resolved value (recursively through + inputs/outputs; the tables themselves excluded) that still holds one + (e.g. `$run` with no `run:` set, or a misspelt name).""" missing = [k for k in REQUIRED if get(config, k) in (None, "", PLACEHOLDER) or "$" in str(get(config, k))] - dollar = [k for key in MACHINE_KEYS for k in _dollar_keys(config.get(key), key)] + dollar = [k for key in config if key not in TABLES + for k in _dollar_keys(config[key], key)] return missing + [k for k in dollar if k not in missing] +def catalogue_source(config): + """`inputs.catalogues` and the retrieve mode its prefix implies. + + The source is a local directory (symlinked in) or a vos: URL (downloaded). + Returns ("", "symlink") when unset or the placeholder; raises ValueError + if `tile_detection: unions_catalogue` then has nothing to fetch. + """ + source = get(config, "inputs.catalogues") or "" + if source == PLACEHOLDER: + source = "" + if config.get("tile_detection") == "unions_catalogue" and not source: + raise ValueError( + "tile_detection=unions_catalogue needs `inputs.catalogues:` (a local " + "directory or vos: URL holding the per-tile CFIS..r.cat files).") + return source, "vos" if source.startswith("vos:") else "symlink" + + def load(config_yaml, run_config=None): with open(config_yaml) as f: config = yaml.safe_load(f) or {} if run_config: with open(run_config) as f: config = merge(config, yaml.safe_load(f) or {}) - return apply_machine_defaults(config) + return apply_defaults(config) if __name__ == "__main__": diff --git a/workflow/scripts/run_report.py b/workflow/scripts/run_report.py index 11cdcffe4..a8efe82ae 100644 --- a/workflow/scripts/run_report.py +++ b/workflow/scripts/run_report.py @@ -49,6 +49,7 @@ import argparse import json +import os import sqlite3 import sys from collections import defaultdict @@ -59,6 +60,11 @@ TILE_STAGES = ["tile_get_images", "tile_uncompress", "tile_find_exposures", "tile_merge_headers", "tile_detect", "tile_vignets", "tile_ngmix", "tile_merge_cats", "tile_make_cat"] +# The catalogue fetch exists only when the run converts the UNIONS catalogue +# (config.yaml's tile_detection); the callers pass the mode in the environment +# so a SExtractor run does not report the stage as not run. +if os.environ.get("SP_TILE_DETECTION") == "unions_catalogue": + TILE_STAGES.insert(TILE_STAGES.index("tile_detect"), "tile_get_catalogue") EXP_STAGES = ["exp_get_images", "exp_split", "exp_psf"] # exp_persist is DELIBERATELY NOT in that list. This report disk-scans the