diff --git a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py index 7bc0cb76b..2b1bf8c34 100644 --- a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py +++ b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py @@ -239,10 +239,15 @@ def process(self): my_mask[inside_circle] = True for name in self._input_file_list: + # The source to read and the NAME to report it by; identical for a + # plain [path] entry (see MergeStarCatPSFEX's docstring on the + # [fileobj, name] form). This class takes its CCD numbers from the + # data's own CCD_ID_LIST, so the name is only ever used in messages. + source, label = name[0], name[-1] try: - starcat_j = fits.open(name[0], memmap=False, ignore_missing_simple=True) + starcat_j = fits.open(source, memmap=False, ignore_missing_simple=True) except ValueError: - print(f"Error for file {name[0]}, check FITS file integrity") + print(f"Error for file {label}, check FITS file integrity") #raise continue @@ -524,7 +529,15 @@ class MergeStarCatPSFEX(object): Parameters ---------- input_file_list : list - Input files + Input entries. Each entry is a list, as the module runner builds them: + ``[path]`` from the file handler. An entry may also carry a name + alongside an already-open source, ``[fileobj, name]`` — ``fits.open`` + takes the first element and the CCD number is parsed from the LAST, + which is the same string in the one-element case. That is what lets a + caller merge catalogues it never wrote to disk (the Snakemake + workflow's ``star_cat_merge`` reads them out of the per-exposure tars + with ``tarfile`` + ``BytesIO``), without this class learning anything + about where they came from. output_dir : str Output directory w_log : logging.Logger @@ -569,10 +582,15 @@ def process(self): ) for name in self._input_file_list: + # The source to read and the NAME to parse the CCD number out of. + # Identical for a plain [path] entry; different only when the caller + # hands over an open file-like object plus the member name it came + # under (see the class docstring). + source, label = name[0], name[-1] try: - starcat_j = fits.open(name[0], memmap=False, ignore_missing_simple=True) + starcat_j = fits.open(source, memmap=False, ignore_missing_simple=True) except OSError as e: - print(f"Error while opening file '{name[0]}'") + print(f"Error while opening file '{label}'") #raise continue @@ -614,7 +632,7 @@ def process(self): psfex_acc += list(np.zeros_like(data_j["X"])) # CCD number - ccd_nb += [re.split(r"\-([0-9]*)\-([0-9]+)\.", name[0])[-2]] * len( + ccd_nb += [re.split(r"\-([0-9]*)\-([0-9]+)\.", label)[-2]] * len( data_j["RA"] ) @@ -786,7 +804,11 @@ def process(self): ) for name in self._input_file_list: - starcat_j = fits.open(name[0], memmap=False) + # The source to read and the NAME to parse the CCD number out of; + # identical for a plain [path] entry (see MergeStarCatPSFEX's + # docstring on the [fileobj, name] form). + source, label = name[0], name[-1] + starcat_j = fits.open(source, memmap=False) data_j = starcat_j[self._hdu_table].data @@ -811,7 +833,7 @@ def process(self): snr += list(data_j["SNR_WIN"]) # CCD number - ccd_nb += [re.split(r"\-([0-9]*)\-([0-9]+)\.", name[0])[-2]] * len( + ccd_nb += [re.split(r"\-([0-9]*)\-([0-9]+)\.", label)[-2]] * len( data_j["XWIN_IMAGE"] ) diff --git a/workflow/README.md b/workflow/README.md index 55611cedc..eb801da6a 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -107,7 +107,9 @@ and the run fails if either phase failed. ## The launch code snapshot `sp run` copies the code it is about to launch — `workflow/` (config symlinks -dereferenced), `src/` and the profile — into `/code`, records HEAD +dereferenced), `src/`, the repo's `scripts/` (`final_cat_merge` loads +`scripts/python/create_final_cat.py` by path) and the profile — into +`/code`, records HEAD plus a dirty flag in `/code/snapshot.json`, and runs the campaign entirely out of that copy. It matters because a campaign is not one process: the SLURM executor re-invokes snakemake on every job's node, so jobs re-parse the @@ -156,15 +158,17 @@ 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 (no temp()) - tile.smk per-tile: exp forest, merge_headers, detect, vignets, ngmix, merge, make_cat + exposure.smk per-exposure: get_images, split, psf, persist (no temp()); campaign star_cat_merge + tile.smk per-tile: exp forest, merge_headers, detect, vignets, ngmix, merge, make_cat; campaign final_cat_merge scripts/ - sp_rule.py the thin per-unit wrapper (isolation furniture, config copy, log-sync, count check) build_index.py prepare-phase run_index.sqlite builder (plain script) build_forest.py per-tile exposure symlink forest (group-compatible shell) completeness.py the ported count table (shared by sp_rule + run_report) run_report.py standalone report (NOT a DAG node; run_report hooks call it) container.py image layers + the resolution order behind `sp container` (stdlib-only) + persist_exp.py ONE exposure's keepable PSF products -> one tar on products_dir (the exp_persist rule) + merge_star_cat.py ALL exposures' validation_psf, read out of the tars -> full_starcat (the star_cat_merge rule) + 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) profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; keep-going ``` @@ -240,6 +244,54 @@ profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; kee trigger reads that cut as a reason to rerun the very tiles it protects. Know the consequence — `--forcerun` on a tile whose `final_cat` exists will not rebuild its reclaimed exposures. Delete the `final_cat` first. +- **PSF products leave scratch before the purge does.** `exp_persist` packs + the files named by `persist_exp:` in `config.yaml` (default: the psfex_interp + `validation_psf-*.fits`, the rho/tau statistics input) from the exposure's + scratch store into ONE uncompressed tar, + `/exp///psf/.tar` (inodes, not bytes, bind + on /project), and writes ONE manifest beside it recording the patterns, the + members and their sizes. The + threat it answers is the /scratch purge, not `clean_exposure` — the store goes + in 60 days whether or not the workflow reclaimed it — so it runs even with + `clean: false`, requested directly by `rule all`. `clean_exposure` takes its + manifest as an input, so reclamation can never overtake the copy. It is a + rule of its own rather than a `cp` on the end of `exp_psf` because the keep + list rides on `params`: adding a pattern reruns seconds of packing, not four + hours of PSF fitting per exposure. A pattern that matches nothing is a + recorded warning (setools rejects sparse CCDs); matching nothing at all is a + failure. A `localrule`, by the same arithmetic as `clean_exposure`. +- **The campaign ends in two merged catalogues, and the workflow now makes + both.** 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. + `star_cat_merge` stacks every exposure's every CCD's `validation_psf-*.fits` + into one `/full_starcat-0000000.fits` — the rho/tau statistics + input, at the path sp_validation hardcodes. It reads the members straight out + of the per-exposure tars (`tarfile` + `BytesIO`; unpacking ~800k files to + merge them would defeat the tar's whole purpose) and stacks them with + `MergeStarCatPSFEX`, the same class the old `merge_starcat_runner` called, so + the column list has exactly one definition. Its input is the same + `exp_persist` manifest set `rule all` already requests, so it pulls nothing + new into the DAG, and it exists only when `persist_exp:` keeps a + `validation_psf-*.fits`-shaped file — otherwise no job, and a warning at parse + time rather than a failure on a node. + `final_cat_merge` collects every ready tile's `final_cat-.fits` into + `/final_cat_.hdf5`: one dataset per tile under a group + named for the campaign, the `final_cat.param` columns, an `n_tiles` attribute. + That schema is what sp_validation's reader opens, so it is fixed; the column + extraction reuses `scripts/python/create_final_cat.py` while the file is + written here, because that script's own discovery walks a directory layout + this workflow does not have. `campaign:` in `config.yaml` names the group and + defaults to the persistent root's basename. + Both rebuild from the whole persistent root rather than appending, so the + output is a function of its input set: byte-stable on a no-op rerun + (tmp-then-`cmp`-then-`mv`), and rebuilt when a tile or exposure is appended + (the input list's fingerprint rides on `params`). Neither is a `localrule` — + one job over ~20k units is real work — and neither puts its input paths in its + shell, which is not fastidiousness: ~20k paths is an order of magnitude over + Linux's 128 KiB `MAX_ARG_STRLEN` for a single argv entry, so each script + rediscovers the set under `products_dir` while the fingerprint travels on + `params`. - **A dead tile can be told to stop pinning exposures.** An exposure is cleanable only once every consuming tile has its vignets, so one permanently-failed tile holds its ~80 exposures for the life of the diff --git a/workflow/Snakefile b/workflow/Snakefile index bbdfeefb6..587f271b2 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 fnmatch import functools import hashlib import json @@ -40,6 +41,10 @@ import sys from pathlib import Path from snakemake.exceptions import WorkflowError +# Explicit rather than relying on the name snakemake injects into this +# namespace: the one place we log at parse time is a branch that only a +# non-default keep list reaches, and a NameError there would be found by a user. +from snakemake.logging import logger # Resolved relative to THIS file, not the working directory: snakemake runs with # --directory on /scratch (bin/sp) so .snakemake/ state never lands on /project @@ -92,6 +97,14 @@ RUN_DIR = Path(OUTPUTS["run_dir"]) # second path: one root, exactly the pre-D5 layout. PRODUCTS_DIR = Path(OUTPUTS.get("products_dir") or RUN_DIR) INDEX_DB = Path(OUTPUTS["index_db"]) +# The campaign's NAME — what the two campaign-level merges label their output +# with (`final_cat_.hdf5`, and the group inside it that holds the +# campaign's per-tile datasets). It +# defaults to the persistent root's basename, which is already how every +# campaign here is named (smk-g4, smk-g5, smk-g6: run_dir, products_dir and +# index all end in it), so the common case needs no key at all. Set `campaign:` +# in config.yaml when the two must differ. +CAMPAIGN = config.get("campaign") or PRODUCTS_DIR.name SCRIPTS = Path(workflow.basedir) / "scripts" # The config chain is the repo's committed directory (D2). The configs and # rules that set their environment variables must be versioned together. There @@ -123,8 +136,14 @@ if PHASE not in ("prepare", "compute", "passthrough"): raise WorkflowError( f"SP_PHASE={PHASE!r} is not one of prepare, compute, passthrough.") +# DEDUPED, order preserved. The tile list is appended to by hand across a +# campaign, so a tile can appear twice — harmless for a per-tile target, which +# is the same path requested twice, but not for the campaign-level merges: their +# fingerprint counts what is in this list and the job derives a deduped set from +# the same file (build_index.campaign_tiles), so a duplicate line would make the +# two disagree about a set they must name identically. with open(config["tile_list"]) as f: - TILES = [ln.strip() for ln in f if ln.strip()] + TILES = list(dict.fromkeys(ln.strip() for ln in f if ln.strip())) # --- ngmix scatter (D4) ---------------------------------------------------- # Native directive: `--set-scatter ngmix=N` overrides it, N=1 degenerates to one @@ -263,6 +282,7 @@ EXP_DIR = str(RUN_DIR / "exp" / "{shard}" / "{exp}") # The persistent root mirrors the scratch one, shard for shard, so the two trees # read as the same campaign seen from two filesystems. PROD_TILE_DIR = str(PRODUCTS_DIR / "tiles" / "{shard}" / "{tile}") +PROD_EXP_DIR = str(PRODUCTS_DIR / "exp" / "{shard}" / "{exp}") def tile_dir(tile): return f"{RUN_DIR}/tiles/{tile[:2]}/{tile}" @@ -276,6 +296,18 @@ def tile_manifest(tile, stage): def exp_manifest(exp, stage): return f"{exp_dir(exp)}/manifests/{stage}.json" +def prod_exp_dir(exp): + """The exposure's dir on the PERSISTENT root — where exp_persist writes. + + Sharded identically to the scratch one, so the two trees read as the same + campaign seen from two filesystems, exposure side as well as tile side.""" + return f"{PRODUCTS_DIR}/exp/{exp[:2]}/{exp}" + +def prod_exp_manifest(exp, stage): + """A manifest that must SURVIVE reclamation, so it is not in the exposure's + scratch manifests/ dir (clean_exposure deletes that wholesale).""" + return f"{prod_exp_dir(exp)}/manifests/{stage}.json" + def forest_dir(tile): return f"{tile_dir(tile)}/exp_forest" @@ -314,6 +346,9 @@ SCRIPT_HASH = script_hash("completeness.py") FOREST_HASH = script_hash("build_forest.py") CLEAN_HASH = script_hash("clean_exposure.py") CLEAN_TILE_HASH = script_hash("clean_tile.py") +PERSIST_HASH = script_hash("persist_exp.py") +MERGE_STAR_HASH = script_hash("merge_star_cat.py") +MERGE_FINAL_HASH = script_hash("merge_final_cat.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 @@ -431,6 +466,266 @@ def clean_targets(): out.append(tombstone(exp)) return sorted(out) +# --- persisted exposure products (D5) -------------------------------------- +# The keep list is config, not a rule input, and it is READ HERE so that exactly +# one place converts it into the form the rule carries. An empty list is a +# deliberate "keep nothing" and produces no jobs at all. +PERSIST_EXP = list(config.get("persist_exp") or []) + + +def persist_targets(): + """Which exposures this invocation must pack PSF products off scratch for. + + `rule all` requests these DIRECTLY rather than reaching them only through + clean_exposure. Persistence and reclamation are different concerns — the + /scratch purge takes the store whether or not `clean:` is on — and hanging + the copy off the clean rule alone would mean a campaign run with clean:false + persists nothing and loses everything at the purge. + + Scope is the ready tiles' exposures, which `all` already builds through the + tile chain, so nothing new is pulled into the DAG by asking. + + EXCEPT AN EXPOSURE WHOSE STORE IS GONE. Its exp_psf manifest is not there, + so requesting its persist manifest would make the DAG rebuild the whole + exposure chain from VOS — the avalanche tile.smk's reclaimed-edge cut exists + to prevent, arriving through a new target instead. exp_store_reclaimed() + below is that test, and what it does NOT test is the tombstone. + + HEAD PROCESS ONLY, for the same reason as clean_targets() above. + """ + if not workflow.is_main_process: + return [] + return persist_manifests() + + +@functools.lru_cache(maxsize=1) +def persist_manifests(): + """persist_targets() without the head-process guard, memoised. + + The guard on persist_targets() is a cost decision, not a correctness one: + `rule all` reads it at MODULE level, so a job parse would pay a whole + campaign's index walk for a target it can never schedule. star_cat_merge + reads the same list through an INPUT FUNCTION, which snakemake evaluates + only for parses that actually build that job — the head process, and the one + merge job's own re-parse under the slurm executor, which genuinely needs it. + So this half carries no guard and the memo keeps either parse to one walk. + """ + if not PERSIST_EXP: + return [] + exps = {e for t in TILES_READY for e in tile_exposures(t)} + return sorted(prod_exp_manifest(e, "exp_persist") for e in exps + if not exp_store_reclaimed(e)) + + +def exp_store_reclaimed(exp): + """True when this exposure's PSF products exist ONLY on the persistent root. + + The one condition both the persist target list and star_cat_merge's edge + choice turn on, and it deliberately does NOT read the tombstone. + + THE TOMBSTONE ALONE IS NOT THE EVIDENCE. It says clean_exposure ran, and + clean_exposure is only one of the two ways a scratch store disappears. + + WHAT THE TEST HAS TO SEPARATE is a store that is GONE from one that has not + been BUILT yet, and no single file says that. This runs at parse time, before + any exp_psf job of a fresh campaign has run, so "the exp_psf manifest is + missing" alone would skip every exposure of a new campaign and persist + nothing at all. The question is therefore: is there evidence this exposure + once had a store? Two files carry it, and either will do: + + * the TOMBSTONE — clean_exposure ran, so the store was built and reclaimed, + and whatever was going to be packed was packed before it went; + * a PERSISTED MANIFEST with no exp_psf manifest beside it — exp_persist + ran, so the store existed, and it is not there now. This is the purge + case, and it is the one the tombstone cannot see: /scratch is purged on + a 60-day window whether or not this workflow reclaimed anything, and it + leaves nothing behind. Keying on the tombstone alone meant that after a + purge, or on any campaign run with `clean: false`, every exposure looked + live, its persist manifest was requested, its exp_psf manifest was not + there, and snakemake rebuilt the entire exposure chain from VOS. + + An exposure with a LIVE store and a manifest is not reclaimed and is still + asked for, which is what lets an edit to `persist_exp:` re-pack in seconds + rather than be silently ignored — the whole reason exp_persist is a rule of + its own. An exposure with neither file is asked for too: either it has not + run yet, or it was purged having saved nothing, and only the DAG can tell + those apart by trying. + """ + return (Path(tombstone(exp)).exists() + or (Path(prod_exp_manifest(exp, "exp_persist")).exists() + and not Path(exp_manifest(exp, "exp_psf")).exists())) + +# --- the campaign-level merges --------------------------------------------- +# Two rules, one job each per campaign, both writing to the persistent root, and +# both the LAST link of a chain whose per-unit half the workflow already had: +# the exposure side ends in one `full_starcat-0000000.fits` (every CCD's PSF +# validation catalogue, stacked — the rho/tau statistics input) and the tile side +# in one `final_cat_.hdf5` (every tile's final catalogue — the shear +# catalogue sp_validation reads). Until they existed the workflow's product set +# was two files short of what the old `combine_runs.bash` + `create_final_cat.py` +# chain delivered, and every campaign ended with a manual merge. +# +# NEITHER IS A LOCALRULE, and the arithmetic runs the opposite way from +# exp_persist's. Those rules are ~20k jobs of seconds each, so submitting them +# costs more in scheduling latency than the work; these are ONE job each over the +# whole campaign — ~800k catalogues stacked in memory, or ~20k catalogues read +# end to end at DR6 scale. That is a compute job, and it belongs on a node. +# +# NEITHER PUTS ITS INPUT PATHS IN ITS SHELL. `{input}` at DR6 scale is ~20k paths +# in a single argv entry, an order of magnitude over Linux's 128 KiB +# MAX_ARG_STRLEN, and the job would die on exec. So each rule's `input` is the +# DAG EDGE (what must exist first) and each script rediscovers the same set from +# the tile list and the index; what travels is a FINGERPRINT of that set's unit +# ids, on `params`, which is what makes the merge rerun when the set changes and +# not otherwise. +# The scripts' docstrings argue the rediscovery — it is also what lets the merges +# cover exposures whose scratch stores reclamation has since taken. + + +def unit_fingerprint(units): + """A short digest of a set of UNIT IDS, for a merge rule's `params`. + + `params` is a rerun trigger and a set of ids is not: appending a tile grows + the set, moves the digest and reruns the merge, while a rerun over the same + set leaves it alone. Sorted before hashing because the ORDER is not part of + what changed. + + IDS RATHER THAN THE RULE'S `input` PATHS. A path can change while the set + does not — star_cat_merge's edge for one exposure flips from its manifest to + its tar when the store is reclaimed — and a merge that reruns over identical + content on every reclamation pass is a rerun trigger firing on bookkeeping. + The ids are also exactly what the job derives on its own side, so both + halves agree on the set and on how it is named. + """ + joined = "\n".join(sorted(str(u) for u in units)) + return f"{len(units)}:{hashlib.md5(joined.encode()).hexdigest()[:12]}" + + +# The tar member the star merge consumes. The keep list is globs, so the test is +# "would this member be kept", not a string comparison — `validation_psf-*.fits`, +# `validation_psf*`, `*.fits` and a bare `*` all say yes, and all are things a +# user might reasonably write. +_STAR_CAT_MEMBER = "validation_psf-2605805-12.fits" +STAR_CAT_MERGE = any(fnmatch.fnmatch(_STAR_CAT_MEMBER, p) for p in PERSIST_EXP) + +# LOUD AT PARSE TIME, once, and only where it can be acted on. A keep list +# without the validation catalogues is a legitimate configuration (persist the +# PSF models alone, say) — it is not an error, so it must not become a job that +# fails on a node an hour later. It is worth SAYING, because the omission is +# silent otherwise and the missing product only surfaces when a rho-statistics +# run cannot find its input. +if PERSIST_EXP and not STAR_CAT_MERGE and workflow.is_main_process \ + and PHASE == "compute": + logger.warning( + f"star_cat_merge: no job — persist_exp {PERSIST_EXP} keeps no " + f"'{_STAR_CAT_MEMBER}'-shaped file, so there is nothing to stack into " + f"{PRODUCTS_DIR}/full_starcat-0000000.fits (the rho/tau statistics " + f"input). Add 'validation_psf-*.fits' to persist_exp: to get it.") + + +def full_starcat(): + """The campaign's merged star catalogue. The NAME is not ours to choose: + sp_validation hardcodes `full_starcat-0000000.fits` beside its data dir.""" + return f"{PRODUCTS_DIR}/full_starcat-0000000.fits" + + +def final_cat_hdf5(): + """The campaign's merged shear catalogue — sp_validation's galaxy_cat_path.""" + return f"{PRODUCTS_DIR}/final_cat_{CAMPAIGN}.hdf5" + + +def prod_exp_tar(exp): + """The tar exp_persist writes. Not a declared output of anything — see + star_cat_inputs().""" + return f"{prod_exp_dir(exp)}/psf/{exp}.tar" + + +@functools.lru_cache(maxsize=1) +def star_cat_inputs(): + """What star_cat_merge waits for: every exposure of TILES_READY whose PSF + products are on the persistent root, live and reclaimed alike. + + THE SAME SET merge_star_cat.py derives at job time, and that equality is + load-bearing — the fingerprint on `params` is taken over THIS list, so + anything the job stacked that was not in it would be rows no rerun trigger + could see. The job states the rule from its own side: same tile list, same + index, exp_persist manifest present on the persistent root. By the time it + runs, every exposure below has one. + + RECLAIMED EXPOSURES BELONG IN THE STAR CATALOGUE. Carrying their PSF + products off scratch is exactly what exp_persist is for, and a merge that + dropped them would shrink the campaign's star catalogue every time + reclamation ran. But their exp_psf manifest is gone, so REQUESTING their + exp_persist manifest rebuilds the whole exposure chain from VOS — the + avalanche persist_targets() drops them to avoid. + ancient() DOES NOT HELP: it suppresses the timestamp comparison, not the + missing input, and snakemake schedules the chain anyway. Measured on smk-g6 + with one reclaimed exposure given a manifest by hand: the dry run grew + exp_get_images, exp_split, exp_psf and exp_persist jobs. + + So a reclaimed exposure is depended on through its TAR instead. The tar is + not a declared output of any rule (exp_persist declares only its manifest, + deliberately — persist_exp.py says why), so a tar that exists is a DAG leaf: + snakemake requires it and builds nothing. A live exposure keeps its manifest + edge, which is what orders the merge after the packing; its tar does not + exist yet, so it could not serve as the edge. + + An exposure reclaimed by a workflow PREDATING exp_persist has neither tar nor + manifest and is in no set at all. Nothing short of rebuilding its chain from + VOS recovers it; the merge reports how many exposures it found. + """ + if not PERSIST_EXP: + return [] + live, reclaimed = [], [] + for exp in sorted({e for t in TILES_READY for e in tile_exposures(t)}): + if not exp_store_reclaimed(exp): + live.append(prod_exp_manifest(exp, "exp_persist")) + elif Path(prod_exp_tar(exp)).exists(): + reclaimed.append(prod_exp_tar(exp)) + return live + reclaimed + + +@functools.lru_cache(maxsize=1) +def star_cat_exposures(): + """The exposure IDs star_cat_merge stacks — what its fingerprint is taken + over. + + THE IDS, NOT THE PATHS, and the difference is a rerun. An exposure's edge + FLIPS from its manifest to its tar the moment its store is reclaimed, so a + fingerprint over paths moves on every reclamation pass and reruns the merge + over content that did not change. The ids move only when the set does, which + is what the trigger is for. It is also what merge_star_cat.py derives on the + job side, so the two agree on the set AND on how it is named. + """ + if not PERSIST_EXP: + return [] + return sorted(e for e in {e for t in TILES_READY for e in tile_exposures(t)} + if Path(prod_exp_manifest(e, "exp_persist")).exists() + or not exp_store_reclaimed(e)) + + +def star_cat_targets(): + """`full_starcat` when there is anything to stack into it, else nothing. + + Three ways to get nothing, and all three are states rather than errors: the + keep list holds no validation catalogue (warned about above), `persist_exp:` + is empty at all, or every exposure in scope is already tombstoned — a + campaign resumed after reclamation, whose exposures were cleaned by a + workflow that predates exp_persist and therefore left neither tar nor + manifest to read. A rule with an empty input list would still be a JOB, and + it would write an empty star catalogue over a good one. + """ + if not STAR_CAT_MERGE or not workflow.is_main_process: + return [] + return [full_starcat()] if star_cat_inputs() else [] + + +def final_cat_targets(): + """The merged hdf5, whenever this campaign has a tile to put in it.""" + if not workflow.is_main_process or not TILES_READY: + return [] + return [final_cat_hdf5()] + # --- tile reclamation (D5) -------------------------------------------------- # A separate flag from `clean:` (config.yaml carries the full # argument): exposure reclamation costs nothing but a rebuild if a tile is @@ -601,11 +896,21 @@ include: "rules/tile.smk" # localrule would: a local job cannot be fused into a submitted group. The old # star-catalogue rules were exactly that, and they are gone with the internal # mask generation.) -localrules: all, prepare_all_tiles, clean_exposure, clean_tile +# +# exp_persist joins them for the same arithmetic — one tar of a few MB per +# exposure, ~20k of them at DR6 scale, each far shorter than the scheduling +# latency that would submit it (exposure.smk argues the placement in full). It +# sits mid-chain between exp_psf and clean_exposure, but both of those are +# outside every group already (exp_psf is heavy, clean_exposure is local), so it +# adds no new grouping constraint. +localrules: all, prepare_all_tiles, clean_exposure, clean_tile, exp_persist rule all: input: [final_cat(t) for t in TILES_READY], + persist_targets(), + star_cat_targets(), + final_cat_targets(), clean_targets(), clean_tile_targets(), diff --git a/workflow/bin/sp b/workflow/bin/sp index fdc3d4f59..e3aec91a4 100755 --- a/workflow/bin/sp +++ b/workflow/bin/sp @@ -72,7 +72,10 @@ STATE_DIR="${SP_STATE_DIR:-${RUN_DIR}-state}"; mkdir -p "$STATE_DIR" # WHAT. `sp run` copies the code it is about to launch into $STATE_DIR/code and # runs the campaign entirely out of that copy: the Snakefile, the rules, the # scripts, the ini chain (symlinks DEREFERENCED -- workflow/config/cfis points -# into example/, and the copy must be self-contained), src/, and the profile. +# into example/, and the copy must be self-contained), src/, the repo's own +# scripts/ (final_cat_merge loads scripts/python/create_final_cat.py by path -- +# it is a script, not an installed module, and the hdf5 layout it defines must +# be pinned to the campaign like everything else here), and the profile. # Every workflow-internal path hangs off `workflow.basedir`, which IS the # snapshot, so they all follow it for free; the profile's PYTHONPATH pin is the # one that cannot (YAML splices nothing) and is rewritten below. @@ -91,10 +94,10 @@ snapshot_code() { mkdir -p "$SNAPSHOT" if command -v rsync >/dev/null 2>&1; then rsync -a --delete --copy-links --exclude '__pycache__' --exclude '*.egg-info' \ - "$HERE" "$REPO/src" "$REPO/profiles" "$SNAPSHOT/" + "$HERE" "$REPO/src" "$REPO/scripts" "$REPO/profiles" "$SNAPSHOT/" else rm -rf "$SNAPSHOT"; mkdir -p "$SNAPSHOT" - cp -rL "$HERE" "$REPO/src" "$REPO/profiles" "$SNAPSHOT/" + cp -rL "$HERE" "$REPO/src" "$REPO/scripts" "$REPO/profiles" "$SNAPSHOT/" find "$SNAPSHOT" -name __pycache__ -type d -prune -exec rm -rf {} + fi diff --git a/workflow/config.yaml b/workflow/config.yaml index 63b29413c..70ace906e 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -57,6 +57,19 @@ outputs: # (-state; bin/sp explains why). products_dir: /project/def-mjhudson/cdaley/sp-products/smk-g6 +# The campaign's NAME. It labels the two campaign-level merges' output — +# /final_cat_.hdf5 and the group inside it holding that +# campaign's per-tile datasets — +# and nothing else; the per-unit stores are named by their own IDs. UNSET means +# the persistent root's basename, which is already how every campaign here is +# named (run_dir, products_dir and index_db all end in smk-g6), so this key only +# earns its place when the two must differ. +# +# It is not a rule input, so renaming a campaign mid-flight changes the merged +# catalogue's PATH and therefore builds a new one; the per-tile catalogues it +# reads are untouched. +# campaign: smk-g6 + # There is no config_src knob: the config chain is workflow/config/cfis, resolved # relative to the Snakefile. The configs interpolate $SP_RUN / $SP_UNIT_NUM / # $SP_CONFIG / $SP_EXP / $NGMIX_* and the rules export them -- configs and rules @@ -68,6 +81,88 @@ outputs: # would otherwise have to rebuild from tile headers. index_db: /project/def-mjhudson/cdaley/sp-products/smk-g6/index/run_index.sqlite +# Per-exposure PSF products to carry onto the persistent root before the scratch +# store goes (`exp_persist`, exposure.smk). A list of plain file-name globs, +# matched recursively under the PSF chain's four module output dirs +# (/exp///output/run_sp_exp_SxSePsfPi/*/output/ — +# sextractor_runner, setools_runner, psfex_runner, psfex_interp_runner). +# Matches are packed, flat, into ONE uncompressed tar per exposure: +# /exp///psf/.tar, with a manifest listing the +# members beside it. One tar rather than loose copies because inodes, not bytes, +# bind on /project (~1 M-file group quota; loose copies would be ~200 files per +# exposure, ~2 M at DR6 scale). FITS members read straight from the tar: +# fits.open(io.BytesIO(tarfile.open(t).extractfile(m).read())). +# +# WHY COPY RATHER THAN EXEMPT THESE FROM CLEANUP. Reclamation is not the threat. +# run_dir is /scratch and is PURGED on a 60-day window whether or not +# clean_exposure ever ran; products_dir is /project, backed up and not purged. +# The only way a per-exposure product outlives its campaign is to leave the +# filesystem. (Ordering is free: clean_exposure takes the exp_persist manifest +# as an input, so a store is never reclaimed before its keepers are written.) +# +# EDITING THIS LIST IS CHEAP. It rides on exp_persist's `params`, so a change +# reruns the packing (seconds) and NOT exp_psf (four hours per exposure). That +# separation is the whole reason exp_persist is a rule of its own. +# +# The default is the minimum: the psfex_interp VALIDATION catalogue, one per +# CCD, which is the input to the rho/tau statistics. Without it the PSF +# diagnostics cannot be recomputed after a purge without rebuilding the exposure +# chain from VOS. +# +# THIS LIST GATES `star_cat_merge`. That campaign-level rule stacks every +# exposure's every CCD's validation_psf into ONE +# /full_starcat-0000000.fits, reading the members straight out of +# the tars. A keep list that matches no `validation_psf-*.fits` is a legitimate +# configuration (keep the PSF models alone, say) and produces NO merge job and a +# warning at parse time — not a failure on a node an hour later. Note the +# corollary: an exposure already reclaimed by a workflow that predates +# exp_persist left no tar, so it contributes nothing and cannot be recovered +# short of rebuilding its chain from VOS. +# +# OPT-IN CANDIDATES, and what each buys. Sizes are per exposure (40 CCDs), +# measured on smk-m2 (127 exposures, 64 tiles); a 64-tile campaign with all of +# the measured ones on came to 7.2 GB: +# validation_psf-*.fits (the default) 2.0 MB +# *.psf the PSFEx model itself. Keeping it means the PSF +# can be re-interpolated at ANY position later +# without rebuilding the exposure chain — the +# single most capability-adding entry here. +# 2.8 MB +# psfex_cat-*.cat PSFEx's own output catalogue (FITS_LDAC): the +# per-star FLAGS_PSF / CHI2_PSF, i.e. WHICH stars +# outlier rejection clipped. Not recoverable from +# anything else (the .psf header keeps only the +# LOADED/ACCEPTED counts). unmeasured +# star_selection-*.fits the PRE-SPLIT selection (setools writes it under +# mask/). The only file that can answer "which +# stars were rejected by the selection cuts, and +# why" — the split samples have already lost the +# rejects. 24.5 MB +# star_split_ratio_80-*.fits setools' 80% TRAINING star sample, the set PSFEx +# actually fitted. Rows duplicate star_selection. +# 19.9 MB +# star_split_ratio_20-*.fits the 20% VALIDATION sample — the positions the +# validation_psf rows correspond to. Rows +# duplicate star_selection. 7.1 MB +# star_stat-*.txt setools' per-CCD STAT block (star counts, +# stars/deg^2, FWHM mode and cuts, under stat/): +# the selection's summary without its catalogue. +# unmeasured +# A production keep list is `validation_psf` + `*.psf` + `psfex_cat` (~5 MB per +# exposure); the star_split files are only worth it if star_selection is off. +# PSFEx residual/check images and its XML diagnostics are NOT candidates as the +# chain stands: the committed default.psfex sets CHECKIMAGE_TYPE NONE and +# WRITE_XML N, so nothing is emitted to match. They are a config change first, +# a pattern second. +# +# NOTE ON products_dir DEFAULTING TO run_dir (a fixture or smoke test): the tar +# then lands beside the store on the same filesystem and buys nothing, and the +# manifest sits in the exposure's own manifests/ dir, which clean_exposure +# deletes wholesale — so a one-root run re-persists after every reclamation. +# Harmless, and exactly the pre-D5 behaviour a one-root run asks for. +persist_exp: + - validation_psf-*.fits + # 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 4e1ef45a4..7b6d33882 100644 --- a/workflow/rules/exposure.smk +++ b/workflow/rules/exposure.smk @@ -1,6 +1,6 @@ """Exposure chain — per exposure, keyed by exp base id (dedup is structural). - exp_get_images -> exp_split -> exp_psf + exp_get_images -> exp_split -> exp_psf -> exp_persist 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 @@ -17,6 +17,13 @@ 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 +/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 +argument is in workflow/scripts/persist_exp.py. + NO temp() anywhere in this file, ever (D5). Exposures overlap tiles by construction (~7-10 tiles each), so their consumer set closes over the CAMPAIGN, not over one invocation — reclamation here is clean_exposure's job (S5), driven @@ -106,6 +113,59 @@ rule exp_psf: sp_shell("exp_psf", f"config_exp_{PSF_MODEL}.ini") +# --- persistence (D5) ------------------------------------------------------- +# The counterpart of reclamation, and it must come first in the DAG: this packs +# the exposure's keepable PSF products into one tar on the persistent root, and +# clean_exposure below takes its manifest as an input so the store is never +# reclaimed before the keepers have left /scratch. The purge would take them +# anyway — that, not clean_exposure, is what this rule exists for +# (persist_exp.py's docstring argues both halves, and config.yaml's +# `persist_exp:` block carries the keep list and its candidates). +# +# A LOCALRULE (declared in the Snakefile), by exactly the arithmetic that made +# clean_exposure one: the body is a `tar` of a few MB from one shared filesystem +# to another, seconds of work, and one sbatch per exposure would be ~20k +# submissions at DR6 scale for jobs shorter than the scheduling latency. The +# grouping constraint that binds mid-chain localrules (this file's docstring) +# does not bite here: exp_persist's only neighbours are exp_psf, which is too +# heavy to ever fuse, and clean_exposure, which is local itself. +# +# ONE DECLARED OUTPUT, AND IT IS A MANIFEST, NOT THE TAR OR A directory(). The +# tar is not declared: a directory output would attest that a directory exists, +# where what we want written down is WHICH files were packed and how big each was — +# the provenance a rho-statistics run months from now needs in order to know +# what it is reading. The manifest is byte-stable, so a no-op rerun does not +# move its mtime and does not make clean_exposure look out of date. +# +# THE KEEP LIST RIDES ON params. That is the entire reason this is not three +# lines of tar appended to exp_psf's shell: `params` is a rerun trigger, so +# adding a pattern reruns the packing and leaves the PSF chain alone. +rule exp_persist: + input: + rules.exp_psf.output.manifest + output: + manifest = f"{PROD_EXP_DIR}/manifests/exp_persist.json" + # No `log:`: the script's only failure modes are "nothing matched" and a + # name collision, both of which it reports on stderr and neither of which + # has a per-CCD verdict worth a completeness record. + params: + patterns = " ".join(f"--pattern '{p}'" for p in PERSIST_EXP), + exp_dir = lambda wc: exp_dir(wc.exp), + dest = lambda wc: f"{prod_exp_dir(wc.exp)}/psf", + script_hash = PERSIST_HASH + threads: 1 + retries: 2 + resources: + mem_mb = 2000, + runtime = 10 + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/persist_exp.py" + " --exp-dir '{params.exp_dir}' --exp {wildcards.exp}" + " --dest '{params.dest}' --manifest {output.manifest}" + " {params.patterns}" + + # --- 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 @@ -137,7 +197,16 @@ rule clean_exposure: # spatial neighbours. In-scope consumers keep their edge: they may run in # this DAG, so the clean must be ordered after them. lambda wc: [tile_manifest(t, "tile_vignets") - for t in clean_consumers(wc.exp) if t in READY_SET] + for t in clean_consumers(wc.exp) if t in READY_SET], + # The keepers must be off /scratch before the store goes. Unlike the + # consumer edges above, this edge does not depend on scope: it is the + # same exposure's own rule, so it drags nothing into the DAG that this + # exposure's chain did not already put there. It is conditional only on + # there being a keep list at all — with `persist_exp:` empty, "keep + # nothing" is a coherent instruction and must not become a dependency on + # a rule that would fail for having nothing to copy. + lambda wc: ([prod_exp_manifest(wc.exp, "exp_persist")] + if PERSIST_EXP else []) output: tombstone = f"{EXP_DIR}/cleaned.json" params: @@ -151,3 +220,70 @@ rule clean_exposure: f"python {SCRIPTS}/clean_exposure.py" " --exp-dir $(dirname {output.tombstone}) --exp {wildcards.exp}" " --tombstone {output.tombstone} --consumers '{params.consumers}'" + + +# --- the campaign's star catalogue ------------------------------------------ +# ONE job per campaign: every exposure's every CCD's `validation_psf--.fits`, +# stacked into `/full_starcat-0000000.fits`. That file is the +# rho/tau statistics input and sp_validation reads it at exactly that path, +# doing no merging of its own; the old bash chain built it with +# `combine_runs.bash psf` + a `merge_starcat_runner` pass, and the workflow +# emitted neither. The stacking itself is `MergeStarCatPSFEX` — the same class +# the old runner called, reused rather than restated, so a column added to the +# module is a column added here (merge_star_cat.py argues the reuse and the +# tar-member reading). +# +# THE INPUT IS star_cat_inputs() (Snakefile): every exposure of TILES_READY whose +# PSF products are on the persistent root — the live ones through the exp_persist +# manifest edge `rule all` already requests, the RECLAIMED ones through their TAR, +# which no rule declares and which therefore requires nothing to be built. That +# asymmetry is not a flourish; requesting a reclaimed exposure's manifest +# rebuilds its whole chain from VOS, and ancient() does not prevent it (measured +# — the Snakefile carries the numbers). Nothing new enters the DAG either way. It +# is read through an INPUT FUNCTION rather than at module level so that only a +# parse which actually builds this job pays for the walk. +# +# THE PATHS DO NOT REACH THE SHELL, and that is not a style choice: ~20k manifest +# paths is an order of magnitude over Linux's 128 KiB MAX_ARG_STRLEN for a single +# argv entry, so `{input}` here would be a job that dies on exec at DR6 scale. +# The job is handed the two small files the Snakefile itself started from — the +# tile list and the index — and derives THE SAME SET from them; `params.inputs` +# carries that set's FINGERPRINT, which is the rerun trigger. The equality is +# the point: a job that stacked anything the fingerprint did not see would be +# rows no rerun trigger could notice, which is what a glob over products_dir +# would have given on a root shared with an earlier, larger tile list. +# Byte-stable output otherwise (tmp-then-cmp-then-mv), so a no-op rerun does not +# move its mtime. +# +# NOT A LOCALRULE. exp_persist is local because it is 20k jobs of seconds; this +# is one job that holds a campaign's stars in memory (~800k catalogues at DR6 +# scale). mem_mb is a guess scaled by attempt, not a measurement — the campaigns +# run so far are 127 exposures, three orders of magnitude short of the case this +# sizing is for, and the first DR6-scale run should replace this number with a +# benchmark. +# +# NO JOB AT ALL when `persist_exp:` keeps no validation catalogue, or when every +# exposure in scope is tombstoned: star_cat_targets() (Snakefile) simply does not +# request the output, and the parse says so rather than a node failing later. +rule star_cat_merge: + input: + lambda wc: star_cat_inputs() + output: + star_cat = full_starcat() + params: + products_dir = str(PRODUCTS_DIR), + tile_list = str(config["tile_list"]), + index_db = str(INDEX_DB), + inputs = unit_fingerprint(star_cat_exposures()), + script_hash = MERGE_STAR_HASH + threads: 1 + resources: + mem_mb = lambda wc, attempt: 16000 * attempt, + runtime = 120 + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/merge_star_cat.py" + " --products-dir '{params.products_dir}'" + " --tile-list '{params.tile_list}' --index-db '{params.index_db}'" + " --output {output.star_cat}" + f" --psf-model {PSF_MODEL}" diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index 68d6da9ba..ed5daceb4 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -898,3 +898,63 @@ rule clean_tile: f"python {SCRIPTS}/clean_tile.py" " --tile-dir $(dirname {output.tombstone}) --tile {wildcards.tile}" " --tombstone {output.tombstone}" + + +# --- the campaign's shear catalogue ----------------------------------------- +# ONE job per campaign, the tile-side twin of exposure.smk's star_cat_merge, and +# the same three design calls hold: the input is the list `rule all` already +# requests (every ready tile's final_cat), the paths never reach the shell +# (MAX_ARG_STRLEN), and a fingerprint on `params` is what makes it rerun when a +# tile is appended. The job derives the same set the fingerprint was taken over +# from the tile list and the index rather than globbing products_dir — on a +# products root shared with an earlier, larger tile list a glob would merge tiles +# no rerun trigger ever saw. +# +# THE OUTPUT SCHEMA IS AN INTERFACE, NOT A CHOICE. sp_validation opens this file +# as its `galaxy_cat_path`: one dataset per tile under a named group, the +# columns of workflow/config/cfis/final_cat.param, an `n_tiles` attribute on the +# root. The group is named for the CAMPAIGN, which is the only unit this +# workflow has above the tile. So the rule reuses +# scripts/python/create_final_cat.py's column extraction rather than restating +# it, and writes the file itself — merge_final_cat.py argues that split, the one +# legacy literal in the schema, and the two places where the reference +# implementation had to be pinned down to be reproducible. +# +# THE INPUT IS final_cat, NOT the tile_make_cat manifest, for the same reason +# clean_tile's is: final_cat on the persistent root IS the campaign's +# tile-finished marker (see final_cat() in the Snakefile), and it is the file +# this rule actually reads. +# +# NOT A LOCALRULE, and here the reason is IO rather than memory: the job reads +# every tile's catalogue end to end on every run — ~32-46 MB per tile, so ~2 GB +# for a 64-tile campaign and ~800 GB at DR6's 23k tiles. It rebuilds rather than +# appends because a DAG output must be a function of its input set +# (merge_final_cat.py); incremental update by hand is what +# `create_final_cat.py -s add` remains for. Memory is one tile's catalogue at a +# time plus the hdf5 write buffer, which is why mem_mb is modest where +# star_cat_merge's is not. +rule final_cat_merge: + input: + lambda wc: [final_cat(t) for t in TILES_READY] + output: + merged = final_cat_hdf5() + params: + products_dir = str(PRODUCTS_DIR), + tile_list = str(config["tile_list"]), + index_db = str(INDEX_DB), + param_file = str(CONFIG_DIR / "final_cat.param"), + campaign = CAMPAIGN, + inputs = unit_fingerprint(TILES_READY), + script_hash = MERGE_FINAL_HASH + threads: 1 + resources: + mem_mb = lambda wc, attempt: 8000 * attempt, + runtime = 120 + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/merge_final_cat.py" + " --products-dir '{params.products_dir}'" + " --tile-list '{params.tile_list}' --index-db '{params.index_db}'" + " --output {output.merged}" + " --campaign '{params.campaign}'" + " --param-file '{params.param_file}'" diff --git a/workflow/scripts/build_index.py b/workflow/scripts/build_index.py index 2dd191306..24290b561 100644 --- a/workflow/scripts/build_index.py +++ b/workflow/scripts/build_index.py @@ -152,6 +152,54 @@ def build(tile_ids: list[str], run_dir: Path, db_path: Path, "n_missing": len(missing)} +# --- reading it back, for the campaign-level merges ------------------------- +# The Snakefile loads this index into dicts at parse time and derives the +# campaign's unit sets from them (TILES_READY, and the exposures those tiles +# read). A merge JOB has to derive the same two sets, and cannot be handed them +# on its command line — ~20k paths is an order of magnitude over Linux's 128 KiB +# MAX_ARG_STRLEN for a single argv entry. So it is given the two things the +# Snakefile itself started from, the tile list and this database, and rebuilds +# the sets here. Both halves therefore read the schema through one module rather +# than two hand-written queries that could drift apart. + + +def campaign_tiles(tile_list: Path, db_path: Path) -> list[str]: + """The campaign's ready tiles: declared in the list AND indexed. + + Exactly the Snakefile's TILES_READY, computed the same way from the same two + files — a declared tile with no indexed exposure list cannot have been + computed, so it has no catalogue to merge. + """ + # DEDUPED, order preserved. The tile list is appended to by hand across a + # campaign, so a tile can appear twice; a merge would then try to write that + # tile's dataset twice and die on the second. Deduping here rather than at + # the call sites keeps the answer the same for every reader of the index. + seen, declared = set(), [] + with open(tile_list) as f: + for line in f: + tile = line.strip() + if tile and tile not in seen: + seen.add(tile) + declared.append(tile) + con = sqlite3.connect(db_path, timeout=60) + indexed = {r[0] for r in con.execute("SELECT DISTINCT tile_id FROM tile_exposures")} + con.close() + return [t for t in declared if t in indexed] + + +def campaign_exposures(tile_list: Path, db_path: Path) -> list[str]: + """Every exposure the campaign's ready tiles read, sorted. + + Exactly the set the Snakefile's persist_manifests() builds its manifest + paths from. + """ + tiles = set(campaign_tiles(tile_list, db_path)) + con = sqlite3.connect(db_path, timeout=60) + rows = con.execute("SELECT tile_id, exp_id FROM tile_exposures").fetchall() + con.close() + return sorted({e for t, e in rows if t in tiles}) + + def main() -> None: p = argparse.ArgumentParser(description=__doc__) p.add_argument("--tile-list", required=True, type=Path, diff --git a/workflow/scripts/merge_final_cat.py b/workflow/scripts/merge_final_cat.py new file mode 100644 index 000000000..c12d5fcee --- /dev/null +++ b/workflow/scripts/merge_final_cat.py @@ -0,0 +1,237 @@ +#!/usr/bin/env python3 +"""Collect the campaign's per-tile final catalogues into ONE hdf5 file. + +Run as the shell of the campaign-level ``final_cat_merge`` rule, never by hand. + +WHAT IT PRODUCES, AND FOR WHOM. ``/final_cat_.hdf5``: +one dataset per tile, carrying the columns named by +``workflow/config/cfis/final_cat.param``, plus an ``n_tiles`` attribute on the +file root. sp_validation opens that file as its ``galaxy_cat_path`` +(``sp_validation/catalog.py``), so its SCHEMA is an interface and not a choice — +see ``SPVAL_GROUP`` below for the one legacy literal in it. + +(sp_validation's own ``merge_catalogues`` is a different layer entirely: it +works over already-calibrated ``shape_catalog_comprehensive_*.fits``. It does +not do this merge, and this does not do that one.) + +WHAT IT REUSES, AND WHAT IT DOES NOT. The column extraction is +``create_final_cat.py``'s — ``read_param_file`` for the parameter list, +``read_data`` and ``copy_data`` for pulling those columns out of one catalogue +with their FITS dtypes — so the column grammar keeps exactly one definition. +Its ``process()`` is NOT used and neither is any of its discovery: that function +walks a directory tree the workflow does not have and never will, and it groups +by a unit ShapePipe v2 no longer has. This script walks the workflow's own +products tree instead (``tiles/<2-char prefix>//final_cat-.fits``) and +writes the hdf5 itself. + +WHERE ``create_final_cat.py`` IS FOUND. Beside this workflow, at +``/scripts/python/create_final_cat.py`` — resolved relative to THIS file, +so it follows the launch code snapshot (``bin/sp``) exactly as +``workflow/scripts/*`` does, and a campaign never reads a mid-run edit. It is +loaded by path rather than imported: it is a script, not an installed module, +and the container's ``shapepipe`` install does not carry it. + +IT REBUILDS THE WHOLE FILE, IT DOES NOT APPEND. ``create_final_cat.py``'s own +``process()`` skips tiles already in the file, which is right for a hand-driven +incremental update (``-s add`` / ``-s remove`` are that tool's job). A DAG rule +wants the opposite: the output must be a pure function of the input set, so that +a no-op rerun is byte-stable and a changed set is visibly a different file. +Appending would make the result depend on the order campaigns were run in, and +would silently keep a tile whose catalogue was later rebuilt. The cost is +reading every tile's catalogue on every run of the rule — real work at DR6 scale +(~20k tiles), which is why this is not a localrule. + +BYTE-STABLE ON A NO-OP RERUN: written to a tmp path, compared, moved only if it +differs (the pattern ``persist_exp.py`` and ``clean_exposure.py`` use). Tiles +are visited in sorted ID order so the file is a function of the input set alone. +An unconditional rewrite would move the output's mtime every invocation. + +TWO PLACES WHERE THE REFERENCE IMPLEMENTATION IS NOT DETERMINISTIC, and where +this script therefore pins the behaviour down rather than copying it. Both are +in the DTYPE, and both are invisible when a human runs the tool once by hand: + + * ``copy_data`` allocates ``np.empty`` with the SOURCE catalogue's full + dtype and then fills only the requested columns, so every column NOT in + ``final_cat.param`` reaches the hdf5 file as uninitialised memory — + different bytes on every run, and meaningless data in the file besides. We + hand ``copy_data`` a dtype restricted to the requested columns, so every + field it writes is a field it fills. The file then carries exactly the + ``final_cat.param`` columns, which is what sp_validation reads and what the + parameter file is for. + * ``read_param_file`` returns ``list(set(...))``, whose order varies with the + process's string hash seed. Column ORDER in a structured dtype is part of + the file, so that alone would defeat the byte comparison. We order the + fields by the source catalogue's own column order instead. + +WHICH TILES — AND WHY THE JOB DERIVES THE SET RATHER THAN BEING TOLD IT. The set +is the CAMPAIGN's: every tile both declared in ``tile_list`` and present in the +index, which is exactly the Snakefile's TILES_READY, rebuilt here from the same +two files the Snakefile started from (``--tile-list`` and ``--index-db``, read +through ``build_index.campaign_tiles`` 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 by ``MAX_ARG_STRLEN``; the rule's ``input`` is the DAG edge +and its ``params`` carries a fingerprint of that same list, which is the rerun +trigger. + +THE TWO SETS ARE THE SAME SET, which is the point of deriving it this way rather +than globbing ``/tiles``: a products root shared with an earlier, +larger tile list would hand the job tiles the fingerprint never saw and no rerun +trigger would notice. A tile in the derived set whose catalogue is missing is a +hard error here, not a skip — under the DAG it cannot happen, since every one of +them is a declared input of this job. +""" + +import argparse +import filecmp +import importlib.util +import sys +from pathlib import Path + +import h5py +import numpy as np +from astropy.io import fits + +# Same directory; the rule invokes this file by path, so it is sys.path[0]. +import build_index + +# /scripts/python/create_final_cat.py, from /workflow/scripts/this. +CFC_PATH = (Path(__file__).resolve().parents[2] + / "scripts" / "python" / "create_final_cat.py") + + +def spval_group(campaign: str) -> str: + """The hdf5 group the campaign's per-tile datasets live under. + + ``patches/`` is a LEGACY KEY IN sp_validation's FILE SCHEMA, kept verbatim + only so its reader works unchanged (CosmoStat/sp_validation#340 tracks + removing it); it names nothing in this workflow, which has campaigns and + tiles and no other unit. This is the one place the literal appears — + everything else here says campaign. + """ + return f"patches/{campaign}" + + +def load_create_final_cat(): + """The hdf5 layout's definition, loaded by path (see the module docstring).""" + if not CFC_PATH.exists(): + sys.exit(f"merge_final_cat: {CFC_PATH} is not there — the launch code " + f"snapshot must carry scripts/python/ (see bin/sp).") + spec = importlib.util.spec_from_file_location("create_final_cat", CFC_PATH) + mod = importlib.util.module_from_spec(spec) + spec.loader.exec_module(mod) + return mod + + +def catalogues(products_dir: Path, tile_list: Path, index_db: Path) -> list: + """``(tile ID, path)`` for the campaign's tiles, in ID order. + + Not a glob over the products root: see the module docstring on why the set + is the campaign's and not the filesystem's. + """ + out, missing = [], [] + for tile in sorted(build_index.campaign_tiles(tile_list, index_db)): + path = (products_dir / "tiles" / tile[:2] / tile + / f"final_cat-{tile}.fits") + if path.exists(): + out.append((tile, path)) + else: + missing.append(tile) + if missing: + sys.exit(f"merge_final_cat: {len(missing)} campaign tile(s) have no " + f"final catalogue: {' '.join(missing[:5])}" + f"{' ...' if len(missing) > 5 else ''}") + return out + + +def check_columns(path: Path, hdu: int, wanted: list) -> None: + """Fail loudly, and by name, when a catalogue lacks a requested column.""" + with fits.open(path, memmap=False) as hdu_list: + present = set(hdu_list[hdu].columns.names) + missing = sorted(c for c in wanted if c not in present) + if missing: + sys.exit(f"merge_final_cat: {path} is missing {len(missing)} of the " + f"{len(wanted)} requested column(s): {' '.join(missing)}") + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--products-dir", required=True, type=Path, + help="the persistent root; per-tile catalogues are found " + "beneath it") + p.add_argument("--tile-list", required=True, type=Path, + help="the campaign's tile list (config tile_list)") + p.add_argument("--index-db", required=True, type=Path, + help="the campaign's run index (config outputs.index_db)") + p.add_argument("--output", required=True, type=Path) + p.add_argument("--campaign", required=True, + help="names the campaign's group in the output file") + p.add_argument("--param-file", required=True, type=Path, + help="workflow/config/cfis/final_cat.param — the column list") + p.add_argument("--hdu", type=int, default=1) + args = p.parse_args() + + cfc = load_create_final_cat() + param_list = cfc.read_param_file(str(args.param_file), verbose=False) + if not param_list: + sys.exit(f"merge_final_cat: no columns read from {args.param_file}") + # read_data/copy_data read their knobs out of this dict, exactly as + # create_final_cat.py's own main() builds it. + params = {"hdu_num": args.hdu, "param_list": param_list, "verbose": False} + + tiles = catalogues(args.products_dir, args.tile_list, args.index_db) + if not tiles: + # An empty hdf5 would satisfy every downstream existence check and + # produce an empty shear catalogue. + sys.exit(f"merge_final_cat: no tile in {args.tile_list} is indexed in " + f"{args.index_db}, so there is nothing to merge") + + # tmp-then-cmp-then-mv; the tmp never outlives this process. + args.output.parent.mkdir(parents=True, exist_ok=True) + tmp = args.output.with_name(args.output.name + ".tmp") + try: + tmp.unlink(missing_ok=True) # h5py "a" would reopen a stale one + with h5py.File(tmp, "w") as hdf5_file: + group = hdf5_file.create_group(spval_group(args.campaign)) + columns = None + for tile, path in tiles: + # BEFORE read_data, and not inside it. read_data wraps its + # column selection in a bare `except:` that prints and falls + # through, so a missing column leaves its return values unbound + # and the caller sees UnboundLocalError from the return + # statement — the real name, and every other missing name, never + # reaches the caller at all. Reading the header costs nothing + # next to reading the table. + check_columns(path, args.hdu, params["param_list"]) + extracted, dtype = cfc.read_data(str(path), params) + # Requested columns, in the SOURCE catalogue's order (see the + # module docstring on determinism). Computed from the first + # tile and reused, so a tile whose catalogue is missing a + # column fails loudly on the assignment rather than quietly + # producing a differently-shaped dataset. + if columns is None: + columns = [c for c in dtype.names + if c in set(params["param_list"])] + subset = np.dtype([(c, dtype[c]) for c in columns]) + group.create_dataset( + tile, + data=cfc.copy_data(columns, extracted, subset), + dtype=subset, + ) + # The same attribute create_final_cat.py's print_list() writes, and + # what sp_validation reads to know how many tiles it is holding. + hdf5_file.attrs["n_tiles"] = len(tiles) + + if args.output.exists() and filecmp.cmp(tmp, args.output, shallow=False): + print(f"[merge_final_cat] unchanged: {args.output}") + else: + tmp.replace(args.output) # atomic: same filesystem + print(f"[merge_final_cat] {len(tiles)} tile(s), " + f"{len(param_list)} column(s) -> {args.output} " + f"(group {spval_group(args.campaign)})") + finally: + tmp.unlink(missing_ok=True) + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/merge_star_cat.py b/workflow/scripts/merge_star_cat.py new file mode 100644 index 000000000..25c6e4c22 --- /dev/null +++ b/workflow/scripts/merge_star_cat.py @@ -0,0 +1,265 @@ +#!/usr/bin/env python3 +"""Concatenate the campaign's per-CCD PSF validation catalogues into ONE full_starcat. + +Run as the shell of the campaign-level ``star_cat_merge`` rule, never by hand. + +WHAT IT PRODUCES, AND FOR WHOM. ``/full_starcat-0000000.fits``: +every exposure's every CCD's ``validation_psf--.fits`` row, stacked, +with a ``CCD_NB`` column recording which CCD each row came from. It is the input +to the rho/tau statistics — sp_validation reads exactly this path +(``star_cat_path`` in its ``scripts/calibration/params.py``) and does no merging +of its own. Historically it was ``combine_runs.bash psf`` + a +``merge_starcat_runner`` pass; the workflow emitted neither, so the product set +was short one file. This script is that pass, driven by the DAG instead of by +bash. + +IT DOES NOT REIMPLEMENT THE COLUMN LIST. The stacking, the column names and the +CCD_NB parse all live in ``MergeStarCatPSFEX`` +(``shapepipe.modules.merge_starcat_package.merge_starcat``), which is what the +old runner called. This script only decides WHICH catalogues that class is +handed, and where the result lands. A column added to the module is a column +added here for free — which is the entire reason for the indirection. + +IT READS THE TARS, IT DOES NOT UNPACK THEM, AND IT STREAMS. ``exp_persist`` +packs each exposure's keepers into one uncompressed tar on the persistent root +(``/exp///psf/.tar``) precisely because inodes, +not bytes, bind on /project. Unpacking ~20k tars × ~40 members to merge them +would materialise ~800k files on the filesystem that design exists to protect, +and then delete them. So members are read out of the tars in memory +(``tarfile.extractfile(m).read()`` -> ``io.BytesIO``) and handed to the merge +class as ``[fileobj, member_name]`` pairs — ONE AT A TIME, lazily, through +``TarMembers`` below, because materialising them all first is ~40 GB at DR6 +scale. The member NAME is what the CCD_NB regex parses, which is why the pair +carries it; the class takes the name from the last element of the entry, so a +plain ``[path]`` entry behaves exactly as it always did. + +WHICH EXPOSURES — AND WHY THE JOB DERIVES THE SET RATHER THAN BEING TOLD IT. +The set is the CAMPAIGN's: every exposure read by a tile that is both declared +in ``tile_list`` and present in the index, which is the Snakefile's TILES_READY +walked one edge further. This script rebuilds it from the same two files the +Snakefile started from (``--tile-list`` and ``--index-db``, both small, both on +the persistent root, both read through ``build_index.campaign_exposures`` so +there is one query and not two that can drift), and then takes the exposures +whose ``exp_persist`` manifest is on the persistent root. + +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 by +``MAX_ARG_STRLEN``. Passing them would be a job that dies before it starts. So +the rule's ``input`` is the DAG EDGE — what must exist before this runs — and +the rule's ``params`` carries a FINGERPRINT of that same list, which is what +makes the merge rerun when the set changes. + +THE TWO SETS ARE THE SAME SET, and that equality is the point of deriving it +this way rather than globbing the tree. The rule's input is ``star_cat_inputs()`` +(Snakefile): for each exposure of TILES_READY whose PSF products are on the +persistent root, an edge — the ``exp_persist`` manifest for a live exposure, the +TAR for one whose scratch store reclamation already took (that function argues +the asymmetry, which is about not rebuilding a reclaimed chain from VOS). +Nothing at all for an exposure reclaimed before ``exp_persist`` existed, which +left neither and is unrecoverable short of that rebuild. What this script +selects is the same rule stated from the job's side: same tiles, same index, +manifest present — and by the time the job runs, every exposure with an edge has +one. A glob over ``/exp`` would NOT be the same set: it would +sweep in exposures of an earlier, larger tile list sharing the products root, +stacking rows the fingerprint never saw and no rerun trigger would notice. + +THE MANIFEST, NOT THE TAR, IS WHAT IT READS FIRST: the manifest records what was +actually packed, pattern by pattern, member by member, with sizes. Selecting +members from it means this script never guesses at tar contents, and an exposure +whose keep list did not include the validation catalogues contributes nothing +visibly rather than silently. + +BYTE-STABLE ON A NO-OP RERUN: written to a tmp path, compared, and moved only +if it differs (the pattern ``persist_exp.py`` and ``clean_exposure.py`` use). +An unconditional rewrite would move the output's mtime on every invocation. +Members are visited in sorted (exposure, member) order so the row order is a +function of the input set alone. + +PSFEX ONLY, DELIBERATELY. ``PSF_MODEL`` is ``psfex`` in every campaign the +workflow has run; ``MergeStarCatMCCD`` and ``MergeStarCatSetools`` exist beside +it and take the same constructor, so the hook is the one-line class choice in +``merge_class()`` below — an implementation, not a design, away. +""" + +import argparse +import filecmp +import io +import json +import logging +import shutil +import sys +import tarfile +import tempfile +from fnmatch import fnmatch +from pathlib import Path + +from shapepipe.modules.merge_starcat_package import merge_starcat + +# Same directory; the rule invokes this file by path, so it is sys.path[0]. +import build_index + +# The output name is not ours to choose: sp_validation hardcodes it +# (`star_cat_path = f"{data_dir}/full_starcat-0000000.fits"`), and +# MergeStarCatPSFEX writes exactly this basename into the output dir it is +# given. Kept here as the name this script promises to produce. +OUT_NAME = "full_starcat-0000000.fits" + +# The keep-list pattern whose members this merge consumes. The rule refuses to +# exist unless `persist_exp:` contains a pattern matching this shape (the +# Snakefile does that check at parse time), so by the time we get here the +# members are expected to be present. +MEMBER_PATTERN = "validation_psf-*.fits" + + +def merge_class(psf_model: str): + """The merge class for this PSF model — the one-line MCCD/setools hook. + + Only psfex is exercised: it is what every campaign has run. MCCD reaches the + tars unchanged (it takes its CCD numbers from the data, and it now reports + by the entry's name like the others). SETOOLS would need one more thing — + it passes ``input_file_list[0][0]`` to file_io as a template path, which a + streamed entry is not — so wiring setools to this path is a change to that + class, not a change here. + """ + try: + return {"psfex": merge_starcat.MergeStarCatPSFEX, + "mccd": merge_starcat.MergeStarCatMCCD, + "setools": merge_starcat.MergeStarCatSetools}[psf_model] + except KeyError: + sys.exit(f"merge_star_cat: unknown psf_model {psf_model!r}") + + +def manifests(products_dir: Path, tile_list: Path, index_db: Path) -> list: + """The campaign's exp_persist manifests that are on disk, in exposure order. + + Not a glob over the products root: see the module docstring on why the set + is the campaign's and not the filesystem's. + """ + out = [] + for exp in build_index.campaign_exposures(tile_list, index_db): + path = (products_dir / "exp" / exp[:2] / exp / "manifests" + / "exp_persist.json") + if path.exists(): + out.append(path) + return out + + +def selection(manifest_paths: list, pattern: str) -> tuple: + """``[(tar path, [member names])]`` for the merge, and the empty exposures. + + Reads the manifests only. Every tar is checked for existence HERE, so a + products root missing a file fails before a single row is stacked rather + than an hour in. + """ + chosen, empty = [], [] + for man_path in manifest_paths: + man = json.loads(man_path.read_text()) + wanted = sorted(f["name"] for f in man["files"] + if fnmatch(f["name"], pattern)) + if not wanted: + empty.append(man["unit"]) + continue + tar_path = Path(man["tar"]) + if not tar_path.exists(): + sys.exit(f"merge_star_cat: {man_path} names a tar that is not " + f"there: {tar_path}") + chosen.append((tar_path, wanted)) + return chosen, empty + + +class TarMembers: + """The merge class's input list, materialised ONE TAR AT A TIME. + + ``MergeStarCatPSFEX`` wants something it can take the length of and iterate + once, handing it ``[fileobj, name]`` entries; it never indexes and never + rewinds. So it does not need a list, and a list is the one thing we cannot + afford: reading every member up front is the whole campaign in memory at + once — ~2 MB per exposure, so ~40 GB at DR6's ~20k exposures, against a + rule asking for 16 GB. Read lazily, peak memory is ONE member's bytes plus + the merge class's own accumulators, which are the real and unavoidable term. + + ``__len__`` comes from the manifests, so the class can log the count before + a single tar is opened. + """ + + def __init__(self, chosen): + self._chosen = chosen + + def __len__(self): + return sum(len(names) for _, names in self._chosen) + + def __iter__(self): + for tar_path, names in self._chosen: + with tarfile.open(tar_path) as tf: + for name in names: + member = tf.extractfile(name) + if member is None: + sys.exit(f"merge_star_cat: {tar_path} has no member " + f"{name}, which its manifest lists") + yield [io.BytesIO(member.read()), name] + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--products-dir", required=True, type=Path, + help="the persistent root; exp_persist manifests and tars " + "are found beneath it") + p.add_argument("--tile-list", required=True, type=Path, + help="the campaign's tile list (config tile_list)") + p.add_argument("--index-db", required=True, type=Path, + help="the campaign's run index (config outputs.index_db)") + p.add_argument("--output", required=True, type=Path, + help=f"the merged catalogue; its basename is {OUT_NAME}") + p.add_argument("--psf-model", default="psfex") + p.add_argument("--pattern", default=MEMBER_PATTERN, + help="tar-member glob to merge; default %(default)s") + args = p.parse_args() + + if args.output.name != OUT_NAME: + # The merge class writes OUT_NAME into a directory it is handed; a + # differently-named declared output would silently never be produced. + sys.exit(f"merge_star_cat: --output must be named {OUT_NAME} " + f"(got {args.output.name})") + + log = logging.getLogger("merge_star_cat") + logging.basicConfig(format="[merge_star_cat] %(message)s", + level=logging.INFO, stream=sys.stdout) + + manifest_paths = manifests(args.products_dir, args.tile_list, args.index_db) + chosen, empty = selection(manifest_paths, args.pattern) + file_list = TarMembers(chosen) + if not len(file_list): + # Not a no-op: an empty star catalogue would pass every downstream + # existence check and produce meaningless rho statistics. + sys.exit(f"merge_star_cat: no member matched {args.pattern!r} in any " + f"of {len(manifest_paths)} exp_persist manifest(s) for this " + f"campaign — is '{args.pattern}' in the persist_exp keep list?") + if empty: + log.info(f"{len(empty)} exposure(s) persisted no {args.pattern}: " + f"{', '.join(sorted(empty)[:5])}" + f"{' ...' if len(empty) > 5 else ''}") + + # tmp-then-cmp-then-mv. The merge class chooses its own basename inside the + # directory it is given, so the tmp is a DIRECTORY, not a file path, and it + # never outlives this process — an orphan on /project is an inode nothing + # revisits. + args.output.parent.mkdir(parents=True, exist_ok=True) + tmp_dir = Path(tempfile.mkdtemp(dir=args.output.parent, + prefix=".star_cat_merge.")) + try: + merge_class(args.psf_model)(file_list, str(tmp_dir), log).process() + tmp = tmp_dir / OUT_NAME + if not tmp.exists(): + sys.exit(f"merge_star_cat: the merge wrote no {OUT_NAME}") + if args.output.exists() and filecmp.cmp(tmp, args.output, shallow=False): + log.info(f"unchanged: {args.output}") + else: + tmp.replace(args.output) # atomic: same filesystem + log.info(f"{len(file_list)} catalogue(s) from {len(chosen)} " + f"exposure(s) -> {args.output}") + finally: + shutil.rmtree(tmp_dir, ignore_errors=True) + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/persist_exp.py b/workflow/scripts/persist_exp.py new file mode 100644 index 000000000..40e8c9354 --- /dev/null +++ b/workflow/scripts/persist_exp.py @@ -0,0 +1,205 @@ +#!/usr/bin/env python3 +"""Pack ONE exposure's keepable PSF products into a tar off scratch, and record what went. + +Run as the shell of the in-DAG ``exp_persist`` rule, never by hand. + +WHY A COPY AND NOT AN EXEMPTION FROM CLEANUP. The obvious alternative — teach +``clean_exposure`` to spare these files — does not work, because reclamation is +not what threatens them. The exposure store lives on ``run_dir``, which is +/scratch: a 60-day purge takes everything there whether or not this workflow +ever cleaned it. ``products_dir`` is /project, backed up and not purged. So the +only way a per-exposure product outlives its campaign is to LEAVE THE +FILESYSTEM, and that is a copy. Reclamation ordering then falls out for free: +``clean_exposure`` takes this rule's manifest as an input, so the store is never +deleted before its keepers have been written elsewhere. + +WHY A SEPARATE RULE AND NOT A ``cp`` APPENDED TO ``exp_psf``. The list of what +to keep is a decision that will be revisited — rho statistics want one file +today, a residual study may want three tomorrow — and ``exp_psf`` is four hours +per exposure. The list rides on this rule's ``params``, so editing it makes +snakemake rerun THIS rule (seconds of cp) and leaves the PSF chain alone. Folded +into ``exp_psf``, the same edit would re-derive every PSF model in the campaign. + +WHAT IT SEARCHES. ``/output/run_sp_exp_SxSePsfPi/*/output/`` — the four +module output dirs of the PSF config (sextractor, setools, psfex, psfex_interp) +— RECURSIVELY. The recursion is not laziness: setools does not write flat, it +writes into ``mask/``, ``rand_split/``, ``new_cat/``, ``plot/`` and ``stat/`` +beneath its own output dir, so a caller who wrote ``star_split_ratio_80-*.fits`` +meaning "the training star sample" would match nothing under a non-recursive +glob. Patterns are therefore plain FILE names and the layout is ours to know, +not the config author's. + +ZERO MATCHES FOR ONE PATTERN IS A WARNING, NOT A FAILURE. setools rejects sparse +CCDs (~0.2% attrition, tolerated by exp_psf's own count floor), so per-CCD +counts are not fixed, and a pattern naming an optional diagnostic may legitimately +find nothing. ZERO FILES IN TOTAL IS A FAILURE: it means the store was not what +we think it is, and writing a green manifest over that would let +``clean_exposure`` delete an exposure whose products were never saved. + +The manifest lists every member (name, pattern, source path, bytes), so a reader +knows what the tar holds without opening it. + +ONE UNCOMPRESSED TAR PER EXPOSURE, ``/.tar``, NOT LOOSE COPIES. +Inodes, not bytes, are what bind on /project: the group quota is ~1 M files, +and loose per-CCD copies are ~200 per exposure with all candidates on — ~25k for +a 64-tile campaign, ~2 M at DR6 scale, against ~7 GB of bytes. A tar collapses +that to one inode per exposure and costs nothing to read: FITS members go +``tarfile.open(t).extractfile(m).read()`` -> ``fits.open(io.BytesIO(...))``, +which is why a tar rather than a multi-HDU FITS bundle (the keep list mixes +FITS, ``.psf`` and ``.txt``; a FITS container could not hold the last two). +Uncompressed because FITS barely compresses and a plain tar is seekable. + +Members are FLAT — file name only, no module subtree — because the module a +file came from is already in its name and the consumer globs member names. A +name collision between two modules is therefore a hard error rather than a +silent overwrite; nothing in the current config can produce one, and if a +future one can we want to hear about it. + +The tar is written DETERMINISTICALLY (ownership zeroed, members in sorted +order, source mtimes kept), tmp-then-``cmp``-then-``mv``: a rerun over an +unchanged store produces a byte-identical tar and leaves the existing one's +mtime alone. + +The manifest is the rule's ONLY declared output, and it lives on the persistent +root beside the tar (``/exp///manifests/``, beside the tar's ``psf/``), NOT in +the exposure's scratch ``manifests/`` dir which ``clean_exposure`` deletes +wholesale. It is deliberately NOT a ``directory()`` output: what was copied, and +how big each file was, is provenance we want written down, and a directory +output attests only that some directory exists. + +It carries no timestamp and is written tmp-then-``cmp``-then-``mv`` (the pattern +``clean_exposure`` uses), so a rerun that packs the same files leaves the mtime +alone — mtime is a rerun trigger, and an unconditional rewrite would make every +downstream ``clean_exposure`` look out of date once per invocation. +""" + +import argparse +import filecmp +import json +import sys +import tarfile +from pathlib import Path + +# The PSF chain's run dir (RUN_NAME in config_exp_psfex.ini). Hardcoded rather +# than passed: this rule persists the PSF stage's products and nothing else, and +# a knob here would be a knob for "persist some other stage", which is a +# different rule. +RUN_NAME = "run_sp_exp_SxSePsfPi" + + +def collect(exp_dir: Path, patterns: list) -> tuple: + """Matched files per pattern, in a stable order, plus the empty patterns.""" + root = exp_dir / "output" / RUN_NAME + found, empty = {}, [] + for pat in patterns: + # One glob per module output dir, recursive beneath it (see the module + # docstring on setools' subdirectories). sorted() over the union keeps + # the manifest byte-stable across filesystem readdir order. + hits = sorted({p for mod in sorted(root.glob("*/output")) + for p in mod.rglob(pat) if p.is_file()}) + if hits: + found[pat] = hits + else: + empty.append(pat) + return found, empty + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--exp-dir", required=True, type=Path, + help="the exposure's scratch store") + p.add_argument("--exp", required=True) + p.add_argument("--dest", required=True, type=Path, + help="/exp///psf; the tar is " + "/.tar") + p.add_argument("--manifest", required=True, type=Path) + p.add_argument("--pattern", action="append", default=[], + help="repeatable; a plain file-name glob") + args = p.parse_args() + + if not args.pattern: + sys.exit("persist_exp: no --pattern given (config persist_exp is empty)") + + found, empty = collect(args.exp_dir, args.pattern) + if not found: + sys.exit(f"persist_exp: {args.exp}: no file matched any of " + f"{args.pattern} under {args.exp_dir}/output/{RUN_NAME}") + + args.dest.mkdir(parents=True, exist_ok=True) + tar_path = args.dest / f"{args.exp}.tar" + # A file matched by TWO patterns is one file, not a collision. Keep lists + # overlap on purpose — `validation_psf-*.fits` alongside `*.fits` is a + # perfectly ordinary way to say "the validation catalogues, and everything + # else FITS while we are here" — and treating the second match as a name + # clash failed every exposure in the campaign. What must still be fatal is + # two DIFFERENT paths landing on one flat member name, which would silently + # overwrite; that is a same-name/different-source test, and the first + # pattern to match a file is the one recorded for it. + seen, files = {}, [] + for pat, hits in found.items(): + for src in hits: + if src.name in seen: + if seen[src.name][0] == src: + continue # same file, a second matching pattern + sys.exit(f"persist_exp: {args.exp}: two source files are both " + f"named {src.name} ({seen[src.name][0]} and {src}); tar " + f"members are flat, so this would silently overwrite") + seen[src.name] = (src, pat) + files.append({"name": src.name, "pattern": pat, + "src": str(src), "bytes": src.stat().st_size}) + files.sort(key=lambda f: f["name"]) + + def anonymous(ti: tarfile.TarInfo) -> tarfile.TarInfo: + # Ownership is the one thing that would differ between two writes of + # the same files from different accounts/nodes; drop it. mtime stays: + # it is the product's, and it is stable while the store is. + ti.uid = ti.gid = 0 + ti.uname = ti.gname = "" + return ti + + # tmp-then-cmp-then-mv, and the tmp NEVER outlives a failure: an orphaned + # .tmp on /project is an inode nothing revisits — the leak this whole tar + # design exists to avoid, one per failed attempt at DR6 scale. + tmp = tar_path.with_name(tar_path.name + ".tmp") + try: + with tarfile.open(tmp, "w", format=tarfile.PAX_FORMAT) as tf: + for f in files: + tf.add(seen[f["name"]][0], arcname=f["name"], filter=anonymous) + if tar_path.exists() and filecmp.cmp(tmp, tar_path, shallow=False): + tmp.unlink() # unchanged: leave the mtime alone + else: + tmp.replace(tar_path) # atomic: no half-written archive + finally: + tmp.unlink(missing_ok=True) + + body = { + "stage": "exp_persist", "level": "exp", "unit": args.exp, + "status": "complete", + "tar": str(tar_path), + "patterns": list(args.pattern), + # The warning the docstring argues for: named patterns that matched + # nothing. Present as a key even when empty, so a reader never has to + # wonder whether an old manifest predates the field. + "patterns_unmatched": empty, + "n_files": len(files), + "bytes": sum(f["bytes"] for f in files), + "files": files, + } + 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") + if args.manifest.exists() and filecmp.cmp(tmp, args.manifest, shallow=False): + tmp.unlink() # unchanged: leave the mtime alone + else: + tmp.replace(args.manifest) + finally: + tmp.unlink(missing_ok=True) + + warn = f" ({len(empty)} pattern(s) matched nothing: {empty})" if empty else "" + print(f"[persist_exp] {args.exp}: {len(files)} file(s), " + f"{body['bytes'] / 1e6:.1f} MB -> {tar_path}{warn}") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/run_report.py b/workflow/scripts/run_report.py index 24c07670b..11cdcffe4 100644 --- a/workflow/scripts/run_report.py +++ b/workflow/scripts/run_report.py @@ -61,6 +61,13 @@ "tile_ngmix", "tile_merge_cats", "tile_make_cat"] EXP_STAGES = ["exp_get_images", "exp_split", "exp_psf"] +# exp_persist is DELIBERATELY NOT in that list. This report disk-scans the +# scratch run_dir, and exp_persist's manifest is the one exposure manifest that +# lives on products_dir instead — that placement is what makes it survive +# clean_exposure. Listed here it would read as "not run" for every exposure in +# the campaign. Reporting on the persisted products means scanning the second +# root, which is a report this one does not yet do. + # The manifests clean_tile leaves on disk (workflow/scripts/clean_tile.py names # the mechanism that owns each). Their presence is therefore NOT evidence that a # tile's chain was rebuilt, which absorb_tombstones needs to know