From f5c76a407a81f0872a70baa55ecb892e52ea3f5a Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 20:33:52 -0400 Subject: [PATCH 1/7] feat(workflow): emit the instrument flags as a healsparse defect map (#878) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The exposures' instrument flag images — bad columns, saturated pixels, bleed trails — are the one masking input this campaign has that never leaves the pixel domain: exp_split splits them per CCD, SExtractor reads them as IMAFLAGS_ISO, and that is the end of it. The survey footprint is built from the CCD corner WCS in the headers, so it cannot subtract them and would silently include defective pixels. The lost area is percent-level, but it carries exactly the thin small-scale geometry an accurate window function needs, and only ShapePipe ever opens these files. Two rules, both in exposure.smk. exp_defect_map rasterizes ONE exposure's flag splits into a boolean healsparse fragment on the persistent root — nside_sparse 131072 over nside_coverage 128, bit-packed, True = masked, the same form as every map in the UNIONS ladder, so it drops into that ladder without a resolution change. It hangs off exp_split rather than exp_psf, so re-rasterizing at a different fidelity never touches the four-hour PSF chain, and clean_exposure takes its manifest as an input, so reclamation cannot overtake the copy. The WCS comes from the image split beside each flag file (headers only, never pixels): the flag mosaic's HDUs carry no WCS at all — checked on a real exposure — and headers-.npy is a pickled array of astropy WCS instances a second consumer should not inherit. Rasterization samples each flagged CCD pixel across its full extent, corners included, and is CONSERVATIVE: a healpix pixel is masked if any part of the flagged region touches it, because the centre-based alternative erases a one-pixel bad column, which is the geometry #878 exists to keep. `oversample` rides on params with the measured convergence table behind its default of 3. defect_map_merge unions the campaign's fragments into /defect_map_.hsp. It RECONCILES like final_cat_merge — a new exposure is OR-ed in on the spot, an exposure that left the campaign or a fragment that changed forces a rebuild (a union cannot be un-OR-ed), a no-op leaves the file untouched — against a sidecar that records which exposures are already in it, declared as a second output because the two are only meaningful together. Memory is flat in the exposure count: each fragment is read, reduced to its pixel ids and dropped, so the job holds one accumulator (the campaign's footprint) and one 2 MB fragment whether the campaign is 127 exposures or 20k. The map is a campaign PRODUCT, not an input: nothing here reads it back. Measured on 2079612p inside the container, 40 CCDs, 4.4% of pixels flagged: 37 s and 0.62 GB peak RSS at oversample 3, a 2.0 MB fragment of 359858 healpix pixels over 13 coverage pixels (95 s / 1.29 GB / 360840 pixels at oversample 5, so the default is within 0.3% of it at a third of the cost). Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- workflow/Snakefile | 132 +++++++++++++ workflow/config.yaml | 51 +++++ workflow/rules/exposure.smk | 159 ++++++++++++++- workflow/scripts/defect_map_exp.py | 279 +++++++++++++++++++++++++++ workflow/scripts/merge_defect_map.py | 279 +++++++++++++++++++++++++++ 5 files changed, 893 insertions(+), 7 deletions(-) create mode 100644 workflow/scripts/defect_map_exp.py create mode 100644 workflow/scripts/merge_defect_map.py diff --git a/workflow/Snakefile b/workflow/Snakefile index 64640f615..0cde94c5a 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -386,6 +386,8 @@ MERGE_FINAL_HASH = ":".join(( path_hash(Path(workflow.basedir).parent / "scripts" / "python" / "create_final_cat.py"), path_hash(CONFIG_DIR / "final_cat.param"))) +DEFECT_HASH = script_hash("defect_map_exp.py") +MERGE_DEFECT_HASH = script_hash("merge_defect_map.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 @@ -529,6 +531,109 @@ for _entry in PERSIST_EXP: except KeyError as _exc: raise WorkflowError(f"config persist_exp: {_exc.args[0]}") +# --- the instrument-flag defect map (#878) ---------------------------------- +# The exposures' flag images are the campaign's one masking input that never +# leaves the pixel domain, so the survey footprint — built from CCD corner WCS — +# cannot subtract them. exp_defect_map rasterizes them per exposure and +# defect_map_merge unions the fragments; both are argued in +# workflow/scripts/defect_map_exp.py and merge_defect_map.py. +# +# THE RESOLUTION IS THE LADDER'S, not a choice made here: nside_sparse 131072 +# over nside_coverage 128, boolean, True = masked, matching every map under +# `inputs.masks` and every entry in config_tile_Mc.ini's MASK_EXT_PATHS. It is +# config so that a campaign against a different ladder can say so, not so that +# this one can drift. +_DEFECT = dict(config.get("defect_map") or {}) +DEFECT_NSIDE = int(_DEFECT.get("nside", 131072)) +DEFECT_NSIDE_COVERAGE = int(_DEFECT.get("nside_coverage", 128)) +DEFECT_OVERSAMPLE = int(_DEFECT.get("oversample", 3)) +if DEFECT_OVERSAMPLE < 2: + raise WorkflowError( + f"config defect_map.oversample={DEFECT_OVERSAMPLE}: 2 is the geometric " + f"floor (the four corners of a CCD pixel bound its footprint); 1 would " + f"sample pixel centres alone and lose thin defects entirely.") +if DEFECT_NSIDE_COVERAGE >= DEFECT_NSIDE: + raise WorkflowError( + f"config defect_map: nside_coverage={DEFECT_NSIDE_COVERAGE} must be " + f"below nside={DEFECT_NSIDE}.") + + +def defect_map(): + """The campaign's merged defect map. A PRODUCT, not an input: nothing in + this workflow opens it, and config_tile_Mc.ini deliberately does not name + it (that file's header carries the recipe for adding it after a campaign + has produced one).""" + return f"{PRODUCTS_DIR}/defect_map_{CAMPAIGN}.hsp" + + +def defect_map_sidecar(): + """The record of which exposures are already in the map — what makes an + append cheap and a removal correct (merge_defect_map.py argues it).""" + return f"{PRODUCTS_DIR}/defect_map_{CAMPAIGN}.json" + + +def prod_exp_fragment(exp): + """The healsparse fragment exp_defect_map writes. Not a declared output of + anything — the rule declares its manifest, exactly as exp_persist does.""" + return f"{prod_exp_dir(exp)}/defect/defect-{exp}.hsp" + + +@functools.lru_cache(maxsize=1) +def defect_map_exposures(): + """The exposures the union covers: every exposure of TILES_READY. + + SIMPLER THAN THE STAR CATALOGUE'S SET, and the difference is where the two + rules hang. exp_persist waits on exp_psf, so a reclaimed exposure's manifest + cannot be requested without rebuilding a four-hour chain from VOS — hence + star_cat_inputs()'s manifest-or-tar split. exp_defect_map waits on + exp_split, and its manifest lives on the persistent root, so an exposure + that already has one is a DAG leaf and an exposure that does not is one + split away. No special case is needed and none is invented. + + THE IDS, NOT THE PATHS: the fingerprint on `params` must move when the SET + moves and not when a path does. + """ + return sorted({e for t in TILES_READY for e in tile_exposures(t)}) + + +@functools.lru_cache(maxsize=1) +def defect_map_inputs(): + """What defect_map_merge waits for: the fragment manifests of the exposures + whose stores this invocation can still reach. + + An exposure whose scratch store is gone (reclaimed, or purged) and which has + no fragment manifest yet cannot get one without rebuilding its chain from + VOS — the avalanche persist_targets() drops such exposures to avoid. It is + left out of the edge; the job skips it and records it as missing. Every + other exposure's manifest is either already there or one exp_split away. + """ + return [prod_exp_manifest(e, "exp_defect_map") + for e in defect_map_exposures() + if not exp_store_reclaimed(e) + or Path(prod_exp_manifest(e, "exp_defect_map")).exists()] + + +def defect_map_targets(): + """The merged map, whenever this campaign has an exposure to rasterize.""" + if not workflow.is_main_process: + return [] + return [defect_map(), defect_map_sidecar()] if defect_map_inputs() else [] + + +def defect_map_fragment_targets(): + """The per-exposure fragments `rule all` requests DIRECTLY. + + Same argument as persist_targets(): the fragment must leave /scratch whether + or not the campaign ever merges or reclaims, because the flag splits go with + the store at the purge. Reaching them only through defect_map_merge would + make a campaign that never merged lose them. + + HEAD PROCESS ONLY, for the same cost reason as persist_targets(). + """ + if not workflow.is_main_process: + return [] + return defect_map_inputs() + def persist_targets(): """Which exposures this invocation must pack PSF products off scratch for. @@ -833,6 +938,18 @@ STAR_CAT_PRODUCT = _persist.ALWAYS STAR_CAT_PATTERN = _persist.resolve(_persist.ALWAYS) TILE_BYTES_DEFAULT = 46_000_000 +# DEFECT SIDE. The union's footprint is what sizes it, not the exposure count: +# the accumulator is one BIT per sparse pixel of every coverage pixel the +# campaign touches, so a coverage pixel costs (nside/nside_coverage)^2 / 8 bytes +# — 128 KiB at the ladder's 131072/128 — and the loop holds one fragment +# (~2 MB) at a time. Two ways to know the count, in order: the sidecar the last +# merge wrote, which is the campaign's own measured footprint; and before there +# is one, the exposures times the 13 coverage pixels ONE exposure touched +# (measured, 2079612p), which ignores overlap and is therefore an overestimate +# in the safe direction, capped at the whole sky. +DEFECT_MEM_BASE_MB = 800 # interpreter + astropy + healpy + healsparse +DEFECT_COV_PER_EXP = 13 + def _size(path, default): """Bytes on disk, or the documented per-unit default if it is not there.""" @@ -842,6 +959,19 @@ def _size(path, default): return default +def defect_map_cov_bytes(): + """Bytes the union's accumulator occupies: coverage pixels x the bit-packed + block size. See the sizing block above for where the count comes from.""" + block = (DEFECT_NSIDE // DEFECT_NSIDE_COVERAGE) ** 2 // 8 + try: + record = json.loads(Path(defect_map_sidecar()).read_text()) + n_cov = int(record["n_coverage_pixels"]) + except (OSError, ValueError, KeyError): + n_cov = min(DEFECT_COV_PER_EXP * len(defect_map_exposures()), + 12 * DEFECT_NSIDE_COVERAGE ** 2) + return n_cov * block + + def star_cat_max_bytes(): """The LARGEST exposure's psf_validation members — what sizes the merge. @@ -1107,8 +1237,10 @@ rule all: input: [final_cat(t) for t in TILES_READY], persist_targets(), + defect_map_fragment_targets(), star_cat_targets(), final_cat_targets(), + defect_map_targets(), clean_targets(), clean_tile_targets(), diff --git a/workflow/config.yaml b/workflow/config.yaml index f119139be..86581604a 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -180,6 +180,57 @@ outputs: persist_exp: - psf_model +# The instrument-flag DEFECT MAP (CosmoStat/shapepipe#878). The exposures' flag +# images -- bad columns, saturated pixels, bleed trails -- are the one masking +# input this campaign has that never leaves the pixel domain: exp_split splits +# them per CCD and SExtractor reads them as IMAFLAGS_ISO, and that is the end of +# it. The survey footprint is built from the CCD corner WCS in the headers, so +# it cannot subtract them and would silently include defective pixels. So +# `exp_defect_map` rasterizes each exposure's flags into a boolean healsparse +# fragment beside its PSF tar on products_dir, and `defect_map_merge` unions the +# campaign's fragments into /defect_map_.hsp. +# +# IT IS A PRODUCT, NOT AN INPUT. Nothing in the workflow reads it back; it is +# for the footprint, downstream. In particular config_tile_Mc.ini's +# MASK_EXT_PATHS does NOT name it -- that ladder names maps that exist before +# the run starts, and this one exists only after it. That file's header carries +# the recipe for adding it once a campaign has produced one. +# +# Measured on one real exposure (2079612p, 40 CCDs, 4.4% of pixels flagged) on +# this login node inside the campaign container: 37 s at oversample 3, peak RSS +# 0.62 GB, 359858 healpix pixels, a 2.0 MB fragment over 13 coverage pixels +# (21 s / 0.41 GB / 357896 pixels at oversample 2; 95 s / 1.29 GB / 360840 at +# oversample 5 -- so 3 is within 0.3% of 5 over the whole exposure, at a third +# of the cost). ~40 GB and ~20k inodes at DR6 scale, against a ~1 M-inode group +# quota. +# +# EVERY KEY BELOW DEFAULTS TO THE MASK LADDER'S OWN CONVENTION and the block can +# be omitted entirely. nside/nside_coverage match every map under `inputs.masks` +# (boolean, True = masked) so the fragments drop into that ladder without a +# resolution change; they are config so a campaign against a different ladder +# can say so, not so this one can drift. +# +# `oversample` is samples per CCD pixel per axis, endpoints included, so 2 means +# the four corners -- the geometric floor, since a healpix pixel at nside 131072 +# (1.61") is ~74 times a MegaCam pixel (0.187") and cannot hide inside one. +# Raising it only fills boundary and rounding gaps, at linear cost. Measured on +# 2079612p CCD 0 (430886 flagged pixels) against a 12x12 interior reference of +# 14229 healpix pixels: +# +# oversample samples/pixel healpix pixels missed vs reference +# 2 (corners) 4 14081 169 (1.2%) +# 3 9 14198 54 (0.4%) +# 5 25 14254 5 (0.04%) +# +# The rasterization is CONSERVATIVE by design -- a healpix pixel is masked if +# any part of the flagged region touches it -- because the centre-based +# alternative erases a one-pixel bad column, which is the geometry #878 exists +# to keep. +defect_map: + nside: 131072 + nside_coverage: 128 + oversample: 3 + # Rolling exposure-store reclamation (D5). When true, the COMPUTE DAG grows one # `clean_exposure` job per exposure. It fires once every campaign tile that reads # that exposure has its vignets, deletes the exposure's store AND its manifests, diff --git a/workflow/rules/exposure.smk b/workflow/rules/exposure.smk index 5845adb8d..33fd71763 100644 --- a/workflow/rules/exposure.smk +++ b/workflow/rules/exposure.smk @@ -1,24 +1,31 @@ """Exposure chain — per exposure, keyed by exp base id (dedup is structural). exp_get_images -> exp_split -> exp_psf -> exp_persist + `-> exp_defect_map Each in the exposure's own sharded work dir, chained by manifests; every config reads fixed ``$SP_RUN/output/run_sp_exp_*`` INPUT_DIRs, so nothing resolves a run log. There is no `prepare_exposures` aggregation target: these chains hang off the compute DAG (`all` <- final_cat <- tile chain <- exposure manifests). -NO MASK RULE, and that is the design (PR #847). ShapePipe generates no masks. -The only mask that reaches pixels is the instrument flag image delivered with -the exposure, which ``exp_split`` splits per CCD alongside image and weight and -SExtractor reads directly. Sky-fixed masks are healsparse maps, queried once per +NO MASK-GENERATION RULE, and that is the design (PR #847). ShapePipe generates +no masks. The only mask that reaches pixels is the instrument flag image +delivered with the exposure, which ``exp_split`` splits per CCD alongside image +and weight and SExtractor reads directly. ``exp_defect_map`` does not generate +that mask, it EXPORTS it: the flag image is the one masking input the campaign +has that never becomes sky-fixed, so the survey footprint cannot subtract it +(#878), and the rule rasterizes it into the same healsparse form as every other +mask here. Nothing in this workflow reads the result back. + +Sky-fixed masks are healsparse maps, queried once per object: ``mask_query`` (inside exp_psf's config chain) writes ``FLAG_EXT`` onto each CCD's SExtractor catalogue for setools' star cut, and ``make_cat`` writes the per-band ``MASK_`` columns on the tile side. Neither needs a rule, a star catalogue, or a network fetch — hence no ``star_catalogue`` / ``exp_star_cat`` here, and no ``exp_mask``. -``exp_persist`` is the one rule here that writes to the PERSISTENT root: it -packs the PSF products named by `persist_exp:` into one tar per exposure off +``exp_persist`` is one of the two rules here that write to the PERSISTENT root +(``exp_defect_map`` is the other): it packs the PSF products named by `persist_exp:` into one tar per exposure off /scratch before the purge (or clean_exposure) can take them. It is a separate rule from exp_psf precisely so that editing that list costs a re-pack and not a four-hour refit; the full @@ -169,6 +176,75 @@ rule exp_persist: " {params.patterns}" +# --- the pixel-domain mask leaves the pixel domain (#878) ------------------- +# The second rule here that writes to the PERSISTENT root, and the second one +# clean_exposure must wait for. It rasterizes this exposure's per-CCD instrument +# flag splits — the ONE masking input the campaign has that never becomes a +# sky-fixed map — into a boolean healsparse fragment at the mask ladder's own +# resolution, so the footprint can finally subtract bad columns, saturated +# pixels and bleed trails. The full argument, the WCS source and the measured +# oversampling table are in workflow/scripts/defect_map_exp.py. +# +# AFTER exp_split, NOT AFTER exp_psf, and it is deliberately not chained behind +# the PSF work: the flag splits exist the moment the split finishes, and hanging +# a two-minute rasterization off a four-hour rule would make re-rasterizing the +# campaign at a different oversampling cost the PSF chain. It runs in parallel +# with exp_psf, and both are ordered before clean_exposure. +# +# ONE DECLARED OUTPUT, AND IT IS A MANIFEST, exactly as exp_persist: the +# fragment is written beside it on the persistent root and the manifest records +# the per-CCD healpix counts. Byte-stable, so a no-op rerun does not move the +# mtime clean_exposure reads. +# +# NOT A LOCALRULE, and this is where it parts company with exp_persist. That +# rule is a tar of a few MB — seconds, far shorter than the scheduling latency. +# This one is 40 CCDs of WCS transforms and ang2pix over ~17 M flagged pixels: +# 37 s measured end to end on a real exposure at oversample 3 (2079612p, on this +# login node, inside the campaign container; 21 s at oversample 2). That is real +# work, it is CPU-bound, and running ~20k of them under local-cores in the head +# process would serialise the campaign behind them. +# +# THE RESOLUTION AND THE OVERSAMPLING RIDE ON params. Both are what the fragment +# IS, and both are decisions that may be revisited; on params they re-rasterize +# (a minute) and leave the split and the PSF chain alone. +rule exp_defect_map: + input: + rules.exp_split.output.manifest + output: + manifest = f"{PROD_EXP_DIR}/manifests/exp_defect_map.json" + # No `log:`, for exp_persist's reason: the failure modes are "no flag split + # under the store" and "a flag split with no image beside it", both reported + # on stderr, neither with a per-CCD verdict worth a completeness record. + params: + exp_dir = lambda wc: exp_dir(wc.exp), + dest = lambda wc: f"{prod_exp_dir(wc.exp)}/defect", + nside = DEFECT_NSIDE, + nside_cov = DEFECT_NSIDE_COVERAGE, + oversample = DEFECT_OVERSAMPLE, + script_hash = DEFECT_HASH + threads: 1 + retries: 2 + resources: + # Measured on 2079612p inside the container: peak RSS 0.62 GB at + # oversample 3 (0.41 GB at 2), dominated by one CCD's sample arrays + # (430k flagged pixels x 9 samples x two float64 coordinate arrays) + # plus the bit-packed fragment's 13 coverage pixels. Flat in the number + # of CCDs — they are rasterized one at a time — so this is sized on the + # worst CCD, not on the exposure, at ~3x it. It scales with + # `oversample`, which is why that knob is not free. + mem_mb = lambda wc, attempt: 2000 * attempt, + # 37 s measured; a factor of ~30 for a dirtier exposure, a higher + # oversampling and a busy filesystem. + runtime = 20 + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/defect_map_exp.py" + " --exp-dir '{params.exp_dir}' --exp {wildcards.exp}" + " --dest '{params.dest}' --manifest {output.manifest}" + " --nside {params.nside} --nside-coverage {params.nside_cov}" + " --oversample {params.oversample}" + + # --- reclamation (D5) ------------------------------------------------------- # The one exception to "no reclamation in this file": clean_exposure OWNS # exposure-level deletion, and it is a real job, not temp() bookkeeping, because @@ -207,7 +283,8 @@ rule clean_exposure: # exposure's chain did not already put there. It is UNCONDITIONAL now: # exp_persist always packs the star catalogue's inputs, so there is no # keep list under which this rule has nothing to wait for. - lambda wc: [prod_exp_manifest(wc.exp, "exp_persist")] + lambda wc: [prod_exp_manifest(wc.exp, "exp_persist"), + prod_exp_manifest(wc.exp, "exp_defect_map")] output: tombstone = f"{EXP_DIR}/cleaned.json" params: @@ -302,3 +379,71 @@ rule star_cat_merge: " --tile-list '{params.tile_list}' --index-db '{params.index_db}'" " --output {output.star_cat}" " --campaign '{params.campaign}'" + + +# --- the campaign's defect map (#878) --------------------------------------- +# ONE job per campaign: every exposure's fragment, OR-ed into +# `/defect_map_.hsp`. It is the exposure side's third +# campaign product, and the only one nothing downstream in this workflow opens: +# it exists so the survey footprint can subtract the pixel-domain masking the +# CCD corner WCS cannot see. merge_defect_map.py argues the reconcile, the +# memory-flat accumulation and why the map is a campaign PRODUCT rather than an +# entry in config_tile_Mc.ini's MASK_EXT_PATHS. +# +# TWO DECLARED OUTPUTS, the map and the SIDECAR that records which exposures are +# already in it. They are written together and they are only meaningful +# together: a union cannot be un-OR-ed, so the record of what went in is what +# makes an append cheap and a removal correct. Declaring both means a failed job +# takes both, and the next invocation rebuilds rather than reconciling against a +# record that does not describe the map beside it. +# +# THE INPUT IS THE FRAGMENT MANIFESTS, not the fragments: the manifest is what +# exp_defect_map declares, so it is the edge that orders this after the +# rasterizations. Reclaimed exposures need no special case here — unlike the +# star catalogue's tars, the fragment manifest lives on the persistent root and +# survives reclamation, and the rule that writes it hangs off exp_split, not off +# a store that reclamation took. An exposure cleaned by a workflow that predates +# this rule simply has no fragment; the job skips it and says so. +# +# THE PATHS DO NOT REACH THE SHELL (~20k of them at DR6 scale, an order of +# magnitude over MAX_ARG_STRLEN): the job is handed the tile list and the index +# and derives the same set, with the set's FINGERPRINT on `params` as the rerun +# trigger. Same discipline as the two merges above. +# +# NOT A LOCALRULE: the accumulator is the campaign's footprint at nside 131072, +# ~3 GB resident at DR6 scale. +rule defect_map_merge: + input: + lambda wc: defect_map_inputs() + output: + defect_map = defect_map(), + sidecar = defect_map_sidecar() + params: + products_dir = str(PRODUCTS_DIR), + tile_list = str(config["tile_list"]), + index_db = str(INDEX_DB), + nside = DEFECT_NSIDE, + nside_cov = DEFECT_NSIDE_COVERAGE, + inputs = unit_fingerprint(defect_map_exposures()), + script_hash = MERGE_DEFECT_HASH + threads: 1 + resources: + # Sized on the FOOTPRINT, not on the exposure count: the accumulator is + # one bit per sparse pixel of every touched coverage pixel, and the loop + # holds one fragment at a time (merge_defect_map.py's memory argument). + # DEFECT_COV_BYTES is that arithmetic, from the sidecar's own recorded + # coverage count once there is one and from a measured per-exposure + # figure before that; the Snakefile's sizing block carries both. + mem_mb = lambda wc, attempt: attempt * ( + DEFECT_MEM_BASE_MB + 2 * defect_map_cov_bytes() // 1_000_000), + # Dominated by reading fragments (~2 MB each) and setting their pixels; + # ~1 s per exposure measured, over a floor that covers writing the map. + runtime = lambda wc, attempt: attempt * ( + 20 + len(defect_map_exposures()) // 30) + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/merge_defect_map.py" + " --products-dir '{params.products_dir}'" + " --tile-list '{params.tile_list}' --index-db '{params.index_db}'" + " --output {output.defect_map} --sidecar {output.sidecar}" + " --nside {params.nside} --nside-coverage {params.nside_cov}" diff --git a/workflow/scripts/defect_map_exp.py b/workflow/scripts/defect_map_exp.py new file mode 100644 index 000000000..17a9c767e --- /dev/null +++ b/workflow/scripts/defect_map_exp.py @@ -0,0 +1,279 @@ +#!/usr/bin/env python3 +"""Rasterize ONE exposure's instrument flags into a boolean healsparse fragment. + +Run as the shell of the in-DAG ``exp_defect_map`` rule, never by hand. + +WHY THIS EXISTS (CosmoStat/shapepipe#878). The survey footprint is built as a +positive coverage map minus the healsparse mask bits, and every masking input +the campaign has is already sky-fixed and queried per object — except one. The +instrument flag image delivered with each exposure (bad columns, saturated +pixels, bleed trails) never leaves the PIXEL domain: ``exp_split`` splits it per +CCD, SExtractor reads it as ``IMAFLAGS_ISO``, and that is the end of it. The +coverage map is built from the CCD corner WCS in the headers, so it cannot +subtract those pixels and the footprint silently includes them. The lost area is +percent-level, but it is exactly the thin, small-scale geometry an accurate +window function needs. Only ShapePipe ever opens these files, so the map has to +come from here. + +WHAT IT WRITES. ``/defect-.hsp``: a boolean ``HealSparseMap``, +``nside_sparse`` 131072 over ``nside_coverage`` 128, ``bit_packed``, ``True`` +where any flag bit is set. Those are not free choices — they are the +convention every other map in this campaign's ladder already follows (the +UNIONS ugriz bit maps under ``inputs.masks``, and ``config_tile_Mc.ini``'s +``MASK_EXT_PATHS``), so the fragments and the map they union into drop into that +ladder without a resolution change. ``True = masked`` likewise. + +ANY BIT, NOT A BIT TABLE. The flag image is a bitmask, but the CFIS flags are +"this pixel is not to be trusted" in several flavours (the campaign's own +exposures carry values 1, 2, 3, 8 and 11), and nothing downstream distinguishes +them: SExtractor's ``IMAFLAGS_ISO`` cut is nonzero-vs-zero. So the fragment is +``!= 0`` and the map is boolean. A per-bit ladder would be a different product +answering a question nobody has asked yet. + +WHERE THE WCS COMES FROM, AND WHY NOT FROM THE FLAG FILE. The flag mosaic's +per-CCD HDUs carry NO WCS at all — checked on a real exposure: ``CTYPE`` empty, +``CRVAL`` 0, the identity transform. Only the image mosaic is astrometric, so +``split_exp`` saves headers for the image suffix alone. The fragment therefore +takes each CCD's WCS from its ``image--.fits`` split, which sits +beside the flag split and carries the full SCAMP header (``RA---TAN`` with PV +distortion). The image PIXELS are never read: only ``fits.getheader``. + +NOT ``headers-.npy``, which is the other thing ``exp_split`` writes and +would be one small file instead of forty header reads. It is a pickled object +array of ``astropy.wcs.WCS`` INSTANCES, so reading it unpickles astropy objects +across whatever container rebuild happens next; the FITS headers are +self-describing text and cost milliseconds. ``merge_headers`` may live with the +pickle because it is the file's author's own consumer; a second consumer should +not inherit the coupling. + +HOW A PIXEL IS RASTERIZED, AND WHY IT IS SAMPLED RATHER THAN INTEGRATED. A +healpix pixel at nside 131072 is 1.61 arcsec across; a MegaCam pixel is 0.187 +arcsec. So one healpix pixel covers ~74 CCD pixels, and the question is never +"which healpix pixels does this CCD pixel cover" but "which healpix pixels does +the flagged REGION touch". The answer is taken by sampling each flagged CCD +pixel on an ``oversample`` x ``oversample`` grid spanning its full extent +(``linspace(-0.5, 0.5)``, so the four corners are always sampled), converting to +sky through that CCD's WCS and binning with ``ang2pix``. + +THE RASTERIZATION IS CONSERVATIVE, deliberately: a healpix pixel is masked if +any part of the flagged region falls in it. The alternative — mask when the +healpix pixel's CENTRE is flagged — erases a one-pixel bad column entirely, +which is precisely the geometry #878 exists to keep. The cost is that a thin +defect is widened to the healpix resolution; at 1.61 arcsec that is the price +of the ladder's nside. + +WHAT ``oversample`` BUYS, MEASURED (exposure 2079612p, CCD 0, 430886 flagged +pixels; the reference is a 12x12 interior grid, 14229 healpix pixels): + + oversample samples/pixel healpix pixels missed vs reference + 2 (corners) 4 14081 169 (1.2%) + 3 9 14198 54 (0.4%) + 5 25 14254 5 (0.04%) + +and the cost is linear in the sample count. 3 is the default in ``config.yaml`` +for that reason, and it rides on the rule's ``params`` so raising it re-rasterizes +without touching anything upstream. 2 is the geometric floor: the four corners +of a CCD pixel bound its footprint, and a healpix pixel 74 times its area cannot +sit inside it — the residual 1.2% is boundary and rounding, not a class of +missed defect. + +BYTE-STABLE, tmp-then-``cmp``-then-``mv``, the pattern ``persist_exp`` uses: +healsparse's FITS output carries no timestamp (checked), so re-rasterizing an +unchanged store produces an identical file and leaves its mtime alone. mtime is +a rerun trigger and ``clean_exposure`` waits on this rule's manifest, so an +unconditional rewrite would make every reclamation look out of date once per +invocation. + +THE MANIFEST IS THE ONLY DECLARED OUTPUT and it lives on the PERSISTENT root +beside the fragment (``/exp///manifests/``), not in +the exposure's scratch ``manifests/`` which ``clean_exposure`` deletes wholesale +— same placement, and same reason, as ``exp_persist``. It records the per-CCD +healpix counts, so a reader can see which CCD contributed what without opening +the map. + +ORDERED BEFORE RECLAMATION. ``clean_exposure`` takes this manifest as an input, +exactly as it takes ``exp_persist``'s: the flag splits live on /scratch and go +with the store, so the fragment must be on /project before anything is deleted. +""" + +import argparse +import filecmp +import json +import re +import sys +import warnings +from pathlib import Path + +import numpy as np +from astropy.io import fits +from astropy.wcs import WCS + +import healpy as hp +import healsparse as hsp + +# The split stage's run dir (RUN_NAME in config_exp_Sp.ini) and its module. +# Hardcoded for the same reason persist_exp.py hardcodes its own: this rule +# rasterizes the SPLIT stage's flag images and nothing else, and a knob here +# would be a knob for "rasterize some other stage". +RUN_NAME = "run_sp_exp_Sp" +MODULE = "split_exp_runner" + +# `flag-2079612-13.fits` -> ccd 13. The number string is the exposure's +# ($SP_UNIT_NUM, `-2079612`), so the CCD is what follows the last dash. +_CCD = re.compile(r"-(\d+)\.fits$") + + +def split_dir(exp_dir: Path) -> Path: + return exp_dir / "output" / RUN_NAME / MODULE / "output" + + +def ccd_files(exp_dir: Path) -> list: + """``(ccd, flag path, image path)`` for every CCD this exposure split, in + CCD order. + + Driven by the FLAG files, which are what we rasterize; the image split is + looked up beside each one for its header alone. A flag split with no image + beside it is a hard error rather than a skip — ``exp_split`` writes the + three suffixes in one pass, so a missing image means the store is not what + we think it is, and a fragment quietly short of a CCD would be a hole in the + footprint nothing downstream could notice. + """ + root = split_dir(exp_dir) + out = [] + for flag in sorted(root.glob("flag-*.fits")): + match = _CCD.search(flag.name) + if not match: + continue + image = flag.with_name(flag.name.replace("flag-", "image-", 1)) + if not image.exists(): + sys.exit(f"defect_map_exp: {flag.name} has no {image.name} beside " + f"it in {root}; the WCS lives on the image split (see the " + f"module docstring)") + out.append((int(match.group(1)), flag, image)) + return sorted(out) + + +def offsets(oversample: int) -> tuple: + """Sample offsets within one CCD pixel, in pixel units. + + ``linspace`` with both endpoints, so the CORNERS are always sampled: they + are what bounds the pixel's footprint, and the interior samples only fill + boundary and rounding gaps (the docstring's table measures how many). + """ + if oversample < 2: + sys.exit(f"defect_map_exp: oversample={oversample} would sample the " + f"pixel centre alone and lose the pixel's extent; 2 is the " + f"geometric floor (its four corners)") + step = np.linspace(-0.5, 0.5, oversample) + grid_x, grid_y = np.meshgrid(step, step) + return grid_x.ravel(), grid_y.ravel() + + +def rasterize_ccd(flag_path: Path, image_path: Path, nside: int, + off_x, off_y) -> np.ndarray: + """The healpix pixel ids (NEST, ``nside``) this CCD's flags touch. + + One CCD at a time and one array at a time: the whole point of the loop in + ``main`` is that the job's footprint is one CCD's samples, not the + exposure's. + """ + with warnings.catch_warnings(): + # SCAMP headers carry a deprecated RADECSYS and a redundant SIP block + # beside the PV distortion astropy actually uses; both are FITSFixedWarning + # noise on every one of 40 CCDs and neither changes the transform. + warnings.simplefilter("ignore") + wcs = WCS(fits.getheader(image_path)) + data = fits.getdata(flag_path) + rows, cols = np.nonzero(data) + if rows.size == 0: + return np.empty(0, dtype=np.int64) + # 1-based FITS pixel coordinates, sampled across each flagged pixel's extent. + x = (cols[:, None] + 1.0 + off_x[None, :]).ravel() + y = (rows[:, None] + 1.0 + off_y[None, :]).ravel() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + ra, dec = wcs.all_pix2world(x, y, 1) + return np.unique(hp.ang2pix(nside, ra, dec, lonlat=True, nest=True)) + + +def write_stable(tmp: Path, dest: Path) -> None: + """Move ``tmp`` onto ``dest``, or drop it when the bytes already match.""" + if dest.exists() and filecmp.cmp(tmp, dest, shallow=False): + tmp.unlink() + else: + tmp.replace(dest) + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--exp-dir", required=True, type=Path, + help="the exposure's scratch store") + parser.add_argument("--exp", required=True) + parser.add_argument("--dest", required=True, type=Path, + help="/exp///defect; the " + "fragment is /defect-.hsp") + parser.add_argument("--manifest", required=True, type=Path) + parser.add_argument("--nside", type=int, default=131072, + help="nside_sparse; the mask ladder's resolution") + parser.add_argument("--nside-coverage", type=int, default=128) + parser.add_argument("--oversample", type=int, default=3, + help="samples per CCD pixel per axis (see the module " + "docstring's measured table)") + args = parser.parse_args() + + ccds = ccd_files(args.exp_dir) + if not ccds: + sys.exit(f"defect_map_exp: {args.exp}: no flag split under " + f"{split_dir(args.exp_dir)}") + + off_x, off_y = offsets(args.oversample) + fragment = hsp.HealSparseMap.make_empty( + args.nside_coverage, args.nside, np.bool_, bit_packed=True) + per_ccd = {} + for ccd, flag_path, image_path in ccds: + pixels = rasterize_ccd(flag_path, image_path, args.nside, off_x, off_y) + per_ccd[str(ccd)] = int(pixels.size) + if pixels.size: + fragment[pixels] = True + + args.dest.mkdir(parents=True, exist_ok=True) + frag_path = args.dest / f"defect-{args.exp}.hsp" + tmp = frag_path.with_name(frag_path.name + ".tmp") + try: + fragment.write(str(tmp), clobber=True) + write_stable(tmp, frag_path) + finally: + tmp.unlink(missing_ok=True) + + body = { + "stage": "exp_defect_map", "level": "exp", "unit": args.exp, + "status": "complete", + "map": str(frag_path), + "nside": args.nside, + "nside_coverage": args.nside_coverage, + "oversample": args.oversample, + "n_ccds": len(ccds), + # Per CCD, so a reader can see WHICH CCD contributed what without + # opening the map — a CCD at zero is a real thing (a clean chip) and a + # whole exposure at zero is not. + "pixels_per_ccd": per_ccd, + "n_pixels": int(fragment.n_valid), + "n_coverage_pixels": int(fragment.coverage_mask.sum()), + "bytes": frag_path.stat().st_size, + } + args.manifest.parent.mkdir(parents=True, exist_ok=True) + tmp = args.manifest.with_name(args.manifest.name + ".tmp") + try: + tmp.write_text(json.dumps(body, indent=2, sort_keys=True) + "\n") + write_stable(tmp, args.manifest) + finally: + tmp.unlink(missing_ok=True) + + print(f"[defect_map_exp] {args.exp}: {len(ccds)} CCD(s), " + f"{body['n_pixels']} healpix pixel(s) over " + f"{body['n_coverage_pixels']} coverage pixel(s), " + f"{body['bytes'] / 1e6:.1f} MB -> {frag_path}") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/merge_defect_map.py b/workflow/scripts/merge_defect_map.py new file mode 100644 index 000000000..dab1dd219 --- /dev/null +++ b/workflow/scripts/merge_defect_map.py @@ -0,0 +1,279 @@ +#!/usr/bin/env python3 +"""Union the campaign's per-exposure defect fragments into ONE healsparse map. + +Run as the shell of the campaign-level ``defect_map_merge`` rule, never by hand. + +WHAT IT PRODUCES, AND FOR WHOM. ``/defect_map_.hsp``: +a boolean ``HealSparseMap``, ``True`` where any exposure of the campaign flagged +the sky, at the mask ladder's own resolution (nside_sparse 131072 over +nside_coverage 128, ``bit_packed``). It is the pixel-domain half of the masking +the footprint could not otherwise see (CosmoStat/shapepipe#878; +``defect_map_exp.py`` argues why the map has to come from here at all). + +IT IS A CAMPAIGN PRODUCT, NOT AN INPUT. Nothing in this workflow consumes it: +``config_tile_Mc.ini``'s ``MASK_EXT_PATHS`` names maps that exist before the run +starts, and this one exists only after it. Adding it to that ladder is a +deliberate, later config edit against a path a campaign has actually produced — +that file's header carries the recipe, and deliberately not the entry. + +TRUE MEANS MASKED, matching every other map in the ladder. A position outside +the campaign's coverage reads ``False``, indistinguishable from clean, exactly +as the UNIONS bit maps behave outside theirs; the union's coverage is the +campaign's exposures and nothing else says where that is. + +IT RECONCILES, AND THE ASYMMETRY IS THE WHOLE DESIGN. The output must be a +function of the input set — that is what makes the rule's fingerprint mean +anything — but a union is not invertible: + + * an exposure whose fragment is NEW is OR-ed into the map on the spot. This is + the common case (a campaign grows by appending tiles) and it reads exactly + the appended exposures; + * an exposure that LEFT the campaign, or whose fragment CHANGED, forces a + REBUILD from every fragment, because nothing can un-OR a pixel that two + exposures both set. Rebuilding is honest about that rather than leaving a + stale bit nothing would ever notice; + * a plan with neither leaves the file UNTOUCHED — not rewritten identically, + untouched, so its mtime cannot move. mtime is a rerun trigger. + +Which exposures are already in the map is recorded in a SIDECAR beside it +(``defect_map_.json``), each with its fragment's size and mtime — the +same "did this change since we read it" stamp ``merge_final_cat.py`` records on +each hdf5 dataset, and for the same reason. The map itself cannot carry that +record: a healsparse FITS header is no place for twenty thousand exposures. The +sidecar is therefore a DECLARED OUTPUT of the rule alongside the map; losing one +without the other would be a map nobody can reconcile, and snakemake removing +both on a failure is the correct recovery (the next run rebuilds). + +MEMORY IS FLAT IN THE NUMBER OF EXPOSURES, which is the reason for the loop +below and not for a list comprehension over ``HealSparseMap.read``. Fragments +are never held together and are never OR-ed as maps: each is read, reduced to +its ``valid_pixels`` (~360k int64, ~3 MB on a real exposure), set into the +accumulator, and dropped. The job's footprint is therefore ONE accumulator plus +ONE fragment, whether the campaign is 127 exposures or 20k. The accumulator is +a function of the campaign's FOOTPRINT, not of its exposure count: a coverage +pixel costs (nside/nside_coverage)^2 / 8 bytes bit-packed — 128 KiB at the +ladder's resolution — so a DR6-scale footprint (~23k coverage pixels, measured +on the UNIONS ugriz maps) is ~3 GB resident and the rule is sized on exactly +that count. + +WHICH EXPOSURES — AND WHY THE JOB DERIVES THE SET. The campaign's: every +exposure of every tile both declared in ``tile_list`` and present in the index, +through ``build_index.campaign_exposures`` so there is one definition and not +two that can drift. It is derived rather than passed because at DR6 scale the +set is ~20k paths and a shell command reaches ``execve`` as a single argv entry +capped at 128 KiB; the rule's ``input`` is the DAG edge and its ``params`` +carries a fingerprint of the same ids. + +A CAMPAIGN EXPOSURE WITH NO FRAGMENT IS SKIPPED, NOT AN ERROR, and that is the +one place this differs from ``merge_final_cat``. Fragments accumulate on the +persistent root across campaigns and survive reclamation, but an exposure +reclaimed by a workflow PREDATING this rule has none and never will without a +rebuild from VOS. Failing would make the map unbuildable for exactly the +campaigns that most want it; the count of skipped exposures is reported and +recorded in the sidecar instead. +""" + +import argparse +import filecmp +import json +import sys +from pathlib import Path + +import numpy as np + +import healsparse as hsp + +# Same directory; the rule invokes this file by path, so it is sys.path[0]. +import build_index + + +def fragment_path(products_dir: Path, exp: str) -> Path: + """Where ``exp_defect_map`` wrote this exposure's fragment.""" + return (products_dir / "exp" / exp[:2] / exp / "defect" + / f"defect-{exp}.hsp") + + +def stamp(path: Path) -> list: + """The fragment's identity, as recorded in the sidecar. + + Size and mtime, not a checksum: the fragment is ~2 MB, it is written + byte-stably (so a no-op re-rasterization does not move its mtime), and the + question is only "did this change since we read it". + """ + st = path.stat() + return [st.st_size, st.st_mtime_ns] + + +def fragments(products_dir: Path, tile_list: Path, index_db: Path) -> tuple: + """``({exp: fragment path}, [exposures with no fragment])``, in exposure order.""" + have, missing = {}, [] + for exp in build_index.campaign_exposures(tile_list, index_db): + path = fragment_path(products_dir, exp) + if path.exists(): + have[exp] = path + else: + missing.append(exp) + return have, missing + + +class Plan: + """What reconciling this campaign into this map requires. + + ``append`` is the cheap path — OR these fragments into the map on disk. + ``rebuild`` is the honest one: a union cannot drop a pixel, so a removal or + a changed fragment means reading every fragment again. + """ + + def __init__(self, append, rebuild, reason): + self.append, self.rebuild, self.reason = append, rebuild, reason + + def empty(self): + return not (self.append or self.rebuild) + + def describe(self): + if self.rebuild: + return f"rebuilt from {len(self.rebuild)} fragment(s) ({self.reason})" + return f"{len(self.append)} fragment(s) appended" + + +def read_sidecar(path: Path) -> dict: + try: + return json.loads(path.read_text()) + except (OSError, ValueError): + return {} + + +def reconcile_plan(output: Path, sidecar: Path, have: dict) -> Plan: + """Compare what is on disk with the campaign, WITHOUT writing anything. + + A missing map, or a sidecar that does not describe it, is a rebuild: the two + are written together and either one alone is not evidence about the other. + """ + record = read_sidecar(sidecar) + known = record.get("exposures") or {} + if not output.exists() or not known: + return Plan([], sorted(have), "no map on disk") + + gone = sorted(set(known) - set(have)) + if gone: + return Plan([], sorted(have), + f"{len(gone)} exposure(s) left the campaign") + changed = sorted(exp for exp, path in have.items() + if exp in known and list(known[exp]) != stamp(path)) + if changed: + return Plan([], sorted(have), + f"{len(changed)} fragment(s) changed on disk") + return Plan(sorted(set(have) - set(known)), [], "") + + +def accumulate(target, paths, nside_coverage, nside) -> None: + """OR each fragment into ``target``, ONE AT A TIME (see the docstring). + + A fragment at the wrong resolution is a hard error rather than a silent + upgrade: the ladder's nside is a campaign-wide decision, and a fragment that + disagrees with it was rasterized by a differently-configured run. + """ + for exp, path in paths: + frag = hsp.HealSparseMap.read(str(path)) + if (frag.nside_sparse, frag.nside_coverage) != (nside, nside_coverage): + sys.exit(f"merge_defect_map: {path} is nside_sparse " + f"{frag.nside_sparse} / nside_coverage " + f"{frag.nside_coverage}, not {nside} / {nside_coverage}; " + f"re-rasterize {exp} before merging") + pixels = frag.valid_pixels + del frag + if pixels.size: + target[pixels] = True + del pixels + + +def write_stable(tmp: Path, dest: Path) -> None: + if dest.exists() and filecmp.cmp(tmp, dest, shallow=False): + tmp.unlink() + else: + tmp.replace(dest) + + +def apply_plan(output: Path, sidecar: Path, plan: Plan, have: dict, + missing: list, nside_coverage: int, nside: int) -> dict: + """Carry the plan out on tmp copies, then move both files into place. + + Map and sidecar are moved together at the end, so a crash mid-merge leaves + the previous PAIR intact rather than a map the record no longer describes. + """ + if plan.rebuild: + target = hsp.HealSparseMap.make_empty( + nside_coverage, nside, np.bool_, bit_packed=True) + todo = plan.rebuild + else: + target = hsp.HealSparseMap.read(str(output)) + todo = plan.append + accumulate(target, [(exp, have[exp]) for exp in todo], + nside_coverage, nside) + + record = { + "campaign_exposures": len(have) + len(missing), + "nside": nside, + "nside_coverage": nside_coverage, + "n_pixels": int(target.n_valid), + "n_coverage_pixels": int(target.coverage_mask.sum()), + # What the NEXT invocation reconciles against; sorted so the sidecar is + # byte-stable for a given campaign state. + "exposures": {exp: stamp(path) for exp, path in sorted(have.items())}, + # Recorded rather than merely printed: a map short of exposures should + # say so on disk, not only in a job log nobody keeps. + "exposures_without_fragment": missing, + } + + map_tmp = output.with_name(output.name + ".tmp") + side_tmp = sidecar.with_name(sidecar.name + ".tmp") + try: + target.write(str(map_tmp), clobber=True) + side_tmp.write_text(json.dumps(record, indent=2, sort_keys=True) + "\n") + write_stable(map_tmp, output) + write_stable(side_tmp, sidecar) + finally: + map_tmp.unlink(missing_ok=True) + side_tmp.unlink(missing_ok=True) + return record + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--products-dir", required=True, type=Path, + help="the persistent root; fragments are found beneath it") + parser.add_argument("--tile-list", required=True, type=Path, + help="the campaign's tile list (config tile_list)") + parser.add_argument("--index-db", required=True, type=Path, + help="the campaign's run index (config outputs.index_db)") + parser.add_argument("--output", required=True, type=Path) + parser.add_argument("--sidecar", required=True, type=Path, + help="the reconciliation record written beside the map") + parser.add_argument("--nside", type=int, default=131072) + parser.add_argument("--nside-coverage", type=int, default=128) + args = parser.parse_args() + + have, missing = fragments(args.products_dir, args.tile_list, args.index_db) + if not have: + # An empty map would satisfy every downstream existence check and mask + # nothing anywhere. + sys.exit(f"merge_defect_map: no campaign exposure in {args.tile_list} " + f"has a defect fragment under {args.products_dir}") + + args.output.parent.mkdir(parents=True, exist_ok=True) + plan = reconcile_plan(args.output, args.sidecar, have) + if plan.empty(): + print(f"[merge_defect_map] unchanged: {args.output} " + f"({len(have)} exposure(s))") + return + record = apply_plan(args.output, args.sidecar, plan, have, missing, + args.nside_coverage, args.nside) + warn = (f"; {len(missing)} campaign exposure(s) have no fragment" + if missing else "") + print(f"[merge_defect_map] {plan.describe()} -> {args.output} " + f"({record['n_pixels']} healpix pixel(s) over " + f"{record['n_coverage_pixels']} coverage pixel(s)){warn}") + + +if __name__ == "__main__": + main() From 6783f4de8ebb143c36cb30e898655f5fea9a798d Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 20:34:01 -0400 Subject: [PATCH 2/7] docs(workflow): the defect map in the README, and how to join it to the ladder MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit config_tile_Mc.ini does NOT gain a MASK_EXT_PATHS entry, and the header now says why: every path in that ladder exists before the run starts, and the defect map exists only after it, so wiring it would name a file the first run of a fresh campaign cannot have and every tile would fail on a missing map. What the header gains instead is the recipe — where the map lands, copying or symlinking it into the inputs.masks root, the one MASK_EXT_PATHS entry and the matching final_cat.param line — plus the two things to know before doing it: the map is the producing campaign's own footprint, so a catalogue built from a different tile list reads False (unknown, not clean) outside it, and the rasterization is conservative, widening a one-pixel bad column to the 1.61" healpix resolution. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- workflow/README.md | 37 ++++++++++++++++++++++++- workflow/config/cfis/config_tile_Mc.ini | 28 +++++++++++++++++++ 2 files changed, 64 insertions(+), 1 deletion(-) diff --git a/workflow/README.md b/workflow/README.md index 6b21ef097..c9b176cd4 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -158,7 +158,7 @@ workflow/ bin/sp committed launcher (module load + /project venv + launch code snapshot + run/report/container/cancel) rules/ prepare.smk tile get_images/uncompress/find_exposures - exposure.smk per-exposure: get_images, split, psf, persist (no temp()); campaign star_cat_merge + exposure.smk per-exposure: get_images, split, psf, persist, defect_map (no temp()); campaign star_cat_merge, defect_map_merge tile.smk per-tile: exp forest, merge_headers, detect, vignets, ngmix, merge, make_cat; campaign final_cat_merge scripts/ build_index.py prepare-phase run_index.sqlite builder (plain script) @@ -171,6 +171,8 @@ workflow/ merge_star_cat.py ALL exposures' validation_psf, out of the tars -> full_starcat_.hdf5 merge_final_cat.py ALL tiles' final_cat -> final_cat_.hdf5 (the final_cat_merge rule) clean_exposure.py ONE exposure's store + manifests + logs -> tombstone (the clean_exposure rule) + defect_map_exp.py ONE exposure's per-CCD instrument flags -> a boolean healsparse fragment (the exp_defect_map rule) + merge_defect_map.py ALL exposures' fragments -> defect_map_.hsp (the defect_map_merge rule) profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; keep-going ``` @@ -221,6 +223,36 @@ profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; kee catalogue server, staged, or rasterized, which is why the old `star_catalogue` / `exp_star_cat` / `exp_mask` rules and their cache root are gone. +- **The one pixel-domain mask now leaves the pixel domain.** The instrument + flag image is the exception to everything above: bad columns, saturated + pixels and bleed trails, split per CCD by `exp_split` and read by SExtractor + as `IMAFLAGS_ISO`, and never sky-fixed. The survey footprint is built from the + CCD corner WCS in the headers, so it cannot subtract them — the footprint + would silently include defective pixels, and the lost area, though only + percent-level, carries exactly the thin small-scale geometry an accurate + window function needs ([#878](https://github.com/CosmoStat/shapepipe/issues/878)). + So `exp_defect_map` rasterizes each exposure's flags into a boolean healsparse + fragment on the persistent root, and `defect_map_merge` unions the campaign's + fragments into `/defect_map_.hsp`. Same form as every + other map here — nside 131072 over coverage 128, `True` = masked — so it drops + into the ladder unchanged. **It is a product, not an input:** nothing in the + workflow reads it back, and `config_tile_Mc.ini` deliberately does not name it + in `MASK_EXT_PATHS`, because that ladder names maps that exist before the run + and this one exists only after it; that file's header carries the recipe for + adding it once a campaign has produced one. The rule hangs off `exp_split`, + not off `exp_psf`, so re-rasterizing the campaign at a different fidelity + never touches the PSF chain, and `clean_exposure` takes its manifest as an + input, so reclamation cannot overtake the copy. Its resolution and its + oversampling ride on `params`; `config.yaml`'s `defect_map:` block carries + both, the measured convergence table behind the default, and the measurement + on one real exposure (37 s, 0.62 GB, a 2.0 MB fragment). The union + RECONCILES like `final_cat_merge` — a new exposure is OR-ed in on the spot, an + exposure that left the campaign or a fragment that changed forces a rebuild + (a union cannot be un-OR-ed), and a no-op leaves the file untouched — against + a sidecar `defect_map_.json` that records which exposures are + already in it. Memory is flat in the exposure count: fragments are read one at + a time and reduced to their pixel ids, so the job holds one accumulator (the + campaign's footprint, ~3 GB at DR6 scale) and one 2 MB fragment. - **External masks are wired, on the tile side only.** `inputs.masks` is a third input root beside tiles and exposures, exported as `$SP_INPUT_MASKS` and pointing at the UNIONS DR6 ugriz bit ladder: one boolean healsparse map per @@ -305,6 +337,9 @@ profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; kee and an unknown *name* is a parse-time error listing the valid ones. The list is exposure-side only; tile-side retention is #844 follow-up. - **The campaign ends in two merged catalogues, and the workflow makes both.** + (Three campaign products, counting the defect map above — but that one is a + map for the footprint, not a catalogue, and nothing downstream of it lives + here.) Everything above is per unit; the two products downstream analysis actually opens are per *campaign*, and until these rules existed each was a manual pass after the run. diff --git a/workflow/config/cfis/config_tile_Mc.ini b/workflow/config/cfis/config_tile_Mc.ini index a8316f8b7..05cb6f59f 100644 --- a/workflow/config/cfis/config_tile_Mc.ini +++ b/workflow/config/cfis/config_tile_Mc.ini @@ -117,4 +117,32 @@ SHAPE_MEASUREMENT_TYPE = ngmix # PSF-star selection deliberately does NOT consume these: MASK_PATHS in # config_exp_psfex.ini stays commented out, keeping the star diet narrow # (instrument flags only). See that file's header. +# +# THE CAMPAIGN'S OWN DEFECT MAP IS NOT IN THE LADDER BELOW, AND THAT IS +# DELIBERATE (CosmoStat/shapepipe#878). The workflow now emits one more boolean +# healsparse map in exactly this form -- nside 131072, coverage 128, True = +# masked -- from the exposures' instrument flag images: bad columns, saturated +# pixels and bleed trails, the one masking input that otherwise never leaves the +# pixel domain. But it is a campaign PRODUCT, not an input. Every path above +# exists before the run starts; that one exists only after it, so wiring it here +# would name a file the first run of a fresh campaign cannot have, and every +# tile would fail on a missing map. +# +# TO ADD IT, AFTER A CAMPAIGN HAS PRODUCED ONE. It lands at +# /defect_map_.hsp (`campaign:` in workflow/config.yaml, +# defaulting to the persistent root's basename). Copy or symlink it into the +# `inputs.masks` root -- the same root everything above resolves through, so the +# ladder keeps one location -- and append one entry: +# +# MASK_EXT_PATHS = ..., defect:$SP_INPUT_MASKS/defect_map_.hsp +# +# then add the matching `MASK_defect` line to final_cat.param, which is what +# turns the query into a column. The label is free; `defect` reads better than +# an `n` name, since this map is not one bit of the UNIONS ladder. +# +# TWO THINGS TO KNOW BEFORE YOU DO. It is the campaign's OWN footprint, so a +# catalogue built from a different tile list reads False -- not clean, unknown -- +# wherever that campaign had no exposure; and the rasterization is conservative, +# widening a one-pixel bad column to the 1.61" healpix resolution, which is the +# price of keeping thin defects at all. Nothing here cuts on it either way. MASK_EXT_PATHS = n1:$SP_INPUT_MASKS/mask_ugriz_nside131072_n1.hsp, n2:$SP_INPUT_MASKS/mask_ugriz_nside131072_n2.hsp, n4:$SP_INPUT_MASKS/mask_ugriz_nside131072_n4.hsp, n8:$SP_INPUT_MASKS/mask_ugriz_nside131072_n8.hsp, n16:$SP_INPUT_MASKS/mask_ugriz_nside131072_n16.hsp, n32:$SP_INPUT_MASKS/mask_ugriz_nside131072_n32.hsp, n64:$SP_INPUT_MASKS/mask_ugriz_nside131072_n64.hsp, n128:$SP_INPUT_MASKS/mask_ugriz_nside131072_n128.hsp, n256:$SP_INPUT_MASKS/mask_ugriz_nside131072_n256.hsp, n1024:$SP_INPUT_MASKS/mask_ugriz_nside131072_n1024.hsp, n2048:$SP_INPUT_MASKS/mask_ugriz_nside131072_n2048.hsp From 958b56230c9864e828d5ab5b82a5fcf4d60901c9 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 21:13:25 -0400 Subject: [PATCH 3/7] fix(defect-map): a short split is an error, a bad chip is not an OOM MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three defects in the two scripts, found reviewing them. A CCD WHOSE FLAG SPLIT IS MISSING was silently dropped. ccd_files globbed flag-*.fits and rasterized whatever came back, so an incompletely materialised split dir — an age-based scratch purge deleting files one at a time, a store copied or restored half way, a truncated rsync — produced a fragment covering half the exposure, written with "status": "complete", and nothing downstream able to notice. Reproduced on a fixture: remove one flag split of two and the script exited 0 with n_ccds 1. The guard it had was on the opposite, unreachable direction (a flag without its image), which split_exp cannot produce. The loop is now driven by the EXPECTED CCD count, which comes in on params from config_exp_Sp.ini's own N_HDU, and a short split dir is a hard error. A FULLY FLAGGED CCD blew the rule's memory budget. A MegaCam exposure routinely carries a dead or saturated chip: 2048 x 4612 = 9.4M flagged pixels, 85M samples at oversample 3, held in one shot several GB of coordinate arrays and astropy temporaries against a 2 GB request — so the job OOMed on all three attempts, the exposure never got a fragment, and because clean_exposure waits on that manifest it could never be reclaimed either. The 0.62 GB measured on a 4.1%-flagged exposure said nothing about that case. The flagged-pixel loop is now batched at CHUNK source pixels and np.unique-accumulated, so peak RSS is set by the batch and not by how bad the CCD is. Measured on exactly that case, a fully flagged chip against 2079612p's real WCS: 0.74 GB and 18.3 s at 500k, 1.13 GB and 18.9 s at 1M — time is flat in the batch size and memory is linear in it, so 500k is free. The real exposure is unchanged bit for bit (359858 healpix pixels, 2.0 MB) at 34 s. THE SIDECAR WENT STALE ON A NO-OP. main() returned on an empty plan without touching it, but two of its fields describe the CAMPAIGN and not the map: add tiles whose exposures all lack fragments and campaign_exposures and exposures_without_fragment are wrong on disk while the union is untouched and correct. That is exactly the record the docstring promises is "recorded rather than merely printed". The record is now built apart from writing the map and rewritten alone when it differs; write_stable keeps the map's mtime where it is. Also: the stamp is IMPORTED from hdf5_reconcile rather than re-derived, so the two halves of a campaign's reconciliation answer "did this source change" the same way; and the docstring now says plainly, in hdf5_reconcile's own terms, that the map's CONTENT is a function of the input set and its BYTES are not (measured: 1,586,880 B rebuilt vs 1,589,760 B appended, identical valid_pixels). tests/unit/test_defect_map_reconcile.py pins the four reconcile branches, the sidecar refresh and the nside mismatch — fragments made with make_empty, since the merge half never opens a FITS image or a WCS. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- tests/unit/test_defect_map_reconcile.py | 208 ++++++++++++++++++++++++ workflow/scripts/defect_map_exp.py | 99 ++++++++--- workflow/scripts/merge_defect_map.py | 111 +++++++++---- 3 files changed, 370 insertions(+), 48 deletions(-) create mode 100644 tests/unit/test_defect_map_reconcile.py diff --git a/tests/unit/test_defect_map_reconcile.py b/tests/unit/test_defect_map_reconcile.py new file mode 100644 index 000000000..ddc8fb133 --- /dev/null +++ b/tests/unit/test_defect_map_reconcile.py @@ -0,0 +1,208 @@ +"""``merge_defect_map`` agrees the campaign's defect map, or rebuilds it. + +The union of the per-exposure defect fragments (CosmoStat/shapepipe#878) is the +one campaign product whose reconciliation is ASYMMETRIC, and that asymmetry is +the reason this file exists. A union cannot be un-OR-ed: two exposures both set +a pixel and nothing in the map records which. So + + * a NEW fragment is OR-ed in on the spot — the cheap, common path; + * a fragment that LEFT the campaign, or one that CHANGED on disk, forces a + REBUILD from every fragment; + * neither leaves the map UNTOUCHED, mtime included, because mtime is a + rerun trigger and an unconditional rewrite makes every invocation look + like a change. + +Those are four branches of ``reconcile_plan`` whose failure mode is silent: an +edit that stopped treating a removal as a rebuild leaves the map carrying bits +from exposures the campaign no longer has, and nothing downstream — nothing in +this workflow reads the map at all — would ever notice. Hence the pins here. + +The SIDECAR is pinned alongside, for the fifth case the plan cannot see: the +record carries two fields about the CAMPAIGN (how many exposures it has, which +of them have no fragment) that can move while the map cannot. The docstring +promises a short map says so on disk; that only holds if a no-op still refreshes +the record. + +Fragments are made with ``HealSparseMap.make_empty`` and a handful of pixel ids +— the merge half never opens a FITS image or a WCS, so nothing here needs one. +Needs healsparse, so it runs inside the container and skips outside. +""" + +import importlib.util +import json +import sys +from pathlib import Path + +import pytest + + +pytestmark = pytest.mark.unions + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPTS = REPO_ROOT / "workflow" / "scripts" +SCRIPT = SCRIPTS / "merge_defect_map.py" + +NSIDE = 4096 +NSIDE_COV = 32 + + +def _load(): + """Import the rule's script by path — scripts/ is not a package.""" + assert SCRIPT.exists(), f"{SCRIPT} not found; the rule calls it by path" + sys.path.insert(0, str(SCRIPTS)) + try: + spec = importlib.util.spec_from_file_location("_merge_defect_map", + SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + finally: + sys.path.remove(str(SCRIPTS)) + return module + + +@pytest.fixture(scope="module") +def merge(): + pytest.importorskip("healsparse") + pytest.importorskip("h5py") # hdf5_reconcile, where the stamp lives + return _load() + + +def _fragment(merge, root: Path, exp: str, pixels) -> Path: + """Write one exposure's fragment where ``fragment_path`` expects it.""" + import numpy as np + import healsparse as hsp + + path = merge.fragment_path(root, exp) + path.parent.mkdir(parents=True, exist_ok=True) + frag = hsp.HealSparseMap.make_empty(NSIDE_COV, NSIDE, np.bool_, + bit_packed=True) + frag[np.asarray(pixels, dtype=np.int64)] = True + frag.write(str(path), clobber=True) + return path + + +def _valid(path: Path): + import healsparse as hsp + return set(int(p) for p in hsp.HealSparseMap.read(str(path)).valid_pixels) + + +def _run(merge, root: Path, have: dict, missing=()): + """plan + apply, the way ``main`` does, returning (plan, record).""" + output = root / "defect_map_test.hsp" + sidecar = root / "defect_map_test.json" + plan = merge.reconcile_plan(output, sidecar, have) + if plan.empty(): + return plan, json.loads(sidecar.read_text()) + record = merge.apply_plan(output, sidecar, plan, have, list(missing), + NSIDE_COV, NSIDE) + return plan, record + + +@pytest.fixture +def campaign(merge, tmp_path): + """Two exposures, disjoint pixels, merged once. The starting state.""" + have = { + "2079612p": _fragment(merge, tmp_path, "2079612p", [10, 11, 12]), + "2079613p": _fragment(merge, tmp_path, "2079613p", [20, 21]), + } + plan, _ = _run(merge, tmp_path, have) + assert plan.rebuild, "the first merge has no map to append to" + return tmp_path, have + + +def test_first_merge_is_the_union(merge, campaign): + root, _ = campaign + assert _valid(root / "defect_map_test.hsp") == {10, 11, 12, 20, 21} + + +def test_append_reads_only_the_new_fragment(merge, campaign): + """A grown campaign is an APPEND, not a rebuild — that is the cheap path.""" + root, have = campaign + have = dict(have) + have["2079614p"] = _fragment(merge, root, "2079614p", [30]) + plan, record = _run(merge, root, have) + assert plan.rebuild == [], plan.describe() + assert plan.append == ["2079614p"] + assert _valid(root / "defect_map_test.hsp") == {10, 11, 12, 20, 21, 30} + assert set(record["exposures"]) == set(have) + + +def test_removal_forces_a_rebuild_and_drops_the_pixels(merge, campaign): + """The case a union cannot do incrementally, and the reason for rebuild.""" + root, have = campaign + have = {k: v for k, v in have.items() if k != "2079613p"} + plan, _ = _run(merge, root, have) + assert plan.append == [] and plan.rebuild == ["2079612p"], plan.describe() + assert "left the campaign" in plan.reason + assert _valid(root / "defect_map_test.hsp") == {10, 11, 12} + + +def test_changed_fragment_forces_a_rebuild(merge, campaign): + """A restamped fragment is not trusted to be a superset of what went in.""" + root, have = campaign + _fragment(merge, root, "2079613p", [20, 21, 22]) + import os + os.utime(have["2079613p"], (0, 0)) + plan, _ = _run(merge, root, have) + assert plan.rebuild == ["2079612p", "2079613p"], plan.describe() + assert "changed on disk" in plan.reason + assert _valid(root / "defect_map_test.hsp") == {10, 11, 12, 20, 21, 22} + + +def test_no_op_leaves_the_map_untouched(merge, campaign): + """UNTOUCHED, not rewritten identically: mtime is a rerun trigger.""" + root, have = campaign + output = root / "defect_map_test.hsp" + before = output.stat().st_mtime_ns + plan, _ = _run(merge, root, have) + assert plan.empty(), plan.describe() + assert output.stat().st_mtime_ns == before + + +def test_no_op_still_refreshes_a_stale_sidecar(merge, campaign): + """The campaign moved, the map could not: the RECORD must still say so. + + Tiles whose exposures were all reclaimed by a workflow predating this rule + add nothing to merge and nothing to remove — an empty plan — but they change + what the campaign asked for. A sidecar that kept reporting the old counts + would make a short map look complete on disk. + """ + root, have = campaign + sidecar = root / "defect_map_test.json" + output = root / "defect_map_test.hsp" + before_map = output.stat().st_mtime_ns + plan, record = _run(merge, root, have) + assert plan.empty() + + missing = ["2079999p"] + plan = merge.reconcile_plan(output, sidecar, have) + assert plan.empty(), "an exposure with no fragment is not in the plan" + fresh = merge.build_record( + have, missing, NSIDE_COV, NSIDE, + record["n_pixels"], record["n_coverage_pixels"]) + merge.write_sidecar(sidecar, fresh) + + after = json.loads(sidecar.read_text()) + assert after["exposures_without_fragment"] == missing + assert after["campaign_exposures"] == len(have) + 1 + assert output.stat().st_mtime_ns == before_map, "map must not move" + + +def test_nside_mismatch_is_an_error_not_an_upgrade(merge, campaign): + """The ladder's resolution is a campaign decision, not a per-fragment one.""" + import numpy as np + import healsparse as hsp + + root, have = campaign + odd = merge.fragment_path(root, "2079615p") + odd.parent.mkdir(parents=True, exist_ok=True) + frag = hsp.HealSparseMap.make_empty(NSIDE_COV, NSIDE // 2, np.bool_, + bit_packed=True) + frag[np.asarray([5], dtype=np.int64)] = True + frag.write(str(odd), clobber=True) + + have = dict(have) + have["2079615p"] = odd + with pytest.raises(SystemExit) as exc: + _run(merge, root, have) + assert "re-rasterize 2079615p" in str(exc.value) diff --git a/workflow/scripts/defect_map_exp.py b/workflow/scripts/defect_map_exp.py index 17a9c767e..7473f42b3 100644 --- a/workflow/scripts/defect_map_exp.py +++ b/workflow/scripts/defect_map_exp.py @@ -77,6 +77,13 @@ sit inside it — the residual 1.2% is boundary and rounding, not a class of missed defect. +MEMORY IS BOUNDED BY THE BATCH, NOT BY THE EXPOSURE. Peak RSS is one CCD's +flag image plus one batch of ``CHUNK`` source pixels' worth of coordinates — +measured 0.62 GB on a real 4.1%-flagged exposure at oversample 3, and 0.74 GB on +a FULLY flagged CCD, which is the case the batching exists for (a dead or +saturated MegaCam chip is 9.4M flagged pixels and, unbatched, several GB; +``rasterize_ccd`` argues it). The bound is the batch, not the exposure. + BYTE-STABLE, tmp-then-``cmp``-then-``mv``, the pattern ``persist_exp`` uses: healsparse's FITS output carries no timestamp (checked), so re-rasterizing an unchanged store produces an identical file and leaves its mtime alone. mtime is @@ -127,16 +134,25 @@ def split_dir(exp_dir: Path) -> Path: return exp_dir / "output" / RUN_NAME / MODULE / "output" -def ccd_files(exp_dir: Path) -> list: +def ccd_files(exp_dir: Path, n_ccds: int) -> list: """``(ccd, flag path, image path)`` for every CCD this exposure split, in - CCD order. - - Driven by the FLAG files, which are what we rasterize; the image split is - looked up beside each one for its header alone. A flag split with no image - beside it is a hard error rather than a skip — ``exp_split`` writes the - three suffixes in one pass, so a missing image means the store is not what - we think it is, and a fragment quietly short of a CCD would be a hole in the - footprint nothing downstream could notice. + CCD order — all ``n_ccds`` of them or none at all. + + COUNTED AGAINST THE EXPECTED CCD COUNT, not against whatever is on disk, and + that is the whole guard. ``split_exp`` writes image, weight and flag for + each of ``N_HDU`` CCDs in one pass, so the reachable failure is not "a flag + without its image" — it is an incompletely MATERIALISED split dir: an + age-based scratch purge deleting files one at a time, a store copied or + restored half way, a truncated rsync. Globbing for ``flag-*.fits`` and + rasterizing whatever comes back turns that into a fragment covering half the + exposure, written with ``"status": "complete"`` and with nothing downstream + able to notice — a hole in the footprint, which is precisely what this rule + exists to prevent. So the expected count comes in on ``params`` (the + Snakefile reads ``N_HDU`` from ``config_exp_Sp.ini``) and a short split is a + hard error. + + The image split is looked up beside each flag for its header alone, and a + missing one is the same hard error for the same reason. """ root = split_dir(exp_dir) out = [] @@ -150,6 +166,12 @@ def ccd_files(exp_dir: Path) -> list: f"it in {root}; the WCS lives on the image split (see the " f"module docstring)") out.append((int(match.group(1)), flag, image)) + if len(out) != n_ccds: + sys.exit(f"defect_map_exp: {root} holds {len(out)} flag split(s), not " + f"the {n_ccds} this exposure was split into; the split dir is " + f"incomplete and a fragment built from it would be a hole in " + f"the footprint marked complete (see ccd_files' docstring). " + f"Re-run exp_split for this exposure.") return sorted(out) @@ -169,13 +191,36 @@ def offsets(oversample: int) -> tuple: return grid_x.ravel(), grid_y.ravel() +# Flagged CCD pixels converted per batch. Peak RSS is set by THIS, not by how +# bad the CCD is: one batch at oversample 3 is 500k x 9 samples x two float64 +# coordinate arrays in and two out, plus astropy's PV/SIP temporaries, and the +# accumulated result is a deduplicated int64 pixel list bounded by the CCD's +# healpix footprint (124601 ids for a WHOLE CCD at nside 131072), not by the +# sample count. Measured on a fully flagged MegaCam chip (2048 x 4612 = 9.4M +# pixels, the worst case there is) against 2079612p CCD 1's real WCS: 0.74 GB +# peak and 18.3 s at 500k, 1.13 GB and 18.9 s at 1M. Time is flat in the batch +# size and memory is linear in it, so the smaller batch is free. +CHUNK = 500_000 + + def rasterize_ccd(flag_path: Path, image_path: Path, nside: int, off_x, off_y) -> np.ndarray: """The healpix pixel ids (NEST, ``nside``) this CCD's flags touch. - One CCD at a time and one array at a time: the whole point of the loop in - ``main`` is that the job's footprint is one CCD's samples, not the - exposure's. + ONE CCD AT A TIME AND, WITHIN IT, ONE BATCH AT A TIME. The first is why the + loop in ``main`` is a loop; the second is why this one is. A MegaCam + exposure routinely carries a dead or saturated chip, and a FULLY flagged CCD + is 2048 x 4612 = 9.4M nonzero pixels — 85M samples at oversample 3. Held in + one shot that is ~680 MB per coordinate array in and the same again out of + ``all_pix2world``, plus astropy's own PV/SIP temporaries: several GB, well + over the rule's request, on all three attempts. Batched it is 0.74 GB + (measured), inside the request with room to spare. The exposure would then + never get a fragment AND, because ``clean_exposure`` waits on this rule's + manifest, never be reclaimed either. The 0.62 GB measured on a 4.4%-flagged + exposure says nothing about that case; ``CHUNK`` does. + + The batches are ``np.unique``-reduced as they go, so what survives across + them is the CCD's healpix footprint and not its samples. """ with warnings.catch_warnings(): # SCAMP headers carry a deprecated RADECSYS and a redundant SIP block @@ -185,15 +230,25 @@ def rasterize_ccd(flag_path: Path, image_path: Path, nside: int, wcs = WCS(fits.getheader(image_path)) data = fits.getdata(flag_path) rows, cols = np.nonzero(data) + del data if rows.size == 0: return np.empty(0, dtype=np.int64) - # 1-based FITS pixel coordinates, sampled across each flagged pixel's extent. - x = (cols[:, None] + 1.0 + off_x[None, :]).ravel() - y = (rows[:, None] + 1.0 + off_y[None, :]).ravel() - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - ra, dec = wcs.all_pix2world(x, y, 1) - return np.unique(hp.ang2pix(nside, ra, dec, lonlat=True, nest=True)) + found = np.empty(0, dtype=np.int64) + for start in range(0, rows.size, CHUNK): + stop = start + CHUNK + # 1-based FITS pixel coordinates, sampled across each flagged pixel's + # extent. + x = (cols[start:stop, None] + 1.0 + off_x[None, :]).ravel() + y = (rows[start:stop, None] + 1.0 + off_y[None, :]).ravel() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + ra, dec = wcs.all_pix2world(x, y, 1) + del x, y + pixels = hp.ang2pix(nside, ra, dec, lonlat=True, nest=True) + del ra, dec + found = np.union1d(found, pixels) + del pixels + return found def write_stable(tmp: Path, dest: Path) -> None: @@ -213,6 +268,10 @@ def main() -> None: help="/exp///defect; the " "fragment is /defect-.hsp") parser.add_argument("--manifest", required=True, type=Path) + parser.add_argument("--n-ccds", required=True, type=int, + help="how many CCDs exp_split wrote (N_HDU in " + "config_exp_Sp.ini); a split dir short of this is " + "an error, not a smaller fragment") parser.add_argument("--nside", type=int, default=131072, help="nside_sparse; the mask ladder's resolution") parser.add_argument("--nside-coverage", type=int, default=128) @@ -221,7 +280,7 @@ def main() -> None: "docstring's measured table)") args = parser.parse_args() - ccds = ccd_files(args.exp_dir) + ccds = ccd_files(args.exp_dir, args.n_ccds) if not ccds: sys.exit(f"defect_map_exp: {args.exp}: no flag split under " f"{split_dir(args.exp_dir)}") diff --git a/workflow/scripts/merge_defect_map.py b/workflow/scripts/merge_defect_map.py index dab1dd219..2d7e13c58 100644 --- a/workflow/scripts/merge_defect_map.py +++ b/workflow/scripts/merge_defect_map.py @@ -35,6 +35,20 @@ * a plan with neither leaves the file UNTOUCHED — not rewritten identically, untouched, so its mtime cannot move. mtime is a rerun trigger. +WHAT IS AND IS NOT A FUNCTION OF THE INPUT SET, in the same terms +``hdf5_reconcile.py`` sets them. The map's CONTENT is: the same fragments give +the same valid pixels, the same counts and the same sidecar, whether they +arrived at once or one append at a time. Its BYTES are not — reaching a state by +append rather than by rebuild round-trips the map through healsparse's reader +and writer, which can lay the same pixels out in a different number of 2880-byte +FITS blocks (measured: 1,586,880 B rebuilt vs 1,589,760 B appended, identical +``valid_pixels``). That is the trade for not re-reading the campaign. It means +``write_stable`` below can move the map's mtime on a rebuild that changed +nothing — cheap while nothing consumes the map, and the thing to fix (by +rebuilding whenever the append path would rewrite anyway) if something ever +does. The no-op case is unaffected: it compares the PLAN, not the bytes, and +never opens the map at all. + Which exposures are already in the map is recorded in a SIDECAR beside it (``defect_map_.json``), each with its fragment's size and mtime — the same "did this change since we read it" stamp ``merge_final_cat.py`` records on @@ -85,6 +99,14 @@ # Same directory; the rule invokes this file by path, so it is sys.path[0]. import build_index +# The two hdf5 merges' reconciler. This file cannot use its `plan`/`apply` — a +# union is not a group of independent datasets, so removing an exposure is a +# rebuild here and a `del` there — but "did this source change since we read +# it" is the SAME question, and answering it twice in two ways is how the two +# halves of a campaign's reconciliation drift apart. So the stamp is imported, +# not re-derived, and the vocabulary below (`stamp`, `Plan`, `empty`, +# `describe`, plan-then-apply, untouched-on-a-no-op) is deliberately theirs. +from hdf5_reconcile import stamp def fragment_path(products_dir: Path, exp: str) -> Path: @@ -93,15 +115,9 @@ def fragment_path(products_dir: Path, exp: str) -> Path: / f"defect-{exp}.hsp") -def stamp(path: Path) -> list: - """The fragment's identity, as recorded in the sidecar. - - Size and mtime, not a checksum: the fragment is ~2 MB, it is written - byte-stably (so a no-op re-rasterization does not move its mtime), and the - question is only "did this change since we read it". - """ - st = path.stat() - return [st.st_size, st.st_mtime_ns] +def sidecar_stamp(path: Path) -> list: + """``stamp`` as JSON round-trips it: a list, so a read record compares.""" + return list(stamp(path)) def fragments(products_dir: Path, tile_list: Path, index_db: Path) -> tuple: @@ -159,7 +175,7 @@ def reconcile_plan(output: Path, sidecar: Path, have: dict) -> Plan: return Plan([], sorted(have), f"{len(gone)} exposure(s) left the campaign") changed = sorted(exp for exp, path in have.items() - if exp in known and list(known[exp]) != stamp(path)) + if exp in known and list(known[exp]) != sidecar_stamp(path)) if changed: return Plan([], sorted(have), f"{len(changed)} fragment(s) changed on disk") @@ -194,6 +210,44 @@ def write_stable(tmp: Path, dest: Path) -> None: tmp.replace(dest) +def build_record(have: dict, missing: list, nside_coverage: int, nside: int, + n_pixels: int, n_coverage: int) -> dict: + """The sidecar: what is in the map, and what the campaign wanted but lacked. + + Built apart from writing the map because it is not only the map's record. + Two of its fields describe the CAMPAIGN — how many exposures it has and + which of them have no fragment — and those can move while the map itself + cannot: add tiles whose exposures were all reclaimed by a workflow + predating this rule and there is nothing to append, nothing to rebuild, and + a sidecar still reporting the previous campaign's counts. The docstring says + a short map should say so ON DISK; that means the record has to be rewritten + even when the map is untouched. + """ + return { + "campaign_exposures": len(have) + len(missing), + "nside": nside, + "nside_coverage": nside_coverage, + "n_pixels": n_pixels, + "n_coverage_pixels": n_coverage, + # What the NEXT invocation reconciles against; sorted so the sidecar is + # byte-stable for a given campaign state. + "exposures": {exp: sidecar_stamp(path) + for exp, path in sorted(have.items())}, + # Recorded rather than merely printed: a map short of exposures should + # say so on disk, not only in a job log nobody keeps. + "exposures_without_fragment": sorted(missing), + } + + +def write_sidecar(sidecar: Path, record: dict) -> None: + tmp = sidecar.with_name(sidecar.name + ".tmp") + try: + tmp.write_text(json.dumps(record, indent=2, sort_keys=True) + "\n") + write_stable(tmp, sidecar) + finally: + tmp.unlink(missing_ok=True) + + def apply_plan(output: Path, sidecar: Path, plan: Plan, have: dict, missing: list, nside_coverage: int, nside: int) -> dict: """Carry the plan out on tmp copies, then move both files into place. @@ -211,30 +265,17 @@ def apply_plan(output: Path, sidecar: Path, plan: Plan, have: dict, accumulate(target, [(exp, have[exp]) for exp in todo], nside_coverage, nside) - record = { - "campaign_exposures": len(have) + len(missing), - "nside": nside, - "nside_coverage": nside_coverage, - "n_pixels": int(target.n_valid), - "n_coverage_pixels": int(target.coverage_mask.sum()), - # What the NEXT invocation reconciles against; sorted so the sidecar is - # byte-stable for a given campaign state. - "exposures": {exp: stamp(path) for exp, path in sorted(have.items())}, - # Recorded rather than merely printed: a map short of exposures should - # say so on disk, not only in a job log nobody keeps. - "exposures_without_fragment": missing, - } + record = build_record(have, missing, nside_coverage, nside, + int(target.n_valid), + int(target.coverage_mask.sum())) map_tmp = output.with_name(output.name + ".tmp") - side_tmp = sidecar.with_name(sidecar.name + ".tmp") try: target.write(str(map_tmp), clobber=True) - side_tmp.write_text(json.dumps(record, indent=2, sort_keys=True) + "\n") write_stable(map_tmp, output) - write_stable(side_tmp, sidecar) finally: map_tmp.unlink(missing_ok=True) - side_tmp.unlink(missing_ok=True) + write_sidecar(sidecar, record) return record @@ -263,8 +304,22 @@ def main() -> None: args.output.parent.mkdir(parents=True, exist_ok=True) plan = reconcile_plan(args.output, args.sidecar, have) if plan.empty(): + # The MAP is untouched — that is what an empty plan means, and its mtime + # must not move. The RECORD still can be stale: campaign_exposures and + # exposures_without_fragment describe the campaign, not the map, so + # tiles whose exposures all lack fragments change them without changing + # a single bit of the union. Rewrite it alone when it differs; + # write_stable drops the tmp when it does not. + old_record = read_sidecar(args.sidecar) + record = build_record(have, missing, args.nside_coverage, args.nside, + int(old_record.get("n_pixels", 0)), + int(old_record.get("n_coverage_pixels", 0))) + stale = record != old_record + if stale: + write_sidecar(args.sidecar, record) print(f"[merge_defect_map] unchanged: {args.output} " - f"({len(have)} exposure(s))") + f"({len(have)} exposure(s)" + f"{'; sidecar refreshed' if stale else ''})") return record = apply_plan(args.output, args.sidecar, plan, have, missing, args.nside_coverage, args.nside) From a9daa4f3bcff75863e3a38ad2ca776c573b6ee89 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 21:14:30 -0400 Subject: [PATCH 4/7] fix(defect-map): a reclaimed exposure is a leaf, and the merge asks for 7 GB MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Four defects in the DAG wiring and the sizing. A RERUN TRIGGER RE-AVALANCHED THE CAMPAIGN. defect_map_inputs() named the manifest of any exposure that had one, reclaimed or not, but exp_defect_map's input is exp_split's manifest and that went with the scratch store. So any trigger on the rule — and params is live, carrying the script hash, the nside and the oversampling this file advertises as a minute's work to change — would schedule exp_get_images and exp_split from VOS for every reclaimed exposure: a four-hour chain each, campaign-wide, on a one-character edit. It now mirrors star_cat_inputs(): a live exposure through its manifest, a reclaimed one through its FRAGMENT, which is not a declared output of any rule and is therefore a true leaf. That also gives prod_exp_fragment() its only caller, so the fragment's location is spelled once here and once in merge_defect_map.fragment_path() rather than twice with nothing joining them. Verified on a two-exposure compute fixture: with 2079612 reclaimed and holding a fragment, requesting the map builds 4 jobs (one exp chain, the live exposure's) where both-live builds 7. THE FIRST MERGE OF A LARGE CAMPAIGN ASKED FOR ~52 GB. Before there is a sidecar the coverage-pixel count was estimated as 13 per exposure capped at the FULL SKY, and 12 * 128^2 = 196608 is reached at ~15k exposures — so a DR6-scale campaign estimated the whole sky, 25.8 GB of accumulator and ~52 GB of request (~104 GB on the retry) for a job this file's own comment measures at ~3 GB. Exposures overlap almost completely; the prior has to be a FOOTPRINT, so it is now capped at the ~23k coverage pixels the UNIONS ugriz maps measure over the same sky. A small campaign is still sized on its own exposures, where the per-exposure figure is the honest one: the two-exposure fixture asks for 806 MB. THE MERGE'S RUNTIME REACHED 11.4 H WITH NO RETRIES AND NO PARTIAL STATE, so its mem_mb/runtime attempt scaling was dead code and one walltime overrun eleven hours in lost the whole thing. retries: 1 is declared, the formula is corrected to the ~1 s per exposure actually measured, and it is capped below the walltime Alliance policy lets a job run without checkpointing — with the comment saying that a campaign needing longer needs a resumable accumulator, not a bigger number. SILENCE WAS THE ANSWER when no exposure can contribute. A campaign whose exposures were all reclaimed by a workflow predating this rule builds no defect job at all, and the operator saw nothing rather than "no map is possible here". It now says so once at parse time, like the missing-index note. Also: a comment pointed at DEFECT_COV_BYTES, which does not exist. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- workflow/Snakefile | 119 +++++++++++++++++++++++++++++++----- workflow/rules/exposure.smk | 67 ++++++++++++++------ 2 files changed, 154 insertions(+), 32 deletions(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index 0cde94c5a..296b38fe7 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -31,6 +31,7 @@ the manifest says "this stage succeeded", the log says "here is what happened" (the contract is argued in completeness.py's docstring). """ +import configparser import fnmatch import functools import hashlib @@ -558,6 +559,28 @@ if DEFECT_NSIDE_COVERAGE >= DEFECT_NSIDE: f"below nside={DEFECT_NSIDE}.") +def _split_n_hdu(): + """How many CCDs exp_split writes, from ITS OWN config. + + Read rather than repeated, because it is the number defect_map_exp.py + checks the split dir against: a constant here that drifted from + config_exp_Sp.ini would either fail every exposure or restore the silent + half-exposure fragment the check exists to stop. It rides on the rule's + params, so changing N_HDU re-rasterizes. + """ + parser = configparser.ConfigParser() + parser.read(CONFIG_DIR / "config_exp_Sp.ini") + try: + return int(parser["SPLIT_EXP_RUNNER"]["N_HDU"]) + except (KeyError, ValueError) as exc: + raise WorkflowError( + f"config_exp_Sp.ini: no readable N_HDU ({exc}); exp_defect_map " + f"needs it to tell a complete split dir from a truncated one.") + + +DEFECT_N_CCDS = _split_n_hdu() + + def defect_map(): """The campaign's merged defect map. A PRODUCT, not an input: nothing in this workflow opens it, and config_tile_Mc.ini deliberately does not name @@ -601,23 +624,58 @@ def defect_map_inputs(): """What defect_map_merge waits for: the fragment manifests of the exposures whose stores this invocation can still reach. - An exposure whose scratch store is gone (reclaimed, or purged) and which has - no fragment manifest yet cannot get one without rebuilding its chain from - VOS — the avalanche persist_targets() drops such exposures to avoid. It is - left out of the edge; the job skips it and records it as missing. Every - other exposure's manifest is either already there or one exp_split away. + THE SAME SPLIT AS star_cat_inputs(), AND FOR THE SAME REASON. A live + exposure is depended on through its MANIFEST, which is what orders the merge + after the rasterization. A reclaimed one is depended on through its + FRAGMENT — the .hsp, which is not a declared output of any rule, so a + fragment that exists is a DAG leaf snakemake requires and builds nothing + for. + + Naming the manifest for a reclaimed exposure would be the avalanche + persist_targets() drops such exposures to avoid: this rule's input is + exp_split's manifest, which went with the scratch store, so ANY rerun + trigger on exp_defect_map — and `params` is a live trigger, carrying the + script hash, the nside and the oversampling this file advertises as cheap to + change — would schedule exp_get_images and exp_split from VOS, a four-hour + chain per exposure, campaign-wide, on a one-character edit. The fragment is + already on the persistent root and cannot be rebuilt in place; requiring it + asks for nothing. + + An exposure reclaimed by a workflow PREDATING this rule has neither, and is + left out of the edge entirely; the job skips it and records it as missing. """ - return [prod_exp_manifest(e, "exp_defect_map") - for e in defect_map_exposures() - if not exp_store_reclaimed(e) - or Path(prod_exp_manifest(e, "exp_defect_map")).exists()] + live, reclaimed = [], [] + for exp in defect_map_exposures(): + if not exp_store_reclaimed(exp): + live.append(prod_exp_manifest(exp, "exp_defect_map")) + elif Path(prod_exp_fragment(exp)).exists(): + reclaimed.append(prod_exp_fragment(exp)) + return live + reclaimed def defect_map_targets(): - """The merged map, whenever this campaign has an exposure to rasterize.""" + """The merged map, whenever this campaign has an exposure to rasterize. + + SILENCE IS THE WRONG ANSWER when there is none. Every exposure of a campaign + reclaimed by a workflow predating this rule is out of defect_map_inputs(), + so the DAG carries no defect job at all and an operator sees nothing — not + "no map is possible here". Say it once, at parse time, like the missing + index note above. + """ if not workflow.is_main_process: return [] - return [defect_map(), defect_map_sidecar()] if defect_map_inputs() else [] + if not defect_map_inputs(): + if defect_map_exposures(): + logger.warning( + f"defect map: none of this campaign's " + f"{len(defect_map_exposures())} exposure(s) can contribute — " + f"each is reclaimed and has no fragment on the persistent root, " + f"so no fragment can be rasterized without rebuilding its chain " + f"from VOS. No defect map will be produced for this campaign " + f"(workflow/scripts/merge_defect_map.py argues why that is not " + f"an error).") + return [] + return [defect_map(), defect_map_sidecar()] def defect_map_fragment_targets(): @@ -628,6 +686,11 @@ def defect_map_fragment_targets(): the store at the purge. Reaching them only through defect_map_merge would make a campaign that never merged lose them. + It is defect_map_inputs() verbatim, so it inherits that function's split: + a live exposure's manifest, which is the thing to BUILD, and a reclaimed + one's existing fragment, which is already where it needs to be and asks for + nothing. + HEAD PROCESS ONLY, for the same cost reason as persist_targets(). """ if not workflow.is_main_process: @@ -945,10 +1008,16 @@ TILE_BYTES_DEFAULT = 46_000_000 # (~2 MB) at a time. Two ways to know the count, in order: the sidecar the last # merge wrote, which is the campaign's own measured footprint; and before there # is one, the exposures times the 13 coverage pixels ONE exposure touched -# (measured, 2079612p), which ignores overlap and is therefore an overestimate -# in the safe direction, capped at the whole sky. +# (measured, 2079612p), capped at the survey FOOTPRINT — exposures overlap +# almost completely, so the per-exposure sum is only honest while the campaign +# is small, and unbounded it would ask for the whole sky +# (defect_map_cov_bytes() carries the arithmetic). DEFECT_MEM_BASE_MB = 800 # interpreter + astropy + healpy + healsparse DEFECT_COV_PER_EXP = 13 +# The merge holds no partial state, so it must stay inside the walltime that +# Alliance policy lets a job run without checkpointing. A campaign that needs +# longer needs a resumable accumulator, not a bigger number here. +DEFECT_RUNTIME_CAP_MIN = 660 # 11 h def _size(path, default): @@ -959,16 +1028,36 @@ def _size(path, default): return default +# The measured DR6 FOOTPRINT in nside-128 coverage pixels: ~23k, taken from the +# UNIONS ugriz maps staged under `inputs.masks`, which cover the same sky the +# exposures do. It is the ceiling on the pre-sidecar estimate below, and it is a +# FOOTPRINT rather than a whole sky: the survey is ~5000 deg^2, an eighth of the +# 12 * 128^2 = 196608 coverage pixels a full-sky cap would allow. +DEFECT_COV_FOOTPRINT = 23_000 + + def defect_map_cov_bytes(): """Bytes the union's accumulator occupies: coverage pixels x the bit-packed - block size. See the sizing block above for where the count comes from.""" + block size. + + ONCE THERE IS A SIDECAR the count is the map's own, measured. Before that it + is a prior, and the prior must be a FOOTPRINT estimate, not a sum over + exposures. Exposures overlap almost completely — each tile is covered by 7-10 + of them and they tile the same sky — so DEFECT_COV_PER_EXP * n_exposures + passes the survey footprint at ~1800 exposures and, uncapped at the full sky, + asks for 25.8 GB of accumulator (and ~52 GB of rule, doubled again on a + retry) for a job the merge measures at ~3 GB. Capping at the measured DR6 + footprint keeps the first-ever merge of a large campaign asking for what it + needs; a small campaign is still sized on its own exposures, where the + per-exposure figure is the honest one. + """ block = (DEFECT_NSIDE // DEFECT_NSIDE_COVERAGE) ** 2 // 8 try: record = json.loads(Path(defect_map_sidecar()).read_text()) n_cov = int(record["n_coverage_pixels"]) except (OSError, ValueError, KeyError): n_cov = min(DEFECT_COV_PER_EXP * len(defect_map_exposures()), - 12 * DEFECT_NSIDE_COVERAGE ** 2) + DEFECT_COV_FOOTPRINT) return n_cov * block diff --git a/workflow/rules/exposure.smk b/workflow/rules/exposure.smk index 33fd71763..f37c3f056 100644 --- a/workflow/rules/exposure.smk +++ b/workflow/rules/exposure.smk @@ -198,8 +198,8 @@ rule exp_persist: # # NOT A LOCALRULE, and this is where it parts company with exp_persist. That # rule is a tar of a few MB — seconds, far shorter than the scheduling latency. -# This one is 40 CCDs of WCS transforms and ang2pix over ~17 M flagged pixels: -# 37 s measured end to end on a real exposure at oversample 3 (2079612p, on this +# This one is 40 CCDs of WCS transforms and ang2pix over ~16 M flagged pixels: +# 34 s measured end to end on a real exposure at oversample 3 (2079612p, on this # login node, inside the campaign container; 21 s at oversample 2). That is real # work, it is CPU-bound, and running ~20k of them under local-cores in the head # process would serialise the campaign behind them. @@ -213,11 +213,17 @@ rule exp_defect_map: output: manifest = f"{PROD_EXP_DIR}/manifests/exp_defect_map.json" # No `log:`, for exp_persist's reason: the failure modes are "no flag split - # under the store" and "a flag split with no image beside it", both reported - # on stderr, neither with a per-CCD verdict worth a completeness record. + # under the store", "fewer flag splits than N_HDU" and "a flag split with no + # image beside it", all reported on stderr, none with a per-CCD verdict + # worth a completeness record. params: exp_dir = lambda wc: exp_dir(wc.exp), dest = lambda wc: f"{prod_exp_dir(wc.exp)}/defect", + # config_exp_Sp.ini's own N_HDU, read at parse time. A split dir short + # of it is a hard error, not a smaller fragment: half an exposure's + # defects, written "complete", is a hole in the footprint nothing + # downstream can see (defect_map_exp.py's ccd_files argues it). + n_ccds = DEFECT_N_CCDS, nside = DEFECT_NSIDE, nside_cov = DEFECT_NSIDE_COVERAGE, oversample = DEFECT_OVERSAMPLE, @@ -226,21 +232,31 @@ rule exp_defect_map: retries: 2 resources: # Measured on 2079612p inside the container: peak RSS 0.62 GB at - # oversample 3 (0.41 GB at 2), dominated by one CCD's sample arrays - # (430k flagged pixels x 9 samples x two float64 coordinate arrays) - # plus the bit-packed fragment's 13 coverage pixels. Flat in the number - # of CCDs — they are rasterized one at a time — so this is sized on the - # worst CCD, not on the exposure, at ~3x it. It scales with - # `oversample`, which is why that knob is not free. + # oversample 3 (0.41 GB at 2), dominated by one BATCH of a CCD's sample + # arrays plus the bit-packed fragment's 13 coverage pixels. + # + # FLAT IN BOTH DIRECTIONS THAT COULD BLOW IT: in the number of CCDs, + # which are rasterized one at a time, and in how badly any one of them + # is flagged, which is batched at defect_map_exp.CHUNK source pixels. + # The second is the one worth requesting for — a MegaCam exposure + # routinely carries a dead or saturated chip, 9.4M flagged pixels, and + # unbatched that is several GB and three OOMs, after which the exposure + # has no fragment AND cannot be reclaimed (clean_exposure waits on this + # manifest). Measured on exactly that case, a fully flagged chip against + # a real WCS: 0.74 GB. So 2000 covers the worst CCD at ~2.7x, not the + # measured average at ~3x. It scales with `oversample`, which is why + # that knob is not free. mem_mb = lambda wc, attempt: 2000 * attempt, - # 37 s measured; a factor of ~30 for a dirtier exposure, a higher - # oversampling and a busy filesystem. + # 34 s measured end to end on 2079612p at oversample 3; a factor of + # ~35 for a dirtier exposure (a fully flagged chip is 18 s on its own), + # a higher oversampling and a busy filesystem. runtime = 20 shell: "set -euo pipefail\n" f"python {SCRIPTS}/defect_map_exp.py" " --exp-dir '{params.exp_dir}' --exp {wildcards.exp}" " --dest '{params.dest}' --manifest {output.manifest}" + " --n-ccds {params.n_ccds}" " --nside {params.nside} --nside-coverage {params.nside_cov}" " --oversample {params.oversample}" @@ -427,19 +443,36 @@ rule defect_map_merge: inputs = unit_fingerprint(defect_map_exposures()), script_hash = MERGE_DEFECT_HASH threads: 1 + # Declared so the attempt scaling above is not dead code. One retry, not the + # two the exposure rules take: a failed attempt here has already cost hours, + # both declared outputs go with it, and the retry starts from an empty + # accumulator — there is nothing to salvage and little to gain from a third. + retries: 1 resources: # Sized on the FOOTPRINT, not on the exposure count: the accumulator is # one bit per sparse pixel of every touched coverage pixel, and the loop # holds one fragment at a time (merge_defect_map.py's memory argument). - # DEFECT_COV_BYTES is that arithmetic, from the sidecar's own recorded - # coverage count once there is one and from a measured per-exposure - # figure before that; the Snakefile's sizing block carries both. + # defect_map_cov_bytes() is that arithmetic, from the sidecar's own + # recorded coverage count once there is one and, before that, from a + # per-exposure figure capped at the MEASURED DR6 footprint; the + # Snakefile's sizing block carries both and argues the cap. mem_mb = lambda wc, attempt: attempt * ( DEFECT_MEM_BASE_MB + 2 * defect_map_cov_bytes() // 1_000_000), # Dominated by reading fragments (~2 MB each) and setting their pixels; # ~1 s per exposure measured, over a floor that covers writing the map. - runtime = lambda wc, attempt: attempt * ( - 20 + len(defect_map_exposures()) // 30) + # + # AND IT IS THE EXPENSIVE CASE THAT SETS IT. An append reads exactly the + # appended exposures and finishes in minutes; a REBUILD — any exposure + # leaving the campaign, any fragment restamped — reads all of them, ~40 + # GB at DR6 scale, and that is the ~6 h this formula sizes for at 20k + # exposures. Capped below 12 h because the job holds no partial state + # and Alliance policy asks anything longer to checkpoint: it cannot, so + # it must not ask. If a campaign ever needs longer than DEFECT_RUNTIME_ + # CAP_MIN, the accumulator has to become resumable (write map and + # sidecar every N fragments) rather than the cap being raised. + runtime = lambda wc, attempt: min( + attempt * (20 + len(defect_map_exposures()) // 60), + DEFECT_RUNTIME_CAP_MIN) shell: "set -euo pipefail\n" f"python {SCRIPTS}/merge_defect_map.py" From 1f6e8327fed66ec7559887a1296d0ae153a861db Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 21:14:30 -0400 Subject: [PATCH 5/7] docs(defect-map): the measurements, and three per exposure rather than one MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The README and config.yaml carried the pre-batching numbers and an inode estimate a third of the real cost. The DR6 inode figure was ~20k against a ~1M group quota; it is ~60k, because each exposure gets a defect/ DIRECTORY, the fragment in it and a manifest under manifests/ — three inodes, not one. The fragment was also described as sitting "beside its PSF tar", which is not the layout; the path is spelled out. The timing is the re-measured 34 s (2079612p, 40 CCDs, 4.1% of pixels flagged, oversample 3, inside the campaign container), and both files now carry the worst-case memory bound the batching exists for — a fully flagged chip at 0.74 GB — rather than only the average exposure's 0.62 GB. The README also records what the reconcile does and does not promise (content, not bytes) and points at the new unit test. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- workflow/README.md | 12 ++++++++++-- workflow/config.yaml | 20 ++++++++++++++------ 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/workflow/README.md b/workflow/README.md index c9b176cd4..2276f963c 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -245,14 +245,22 @@ profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; kee input, so reclamation cannot overtake the copy. Its resolution and its oversampling ride on `params`; `config.yaml`'s `defect_map:` block carries both, the measured convergence table behind the default, and the measurement - on one real exposure (37 s, 0.62 GB, a 2.0 MB fragment). The union + on one real exposure (34 s, 0.62 GB, a 2.0 MB fragment; the rasterization is + batched, so a fully flagged chip — the worst case, and one a real exposure + carries whenever a chip is dead — is 0.74 GB rather than several). The union RECONCILES like `final_cat_merge` — a new exposure is OR-ed in on the spot, an exposure that left the campaign or a fragment that changed forces a rebuild (a union cannot be un-OR-ed), and a no-op leaves the file untouched — against a sidecar `defect_map_.json` that records which exposures are already in it. Memory is flat in the exposure count: fragments are read one at a time and reduced to their pixel ids, so the job holds one accumulator (the - campaign's footprint, ~3 GB at DR6 scale) and one 2 MB fragment. + campaign's footprint, ~3 GB at DR6 scale) and one 2 MB fragment. What the map + and the sidecar say is a function of the input set; the map's BYTES are not, + because reaching a state by append rather than by rebuild round-trips it + through healsparse's writer (`merge_defect_map.py` measures the difference and + says what would have to change if anything ever consumed the map). + `tests/unit/test_defect_map_reconcile.py` pins the four reconcile branches and + the sidecar refresh. - **External masks are wired, on the tile side only.** `inputs.masks` is a third input root beside tiles and exposures, exported as `$SP_INPUT_MASKS` and pointing at the UNIONS DR6 ugriz bit ladder: one boolean healsparse map per diff --git a/workflow/config.yaml b/workflow/config.yaml index 86581604a..f8a0687dd 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -187,8 +187,11 @@ persist_exp: # it. The survey footprint is built from the CCD corner WCS in the headers, so # it cannot subtract them and would silently include defective pixels. So # `exp_defect_map` rasterizes each exposure's flags into a boolean healsparse -# fragment beside its PSF tar on products_dir, and `defect_map_merge` unions the -# campaign's fragments into /defect_map_.hsp. +# fragment on products_dir, at +# /exp///defect/defect-.hsp -- its own subdir +# beside the exposure's manifests/ and its PSF tar, not beside the tar itself -- +# and `defect_map_merge` unions the campaign's fragments into +# /defect_map_.hsp. # # IT IS A PRODUCT, NOT AN INPUT. Nothing in the workflow reads it back; it is # for the footprint, downstream. In particular config_tile_Mc.ini's @@ -196,13 +199,18 @@ persist_exp: # the run starts, and this one exists only after it. That file's header carries # the recipe for adding it once a campaign has produced one. # -# Measured on one real exposure (2079612p, 40 CCDs, 4.4% of pixels flagged) on -# this login node inside the campaign container: 37 s at oversample 3, peak RSS +# Measured on one real exposure (2079612p, 40 CCDs, 4.1% of pixels flagged) on +# this login node inside the campaign container: 34 s at oversample 3, peak RSS # 0.62 GB, 359858 healpix pixels, a 2.0 MB fragment over 13 coverage pixels # (21 s / 0.41 GB / 357896 pixels at oversample 2; 95 s / 1.29 GB / 360840 at # oversample 5 -- so 3 is within 0.3% of 5 over the whole exposure, at a third -# of the cost). ~40 GB and ~20k inodes at DR6 scale, against a ~1 M-inode group -# quota. +# of the cost). The rasterization is BATCHED, so peak memory is the batch and +# not the exposure: a fully flagged MegaCam chip -- 9.4M pixels, the worst case +# there is, and one a real exposure carries whenever a chip is dead -- measures +# 0.74 GB and 18 s on its own. At DR6 scale (~20k exposures) the products are +# ~40 GB and ~60k inodes, +# against a ~1 M-inode group quota: THREE per exposure, not one -- the defect/ +# directory, the fragment in it, and the manifest under manifests/. # # EVERY KEY BELOW DEFAULTS TO THE MASK LADDER'S OWN CONVENTION and the block can # be omitted entirely. nside/nside_coverage match every map under `inputs.masks` From db4f426a5e476d8286f886e076a33cf1378d8122 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 21:46:32 -0400 Subject: [PATCH 6/7] fix(defect-map): reclamation never waits on a rasterization it cannot run MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit clean_exposure's exp_defect_map edge was unconditional, which reopened the avalanche defect_map_inputs() is written to avoid. An exposure whose scratch store went to the 60-day /scratch purge, or to a `clean: false` run, has NO tombstone — exp_store_reclaimed()'s docstring names that case — so clean_targets() still asks for its tombstone, the missing exp_defect_map manifest is scheduled, and its own input (exp_split's manifest) went with the store: exp_get_images and exp_split from VOS, four hours per exposure, campaign-wide, on the first run of this branch. Reproduced on a fixture campaign (tile 000.000, exposures 2079612 live and 2079613 purged-without-tombstone, tile_vignets present): a dry run of 2079613's tombstone scheduled 13 jobs including exp_get_images and exp_split; with the edge made conditional it schedules one, clean_exposure. The live exposure still gets exp_defect_map ordered before its clean, and a reclaimed exposure that already has a fragment depends on the fragment, a leaf that builds nothing. ccd_files' short-split error now says what to do about the store it pins: a damaged split dir is still a hard error rather than a half-footprint fragment marked complete, but the operator is told that deleting the scratch exp_split manifest is what rebuilds the split and unblocks the clean. defect_map_merge's mem_mb also goes through capped_mem(), as star_cat_merge and final_cat_merge already do: defect_map_cov_bytes() scales as nside^2 and defect_map.nside is an advertised knob, so one ladder change turns the request into a job SLURM never schedules and snakemake never diagnoses. Verified on the fixture with max_mem_mb=700: "defect_map_merge: sized at 806 MB, capped to max_mem_mb=700" at parse time, mem_mb=700 on the job. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- workflow/README.md | 7 +++++- workflow/rules/exposure.smk | 39 +++++++++++++++++++++++++++--- workflow/scripts/defect_map_exp.py | 6 ++++- 3 files changed, 46 insertions(+), 6 deletions(-) diff --git a/workflow/README.md b/workflow/README.md index 2276f963c..db01d8207 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -242,7 +242,12 @@ profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; kee adding it once a campaign has produced one. The rule hangs off `exp_split`, not off `exp_psf`, so re-rasterizing the campaign at a different fidelity never touches the PSF chain, and `clean_exposure` takes its manifest as an - input, so reclamation cannot overtake the copy. Its resolution and its + input for a LIVE exposure, so reclamation cannot overtake the copy — and only + for a live one: an exposure whose store went to the /scratch purge (no + tombstone, nothing left to rasterize) is asked for its existing fragment if it + has one and for nothing if it does not, the same split `defect_map_merge`'s + input makes, because requiring a manifest behind a vanished split dir would + rebuild the whole exposure chain from VOS to reclaim it. Its resolution and its oversampling ride on `params`; `config.yaml`'s `defect_map:` block carries both, the measured convergence table behind the default, and the measurement on one real exposure (34 s, 0.62 GB, a 2.0 MB fragment; the rasterization is diff --git a/workflow/rules/exposure.smk b/workflow/rules/exposure.smk index f37c3f056..218708f9c 100644 --- a/workflow/rules/exposure.smk +++ b/workflow/rules/exposure.smk @@ -299,8 +299,32 @@ rule clean_exposure: # exposure's chain did not already put there. It is UNCONDITIONAL now: # exp_persist always packs the star catalogue's inputs, so there is no # keep list under which this rule has nothing to wait for. - lambda wc: [prod_exp_manifest(wc.exp, "exp_persist"), - prod_exp_manifest(wc.exp, "exp_defect_map")] + lambda wc: [prod_exp_manifest(wc.exp, "exp_persist")], + # The fragment must be off /scratch before the store goes too — the flag + # splits go with it — but this edge is CONDITIONAL, and it is exactly + # defect_map_inputs()' split (Snakefile), for exactly its reason. + # + # Naming the exp_defect_map manifest unconditionally reopens the + # avalanche that function is written to avoid. An exposure whose store + # went to the 60-day /scratch purge, or to a `clean: false` run, has NO + # TOMBSTONE — exp_store_reclaimed()'s docstring names that case — so + # clean_targets() still asks for one, and the manifest it would then + # require sits behind exp_split's manifest, which went with the store: + # snakemake schedules exp_get_images and exp_split from VOS, four hours + # per exposure, campaign-wide, on the first run of this branch. + # + # So: a LIVE exposure is asked for its manifest (the thing to build, and + # the thing that orders this rule after the rasterization); a RECLAIMED + # one is asked for its FRAGMENT if it has one — already on the + # persistent root, no rule's declared output, hence a leaf that requires + # nothing — and for nothing at all if it has neither, which is an + # exposure reclaimed by a workflow predating this rule and whose flags + # are gone either way. Blocking its tombstone would pin its scratch + # store forever without recovering a single flag. + lambda wc: ([prod_exp_manifest(wc.exp, "exp_defect_map")] + if not exp_store_reclaimed(wc.exp) + else [prod_exp_fragment(wc.exp)] + if Path(prod_exp_fragment(wc.exp)).exists() else []) output: tombstone = f"{EXP_DIR}/cleaned.json" params: @@ -456,8 +480,15 @@ rule defect_map_merge: # recorded coverage count once there is one and, before that, from a # per-exposure figure capped at the MEASURED DR6 footprint; the # Snakefile's sizing block carries both and argues the cap. - mem_mb = lambda wc, attempt: attempt * ( - DEFECT_MEM_BASE_MB + 2 * defect_map_cov_bytes() // 1_000_000), + # capped_mem() for the same reason star_cat_merge and final_cat_merge + # take it: defect_map_cov_bytes() scales as nside^2 and defect_map.nside + # is an advertised knob, so one ladder change turns this into a request + # no partition can schedule — a job that sits PENDING while the campaign + # looks alive, instead of a diagnosable OOM and a parse-time warning. + mem_mb = lambda wc, attempt: capped_mem( + attempt * (DEFECT_MEM_BASE_MB + + 2 * defect_map_cov_bytes() // 1_000_000), + "defect_map_merge"), # Dominated by reading fragments (~2 MB each) and setting their pixels; # ~1 s per exposure measured, over a floor that covers writing the map. # diff --git a/workflow/scripts/defect_map_exp.py b/workflow/scripts/defect_map_exp.py index 7473f42b3..37517d7ef 100644 --- a/workflow/scripts/defect_map_exp.py +++ b/workflow/scripts/defect_map_exp.py @@ -171,7 +171,11 @@ def ccd_files(exp_dir: Path, n_ccds: int) -> list: f"the {n_ccds} this exposure was split into; the split dir is " f"incomplete and a fragment built from it would be a hole in " f"the footprint marked complete (see ccd_files' docstring). " - f"Re-run exp_split for this exposure.") + f"Re-run exp_split for this exposure: delete its scratch " + f"manifests/exp_split.json and the workflow rebuilds the " + f"split. Until it does, clean_exposure waits on this rule " + f"and the store stays — a damaged store is not reclaimed " + f"silently.") return sorted(out) From f760687b900766abd1768137d739b0e5a98c03d4 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 9 Sep 2026 21:50:20 -0400 Subject: [PATCH 7/7] perf(defect-map): the merge allocates its accumulator once MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A rebuild set pixels into a fresh accumulator fragment by fragment, so every fragment that touched an unseen coverage pixel made healsparse grow — and copy — the sparse array. union_coverage() now reads each fragment's COVERAGE TABLE first (HealSparseCoverage.read, kilobytes, never the sparse array) and seeds make_empty(cov_pixels=...) with the union, so the final array is allocated once. MEASURED rather than assumed, on synthetic fragments at the campaign's own resolution (nside 131072 / coverage 128), 400 fragments over 5131 coverage pixels, a 656 MB accumulator: accumulation 8.4 s unseeded, 7.1 s seeded, with a 0.7 s coverage pre-pass. The reallocations are ~16% of the accumulation, not the dominant term a naive reading predicts — healsparse grows the array in blocks — so the review's tens-of-TB estimate does not hold, but the pre-pass is linear where the copying is not, and the win grows with the accumulator. Only the rebuild path seeds: an append starts from the map on disk, whose coverage is already most of the footprint. A fragment at the wrong coverage resolution abandons the seed so accumulate()'s named error still stands; the existing nside-mismatch test pins that. tests/unit/test_defect_map_reconcile.py passes (7/7, run in the container), and an end-to-end fixture merge reproduces rebuild == OR of fragments, no-op untouched, append, and removal-rebuild. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QbnPCyzuDNTgkg715pHhar --- workflow/scripts/merge_defect_map.py | 46 ++++++++++++++++++++++++++-- 1 file changed, 44 insertions(+), 2 deletions(-) diff --git a/workflow/scripts/merge_defect_map.py b/workflow/scripts/merge_defect_map.py index 2d7e13c58..385970cf6 100644 --- a/workflow/scripts/merge_defect_map.py +++ b/workflow/scripts/merge_defect_map.py @@ -182,6 +182,46 @@ def reconcile_plan(output: Path, sidecar: Path, have: dict) -> Plan: return Plan(sorted(set(have) - set(known)), [], "") +def union_coverage(paths, nside_coverage) -> np.ndarray: + """The coverage pixels every fragment in ``paths`` touches, together. + + A PRE-PASS, so the accumulator is allocated ONCE. ``target[pixels] = True`` + into a coverage pixel the map has not seen yet makes healsparse GROW its + sparse array, which copies it; at DR6 scale that array is gigabytes and a + rebuild discovers coverage pixels all the way through the campaign, so the + copies dominate everything the "reading fragments dominates" comment below + models. Seeding the coverage up front turns O(fragments) reallocations of a + growing array into one allocation of the final one. + + It costs a second read of each fragment's COVERAGE TABLE only — + ``HealSparseCoverage.read`` never touches the sparse array — which is + kilobytes against the megabytes the accumulation itself reads. + + MEASURED, on synthetic fragments at the campaign's own resolution (nside + 131072 / coverage 128), 400 fragments discovering 5131 coverage pixels — a + 656 MB accumulator: accumulation 8.4 s unseeded, 7.1 s seeded, with a 0.7 s + coverage pre-pass. So the reallocations are ~16% of the accumulation here, + not the dominant term a naive "copy the array once per new coverage pixel" + reading predicts (healsparse grows the sparse array in blocks). The win + grows with the accumulator; the pre-pass does not. Both paths produced + identical maps. + + Only the rebuild path uses it: an append starts from the map on disk, whose + coverage is already most of the footprint, and reads a handful of fragments. + """ + mask = None + for path in paths: + cov = hsp.HealSparseCoverage.read(str(path)) + if cov.nside_coverage != nside_coverage: + # accumulate() is the one place that reports a fragment built at the + # wrong resolution, with the exposure id and what to do about it. + # Here it is only a seed: give up on it and let that error stand. + return None + mask = (cov.coverage_mask.copy() if mask is None + else mask | cov.coverage_mask) + return None if mask is None else np.where(mask)[0] + + def accumulate(target, paths, nside_coverage, nside) -> None: """OR each fragment into ``target``, ONE AT A TIME (see the docstring). @@ -256,9 +296,11 @@ def apply_plan(output: Path, sidecar: Path, plan: Plan, have: dict, the previous PAIR intact rather than a map the record no longer describes. """ if plan.rebuild: - target = hsp.HealSparseMap.make_empty( - nside_coverage, nside, np.bool_, bit_packed=True) todo = plan.rebuild + target = hsp.HealSparseMap.make_empty( + nside_coverage, nside, np.bool_, bit_packed=True, + cov_pixels=union_coverage( + [have[exp] for exp in todo], nside_coverage)) else: target = hsp.HealSparseMap.read(str(output)) todo = plan.append