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/README.md b/workflow/README.md index 6b21ef097..db01d8207 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,49 @@ 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 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 + 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. 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 @@ -305,6 +350,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/Snakefile b/workflow/Snakefile index 64640f615..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 @@ -386,6 +387,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 +532,171 @@ 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 _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 + 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. + + 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. + """ + 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. + + 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 [] + 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(): + """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. + + 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: + return [] + return defect_map_inputs() + def persist_targets(): """Which exposures this invocation must pack PSF products off scratch for. @@ -833,6 +1001,24 @@ 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), 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): """Bytes on disk, or the documented per-unit default if it is not there.""" @@ -842,6 +1028,39 @@ 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. + + 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()), + DEFECT_COV_FOOTPRINT) + return n_cov * block + + def star_cat_max_bytes(): """The LARGEST exposure's psf_validation members — what sizes the merge. @@ -1107,8 +1326,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..f8a0687dd 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -180,6 +180,65 @@ 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 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 +# 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.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). 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` +# (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/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 diff --git a/workflow/rules/exposure.smk b/workflow/rules/exposure.smk index 5845adb8d..218708f9c 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,91 @@ 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 ~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. +# +# 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", "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, + 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 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, + # 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}" + + # --- 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 +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")] + 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: @@ -302,3 +419,95 @@ 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 + # 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_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. + # 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. + # + # 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" + " --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..37517d7ef --- /dev/null +++ b/workflow/scripts/defect_map_exp.py @@ -0,0 +1,342 @@ +#!/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. + +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 +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, n_ccds: int) -> list: + """``(ccd, flag path, image path)`` for every CCD this exposure split, in + 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 = [] + 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)) + 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: 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) + + +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() + + +# 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, 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 + # 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) + del data + if rows.size == 0: + return np.empty(0, dtype=np.int64) + 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: + """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("--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) + 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, args.n_ccds) + 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..385970cf6 --- /dev/null +++ b/workflow/scripts/merge_defect_map.py @@ -0,0 +1,376 @@ +#!/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. + +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 +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 +# 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: + """Where ``exp_defect_map`` wrote this exposure's fragment.""" + return (products_dir / "exp" / exp[:2] / exp / "defect" + / f"defect-{exp}.hsp") + + +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: + """``({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]) != sidecar_stamp(path)) + if changed: + return Plan([], sorted(have), + f"{len(changed)} fragment(s) changed on disk") + 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). + + 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 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. + + 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: + 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 + accumulate(target, [(exp, have[exp]) for exp in todo], + nside_coverage, nside) + + 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") + try: + target.write(str(map_tmp), clobber=True) + write_stable(map_tmp, output) + finally: + map_tmp.unlink(missing_ok=True) + write_sidecar(sidecar, record) + 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(): + # 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"{'; sidecar refreshed' if stale else ''})") + 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()