From 598a87db31e412a4e4157b0f56be325d31a53b2a Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 17 Sep 2026 13:11:04 +0200 Subject: [PATCH 01/21] Add the UNIONS tile-catalogue inis to the workflow config chain Move the two v2.0-era configs (fetch Steven Gwyn's per-tile .cat, convert it with read_ext_sexcat_runner) into workflow/config/cfis and adapt them to the workflow's conventions: fixed run names, NUMBER_LIST from SP_UNIT_NUM, explicit chained INPUT_DIRs, and the catalogue source and retrieve mode taken from the environment the rules export. The converter writes into run_sp_tile_Sx so that the downstream tile chain reads one path in either detection mode. get_images_runner expands RETRIEVE so it can be set per run. Co-Authored-By: Claude Fable 5.1 --- src/shapepipe/modules/get_images_runner.py | 4 ++-- .../config/cfis/config_tile_Gic.ini | 17 ++++++++++------- .../config}/cfis/config_tile_Uc.ini | 17 +++++++++++------ 3 files changed, 23 insertions(+), 15 deletions(-) rename example/cfis/config_tile_Git_cat_vos.ini => workflow/config/cfis/config_tile_Gic.ini (76%) rename {example => workflow/config}/cfis/config_tile_Uc.ini (73%) diff --git a/src/shapepipe/modules/get_images_runner.py b/src/shapepipe/modules/get_images_runner.py index b1a03b15a..f38a9b708 100644 --- a/src/shapepipe/modules/get_images_runner.py +++ b/src/shapepipe/modules/get_images_runner.py @@ -27,8 +27,8 @@ def get_images_runner( """Define The Get Images Runner.""" # Read config file section - # Copy/download method - retrieve_method = config.get(module_config_sec, "RETRIEVE") + # Copy/download method; may be an environment variable + retrieve_method = config.getexpanded(module_config_sec, "RETRIEVE") retrieve_ok = ["vos", "symlink"] if retrieve_method not in retrieve_ok: raise ValueError( diff --git a/example/cfis/config_tile_Git_cat_vos.ini b/workflow/config/cfis/config_tile_Gic.ini similarity index 76% rename from example/cfis/config_tile_Git_cat_vos.ini rename to workflow/config/cfis/config_tile_Gic.ini index d2a18d816..ab38f38fe 100644 --- a/example/cfis/config_tile_Git_cat_vos.ini +++ b/workflow/config/cfis/config_tile_Gic.ini @@ -1,4 +1,5 @@ -# ShapePipe configuration file for: get external tile catalogue +# ShapePipe configuration file for: get the UNIONS per-tile object catalogue +# (tile_detection: unions_catalogue in the workflow run config) ## Default ShapePipe options @@ -11,7 +12,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 +53,7 @@ TIMEOUT = 96:00:00 ## Module options -# Get external tile catalogue +# Get the external tile catalogue [GET_IMAGES_RUNNER] FILE_PATTERN = tile_numbers @@ -64,8 +65,9 @@ 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 # Input file pattern including tile number as dummy template INPUT_FILE_PATTERN = CFIS.000.000.r @@ -79,8 +81,9 @@ INPUT_NUMBERING = \d{3}\.\d{3} # Output file pattern without number OUTPUT_FILE_PATTERN = CFIS_cat- -# 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/example/cfis/config_tile_Uc.ini b/workflow/config/cfis/config_tile_Uc.ini similarity index 73% rename from example/cfis/config_tile_Uc.ini rename to workflow/config/cfis/config_tile_Uc.ini index 628e81a0c..a13def13f 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,7 +61,7 @@ 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_Mh_exp/merge_headers_runner/output FILE_PATTERN = CFIS_cat, CFIS_image, log_exp_headers From 8702483d8b9e9680dc6ec9fabb2d90f97eee2a69 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 17 Sep 2026 13:15:56 +0200 Subject: [PATCH 02/21] Add tile_detection: unions_catalogue to the workflow A run config key, validated at parse time, that picks where the tile's galaxy sample comes from. `sextractor` is the existing SExtractor rule, byte for byte. `unions_catalogue` replaces it with two rules: tile_get_catalogue fetches the UNIONS per-tile .cat from inputs.catalogues (a local mirror or a vos: URL), and tile_detect converts it with read_ext_sexcat_runner into run_sp_tile_Sx, linking the module output as sextractor_runner so that tile_vignets, ngmix and make_cat read the one path they always have; the manifest is tile_detect.json in both modes. completeness.py checks tile_detect per mode through SP_TILE_DETECTION, exported on that rule alone so a SExtractor run's prologue is unchanged; run_report lists the fetch stage only for catalogue runs. The catalogue's per-tile unique object ID reaches the final catalogue as TILE_UNIQUE_ID through make_cat's copy of every sexcat column. Co-Authored-By: Claude Fable 5.1 --- workflow/README.md | 6 +- workflow/Snakefile | 22 +++++- workflow/bin/sp | 3 +- workflow/config.yaml | 11 +++ workflow/config/cfis/final_cat.param | 5 ++ workflow/rules/tile.smk | 113 ++++++++++++++++++++++----- workflow/scripts/completeness.py | 27 ++++++- workflow/scripts/run_report.py | 6 ++ 8 files changed, 169 insertions(+), 24 deletions(-) diff --git a/workflow/README.md b/workflow/README.md index 55611cedc..aff3cee0a 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -26,6 +26,10 @@ uv pip install 'snakemake>=9,<10' 'snakemake-executor-plugin-slurm>=2.7,<3' # outputs.products_dir/index_db, and container. # `psf_model` is `psfex` or `mccd`; mccd is wired but unvalidated here, while psfex is exercised by smk-g4 through smk-g6. +# `tile_detection` is `sextractor` (the tile is detected with SExtractor) or +# `unions_catalogue` (the UNIONS per-tile catalogue at `inputs.catalogues` is +# fetched and converted in place, and its object ID rides into the final +# catalogue as TILE_UNIQUE_ID; no segmentation map, so not with ngmix's uberseg). # The committed launcher loads apptainer/1.4.5 + the /project venv, so a # fresh shell always has the right state. @@ -157,7 +161,7 @@ workflow/ rules/ prepare.smk tile get_images/uncompress/find_exposures exposure.smk per-exposure: get_images, split, psf (no temp()) - tile.smk per-tile: exp forest, merge_headers, detect, vignets, ngmix, merge, make_cat + tile.smk per-tile: exp forest, merge_headers, detect (SExtractor, or fetch + convert the UNIONS catalogue), vignets, ngmix, merge, make_cat scripts/ sp_rule.py the thin per-unit wrapper (isolation furniture, config copy, log-sync, count check) build_index.py prepare-phase run_index.sqlite builder (plain script) diff --git a/workflow/Snakefile b/workflow/Snakefile index bbdfeefb6..d7091451e 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -100,7 +100,24 @@ CONFIG_DIR = Path(workflow.basedir) / "config" / "cfis" 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)}.") +# The catalogue source: a local directory (symlinked in) or a vos: URL +# (downloaded); the retrieve mode follows from the prefix. +CATALOGUES = INPUTS.get("catalogues") or "" +if TILE_DETECTION == "unions_catalogue" and not CATALOGUES: + raise WorkflowError( + "tile_detection=unions_catalogue needs `inputs.catalogues:` (a local " + "directory or vos: URL holding the per-tile CFIS..r.cat files).") +CATALOGUE_RETRIEVE = "vos" if CATALOGUES.startswith("vos:") else "symlink" # 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 @@ -622,7 +639,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 fdc3d4f59..bb3565fa5 100755 --- a/workflow/bin/sp +++ b/workflow/bin/sp @@ -42,7 +42,7 @@ source "$VENV/bin/activate" cfg() { python -c 'import sys, yaml value = yaml.safe_load(open(sys.argv[1])) for key in sys.argv[2].split("."): - value = value[key] + value = (value or {}).get(key) print(value or "")' "$CONFIG" "$1"; } RUN_DIR="$(cfg outputs.run_dir)"; INDEX_DB="$(cfg outputs.index_db)" @@ -206,6 +206,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 63b29413c..8c3a0e58e 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -24,6 +24,17 @@ inputs: # Pre-staged P3 data (get_images RETRIEVE=symlink). tiles: /project/def-mjhudson/unions-wl/tiles exposures: /project/def-mjhudson/unions-wl/exposures + # The UNIONS per-tile object catalogues (CFIS..r.cat), read only when + # tile_detection is unions_catalogue: a local directory (symlinked in) or a + # vos: URL (downloaded, e.g. vos:cfis/tiles_DR6). + # catalogues: vos:cfis/tiles_DR6 + +# Where the tile's galaxy sample comes from: `sextractor` runs SExtractor on +# the tile image; `unions_catalogue` fetches the UNIONS per-tile catalogue +# (inputs.catalogues) and converts it in place, carrying its unique object ID +# into the final catalogue as TILE_UNIQUE_ID. The catalogue path has no +# segmentation map, so ngmix's BLEND_HANDLING = uberseg needs sextractor. +tile_detection: sextractor # The container every job runs inside (apptainer software-deployment in the profile). container: /project/def-mjhudson/cdaley/containers/shapepipe-develop-runtime.sif diff --git a/workflow/config/cfis/final_cat.param b/workflow/config/cfis/final_cat.param index 00bcb3f73..9bb1cc285 100644 --- a/workflow/config/cfis/final_cat.param +++ b/workflow/config/cfis/final_cat.param @@ -6,6 +6,11 @@ YWIN_WORLD # Can maybe be removed. TILE_ID +# UNIONS per-tile unique object ID: present only in catalogues produced with +# tile_detection: unions_catalogue (workflow/config.yaml). create_final_cat.py +# needs every listed column in every tile, so list it for such runs only. +#TILE_UNIQUE_ID + # flags FLAGS IMAFLAGS_ISO diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index 68d6da9ba..be00579e3 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -41,6 +41,20 @@ 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 +(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, the multi-epoch post-processing, and the catalogue's per-tile +unique object ID as ``TILE_UNIQUE_ID`` -- a column make_cat carries into the +final catalogue with every other sexcat column. 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. What +this path does not have is a segmentation map: ngmix's ``BLEND_HANDLING = +uberseg`` needs the SExtractor mode. + 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 @@ -444,24 +458,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), 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 diff --git a/workflow/scripts/completeness.py b/workflow/scripts/completeness.py index ff776e6e7..061499f2b 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. COMPLETENESS = { # --- tile prepare (phase A) --- # get_images counts are CONFIG-FLAVOR-DEPENDENT: the v2.0 bash table said 4/6 @@ -122,7 +129,13 @@ # --- tile post --- "tile_merge_headers": {"merge_headers_runner": dict(expect=1)}, - "tile_detect": {"sextractor_runner": dict(expect=2)}, + # The fetched UNIONS catalogue: one .cat per tile. + "tile_get_catalogue": {"get_images_runner": dict(expect=1)}, + "tile_detect": { + "sextractor": {"sextractor_runner": dict(expect=2)}, + # One 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), @@ -190,6 +203,13 @@ def check_counts(stage, run_dir): raise ValueError( f"Invalid SP_PSF={psf_model!r}; expected one of psfex, mccd." ) 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) @@ -221,6 +241,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/run_report.py b/workflow/scripts/run_report.py index 24c07670b..2462cdae0 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"] # The manifests clean_tile leaves on disk (workflow/scripts/clean_tile.py names From 379be4850b85f690653fd2a75a9398dc4e2884ca Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 17 Sep 2026 13:15:56 +0200 Subject: [PATCH 03/21] Test the UNIONS catalogue conversion and its workflow wiring read_ext_sexcat on a synthetic ASCII catalogue and image: LDAC layout, the tile header in LDAC_IMHEAD, the SExtractor aliases, zero-padded VIGNET stamps cut at the right pixels, TILE_UNIQUE_ID = tile_id * 1e6 + NUMBER, and that column surviving make_cat.save_sextractor_data into the final catalogue. On the workflow side: the two detection inis share the stage's run dir, the fetch ini feeds the converter, the tile_detect count table is selected per mode with sextractor as the unset default, and run_report lists the fetch stage only for catalogue runs. Co-Authored-By: Claude Fable 5.1 --- tests/module/test_read_ext_sexcat.py | 125 +++++++++++++++++++++ tests/unit/test_workflow_tile_detection.py | 118 +++++++++++++++++++ 2 files changed, 243 insertions(+) create mode 100644 tests/module/test_read_ext_sexcat.py create mode 100644 tests/unit/test_workflow_tile_detection.py diff --git a/tests/module/test_read_ext_sexcat.py b/tests/module/test_read_ext_sexcat.py new file mode 100644 index 000000000..0377755ec --- /dev/null +++ b/tests/module/test_read_ext_sexcat.py @@ -0,0 +1,125 @@ +"""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 ``TILE_UNIQUE_ID``. The +last test follows that column through ``make_cat.save_sextractor_data`` into +the final catalogue. +""" + +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), "-301-279", 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_tile_unique_id_and_aliases(ldac): + with fits.open(ldac) as hdul: + data = hdul["LDAC_OBJECTS"].data + numbers = np.array([o[0] for o in OBJECTS]) + npt.assert_array_equal(data["NUMBER"], numbers) + npt.assert_array_equal(data["TILE_UNIQUE_ID"], 301279 * 10**6 + numbers) + assert data["TILE_UNIQUE_ID"].dtype == np.int64 + 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_zero_padded(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 zero. + assert (vignets[1, :, :2] == 0).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:, :] == 0).all() + assert (vignets[2, :, STAMP // 2 + 1:] == 0).all() + assert vignets[2, STAMP // 2, STAMP // 2] == 29 * 1000 + 39 + + +@pytest.mark.parametrize("number", ["-1301-279", "-301-1279"]) +def test_four_digit_grid_index_is_rejected(number): + with pytest.raises(ValueError): + rs._tile_id_from_file_number_string(number) + + +def test_tile_unique_id_reaches_the_final_catalogue(ldac, tmp_path): + """make_cat copies every sexcat column, so the ID needs no extra wiring.""" + 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) diff --git a/tests/unit/test_workflow_tile_detection.py b/tests/unit/test_workflow_tile_detection.py new file mode 100644 index 000000000..744dcb595 --- /dev/null +++ b/tests/unit/test_workflow_tile_detection.py @@ -0,0 +1,118 @@ +"""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 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 From affcdc099c88018b569870fd9bd1ede5ab01d8ce Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 17 Sep 2026 13:29:12 +0200 Subject: [PATCH 04/21] Fix TILE_UNIQUE_ID merge, NUMBER contiguity, and FITS int64 dtype test final_cat.param lists TILE_UNIQUE_ID again; merge_final_cat.py and create_final_cat.py now drop a requested column that is absent from a given catalogue (logging once) instead of failing, so sextractor-mode runs (without the column) and unions_catalogue-mode runs (with it) both merge cleanly. read_ext_sexcat's NUMBER passthrough is not guaranteed contiguous or in order in the external catalogue, which ngmix_range requires; TILE_UNIQUE_ID already preserves the original NUMBER, so the converter now renumbers NUMBER to 1..n_obj in output row order and says so in its docstring. test_read_ext_sexcat's TILE_UNIQUE_ID dtype check compared against np.int64 directly; a FITS 'K' column reads back big-endian, so the check now compares dtype.kind and itemsize instead. Co-Authored-By: Claude Fable 5.1 --- scripts/python/create_final_cat.py | 38 +++++--- scripts/python/merge_final_cat.py | 45 +++++++++ .../read_ext_sexcat.py | 13 ++- tests/module/test_read_ext_sexcat.py | 9 +- tests/unit/test_merge_final_cat.py | 92 +++++++++++++++++++ workflow/config/cfis/final_cat.param | 7 +- 6 files changed, 185 insertions(+), 19 deletions(-) create mode 100644 tests/unit/test_merge_final_cat.py diff --git a/scripts/python/create_final_cat.py b/scripts/python/create_final_cat.py index 2b583b857..577c2b559 100755 --- a/scripts/python/create_final_cat.py +++ b/scripts/python/create_final_cat.py @@ -297,9 +297,17 @@ def get_patch_group(hdf5_file, patch, verbose=False): return patch_group +_warned_missing_columns = set() + + def read_data(fits_file, params): """Read Data. + Some requested columns (e.g. ``TILE_UNIQUE_ID``, present only for tiles + produced with ``tile_detection: unions_catalogue``) may be absent from a + given tile's catalogue; such columns are skipped for that tile, with one + log line the first time each is seen missing. + """ with fits.open(fits_file) as hdu_list: try: @@ -312,18 +320,22 @@ def read_data(fits_file, params): if params["param_list"] is None: params["param_list"] = [col for col in data.keys()] - try: - extracted_data = {col: data[col] for col in params["param_list"]} - dtype = data.dtype - except: - print(f"Error for ID {id}, path {fits_file}") - for col in params["param_list"]: - if col not in data: - print(col, end=" ") - print() - continue + available = set(data.dtype.names) + param_list = [] + for col in params["param_list"]: + if col in available: + param_list.append(col) + elif col not in _warned_missing_columns: + print( + f"Column '{col}' not found in input catalogue " + f"(e.g. {fits_file}), skipping in merged catalogue" + ) + _warned_missing_columns.add(col) + + extracted_data = {col: data[col] for col in param_list} + dtype = np.dtype([(col, data.dtype[col]) for col in param_list]) - return extracted_data, dtype + return extracted_data, dtype, param_list def copy_data(param_list, extracted_data, dtype): @@ -463,9 +475,9 @@ def process(params): print(f"Run without output file found for {id}, skipping") continue - extracted_data, dtype = read_data(fits_file, params) + extracted_data, dtype, param_list = read_data(fits_file, params) - structured_data = copy_data(params["param_list"], extracted_data, dtype) + structured_data = copy_data(param_list, extracted_data, dtype) # Create a new dataset try: diff --git a/scripts/python/merge_final_cat.py b/scripts/python/merge_final_cat.py index 9d9c10d21..422df61a8 100755 --- a/scripts/python/merge_final_cat.py +++ b/scripts/python/merge_final_cat.py @@ -247,6 +247,43 @@ def read_param_file(path, verbose=False): return param_list +def filter_available_columns(param_list, available_columns): + """Filter Available Columns. + + Return the subset of ``param_list`` present in ``available_columns``, + printing one line for each requested column that is missing (e.g. + ``TILE_UNIQUE_ID``, present only in catalogues produced with + ``tile_detection: unions_catalogue``) so the merge still proceeds for + catalogues produced with other tile-detection settings. + + Parameters + ---------- + param_list: list of str + requested column names + available_columns: iterable of str + column names present in the catalogue + + Returns + ------- + list of str + subset of ``param_list`` present in ``available_columns`` + + """ + if not param_list: + return param_list + + available_columns = set(available_columns) + missing = [p for p in param_list if p not in available_columns] + + for p in missing: + print( + f"Column '{p}' not found in input catalogue, skipping in " + "merged catalogue" + ) + + return [p for p in param_list if p in available_columns] + + def get_data(path, hdu_num, param_list): """Get Data. @@ -345,6 +382,14 @@ def main(argv=None): if param.verbose: print(f"{len(lpath)} files files to merge found") + # Drop requested columns absent from the catalogues (e.g. TILE_UNIQUE_ID + # for tile_detection: sextractor runs) so the merge still proceeds. + with fits.open(lpath[0]) as hdu_list: + available_columns = hdu_list[param.hdu_num].columns.names + param.param_list = filter_available_columns( + param.param_list, available_columns + ) + count = 0 # Determine number of columns and keys from first catalogue file 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..c117ef243 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 @@ -219,7 +219,13 @@ def make_ldac_from_ascii( 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``). + ``file_number_string`` (e.g. ``'-301-279'`` → ``tile_id = 301279``) and + ``NUMBER`` is the input catalogue's original object number, so + ``TILE_UNIQUE_ID`` preserves the identity assigned upstream. The + ``NUMBER`` column written to ``LDAC_OBJECTS`` is then overwritten with a + running index ``1..n_obj`` in output row order: downstream ShapePipe + (e.g. ``ngmix_range``) assumes a contiguous, in-order ``NUMBER``, which + the input catalogue is not guaranteed to have. Parameters ---------- @@ -247,6 +253,11 @@ def make_ldac_from_ascii( tile_id * 10**6 + np.array(cat_data["NUMBER"], dtype=np.int64) ) + # NUMBER is not guaranteed to be a contiguous, in-order running index in + # the input catalogue; downstream code (ngmix_range) requires exactly + # that, and the original identity is preserved above in TILE_UNIQUE_ID. + cat_data["NUMBER"] = np.arange(1, n_obj + 1, dtype=np.int64) + with fits.open(image_path) as hdul: img_header = hdul[0].header image_data = hdul[0].data.astype(np.float32) diff --git a/tests/module/test_read_ext_sexcat.py b/tests/module/test_read_ext_sexcat.py index 0377755ec..a17527d6b 100644 --- a/tests/module/test_read_ext_sexcat.py +++ b/tests/module/test_read_ext_sexcat.py @@ -76,9 +76,14 @@ def test_tile_unique_id_and_aliases(ldac): with fits.open(ldac) as hdul: data = hdul["LDAC_OBJECTS"].data numbers = np.array([o[0] for o in OBJECTS]) - npt.assert_array_equal(data["NUMBER"], numbers) + # NUMBER is renumbered to a contiguous 1..n_obj running index in output + # row order; the original NUMBER survives only inside TILE_UNIQUE_ID. + npt.assert_array_equal(data["NUMBER"], np.arange(1, len(OBJECTS) + 1)) npt.assert_array_equal(data["TILE_UNIQUE_ID"], 301279 * 10**6 + numbers) - assert data["TILE_UNIQUE_ID"].dtype == np.int64 + # A FITS 'K' column reads back as big-endian ('>i8'), not np.int64's + # native byte order, so compare kind and width rather than dtype. + assert data["TILE_UNIQUE_ID"].dtype.kind == "i" + assert data["TILE_UNIQUE_ID"].dtype.itemsize == 8 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"]) diff --git a/tests/unit/test_merge_final_cat.py b/tests/unit/test_merge_final_cat.py new file mode 100644 index 000000000..981c38a52 --- /dev/null +++ b/tests/unit/test_merge_final_cat.py @@ -0,0 +1,92 @@ +"""``merge_final_cat.py`` and ``create_final_cat.py`` drop catalogue-param +columns absent from a given input catalogue. + +``TILE_UNIQUE_ID`` is only present in per-tile catalogues produced with +``tile_detection: unions_catalogue``; a run using ``tile_detection: +sextractor`` never has it. Listing it in ``final_cat.param`` must not break +the merge for such runs: ``merge_final_cat.filter_available_columns`` and +``create_final_cat.read_data`` are what make that column (and any other +requested-but-absent column) optional, logging one line per drop instead of +failing partway through the merge. +""" + +import importlib.util +from pathlib import Path + +import numpy as np +from astropy.io import fits + +REPO_ROOT = Path(__file__).resolve().parents[2] +MERGE_SCRIPT = REPO_ROOT / "scripts" / "python" / "merge_final_cat.py" +CREATE_SCRIPT = REPO_ROOT / "scripts" / "python" / "create_final_cat.py" + + +def _load(script): + """Import a script by path — ``scripts/python`` is not a package.""" + assert script.exists(), f"{script} not found; the rule calls it by path" + spec = importlib.util.spec_from_file_location(f"_{script.stem}", script) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def test_present_columns_are_kept(): + m = _load(MERGE_SCRIPT) + param_list = ["XWIN_WORLD", "TILE_ID"] + assert m.filter_available_columns( + param_list, ["XWIN_WORLD", "TILE_ID", "FLAGS"] + ) == param_list + + +def test_missing_column_is_dropped_and_logged(capsys): + m = _load(MERGE_SCRIPT) + kept = m.filter_available_columns( + ["XWIN_WORLD", "TILE_UNIQUE_ID"], ["XWIN_WORLD", "FLAGS"] + ) + assert kept == ["XWIN_WORLD"] + out = capsys.readouterr().out + assert "TILE_UNIQUE_ID" in out + + +def test_empty_param_list_means_copy_all_columns(): + m = _load(MERGE_SCRIPT) + assert m.filter_available_columns([], ["XWIN_WORLD"]) == [] + + +def _write_cat(path, columns): + cols = [ + fits.Column(name=name, format="K", array=np.array(values, dtype=np.int64)) + for name, values in columns.items() + ] + hdu = fits.BinTableHDU.from_columns(cols) + fits.HDUList([fits.PrimaryHDU(), hdu]).writeto(path, overwrite=True) + + +def test_create_final_cat_read_data_skips_missing_column(tmp_path, capsys): + m = _load(CREATE_SCRIPT) + + cat_path = tmp_path / "final_cat-100-100.fits" + _write_cat(cat_path, {"XWIN_WORLD": [1, 2], "TILE_ID": [100, 100]}) + + params = {"hdu_num": 1, "param_list": ["XWIN_WORLD", "TILE_UNIQUE_ID"]} + extracted_data, dtype, param_list = m.read_data(str(cat_path), params) + + assert param_list == ["XWIN_WORLD"] + assert set(extracted_data) == {"XWIN_WORLD"} + assert "TILE_UNIQUE_ID" not in dtype.names + assert "TILE_UNIQUE_ID" in capsys.readouterr().out + + +def test_create_final_cat_read_data_keeps_present_column(tmp_path): + m = _load(CREATE_SCRIPT) + + cat_path = tmp_path / "final_cat-100-100.fits" + _write_cat( + cat_path, + {"XWIN_WORLD": [1, 2], "TILE_UNIQUE_ID": [100000001, 100000002]}, + ) + params = {"hdu_num": 1, "param_list": ["XWIN_WORLD", "TILE_UNIQUE_ID"]} + extracted_data, dtype, param_list = m.read_data(str(cat_path), params) + + assert param_list == ["XWIN_WORLD", "TILE_UNIQUE_ID"] + assert set(dtype.names) == {"XWIN_WORLD", "TILE_UNIQUE_ID"} diff --git a/workflow/config/cfis/final_cat.param b/workflow/config/cfis/final_cat.param index 9bb1cc285..f51f0e767 100644 --- a/workflow/config/cfis/final_cat.param +++ b/workflow/config/cfis/final_cat.param @@ -7,9 +7,10 @@ YWIN_WORLD TILE_ID # UNIONS per-tile unique object ID: present only in catalogues produced with -# tile_detection: unions_catalogue (workflow/config.yaml). create_final_cat.py -# needs every listed column in every tile, so list it for such runs only. -#TILE_UNIQUE_ID +# tile_detection: unions_catalogue (workflow/config.yaml). merge_final_cat.py +# drops this column (with a log line) when merging catalogues produced with +# tile_detection: sextractor, where it is absent. +TILE_UNIQUE_ID # flags FLAGS From 930b5afc970ad9b77bd975d2761c559940b4724f Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 17 Sep 2026 15:04:58 +0200 Subject: [PATCH 05/21] Build the UberSeg segmentation map for the UNIONS catalogue path UberSeg identifies a stamp's central object BY LABEL: uberseg_weight keeps the pixels whose nearest segmentation footprint carries the object's own catalogue NUMBER. The UNIONS per-tile catalogue carries no segmentation map, and a SExtractor run's labels are its own detection numbers, so the two numberings have to be bridged before ngmix can mask a neighbour. `blend_handling:` (noisefill | uberseg) joins `tile_detection:` in the run config. The pair unions_catalogue + uberseg adds tile_segmentation: one SExtractor run on the tile image whose only product is the SEGMENTATION check image (config_tile_Sg.ini), followed by seg_relabel.py, which gives each catalogue object the footprint its position falls in, relabels it with that object's NUMBER, marks every unclaimed footprint a neighbour, and paints a small disc for an object that falls on sky or on a footprint already claimed. The result lands beside the sexcat, where vignetmaker's segmentation run looks whichever mode wrote the catalogue. Steven's catalogue stays the sample; SExtractor runs once, for the one product the catalogue cannot give. Co-Authored-By: Claude Fable 5.1 --- tests/unit/test_seg_relabel.py | 145 ++++++++++++++++ tests/unit/test_workflow_tile_detection.py | 83 +++++++++ workflow/Snakefile | 21 ++- workflow/bin/sp | 1 + workflow/config/cfis/config_tile_Sg.ini | 131 ++++++++++++++ workflow/rules/tile.smk | 72 +++++++- workflow/scripts/completeness.py | 11 ++ workflow/scripts/run_report.py | 6 + workflow/scripts/seg_relabel.py | 191 +++++++++++++++++++++ 9 files changed, 658 insertions(+), 3 deletions(-) create mode 100644 tests/unit/test_seg_relabel.py create mode 100644 workflow/config/cfis/config_tile_Sg.ini create mode 100644 workflow/scripts/seg_relabel.py diff --git a/tests/unit/test_seg_relabel.py b/tests/unit/test_seg_relabel.py new file mode 100644 index 000000000..56d9600d1 --- /dev/null +++ b/tests/unit/test_seg_relabel.py @@ -0,0 +1,145 @@ +"""The segmentation map UberSeg sees must call THIS object self. + +``uberseg_weight`` keeps the pixels whose nearest footprint carries the +object's catalogue ``NUMBER`` and zeros the rest; it never looks at the label +under the object's position. So the one thing that has to hold, for every +object in the UNIONS catalogue, is that the relabelled map carries that +object's ``NUMBER`` on its own pixels and something else everywhere else. +These tests assert exactly that, on maps small enough to read by eye — +SExtractor is not involved. +""" + +import importlib.util +import sys +from pathlib import Path + +import numpy as np +import pytest + +SCRIPTS = Path(__file__).resolve().parents[2] / "workflow" / "scripts" + + +def _load(): + sys.path.insert(0, str(SCRIPTS)) + try: + spec = importlib.util.spec_from_file_location( + "_seg_relabel", SCRIPTS / "seg_relabel.py") + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + finally: + sys.path.remove(str(SCRIPTS)) + return module + + +seg_relabel = _load() + + +def _map(): + """Two SExtractor 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 + + +def test_each_object_owns_the_footprint_it_sits_in(): + """The claim is by position; the label that comes out is the NUMBER.""" + seg = _map() + # FITS 1-indexed centres of the two footprints, in an order that is NOT + # the label order, so a relabelling that merely renumbered would fail. + out, counts = seg_relabel.relabel( + 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(): + """A detection no catalogue object sits in is masked, not measured.""" + seg = _map() + out, counts = seg_relabel.relabel( + seg, number=np.array([1]), + x_image=np.array([4.0]), y_image=np.array([4.0])) + assert counts["matched"] == 1 + assert set(np.unique(out[seg == 9])) == {seg_relabel.NEIGHBOUR_LABEL} + assert seg_relabel.NEIGHBOUR_LABEL not in np.unique(out[seg == 7]) + + +def test_an_object_on_sky_still_gets_a_self(): + """No footprint at its position -> a disc of its own NUMBER. + + Without one, uberseg_weight would zero the whole stamp and + Ngmix._check_central_seg_label would raise on a map lacking the label. + """ + seg = _map() + out, counts = seg_relabel.relabel( + 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 + # Centred on the object, and nowhere else. + 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(): + """The first claimant keeps the blend; the second takes back its centre.""" + seg = _map() + out, counts = seg_relabel.relabel( + 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 + assert out[15, 15] == 2 + # The claimant still holds the bulk of the footprint. + assert np.count_nonzero(out == 1) > np.count_nonzero(out == 2) + + +def test_every_object_is_self_somewhere(): + """The invariant, over a randomised map: every object on the image keeps + pixels of its own — matched, unmatched, blended or doubled.""" + 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 = seg_relabel.relabel(seg, number, x_image, y_image) + # The four claim outcomes partition the catalogue; shared_pixel is a + # separate axis, counted on top. + assert sum(counts[k] for k in + ("matched", "unclaimed", "shared", "off_image")) == len(number) + assert counts["off_image"] == 0 + for num in number: + assert np.any(out == num), f"object {num} has no self pixels" + # And nothing outside a footprint or a disc was invented. + assert set(np.unique(out)) <= set(number) | {0, seg_relabel.NEIGHBOUR_LABEL} + + +def test_two_objects_on_one_pixel_split_the_contest(): + """A pixel has one label: the lower NUMBER keeps it, the other keeps its + disc, and the collision is counted rather than hidden.""" + seg = _map() + out, counts = seg_relabel.relabel( + 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, 61.0)]) +def test_a_position_off_the_image_is_counted_not_crashed(x, y): + seg = _map() + out, counts = seg_relabel.relabel( + seg, 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) diff --git a/tests/unit/test_workflow_tile_detection.py b/tests/unit/test_workflow_tile_detection.py index 744dcb595..eb7695128 100644 --- a/tests/unit/test_workflow_tile_detection.py +++ b/tests/unit/test_workflow_tile_detection.py @@ -12,6 +12,7 @@ import configparser import importlib.util +import re import sys from pathlib import Path @@ -116,3 +117,85 @@ def test_report_lists_the_fetch_stage_only_for_catalogue_runs(monkeypatch, mode, assert ("tile_get_catalogue" in stages) is present if present: assert stages.index("tile_get_catalogue") == stages.index("tile_detect") - 1 + + +# --- the segmentation half: unions_catalogue + uberseg ---------------------- + + +def test_blend_handlings_mirror_the_ngmix_module(): + """completeness.BLEND_HANDLINGS is a copy; the copy must stay true. + + The Snakefile validates `blend_handling:` against it in the launcher venv, + outside the container where shapepipe is importable, so the tuple is + mirrored rather than imported. Read out of the source text here for the + same reason: this file is container-free. + """ + src = (REPO_ROOT / "src" / "shapepipe" / "modules" / "ngmix_package" + / "ngmix.py").read_text() + match = re.search(r"^BLEND_HANDLINGS = \(([^)]*)\)", src, re.M) + assert match, "ngmix.py no longer defines BLEND_HANDLINGS" + assert completeness.BLEND_HANDLINGS == tuple( + part.strip().strip('"\'') for part in match.group(1).split(",") + if part.strip()) + + +def test_segmentation_ini_matches_its_stage_and_yields_the_map(): + level, sg = completeness.STAGE_DIR["tile_segmentation"] + assert level == "tile" + ini = _ini("config_tile_Sg.ini") + assert ini["DEFAULT"]["RUN_NAME"].strip() == sg + assert ini["DEFAULT"]["RUN_DATETIME"].strip() == "False" + sx = ini["SEXTRACTOR_RUNNER"] + # The check image is the whole point, and it is the only one asked for. + assert [c.strip() for c in sx["CHECKIMAGE"].split(",")] == ["SEGMENTATION"] + # No multi-epoch post-processing: those extensions ride on the sexcat + # config_tile_Uc.ini writes, and without it the WCS sqlite is not an input. + assert sx["MAKE_POST_PROCESS"].strip() == "False" + assert "log_exp_headers" not in sx["FILE_PATTERN"] + # Same detection settings as the SExtractor mode, so the footprints are + # the ones that mode would have drawn. + detect = _ini("config_tile_Sx.ini")["SEXTRACTOR_RUNNER"] + for key in ("DOT_SEX_FILE", "DOT_PARAM_FILE", "DOT_CONV_FILE", + "WEIGHT_IMAGE", "FLAG_IMAGE"): + assert sx[key].strip() == detect[key].strip(), key + + +def test_segmentation_stage_is_checked(tmp_path): + """Two products: the SEGMENTATION check image and SExtractor's own sexcat.""" + ok, details = completeness.check_counts( + "tile_segmentation", _stage_dir(tmp_path, "sextractor_runner", 2)) + assert ok and details == [("sextractor_runner", 2, 2, False)] + + +@pytest.mark.parametrize("detection, blend, present", [ + ("unions_catalogue", "uberseg", True), + ("unions_catalogue", "noisefill", False), + ("sextractor", "uberseg", False), +]) +def test_report_lists_the_segmentation_stage_only_for_the_pair( + monkeypatch, detection, blend, present): + monkeypatch.setenv("SP_TILE_DETECTION", detection) + monkeypatch.setenv("SP_BLEND_HANDLING", blend) + stages = _load("run_report").TILE_STAGES + assert ("tile_segmentation" in stages) is present + if present: + assert stages.index("tile_segmentation") == stages.index( + "tile_detect") + 1 + + +def test_run_config_defaults_to_the_catalogue_and_declares_its_source(): + """The committed run config must parse under its own defaults. + + `tile_detection: unions_catalogue` is refused by the Snakefile without + `inputs.catalogues`, so the two settings travel together. + """ + text = (REPO_ROOT / "workflow" / "config.yaml").read_text() + config = {} + for line in text.splitlines(): + if line.startswith("tile_detection:") or line.startswith( + "blend_handling:"): + key, _, value = line.partition(":") + config[key] = value.strip() + assert config["tile_detection"] == "unions_catalogue" + assert config["blend_handling"] in completeness.BLEND_HANDLINGS + assert re.search(r"^ catalogues: \S+", text, re.M) diff --git a/workflow/Snakefile b/workflow/Snakefile index d7091451e..ac86c5580 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -100,7 +100,8 @@ CONFIG_DIR = Path(workflow.basedir) / "config" / "cfis" sys.path.insert(0, str(SCRIPTS)) import build_index # noqa: E402 -from completeness import STAGE_DIR, TILE_DETECTIONS # noqa: E402 +from completeness import (BLEND_HANDLINGS, STAGE_DIR, # noqa: E402 + TILE_DETECTIONS) # 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). @@ -119,6 +120,20 @@ if TILE_DETECTION == "unions_catalogue" and not CATALOGUES: "directory or vos: URL holding the per-tile CFIS..r.cat files).") CATALOGUE_RETRIEVE = "vos" if CATALOGUES.startswith("vos:") else "symlink" +# How ngmix treats a stamp's neighbours (its BLEND_HANDLING option): the +# noise-fill treatment needs nothing of the workflow, while uberseg needs a +# segmentation map labelled in the numbering of the catalogue being measured. +# SExtractor detection writes that map as a by-product (config_tile_Sx.ini's +# CHECKIMAGE); the UNIONS catalogue has none, so the pair below is what adds +# tile.smk's tile_segmentation rule. +BLEND_HANDLING = config.get("blend_handling", "noisefill") +if BLEND_HANDLING not in BLEND_HANDLINGS: + raise WorkflowError( + f"Invalid blend_handling={BLEND_HANDLING!r}; expected one of " + f"{sorted(BLEND_HANDLINGS)}.") +TILE_SEGMENTATION = (TILE_DETECTION == "unions_catalogue" + and BLEND_HANDLING == "uberseg") + # 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 # --dag`, `sp exp_psf ...`). It gates the two parse-time side effects — the index @@ -353,6 +368,9 @@ CLEAN_TILE_HASH = script_hash("clean_tile.py") # fresh root: a resume across it re-measures the tile, and the failure modes a # mid-campaign params change can reach are catalogued at tile.smk's range_hash. NGMIX_RANGE_HASH = script_hash("ngmix_range.py") +# seg_relabel.py decides which footprint is "self" for every object UberSeg +# masks, so an edit to it changes every mask; it rides on tile_segmentation. +SEG_RELABEL_HASH = script_hash("seg_relabel.py") # --- exposure reclamation (D5, S5) ----------------------------------------- @@ -640,6 +658,7 @@ rule prepare_all_tiles: # via `sp report`. def _report(status): shell(f"SP_TILE_DETECTION='{TILE_DETECTION}' " + f"SP_BLEND_HANDLING='{BLEND_HANDLING}' " f"python {SCRIPTS}/run_report.py --run-dir {RUN_DIR} " f"--index {INDEX_DB} --status {status} || true") diff --git a/workflow/bin/sp b/workflow/bin/sp index bb3565fa5..d88138db6 100755 --- a/workflow/bin/sp +++ b/workflow/bin/sp @@ -207,6 +207,7 @@ case "$cmd" in # either way -- this is consistency, not safety. shift SP_TILE_DETECTION="$(cfg tile_detection)" \ + SP_BLEND_HANDLING="$(cfg blend_handling)" \ python "$(code_root)/workflow/scripts/run_report.py" \ --run-dir "$RUN_DIR" --index "$INDEX_DB" \ --status manual "$@" diff --git a/workflow/config/cfis/config_tile_Sg.ini b/workflow/config/cfis/config_tile_Sg.ini new file mode 100644 index 000000000..b02592fd5 --- /dev/null +++ b/workflow/config/cfis/config_tile_Sg.ini @@ -0,0 +1,131 @@ +# ShapePipe configuration file for the tile SEGMENTATION map, used when +# `tile_detection: unions_catalogue` is combined with `blend_handling: uberseg` +# in the workflow run config. +# +# The UNIONS per-tile catalogue is the galaxy sample; this run exists only for +# the one product that catalogue cannot carry, the segmentation map UberSeg +# needs. Its own sexcat is a by-product and nothing reads it. The detection +# settings are those of config_tile_Sx.ini, so the footprints are the ones a +# SExtractor-detection run would have drawn. +# +# The labels of the SEGMENTATION check image are this run's own NUMBERs, which +# are not the converted catalogue's; workflow/scripts/seg_relabel.py maps them +# into the catalogue's NUMBER space and writes the result beside the sexcat. + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Sg + +# Add date and time to RUN_NAME, optional, default: True +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = sextractor_runner + +# Run mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +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 + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 16 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[SEXTRACTOR_RUNNER] + +# The tile image and its weight. The merged WCS headers are absent because +# MAKE_POST_PROCESS is False: the multi-epoch extensions belong to the sexcat +# the chain reads, which config_tile_Uc.ini writes. +INPUT_DIR = $SP_RUN/output/run_sp_tile_Git/get_images_runner/output, $SP_RUN/output/run_sp_tile_Uz/uncompress_fits_runner/output + +FILE_PATTERN = CFIS_image, CFIS_weight + +FILE_EXT = .fits, .fits + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +# SExtractor executable path +EXEC_PATH = source-extractor + +# SExtractor configuration files +DOT_SEX_FILE = $SP_CONFIG/default_tile.sex +DOT_PARAM_FILE = $SP_CONFIG/default_noimaflags.param +DOT_CONV_FILE = $SP_CONFIG/default.conv + +# Use input weight image if True +WEIGHT_IMAGE = True + +# Use input flag image if True +FLAG_IMAGE = False + +# Use input PSF file if True +PSF_FILE = False + +# Use distinct image for detection (SExtractor in +# dual-image mode) if True +DETECTION_IMAGE = False + +# Distinct weight image for detection (SExtractor +# in dual-image mode) +DETECTION_WEIGHT = False + +ZP_FROM_HEADER = False + +BKG_FROM_HEADER = False + +# Type of image check (optional), default not used, can be a list of +# BACKGROUND, BACKGROUND_RMS, INIBACKGROUND, +# MINIBACK_RMS, -BACKGROUND, #FILTERED, +# OBJECTS, -OBJECTS, SEGMENTATION, APERTURES +# +# SEGMENTATION alone: this run's reason to exist. BACKGROUND is the tile +# background config_tile_Sx.ini writes for the detection chain, and nothing +# reads it from here. +CHECKIMAGE = SEGMENTATION + +# File name suffix for the output sextractor files (optional) +SUFFIX = sexcat + +## Post-processing + +# The multi-epoch extensions ride on the sexcat the chain reads, not on this +# run's by-product catalogue. +MAKE_POST_PROCESS = False diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index be00579e3..07253e52b 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -52,8 +52,14 @@ final catalogue with every other sexcat column. 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. What -this path does not have is a segmentation map: ngmix's ``BLEND_HANDLING = -uberseg`` needs the SExtractor mode. +Steven's catalogue does not carry is a segmentation map, and ngmix's +``BLEND_HANDLING = uberseg`` needs one. ``blend_handling: uberseg`` therefore +adds a third rule, tile_segmentation: a SExtractor run on the tile image whose +only product is the SEGMENTATION check image (config_tile_Sg.ini), followed by +seg_relabel.py, which maps that image's labels into the converted catalogue's +NUMBER space and writes it beside the sexcat -- the path vignetmaker's +segmentation run reads. SExtractor runs once, for the one thing the catalogue +cannot give. 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 @@ -540,6 +546,67 @@ else: '"$SP_RUN/output/run_sp_tile_Sx/sextractor_runner"\n' "fi\n") +# Where the tile catalogue lands, whichever mode wrote it: the one path every +# downstream config reads, and where the relabelled segmentation map goes. +SEXCAT_DIR = "$SP_RUN/output/run_sp_tile_Sx/sextractor_runner/output" + +if TILE_SEGMENTATION: + + # The segmentation map UberSeg needs, and nothing else. + # + # THE SAMPLE IS STEVEN'S CATALOGUE, ALWAYS. This SExtractor run detects on + # the tile image with config_tile_Sx.ini's settings, but its catalogue is + # discarded: what leaves the rule is the SEGMENTATION check image. That + # image's labels are this run's own NUMBERs, which have nothing to do with + # the converted catalogue's, so the post step relabels it -- UberSeg finds + # the central object BY LABEL (uberseg_weight's `object_number` is the + # catalogue NUMBER), so an unmapped map would invert every mask rather than + # fail. seg_relabel.py's docstring owns the mapping and its fallback. + # + # It writes into tile_detect's run dir, which is where vignetmaker's + # segmentation run looks (`run_sp_tile_Sx:sextractor_runner`, FILE_PATTERN + # `segmentation`) whichever mode produced the sexcat. Safe because unit_pre + # clears only THIS stage's run dir (run_sp_tile_Sg), and a tile_detect + # rerun re-schedules this rule through the manifest edge below. + rule tile_segmentation: + input: + sx = rules.tile_detect.output.manifest, + uz = f"{TILE_DIR}/manifests/tile_uncompress.json", + output: + manifest = f"{TILE_DIR}/manifests/tile_segmentation.json" + log: + f"{TILE_DIR}/logs/tile_segmentation.json" + params: + pre = lambda wc: unit_pre("tile_segmentation", wc.tile), + script_hash = SCRIPT_HASH, + relabel_hash = SEG_RELABEL_HASH + threads: 8 + resources: + # The same SExtractor run as tile_detect's sextractor mode. + mem_mb = lambda wc, attempt: 16000 * attempt, + runtime = 180 + shell: + sp_shell( + "tile_segmentation", "config_tile_Sg.ini", + post=( + "if [ $rc -eq 0 ]; then\n" + f" python {SCRIPTS}/seg_relabel.py" + f' --sexcat "{SEXCAT_DIR}/sexcat$SP_UNIT_NUM.fits"' + ' --segmentation "$SP_RUN/output/run_sp_tile_Sg' + '/sextractor_runner/output/segmentation$SP_UNIT_NUM.fits"' + f' --output "{SEXCAT_DIR}/segmentation' + '$SP_UNIT_NUM.fits" || rc=1\n' + "fi\n")) + + +# The segmentation edge: present only when tile_segmentation is defined, so the +# DAG carries it exactly when the map has to be built. +def tile_seg(wc): + if not TILE_SEGMENTATION: + return [] + return [f"{tile_dir(wc.tile)}/manifests/tile_segmentation.json"] + + # Configured PSF interpolation to galaxies + vignet postage stamps: the last # stage that reads exposure products, and the bulk intra-tile intermediate. The store it # writes is node-local (see TILE_LOCAL above). @@ -547,6 +614,7 @@ rule tile_vignets: group: TILE_GROUP input: sx = rules.tile_detect.output.manifest, + seg = tile_seg, forest = rules.tile_exp_forest.output.forest, split = tile_exp_split, psf = tile_exp_psf, diff --git a/workflow/scripts/completeness.py b/workflow/scripts/completeness.py index 061499f2b..c4a630dd3 100644 --- a/workflow/scripts/completeness.py +++ b/workflow/scripts/completeness.py @@ -78,6 +78,13 @@ # when it is not the default, so a SExtractor run's prologue is unchanged. TILE_DETECTIONS = ("sextractor", "unions_catalogue") +# ngmix's neighbour treatments, mirroring BLEND_HANDLINGS in +# shapepipe.modules.ngmix_package.ngmix. Mirrored rather than imported: the +# Snakefile parses this module in the launcher venv, outside the container +# where shapepipe lives. tests/unit/test_workflow_tile_detection.py asserts +# the two tuples agree. +BLEND_HANDLINGS = ("noisefill", "uberseg") + # stage -> {runner_subdir: {expect, [warn], [subpath]}} # exp_psf and tile_vignets are selected by $SP_PSF at check time, tile_detect # by $SP_TILE_DETECTION. @@ -131,6 +138,9 @@ "tile_merge_headers": {"merge_headers_runner": dict(expect=1)}, # The fetched UNIONS catalogue: one .cat per tile. "tile_get_catalogue": {"get_images_runner": dict(expect=1)}, + # The segmentation-only SExtractor run: its SEGMENTATION check image, plus + # the sexcat SExtractor always writes and nothing here reads. + "tile_segmentation": {"sextractor_runner": dict(expect=2)}, "tile_detect": { "sextractor": {"sextractor_runner": dict(expect=2)}, # One FITS-LDAC sexcat, converted from the fetched catalogue. @@ -242,6 +252,7 @@ def check_counts(stage, run_dir): "exp_psf": ("exp", "run_sp_exp_SxSePsf"), "tile_merge_headers": ("tile", "run_sp_tile_Mh_exp"), "tile_get_catalogue": ("tile", "run_sp_tile_Gic"), + "tile_segmentation": ("tile", "run_sp_tile_Sg"), # 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"), diff --git a/workflow/scripts/run_report.py b/workflow/scripts/run_report.py index 2462cdae0..2b9e82f2a 100644 --- a/workflow/scripts/run_report.py +++ b/workflow/scripts/run_report.py @@ -65,6 +65,12 @@ # 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") + # The segmentation run is the second half of that pair: it exists only when + # the converted catalogue has to be measured with ngmix's uberseg blend + # handling, which needs a segmentation map the catalogue does not carry. + if os.environ.get("SP_BLEND_HANDLING") == "uberseg": + TILE_STAGES.insert(TILE_STAGES.index("tile_detect") + 1, + "tile_segmentation") EXP_STAGES = ["exp_get_images", "exp_split", "exp_psf"] # The manifests clean_tile leaves on disk (workflow/scripts/clean_tile.py names diff --git a/workflow/scripts/seg_relabel.py b/workflow/scripts/seg_relabel.py new file mode 100644 index 000000000..bf620a874 --- /dev/null +++ b/workflow/scripts/seg_relabel.py @@ -0,0 +1,191 @@ +"""Map a SExtractor SEGMENTATION map into the tile catalogue's NUMBER space. + +WHY THIS EXISTS. UberSeg identifies the central object of a stamp BY LABEL: it +keeps the pixels whose nearest segmentation footprint carries the object's own +``NUMBER`` and zeros every other pixel's weight +(``ngmix_package/ngmix.py::uberseg_weight``, called with ``object_number = +obj_id``, the catalogue's ``NUMBER``). It never reads the label under the +object's position. So the seg stamp ngmix overlays must be labelled in the +SAME numbering as the catalogue it measures, or every mask is inverted. + +Under ``tile_detection: unions_catalogue`` the two numberings are different by +construction: the sample is Steven Gwyn's per-tile catalogue, renumbered 1..N +by ``read_ext_sexcat``, while the segmentation map comes from a separate +SExtractor run whose labels are ITS OWN detection numbers. This script is the +bridge, and it is the whole of the guarantee: + + * each catalogue object claims the footprint its position falls in -- the one + lookup that is meaningful, since the two runs share the tile's pixel grid; + * that footprint is relabelled with the object's catalogue ``NUMBER``; + * every footprint no catalogue object claims is relabelled ``NEIGHBOUR_LABEL`` + (negative, so it can never collide with a ``NUMBER``): UberSeg's partition + only ever asks "self or not self", so neighbours need no identity; + * an object whose position falls on sky, or on a footprint another object + claimed first, is given a small disc of its own ``NUMBER`` at its position. + +THE FALLBACK IS NOT COSMETIC. Without a footprint of its own an object has no +"self" in the Voronoi partition: ``uberseg_weight`` would hand it an all-zero +weight, and ``Ngmix._check_central_seg_label`` raises on a stamp that does not +contain the object's label at all. The disc gives such an object a defined, +conservative core -- small, centred, and (in the shared-footprint case) taking +its pixels from the neighbour that claimed the blend, which is the right way +round: the claimant keeps the bulk, the unclaimed object keeps a centre. + +Claiming is in ``NUMBER`` order and first come first served, so the mapping is +deterministic and reproducible from the same two inputs. + +Stdlib + numpy/astropy (both in the container). Writes ``int32``, preserving +the check image's header so the seg map stays on the tile's WCS -- vignetmaker +run 3 cuts the stamps by position, not by row. +""" + +import argparse +import sys + +import numpy as np +from astropy.io import fits + +# The label every unclaimed footprint carries. Negative so that no catalogue +# NUMBER (1..N) can ever collide with it; UberSeg tests `seg != object_number`, +# so one shared label for all neighbours is enough. +NEIGHBOUR_LABEL = -1 + +# Radius, in pixels, of the disc given to an object with no footprint of its +# own. ~0.56" at the CFIS pixel scale: a plausible minimum galaxy core, small +# enough not to take a blend away from the object that owns its footprint. +FALLBACK_RADIUS = 3 + + +def relabel(seg, number, x_image, y_image, fallback_radius=FALLBACK_RADIUS): + """Relabel ``seg`` into the numbering of a catalogue. + + Parameters + ---------- + seg : numpy.ndarray + Segmentation map, 0 for sky and one positive label per detection. + number : numpy.ndarray + The catalogue's ``NUMBER`` column. + x_image, y_image : numpy.ndarray + The catalogue's ``XWIN_IMAGE``/``YWIN_IMAGE`` columns, FITS 1-indexed + pixel coordinates on the same grid as ``seg``. + fallback_radius : int, optional + Radius of the disc painted for an object with no footprint of its own. + + Returns + ------- + (numpy.ndarray, dict) + The relabelled map (``int32``) and a count of each outcome: + the four claim outcomes, which partition the catalogue -- + ``matched``, ``unclaimed`` (fell on sky), ``shared`` (fell on a + footprint already claimed) and ``off_image`` -- plus + ``shared_pixel``, a separate axis counting objects that round to the + same pixel as a lower-numbered one. + """ + number = np.asarray(number) + # FITS pixel centres are 1-based; round to the pixel the centroid is in. + col = np.rint(np.asarray(x_image, dtype=float)).astype(np.int64) - 1 + row = np.rint(np.asarray(y_image, dtype=float)).astype(np.int64) - 1 + n_row, n_col = seg.shape + inside = (col >= 0) & (col < n_col) & (row >= 0) & (row < n_row) + + # Claim, in NUMBER order, the footprint each object's centre falls in. + owner = {} + counts = dict(matched=0, unclaimed=0, shared=0, off_image=0, + shared_pixel=0) + fallback = [] + for i in np.argsort(number, kind="stable"): + if not inside[i]: + counts["off_image"] += 1 + continue + label = int(seg[row[i], col[i]]) + if label == 0: + counts["unclaimed"] += 1 + fallback.append(i) + elif label in owner: + counts["shared"] += 1 + fallback.append(i) + else: + owner[label] = int(number[i]) + counts["matched"] += 1 + + # Every footprint becomes a neighbour unless an object claimed it. A lookup + # table over the labels present is O(pixels) once, rather than one pass per + # object. + labels = np.unique(seg) + table = np.full(int(labels.max()) + 1, NEIGHBOUR_LABEL, dtype=np.int32) + table[0] = 0 + for label, num in owner.items(): + table[label] = num + out = table[np.clip(seg, 0, None)] + + # The fallback discs go on before the centre pass below, so an object that + # shares a claimed footprint takes its centre back from the claimant. + if fallback: + offsets = np.argwhere( + np.add.outer( + np.arange(-fallback_radius, fallback_radius + 1) ** 2, + np.arange(-fallback_radius, fallback_radius + 1) ** 2, + ) <= fallback_radius ** 2 + ) - fallback_radius + for i in fallback: + rows = np.clip(row[i] + offsets[:, 0], 0, n_row - 1) + cols = np.clip(col[i] + offsets[:, 1], 0, n_col - 1) + out[rows, cols] = int(number[i]) + + # EVERY OBJECT ENDS WITH ITS CENTRE PIXEL, and this last pass is what makes + # that true rather than usually true. A disc painted for one object can lie + # over another's pixels -- over a small footprint, or over an earlier disc + # -- and an object with no pixel of its own is not a soft failure: UberSeg + # would hand it an all-zero weight and Ngmix._check_central_seg_label + # raises on a stamp that does not contain the label. Re-stamping the centre + # costs the overlapping object one pixel and nothing else. + # + # Two catalogue objects rounding to the SAME pixel is the one contest a + # pixel cannot settle twice: the lower NUMBER keeps it and the other is + # counted as ``shared_pixel``. The loser still holds the rest of its disc, + # so it keeps a self of its own; the count is there because two objects + # that close are worth seeing in the stage log. + claimed_pixel = {} + for i in np.argsort(number, kind="stable"): + if not inside[i]: + continue + pixel = (int(row[i]), int(col[i])) + if pixel in claimed_pixel: + counts["shared_pixel"] += 1 + continue + claimed_pixel[pixel] = True + out[row[i], col[i]] = int(number[i]) + + return out, counts + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--sexcat", required=True, + help="the tile catalogue (FITS-LDAC) being measured") + parser.add_argument("--segmentation", required=True, + help="the SExtractor SEGMENTATION check image") + parser.add_argument("--output", required=True, + help="where to write the relabelled map") + parser.add_argument("--fallback-radius", type=int, default=FALLBACK_RADIUS) + args = parser.parse_args(argv) + + with fits.open(args.sexcat) as hdus: + cat = hdus["LDAC_OBJECTS"].data + number = np.array(cat["NUMBER"]) + x_image = np.array(cat["XWIN_IMAGE"]) + y_image = np.array(cat["YWIN_IMAGE"]) + with fits.open(args.segmentation) as hdus: + seg = hdus[0].data + header = hdus[0].header + + out, counts = relabel(seg, number, x_image, y_image, args.fallback_radius) + fits.PrimaryHDU(data=out, header=header).writeto(args.output, + overwrite=True) + print(f"seg_relabel: {len(number)} objects, " + ", ".join( + f"{k}={v}" for k, v in counts.items())) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) From 57b75258655289c7349a40132f7c025351e4ca9b Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 17 Sep 2026 15:05:04 +0200 Subject: [PATCH 06/21] Default tile_detection to the UNIONS catalogue The per-tile catalogue is the sample; SExtractor detection is the explicit alternative. `inputs.catalogues` travels with it -- the Snakefile refuses unions_catalogue without one -- and is committed as TBD for the run to set. Co-Authored-By: Claude Fable 5.1 --- workflow/README.md | 11 +++++++---- workflow/config.yaml | 26 +++++++++++++++++--------- 2 files changed, 24 insertions(+), 13 deletions(-) diff --git a/workflow/README.md b/workflow/README.md index aff3cee0a..282c35734 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -26,10 +26,13 @@ uv pip install 'snakemake>=9,<10' 'snakemake-executor-plugin-slurm>=2.7,<3' # outputs.products_dir/index_db, and container. # `psf_model` is `psfex` or `mccd`; mccd is wired but unvalidated here, while psfex is exercised by smk-g4 through smk-g6. -# `tile_detection` is `sextractor` (the tile is detected with SExtractor) or -# `unions_catalogue` (the UNIONS per-tile catalogue at `inputs.catalogues` is -# fetched and converted in place, and its object ID rides into the final -# catalogue as TILE_UNIQUE_ID; no segmentation map, so not with ngmix's uberseg). +# `tile_detection` is `unions_catalogue` (the default: the UNIONS per-tile +# catalogue at `inputs.catalogues` is fetched and converted in place, and its +# object ID rides into the final catalogue as TILE_UNIQUE_ID) or `sextractor` +# (the tile is detected with SExtractor). +# `blend_handling` is ngmix's neighbour treatment, `noisefill` or `uberseg`. +# With `unions_catalogue`, `uberseg` adds one SExtractor run per tile for the +# segmentation map alone, relabelled into the catalogue's numbering. # The committed launcher loads apptainer/1.4.5 + the /project venv, so a # fresh shell always has the right state. diff --git a/workflow/config.yaml b/workflow/config.yaml index 8c3a0e58e..5671c2650 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -24,17 +24,25 @@ inputs: # Pre-staged P3 data (get_images RETRIEVE=symlink). tiles: /project/def-mjhudson/unions-wl/tiles exposures: /project/def-mjhudson/unions-wl/exposures - # The UNIONS per-tile object catalogues (CFIS..r.cat), read only when + # The UNIONS per-tile object catalogues (CFIS..r.cat), read when # tile_detection is unions_catalogue: a local directory (symlinked in) or a - # vos: URL (downloaded, e.g. vos:cfis/tiles_DR6). - # catalogues: vos:cfis/tiles_DR6 + # vos: URL (downloaded, e.g. vos:cfis/tiles_DR6). TBD: set this to the + # catalogue store for the run. + catalogues: TBD -# Where the tile's galaxy sample comes from: `sextractor` runs SExtractor on -# the tile image; `unions_catalogue` fetches the UNIONS per-tile catalogue -# (inputs.catalogues) and converts it in place, carrying its unique object ID -# into the final catalogue as TILE_UNIQUE_ID. The catalogue path has no -# segmentation map, so ngmix's BLEND_HANDLING = uberseg needs sextractor. -tile_detection: sextractor +# Where the tile's galaxy sample comes from. `unions_catalogue` fetches the +# UNIONS per-tile catalogue (inputs.catalogues) and converts it in place, +# carrying its unique object ID into the final catalogue as TILE_UNIQUE_ID; +# `sextractor` runs SExtractor on the tile image and takes its detections. +tile_detection: unions_catalogue + +# How ngmix treats a stamp's neighbours (its BLEND_HANDLING option): +# `noisefill` replaces neighbour pixels with a noise realisation, `uberseg` +# hard-masks them using the tile segmentation map. Under +# tile_detection: unions_catalogue, `uberseg` adds the tile_segmentation rule +# -- one SExtractor run on the tile image for the segmentation map alone, +# relabelled into the catalogue's numbering (workflow/rules/tile.smk). +blend_handling: noisefill # The container every job runs inside (apptainer software-deployment in the profile). container: /project/def-mjhudson/cdaley/containers/shapepipe-develop-runtime.sif From 36e96b3906a3f962287822acbe0897c811aad328 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 28 Sep 2026 16:11:45 +0200 Subject: [PATCH 07/21] Build TILE_UNIQUE_ID once, in make_cat, from one helper cfis.get_tile_unique_id defines the survey-wide object ID, tile_id * 10**6 + NUMBER with tile_id = RRR * 1000 + DDD, and raises when NUMBER (or tile_id) falls outside [0, 10**6). cfis.get_tile_id parses the tile, split_tile_unique_id inverts the encoding, and get_tile_number now refuses tile components with more than three digits. make_cat writes TILE_UNIQUE_ID next to TILE_ID for every final catalogue, whatever produced the detection catalogue. read_ext_sexcat is a plain format converter again: it copies NUMBER unchanged and adds no ID. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Uemfjv9ybCwtKtprksZbVY --- .../modules/make_cat_package/__init__.py | 5 +- .../modules/make_cat_package/make_cat.py | 25 ++--- .../read_ext_sexcat.py | 62 +------------ .../modules/read_ext_sexcat_runner.py | 1 - src/shapepipe/utilities/cfis.py | 92 ++++++++++++++++++- tests/module/test_make_cat.py | 61 ++++++++++++ tests/unit/test_tile_unique_id.py | 49 ++++++++++ 7 files changed, 222 insertions(+), 73 deletions(-) create mode 100644 tests/unit/test_tile_unique_id.py diff --git a/src/shapepipe/modules/make_cat_package/__init__.py b/src/shapepipe/modules/make_cat_package/__init__.py index 838a91c4d..3b37dccc7 100644 --- a/src/shapepipe/modules/make_cat_package/__init__.py +++ b/src/shapepipe/modules/make_cat_package/__init__.py @@ -23,7 +23,10 @@ galaxies for weak-lensing post-processing. This includes galaxy detection and basic measurement parameters, the PSF model at galaxy positions, the spread-model classification, and the shape measurement. 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 f9f3b69d3..7c21dc966 100644 --- a/src/shapepipe/modules/make_cat_package/make_cat.py +++ b/src/shapepipe/modules/make_cat_package/make_cat.py @@ -16,7 +16,7 @@ from sqlitedict import SqliteDict 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): @@ -95,7 +95,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 ---------- @@ -117,22 +121,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/read_ext_sexcat_package/read_ext_sexcat.py b/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py index f1224d791..3bbccdcbd 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 @@ -94,15 +94,13 @@ def _extract_vignets(image_data, x_pos, y_pos, stamp_size): 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 +125,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,46 +156,10 @@ 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, w_log=None, ): @@ -214,12 +172,9 @@ 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``. Parameters ---------- @@ -229,8 +184,6 @@ 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 w_log : logging.Logger, optional @@ -242,11 +195,6 @@ 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) @@ -263,7 +211,7 @@ def make_ldac_from_ascii( ) 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..5c6e7b1db 100644 --- a/src/shapepipe/modules/read_ext_sexcat_runner.py +++ b/src/shapepipe/modules/read_ext_sexcat_runner.py @@ -55,7 +55,6 @@ def read_ext_sexcat_runner( cat_path, image_path, output_path, - file_number_string, stamp_size=stamp_size, w_log=w_log, ) 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/test_make_cat.py b/tests/module/test_make_cat.py index 576ad2c95..4fd9c9078 100644 --- a/tests/module/test_make_cat.py +++ b/tests/module/test_make_cat.py @@ -20,8 +20,10 @@ from astropy.io import fits from sqlitedict import SqliteDict +from shapepipe.modules.make_cat_package import make_cat from shapepipe.modules.make_cat_package.make_cat import SaveCatalogue from shapepipe.modules.ngmix_package.ngmix import Ngmix +from shapepipe.utilities import cfis class _NullLogger: @@ -670,3 +672,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/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) From 4e76ccc3c6ffebf5e9b4b408e4afb0651ce57944 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 28 Sep 2026 16:11:45 +0200 Subject: [PATCH 08/21] Select ngmix chunks by catalogue row, not NUMBER value ID_OBJ_MIN/MAX are now 1-based closed row ranges (ngmix.chunk_rows), so the partition ngmix_range.py writes covers every object once however NUMBER is ordered or spaced. ngmix_range emits NGMIX_ROW_MIN/MAX, and checks that the EPOCH extensions are row-aligned instead of requiring NUMBER = 1..N. NUMBER stays the object's identity: vignet/PSF store keys, the uberseg label check and position seeding are unchanged. fake_psf keys its PSF store by NUMBER rather than row + 1. For SExtractor catalogues (NUMBER = row) chunk contents are unchanged, but ngmix_range.py's hash changes, so a resume across this commit re-measures. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Uemfjv9ybCwtKtprksZbVY --- .../modules/fake_psf_package/fake_psf.py | 5 +- .../modules/ngmix_package/__init__.py | 13 ++- src/shapepipe/modules/ngmix_package/ngmix.py | 47 ++++++-- src/shapepipe/modules/ngmix_runner.py | 8 +- tests/module/test_ngmix_chunk_rows.py | 110 ++++++++++++++++++ tests/unit/test_ngmix_range.py | 30 ++--- workflow/Snakefile | 2 +- .../config/cfis/config_tile_Ng_template.ini | 11 +- workflow/rules/tile.smk | 2 +- workflow/scripts/ngmix_range.py | 52 +++++---- 10 files changed, 211 insertions(+), 69 deletions(-) create mode 100644 tests/module/test_ngmix_chunk_rows.py 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/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 b2e092240..147970798 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -125,13 +125,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):: @@ -414,11 +441,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 @@ -987,11 +1014,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 27648947f..0275ba357 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -109,10 +109,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/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/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/workflow/Snakefile b/workflow/Snakefile index 2befdb5ec..9abeb59f8 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -369,7 +369,7 @@ CLEAN_HASH = script_hash("clean_exposure.py") CLEAN_TILE_HASH = script_hash("clean_tile.py") # 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 diff --git a/workflow/config/cfis/config_tile_Ng_template.ini b/workflow/config/cfis/config_tile_Ng_template.ini index 04606636d..640ab5e41 100644 --- a/workflow/config/cfis/config_tile_Ng_template.ini +++ b/workflow/config/cfis/config_tile_Ng_template.ini @@ -117,9 +117,8 @@ MAG_ZP = 30.0 # Pixel scale in arcsec 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/workflow/rules/tile.smk b/workflow/rules/tile.smk index ccbc83fbf..31bbc022c 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -535,7 +535,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 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__": From e34b27598d0bce34c8c9a0a5bf9d5994ae618be4 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 28 Sep 2026 16:30:23 +0200 Subject: [PATCH 09/21] Segmentation run uses the tile detection's MegaPipe filter config_tile_Sg.ini takes DOT_CONV_FILE from config_tile_Sx.ini (#896), so the uberseg footprints are the ones SExtractor mode would draw. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Uemfjv9ybCwtKtprksZbVY --- workflow/config/cfis/config_tile_Sg.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/config/cfis/config_tile_Sg.ini b/workflow/config/cfis/config_tile_Sg.ini index b02592fd5..9ad3931c3 100644 --- a/workflow/config/cfis/config_tile_Sg.ini +++ b/workflow/config/cfis/config_tile_Sg.ini @@ -88,7 +88,7 @@ EXEC_PATH = source-extractor # SExtractor configuration files DOT_SEX_FILE = $SP_CONFIG/default_tile.sex DOT_PARAM_FILE = $SP_CONFIG/default_noimaflags.param -DOT_CONV_FILE = $SP_CONFIG/default.conv +DOT_CONV_FILE = $SP_CONFIG/gauss_3.0_7x7.conv # Use input weight image if True WEIGHT_IMAGE = True From 2e0902791d360de4265bff532766a8fbf613800e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 28 Sep 2026 17:18:27 +0200 Subject: [PATCH 10/21] Move UberSeg segmentation to its own branch The tile_segmentation rule, seg_relabel.py, config_tile_Sg.ini and the blend_handling setting move to feat/uberseg-seg-map, which feeds the consolidated blend/defect PR. This branch keeps the Gwyn catalogue path, TILE_UNIQUE_ID and row-based ngmix chunking; the catalogue path has no segmentation map, so uberseg needs the SExtractor mode. Co-Authored-By: Claude Opus 5.5 --- tests/unit/test_seg_relabel.py | 145 ---------------- tests/unit/test_workflow_tile_detection.py | 66 ------- workflow/README.md | 3 - workflow/Snakefile | 21 +-- workflow/bin/sp | 1 - workflow/config.yaml | 8 - workflow/config/cfis/config_tile_Sg.ini | 131 -------------- workflow/rules/tile.smk | 72 +------- workflow/scripts/completeness.py | 11 -- workflow/scripts/run_report.py | 6 - workflow/scripts/seg_relabel.py | 191 --------------------- 11 files changed, 3 insertions(+), 652 deletions(-) delete mode 100644 tests/unit/test_seg_relabel.py delete mode 100644 workflow/config/cfis/config_tile_Sg.ini delete mode 100644 workflow/scripts/seg_relabel.py diff --git a/tests/unit/test_seg_relabel.py b/tests/unit/test_seg_relabel.py deleted file mode 100644 index 56d9600d1..000000000 --- a/tests/unit/test_seg_relabel.py +++ /dev/null @@ -1,145 +0,0 @@ -"""The segmentation map UberSeg sees must call THIS object self. - -``uberseg_weight`` keeps the pixels whose nearest footprint carries the -object's catalogue ``NUMBER`` and zeros the rest; it never looks at the label -under the object's position. So the one thing that has to hold, for every -object in the UNIONS catalogue, is that the relabelled map carries that -object's ``NUMBER`` on its own pixels and something else everywhere else. -These tests assert exactly that, on maps small enough to read by eye — -SExtractor is not involved. -""" - -import importlib.util -import sys -from pathlib import Path - -import numpy as np -import pytest - -SCRIPTS = Path(__file__).resolve().parents[2] / "workflow" / "scripts" - - -def _load(): - sys.path.insert(0, str(SCRIPTS)) - try: - spec = importlib.util.spec_from_file_location( - "_seg_relabel", SCRIPTS / "seg_relabel.py") - module = importlib.util.module_from_spec(spec) - spec.loader.exec_module(module) - finally: - sys.path.remove(str(SCRIPTS)) - return module - - -seg_relabel = _load() - - -def _map(): - """Two SExtractor 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 - - -def test_each_object_owns_the_footprint_it_sits_in(): - """The claim is by position; the label that comes out is the NUMBER.""" - seg = _map() - # FITS 1-indexed centres of the two footprints, in an order that is NOT - # the label order, so a relabelling that merely renumbered would fail. - out, counts = seg_relabel.relabel( - 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(): - """A detection no catalogue object sits in is masked, not measured.""" - seg = _map() - out, counts = seg_relabel.relabel( - seg, number=np.array([1]), - x_image=np.array([4.0]), y_image=np.array([4.0])) - assert counts["matched"] == 1 - assert set(np.unique(out[seg == 9])) == {seg_relabel.NEIGHBOUR_LABEL} - assert seg_relabel.NEIGHBOUR_LABEL not in np.unique(out[seg == 7]) - - -def test_an_object_on_sky_still_gets_a_self(): - """No footprint at its position -> a disc of its own NUMBER. - - Without one, uberseg_weight would zero the whole stamp and - Ngmix._check_central_seg_label would raise on a map lacking the label. - """ - seg = _map() - out, counts = seg_relabel.relabel( - 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 - # Centred on the object, and nowhere else. - 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(): - """The first claimant keeps the blend; the second takes back its centre.""" - seg = _map() - out, counts = seg_relabel.relabel( - 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 - assert out[15, 15] == 2 - # The claimant still holds the bulk of the footprint. - assert np.count_nonzero(out == 1) > np.count_nonzero(out == 2) - - -def test_every_object_is_self_somewhere(): - """The invariant, over a randomised map: every object on the image keeps - pixels of its own — matched, unmatched, blended or doubled.""" - 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 = seg_relabel.relabel(seg, number, x_image, y_image) - # The four claim outcomes partition the catalogue; shared_pixel is a - # separate axis, counted on top. - assert sum(counts[k] for k in - ("matched", "unclaimed", "shared", "off_image")) == len(number) - assert counts["off_image"] == 0 - for num in number: - assert np.any(out == num), f"object {num} has no self pixels" - # And nothing outside a footprint or a disc was invented. - assert set(np.unique(out)) <= set(number) | {0, seg_relabel.NEIGHBOUR_LABEL} - - -def test_two_objects_on_one_pixel_split_the_contest(): - """A pixel has one label: the lower NUMBER keeps it, the other keeps its - disc, and the collision is counted rather than hidden.""" - seg = _map() - out, counts = seg_relabel.relabel( - 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, 61.0)]) -def test_a_position_off_the_image_is_counted_not_crashed(x, y): - seg = _map() - out, counts = seg_relabel.relabel( - seg, 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) diff --git a/tests/unit/test_workflow_tile_detection.py b/tests/unit/test_workflow_tile_detection.py index 98565a6c7..c2d5767a9 100644 --- a/tests/unit/test_workflow_tile_detection.py +++ b/tests/unit/test_workflow_tile_detection.py @@ -12,7 +12,6 @@ import configparser import importlib.util -import re import sys from pathlib import Path @@ -119,70 +118,6 @@ def test_report_lists_the_fetch_stage_only_for_catalogue_runs(monkeypatch, mode, assert stages.index("tile_get_catalogue") == stages.index("tile_detect") - 1 -# --- the segmentation half: unions_catalogue + uberseg ---------------------- - - -def test_blend_handlings_mirror_the_ngmix_module(): - """completeness.BLEND_HANDLINGS is a copy; the copy must stay true. - - The Snakefile validates `blend_handling:` against it in the launcher venv, - outside the container where shapepipe is importable, so the tuple is - mirrored rather than imported. Read out of the source text here for the - same reason: this file is container-free. - """ - src = (REPO_ROOT / "src" / "shapepipe" / "modules" / "ngmix_package" - / "ngmix.py").read_text() - match = re.search(r"^BLEND_HANDLINGS = \(([^)]*)\)", src, re.M) - assert match, "ngmix.py no longer defines BLEND_HANDLINGS" - assert completeness.BLEND_HANDLINGS == tuple( - part.strip().strip('"\'') for part in match.group(1).split(",") - if part.strip()) - - -def test_segmentation_ini_matches_its_stage_and_yields_the_map(): - level, sg = completeness.STAGE_DIR["tile_segmentation"] - assert level == "tile" - ini = _ini("config_tile_Sg.ini") - assert ini["DEFAULT"]["RUN_NAME"].strip() == sg - assert ini["DEFAULT"]["RUN_DATETIME"].strip() == "False" - sx = ini["SEXTRACTOR_RUNNER"] - # The check image is the whole point, and it is the only one asked for. - assert [c.strip() for c in sx["CHECKIMAGE"].split(",")] == ["SEGMENTATION"] - # No multi-epoch post-processing: those extensions ride on the sexcat - # config_tile_Uc.ini writes, and without it the WCS sqlite is not an input. - assert sx["MAKE_POST_PROCESS"].strip() == "False" - assert "log_exp_headers" not in sx["FILE_PATTERN"] - # Same detection settings as the SExtractor mode, so the footprints are - # the ones that mode would have drawn. - detect = _ini("config_tile_Sx.ini")["SEXTRACTOR_RUNNER"] - for key in ("DOT_SEX_FILE", "DOT_PARAM_FILE", "DOT_CONV_FILE", - "WEIGHT_IMAGE", "FLAG_IMAGE"): - assert sx[key].strip() == detect[key].strip(), key - - -def test_segmentation_stage_is_checked(tmp_path): - """Two products: the SEGMENTATION check image and SExtractor's own sexcat.""" - ok, details = completeness.check_counts( - "tile_segmentation", _stage_dir(tmp_path, "sextractor_runner", 2)) - assert ok and details == [("sextractor_runner", 2, 2, False)] - - -@pytest.mark.parametrize("detection, blend, present", [ - ("unions_catalogue", "uberseg", True), - ("unions_catalogue", "noisefill", False), - ("sextractor", "uberseg", False), -]) -def test_report_lists_the_segmentation_stage_only_for_the_pair( - monkeypatch, detection, blend, present): - monkeypatch.setenv("SP_TILE_DETECTION", detection) - monkeypatch.setenv("SP_BLEND_HANDLING", blend) - stages = _load("run_report").TILE_STAGES - assert ("tile_segmentation" in stages) is present - if present: - assert stages.index("tile_segmentation") == stages.index( - "tile_detect") + 1 - - @pytest.mark.parametrize("machine", ["nibi", "candide"]) @pytest.mark.parametrize("input_type", ["data", "image_sims"]) def test_machine_defaults_pair_the_catalogue_with_its_source( @@ -200,7 +135,6 @@ def test_machine_defaults_pair_the_catalogue_with_its_source( over.write_text(f"input_type: {input_type}\n") config = run_config.load( str(REPO_ROOT / "workflow" / "config.yaml"), str(over)) - assert config["blend_handling"] in completeness.BLEND_HANDLINGS if input_type == "data": assert config["tile_detection"] == "unions_catalogue" assert config["inputs"]["catalogues"] diff --git a/workflow/README.md b/workflow/README.md index 31137bd3d..2c4a4a181 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -32,9 +32,6 @@ uv pip install 'snakemake>=9,<10' 'snakemake-executor-plugin-slurm>=2.7,<3' # fetched and converted in place, keeping its NUMBER) 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. -# `blend_handling` is ngmix's neighbour treatment, `noisefill` or `uberseg`. -# With `unions_catalogue`, `uberseg` adds one SExtractor run per tile for the -# segmentation map alone, relabelled into the catalogue's numbering. # The committed launcher loads apptainer/1.4.5 + the /project venv, so a # fresh shell always has the right state. diff --git a/workflow/Snakefile b/workflow/Snakefile index 0c8a1039d..4cafdf88a 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -153,8 +153,7 @@ CONFIG_DIR = Path(workflow.basedir) / "config" / { sys.path.insert(0, str(SCRIPTS)) import build_index # noqa: E402 -from completeness import (BLEND_HANDLINGS, STAGE_DIR, # noqa: E402 - TILE_DETECTIONS) +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). @@ -175,20 +174,6 @@ if TILE_DETECTION == "unions_catalogue" and not CATALOGUES: "directory or vos: URL holding the per-tile CFIS..r.cat files).") CATALOGUE_RETRIEVE = "vos" if CATALOGUES.startswith("vos:") else "symlink" -# How ngmix treats a stamp's neighbours (its BLEND_HANDLING option): the -# noise-fill treatment needs nothing of the workflow, while uberseg needs a -# segmentation map labelled in the numbering of the catalogue being measured. -# SExtractor detection writes that map as a by-product (config_tile_Sx.ini's -# CHECKIMAGE); the UNIONS catalogue has none, so the pair below is what adds -# tile.smk's tile_segmentation rule. -BLEND_HANDLING = config.get("blend_handling", "noisefill") -if BLEND_HANDLING not in BLEND_HANDLINGS: - raise WorkflowError( - f"Invalid blend_handling={BLEND_HANDLING!r}; expected one of " - f"{sorted(BLEND_HANDLINGS)}.") -TILE_SEGMENTATION = (TILE_DETECTION == "unions_catalogue" - and BLEND_HANDLING == "uberseg") - # 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 # --dag`, `sp exp_psf ...`). It gates the two parse-time side effects — the index @@ -423,9 +408,6 @@ CLEAN_TILE_HASH = script_hash("clean_tile.py") # fresh root: a resume across it re-measures the tile, and the failure modes a # mid-campaign params change can reach are catalogued at tile.smk's range_hash. NGMIX_RANGE_HASH = script_hash("ngmix_range.py") -# seg_relabel.py decides which footprint is "self" for every object UberSeg -# masks, so an edit to it changes every mask; it rides on tile_segmentation. -SEG_RELABEL_HASH = script_hash("seg_relabel.py") # --- exposure reclamation (D5, S5) ----------------------------------------- @@ -779,7 +761,6 @@ rule prepare_all_tiles: # via `sp report`. def _report(status): shell(f"SP_TILE_DETECTION='{TILE_DETECTION}' " - f"SP_BLEND_HANDLING='{BLEND_HANDLING}' " f"python {SCRIPTS}/run_report.py --run-dir {RUN_DIR} " f"--index {INDEX_DB} --status {status} || true") diff --git a/workflow/bin/sp b/workflow/bin/sp index 86b669088..bcd5295d6 100755 --- a/workflow/bin/sp +++ b/workflow/bin/sp @@ -338,7 +338,6 @@ case "$cmd" in # either way -- this is consistency, not safety. shift SP_TILE_DETECTION="$(cfg tile_detection)" \ - SP_BLEND_HANDLING="$(cfg blend_handling)" \ 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 b6e7007d7..266413d30 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -24,14 +24,6 @@ input_type: data # Run name, available as `$run` in the paths below. run: smk-g6 -# How ngmix treats a stamp's neighbours (its BLEND_HANDLING option): -# `noisefill` replaces neighbour pixels with a noise realisation, `uberseg` -# hard-masks them using the tile segmentation map. Under -# tile_detection: unions_catalogue, `uberseg` adds the tile_segmentation rule -# -- one SExtractor run on the tile image for the segmentation map alone, -# relabelled into the catalogue's numbering (workflow/rules/tile.smk). -blend_handling: noisefill - # tile_detection is set per (machine, input_type) in the table below: # `unions_catalogue` fetches the UNIONS per-tile catalogue (inputs.catalogues: # a local directory, symlinked in, or a vos: URL) and converts it in place; diff --git a/workflow/config/cfis/config_tile_Sg.ini b/workflow/config/cfis/config_tile_Sg.ini deleted file mode 100644 index 9ad3931c3..000000000 --- a/workflow/config/cfis/config_tile_Sg.ini +++ /dev/null @@ -1,131 +0,0 @@ -# ShapePipe configuration file for the tile SEGMENTATION map, used when -# `tile_detection: unions_catalogue` is combined with `blend_handling: uberseg` -# in the workflow run config. -# -# The UNIONS per-tile catalogue is the galaxy sample; this run exists only for -# the one product that catalogue cannot carry, the segmentation map UberSeg -# needs. Its own sexcat is a by-product and nothing reads it. The detection -# settings are those of config_tile_Sx.ini, so the footprints are the ones a -# SExtractor-detection run would have drawn. -# -# The labels of the SEGMENTATION check image are this run's own NUMBERs, which -# are not the converted catalogue's; workflow/scripts/seg_relabel.py maps them -# into the catalogue's NUMBER space and writes the result beside the sexcat. - - -## Default ShapePipe options -[DEFAULT] - -# verbose mode (optional), default: True, print messages on terminal -VERBOSE = True - -# Name of run (optional) default: shapepipe_run -RUN_NAME = run_sp_tile_Sg - -# Add date and time to RUN_NAME, optional, default: True -RUN_DATETIME = False - - -## ShapePipe execution options -[EXECUTION] - -# Module name, single string or comma-separated list of valid module runner names -MODULE = sextractor_runner - -# Run mode, SMP or MPI -MODE = SMP - - -## ShapePipe file handling options -[FILE] - -# Log file master name, optional, default: shapepipe -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 - -# Output directory -OUTPUT_DIR = $SP_RUN/output - - -## ShapePipe job handling options -[JOB] - -# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial -SMP_BATCH_SIZE = 16 - -# Timeout value (optional), default is None, i.e. no timeout limit applied -TIMEOUT = 96:00:00 - - -## Module options - -[SEXTRACTOR_RUNNER] - -# The tile image and its weight. The merged WCS headers are absent because -# MAKE_POST_PROCESS is False: the multi-epoch extensions belong to the sexcat -# the chain reads, which config_tile_Uc.ini writes. -INPUT_DIR = $SP_RUN/output/run_sp_tile_Git/get_images_runner/output, $SP_RUN/output/run_sp_tile_Uz/uncompress_fits_runner/output - -FILE_PATTERN = CFIS_image, CFIS_weight - -FILE_EXT = .fits, .fits - -# NUMBERING_SCHEME (optional) string with numbering pattern for input files -NUMBERING_SCHEME = -000-000 - -# SExtractor executable path -EXEC_PATH = source-extractor - -# SExtractor configuration files -DOT_SEX_FILE = $SP_CONFIG/default_tile.sex -DOT_PARAM_FILE = $SP_CONFIG/default_noimaflags.param -DOT_CONV_FILE = $SP_CONFIG/gauss_3.0_7x7.conv - -# Use input weight image if True -WEIGHT_IMAGE = True - -# Use input flag image if True -FLAG_IMAGE = False - -# Use input PSF file if True -PSF_FILE = False - -# Use distinct image for detection (SExtractor in -# dual-image mode) if True -DETECTION_IMAGE = False - -# Distinct weight image for detection (SExtractor -# in dual-image mode) -DETECTION_WEIGHT = False - -ZP_FROM_HEADER = False - -BKG_FROM_HEADER = False - -# Type of image check (optional), default not used, can be a list of -# BACKGROUND, BACKGROUND_RMS, INIBACKGROUND, -# MINIBACK_RMS, -BACKGROUND, #FILTERED, -# OBJECTS, -OBJECTS, SEGMENTATION, APERTURES -# -# SEGMENTATION alone: this run's reason to exist. BACKGROUND is the tile -# background config_tile_Sx.ini writes for the detection chain, and nothing -# reads it from here. -CHECKIMAGE = SEGMENTATION - -# File name suffix for the output sextractor files (optional) -SUFFIX = sexcat - -## Post-processing - -# The multi-epoch extensions ride on the sexcat the chain reads, not on this -# run's by-product catalogue. -MAKE_POST_PROCESS = False diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index 9996b4655..e37984682 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -51,14 +51,8 @@ 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. What -Steven's catalogue does not carry is a segmentation map, and ngmix's -``BLEND_HANDLING = uberseg`` needs one. ``blend_handling: uberseg`` therefore -adds a third rule, tile_segmentation: a SExtractor run on the tile image whose -only product is the SEGMENTATION check image (config_tile_Sg.ini), followed by -seg_relabel.py, which maps that image's labels into the converted catalogue's -NUMBER space and writes it beside the sexcat -- the path vignetmaker's -segmentation run reads. SExtractor runs once, for the one thing the catalogue -cannot give. +this path does not have is a segmentation map: ngmix's ``BLEND_HANDLING = +uberseg`` needs the SExtractor mode. 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 @@ -560,67 +554,6 @@ else: '"$SP_RUN/output/run_sp_tile_Sx/sextractor_runner"\n' "fi\n") -# Where the tile catalogue lands, whichever mode wrote it: the one path every -# downstream config reads, and where the relabelled segmentation map goes. -SEXCAT_DIR = "$SP_RUN/output/run_sp_tile_Sx/sextractor_runner/output" - -if TILE_SEGMENTATION: - - # The segmentation map UberSeg needs, and nothing else. - # - # THE SAMPLE IS STEVEN'S CATALOGUE, ALWAYS. This SExtractor run detects on - # the tile image with config_tile_Sx.ini's settings, but its catalogue is - # discarded: what leaves the rule is the SEGMENTATION check image. That - # image's labels are this run's own NUMBERs, which have nothing to do with - # the converted catalogue's, so the post step relabels it -- UberSeg finds - # the central object BY LABEL (uberseg_weight's `object_number` is the - # catalogue NUMBER), so an unmapped map would invert every mask rather than - # fail. seg_relabel.py's docstring owns the mapping and its fallback. - # - # It writes into tile_detect's run dir, which is where vignetmaker's - # segmentation run looks (`run_sp_tile_Sx:sextractor_runner`, FILE_PATTERN - # `segmentation`) whichever mode produced the sexcat. Safe because unit_pre - # clears only THIS stage's run dir (run_sp_tile_Sg), and a tile_detect - # rerun re-schedules this rule through the manifest edge below. - rule tile_segmentation: - input: - sx = rules.tile_detect.output.manifest, - uz = f"{TILE_DIR}/manifests/tile_uncompress.json", - output: - manifest = f"{TILE_DIR}/manifests/tile_segmentation.json" - log: - f"{TILE_DIR}/logs/tile_segmentation.json" - params: - pre = lambda wc: unit_pre("tile_segmentation", wc.tile), - script_hash = SCRIPT_HASH, - relabel_hash = SEG_RELABEL_HASH - threads: 8 - resources: - # The same SExtractor run as tile_detect's sextractor mode. - mem_mb = lambda wc, attempt: 16000 * attempt, - runtime = 180 - shell: - sp_shell( - "tile_segmentation", "config_tile_Sg.ini", - post=( - "if [ $rc -eq 0 ]; then\n" - f" python {SCRIPTS}/seg_relabel.py" - f' --sexcat "{SEXCAT_DIR}/sexcat$SP_UNIT_NUM.fits"' - ' --segmentation "$SP_RUN/output/run_sp_tile_Sg' - '/sextractor_runner/output/segmentation$SP_UNIT_NUM.fits"' - f' --output "{SEXCAT_DIR}/segmentation' - '$SP_UNIT_NUM.fits" || rc=1\n' - "fi\n")) - - -# The segmentation edge: present only when tile_segmentation is defined, so the -# DAG carries it exactly when the map has to be built. -def tile_seg(wc): - if not TILE_SEGMENTATION: - return [] - return [f"{tile_dir(wc.tile)}/manifests/tile_segmentation.json"] - - # Configured PSF interpolation to galaxies + vignet postage stamps: the last # stage that reads exposure products, and the bulk intra-tile intermediate. The store it # writes is node-local (see TILE_LOCAL above). @@ -628,7 +561,6 @@ rule tile_vignets: group: TILE_GROUP input: sx = rules.tile_detect.output.manifest, - seg = tile_seg, forest = rules.tile_exp_forest.output.forest, split = tile_exp_split, psf = tile_exp_psf, diff --git a/workflow/scripts/completeness.py b/workflow/scripts/completeness.py index 826f8f5da..6620ca015 100644 --- a/workflow/scripts/completeness.py +++ b/workflow/scripts/completeness.py @@ -78,13 +78,6 @@ # when it is not the default, so a SExtractor run's prologue is unchanged. TILE_DETECTIONS = ("sextractor", "unions_catalogue") -# ngmix's neighbour treatments, mirroring BLEND_HANDLINGS in -# shapepipe.modules.ngmix_package.ngmix. Mirrored rather than imported: the -# Snakefile parses this module in the launcher venv, outside the container -# where shapepipe lives. tests/unit/test_workflow_tile_detection.py asserts -# the two tuples agree. -BLEND_HANDLINGS = ("noisefill", "uberseg") - # stage -> {runner_subdir: {expect, [warn], [subpath]}} # exp_psf and tile_vignets are selected by $SP_PSF at check time, tile_detect # by $SP_TILE_DETECTION. @@ -151,9 +144,6 @@ "tile_merge_headers": {"merge_headers_runner": dict(expect=1)}, # The fetched UNIONS catalogue: one .cat per tile. "tile_get_catalogue": {"get_images_runner": dict(expect=1)}, - # The segmentation-only SExtractor run: its SEGMENTATION check image, plus - # the sexcat SExtractor always writes and nothing here reads. - "tile_segmentation": {"sextractor_runner": dict(expect=2)}, "tile_detect": { "sextractor": {"sextractor_runner": dict(expect=2)}, # One FITS-LDAC sexcat, converted from the fetched catalogue. @@ -272,7 +262,6 @@ def check_counts(stage, run_dir): "exp_psf": ("exp", "run_sp_exp_SxSePsf"), "tile_merge_headers": ("tile", "run_sp_tile_Mh_exp"), "tile_get_catalogue": ("tile", "run_sp_tile_Gic"), - "tile_segmentation": ("tile", "run_sp_tile_Sg"), # 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"), diff --git a/workflow/scripts/run_report.py b/workflow/scripts/run_report.py index 2b9e82f2a..2462cdae0 100644 --- a/workflow/scripts/run_report.py +++ b/workflow/scripts/run_report.py @@ -65,12 +65,6 @@ # 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") - # The segmentation run is the second half of that pair: it exists only when - # the converted catalogue has to be measured with ngmix's uberseg blend - # handling, which needs a segmentation map the catalogue does not carry. - if os.environ.get("SP_BLEND_HANDLING") == "uberseg": - TILE_STAGES.insert(TILE_STAGES.index("tile_detect") + 1, - "tile_segmentation") EXP_STAGES = ["exp_get_images", "exp_split", "exp_psf"] # The manifests clean_tile leaves on disk (workflow/scripts/clean_tile.py names diff --git a/workflow/scripts/seg_relabel.py b/workflow/scripts/seg_relabel.py deleted file mode 100644 index 7b8b9b127..000000000 --- a/workflow/scripts/seg_relabel.py +++ /dev/null @@ -1,191 +0,0 @@ -"""Map a SExtractor SEGMENTATION map into the tile catalogue's NUMBER space. - -WHY THIS EXISTS. UberSeg identifies the central object of a stamp BY LABEL: it -keeps the pixels whose nearest segmentation footprint carries the object's own -``NUMBER`` and zeros every other pixel's weight -(``ngmix_package/ngmix.py::uberseg_weight``, called with ``object_number = -obj_id``, the catalogue's ``NUMBER``). It never reads the label under the -object's position. So the seg stamp ngmix overlays must be labelled in the -SAME numbering as the catalogue it measures, or every mask is inverted. - -Under ``tile_detection: unions_catalogue`` the two numberings are different by -construction: the sample is Steven Gwyn's per-tile catalogue, whose NUMBER -``read_ext_sexcat`` keeps as is, while the segmentation map comes from a separate -SExtractor run whose labels are ITS OWN detection numbers. This script is the -bridge, and it is the whole of the guarantee: - - * each catalogue object claims the footprint its position falls in -- the one - lookup that is meaningful, since the two runs share the tile's pixel grid; - * that footprint is relabelled with the object's catalogue ``NUMBER``; - * every footprint no catalogue object claims is relabelled ``NEIGHBOUR_LABEL`` - (negative, so it can never collide with a ``NUMBER``): UberSeg's partition - only ever asks "self or not self", so neighbours need no identity; - * an object whose position falls on sky, or on a footprint another object - claimed first, is given a small disc of its own ``NUMBER`` at its position. - -THE FALLBACK IS NOT COSMETIC. Without a footprint of its own an object has no -"self" in the Voronoi partition: ``uberseg_weight`` would hand it an all-zero -weight, and ``Ngmix._check_central_seg_label`` raises on a stamp that does not -contain the object's label at all. The disc gives such an object a defined, -conservative core -- small, centred, and (in the shared-footprint case) taking -its pixels from the neighbour that claimed the blend, which is the right way -round: the claimant keeps the bulk, the unclaimed object keeps a centre. - -Claiming is in ``NUMBER`` order and first come first served, so the mapping is -deterministic and reproducible from the same two inputs. - -Stdlib + numpy/astropy (both in the container). Writes ``int32``, preserving -the check image's header so the seg map stays on the tile's WCS -- vignetmaker -run 3 cuts the stamps by position, not by row. -""" - -import argparse -import sys - -import numpy as np -from astropy.io import fits - -# The label every unclaimed footprint carries. Negative so that no catalogue -# NUMBER (positive) can ever collide with it; UberSeg tests `seg != object_number`, -# so one shared label for all neighbours is enough. -NEIGHBOUR_LABEL = -1 - -# Radius, in pixels, of the disc given to an object with no footprint of its -# own. ~0.56" at the CFIS pixel scale: a plausible minimum galaxy core, small -# enough not to take a blend away from the object that owns its footprint. -FALLBACK_RADIUS = 3 - - -def relabel(seg, number, x_image, y_image, fallback_radius=FALLBACK_RADIUS): - """Relabel ``seg`` into the numbering of a catalogue. - - Parameters - ---------- - seg : numpy.ndarray - Segmentation map, 0 for sky and one positive label per detection. - number : numpy.ndarray - The catalogue's ``NUMBER`` column. - x_image, y_image : numpy.ndarray - The catalogue's ``XWIN_IMAGE``/``YWIN_IMAGE`` columns, FITS 1-indexed - pixel coordinates on the same grid as ``seg``. - fallback_radius : int, optional - Radius of the disc painted for an object with no footprint of its own. - - Returns - ------- - (numpy.ndarray, dict) - The relabelled map (``int32``) and a count of each outcome: - the four claim outcomes, which partition the catalogue -- - ``matched``, ``unclaimed`` (fell on sky), ``shared`` (fell on a - footprint already claimed) and ``off_image`` -- plus - ``shared_pixel``, a separate axis counting objects that round to the - same pixel as a lower-numbered one. - """ - number = np.asarray(number) - # FITS pixel centres are 1-based; round to the pixel the centroid is in. - col = np.rint(np.asarray(x_image, dtype=float)).astype(np.int64) - 1 - row = np.rint(np.asarray(y_image, dtype=float)).astype(np.int64) - 1 - n_row, n_col = seg.shape - inside = (col >= 0) & (col < n_col) & (row >= 0) & (row < n_row) - - # Claim, in NUMBER order, the footprint each object's centre falls in. - owner = {} - counts = dict(matched=0, unclaimed=0, shared=0, off_image=0, - shared_pixel=0) - fallback = [] - for i in np.argsort(number, kind="stable"): - if not inside[i]: - counts["off_image"] += 1 - continue - label = int(seg[row[i], col[i]]) - if label == 0: - counts["unclaimed"] += 1 - fallback.append(i) - elif label in owner: - counts["shared"] += 1 - fallback.append(i) - else: - owner[label] = int(number[i]) - counts["matched"] += 1 - - # Every footprint becomes a neighbour unless an object claimed it. A lookup - # table over the labels present is O(pixels) once, rather than one pass per - # object. - labels = np.unique(seg) - table = np.full(int(labels.max()) + 1, NEIGHBOUR_LABEL, dtype=np.int32) - table[0] = 0 - for label, num in owner.items(): - table[label] = num - out = table[np.clip(seg, 0, None)] - - # The fallback discs go on before the centre pass below, so an object that - # shares a claimed footprint takes its centre back from the claimant. - if fallback: - offsets = np.argwhere( - np.add.outer( - np.arange(-fallback_radius, fallback_radius + 1) ** 2, - np.arange(-fallback_radius, fallback_radius + 1) ** 2, - ) <= fallback_radius ** 2 - ) - fallback_radius - for i in fallback: - rows = np.clip(row[i] + offsets[:, 0], 0, n_row - 1) - cols = np.clip(col[i] + offsets[:, 1], 0, n_col - 1) - out[rows, cols] = int(number[i]) - - # EVERY OBJECT ENDS WITH ITS CENTRE PIXEL, and this last pass is what makes - # that true rather than usually true. A disc painted for one object can lie - # over another's pixels -- over a small footprint, or over an earlier disc - # -- and an object with no pixel of its own is not a soft failure: UberSeg - # would hand it an all-zero weight and Ngmix._check_central_seg_label - # raises on a stamp that does not contain the label. Re-stamping the centre - # costs the overlapping object one pixel and nothing else. - # - # Two catalogue objects rounding to the SAME pixel is the one contest a - # pixel cannot settle twice: the lower NUMBER keeps it and the other is - # counted as ``shared_pixel``. The loser still holds the rest of its disc, - # so it keeps a self of its own; the count is there because two objects - # that close are worth seeing in the stage log. - claimed_pixel = {} - for i in np.argsort(number, kind="stable"): - if not inside[i]: - continue - pixel = (int(row[i]), int(col[i])) - if pixel in claimed_pixel: - counts["shared_pixel"] += 1 - continue - claimed_pixel[pixel] = True - out[row[i], col[i]] = int(number[i]) - - return out, counts - - -def main(argv=None): - parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) - parser.add_argument("--sexcat", required=True, - help="the tile catalogue (FITS-LDAC) being measured") - parser.add_argument("--segmentation", required=True, - help="the SExtractor SEGMENTATION check image") - parser.add_argument("--output", required=True, - help="where to write the relabelled map") - parser.add_argument("--fallback-radius", type=int, default=FALLBACK_RADIUS) - args = parser.parse_args(argv) - - with fits.open(args.sexcat) as hdus: - cat = hdus["LDAC_OBJECTS"].data - number = np.array(cat["NUMBER"]) - x_image = np.array(cat["XWIN_IMAGE"]) - y_image = np.array(cat["YWIN_IMAGE"]) - with fits.open(args.segmentation) as hdus: - seg = hdus[0].data - header = hdus[0].header - - out, counts = relabel(seg, number, x_image, y_image, args.fallback_radius) - fits.PrimaryHDU(data=out, header=header).writeto(args.output, - overwrite=True) - print(f"seg_relabel: {len(number)} objects, " + ", ".join( - f"{k}={v}" for k, v in counts.items())) - return 0 - - -if __name__ == "__main__": - sys.exit(main()) From 5511597550c3047c3faf58be0ddca4fbcdb3b231 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 28 Sep 2026 17:30:36 +0200 Subject: [PATCH 11/21] Default tile_detection and psf_model per input type, not per machine Both follow from the kind of input (the UNIONS catalogue and PSFEx for data, SExtractor and the true simulation PSF for image sims), so they move out of the machine entries into a top-level input_types: table. run_config layers config.yaml < input_types[input_type] < machines[machine][input_type] < run config, and expands $variables in every resolved key. Co-Authored-By: Claude Opus 5.5 --- astra.yaml | 5 +- scripts/python/create_final_cat.py | 2 +- tests/unit/test_run_config.py | 107 +++++++++++++++++++++ tests/unit/test_workflow_tile_detection.py | 12 ++- workflow/README.md | 18 ++-- workflow/Snakefile | 5 +- workflow/bin/sp | 2 +- workflow/config.yaml | 46 +++++---- workflow/scripts/run_config.py | 59 ++++++------ 9 files changed, 182 insertions(+), 74 deletions(-) create mode 100644 tests/unit/test_run_config.py diff --git a/astra.yaml b/astra.yaml index 87c844909..e6902a189 100644 --- a/astra.yaml +++ b/astra.yaml @@ -914,9 +914,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 f9f10564e..387fc9383 100755 --- a/scripts/python/create_final_cat.py +++ b/scripts/python/create_final_cat.py @@ -33,7 +33,7 @@ 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 + 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 diff --git a/tests/unit/test_run_config.py b/tests/unit/test_run_config.py new file mode 100644 index 000000000..8e24a8c15 --- /dev/null +++ b/tests/unit/test_run_config.py @@ -0,0 +1,107 @@ +"""The run-config resolver's layering: config.yaml < input_types < machines < +run config, with $variables expanded across every layer. + +Container-free: run_config.py needs only PyYAML. +""" + +import importlib.util +from pathlib import Path + +import pytest + +yaml = pytest.importorskip("yaml") + +REPO_ROOT = Path(__file__).resolve().parents[2] +CONFIG_YAML = REPO_ROOT / "workflow" / "config.yaml" + +_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) diff --git a/tests/unit/test_workflow_tile_detection.py b/tests/unit/test_workflow_tile_detection.py index c2d5767a9..95a174023 100644 --- a/tests/unit/test_workflow_tile_detection.py +++ b/tests/unit/test_workflow_tile_detection.py @@ -120,13 +120,15 @@ def test_report_lists_the_fetch_stage_only_for_catalogue_runs(monkeypatch, mode, @pytest.mark.parametrize("machine", ["nibi", "candide"]) @pytest.mark.parametrize("input_type", ["data", "image_sims"]) -def test_machine_defaults_pair_the_catalogue_with_its_source( +def test_defaults_pair_the_catalogue_with_its_source( machine, input_type, monkeypatch, tmp_path): - """Real data defaults to the catalogue and declares where it lives. + """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`, so the two settings travel together in the machines - table. Image sims have no UNIONS catalogue and fall back to SExtractor. + `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") @@ -139,4 +141,4 @@ def test_machine_defaults_pair_the_catalogue_with_its_source( assert config["tile_detection"] == "unions_catalogue" assert config["inputs"]["catalogues"] else: - assert config.get("tile_detection", "sextractor") == "sextractor" + assert config["tile_detection"] == "sextractor" diff --git a/workflow/README.md b/workflow/README.md index 2c4a4a181..08213697d 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -27,10 +27,10 @@ uv pip install 'snakemake>=9,<10' 'snakemake-executor-plugin-slurm>=2.7,<3' # `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 machines-table default for -# input_type data: the UNIONS per-tile catalogue at `inputs.catalogues` is -# fetched and converted in place, keeping its NUMBER) or `sextractor` (the tile -# is detected with SExtractor; the default for image sims). Either way make_cat +# `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) 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. # The committed launcher loads apptainer/1.4.5 + the /project venv, so a @@ -80,9 +80,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`, `$run` to the run config's `run:`). A value of `TBD` stops the run at parse time until it is set. A run config therefore only needs what differs, e.g. for one SKiLLS shear branch on candide: @@ -90,8 +92,6 @@ SKiLLS shear branch on candide: ```yaml machine: candide input_type: image_sims -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 diff --git a/workflow/Snakefile b/workflow/Snakefile index 4cafdf88a..76845e592 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -45,7 +45,8 @@ from snakemake.exceptions import WorkflowError # --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 @@ -76,7 +77,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( diff --git a/workflow/bin/sp b/workflow/bin/sp index bcd5295d6..e10628214 100755 --- a/workflow/bin/sp +++ b/workflow/bin/sp @@ -38,7 +38,7 @@ # # 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 (scripts/run_config.py). +# SP_PROFILE and input_type, then the input_types: entry (scripts/run_config.py). # To check a resolved value without running: # # python scripts/run_config.py workflow/config.yaml ~/my_run.yaml outputs.run_dir diff --git a/workflow/config.yaml b/workflow/config.yaml index 266413d30..176ae1d54 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -16,19 +16,26 @@ # 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. run: smk-g6 -# tile_detection is set per (machine, input_type) in the table below: -# `unions_catalogue` fetches the UNIONS per-tile catalogue (inputs.catalogues: -# a local directory, symlinked in, or a vos: URL) and converts it in place; -# `sextractor` (the default when no entry sets it, as for image sims) runs -# SExtractor on the tile image. +# Defaults per input_type: what follows from the kind of input, whatever the +# machine. A run config, or a machines: entry below, overrides them; a +# top-level key here would shadow both tables, so none of these is set above. +# - tile_detection: `unions_catalogue` fetches the UNIONS per-tile catalogue +# (inputs.catalogues) and converts it in place; `sextractor` runs SExtractor +# on the tile image +# - psf_model: psfex, mccd, or fake (image sims only: the true simulation +# PSF, read from psf_dict) +input_types: + data: + 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 @@ -36,10 +43,10 @@ run: smk-g6 # - base_dir # - tile_list: text file with tile IDs # - retrieve: method to get input data, allowed are symlink, vos -# - tile_detection: unions_catalogue or sextractor (default) # - inputs.tiles,, .exposures: path for input tile and exposure -# - inputs.catalogues: UNIONS per-tile catalogues (tile_detection: -# unions_catalogue) +# - 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 # - outputs.index_db: path to bookkeeping index sqlite file @@ -49,9 +56,6 @@ machines: nibi: base_dir: /project/def-mjhudson data: - # @sc [decision:star_selection_psf.psf_modelling_software] - psf_model: psfex - tile_detection: unions_catalogue tile_list: $base_dir/cdaley/sp-products/$run/tiles.txt retrieve: symlink inputs: @@ -69,25 +73,19 @@ machines: candide: base_dir: /n17data/UNIONS/WL data: - # @sc [decision:star_selection_psf.psf_modelling_software] - psf_model: psfex - tile_detection: unions_catalogue retrieve: vos inputs: catalogues: TBD 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/workflow/scripts/run_config.py b/workflow/scripts/run_config.py index 59c99d37e..e0db7c98c 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,16 +21,8 @@ import yaml PLACEHOLDER = "TBD" -# Keys the machines: table may default, per (machine, input_type). -# psf_model/psf_dict/tile_detection belong here because they are -# per-input_type facts (image sims have no UNIONS catalogue), -# 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", "tile_detection") +# The tables themselves: read here, never defaults or $-expanded. +TABLES = ("input_types", "machines") REQUIRED = ("tile_list", "inputs.tiles", "inputs.exposures", "outputs.run_dir", "outputs.index_db") @@ -65,19 +64,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 `$run` to the top-level `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` 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 @@ -95,17 +100,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 @@ -130,7 +131,7 @@ def load(config_yaml, run_config=None): 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__": From 9d1e023adbb343f9818de439bafca754f1e69a38 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 28 Sep 2026 17:30:36 +0200 Subject: [PATCH 12/21] Record tile_detection in astra.yaml Real data now takes its tile sample from Gwyn's UNIONS catalogue. The detection sub-analysis gains a tile_detection decision (committed: unions_catalogue), pinned to the input_types table, and notes that its SExtractor decisions govern only the sextractor path. The converter's multi-epoch post-processing is tagged under epoch_membership_ccd_bounds, which applies to both. Co-Authored-By: Claude Opus 5.5 --- astra.yaml | 37 +++++++++++++++++++++---- universes/committed.yaml | 1 + workflow/config.yaml | 1 + workflow/config/cfis/config_tile_Uc.ini | 2 ++ 4 files changed, 35 insertions(+), 6 deletions(-) diff --git a/astra.yaml b/astra.yaml index e6902a189..16401b79a 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,33 @@ 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, 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: >- @@ -600,9 +622,12 @@ 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: diff --git a/universes/committed.yaml b/universes/committed.yaml index 63b6e795b..05d88e097 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 diff --git a/workflow/config.yaml b/workflow/config.yaml index 176ae1d54..edcd19fd7 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -29,6 +29,7 @@ run: smk-g6 # 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 diff --git a/workflow/config/cfis/config_tile_Uc.ini b/workflow/config/cfis/config_tile_Uc.ini index a13def13f..512eeb39b 100644 --- a/workflow/config/cfis/config_tile_Uc.ini +++ b/workflow/config/cfis/config_tile_Uc.ini @@ -80,10 +80,12 @@ VIGNET_SIZE = 51 ## 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 From 21e7c4e7b82b08be2cd4861317f5c9257824bc37 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 28 Sep 2026 17:43:47 +0200 Subject: [PATCH 13/21] Point inputs.catalogues at vos:cfis/tiles_DR6 on both machines Stephen Gwyn's per-tile catalogues live at vos:cfis/tiles_DR6 (the path develop's config_tile_Git_cat_vos.ini already uses). tile_get_catalogue fetches them with vcp inside the job, using ~/.ssl/cadcproxy.pem; that works from a candide compute node in the runtime container. The unset/placeholder check moves from the Snakefile into run_config.catalogue_source, so the parse-time refusal is unit-tested. Co-Authored-By: Claude Opus 5.5 --- tests/unit/test_workflow_tile_detection.py | 24 +++++++++++++++++++++- workflow/Snakefile | 14 ++++--------- workflow/config.yaml | 8 ++++---- workflow/scripts/run_config.py | 17 +++++++++++++++ 4 files changed, 48 insertions(+), 15 deletions(-) diff --git a/tests/unit/test_workflow_tile_detection.py b/tests/unit/test_workflow_tile_detection.py index 95a174023..8f5d7275e 100644 --- a/tests/unit/test_workflow_tile_detection.py +++ b/tests/unit/test_workflow_tile_detection.py @@ -139,6 +139,28 @@ def test_defaults_pair_the_catalogue_with_its_source( str(REPO_ROOT / "workflow" / "config.yaml"), str(over)) if input_type == "data": assert config["tile_detection"] == "unions_catalogue" - assert config["inputs"]["catalogues"] + 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/workflow/Snakefile b/workflow/Snakefile index 76845e592..a95035943 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -164,16 +164,10 @@ if TILE_DETECTION not in TILE_DETECTIONS: raise WorkflowError( f"Invalid tile_detection={TILE_DETECTION!r}; expected one of " f"{sorted(TILE_DETECTIONS)}.") -# The catalogue source: a local directory (symlinked in) or a vos: URL -# (downloaded); the retrieve mode follows from the prefix. -CATALOGUES = INPUTS.get("catalogues") or "" -if CATALOGUES == run_config.PLACEHOLDER: - CATALOGUES = "" -if TILE_DETECTION == "unions_catalogue" and not CATALOGUES: - raise WorkflowError( - "tile_detection=unions_catalogue needs `inputs.catalogues:` (a local " - "directory or vos: URL holding the per-tile CFIS..r.cat files).") -CATALOGUE_RETRIEVE = "vos" if CATALOGUES.startswith("vos:") else "symlink" +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 diff --git a/workflow/config.yaml b/workflow/config.yaml index edcd19fd7..b0b75d903 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -62,9 +62,9 @@ machines: inputs: tiles: $base_dir/unions-wl/tiles exposures: $base_dir/unions-wl/exposures - # UNIONS per-tile catalogues (CFIS..r.cat): a local directory - # or a vos: URL (e.g. vos:cfis/tiles_DR6). - catalogues: TBD + # UNIONS per-tile catalogues (CFIS..r.cat), 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 @@ -76,7 +76,7 @@ machines: data: retrieve: vos inputs: - catalogues: TBD + catalogues: vos:cfis/tiles_DR6 container: /n17data/cdaley/containers/shapepipe_develop-runtime-20260718.sif image_sims: retrieve: symlink diff --git a/workflow/scripts/run_config.py b/workflow/scripts/run_config.py index e0db7c98c..f6799aa5c 100644 --- a/workflow/scripts/run_config.py +++ b/workflow/scripts/run_config.py @@ -125,6 +125,23 @@ def unresolved(config): if get(config, k) in (None, "", PLACEHOLDER) or "$" in str(get(config, k))] +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 {} From 9f034f823afe90707dc556d8138f35e01836fe0b Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 28 Sep 2026 20:59:23 +0200 Subject: [PATCH 14/21] Leave the retired standalone merge_final_cat.py to develop The workflow merges final catalogues through final_cat_merge; the standalone script is retired, so this branch no longer patches it or tests filter_available_columns. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01DksuyF9YrZAwHQeAMXsxp4 --- scripts/python/merge_final_cat.py | 44 --------------------------- tests/unit/test_merge_final_cat.py | 48 ------------------------------ 2 files changed, 92 deletions(-) delete mode 100644 tests/unit/test_merge_final_cat.py diff --git a/scripts/python/merge_final_cat.py b/scripts/python/merge_final_cat.py index 04fd3bf27..9d9c10d21 100755 --- a/scripts/python/merge_final_cat.py +++ b/scripts/python/merge_final_cat.py @@ -247,42 +247,6 @@ def read_param_file(path, verbose=False): return param_list -def filter_available_columns(param_list, available_columns): - """Filter Available Columns. - - Return the subset of ``param_list`` present in ``available_columns``, - printing one line for each requested column that is missing (e.g. - ``TILE_UNIQUE_ID``, absent from final catalogues made before make_cat - wrote it) so the merge still proceeds for such catalogues. - - Parameters - ---------- - param_list: list of str - requested column names - available_columns: iterable of str - column names present in the catalogue - - Returns - ------- - list of str - subset of ``param_list`` present in ``available_columns`` - - """ - if not param_list: - return param_list - - available_columns = set(available_columns) - missing = [p for p in param_list if p not in available_columns] - - for p in missing: - print( - f"Column '{p}' not found in input catalogue, skipping in " - "merged catalogue" - ) - - return [p for p in param_list if p in available_columns] - - def get_data(path, hdu_num, param_list): """Get Data. @@ -381,14 +345,6 @@ def main(argv=None): if param.verbose: print(f"{len(lpath)} files files to merge found") - # Drop requested columns absent from the catalogues (e.g. TILE_UNIQUE_ID - # in older final catalogues) so the merge still proceeds. - with fits.open(lpath[0]) as hdu_list: - available_columns = hdu_list[param.hdu_num].columns.names - param.param_list = filter_available_columns( - param.param_list, available_columns - ) - count = 0 # Determine number of columns and keys from first catalogue file diff --git a/tests/unit/test_merge_final_cat.py b/tests/unit/test_merge_final_cat.py deleted file mode 100644 index d7de17ace..000000000 --- a/tests/unit/test_merge_final_cat.py +++ /dev/null @@ -1,48 +0,0 @@ -"""The standalone ``scripts/python/merge_final_cat.py`` drops catalogue-param -columns absent from its input catalogues. - -``TILE_UNIQUE_ID`` is absent from final catalogues made before make_cat wrote -it; ``filter_available_columns`` makes that column (and any other -requested-but-absent one) optional for this manual tool, logging one line per -drop. The workflow's campaign merge (``final_cat_merge``, through -``create_final_cat.read_data``) does the opposite and raises on a missing -column; tests/unit/test_final_cat_merge_invariants.py holds it to that. -""" - -import importlib.util -from pathlib import Path - -REPO_ROOT = Path(__file__).resolve().parents[2] -MERGE_SCRIPT = REPO_ROOT / "scripts" / "python" / "merge_final_cat.py" - - -def _load(script): - """Import a script by path — ``scripts/python`` is not a package.""" - assert script.exists(), f"{script} not found; the rule calls it by path" - spec = importlib.util.spec_from_file_location(f"_{script.stem}", script) - module = importlib.util.module_from_spec(spec) - spec.loader.exec_module(module) - return module - - -def test_present_columns_are_kept(): - m = _load(MERGE_SCRIPT) - param_list = ["XWIN_WORLD", "TILE_ID"] - assert m.filter_available_columns( - param_list, ["XWIN_WORLD", "TILE_ID", "FLAGS"] - ) == param_list - - -def test_missing_column_is_dropped_and_logged(capsys): - m = _load(MERGE_SCRIPT) - kept = m.filter_available_columns( - ["XWIN_WORLD", "TILE_UNIQUE_ID"], ["XWIN_WORLD", "FLAGS"] - ) - assert kept == ["XWIN_WORLD"] - out = capsys.readouterr().out - assert "TILE_UNIQUE_ID" in out - - -def test_empty_param_list_means_copy_all_columns(): - m = _load(MERGE_SCRIPT) - assert m.filter_available_columns([], ["XWIN_WORLD"]) == [] From 2905a8829c55e01469ca93f10e303d0cb3ccd9a2 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 29 Sep 2026 13:29:59 +0200 Subject: [PATCH 15/21] final_cat_merge: request only the columns the tile's detection mode can carry Under tile_detection: unions_catalogue the detection columns are the UNIONS DR6 catalogue's, which lacks nine SExtractor columns final_cat.param names (MAG_WIN, MAGERR_WIN, FLUX_AUTO, FLUXERR_AUTO, FLUX_APER, FLUXERR_APER, SNR_WIN, FWHM_IMAGE, FWHM_WORLD), so every catalogue-mode merge raised. merge_final_cat.py now takes --tile-detection and subtracts exactly SEXTRACTOR_ONLY_COLUMNS in that mode; the merge stays strict on the rest. tests/module/test_final_cat_columns.py runs the converter on the DR6 header and asserts that list is exactly the requested SExtractor columns DR6 lacks. The params pin moves for final_cat_merge only. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BeQNTBcjTqu6TPktxysovQ --- tests/module/test_final_cat_columns.py | 141 ++++++++++++++++++ tests/unit/test_final_cat_merge_invariants.py | 2 +- tests/workflow/params_pin.json | 4 +- workflow/config/cfis/final_cat.param | 4 + workflow/rules/tile.smk | 2 + workflow/scripts/merge_final_cat.py | 48 +++++- 6 files changed, 197 insertions(+), 4 deletions(-) create mode 100644 tests/module/test_final_cat_columns.py 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/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/workflow/params_pin.json b/tests/workflow/params_pin.json index 2298fad2c..eeb61d950 100644 --- a/tests/workflow/params_pin.json +++ b/tests/workflow/params_pin.json @@ -8,7 +8,7 @@ "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": "1b4b0871bf28296f8e5d6855ae47ff2cb84542c93593a5ae5c0ffc3c84c45479", @@ -24,7 +24,7 @@ "tile_vignets": "9d4ae0d99c18217f2185f245281312454c8a219ec1628176e08c71a5efc4dc91" }, "schema": 1, - "sha256": "1fc6e355ed162e97d5252038b14ed2ae9e52df19ce83485ec42d395232657f3d", + "sha256": "f9772987e517e9827c8d501d310ec323be251a62b7b89aad0725a03f8f4f120e", "unit_pre": { "exp_get_images": "8dec850af212879f225fcf27a5f1281e1a075264c7b97d38c2214395d360168c", "exp_psf": "f2358ddf7385918dc5033d10b37f6dc97a15d02b071a3ea0a4619a5f7e6f5bec", diff --git a/workflow/config/cfis/final_cat.param b/workflow/config/cfis/final_cat.param index 9a1850b46..cc5feaabf 100644 --- a/workflow/config/cfis/final_cat.param +++ b/workflow/config/cfis/final_cat.param @@ -226,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/rules/tile.smk b/workflow/rules/tile.smk index 15c6918f1..289d636bd 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -1040,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), @@ -1064,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/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 From 8da837be818c68f85bd0f50ac84e8e6bd7c35670 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 29 Sep 2026 16:16:06 +0200 Subject: [PATCH 16/21] Mask neighbours in catalogue-mode VIGNETs from the DR6 segmentation map tile_detection: unions_catalogue cut VIGNET straight from the tile image, so ngmix (which masks where the tile VIGNET is -1e30) masked no neighbours: in smk-g10 NGMIX_MAG_NOSHEAR ran >1 mag bright of MAG_AUTO for 6.1% of objects, against 2.7% with SExtractor. tile_get_catalogue now also fetches CFIS..r.seg.fits.fz. The converter relabels it to the catalogue's NUMBER by the footprint under each object's centre pixel (relabel_seg: unclaimed footprints -1, a small disc for objects without a footprint), writes it beside the sexcat as seg-.fits for UberSeg, and sets VIGNET pixels on other objects' footprints to -1e30. Off-image VIGNET pixels are -1e30 too, as SExtractor writes them. Records detection.catalogue_neighbour_marking. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BeQNTBcjTqu6TPktxysovQ --- astra.yaml | 34 ++- .../read_ext_sexcat.py | 195 +++++++++++-- .../modules/read_ext_sexcat_runner.py | 23 +- tests/module/data/dr6_202.301_seg_patch.fits | Bin 0 -> 34560 bytes tests/module/test_read_ext_sexcat.py | 268 +++++++++++++++++- tests/unit/test_workflow_tile_detection.py | 14 +- universes/committed.yaml | 1 + workflow/README.md | 8 +- workflow/config.yaml | 3 +- workflow/config/cfis/config_tile_Gic.ini | 11 +- workflow/config/cfis/config_tile_Uc.ini | 13 +- workflow/rules/tile.smk | 23 +- workflow/scripts/completeness.py | 9 +- 13 files changed, 539 insertions(+), 63 deletions(-) create mode 100644 tests/module/data/dr6_202.301_seg_patch.fits diff --git a/astra.yaml b/astra.yaml index cf9d37faa..c3810ac2b 100644 --- a/astra.yaml +++ b/astra.yaml @@ -348,7 +348,8 @@ analyses: [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: @@ -633,6 +634,37 @@ analyses: 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 + written beside the sexcat as seg-.fits for an UberSeg seg run; + 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: >- 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 3bbccdcbd..43f5fb61e 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,122 @@ 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; UberSeg 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 _image_header(header): + """Header of an image HDU without its compression or extension cards.""" + header = header.copy() + for key in ("XTENSION", "PCOUNT", "GCOUNT", "EXTNAME"): + header.remove(key, ignore_missing=True, remove_all=True) + return header + + +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,32 +174,36 @@ 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 @@ -161,6 +275,8 @@ def make_ldac_from_ascii( image_path, output_cat_path, stamp_size=51, + seg_path=None, + seg_output_path=None, w_log=None, ): """Convert an external ASCII catalogue to FITS-LDAC format. @@ -174,7 +290,11 @@ def make_ldac_from_ascii( The input columns, ``NUMBER`` included, are copied unchanged; a ``VIGNET`` column (postage stamps extracted from the tile image) is - added to ``LDAC_OBJECTS``. + 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`), written to ``seg_output_path``, and + used to set neighbours' pixels in each ``VIGNET`` to ``BIG``, as + SExtractor does. Parameters ---------- @@ -186,6 +306,11 @@ def make_ldac_from_ascii( Path to the output FITS-LDAC catalogue 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) + seg_output_path : str, optional + Path to write the relabelled segmentation map to, required with + ``seg_path`` w_log : logging.Logger, optional Pipeline logger @@ -199,6 +324,32 @@ def make_ldac_from_ascii( 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_header = hdu.header + 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 + fits.PrimaryHDU(seg, header=_image_header(seg_header)).writeto( + seg_output_path, overwrite=True + ) + if w_log: + w_log.info( + f"Relabelled {seg_path} to NUMBER, written to" + + f" {seg_output_path}: " + + ", ".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}" @@ -208,6 +359,8 @@ 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) diff --git a/src/shapepipe/modules/read_ext_sexcat_runner.py b/src/shapepipe/modules/read_ext_sexcat_runner.py index 5c6e7b1db..daa6a9351 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 and the map, relabelled to the catalogue's + NUMBER, is written as ``seg.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") @@ -56,19 +63,21 @@ def read_ext_sexcat_runner( image_path, output_path, stamp_size=stamp_size, + seg_path=seg_path, + seg_output_path=f"{run_dirs['output']}/seg{file_number_string}.fits", 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/tests/module/data/dr6_202.301_seg_patch.fits b/tests/module/data/dr6_202.301_seg_patch.fits new file mode 100644 index 0000000000000000000000000000000000000000..c465e6446fa67de59ea3da93eb2d059b3befcf17 GIT binary patch literal 34560 zcmeHQ3v^Z0ncnx2n<%eBfp6RP&ZgEL zt!uM^xxKS>TW@>Uwq6T~tyW%NSH8L|2EIDmz$B{@S(0PZ!_pc^Yap$Gv|Q_54}p<3%y@XUgxdZ*S?_Lg~1YxP6<}w6m|(KN*MU5?`?|zT!e>e&5q|eXj<5 z3b9jdW~*lTnk!c$kN2DxiWbhF=ajdtr)zt+*4^0C*xB0G+T+jSidguFVv$3C5}!X7 zlRd4jqO!cItk%n)7vJXg){Yh_acg7q77gNCXv|t&-MX4`X^&aU*HqW6Ew8ODFWXqQ ztZtdz4b9;RD|HoX)>bQgMKz_fOA1To3=a`gYN)JUwj$p7(pnNYY4$a?98Ilu2XxPp`26FQiU+CUTV7p5OBnmd+REkS8?p49{%kA;+-PWx9c^7b z?R{H1-CS91*teSFL21>R<*RBd8|)hr-RG#qMfjRJx|+9Yl4A>lv(|9cvMTW8=Pzn~ zbhqzl?da9IaR;N&Ip@1ul0T`xS*7dh%IU7iDX-cs2%W^N%`VE*iu1I&d0JPGwqTA2 zANT`0U;ebYSidRa!~GCx5cp=*mamZ1ocwLWZXx|;OM7n;n>)E5VFAA4!s7gr!Xk&f zl0YgS^}Oc2>tlXlVM!{O)Yx?~@s%u?cWLVQunAR{V*{>f?8BagG~l|f-bKxu+j|T8 zXn}1kTRX2{c|~PyK5HYtC%?C~Eq`-+UvK{A8}hWmOPU%f@WapBu-Dfc_zLq3PE+T} zJ|lf!XV=DOUl&50g@05C8viD6J%`~;hL7Nf141&9)?4zKeU?Rtqh-1FdSsKgy9Ws*FSfYoR}1PDA?e^T%);`e5#@jXD`JG%=I zeP{o>f}^ff{iB{%a5U)(&qx0wAmknMIY9J{`3B%-z-JX4TgdQEzz!kCo=|Wc@r!Ze z{!PK`MGX5D93NqbIpXhx82xc0JwSKn#|mCRt|iaS7=I;qwZXTnY&MCHE*eH&gY`?Pd51AmqK2;t9c*6`XgidOrVbK+K;{e)Iz1 zYkx(C-_GRN4{ss_uFaNHBOG*KE2(jb{V7U-W$qxwL1-J^3 z#^bp@8xV5q{R%FNFa$5uzYM$(EL{Ky{H5Pka5>f$wU-YnSVsClyFAE{^niB7T!z0; za0PHt{|d-J@~$L)&{s)zg0@zAZ&ldL53eHxGJ9^&Byaq z4>EjH!K&#D8yHfY!MG}t5BRIcGMvm%V>k!!RF=s?3yWxWQ0Q3!Pn8It7(jG>+djfz zC*D+y$@~wRgmV6!Ch}7}l@^&eK`tRq?6*XsHp=>#a`I3bbaM?FfPu&+=miy3lAoJJ z6SU`2n2Yh@sY=10>RgSeNqJM}*_sO+rmlYEbs?o%3<}pr@fo~$!0Bh`OTjW2n^hv+ z6tx)C2x0DF*qSQkI~*)g7nOX|Bb}a)M_G>;_?^7Ud9mlxlk#1P7^-K`T;e?UanIVI zZYt0+O884rdBb{fYUbL5_Ze*im zt7H2nRJRmbd*HAU!Ig!`e#$3c`f?OfI)OE6Jk}86XH%rJ8!q40LQ7qQ77c-h^bH}N z-oonvU5reFf?M4Rb|Mpz>qtBDXEtV{jCwe+OOn_6S|}c6Y%HUTc-bOWi~!}R5Y3PB z>QA0BW+?^g$os!9q%M^PMDPg2n)H%eXW)EP2c})U-UzW_TRXVu3@;!ZsO2d?MWD0Wa7dzOqlkGpu>Or>2c+k*DM~d;C$>}6 zaojr~*XBr4lzJ6o4l3P3uTpJd%wq$jY=>S&&Q?;qzg{`I6x*T;g!Et7L^A-DI%x)l zU>ho1lA%ML-DJ-Zpo!&gsTH<#{3w@&e)qe*8Zw}m9!tq3b=*>8%jsjpK{{%Qv9CD7 zrV)zi(RD~)G^0aRshdj106szoy=NwLFZMqx;E62zlG`b1Y)(r+TqL6I%~jzMo!wz? zcEZoZ$y3P!p9)*xu>n^iG_q8rd_JUxNu}6Yg7yGFb;QIMDFVW~oIy}eg%;LUh*WXrymvO=(F@{#e&zNrGKZQvAD7%LB$AI7RuXCFaFlG& z_^bOd9+luh?`CbEkuARWfln0;eEy$4BeEBOcOl|)fHP+lewPvr27oY4_+N|*u?W^*1wF&+8d}YlP8^d6TkaSoU?c zEc#D3bp^WB1*DZ-fXU3Qj^*Pz*7>Y-hZ3bmEEPogr1jDco-BZ|(n^ASjJVl)1`|0b|IfxiE z&J|S8qwgKXqIW-^*+Zf_syELMlwc-62B;~vs3$thR37M zqcw#xQlb$ql$;jv3||l6V!(bG;46{YISkQ^#R@G_A*&IrKWZK!gHOdzk7&3omp)B7C8d?&Ym#Txt&NWzl4#+fPn~k)F#U=b{z~o200JAj!Ejws58;V$yS= zW*LVZEMoSkC#?9i!D)WEv5yVYo?`rn%)o$gO{<#~{IZGsWjS6`9eYfXJc_v9pQp zX`ekN?@2}Go__g2jGH=OY`cC^CPh;hL;U@BGj6c z$Bqm|k$EvK(j-?UXeO|uby%{+b4&}|sTr9~>kemlAwj4lXG09MR-~NgJ69!9j`$&? zO>%SPk~N~8Pb7#(*EJVIM;)?<)A)1AX;3m=kK!}XW8vhLW}GCQ z%S~Oggl+-Gv7XhM%0URM(Up`c$+?Xog}I5#pkPK#Jf2oe!AJytIz-{Zd=gEHF-@Jt zuzUvhF{zB_levi^l2nGFF?B7jR&C}K6?U$5ML4G@mQx!tf^y&~QUxlZ-L5t7Q4t35 zsb(5U*>ji2zBawts%}oD>Bbh=W0RAX5(d8|e24vtee{D^!dZUTuV$Z-FGs;%=r&zi zWx?YPd)oMtRg6SoO2{r`3C&4Fi3bjMo}tI@DW`jOlgm<;|JnGpzJ><%zzc@m{vod2 zjBkP{io3N8DuijUUr|SW?9DWY3-60)g&vU8ssHFPU*4Ea+=5URBSh0&!d*1U=PM2o zI9g+*TtfU&{>lzvMj19m5t*T@Oy7@aLW&@LUam?_iH&4rOqL@4wquz}>oU5K_=^di1HhwDK9OPcCGSw-2IK36(aJ?<1l}9C!c|VjULW7&s*G)3F zB^KRGQqQ8x;x40r?Y7$A(^M}oSxjCuq(V<{274zxo#Z(pIMxR7f|9}ebnIP8jT7NL z6o zpwnUGE{-QwbWO=&^`J2Y6+xmNR*zzK;i`=;My~ z9SA5*OK@h8CS6Lu&`j8FoFWevEzTe%ul95jK{{scm|I;@vz9UK3pusQJPO9^7j23qmqlal zIoxn-$#&Nh5 z%J~Ol9cdEf`ZW8-*^2T9iky2Nq#-CXx5UW6UHgpTt=86$!;vcC-8dzj_nX;* z+z$9d{+gpF-Z?ThR)&h^s^M5-1HY zCtzQp;<2JOY`)NQ21>d19T5SD=-vk)FdqMVg6*yf4JvMBs@_;}!tZy}Z3in?#J60= zO+{B^hPzUDwt3p-_8CcqgQm7?~AcZnUPW94e~jdsl@dnjN!{L8l9V` zwq=i!X~R>+&uiibJ)UT!ZkEHuV+@nB_Ne$m^G^k1iZ4T1QirO`xJdU|cQl}Ei%i*7 zwtJc|2DZ#=;3OE+Ne&cmJ|5SUYEnuX!<19bm|{zqks=9G(iC}tjaJ>5Un&Oi3tRmu0Cl1`L7mM#AJ9E|8r`Lf&f`U&eT>c?DzJ zarlvcxS)qNjiHsb+ssR1F~)hkLTp7e4`XV8%d@1i5^GH6@P>~OqOy8+Hkda7V>%Ny zylt3_C8FZmZusj*<9d>B*Kq|QZ%fVZbzBu@H6bmhH?TxJ*|~Hr1zRnIZ}hC>(NRmTuM1P+4d$RjO8+K8;rRXYL6$4byO0Fc&4*WX3q? zoO%o6AQDeC{{ z_rpj3VSbDf^Zh9_-uXU<)9?K$CD5PsK1lVR7=DXN3~@Ev`%_k5SzTJb*7^R0Skq$m z+uomY1rM^nKOpA3SPvwBuovDRP`^>V@#N(2(feH#zOoVHuSD{#xA}N+$Vchgfzn`qS?ZFf@^dX$_<`kk&w21D~Y^@V_-3&^MjPD!TC9 zs9x1rpWReFfd2<4s9$;S%1Or}QGL!e-_FUp8%Kb75q;5v_cy-vKSD2hW6W2o=`wEK zIK5)@wWAN*DfIF+KYVxGM^W@m(p$Fn{^vin_3O>i?f?C8E{&U}Z|?r;iGQBaug^cH z@f>|MjjPh@-+1{458bEhmjyrg?(ui1Z-c(%x7yVM_&bhe!Ts9ddj`hY{I8C6MbY5>;DCO*S+z;(zk(ceviK9C*OE7d{b0kzjfB3 zg7;|L*Yx7!muEh+O6VJwV%X0$^!>G7He>LMU%_uuZFv2yuitv@5c=-XBVTyu*P9{7 zHQ$`EaNu{;w@+X9yXf4`CkOP6_pQ9~U%#oN?||Mm72}@Q^^IM3J^kjFsqcM#@~>}- z%zRDgO&kC4tpk_f4t)M`eeK6XV{)#I>dohG%YWo>A^tcypkK4H296-~<|&UZtk2$r zzKC8SWPyFZ)<0`)4f@9EwY#L=^yWk3?;YBH41GEJhMS?%2BEjUctRo(3mWuwpJd$~ l+^*|wA3k!p_*u+Nm$U}b8c1s(t%0-#(i%u>;LOs%{{cH+$~*u7 literal 0 HcmV?d00001 diff --git a/tests/module/test_read_ext_sexcat.py b/tests/module/test_read_ext_sexcat.py index a4067179c..b9afa421d 100644 --- a/tests/module/test_read_ext_sexcat.py +++ b/tests/module/test_read_ext_sexcat.py @@ -5,10 +5,15 @@ 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. The last test follows the catalogue through -``make_cat.save_sextractor_data``, which builds ``TILE_UNIQUE_ID``. +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 @@ -85,7 +90,7 @@ def test_number_is_kept_and_aliases_added(ldac): npt.assert_array_equal(data["YWIN_WORLD"], data["DELTA_J2000"]) -def test_vignets_are_cut_from_the_image_and_zero_padded(ldac): +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) @@ -94,14 +99,15 @@ def test_vignets_are_cut_from_the_image_and_zero_padded(ldac): 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 zero. - assert (vignets[1, :, :2] == 0).all() + # 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:, :] == 0).all() - assert (vignets[2, :, STAMP // 2 + 1:] == 0).all() + 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 @@ -117,3 +123,251 @@ def test_tile_unique_id_reaches_the_final_catalogue(ldac, tmp_path): 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_writes_the_relabelled_map_and_marks(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" + seg_out = tmp_path / "seg-301-279.fits" + 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) + header = fits.Header() + header["CRPIX1"] = 10.0 + fits.HDUList([fits.PrimaryHDU(), + fits.CompImageHDU(_seg_map(), header=header)]).writeto(seg_in) + + rs.make_ldac_from_ascii(str(cat), str(img), str(out), stamp_size=SEG_STAMP, + seg_path=str(seg_in), seg_output_path=str(seg_out)) + + with fits.open(seg_out) as hdul: + relabelled = hdul[0].data + assert hdul[0].header["CRPIX1"] == 10.0 + seg = _seg_map() + 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), + seg_output_path=str(seg_out)) + + +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 diff --git a/tests/unit/test_workflow_tile_detection.py b/tests/unit/test_workflow_tile_detection.py index 8f5d7275e..d6c84f070 100644 --- a/tests/unit/test_workflow_tile_detection.py +++ b/tests/unit/test_workflow_tile_detection.py @@ -66,6 +66,18 @@ def test_fetch_ini_matches_its_stage_and_feeds_the_converter(): 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" @@ -81,7 +93,7 @@ def _stage_dir(tmp_path, runner, n): @pytest.mark.parametrize("mode, runner, expect", [ ("sextractor", "sextractor_runner", 2), - ("unions_catalogue", "read_ext_sexcat_runner", 1), + ("unions_catalogue", "read_ext_sexcat_runner", 2), ]) def test_tile_detect_is_checked_per_mode(tmp_path, monkeypatch, mode, runner, expect): monkeypatch.setenv("SP_TILE_DETECTION", mode) diff --git a/universes/committed.yaml b/universes/committed.yaml index 05d88e097..bc34d7601 100644 --- a/universes/committed.yaml +++ b/universes/committed.yaml @@ -24,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 24c0dbe4c..6e91b344f 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -29,9 +29,11 @@ uv pip install 'snakemake>=9,<10' 'snakemake-executor-plugin-slurm>=2.7,<3' # 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) 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. +# 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. diff --git a/workflow/config.yaml b/workflow/config.yaml index c37df88f5..86d8c5a9a 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -64,7 +64,8 @@ machines: inputs: tiles: $base_dir/unions-wl/tiles exposures: $base_dir/unions-wl/exposures - # UNIONS per-tile catalogues (CFIS..r.cat), fetched by vcp in + # 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: diff --git a/workflow/config/cfis/config_tile_Gic.ini b/workflow/config/cfis/config_tile_Gic.ini index ab38f38fe..1e6bdb7f5 100644 --- a/workflow/config/cfis/config_tile_Gic.ini +++ b/workflow/config/cfis/config_tile_Gic.ini @@ -1,4 +1,5 @@ # 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) @@ -53,7 +54,7 @@ TIMEOUT = 96:00:00 ## Module options -# Get the external tile catalogue +# Get the external tile catalogue and its segmentation map [GET_IMAGES_RUNNER] FILE_PATTERN = tile_numbers @@ -67,19 +68,19 @@ NUMBERING_SCHEME = # 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 +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'; the workflow derives it from # inputs.catalogues (vos: URL -> vos, local directory -> symlink) diff --git a/workflow/config/cfis/config_tile_Uc.ini b/workflow/config/cfis/config_tile_Uc.ini index 512eeb39b..e349ee65a 100644 --- a/workflow/config/cfis/config_tile_Uc.ini +++ b/workflow/config/cfis/config_tile_Uc.ini @@ -61,11 +61,11 @@ TIMEOUT = 96:00:00 [READ_EXT_SEXCAT_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_Mh_exp/merge_headers_runner/output +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 @@ -77,6 +77,13 @@ 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; and write +# the map relabelled to the catalogue's NUMBER as seg-.fits beside the +# sexcat (NUMBER on each object's footprint, -1 on unclaimed footprints, 0 sky) +# @sc [decision:detection.catalogue_neighbour_marking] +SEGMENTATION = True + ## Post-processing # Necessary for tiles, to enable multi-exposure processing diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index 289d636bd..6a36e9fee 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -42,17 +42,20 @@ 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 -(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, and the multi-epoch post-processing. It keeps the catalogue's own +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. What -this path does not have is a segmentation map: ngmix's ``BLEND_HANDLING = -uberseg`` needs the SExtractor mode. +tile_detect.json in both modes, so tile_vignets onwards is the same DAG. Beside +the sexcat it writes seg-.fits, the segmentation map relabelled to the +catalogue's NUMBER (-1 on footprints no object claims, 0 on sky), the map an +UberSeg vignet run cuts its seg stamps from. 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 @@ -498,8 +501,8 @@ if TILE_DETECTION == "sextractor": else: - # Fetch the UNIONS per-tile catalogue (CFIS..r.cat), from a local - # mirror or from vos, the way tile_get_images fetches the image. Reads + # 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: diff --git a/workflow/scripts/completeness.py b/workflow/scripts/completeness.py index 6620ca015..e1ac6fbab 100644 --- a/workflow/scripts/completeness.py +++ b/workflow/scripts/completeness.py @@ -142,12 +142,13 @@ # --- tile post --- "tile_merge_headers": {"merge_headers_runner": dict(expect=1)}, - # The fetched UNIONS catalogue: one .cat per tile. - "tile_get_catalogue": {"get_images_runner": dict(expect=1)}, + # 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)}, - # One FITS-LDAC sexcat, converted from the fetched catalogue. - "unions_catalogue": {"read_ext_sexcat_runner": dict(expect=1)}, + # The FITS-LDAC sexcat converted from the fetched catalogue, and the + # segmentation map relabelled to its NUMBER. + "unions_catalogue": {"read_ext_sexcat_runner": dict(expect=2)}, }, "tile_vignets": { "psfex": { From a4da900ce32d79c50d4543f8cf97e2161280047c Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 29 Sep 2026 16:53:22 +0200 Subject: [PATCH 17/21] completeness: expect the seg-stamp vignetmaker run under blend_handling: uberseg tile_vignets gains a third vignetmaker run, the segmentation stamps ngmix's SEG_VIGNET_PATH reads, when SP_BLEND_HANDLING=uberseg. Noise-fill runs export nothing and expect exactly what they did before. Landing this here keeps SCRIPT_HASH(completeness.py) at one campaign boundary with #897. Co-Authored-By: Claude Opus 5.5 --- workflow/scripts/completeness.py | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/workflow/scripts/completeness.py b/workflow/scripts/completeness.py index e1ac6fbab..3ad8fab1d 100644 --- a/workflow/scripts/completeness.py +++ b/workflow/scripts/completeness.py @@ -78,9 +78,17 @@ # when it is not the default, so a SExtractor run's prologue is unchanged. TILE_DETECTIONS = ("sextractor", "unions_catalogue") +# ngmix's neighbour treatments, mirroring BLEND_HANDLINGS in +# shapepipe.modules.ngmix_package.ngmix. Mirrored rather than imported: the +# Snakefile parses this module in the launcher venv, outside the container +# where shapepipe lives. tests/unit/test_workflow_tile_detection.py asserts +# the two tuples agree. +BLEND_HANDLINGS = ("noisefill", "uberseg") + # stage -> {runner_subdir: {expect, [warn], [subpath]}} # exp_psf and tile_vignets are selected by $SP_PSF at check time, tile_detect -# by $SP_TILE_DETECTION. +# by $SP_TILE_DETECTION; tile_vignets adds SEG_VIGNETS under +# $SP_BLEND_HANDLING=uberseg. # @sc [decision:per_unit_completeness] COMPLETENESS = { # --- tile prepare (phase A) --- @@ -181,6 +189,12 @@ } +# tile_vignets' extra run under blend_handling: uberseg: the one seg_vignet +# file ngmix's SEG_VIGNET_PATH reads. The rules export $SP_BLEND_HANDLING only +# for uberseg, so a noise-fill run's prologue is unchanged. +SEG_VIGNETS = {"vignetmaker_runner_run_3": dict(expect=1)} + + def count_products(run_dir, runner, spec): """Count files in ``run_dir//output[/]/`` (live links only). @@ -224,6 +238,14 @@ def check_counts(stage, run_dir): raise ValueError( f"Invalid SP_PSF={psf_model!r}; expected one of {sorted(table)}." ) from exc + if stage == "tile_vignets": + blend = os.environ.get("SP_BLEND_HANDLING", BLEND_HANDLINGS[0]) + if blend not in BLEND_HANDLINGS: + raise ValueError( + f"Invalid SP_BLEND_HANDLING={blend!r}; expected one of " + f"{', '.join(BLEND_HANDLINGS)}.") + if blend == "uberseg": + table = {**table, **SEG_VIGNETS} elif stage == "tile_detect": detection = os.environ.get("SP_TILE_DETECTION", TILE_DETECTIONS[0]) if detection not in TILE_DETECTIONS: From 97c11a7f464905f8a842c91a48e917cdac77e623 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 30 Sep 2026 00:15:17 +0200 Subject: [PATCH 18/21] Revert "completeness: expect the seg-stamp vignetmaker run under blend_handling: uberseg" UberSeg's seg stamps will ride the tile catalogue as a SEG_VIGNET column, so there is no extra vignetmaker run to expect. completeness.py is again byte-identical to 8da837be, keeping SCRIPT_HASH at the value smk-g11 runs. Co-Authored-By: Claude Opus 5.5 --- workflow/scripts/completeness.py | 24 +----------------------- 1 file changed, 1 insertion(+), 23 deletions(-) diff --git a/workflow/scripts/completeness.py b/workflow/scripts/completeness.py index 3ad8fab1d..e1ac6fbab 100644 --- a/workflow/scripts/completeness.py +++ b/workflow/scripts/completeness.py @@ -78,17 +78,9 @@ # when it is not the default, so a SExtractor run's prologue is unchanged. TILE_DETECTIONS = ("sextractor", "unions_catalogue") -# ngmix's neighbour treatments, mirroring BLEND_HANDLINGS in -# shapepipe.modules.ngmix_package.ngmix. Mirrored rather than imported: the -# Snakefile parses this module in the launcher venv, outside the container -# where shapepipe lives. tests/unit/test_workflow_tile_detection.py asserts -# the two tuples agree. -BLEND_HANDLINGS = ("noisefill", "uberseg") - # stage -> {runner_subdir: {expect, [warn], [subpath]}} # exp_psf and tile_vignets are selected by $SP_PSF at check time, tile_detect -# by $SP_TILE_DETECTION; tile_vignets adds SEG_VIGNETS under -# $SP_BLEND_HANDLING=uberseg. +# by $SP_TILE_DETECTION. # @sc [decision:per_unit_completeness] COMPLETENESS = { # --- tile prepare (phase A) --- @@ -189,12 +181,6 @@ } -# tile_vignets' extra run under blend_handling: uberseg: the one seg_vignet -# file ngmix's SEG_VIGNET_PATH reads. The rules export $SP_BLEND_HANDLING only -# for uberseg, so a noise-fill run's prologue is unchanged. -SEG_VIGNETS = {"vignetmaker_runner_run_3": dict(expect=1)} - - def count_products(run_dir, runner, spec): """Count files in ``run_dir//output[/]/`` (live links only). @@ -238,14 +224,6 @@ def check_counts(stage, run_dir): raise ValueError( f"Invalid SP_PSF={psf_model!r}; expected one of {sorted(table)}." ) from exc - if stage == "tile_vignets": - blend = os.environ.get("SP_BLEND_HANDLING", BLEND_HANDLINGS[0]) - if blend not in BLEND_HANDLINGS: - raise ValueError( - f"Invalid SP_BLEND_HANDLING={blend!r}; expected one of " - f"{', '.join(BLEND_HANDLINGS)}.") - if blend == "uberseg": - table = {**table, **SEG_VIGNETS} elif stage == "tile_detect": detection = os.environ.get("SP_TILE_DETECTION", TILE_DETECTIONS[0]) if detection not in TILE_DETECTIONS: From a8facf5b019866480af66d7b10b178b5b9b66694 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 30 Sep 2026 01:16:27 +0200 Subject: [PATCH 19/21] config.yaml: plainer input_types comment (review suggestion) Co-Authored-By: Claude Opus 5.5 --- workflow/config.yaml | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/workflow/config.yaml b/workflow/config.yaml index 86d8c5a9a..b6f22a01e 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -20,14 +20,12 @@ input_type: data # campaign's merged catalogues. Required, and set by the run config. # run: smk-g6 -# Defaults per input_type: what follows from the kind of input, whatever the -# machine. A run config, or a machines: entry below, overrides them; a -# top-level key here would shadow both tables, so none of these is set above. -# - tile_detection: `unions_catalogue` fetches the UNIONS per-tile catalogue -# (inputs.catalogues) and converts it in place; `sextractor` runs SExtractor -# on the tile image -# - psf_model: psfex, mccd, or fake (image sims only: the true simulation -# PSF, read from psf_dict) +# 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] From bced5a3ce509dffd949093837f0ebddead6036bd Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 30 Sep 2026 01:18:52 +0200 Subject: [PATCH 20/21] Drop the relabelled seg-map file output Its only consumer was an UberSeg seg run that no longer exists; the relabelled map is still used in memory for VIGNET neighbour marking. Co-Authored-By: Claude Opus 5.5 --- astra.yaml | 3 +-- .../read_ext_sexcat.py | 26 +++---------------- .../modules/read_ext_sexcat_runner.py | 1 - tests/module/test_read_ext_sexcat.py | 17 +++++------- workflow/config/cfis/config_tile_Uc.ini | 4 +-- workflow/rules/tile.smk | 5 +--- 6 files changed, 13 insertions(+), 43 deletions(-) diff --git a/astra.yaml b/astra.yaml index c3810ac2b..24dd350d3 100644 --- a/astra.yaml +++ b/astra.yaml @@ -647,8 +647,7 @@ analyses: 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 - written beside the sexcat as seg-.fits for an UberSeg seg run; + 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 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 43f5fb61e..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 @@ -53,8 +53,8 @@ def _build_ldac_imhead(img_header): BIG = -1e30 # The label of a footprint no catalogue object claims in the relabelled -# segmentation map. Negative, so it never collides with a NUMBER; UberSeg only -# asks "self or not self", so neighbours need no identity. +# 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 @@ -145,14 +145,6 @@ def relabel_seg(seg, number, x_image, y_image, fallback_radius=FALLBACK_RADIUS): return out, counts -def _image_header(header): - """Header of an image HDU without its compression or extension cards.""" - header = header.copy() - for key in ("XTENSION", "PCOUNT", "GCOUNT", "EXTNAME"): - header.remove(key, ignore_missing=True, remove_all=True) - return header - - def _extract_vignets(image_data, x_pos, y_pos, stamp_size, seg=None, number=None): """Extract postage stamps from a tile image array. @@ -276,7 +268,6 @@ def make_ldac_from_ascii( output_cat_path, stamp_size=51, seg_path=None, - seg_output_path=None, w_log=None, ): """Convert an external ASCII catalogue to FITS-LDAC format. @@ -292,8 +283,7 @@ def make_ldac_from_ascii( ``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`), written to ``seg_output_path``, and - used to set neighbours' pixels in each ``VIGNET`` to ``BIG``, as + ``NUMBER`` (:func:`relabel_seg`) and used to set neighbours' pixels in each ``VIGNET`` to ``BIG``, as SExtractor does. Parameters @@ -308,9 +298,6 @@ def make_ldac_from_ascii( 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) - seg_output_path : str, optional - Path to write the relabelled segmentation map to, required with - ``seg_path`` w_log : logging.Logger, optional Pipeline logger @@ -328,7 +315,6 @@ def make_ldac_from_ascii( 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_header = hdu.header seg_raw = hdu.data if seg_raw.shape != image_data.shape: raise ValueError( @@ -340,13 +326,9 @@ def make_ldac_from_ascii( cat_data["Y_IMAGE"], ) del seg_raw - fits.PrimaryHDU(seg, header=_image_header(seg_header)).writeto( - seg_output_path, overwrite=True - ) if w_log: w_log.info( - f"Relabelled {seg_path} to NUMBER, written to" - + f" {seg_output_path}: " + f"Relabelled {seg_path} to NUMBER: " + ", ".join(f"{k}={v}" for k, v in counts.items()) ) diff --git a/src/shapepipe/modules/read_ext_sexcat_runner.py b/src/shapepipe/modules/read_ext_sexcat_runner.py index daa6a9351..6a7385378 100644 --- a/src/shapepipe/modules/read_ext_sexcat_runner.py +++ b/src/shapepipe/modules/read_ext_sexcat_runner.py @@ -64,7 +64,6 @@ def read_ext_sexcat_runner( output_path, stamp_size=stamp_size, seg_path=seg_path, - seg_output_path=f"{run_dirs['output']}/seg{file_number_string}.fits", w_log=w_log, ) diff --git a/tests/module/test_read_ext_sexcat.py b/tests/module/test_read_ext_sexcat.py index b9afa421d..d403f5a85 100644 --- a/tests/module/test_read_ext_sexcat.py +++ b/tests/module/test_read_ext_sexcat.py @@ -286,30 +286,26 @@ def test_a_claimed_neighbour_is_masked_by_its_number(): @pytest.mark.decision("detection.catalogue_neighbour_marking") -def test_converter_writes_the_relabelled_map_and_marks(tmp_path): +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" - seg_out = tmp_path / "seg-301-279.fits" 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) - header = fits.Header() - header["CRPIX1"] = 10.0 fits.HDUList([fits.PrimaryHDU(), - fits.CompImageHDU(_seg_map(), header=header)]).writeto(seg_in) + 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_output_path=str(seg_out)) + seg_path=str(seg_in)) - with fits.open(seg_out) as hdul: - relabelled = hdul[0].data - assert hdul[0].header["CRPIX1"] == 10.0 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: @@ -322,8 +318,7 @@ def test_converter_writes_the_relabelled_map_and_marks(tmp_path): 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), - seg_output_path=str(seg_out)) + stamp_size=SEG_STAMP, seg_path=str(seg_in)) DR6_PATCH = Path(__file__).parent / "data" / "dr6_202.301_seg_patch.fits" diff --git a/workflow/config/cfis/config_tile_Uc.ini b/workflow/config/cfis/config_tile_Uc.ini index e349ee65a..ab9d04742 100644 --- a/workflow/config/cfis/config_tile_Uc.ini +++ b/workflow/config/cfis/config_tile_Uc.ini @@ -78,9 +78,7 @@ SUFFIX = sexcat 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; and write -# the map relabelled to the catalogue's NUMBER as seg-.fits beside the -# sexcat (NUMBER on each object's footprint, -1 on unclaimed footprints, 0 sky) +# VIGNET pixels to -1e30, as SExtractor does, so ngmix masks them # @sc [decision:detection.catalogue_neighbour_marking] SEGMENTATION = True diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index 6a36e9fee..f105d23b4 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -52,10 +52,7 @@ 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. Beside -the sexcat it writes seg-.fits, the segmentation map relabelled to the -catalogue's NUMBER (-1 on footprints no object claims, 0 on sky), the map an -UberSeg vignet run cuts its seg stamps from. +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 From 0ce4060cf757c03f31a292247e7664d515f99366 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 30 Sep 2026 01:40:49 +0200 Subject: [PATCH 21/21] completeness: the catalogue converter writes one file The relabelled seg map is no longer written, so expect=2 failed every catalogue-mode tile_detect; a runner-driven test now ties the count to the converter's output. Co-Authored-By: Claude Opus 5.5 --- .../modules/read_ext_sexcat_runner.py | 4 +- tests/module/test_read_ext_sexcat.py | 56 +++++++++++++++++++ tests/unit/test_workflow_tile_detection.py | 2 +- workflow/scripts/completeness.py | 5 +- 4 files changed, 61 insertions(+), 6 deletions(-) diff --git a/src/shapepipe/modules/read_ext_sexcat_runner.py b/src/shapepipe/modules/read_ext_sexcat_runner.py index 6a7385378..17e23163d 100644 --- a/src/shapepipe/modules/read_ext_sexcat_runner.py +++ b/src/shapepipe/modules/read_ext_sexcat_runner.py @@ -33,8 +33,8 @@ def read_ext_sexcat_runner( 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 and the map, relabelled to the catalogue's - NUMBER, is written as ``seg.fits``. MAKE_POST_PROCESS runs the + 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, image_path, *extra_inputs = input_file_list diff --git a/tests/module/test_read_ext_sexcat.py b/tests/module/test_read_ext_sexcat.py index d403f5a85..fa406b9a6 100644 --- a/tests/module/test_read_ext_sexcat.py +++ b/tests/module/test_read_ext_sexcat.py @@ -366,3 +366,59 @@ def test_dr6_marks_every_neighbour_pixel_and_no_own_pixel(): 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_workflow_tile_detection.py b/tests/unit/test_workflow_tile_detection.py index d6c84f070..24775304e 100644 --- a/tests/unit/test_workflow_tile_detection.py +++ b/tests/unit/test_workflow_tile_detection.py @@ -93,7 +93,7 @@ def _stage_dir(tmp_path, runner, n): @pytest.mark.parametrize("mode, runner, expect", [ ("sextractor", "sextractor_runner", 2), - ("unions_catalogue", "read_ext_sexcat_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) diff --git a/workflow/scripts/completeness.py b/workflow/scripts/completeness.py index e1ac6fbab..52720b724 100644 --- a/workflow/scripts/completeness.py +++ b/workflow/scripts/completeness.py @@ -146,9 +146,8 @@ "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, and the - # segmentation map relabelled to its NUMBER. - "unions_catalogue": {"read_ext_sexcat_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": {