diff --git a/.gitignore b/.gitignore index db20f3086..4c4b6d0b0 100644 --- a/.gitignore +++ b/.gitignore @@ -139,3 +139,4 @@ code # .felt here is a machine-local symlink into it. Never track it in this repo. /.felt/ /.felt +.snakemake/ diff --git a/example/cfis/config_tile_PiViVi_canfar_sx.ini b/example/cfis/config_tile_PiViVi_canfar_sx.ini index ca7efd9ec..991153143 100644 --- a/example/cfis/config_tile_PiViVi_canfar_sx.ini +++ b/example/cfis/config_tile_PiViVi_canfar_sx.ini @@ -79,7 +79,7 @@ POSITION_PARAMS = XWIN_WORLD,YWIN_WORLD GET_SHAPES = True # Number of stars threshold -STAR_THRESH = 20 +STAR_THRESH = 22 # chi^2 threshold CHI2_THRESH = 2 diff --git a/example/cfis/config_tile_PiViVi_canfar_uc.ini b/example/cfis/config_tile_PiViVi_canfar_uc.ini index d59af1c4a..b33361b6d 100644 --- a/example/cfis/config_tile_PiViVi_canfar_uc.ini +++ b/example/cfis/config_tile_PiViVi_canfar_uc.ini @@ -77,7 +77,7 @@ POSITION_PARAMS = XWIN_WORLD,YWIN_WORLD GET_SHAPES = True # Number of stars threshold -STAR_THRESH = 20 +STAR_THRESH = 22 # chi^2 threshold CHI2_THRESH = 2 diff --git a/profiles/nibi/config.yaml b/profiles/nibi/config.yaml new file mode 100644 index 000000000..7b75b458f --- /dev/null +++ b/profiles/nibi/config.yaml @@ -0,0 +1,97 @@ +# Snakemake profile for the Nibi cluster (Digital Research Alliance). +# +# SLURM-EXECUTOR mode (PRD #848 D-profile): one SLURM job per rule instance, +# each carrying its rule's own attempt-scaled resources (cpus_per_task, +# mem_mb, runtime). Snakemake feeds the queue as jobs finish, so DR6's ~170k +# total jobs never queue at once, and multi-node scaling is inherent — this +# supersedes the earlier one-allocation/local-scheduler mode. +# +# Launch via `workflow/bin/sp` (loads apptainer/1.4.5, uses the /project venv). +# software-deployment-method wraps every job's shell in `apptainer exec` — the +# user never types apptainer. WHICH image is resolved in the Snakefile through +# workflow/scripts/container.py: this user's writable sandbox if they built one, +# else their cached SIF (`sp container pull`), else the `container:` path in +# workflow/config.yaml — the shared /project .sif, and the default whenever the +# cache is empty. `sp container status` prints which layer is live. `sp container +# pull` needs the network: login node or salloc, never a batch job. + +executor: slurm + +# Per-user submit cap, queried 2026-07-30 on nibi: +# sacctmgr show assoc user=cdaley format=Account,MaxSubmitPU -P +# -> def-mjhudson_cpu 1000, def-mjhudson_gpu 1000 (MaxSubmitPU; no site-wide +# MaxSubmitJobs in `scontrol show config`, so the association limit governs) +# Set to ~80% of that (1000) so this workflow never starves other submissions +# under the same account. Re-query if the association limits change. +jobs: 800 + +default-resources: + mem_mb: 2000 + runtime: 120 # minutes + slurm_account: def-mjhudson_cpu + +software-deployment-method: [apptainer] +# Explicit container environment (finding 2/3): the apptainer SDM otherwise drops +# the proven-recipe env and hardcodes --home , hiding ~/.ssl/cadcproxy.pem. +# --cleanenv strict host-env isolation (APPTAINERENV_*/SINGULARITYENV_* survive) +# OMP_NUM_THREADS=1 caps OpenBLAS fork-explosion (verified: pool 32->1) +# MALLOC_ARENA_MAX=2 bounded allocator (im_sims/nibi lesson) +# --home /home/cdaley wins over the SDM's --home ; restores cadcproxy.pem for vos/vcp +# PYTHONPATH SETTLED CALL 3 — pins THIS branch's src/, which is +# develop@97e16d50 plus the four commits that genuinely need +# module code: ngmix chunk fields + position-seeded RNG, +# merge_sep_cats chunk paths, the vizier star-cat helpers, +# and the cherry-picked #873 (seeded setools split). +# NOT shapepipe-prod (drifted to a PR branch mid-run, and the +# live p3-batch1 job reads it) and not the sif default +# (frozen pre-#843, and pre-#873). Production later rebuilds +# the sif at the validated commit and DROPS this --env line. +apptainer-args: "--cleanenv --env OMP_NUM_THREADS=1 --env MALLOC_ARENA_MAX=2 --env PYTHONPATH=/project/def-mjhudson/cdaley/shapepipe-snakemake/src --home /home/cdaley --bind /project --bind /scratch" + +latency-wait: 60 # NFS: wait for outputs to appear after a job +keep-going: true # a failed job poisons only its cone; siblings run on +rerun-incomplete: true # re-do jobs left incomplete by an unclean death +show-failed-logs: true +printshellcmds: true + +# No `keep-incomplete` here. The declared output IS the manifest, so letting +# snakemake delete a failed job's output is exactly the semantics this workflow +# wants: a `.json` on disk means the stage succeeded, and a resume cannot +# schedule downstream work on top of a failure. The post-mortem evidence is the +# rule's `log:` (`/logs/.json`), which completeness.py writes on +# every run and which snakemake preserves through the failure that produced it — +# which is also what `show-failed-logs` above surfaces on the console. + +# rerun-triggers: the v9 default MINUS `input`. params/code/mtime fixes still +# propagate — completeness.py writes the manifest only on change, so mtimes move +# only when reality moves. +# +# `input` is dropped because reclamation needs the tile->exposure edge to be +# CONDITIONAL: a tile whose final_cat is on disk declares no exposure inputs, so +# a neighbour rebuilding a shared exposure cannot drag it along (see tile.smk). +# With the `input` trigger on, that same conditional reads as "set of input files +# has changed" and reruns every finished tile — against an exposure store that +# reclamation has already deleted. Measured on fixture t4: 70 jobs with the +# trigger, 28 without, for one damaged tile in a four-tile chain. +# +# Nothing this workflow relied on is lost. The two genuinely data-derived input +# sets are covered another way: a changed exposure list arrives through the +# (non-ancient) find_exposures manifest's mtime, and the ngmix chunk count rides +# in params. And clean scheduling is structurally gated to bin/sp (SP_PHASE=compute +# + this profile), so a bare snakemake invocation cannot quietly recombine +# reclamation with the `input` trigger; a runtime assertion is impossible +# (the trigger set is unreadable at parse time) — the gate is the launcher. +rerun-triggers: [mtime, params, code, software-env] + +# NO set-threads / set-resources here, deliberately. Profile overrides REPLACE a +# rule's own values (verified snakemake 9.23), which would kill the +# attempt-scaled `mem_mb = lambda wc, attempt: ...` OOM retries and the tuned +# ngmix thread count. The RULES own threads and resources; this profile only +# sets defaults for rules that state nothing (default-resources above). +# +# `group:` fusion of the short rules (PRD D-profile) lives in the rule files, not +# here: the labels are `tile_prep` (prepare.smk), `exp_short` (exp_split + +# exp_mask) and `tile_finish` (tile_merge_cats + tile_make_cat), each documented +# in its file's docstring. Grouping composes resources per toposort level (max +# mem/threads, summed runtime for a linear chain) and leaves the attempt scaling +# above intact, which is the other reason set-resources must stay out of here. diff --git a/scripts/python/create_star_cat.py b/scripts/python/create_star_cat.py index 0e3586f88..404cee4ac 100755 --- a/scripts/python/create_star_cat.py +++ b/scripts/python/create_star_cat.py @@ -118,6 +118,29 @@ def query_vizier(ra, dec, radius_arcmin): return _query_vizier(ra, dec, radius_arcmin, CDS_CAT_ID) +def _write_atomic(table, output_dir, img_number, output_name): + """Write ``table`` to ``output_name`` atomically. + + The catalogue is a run-independent CACHE, and the caller's only test for a + cache hit is ``os.path.isfile``. An in-place ``table.write`` that is killed + part-way (job timeout, OOM, node failure) therefore leaves a truncated FITS + that every later run trusts forever — and ``test -s`` passes on partial + bytes. Writing to a temp and renaming makes the visible file all-or-nothing: + ``os.replace`` is atomic within a directory. + + The temp keeps a ``.fits`` suffix, because astropy picks the writer from the + extension. It is dot-prefixed and PID-tagged so it stays out of the + ``star_cat*`` globs the rules use, and two concurrent writers cannot collide. + """ + tmp = f"{output_dir}/.tmp-{os.getpid()}-star_cat{img_number}.fits" + try: + table.write(tmp, overwrite=True) + os.replace(tmp, output_name) + finally: + if os.path.exists(tmp): + os.remove(tmp) + + def main(input_dir, output_dir, kind): file_list = os.listdir(input_dir) @@ -140,7 +163,7 @@ def main(input_dir, output_dir, kind): f"Focal plane center: ra={ra:.4f}, dec={dec:.4f}, radius={radius:.2f} arcmin" ) table = query_vizier(ra, dec, radius) - table.write(output_name, overwrite=True) + _write_atomic(table, output_dir, img_number, output_name) else: h = fits.getheader(fpath, 0) @@ -156,7 +179,7 @@ def main(input_dir, output_dir, kind): continue table = query_vizier(ra, dec, radius) - table.write(output_name, overwrite=True) + _write_atomic(table, output_dir, img_number, output_name) return 0 diff --git a/src/shapepipe/modules/merge_sep_cats_package/merge_sep_cats.py b/src/shapepipe/modules/merge_sep_cats_package/merge_sep_cats.py index 5ecdb6fb0..b0be88944 100644 --- a/src/shapepipe/modules/merge_sep_cats_package/merge_sep_cats.py +++ b/src/shapepipe/modules/merge_sep_cats_package/merge_sep_cats.py @@ -7,7 +7,6 @@ """ import os -import re import warnings import numpy as np @@ -16,6 +15,48 @@ from shapepipe.pipeline import file_io +def chunk_path(input_file, n): + """Chunk Path. + + Derive chunk ``n``'s input path from chunk 1's. + + The separate catalogues live in ShapePipe run directories whose names differ + only in the chunk number (``run_sp_tile_ngmix_Ng1u`` -> + ``run_sp_tile_ngmix_Ng2u``). The substitution is confined to that + run-directory component: replacing the first "1" found anywhere in the path + breaks for absolute paths whose parent directories carry digits, e.g. a + sharded store ``.../tiles/21/210.282/output/run_..._Ng1u/...``. + + Parameters + ---------- + input_file : str + Path to chunk 1's catalogue + n : int + Chunk number + + Returns + ------- + str + Path to chunk ``n``'s catalogue + + Raises + ------ + ValueError + If no run-directory component of the path carries a chunk number + + """ + parts = input_file.split(os.sep) + for idx in reversed(range(len(parts))): + if parts[idx].startswith("run_") and "1" in parts[idx]: + parts[idx] = parts[idx].replace("1", str(n), 1) + return os.sep.join(parts) + + raise ValueError( + f"Cannot derive chunk {n}'s path from '{input_file}': no 'run_*' " + + "directory component contains a chunk number '1'" + ) + + class MergeSep(object): """Merge Sep. @@ -79,8 +120,7 @@ def process(self): input_path_n = [] input_path_n.append(input_file) for n in range(2, self._n_split_max + 1): - res = re.sub("1", str(n), input_file, 1) - input_path_n.append(res) + input_path_n.append(chunk_path(input_file, n)) # Open first catalogue, read number of extensions and columns cat0 = file_io.FITSCatalogue(input_file, SEx_catalogue=True) diff --git a/src/shapepipe/modules/merge_sep_cats_runner.py b/src/shapepipe/modules/merge_sep_cats_runner.py index 927c79b76..6e99cf3ec 100644 --- a/src/shapepipe/modules/merge_sep_cats_runner.py +++ b/src/shapepipe/modules/merge_sep_cats_runner.py @@ -28,7 +28,7 @@ def merge_sep_cats_runner( ): """Define The Merge SEP Catalogues Runner.""" # Get config entries - n_split_max = config.getint(module_config_sec, "N_SPLIT_MAX") + n_split_max = int(config.getexpanded(module_config_sec, "N_SPLIT_MAX")) file_pattern = config.getlist(module_config_sec, "FILE_PATTERN") file_ext = config.getlist(module_config_sec, "FILE_EXT") diff --git a/src/shapepipe/modules/ngmix_package/__init__.py b/src/shapepipe/modules/ngmix_package/__init__.py index 230606bba..fca6ee0bd 100644 --- a/src/shapepipe/modules/ngmix_package/__init__.py +++ b/src/shapepipe/modules/ngmix_package/__init__.py @@ -40,15 +40,27 @@ (no batch saving) ID_OBJ_MIN : int ID of first galaxy object to be processed; not used if set to ``-1`` - (default) + (default). Environment variables are expanded, so an orchestrator can + set the object range per chunk, for example + ``ID_OBJ_MIN = $SP_NGMIX_ID_OBJ_MIN``. ID_OBJ_MAX : int ID of last galaxy object to be processed; not used if set to ``-1`` - (default) + (default). Environment variables are expanded, as for ``ID_OBJ_MIN``. BKG_RMS_VIGNET_PATH : str, optional Path to a ``background_rms_vignet*.sqlite`` file produced by ``vignetmaker_runner``. The string may contain ``{file_number_string}``, which is replaced by the current tile ID. +Random number generation +======================== + +Each object gets its own random number stream, seeded from its sky position +and CCD (``position_seed``, ngmix#796). The output is therefore identical +whether a tile is processed in one go or split into object chunks with +``ID_OBJ_MIN``/``ID_OBJ_MAX``. The older tile-seeded mode is retired; the +config option ``SEED_FROM_POSITION`` is obsolete. Setting it to ``False`` +raises an error, so a stale config cannot silently change the RNG. + """ __all__ = ["ngmix"] diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index d4f64a74e..fc9779c6a 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -335,13 +335,12 @@ class Ngmix(object): (robust for galaxies); ``"wcs"`` uses the catalog sky position projected through the WCS (better for stars, whose HSM moments are noisy). See :func:`make_ngmix_observation`. - seed_from_position : bool, optional - If ``True``, replace the tile-level RNG with a per-object RNG seeded - from the object's sky position (:func:`position_seed`) inside the - object loop, so metacal's ``fixnoise`` counter-noise and the fit - guesses cancel across Pujol image-simulation branches (ngmix#796). The - default ``False`` leaves the production path byte-identical. See - :func:`position_seed` for the physics and the seed construction. + Notes + ----- + The RNG is always per object and seeded from that object's sky position + (:func:`position_seed`). Results therefore do not depend on how the tile is + split into object chunks, and metacal's ``fixnoise`` counter-noise and the + fit guesses cancel across Pujol image-simulation branches (ngmix#796). Raises ------ @@ -364,7 +363,6 @@ def __init__( id_obj_max=-1, bkg_sub=True, centroid_source="hsm", - seed_from_position=False, metacal_psf="fitgauss", ): @@ -418,20 +416,14 @@ def __init__( self._id_obj_max = id_obj_max self._bkg_sub = bkg_sub self._centroid_source = centroid_source - self._seed_from_position = seed_from_position self._metacal_psf = metacal_psf self._w_log = w_log - # Initiatlise random generator - seed = int(''.join(re.findall(r'\d+', self._file_number_string))) - self._rng = np.random.RandomState(seed) - self._w_log.info(f'Random generator initialisation seed = {seed}') - if self._seed_from_position: - self._w_log.info( - 'SEED_FROM_POSITION on: per-object RNG seeded from sky position' - ' for Pujol noise cancellation (image sims, ngmix#796)' - ) + self._w_log.info( + 'Per-object RNG seeded from sky position (ngmix#796): results are' + ' invariant to how the tile is split into object chunks' + ) @classmethod def MegaCamFlip(self, vign, ccd_nb): @@ -461,18 +453,6 @@ def MegaCamFlip(self, vign, ccd_nb): # swap y axis so origin is on bottom-left return vign - def get_prior(self, T_range=None, F_range=None): - """Get Prior. - - Returns - ------- - ngmix.joint_prior.PriorSimpleSep - """ - return get_prior( - self._pixel_scale, self._rng, - T_range=T_range, F_range=F_range, - ) - def compile_results(self, results): """Compile Results. @@ -809,7 +789,6 @@ def process(self): vignet_cat = self._vignet_cat final_res = [] - prior = self.get_prior() count = 0 n_empty_cat = 0 @@ -843,24 +822,20 @@ def process(self): n_no_epoch += 1 continue - # Position-seeded per-object RNG for Pujol noise cancellation in - # image sims (ngmix#796): the same object gets the same fixnoise - # counter-noise and fit guesses in every shear branch, so both - # cancel in the branch difference. The prior is rebuilt from the - # same per-object RNG because the guesser draws its initial guess - # via prior.sample() (ngmix guessers.py), which consumes the RNG the - # prior was CONSTRUCTED with — so a per-object rng alone would leave - # the guess drawing from the shared tile stream and break - # cancellation. Off in production, where the single tile-level - # self._rng and the tile-level prior carry the whole loop. - if self._seed_from_position: - obj_rng = np.random.RandomState( - position_seed(stamp.ra[0], stamp.dec[0], stamp.ccd) - ) - obj_prior = get_prior(self._pixel_scale, obj_rng) - else: - obj_rng = self._rng - obj_prior = prior + # Position-seeded per-object RNG (ngmix#796). Each object draws from + # a stream fixed by its own (ra, dec, ccd), so the result is + # independent of which chunk the object lands in and of detection + # order, and the same object gets the same fixnoise counter-noise + # and fit guesses in every Pujol shear branch, so both cancel in the + # branch difference. The prior is rebuilt from the same per-object + # RNG because the guesser draws its initial guess via prior.sample() + # (ngmix guessers.py), which consumes the RNG the prior was + # CONSTRUCTED with — a per-object rng alone would leave the guess + # drawing from a shared stream and break both properties. + obj_rng = np.random.RandomState( + position_seed(stamp.ra[0], stamp.dec[0], stamp.ccd) + ) + obj_prior = get_prior(self._pixel_scale, obj_rng) try: flux_guess = ( diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 3d98356c4..fcb50be9f 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -86,9 +86,12 @@ def ngmix_runner( # No batch saving save_batch = -1 - # First and last galaxy ID to process - id_obj_min = config.getint(module_config_sec, "ID_OBJ_MIN") - id_obj_max = config.getint(module_config_sec, "ID_OBJ_MAX") + # First and last galaxy ID to process. Read via ``getexpanded`` so an + # orchestrator can drive the chunk bounds from environment variables + # (``$SP_NGMIX_ID_OBJ_MIN`` and friends); ``getexpanded`` is the only + # accessor in ShapePipe's config that expands ``$VAR``. + id_obj_min = int(config.getexpanded(module_config_sec, "ID_OBJ_MIN")) + id_obj_max = int(config.getexpanded(module_config_sec, "ID_OBJ_MAX")) # Centroid source for the galaxy Jacobian origin: "wcs" (default -- the # catalog sky position projected through the WCS, trusting the astrometry) @@ -99,16 +102,20 @@ def ngmix_runner( else: centroid_source = "wcs" - # Seed the per-object RNG from sky position instead of per tile, so - # metacal's fixnoise counter-noise (and the fit guesses) cancel across - # Pujol image-simulation shear branches (ngmix#796). Default False leaves - # the production path byte-identical. + # Position-seeded RNG is the only mode: every object's RNG comes from its + # own (ra, dec, ccd), so results do not depend on how the tile is split + # into chunks, and metacal's fixnoise counter-noise cancels across Pujol + # image-simulation shear branches (ngmix#796). The retired tile-seed mode + # had neither property. Old configs that disable it must fail loudly. if config.has_option(module_config_sec, "SEED_FROM_POSITION"): - seed_from_position = config.getboolean( - module_config_sec, "SEED_FROM_POSITION" - ) - else: - seed_from_position = False + if not config.getboolean(module_config_sec, "SEED_FROM_POSITION"): + raise ValueError( + "SEED_FROM_POSITION = False is no longer supported: the" + " tile-seeded RNG mode has been retired because it makes" + " results depend on the object chunking. Remove the" + " SEED_FROM_POSITION entry from the ngmix config section" + " (position-seeded RNG is now the only mode)." + ) # Check PSF vignets first: if all are empty dicts {}, the exposures for this # tile are absent from the PSF dictionary and no shape measurement is possible. @@ -161,7 +168,6 @@ def ngmix_runner( id_obj_max=id_obj_max, bkg_sub=bkg_sub, centroid_source=centroid_source, - seed_from_position=seed_from_position, metacal_psf=metacal_psf, ) diff --git a/src/shapepipe/modules/setools_package/setools.py b/src/shapepipe/modules/setools_package/setools.py index c430b52b1..33d8efa8b 100644 --- a/src/shapepipe/modules/setools_package/setools.py +++ b/src/shapepipe/modules/setools_package/setools.py @@ -658,13 +658,18 @@ def _make_rand_split(self): cat_size = len(np.where(mask)[0]) n_keep = int(np.ceil(cat_size * ratio)) - mask_ratio = [] - mask_left = list(range(0, cat_size)) - while len(mask_ratio) != n_keep: - idx = np.random.randint(0, len(mask_left)) - mask_ratio.append(mask_left.pop(idx)) - mask_ratio = np.array(mask_ratio) - mask_left = np.array(mask_left) + # Deterministic split, seeded from the unit's file number: the + # train/validation assignment is a pure function of the input + # catalogue, so the PSF star sample (and everything downstream + # of the PSF model) is reproducible run-to-run. An unseeded + # np.random here made the shear catalogue non-reproducible + # upstream of ngmix's own position seeding. + seed = int( + re.sub(r"\D", "", self._file_number_string) or 0 + ) % (2 ** 32) + perm = np.random.RandomState(seed).permutation(cat_size) + mask_ratio = perm[:n_keep] + mask_left = np.sort(perm[n_keep:]) self.rand_split[key]["mask"] = mask self.rand_split[key][f"ratio_{int(ratio * 100)}"] = mask_ratio self.rand_split[key][f"ratio_{100 - int(ratio * 100)}"] = mask_left diff --git a/src/shapepipe/run.py b/src/shapepipe/run.py index 962953ede..fee0dad3e 100644 --- a/src/shapepipe/run.py +++ b/src/shapepipe/run.py @@ -84,7 +84,7 @@ def _set_run_name(self): Set the name of the current pipeline run. """ - self._run_name = self.config.get("DEFAULT", "RUN_NAME") + self._run_name = self.config.getexpanded("DEFAULT", "RUN_NAME") if self.config.getboolean("DEFAULT", "RUN_DATETIME"): self._run_name += datetime.now().strftime("_%Y-%m-%d_%H-%M-%S") diff --git a/src/shapepipe/utilities/vizier.py b/src/shapepipe/utilities/vizier.py index 1f4fd3472..0bf077a78 100644 --- a/src/shapepipe/utilities/vizier.py +++ b/src/shapepipe/utilities/vizier.py @@ -64,8 +64,16 @@ def query_vizier(ra, dec, radius_arcmin, cat_id): v = Vizier( row_limit=-1, timeout=timeout, vizier_server=server ) + # cache=False: astroquery otherwise pickles every HTTP response into + # $HOME/.astropy/cache/astroquery/Vizier, ~2 MB per query. The + # workflow already caches the RESULT as a FITS catalogue on scratch + # and skips the query when it hits, so the pickle is pure duplicate — + # and at campaign scale (~25k exposures) it is ~50 GB against a + # 50 GB home quota. Home is for source and config, not for a second + # copy of the survey. result = v.query_region( - coord, radius=radius_arcmin * u.arcmin, catalog=cat_id + coord, radius=radius_arcmin * u.arcmin, catalog=cat_id, + cache=False, ) if len(result) > 0: print( diff --git a/tests/unit/test_merge_sep_cats_paths.py b/tests/unit/test_merge_sep_cats_paths.py new file mode 100644 index 000000000..754bc76dd --- /dev/null +++ b/tests/unit/test_merge_sep_cats_paths.py @@ -0,0 +1,40 @@ +"""Chunk-path derivation in the merge_sep_cats module.""" + +import pytest + +from shapepipe.modules.merge_sep_cats_package.merge_sep_cats import chunk_path + + +REL = "./output/run_sp_tile_ngmix_Ng1u/ngmix_runner/output/ngmix-210-282.fits" +ABS = ( + "/scratch/run/tiles/21/210.282/output/run_sp_tile_ngmix_Ng1u" + "/ngmix_runner/output/ngmix-210-282.fits" +) + + +def test_relative_path_unchanged_behaviour(): + """The bash pipeline's relative INPUT_DIR still resolves as before.""" + assert chunk_path(REL, 3) == ( + "./output/run_sp_tile_ngmix_Ng3u/ngmix_runner/output/ngmix-210-282.fits" + ) + + +def test_sharded_absolute_path(): + """Digits in the sharded parent dirs and in the file number are untouched.""" + assert chunk_path(ABS, 2) == ( + "/scratch/run/tiles/21/210.282/output/run_sp_tile_ngmix_Ng2u" + "/ngmix_runner/output/ngmix-210-282.fits" + ) + + +def test_double_digit_chunk(): + assert "Ng12u" in chunk_path(ABS, 12) + + +def test_chunk_one_is_identity(): + assert chunk_path(ABS, 1) == ABS + + +def test_no_run_directory_raises(): + with pytest.raises(ValueError, match="run_"): + chunk_path("/scratch/run/tiles/21/210.282/ngmix-210-282.fits", 2) diff --git a/workflow/README.md b/workflow/README.md new file mode 100644 index 000000000..52e8d58c3 --- /dev/null +++ b/workflow/README.md @@ -0,0 +1,231 @@ +# ShapePipe Snakemake orchestration + +Snakemake workflow that orchestrates real-data ShapePipe runs. It replaces the +`curl_canfar_local.sh → run_job_sp_canfar_v2.0.bash → job_sp_canfar_v2.0.bash` +bash layers and the per-site sbatch reimplementations. **Module code is +untouched**: rules call `shapepipe_run -c ` on the existing config +chains. Design and rationale: +[CosmoStat/shapepipe#848](https://github.com/CosmoStat/shapepipe/issues/848) +(the living PRD). + +Use `workflow/bin/sp` for everything. Bare `snakemake all` outside `sp` is +unsupported: `sp` sets the state directory, the SLURM profile, and +`SP_PHASE`, which the Snakefile needs to build the tile/exposure index at +parse time. Running snakemake directly skips all of that. + +## Quick start (nibi) + +```bash +# One-time: a snakemake env on a SHARED filesystem (/project — NOT /tmp, which +# is node-local; the SLURM executor re-invokes this python inside every job). +uv venv /project/def-mjhudson/cdaley/snakemake-env --python 3.12 +source /project/def-mjhudson/cdaley/snakemake-env/bin/activate +uv pip install 'snakemake>=9,<10' 'snakemake-executor-plugin-slurm>=2.7,<3' + +# Edit workflow/config.yaml: tile_list, run_dir, container, star_cats (the +# star-catalogue cache root). + +# The committed launcher loads apptainer/1.4.5 + the /project venv, so a +# fresh shell always has the right state. +workflow/bin/sp run # bring products on disk up to date with the tile list +workflow/bin/sp report # emit run_report.json now (mid-run is fine) +workflow/bin/sp cancel # scancel this workflow's jobs +workflow/bin/sp container status # which image the jobs will run +``` + +Installed and pinned versions on nibi (`/project/def-mjhudson/cdaley/snakemake-env`, +queried 2026-07-30): `snakemake==9.23.1`, `snakemake-executor-plugin-slurm==2.7.1`. +Pin range: `snakemake>=9,<10`, `snakemake-executor-plugin-slurm>=2.7,<3`. The +v8→v9 breaks matter here: `--use-singularity` became `--sdm`, executors became +plugins, and full `rerun-triggers` became the default. + +Anything other than `run`, `report`, `container`, `cancel` passes straight through to +snakemake with the workflow's profile and state dir — the escape hatch for +`sp --unlock`, `sp --dag`, `sp exp_psf ...`. + +## The container image + +`sp container` owns which image the jobs run inside. Two layers, and the second +only exists if you ask for one: + +* your **cached SIF** (`~/.cache/shapepipe/shapepipe.sif`, `SP_CACHE_DIR` or + `SP_CONTAINER` to move it) — a pristine pull of the published image, private + to you, so nobody else's refresh moves the ground under your running jobs; +* an optional **sandbox** (`~/.cache/shapepipe/sandbox/`, `SP_SANDBOX`) — the + same image unpacked writable, so a `pip install` into it sticks. The escape + hatch for work needing a package the image does not carry yet. + +The Snakefile's `container:` is that resolution, in one order shared by the CLI +and the workflow: **sandbox → cached SIF → the `container:` path in +`config.yaml`**. With an empty cache — the normal case — that lands on the +shared `/project` `.sif` the workflow has always used, so this changes nothing +until you opt in. + +```bash +sp container status # layers present, active one, revision vs HEAD +sp container pull # ghcr.io/cosmostat/shapepipe:develop-runtime +sp container pull --tag docker://... # some other image +sp container sandbox # unpack the SIF writable (opt-in) +sp container exec --writable pip install +sp container exec python -c 'import shapepipe' +sp container resolve # just the path the workflow will run +``` + +`status` reads the image's OCI labels and places its +`org.opencontainers.image.revision` against this checkout's HEAD: +in-sync / behind / ahead / diverged, or unknown when the image carries no label +or the commit was never fetched here. + +**`pull` needs the network.** Compute nodes on Alliance clusters generally have +none, so run it on a login node or inside an `salloc` allocation — never from a +batch job. `pull` and `sandbox` both stage to a sibling path and swap it in, so +an in-flight job never sees a half-written image and a failed rebuild leaves the +one you had intact. + +## Execution: two static invocations + +The exposure job set is data-derived from the tiles' `find_exposures` output, +so it cannot live in the same static DAG that produces it. `sp run` is +therefore two snakemake invocations over one Snakefile: + +1. **PREPARE** — `snakemake prepare_all_tiles`: per-tile static DAG + (`Git_vos → Uz → Fe`), `keep-going` so tile failures are independent. A + nonzero exit here is not fatal to the run — tiles that lost their + exposure list are dropped at the compute parse — but it is a warning: + `SP_MISSING_THRESHOLD` (default 0.0) is the real gate. +2. **COMPUTE** — `snakemake all`: this invocation's *parse* builds the + tile↔exposure index (`build_index.py`, imported at parse time, not a DAG + node) and runs the full tile/exposure compute chain. The index + accumulates across invocations, so appending tiles later changes which + jobs exist without invalidating completed work. + +`sp run` chains both so the UX is one command; both exit codes are checked +and the run fails if either phase failed. + +## Execution mode: one SLURM job per rule + +The profile (`profiles/nibi/config.yaml`) sets `executor: slurm`. Every rule +instance becomes its own SLURM job carrying that rule's own attempt-scaled +resources (`cpus_per_task = threads`, `mem_mb`, `runtime`) and runs inside +`apptainer exec` via the profile's software-deployment method — the workflow +never calls apptainer directly. Snakemake feeds the queue as jobs finish, so +the full campaign's job count never needs to queue at once, and multi-node +scaling is inherent. `jobs:` in the profile caps concurrent submissions at a +fraction of the cluster's per-user submit limit (queried, not invented — see +the comment in the profile file for the query and date). + +The profile intentionally sets no `set-resources` / `set-threads` overrides: +those replace a rule's own values wholesale, which would kill the +attempt-scaled `mem_mb = lambda wc, attempt: ...` OOM retries and the tuned +ngmix thread count. Rules own their own resources; the profile only supplies +defaults for rules that state nothing. + +`group:` labels that fuse short rules (uncompress, merges) into their chunky +neighbours (queue-latency amortization, per the PRD) are **not yet wired**: +they require labels in `workflow/rules/*.smk`, out of scope for this +profile-only pass. + +## Layout + +``` +workflow/ + Snakefile parse-time index load; global container:; onsuccess/onerror report hooks + config.yaml the run: tile list, paths, container, chunk count + bin/sp committed launcher (module load + /project venv + run/report/container/cancel) + rules/ + prepare.smk tile get_images/uncompress/find_exposures + exposure.smk per-exposure: get_images, star_cat, split, mask, psf (no temp()) + tile.smk per-tile: exp forest, merge_headers, mask, detect, vignets, ngmix, merge, make_cat + scripts/ + sp_rule.py the thin per-unit wrapper (isolation furniture, config copy, log-sync, count floor) + 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-floor 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) + 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 +``` + +## How it works + +- **The atom is one rule == one `shapepipe_run` on one unit.** Its single + declared output is that unit's manifest + (`/manifests/.json`), not its product files — a missing + CCD is often legitimate, and at DR6 scale per-CCD declaration means + millions of paths. +- **Manifests are the DAG's currency, and they are success-only.** + `completeness.py check` writes its full verdict — per-runner counts against + floors, scraped failure reasons, the `shapepipe_run` exit status when nonzero + — to the rule's `log:` (`/logs/.json`) on *every* run, and + additionally to the declared `.json` manifest only when that verdict is + a success. So `.json` on disk means "this stage succeeded", and a + resume after an unclean death cannot schedule downstream work on top of a + failure. Nothing unlinks anything: Snakemake deletes a failed job's declared + output natively and never touches its log, which is why the profile runs + *without* `keep-incomplete`. `sp report` reads both dirs — the manifest for + success, the log for failure — and a unit with neither ran nothing. +- **Completeness is a count floor, not a taxonomy.** After a run, + `sp_rule.py` counts products per mandatory runner against + `completeness.py`'s floor and exits nonzero below it. Per-CCD attrition + between floor and `expect` is tolerated. No 3-class taxonomy, no + error-signature whitelist. `--keep-going` isolates a failure to its own + DAG cone. +- **Stores are sharded.** Every tile/exposure runs its own `shapepipe_run` + in `tiles/<2-char prefix>//` or `exp///`. Configs are + committed under `workflow/config/cfis/` and version with the rules that set + the env vars they interpolate — there is no `config_src` knob, and no per-unit + config symlink; `$SP_CONFIG` points straight at the committed directory. +- **Mask star catalogues are built in the DAG.** `exp_star_cat` runs one Vizier + cone query per exposure into the run-independent cache at `star_cats:`, then + fans it out into a real per-unit `star_cat_exp/` directory of 40 per-CCD + symlinks, which `exp_mask` consumes. The directory must be per-unit and real: + the file handler intersects the image numbers it finds across a config's + `INPUT_DIR`s, so a symlink to the whole cache contributes every other + exposure's numbers and the intersection comes out empty. It is a `localrule`, + so the queries run serially in the head process — CDS is never hammered, and + the scheduler never sees a six-second job. The cache makes reruns and later + campaigns free. +- **The index is parse-time data, never a rule input.** Appending tiles + changes which jobs exist without invalidating completed work. +- **Exposure products are not `temp()`.** Exposures overlap tiles, so + `temp()` would cascade destructive reruns when a tile is appended later. + Reclamation is the in-DAG `clean_exposure` rule instead: one job per + exposure, taking every consuming tile's `tile_vignets` manifest as input + (the campaign-wide consumer set comes from the accumulating index), which + deletes the store *and* the exposure's `manifests/` and `logs/`, and leaves a + `cleaned.json` tombstone. Deleting the manifests is what makes a late append + correct: the + appended tile finds an unbuilt chain and regenerates it. The `clean:` flag + in `config.yaml` gates it; flipping it on later reclaims retroactively, + since the missing tombstones schedule exactly the outstanding clean jobs. + The tombstone is written *before* anything is deleted, so a crash can cost + disk but never the record. +- **A finished tile declares no reclaimed exposures.** Deleting an exposure's + manifests would otherwise rerun every other tile that reads it, and those + reruns spread across the exposure-overlap component. So a tile whose + `final_cat` exists drops the exposure manifests that are gone from its input + list, and holds the rest through `ancient()`. This is why the profile runs + with `rerun-triggers: [mtime, params, code, software-env]`: the `input` + 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. +- **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 + campaign. List it under `clean_ignore_tiles:` in `config.yaml` and it leaves + the consumer sets. Retrying an ignored tile later is legal and expensive: + its exposure chains are gone and rebuild from scratch. +- **A reclaimed exposure reports as `cleaned`.** `run_report.py` reads the + absorbed manifests out of `cleaned.json`, so a reclaimed exposure keeps its + per-runner counts and blocks no tile. The logs go with the manifests — a log + claiming `complete` for a store that is gone would contradict the unbuilt + chain the DAG must now see, and its content duplicates the manifest anyway. + The `exp_psf` benchmark tsv lives beside both dirs, not inside either, so + reclamation does not eat the memory-sizing data. +- **Failure is a report, not a gate.** `run_report.py` disk-scans the trees + against the count table and enumerates shortfalls (whole-unit absence vs + per-CCD attrition). It runs standalone — a DAG report node would itself be + poisoned by the failures it must enumerate — and fires automatically from + the COMPUTE invocation's `onsuccess`/`onerror` hooks, or on demand via + `sp report`. diff --git a/workflow/Snakefile b/workflow/Snakefile new file mode 100644 index 000000000..0f2541334 --- /dev/null +++ b/workflow/Snakefile @@ -0,0 +1,499 @@ +"""ShapePipe real-data orchestration — Snakemake workflow. + +Design / rationale: CosmoStat/shapepipe#848 (the living PRD), D1-D5. + +`sp run` is TWO snakemake invocations over this one Snakefile: + + SP_PHASE=prepare snakemake prepare_all_tiles # Git -> Uz -> Fe, per tile + SP_PHASE=compute snakemake all # everything else + +They are two because the exposure job set is *data-derived*: it comes from the +tiles' find_exposures output, and a Snakemake DAG is fixed at parse time. The +join between them is the tile<->exposure index, built HERE at parse time of +invocation 2 (build_index.build(), imported — there is no `sp index` verb) and +loaded into plain dicts. The index is never a rule input, so appending tiles and +rebuilding it changes which jobs exist without invalidating completed work; it +ACCUMULATES across invocations, so a later clean_exposure (S5) sees every +consuming tile of the whole campaign, not just this tile list. + +The atom (D2): one rule == one `shapepipe_run` on one unit; its single declared +output is its MANIFEST (`/manifests/.json`), written by +`completeness.py check`. Product files are not declared — a missing CCD is often +legitimate, and at DR6 scale per-CCD declaration means millions of paths. + +Failure evidence rides on the rule's `log:` (`/logs/.json`, +`unit_log` below): the check writes its full verdict there on every run, and the +manifest only when that verdict is a success. Snakemake deletes a failed job's +declared output and never touches its log, so the two directives say exactly the +two things this workflow needs — "this stage succeeded" and "here is what +happened". +""" + +import hashlib +import json +import os +import sqlite3 +import sys +from pathlib import Path + +from snakemake.exceptions import WorkflowError + +# Resolved relative to THIS file, not the working directory: snakemake runs with +# --directory on /scratch (bin/sp) so .snakemake/ state never lands on /project +# (group quota is a hard 27/27 TiB — a metadata write mid-run died on it live). +configfile: str(Path(workflow.snakefile).parent / "config.yaml") + +# Every job's shell runs inside this container (apptainer software-deployment in +# the profile); the user never types apptainer. WHICH image is the one resolution +# order `sp container` exposes — this user's writable sandbox if they built one, +# else their cached SIF if they pulled one, else `container:` above. An empty +# cache therefore lands on exactly the shared /project image the workflow has +# always run, and a package installed into a sandbox reaches the jobs too. +sys.path.insert(0, str(Path(workflow.snakefile).parent / "scripts")) +import container as _container # noqa: E402 + +# Loud at PARSE time, both ways: a broken SP_CONTAINER override, or nothing to +# resolve at all (empty cache AND no `container:` key). The empty string the +# resolver returns for "none" is a valid-looking container directive that passes +# dry-run and only fails once jobs reach a node. +try: + _image, _kind = _container.resolve_image() +except _container.ContainerError as _exc: + raise WorkflowError(str(_exc)) +if _kind == "none": + raise WorkflowError( + f"No container image resolved: no sandbox, no cached SIF, and no " + f"'{_container.CONFIG_KEY}:' key in {_container.CONFIG_FILE}. " + f"Run `sp container pull`, or restore the key.") + +container: _image + +# --- paths ----------------------------------------------------------------- +# TWO ROOTS (D5). RUN_DIR is the scratch root: bulk intermediates, sized so a +# batch finishes inside the purge window. PRODUCTS_DIR is the persistent root +# for the durable, low-volume products — the final catalogues, the index, the +# report. Snakemake's own state is the one durable-looking thing that stays on +# scratch (bin/sp explains why: tens of thousands of small, hot metadata writes +# against a backed-up, file-count-limited filesystem). +RUN_DIR = Path(config["run_dir"]) +# Defaults to RUN_DIR so a scratch-only run (a fixture, a smoke test) needs no +# second path: one root, exactly the pre-D5 layout. +PRODUCTS_DIR = Path(config.get("products_dir") or RUN_DIR) +# Run-independent root for the mask star catalogues: the HEALPix chunk store the +# star_catalogue rule fills, and the per-exposure cuts exp_star_cat makes from it +# (config.yaml explains the placement). +STAR_CATS = Path(config["star_cats"]) +# This checkout's own tree: the star-cat rules run in the container against its +# src/ (shapepipe.utilities.vizier and .cfis), for the same reason CONFIG_DIR is +# the committed config dir — script, library and rule are one artefact and +# version together. +REPO_DIR = Path(workflow.basedir).parent +INDEX_DB = Path(config["index_db"]) +SCRIPTS = Path(workflow.basedir) / "scripts" +# The config chain is the repo's own committed dir BY CONSTRUCTION (D2): the +# configs and the rules that set the env vars they interpolate are one artefact +# and must version together. Hence no `config_src` knob. +CONFIG_DIR = Path(workflow.basedir) / "config" / "cfis" + +sys.path.insert(0, str(SCRIPTS)) +import build_index # noqa: E402 +from completeness import STAGE_DIR # noqa: E402 + +# SP_PHASE is set by bin/sp and NOWHERE else: `prepare`/`compute` on the two +# invocations of `sp run`, `passthrough` on the escape hatch (`sp --unlock`, `sp +# --dag`, `sp exp_psf ...`). It gates the two parse-time side effects — the index +# build and the report hooks — so that a passthrough parse never mutates durable +# state or dies on a threshold it was not asked about. +# +# Requiring it to be SET is the point of the check below: a bare `snakemake` +# would otherwise parse as a passthrough, find `rule all` gated on an index it +# never built, and exit 0 on an empty DAG. Jobs re-parse this file under the +# slurm executor and inherit the head process's environment, so a submitted job +# always carries the launching phase. +PHASE = os.environ.get("SP_PHASE", "") +if not PHASE: + raise WorkflowError( + "SP_PHASE is not set — run the workflow through workflow/bin/sp " + "(`sp run`), which exports it. A bare snakemake invocation would " + "silently build an empty DAG.") +if PHASE not in ("prepare", "compute", "passthrough"): + raise WorkflowError( + f"SP_PHASE={PHASE!r} is not one of prepare, compute, passthrough.") + +with open(config["tile_list"]) as f: + TILES = [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 +# ngmix job per tile. We take the count and drive our own integer `chunk` +# wildcard rather than snakemake's `{scatteritem}` ("3-of-8") token, because the +# chunk number is not ours alone: it names the run dir +# (`run_sp_tile_ngmix_Ngu`, from the template's RUN_NAME), and +# merge_sep_cats derives chunks 2..N from chunk 1's path by substituting the +# chunk number in that run-directory name — which only works for bare integers. +scattergather: + ngmix=int(config.get("ngmix_chunks", 8)) + +NGMIX_CHUNKS = workflow._scatter["ngmix"] + +# --- parse-time index build + load (D1) ------------------------------------ +# The COMPUTE invocation's parse IS the index build, and it runs UNCONDITIONALLY +# there — no "some Fe output exists" guard. That guard used to make a +# totally-failed prepare produce an empty index, an empty DAG and a green exit 0; +# with it gone, zero Fe outputs means a missing fraction of 1.0, which trips the +# SP_MISSING_THRESHOLD gate and fails loudly, as it should. +# +# In every other phase (prepare, or unset for a passthrough invocation) the parse +# builds NOTHING and only loads whatever index is already on disk. +# +# is_main_process matters as much as PHASE: the SLURM executor re-invokes +# snakemake inside every job, that re-invocation parses this file again, and it +# inherits SP_PHASE=compute from the submitting environment. Without the guard, +# every one of the run's jobs re-runs the build — hundreds of concurrent sqlite +# writers on Lustre, which is exactly the "database is locked" storm that killed +# the first real run (2026-07-31). Job parses only LOAD the index below. +if PHASE == "compute" and workflow.is_main_process: + build_index.build( + TILES, RUN_DIR, INDEX_DB, + missing_threshold=float(os.environ.get("SP_MISSING_THRESHOLD", "0.0"))) + +# EXP: exposure base-id -> original name (2605805 -> 2605805p; the name goes +# verbatim into the fabricated per-unit exp_numbers list so get_images matches +# .fits.fz in the store). TILE_EXP: tile -> [exp_ids]. +# EXP_TILES is the inverse edge — the CAMPAIGN-WIDE consumer set clean_exposure +# is keyed on (D5). +EXP, TILE_EXP, EXP_TILES = {}, {}, {} +if INDEX_DB.exists(): + # timeout=60: Lustre lock handoffs are slow; the default 5 s trips on + # nothing more sinister than a reader in another job's parse. + _con = sqlite3.connect(INDEX_DB, timeout=60) + EXP = dict(_con.execute("SELECT exp_id, name FROM exposures")) + for _tile, _exp in _con.execute("SELECT tile_id, exp_id FROM tile_exposures"): + TILE_EXP.setdefault(_tile, []).append(_exp) + EXP_TILES.setdefault(_exp, []).append(_tile) + _con.close() + +# Tiles this run can actually compute: declared AND indexed. The index spans the +# campaign, so it is intersected with the declared list, not used as it. +TILES_READY = [t for t in TILES if t in TILE_EXP] +READY_SET = set(TILES_READY) + +# A compute invocation with nothing to compute is never a success. Without this, +# an empty intersection yields `rule all` with no inputs, an empty DAG and exit +# 0 — the silent green run the threshold gate exists to prevent. +if PHASE == "compute" and not TILES_READY: + raise WorkflowError( + f"No declared tile has an indexed exposure list: 0 of {len(TILES)} tiles " + f"are ready to compute (index {INDEX_DB}). Run the prepare phase first " + f"(`sp run`), or check {INDEX_DB.parent / 'missing.json'}.") + +wildcard_constraints: + tile = r"\d{3}\.\d{3}", + exp = r"\d{6,7}", + shard = r"\d{2}", + chunk = r"\d+", + +# --- the sharded stores (D2) ---------------------------------------------- +# tiles/<2-char prefix>// and exp/<2-char prefix>// — no directory +# exceeds ~1k entries at full-UNIONS scale. Rules carry the shard as its own +# wildcard because an output pattern cannot compute it; every path the DAG uses +# is built by these helpers, so a mismatched (shard, id) pair is never requested. +TILE_DIR = str(RUN_DIR / "tiles" / "{shard}" / "{tile}") +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}") + +def tile_dir(tile): + return f"{RUN_DIR}/tiles/{tile[:2]}/{tile}" + +def exp_dir(exp): + return f"{RUN_DIR}/exp/{exp[:2]}/{exp}" + +def tile_manifest(tile, stage): + return f"{tile_dir(tile)}/manifests/{stage}.json" + +def exp_manifest(exp, stage): + return f"{exp_dir(exp)}/manifests/{stage}.json" + +def forest_dir(tile): + return f"{tile_dir(tile)}/exp_forest" + +def final_cat(tile): + """The campaign's science product, and its tile-finished marker. + + It lives on the PERSISTENT root, not the scratch one, and both halves of + that sentence are load-bearing. As a product: it is what the campaign is + for, and a 60-day purge must not eat it. As a marker: tile_finished() keys + on this path to cut a finished tile's exposure edges (D5), so if a purge + could remove it, finished tiles would re-declare inputs against exposure + stores that reclamation deleted — the rerun avalanche the cut exists to + prevent, arriving by way of the purge instead. + """ + return f"{PRODUCTS_DIR}/tiles/{tile[:2]}/{tile}/final_cat-{tile}.fits" + +def unit_num(unit): + """$SP_UNIT_NUM: ShapePipe's image-number convention, dot -> dash, leading + dash (tile ``210.282`` -> ``-210-282``; exposure ``2605805`` -> ``-2605805``). + The RULES do this transform; the configs just interpolate $SP_UNIT_NUM into + NUMBER_LIST (the `set_config_number_list` mechanism that replaced the retired + -e/--exclusive flag, #746).""" + return "-" + unit.replace(".", "-") + +# Content hash of completeness.py, computed once at parse time and carried as a +# param on every rule: the default rerun-triggers' `code` trigger hashes only the +# rule's own shell string, NOT external scripts it calls — without this, a fix to +# the count table silently leaves stale manifests in place (bitten live). Scoped +# to completeness.py alone, the one script every shell line runs; build_forest.py +# gets its own hash, on the forest rule only. +SCRIPT_HASH = hashlib.md5((SCRIPTS / "completeness.py").read_bytes()).hexdigest()[:12] +FOREST_HASH = hashlib.md5((SCRIPTS / "build_forest.py").read_bytes()).hexdigest()[:12] +CLEAN_HASH = hashlib.md5((SCRIPTS / "clean_exposure.py").read_bytes()).hexdigest()[:12] +# Same argument for star_cats.py, which both star-cat rules call: their params +# otherwise fingerprint nothing but paths, so an edit to the chunking or the cut +# would never rerun them. ONE hash for both rules because it is one script — and +# that is also why fetch and cut live in one module (they must agree on which +# pixel holds which star). The hash does NOT key the store path — see +# config.yaml's star_cats block on clearing the store after a semantic change. +STAR_CAT_HASH = hashlib.md5( + (SCRIPTS / "star_cats.py").read_bytes()).hexdigest()[:12] + +# --- exposure reclamation (D5, S5) ----------------------------------------- + + +def flag(value, default=False): + """Truthiness for a config value that may arrive as a STRING. + + `--config clean=false` delivers the string "false", and every non-empty + string is truthy in Python — a plain bool() read that as ON and scheduled + the deletions the user had just switched off. YAML booleans pass through + unchanged; only strings are parsed, and an unparseable one is an error, not + a guess. + """ + if value is None: + return default + if isinstance(value, str): + v = value.strip().lower() + if v in ("1", "true", "yes", "on"): + return True + if v in ("0", "false", "no", "off", ""): + return False + raise WorkflowError(f"Cannot read {value!r} as a boolean (use true/false).") + return bool(value) + + +# `clean:` in config.yaml gates the whole mechanism. Off => the rule generates no +# jobs at all (nothing requests a tombstone); flipping it on later reclaims +# RETROACTIVELY, because the exposures already cleaned are exactly the ones with +# a tombstone, so the missing tombstones schedule exactly the clean jobs. +# Only ever active under SP_PHASE=compute: the prepare and passthrough parses +# read no index of their own and must schedule no deletions. +CLEAN = flag(config.get("clean", False)) and PHASE == "compute" + +# Reclamation REQUIRES the `input` rerun-trigger to be off (profiles/nibi sets +# the list; tile.smk explains why). The two are already tied together: CLEAN is +# gated on SP_PHASE=compute, which only workflow/bin/sp sets, and bin/sp always +# launches with that profile. A bare snakemake invocation without SP_PHASE +# schedules no clean job at all, so it cannot meet the incompatible combination. +# There is no runtime assertion because the trigger set is not readable at parse +# time (workflow.dag_settings is None until the DAG is built). + +# Tiles that must not pin an exposure's store. See config.yaml: a permanently +# failed tile otherwise holds every exposure it touches (~80) forever, because +# its vignets manifest will never exist and its exposures are therefore never +# eligible. Listing it here drops it from the consumer sets. +CLEAN_IGNORE_TILES = set(config.get("clean_ignore_tiles") or []) + + +def tombstone(exp): + """The clean_exposure output. Lives BESIDE manifests/ and logs/, not inside + either: the clean job deletes both wholesale, and `sp report` scans it.""" + return f"{exp_dir(exp)}/cleaned.json" + + +def clean_consumers(exp): + """Every tile in the campaign that reads this exposure — the set whose + vignets must all exist before the store may go — minus the ignored tiles.""" + return sorted(t for t in EXP_TILES.get(exp, []) + if t not in CLEAN_IGNORE_TILES) + + +def clean_targets(): + """Which exposures this invocation may clean. + + An exposure is eligible only when every consuming tile is either in this + run's scope or has already produced its vignets on disk. Without that test, + requesting a tombstone for an exposure shared with a LATER batch would drag + that batch's whole tile chain into this DAG through the clean rule's input — + scope expansion by cleanup, which is not a trade anyone asked for. Ineligible + exposures are simply skipped; the invocation that finishes their last + consumer picks them up. Deferral, never loss. + + Consumer sets are the IGNORE-FILTERED ones (clean_consumers), so a tile in + `clean_ignore_tiles` neither gates eligibility nor appears in the job's + input — which is the whole point of that list. + """ + if not CLEAN: + return [] + out = [] + for exp, raw in EXP_TILES.items(): + if not raw: + continue + # An empty set AFTER filtering means every consumer is ignored: nothing + # is left that could ever read this exposure, so it is eligible now. + # That is the whole point of clean_ignore_tiles — all() of an empty set + # is True, and it is true here in the intended sense. + tiles = clean_consumers(exp) + if all(t in READY_SET or Path(tile_manifest(t, "tile_vignets")).exists() + for t in tiles): + out.append(tombstone(exp)) + return sorted(out) + +# --- the shell every rule runs (D2) ---------------------------------------- + +_THREAD_CAPS = " ".join( + f"{k}={v}" for k, v in ( + ("OMP_NUM_THREADS", 1), ("OPENBLAS_NUM_THREADS", 1), + ("MKL_NUM_THREADS", 1), ("NUMEXPR_NUM_THREADS", 1), + ("MALLOC_ARENA_MAX", 2), ("MALLOC_TRIM_THRESHOLD_", 0))) + + +def unit_pre(stage, level, unit, *, exp_name=None, forest=None, env=None, + pre_run=()): + """The unit-furniture + environment prologue, as bash. + + Returned as a rule ``params`` value, NEVER inlined into the ``shell:`` + string: snakemake formats a shell string ONCE, so a ``{output}``/``{threads}`` + placeholder inside a params value would survive literally — and, conversely, + the literal ``${SP_NGMIX_CHUNK}`` braces this prologue needs would blow up + that formatting if they lived in the shell string. Params values are + substituted after formatting, so both hazards go away together. + + What it materialises (the proven v2.0 isolation-by-work-dir-content, NOT + -e/--exclusive): + * ``output/``, ``manifests/`` and ``logs/`` — the last two are a SIBLING + pair, not one dir with two naming conventions: clean_exposure deletes + BOTH wholesale, because a log attesting "complete" for a store that has + been reclaimed contradicts the unbuilt chain the DAG must now see. (The + benchmark tsv is deliberately outside both: it measures the job, it does + not claim anything about the store, so reclamation must not eat it.) + * tile: ``tile_numbers.txt`` (dot format — what get_images reads); + * exposure: a fabricated pseudo-Fe ``exp_numbers-000-000.txt`` holding the + ORIGINAL exposure name from the index (``2605805p``), so get_images + matches ``.fits.fz`` in the store — the bare base id matches + nothing. Written UNCONDITIONALLY: an exists-guard once pinned a stale + pre-fix file with the bare id. + There is no per-unit ``cfis`` symlink any more: $SP_CONFIG points straight at + the committed config dir, and ``star_cat_exp`` is a real per-unit directory + built by the ``exp_star_cat`` rule, not a symlink into a shared pool. + + Finally it ``rm -rf``s this stage's own fixed run dir — ShapePipe's + FileHandler raises on an existing run dir, and it is how a rerun never sees + stale products (D2: the job clears its run dir at start). + """ + work = tile_dir(unit) if level == "tile" else exp_dir(unit) + _, subdir = STAGE_DIR[stage] + lines = [ + "set -euo pipefail", + f"export SP_RUN='{work}'", + f"export SP_UNIT_NUM='{unit_num(unit)}'", + f"export SP_CONFIG='{CONFIG_DIR}'", + # Also set via apptainer-args in the profile; kept here so a hand-run of + # this same line outside snakemake behaves identically. + f"export {_THREAD_CAPS}", + 'mkdir -p "$SP_RUN/output" "$SP_RUN/manifests" "$SP_RUN/logs"', + ] + if forest: + lines.append(f"export SP_EXP='{forest}'") + for k, v in (env or {}).items(): + lines.append(f"export {k}='{v}'") + + if level == "tile": + lines.append(f"printf '%s\\n' '{unit}' > \"$SP_RUN/tile_numbers.txt\"") + else: + fe = "$SP_RUN/output/run_sp_tile_Fe/find_exposures_runner/output" + lines += [f'mkdir -p "{fe}"', + f"printf '%s\\n' '{exp_name or unit}' > \"{fe}/exp_numbers-000-000.txt\""] + + lines += list(pre_run) + lines += [f'rm -rf "$SP_RUN/output/{subdir}"', 'cd "$SP_RUN"'] + return "\n".join(lines) + + +def sp_shell(stage, config_name): + """The rule's shell string: prologue, one shapepipe_run, one completeness check. + + ``{threads}``, ``{output}`` and ``{log}`` are placeholders HERE and nowhere + else (see unit_pre). ``-b {threads}`` makes SMP fork width and cpus_per_task + one number by construction (D4). + + The check runs even when shapepipe_run failed, so a failed unit still leaves + a record of why: the verdict goes to the rule's `log:` on EVERY run, and to + the declared `.json` manifest only when it is a success. Snakemake + deletes a failed job's output and never touches its log, so the manifest + stays success-only currency and the log is the evidence `sp report` reads. + + The verdict is COMPOSED, which is why `--job-rc "$rc"` is passed: the count + floors and shapepipe_run's exit status are two independent statements about + the same job, and the record must reflect both. Without it, a job whose + counts cleared their floors but whose shapepipe_run died (a late runner + raising after the counted ones wrote their files) publishes a SUCCESS + manifest and still exits nonzero — snakemake deletes that manifest as a + failed job's output, and the log is left claiming "complete" for a stage + with nothing to show for it. With the rc composed in, both files agree that + the job failed. + """ + return ( + "{params.pre}\n" + "rc=0\n" + f"shapepipe_run -c \"$SP_CONFIG/{config_name}\" -b {{threads}} || rc=$?\n" + f"python {SCRIPTS}/completeness.py check {stage} {{output.manifest}}" + " --log {log} --job-rc \"$rc\" || rc=1\n" + "exit $rc\n" + ) + + +include: "rules/prepare.smk" +include: "rules/exposure.smk" +include: "rules/tile.smk" + +# --- top-level targets ------------------------------------------------------ +# The aggregation targets, clean_exposure, star_catalogue and exp_star_cat run in +# the head process. clean_exposure is seconds of rmtree and hangs off `all`; +# exp_star_cat is seconds of local FITS work; both would otherwise be ~20k sbatch +# submissions at DR6 scale for work shorter than the scheduling latency. +# star_catalogue is one job either way, and local keeps its CDS concurrency the +# explicit number its thread pool sets (see exposure.smk). +# +# star_catalogue and exp_star_cat are MID-CHAIN localrules, so they must stay out +# of any future `group:` label: a local job cannot be fused into a submitted group. +localrules: all, prepare_all_tiles, clean_exposure, star_catalogue, exp_star_cat + +rule all: + input: + [final_cat(t) for t in TILES_READY], + clean_targets(), + +# Invocation 1 — the static per-tile DAG, known from the tile list alone. +# keep-going makes tile failures independent; the ones that lose their exposure +# list are dropped by the index build at invocation 2's parse. +rule prepare_all_tiles: + input: + [tile_manifest(t, "tile_find_exposures") for t in TILES] + +# --- report hooks ----------------------------------------------------------- +# run_report is NOT a DAG node (a descendant of every job would be poisoned by +# any hard failure — the exact case it exists for). It is a standalone script, +# emitted automatically at the end of the COMPUTE invocation, runnable any time +# via `sp report`. +def _report(status): + shell(f"python {SCRIPTS}/run_report.py --run-dir {RUN_DIR} " + f"--index {INDEX_DB} --status {status} || true") + +if PHASE == "compute": + + onsuccess: + _report("success") + + onerror: + _report("error") diff --git a/workflow/bin/sp b/workflow/bin/sp new file mode 100755 index 000000000..921b95d8d --- /dev/null +++ b/workflow/bin/sp @@ -0,0 +1,124 @@ +#!/usr/bin/env bash +# sp — the committed launcher for the ShapePipe Snakemake workflow (PRD #848 D1). +# +# Three verbs, nothing else: +# +# sp run [ARGS...] bring the products on disk up to date with the tile list. +# Two snakemake invocations over one Snakefile: +# 1. PREPARE snakemake prepare_all_tiles +# 2. COMPUTE snakemake all <- its PARSE builds the index +# ARGS (--jobs, -n, --forcerun, ...) pass through to BOTH. +# sp report [ARGS...] emit run_report.json now (mid-run is fine). +# sp container VERB manage the image every job runs inside: pull, status, +# sandbox, exec, resolve. `sp container --help` documents +# the two layers and the resolution order (sandbox -> your +# cached SIF -> the `container:` path in config.yaml). +# `pull` needs the network: run it on a login node or in an +# salloc allocation, never from a batch job. +# +# Anything else is passed straight through to snakemake with the same profile and +# state dir (the escape hatch: `sp --unlock`, `sp exp_psf ...`, `sp --dag`). +# +# It also loads the apptainer module (snakemake resolves `apptainer` via PATH at +# job runtime) and activates the snakemake venv on the shared /project FS (the +# executor re-invokes it inside jobs, so it cannot live on a node-local path). +# One entry point, so a fresh tmux or a restart after a crash always launches +# with the right state. +set -euo pipefail + +HERE="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" # workflow/ +REPO="$(dirname "$HERE")" +VENV="${SP_SNAKEMAKE_ENV:-/project/def-mjhudson/cdaley/snakemake-env}" +PROFILE="$REPO/profiles/nibi" +SCRIPTS="$HERE/scripts" +CONFIG="$HERE/config.yaml" + +module load apptainer/1.4.5 2>/dev/null || true +# shellcheck disable=SC1091 +source "$VENV/bin/activate" + +# Minimal scalar reader for workflow/config.yaml (key: value, no nesting). +cfg() { sed -n "s/^$1:[[:space:]]*//p" "$CONFIG" | head -1; } +RUN_DIR="$(cfg run_dir)"; INDEX_DB="$(cfg index_db)" + +# Snakemake state (.snakemake: metadata, locks, incomplete markers) lives NEXT TO +# THE RUN on /scratch, never on the persistent root — the one exception to D5's +# "durable, low-volume products go on /project". It is neither: one small, hot +# metadata file per output, rewritten on every job, which at DR6 is >170k files +# of churn against a backed-up filesystem with a 1M-inode group quota. It is also +# reconstructible — losing it costs a re-parse, not a re-run. (The lesson arrived +# the hard way: a hard 27/27 TiB group quota killed a metadata write mid-run, +# live. The quota has since eased; the placement is right on its own merits.) +# --directory only moves state: all data paths are absolute, and the Snakefile +# resolves its own configfile. +STATE_DIR="${SP_STATE_DIR:-${RUN_DIR}-state}"; mkdir -p "$STATE_DIR" + +# SP_MISSING_THRESHOLD gates the compute parse's index build: the fraction of +# declared tiles allowed to be missing their exposure list (default 0.0). +export SP_MISSING_THRESHOLD="${SP_MISSING_THRESHOLD:-0.0}" + +# --snakefile pins the workflow to this checkout: sp must work from any cwd +# (an sbatch head job starts in the submission directory, not the repo). +# +# SP_PHASE is REQUIRED by the Snakefile (a bare `snakemake` would build an empty +# DAG and exit 0). `sp run` sets prepare/compute on its two invocations; every +# other verb and the escape hatch fall through to `passthrough`, which parses the +# index without building it and schedules no side effects. +sm() { + SP_PHASE="${SP_PHASE:-passthrough}" \ + snakemake --snakefile "$HERE/Snakefile" --profile "$PROFILE" \ + --directory "$STATE_DIR" "$@" +} + +cmd="${1:-}" +case "$cmd" in + run) + shift + # PREPARE failing is NOT fatal to the run: keep-going means a failed tile + # poisons only its own cone, and the tiles that lost their exposure list are + # dropped at the compute parse. The real gate is SP_MISSING_THRESHOLD, which + # the compute parse's index build enforces over the WHOLE tile list — so we + # record the failure and go on rather than letting `set -e` abort here. + prep_rc=0 + SP_PHASE=prepare sm prepare_all_tiles "$@" || prep_rc=$? + if [ "$prep_rc" -ne 0 ]; then + echo "" >&2 + echo "############################################################" >&2 + echo "## WARNING: the PREPARE phase exited $prep_rc." >&2 + echo "## Some tiles may be missing their exposure list and will be" >&2 + echo "## dropped from the compute DAG. Continuing to COMPUTE; the" >&2 + echo "## SP_MISSING_THRESHOLD gate (now $SP_MISSING_THRESHOLD) decides" >&2 + echo "## whether that is tolerable." >&2 + echo "############################################################" >&2 + echo "" >&2 + fi + + comp_rc=0 + SP_PHASE=compute sm all "$@" || comp_rc=$? + if [ "$prep_rc" -ne 0 ] || [ "$comp_rc" -ne 0 ]; then + exit "$([ "$comp_rc" -ne 0 ] && echo "$comp_rc" || echo "$prep_rc")" + fi + ;; + report) + shift + python "$SCRIPTS/run_report.py" --run-dir "$RUN_DIR" --index "$INDEX_DB" \ + --status manual "$@" + ;; + container) + # The image layer, outside snakemake entirely: the script is stdlib-only so + # it runs on the bare host, and the apptainer module is already loaded above. + shift + python "$SCRIPTS/container.py" "$@" + ;; + cancel) + # Kept only because it is two lines: scancel this workflow's jobs by name + # before an --unlock. Not part of the design surface. + run="${2:?usage: sp cancel }" + squeue --me --noheader --format='%i %j' \ + | awk -v r="$run" '$2 ~ r {print $1}' | xargs -r scancel + echo "cancelled jobs matching '$run'; safe to --unlock / rerun now" + ;; + *) + sm "$@" + ;; +esac diff --git a/workflow/config.yaml b/workflow/config.yaml new file mode 100644 index 000000000..e67b6867e --- /dev/null +++ b/workflow/config.yaml @@ -0,0 +1,114 @@ +# Run configuration for the ShapePipe Snakemake workflow. +# +# A "run" is declared by a tile list plus the paths below. Everything here is +# read at parse time; none of it is a rule input, so editing it (e.g. appending +# tiles) never invalidates completed work — it only changes which jobs exist. + +# The tile list that scopes this run (one "IDra.IDdec" per line). +# The campaign grows by appending to this file — parse-time config, so +# completed work is never invalidated. Current contents: the 34-tile smk-g4 +# benchmark set, identical to the tile list smk-g3 finished (the 186/187 quad +# plus the first 30 p3-batch1 tiles), so the two campaigns' catalogues are +# comparable object for object. It lives on the persistent root because it +# defines the campaign and must outlive the scratch purge. +# (The 210/211 quad used earlier is unusable: its tile images are symlinks into +# anaennis' moved processed_tiles tree — ~8.5k of the 10.3k staged tiles are +# broken.) +tile_list: /project/def-mjhudson/cdaley/sp-products/smk-g4/tiles34.txt + +# The container every job runs inside (apptainer software-deployment in the profile). +container: /project/def-mjhudson/cdaley/containers/shapepipe-develop-runtime.sif + +# THE TWO ROOTS (D5). +# +# run_dir is the SCRATCH root and the $SP_RUN every config interpolates: bulk +# intermediates, sized so a batch finishes inside the 60-day purge window. The +# sharded per-unit stores live under it: +# /tiles/<2-char prefix>// and /exp/// +run_dir: /scratch/cdaley/shapepipe-output/smk-g4 + +# products_dir is the PERSISTENT root: the durable, low-volume products — the +# final catalogues (/tiles///final_cat-.fits, +# mirroring the scratch tree shard for shard), the index, the report. ~32-46 MB +# per tile, so a full DR6 campaign is a few hundred GB. +# +# The final catalogue is also the tile-finished MARKER that cuts a finished +# tile's exposure edges (D5), which is the second reason it cannot sit on +# scratch: a purge would not merely lose a product, it would make every finished +# tile re-declare inputs against exposure stores reclamation already deleted. +# +# Unset means "one root": products land under run_dir, exactly the pre-D5 +# layout, which is what a fixture or smoke test wants. +# +# Snakemake's own state is the one durable-looking thing that stays on scratch +# (-state; bin/sp explains why). +products_dir: /project/def-mjhudson/cdaley/sp-products/smk-g4 + +# 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 +# are one artefact and must version together, so the dir is fixed by construction. + +# Pre-staged inputs (P3 data already on /project; get_images RETRIEVE=symlink). + +# The mask star-catalogue root — run-independent, shared by every campaign, and +# holding two things: +# /I_305_out/nside32/star_chunk-.fits the SKY store, one GSC 2.3 +# query per HEALPix chunk (~3.4 deg^2, ~25k rows), written by +# `star_catalogue` over the tile list's footprint and never fetched twice; +# /exp/star_cat-.fits the per-exposure cuts +# `exp_star_cat` makes from those chunks, with no network at all. +# Network therefore scales with SKY AREA, not exposure count: exposures overlap +# ~7-10 deep, so a full-UNIONS footprint is ~1.5k queries against ~25k exposures. +# +# On the PERSISTENT root: the sky store is a durable science product bought with +# ~1.5k catalogue-server queries at DR6 scale, and re-buying it after a scratch +# purge is the one cost in this workflow that cannot be paid with local compute. +# (It sat on scratch through smk-g3 only because def-mjhudson /project was then +# hard-full at 27/27 TiB.) +# +# THE STORE IS NOT KEYED BY SCRIPT VERSION. A semantic change to +# workflow/scripts/star_cats.py (padding, catalogue ID, column set) does rerun +# both rules — the script's hash is a param on each — but a chunk already on disk +# is skipped and only re-cut. Clear the store by hand when the change must reach +# the data. Changing NSIDE or the catalogue ID is the exception: those name the +# directory, so a change there fetches into a new one beside the old. +star_cats: /project/def-mjhudson/cdaley/sp-products/star-cat-cache + +# The run index, and — sharing its directory — missing.json and run_report.json. +# On the persistent root with the catalogues (D5): the index is the record of +# which tile reads which exposure, so it is what a post-purge reconstruction +# would otherwise have to rebuild from tile headers. +index_db: /project/def-mjhudson/cdaley/sp-products/smk-g4/index/run_index.sqlite + +# 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, +# and leaves `cleaned.json`, which absorbs the manifests — `sp report` reads them +# back out of the tombstone and reports the exposure as `cleaned`. On for smk-g4: +# reclamation is part of what this campaign measures, and at 34 tiles from +# scratch the un-cleaned exposure high-water is the run's binding disk constraint. +clean: true + +# Tiles that may NOT pin an exposure store (default: empty). +# +# An exposure is eligible for cleaning only once EVERY consuming tile has its +# vignets. One permanently-failed tile therefore holds all its exposures +# (~8 measured on P3: 7.9 exposures/tile) for the life of the campaign. A tile listed here is dropped from the +# consumer sets, and its exposures become eligible. +# +# READ THIS BEFORE ADDING A TILE. Ignoring a tile is a decision to give up its +# exposures' stores. If you later retry that tile, those exposure chains are +# gone and will be REBUILT from scratch — get_images, split, mask, psf, per +# exposure. That is correct, and expensive. Ignore a tile when you have decided +# it is dead, not while you are still debugging it. +clean_ignore_tiles: [] + +# ngmix within-tile chunking: static N chunks (closed ID ranges computed +# per-tile, in-job, from the tile's own sexcat). +ngmix_chunks: 8 + +# The container's installed shapepipe is overridden by the prod worktree via +# --env PYTHONPATH in profiles/nibi (settled call 3); no config knob here. +# build_index.py's missing-tile fraction floor is passed by workflow/bin/sp +# (SP_MISSING_THRESHOLD, default 0.0 = any missing tile is fatal). diff --git a/workflow/config/cfis/config_exp_Gie.ini b/workflow/config/cfis/config_exp_Gie.ini new file mode 100644 index 000000000..7c6297cf0 --- /dev/null +++ b/workflow/config/cfis/config_exp_Gie.ini @@ -0,0 +1,99 @@ +# ShapePipe configuration file for: get images + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = False + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_exp_Gie + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = get_images_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names +INPUT_DIR = $SP_RUN + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +# Get exposures +[GET_IMAGES_RUNNER] + +INPUT_DIR = $SP_RUN/output/run_sp_tile_Fe/find_exposures_runner/output + +FILE_PATTERN = exp_numbers + +FILE_EXT = .txt + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + + +# Paths + +# Output path (optional, default is [FILE]:OUTPUT_DIR +# OUTPUT_PATH = input_images + +# Input path where original images are stored. Can be local path or vos url. +# Single string or list of strings +INPUT_PATH = /project/def-mjhudson/unions-wl/exposures, /project/def-mjhudson/unions-wl/exposures, /project/def-mjhudson/unions-wl/exposures + +# Input file pattern including tile number as dummy template +INPUT_FILE_PATTERN = 000000, 000000.weight, 000000.flag + +# Input file extensions +INPUT_FILE_EXT = .fits.fz, .fits.fz, .fits.fz + +# Input numbering scheme, python regexp +INPUT_NUMBERING = \d{6} + +# Output file pattern without number +OUTPUT_FILE_PATTERN = image-, weight-, flag- + +# Method to retrieve images, one in 'vos', 'symlink' +RETRIEVE = symlink + +# If RETRIEVE=vos, number of attempts to download +# Optional, default=3 +N_TRY = 3 + +# Retrieve command options, optional +RETRIEVE_OPTIONS = --certfile=$HOME/.ssl/cadcproxy.pem + +#CHECK_EXISTING_DIR = $SP_RUN/output/run_sp_Gie_prev diff --git a/workflow/config/cfis/config_exp_Ma.ini b/workflow/config/cfis/config_exp_Ma.ini new file mode 100644 index 000000000..d5b521080 --- /dev/null +++ b/workflow/config/cfis/config_exp_Ma.ini @@ -0,0 +1,86 @@ +# ShapePipe configuration file for masking of exposures + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_exp_Ma + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = mask_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 4 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +### Mask exposures +[MASK_RUNNER] + +# Parent module +INPUT_DIR = $SP_RUN/output/run_sp_exp_Sp/split_exp_runner/output, $SP_RUN/star_cat_exp + +# Update numbering convention, accounting for HDU number of +# single-exposure single-HDU files +NUMBERING_SCHEME = -0000000-0 + +# Input file patterns: image, weight, external flag, external star catalogue +FILE_PATTERN = image, weight, flag, star_cat + +FILE_EXT = .fits, .fits, .fits, .fits + +# Path of mask config file +MASK_CONFIG_PATH = $SP_CONFIG/config_onthefly.mask + +# External mask file flag, use if True, otherwise ignore +USE_EXT_FLAG = True + +# External star catalogue flag, use external cat if True, +# obtain from online catalogue if False +# True: the cat comes from $SP_RUN/star_cat_exp, the per-unit farm the +# exp_star_cat rule builds (40 per-CCD links to this exposure's one cat). +USE_EXT_STAR = True + +# File name suffix for the output flag files (optional) +PREFIX = pipeline + +# Path to check for existing output mask files +CHECK_EXISTING_DIR = $SP_RUN/output/run_sp_exp_Ma/mask_runner/output diff --git a/workflow/config/cfis/config_exp_Sp.ini b/workflow/config/cfis/config_exp_Sp.ini new file mode 100644 index 000000000..dd27d6ccd --- /dev/null +++ b/workflow/config/cfis/config_exp_Sp.ini @@ -0,0 +1,78 @@ +# ShapePipe configuration file for single-exposures, +# split images + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_exp_Sp + +# Add date and time to RUN_NAME, optional, default: True +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = split_exp_runner + +# Run mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed exposure ID (exp split is a "tile-scheme" stage per sp_rule.py). +NUMBER_LIST = $SP_UNIT_NUM + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 8 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[SPLIT_EXP_RUNNER] + +INPUT_DIR = $SP_RUN/output/run_sp_exp_Gie/get_images_runner/output + +FILE_PATTERN = image, weight, flag + +# Matches compressed single-exposure files +FILE_EXT = .fitsfz, .fitsfz, .fitsfz + +NUMBERING_SCHEME = -0000000 + +# OUTPUT_SUFFIX, actually file name prefixes. +# Expected keyword "flag" will lead to a behavior where the data are saved as int. +# The code also expects the image data to use the "image" suffix +# (default value in the pipeline). +OUTPUT_SUFFIX = image, weight, flag + +# Number of HDUs/CCDs of mosaic +N_HDU = 40 diff --git a/workflow/config/cfis/config_exp_psfex.ini b/workflow/config/cfis/config_exp_psfex.ini new file mode 100644 index 000000000..0af871d8e --- /dev/null +++ b/workflow/config/cfis/config_exp_psfex.ini @@ -0,0 +1,181 @@ +# ShapePipe configuration file for single-exposures. PSFex PSF model. +# Process exposures after masking, from star detection to PSF model. + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_exp_SxSePsfPi +#RUN_NAME = run_sp_exp_SxSePsf + +# Add date and time to RUN_NAME, optional, default: True +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = sextractor_runner, setools_runner, psfex_runner, psfex_interp_runner + + +# Run mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = $SP_RUN/output + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 8 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[SEXTRACTOR_RUNNER] + +# Input from two modules +INPUT_DIR = $SP_RUN/output/run_sp_exp_Sp/split_exp_runner/output, $SP_RUN/output/run_sp_exp_Ma/mask_runner/output + +# Read pipeline flag files created by mask module +FILE_PATTERN = image, weight, pipeline_flag + +# Explicit extensions: a 3-entry FILE_PATTERN override must not fall back on +# the decorator's 4-entry FILE_EXT default (length check fails at startup) +FILE_EXT = .fits, .fits, .fits + +NUMBERING_SCHEME = -0000000-0 + +# SExtractor executable path +EXEC_PATH = source-extractor + +# SExtractor configuration files +DOT_SEX_FILE = $SP_CONFIG/default_exp.sex +DOT_PARAM_FILE = $SP_CONFIG//default.param +DOT_CONV_FILE = $SP_CONFIG/default.conv + +# Use input weight image if True +WEIGHT_IMAGE = True + +# Use input flag image if True +FLAG_IMAGE = True + +# Use input PSF file if True +PSF_FILE = False + +# Use distinct image for detection (SExtractor in +# dual-image mode) if True. +DETECTION_IMAGE = False + +# Distinct weight image for detection (SExtractor +# in dual-image mode) +DETECTION_WEIGHT = False + +# True if photometry zero-point is to be read from exposure image header +ZP_FROM_HEADER = True + +# If ZP_FROM_HEADER is True, zero-point key name +ZP_KEY = PHOTZP + +# Background information from image header. +# If BKG_FROM_HEADER is True, background value will be read from header. +# In that case, the value of BACK_TYPE will be set atomatically to MANUAL. +# This is used e.g. for the LSB images. +BKG_FROM_HEADER = False +# LSB images: +# BKG_FROM_HEADER = True + +# If BKG_FROM_HEADER is True, background value key name +# LSB images: +#BKG_KEY = IMMODE + +# Type of image check (optional), default not used, can be a list of +# BACKGROUND, BACKGROUND_RMS, INIBACKGROUND, MINIBACK_RMS, -BACKGROUND, +# FILTERED, OBJECTS, -OBJECTS, SEGMENTATION, APERTURES +CHECKIMAGE = BACKGROUND, BACKGROUND_RMS + +# File name suffix for the output sextractor files (optional) SUFFIX = tile +SUFFIX = sexcat + +## Post-processing + +# Not required for single exposures +MAKE_POST_PROCESS = FALSE + + +[SETOOLS_RUNNER] + +INPUT_DIR = $SP_RUN/output/run_sp_exp_SxSePsfPi/sextractor_runner/output + +# Note: Make sure this doe not match the SExtractor background images +# (sexcat_background*) +FILE_PATTERN = sexcat + +NUMBERING_SCHEME = -0000000-0 + +# SETools config file +SETOOLS_CONFIG_PATH = $SP_CONFIG/star_selection.setools + + +[PSFEX_RUNNER] + +# Use 80% sample for PSF model +FILE_PATTERN = star_split_ratio_80 + +NUMBERING_SCHEME = -0000000-0 + +# Path to executable for the PSF model (optional) +EXEC_PATH = psfex + +# Default psfex configuration file +DOT_PSFEX_FILE = $SP_CONFIG/default.psfex + +[PSFEX_INTERP_RUNNER] + +# Use 20% sample for PSF validation +FILE_PATTERN = star_split_ratio_80, star_split_ratio_20, psfex_cat + +FILE_EXT = .psf, .fits, .cat + +NUMBERING_SCHEME = -0000000-0 + +# Run mode for psfex interpolation: +# CLASSIC: 'classical' run, interpolate to object positions +# MULTI-EPOCH: interpolate for multi-epoch images +# VALIDATION: validation for single-epoch images +MODE = VALIDATION + +# Column names of position parameters +POSITION_PARAMS = XWIN_IMAGE,YWIN_IMAGE + +# If True, measure and store ellipticity of the PSF (using moments) +GET_SHAPES = True + +# Minimum number of stars per CCD for PSF model to be computed +STAR_THRESH = 22 + +# Maximum chi^2 for PSF model to be computed on CCD +CHI2_THRESH = 2 diff --git a/workflow/config/cfis/config_merge_sep_cats.ini b/workflow/config/cfis/config_merge_sep_cats.ini new file mode 100644 index 000000000..9dda5cf0f --- /dev/null +++ b/workflow/config/cfis/config_merge_sep_cats.ini @@ -0,0 +1,85 @@ +# ShapePipe post-run configuration file: merge separated catalogues + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_Ms + +# Add date and time to RUN_NAME, optional, default: True +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = merge_sep_cats_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed tile ID (e.g. -210-282). +NUMBER_LIST = $SP_UNIT_NUM + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +# NOTE: only chunk 1's ngmix output is listed here; merge_sep_cats_runner +# derives the other chunks' paths itself from N_SPLIT_MAX and this pattern +# (chunk dirs are run_sp_tile_ngmix_Ngu, k=1..N_SPLIT_MAX), all under the +# same fixed output dir. +INPUT_DIR = $SP_RUN/output/run_sp_tile_ngmix_Ng1u/ngmix_runner/output + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 8 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[MERGE_SEP_CATS_RUNNER] + +# Input file pattern(s), list of strings with length matching number of expected input file types +# Cannot contain wild cards +FILE_PATTERN = ngmix + +# FILE_EXT (optional) list of string extensions to identify input files +FILE_EXT = .fits + +# Numbering convention, string that exemplifies a numbering pattern. +NUMBERING_SCHEME = -000-000 + +# WARNING (optional, default is 'error'). Use 'always'/'ignore' to +# display/ignore warnings, and not raise error +WARNING = always + +# Maximum number of separated catalogues per input. +# merge_sep_cats_runner.py reads this with getexpanded (runner .py line ~31), so +# $NGMIX_N_CHUNKS DOES expand here, exactly like ID_OBJ_MIN/MAX in the ngmix +# module. The tile_merge_cats rule exports it from the workflow's scattergather +# chunk count, which is the single source of truth; this value must equal that +# count, so do not replace it with a literal. +N_SPLIT_MAX = $NGMIX_N_CHUNKS diff --git a/workflow/config/cfis/config_onthefly.mask b/workflow/config/cfis/config_onthefly.mask new file mode 100644 index 000000000..7c185c602 --- /dev/null +++ b/workflow/config/cfis/config_onthefly.mask @@ -0,0 +1,86 @@ +# Mask module configuration file for single-exposure images + +## Paths to executables +[PROGRAM_PATH] + +WW_PATH = weightwatcher +WW_CONFIG_FILE = $SP_CONFIG/mask_default/default.ww + +# Indicate cds client executable if no external star catalogue is available +# (e.g. no internet access on run nodes) +CDSCLIENT_PATH = findgsc2.2 + + +## Border mask +[BORDER_PARAMETERS] + +BORDER_MAKE = True + +BORDER_WIDTH = 50 +BORDER_FLAG_VALUE = 4 + + +## Halo mask +[HALO_PARAMETERS] + +HALO_MAKE = True + +HALO_MASKMODEL_PATH = $SP_CONFIG/mask_default/halo_mask.reg +HALO_MAG_LIM = 13. +HALO_SCALE_FACTOR = 0.05 +HALO_MAG_PIVOT = 13.8 +HALO_FLAG_VALUE = 2 +HALO_REG_FILE = halo.reg + + +## Diffraction spike mask +[SPIKE_PARAMETERS] + +SPIKE_MAKE = True + +SPIKE_MASKMODEL_PATH = $SP_CONFIG/mask_default/MEGAPRIME_star_i_13.8.reg +SPIKE_MAG_LIM = 18. +SPIKE_SCALE_FACTOR = 0.3 +SPIKE_MAG_PIVOT = 13.8 +SPIKE_FLAG_VALUE = 128 +SPIKE_REG_FILE = spike.reg + + +## Messier mask +[MESSIER_PARAMETERS] + +MESSIER_MAKE = True + +MESSIER_CAT_PATH = $SP_CONFIG/mask_default/Messier_catalog_updated.fits +MESSIER_SIZE_PLUS = 0. +MESSIER_FLAG_VALUE = 16 + + +## NGC mask +[NGC_PARAMETERS] + +NGC_MAKE = True + +NGC_CAT_PATH = $SP_CONFIG/mask_default/ngc_cat.fits +NGC_SIZE_PLUS = 0. +NGC_FLAG_VALUE = 32 + + + +## Missing data parameters +[MD_PARAMETERS] + +MD_MAKE = False + +MD_THRESH_FLAG = 0.3 +MD_THRESH_REMOVE = 0.75 +MD_REMOVE = False + + +## Other parameters +[OTHER] + +TEMP_DIRECTORY = .temp + +KEEP_REG_FILE = False +KEEP_INDIVIDUAL_MASK = False diff --git a/workflow/config/cfis/config_tile_Fe.ini b/workflow/config/cfis/config_tile_Fe.ini new file mode 100644 index 000000000..9546f062e --- /dev/null +++ b/workflow/config/cfis/config_tile_Fe.ini @@ -0,0 +1,76 @@ +# ShapePipe configuration file for: find exposures + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = False + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Fe + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = find_exposures_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed tile ID (e.g. -210-282). +NUMBER_LIST = $SP_UNIT_NUM + +# Input directory, containing input files, single string or list of names +INPUT_DIR = $SP_RUN + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +# Get tiles +[FIND_EXPOSURES_RUNNER] + +INPUT_DIR = $SP_RUN/output/run_sp_tile_Git/get_images_runner/output + +FILE_PATTERN = CFIS_image + +FILE_EXT = .fits + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +# Column number of exposure name in FITS header +COLNUM = 3 + +# Prefix to remove from exposure name +EXP_PREFIX = p + diff --git a/workflow/config/cfis/config_tile_Git.ini b/workflow/config/cfis/config_tile_Git.ini new file mode 100644 index 000000000..72a0f0be7 --- /dev/null +++ b/workflow/config/cfis/config_tile_Git.ini @@ -0,0 +1,93 @@ +# ShapePipe configuration file for: get tile images + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = False + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Git + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = get_images_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names +INPUT_DIR = $SP_RUN + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +# Get tiles +[GET_IMAGES_RUNNER] + +FILE_PATTERN = tile_numbers + +FILE_EXT = .txt + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = + +# Paths + +# Input path where original images are stored. Can be local path or vos url. +# Single string or list of strings +INPUT_PATH = /project/def-mjhudson/unions-wl/tiles, /project/def-mjhudson/unions-wl/tiles + +# Input file pattern including tile number as dummy template +INPUT_FILE_PATTERN = CFIS.000.000.r, CFIS.000.000.r.weight + +# Input file extensions +INPUT_FILE_EXT = .fits, .fits.fz + +# Input numbering scheme, python regexp +INPUT_NUMBERING = \d{3}\.\d{3} + +# Output file pattern without number +OUTPUT_FILE_PATTERN = CFIS_image-, CFIS_weight- + +# Copy/download method, one in 'vos', 'symlink' +RETRIEVE = symlink + +# If RETRIEVE=vos, number of attempts to download +# Optional, default=3 +N_TRY = 3 + +# Copy command options, optional +RETRIEVE_OPTIONS = --certfile=$HOME/.ssl/cadcproxy.pem + +#CHECK_EXISTING_DIR = $SP_RUN/data_tiles diff --git a/workflow/config/cfis/config_tile_Mc.ini b/workflow/config/cfis/config_tile_Mc.ini new file mode 100644 index 000000000..502731431 --- /dev/null +++ b/workflow/config/cfis/config_tile_Mc.ini @@ -0,0 +1,80 @@ +# ShapePipe post-run configuration file: create final catalogs, with +# no spread model on input + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_Mc + +# Add date and time to RUN_NAME, optional, default: True +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = make_cat_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = ./output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 8 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[MAKE_CAT_RUNNER] + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = $SP_RUN/output/run_sp_tile_Sx/sextractor_runner/output, $SP_RUN/output/run_sp_tile_PiViVi/psfex_interp_runner/output, $SP_RUN/output/run_sp_Ms/merge_sep_cats_runner/output + +# Input file pattern(s), list of strings with length matching number of expected input file types +# Cannot contain wild cards +FILE_PATTERN = sexcat, galaxy_psf, ngmix + +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed tile ID (e.g. -210-282) +NUMBER_LIST = $SP_UNIT_NUM + +# FILE_EXT (optional) list of string extensions to identify input files +FILE_EXT = .fits, .sqlite, .fits + +# Numbering convention, string that exemplifies a numbering pattern. +# Matches input single exposures (with 'p' removed) +# Needs to be given in this section, will be updated in module +# sections below +NUMBERING_SCHEME = -000-000 + +SM_DO_CLASSIFICATION = False + +SHAPE_MEASUREMENT_TYPE = ngmix diff --git a/workflow/config/cfis/config_tile_Mh_exp.ini b/workflow/config/cfis/config_tile_Mh_exp.ini new file mode 100644 index 000000000..96512df1a --- /dev/null +++ b/workflow/config/cfis/config_tile_Mh_exp.ini @@ -0,0 +1,76 @@ +# ShapePipe configuration file for merging per-exposure WCS headers +# at the tile level. Input is the exp_numbers file produced by +# find_exposures_runner; EXP_BASE_DIR tells the runner where to find +# the per-exposure split_exp_runner header .npy files. + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Mh_exp + +# Add date and time to RUN_NAME, optional, default: True +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = merge_headers_runner + +# Run mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed tile ID (e.g. -210-282). +NUMBER_LIST = $SP_UNIT_NUM + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 16 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[MERGE_HEADERS_RUNNER] + +# Input: exp_numbers txt file from find_exposures_runner +INPUT_DIR = $SP_RUN/output/run_sp_tile_Fe/find_exposures_runner/output + +FILE_PATTERN = exp_numbers + +FILE_EXT = .txt + +# Tile numbering scheme (RA-Dec, e.g. -301-279) +NUMBERING_SCHEME = -000-000 + +# Root directory containing all per-exposure work directories. +# The runner will walk this tree to collect headers-.npy files. +EXP_BASE_DIR = $SP_EXP diff --git a/workflow/config/cfis/config_tile_Ng_template.ini b/workflow/config/cfis/config_tile_Ng_template.ini new file mode 100644 index 000000000..6e3e33292 --- /dev/null +++ b/workflow/config/cfis/config_tile_Ng_template.ini @@ -0,0 +1,104 @@ +# ShapePipe configuration file for tiles: ngmix + KSB + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +# Per-chunk run dir: the workflow exports SP_NGMIX_CHUNK= for each chunk. +# RUN_NAME is env-expanded (run.py getexpanded); braces keep the trailing "u" +# out of the variable name. +RUN_NAME = run_sp_tile_ngmix_Ng${SP_NGMIX_CHUNK}u + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = ngmix_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed tile ID (e.g. -210-282). +NUMBER_LIST = $SP_UNIT_NUM + +# Input directory, containing input files, single string or list of names +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +# Model-fitting shapes with ngmix +[NGMIX_RUNNER] + +INPUT_DIR = $SP_RUN/output/run_sp_tile_Sx/sextractor_runner/output, $SP_RUN/output/run_sp_tile_PiViVi/psfex_interp_runner/output, $SP_RUN/output/run_sp_tile_PiViVi/vignetmaker_runner_run_2/output, $SP_RUN/output/run_sp_tile_Mh_exp/merge_headers_runner/output + +FILE_PATTERN = sexcat, image_vignet, background_vignet, galaxy_psf, weight_vignet, flag_vignet, log_exp_headers + +FILE_EXT = .fits, .sqlite, .sqlite, .sqlite, .sqlite, .sqlite, .sqlite + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +# BKG_RMS_VIGNET_PATH (optional): per-pixel BACKGROUND_RMS vignets, used as +# 1/RMS^2 inverse-variance ngmix weights. When set, the file must exist for +# every tile (missing file -> error, no per-tile fallback); omit the option +# entirely to fall back to the scalar sigma_mad noise estimate. +BKG_RMS_VIGNET_PATH = $SP_RUN/output/run_sp_tile_PiViVi/vignetmaker_runner_run_2/output/background_rms_vignet{file_number_string}.sqlite + +# Number of objects to batch save during processing, optional. Omit or set +# to -1 for no batch saving. +# 250 (folded from runs/p3-batch1/cfis): worker RSS grows ~4.8 MB/object +# until the flush recycles it -> peak ~1.1 GB + 250 x 4.8 MB ~= 2.3 GB/worker +# (A/B test, job 17607877). Not superseded by -b {threads} (that sets fork +# width, this bounds per-worker memory). +SAVE_BATCH = 250 + +# Magnitude zero-point +MAG_ZP = 30.0 + +# Pixel scale in arcsec +PIXEL_SCALE = 0.186 + +# SEED_FROM_POSITION: per-object RNG seeded from sky position (ra, dec, ccd) +# instead of one ordered per-tile stream, so results are bit-identical under +# any chunking (D4/#796). Required for the chunked ngmix scatter/gather. +SEED_FROM_POSITION = True + +# ID_OBJ_MIN/MAX: this chunk's closed SExtractor NUMBER-column range, +# computed at execution time from the tile's own object count and expanded +# via ShapePipe's getexpanded (ngmix_runner.py verified: env-expanded, not +# plain getint). +ID_OBJ_MIN = $NGMIX_ID_MIN +ID_OBJ_MAX = $NGMIX_ID_MAX diff --git a/workflow/config/cfis/config_tile_PiViVi.ini b/workflow/config/cfis/config_tile_PiViVi.ini new file mode 100644 index 000000000..22294213f --- /dev/null +++ b/workflow/config/cfis/config_tile_PiViVi.ini @@ -0,0 +1,172 @@ +# ShapePipe configuration file for tile, from detection up to shape measurement. +# PSFEx PSF model. + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_PiViVi + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +#MODULE = psfex_interp_runner, + +MODULE = psfex_interp_runner, vignetmaker_runner, vignetmaker_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed tile ID (e.g. -210-282). +NUMBER_LIST = $SP_UNIT_NUM + +# Input directory, containing input files, single string or list of names +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 16 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[PSFEX_INTERP_RUNNER] + +INPUT_DIR = $SP_RUN/output/run_sp_tile_Sx/sextractor_runner/output, $SP_RUN/output/run_sp_tile_Mh_exp/merge_headers_runner/output, $SP_RUN/output/run_sp_tile_Fe/find_exposures_runner/output + +FILE_PATTERN = sexcat, log_exp_headers, exp_numbers + +FILE_EXT = .fits, .sqlite, .txt + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +# Run mode for psfex interpolation: +# CLASSIC: 'classical' run, interpolate to object positions +# MULTI-EPOCH: interpolate for multi-epoch images +# VALIDATION: validation for single-epoch images +MODE = MULTI-EPOCH + +# Column names of position parameters +POSITION_PARAMS = XWIN_WORLD,YWIN_WORLD + +# If True, measure and store ellipticity of the PSF +GET_SHAPES = True + +# Number of stars threshold +STAR_THRESH = 22 + +# chi^2 threshold +CHI2_THRESH = 2 + +# Multi-epoch mode parameters + +# Root directory of per-exposure work directories; replaces ME_DOT_PSF_DIR +# for v2.0 per-exposure pipeline. psfex_runner/output/ dirs are discovered +# by scanning $SP_EXP for the exposures listed in the exp_numbers input file. +ME_DOT_PSF_EXP_DIR = $SP_EXP + +# Input psf file pattern +ME_DOT_PSF_PATTERN = star_split_ratio_80 + + +# Create vignets for tiles weights +[VIGNETMAKER_RUNNER_RUN_1] + +INPUT_DIR = $SP_RUN/output/run_sp_tile_Sx/sextractor_runner/output, $SP_RUN/output/run_sp_tile_Uz/uncompress_fits_runner/output + +FILE_PATTERN = sexcat, CFIS_weight + +FILE_EXT = .fits, .fits + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +MASKING = False +MASK_VALUE = 0 + +# Run mode for psfex interpolation: +# CLASSIC: 'classical' run, interpolate to object positions +# MULTI-EPOCH: interpolate for multi-epoch images +# VALIDATION: validation for single-epoch images +MODE = CLASSIC + +# Coordinate frame type, one in PIX (pixel frame), SPHE (spherical coordinates) +COORD = PIX +POSITION_PARAMS = XWIN_IMAGE,YWIN_IMAGE + +# Vignet size in pixels +STAMP_SIZE = 51 + +# Output file name prefix, file name is _vignet.fits +PREFIX = weight + + +[VIGNETMAKER_RUNNER_RUN_2] + +# Create multi-epoch vignets for tiles corresponding to +# positions on single-exposures + +INPUT_DIR = $SP_RUN/output/run_sp_tile_Sx/sextractor_runner/output, $SP_RUN/output/run_sp_tile_Mh_exp/merge_headers_runner/output, $SP_RUN/output/run_sp_tile_Fe/find_exposures_runner/output + +FILE_PATTERN = sexcat, log_exp_headers, exp_numbers + +FILE_EXT = .fits, .sqlite, .txt + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +MASKING = False +MASK_VALUE = 0 + +# Run mode for psfex interpolation: +# CLASSIC: 'classical' run, interpolate to object positions +# MULTI-EPOCH: interpolate for multi-epoch images +# VALIDATION: validation for single-epoch images +MODE = MULTI-EPOCH + +# Coordinate frame type, one in PIX (pixel frame), SPHE (spherical coordinates) +COORD = SPHE +POSITION_PARAMS = XWIN_WORLD,YWIN_WORLD + +# Vignet size in pixels +STAMP_SIZE = 51 + +# Output file name prefix, file name is vignet.fits +PREFIX = + +# Additional parameters for path and file pattern corresponding to single-exposure +# run outputs. ME_IMAGE_EXP_DIR/ME_IMAGE_EXP_RUNNERS replace ME_IMAGE_DIR for +# the v2.0 per-exposure pipeline; output dirs are discovered by scanning $SP_EXP. +ME_IMAGE_EXP_DIR = $SP_EXP +ME_IMAGE_EXP_RUNNERS = split_exp_runner, split_exp_runner, split_exp_runner, sextractor_runner, sextractor_runner +ME_IMAGE_PATTERN = flag, image, weight, background, background_rms diff --git a/workflow/config/cfis/config_tile_Sx.ini b/workflow/config/cfis/config_tile_Sx.ini new file mode 100644 index 000000000..0ce9ea226 --- /dev/null +++ b/workflow/config/cfis/config_tile_Sx.ini @@ -0,0 +1,118 @@ +# ShapePipe configuration file for tile detection + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Sx + +# Add date and time to RUN_NAME, optional, default: True +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = sextractor_runner + + +# Run mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed tile ID (e.g. -210-282). +NUMBER_LIST = $SP_UNIT_NUM + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = $SP_RUN/output + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 16 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[SEXTRACTOR_RUNNER] + +INPUT_DIR = $SP_RUN/output/run_sp_tile_Git/get_images_runner/output, $SP_RUN/output/run_sp_tile_Uz/uncompress_fits_runner/output, $SP_RUN/output/run_sp_tile_Mh_exp/merge_headers_runner/output + +FILE_PATTERN = CFIS_image, CFIS_weight, log_exp_headers + +FILE_EXT = .fits, .fits, .sqlite + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +# SExtractor executable path +EXEC_PATH = source-extractor + +# SExtractor configuration files +DOT_SEX_FILE = $SP_CONFIG/default_tile.sex +DOT_PARAM_FILE = $SP_CONFIG/default_noimaflags.param +DOT_CONV_FILE = $SP_CONFIG/default.conv + +# Use input weight image if True +WEIGHT_IMAGE = True + +# Use input flag image if True +FLAG_IMAGE = False + +# Use input PSF file if True +PSF_FILE = False + +# Use distinct image for detection (SExtractor in +# dual-image mode) if True +DETECTION_IMAGE = False + +# Distinct weight image for detection (SExtractor +# in dual-image mode) +DETECTION_WEIGHT = False + +ZP_FROM_HEADER = False + +BKG_FROM_HEADER = False + +# Type of image check (optional), default not used, can be a list of +# BACKGROUND, BACKGROUND_RMS, INIBACKGROUND, +# MINIBACK_RMS, -BACKGROUND, #FILTERED, +# OBJECTS, -OBJECTS, SEGMENTATION, APERTURES +CHECKIMAGE = BACKGROUND + +# File name suffix for the output sextractor files (optional) +SUFFIX = sexcat + +## Post-processing + +# Necessary for tiles, to enable multi-exposure processing +MAKE_POST_PROCESS = True + +# World coordinate keywords, SExtractor output. Format: KEY_X,KEY_Y +WORLD_POSITION = XWIN_WORLD,YWIN_WORLD + +# Number of pixels in x,y of a CCD. Format: Nx,Ny +CCD_SIZE = 33,2080,1,4612 diff --git a/workflow/config/cfis/config_tile_Uc.ini b/workflow/config/cfis/config_tile_Uc.ini new file mode 100644 index 000000000..f2ff51e4a --- /dev/null +++ b/workflow/config/cfis/config_tile_Uc.ini @@ -0,0 +1,92 @@ +# ShapePipe configuration file for tile object selection using +# an (external) catalogue + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Uc + +# Add date and time to RUN_NAME, optional, default: True +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = read_ext_sexcat_runner + + +# Run mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed tile ID (e.g. -210-282). +NUMBER_LIST = $SP_UNIT_NUM + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = $SP_RUN/output + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 16 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[READ_EXT_SEXCAT_RUNNER] + +# NOTE(blocker): run_sp_tile_Gic (external-catalog get_images) is not one of +# the 13 configs in this sweep -- no committed config produces it, so this +# path is written by naming-convention analogy, unverified. Uc is otherwise +# fully de-wrapper-ified; do not rely on this input until Gic exists. +INPUT_DIR = $SP_RUN/output/run_sp_tile_Gic/get_images_runner/output, $SP_RUN/output/run_sp_tile_Git/get_images_runner/output, $SP_RUN/output/run_sp_tile_Mh_exp/merge_headers_runner/output + +FILE_PATTERN = CFIS_cat, CFIS_image, log_exp_headers + +FILE_EXT = .cat, .fits, .sqlite + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +# File name suffix for the output sextractor files (optional) +SUFFIX = sexcat + +# Side length of the square postage stamp (vignet) extracted from the tile +# image, in pixels (must be odd). Default: 51 +VIGNET_SIZE = 51 + +## Post-processing + +# Necessary for tiles, to enable multi-exposure processing +MAKE_POST_PROCESS = True + +# World coordinate keywords, SExtractor output. Format: KEY_X,KEY_Y +WORLD_POSITION = ALPHA_J2000,DELTA_J2000 + +# Number of pixels in x,y of a CCD. Format: Nx,Ny +CCD_SIZE = 33,2080,1,4612 diff --git a/workflow/config/cfis/config_tile_Uz.ini b/workflow/config/cfis/config_tile_Uz.ini new file mode 100644 index 000000000..fc8550af4 --- /dev/null +++ b/workflow/config/cfis/config_tile_Uz.ini @@ -0,0 +1,73 @@ +# ShapePipe configuration file for: uncompress FITS image + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Uz + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = uncompress_fits_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the +# dashed tile ID (e.g. -210-282). +NUMBER_LIST = $SP_UNIT_NUM + +# Input directory, containing input files, single string or list of names +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 16 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options +[UNCOMPRESS_FITS_RUNNER] + +INPUT_DIR = $SP_RUN/output/run_sp_tile_Git/get_images_runner/output + +FILE_PATTERN = CFIS_weight + +FILE_EXT = .fitsfz + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +# Input HDU of image data, optional, default=0 +HDU_DATA = 1 + +# Output file pattern +OUTPUT_PATTERN = CFIS_weight diff --git a/workflow/config/cfis/config_tile_onthefly.mask b/workflow/config/cfis/config_tile_onthefly.mask new file mode 100644 index 000000000..18cf5db2d --- /dev/null +++ b/workflow/config/cfis/config_tile_onthefly.mask @@ -0,0 +1,90 @@ +# Mask module config file for tiles + +## Paths to executables +[PROGRAM_PATH] + +WW_PATH = weightwatcher +WW_CONFIG_FILE = $SP_CONFIG/mask_default/default.ww + +# Indicate cds client executable if no external star catalogue is available +# (e.g. no internet access on run nodes) +CDSCLIENT_PATH = findgsc2.2 + +## Border parameters +[BORDER_PARAMETERS] + +BORDER_MAKE = False + +BORDER_WIDTH = 0 +BORDER_FLAG_VALUE = 4 + + +## Halo parameters +[HALO_PARAMETERS] + +HALO_MAKE = True + +HALO_MASKMODEL_PATH = $SP_CONFIG/mask_default/halo_mask.reg +HALO_MAG_LIM = 13. +HALO_SCALE_FACTOR = 0.05 +HALO_MAG_PIVOT = 13.8 +HALO_FLAG_VALUE = 2 +HALO_REG_FILE = halo.reg + + +## Diffraction pike parameters +[SPIKE_PARAMETERS] + +SPIKE_MAKE = True + +SPIKE_MASKMODEL_PATH = $SP_CONFIG/mask_default/MEGAPRIME_star_i_13.8.reg +SPIKE_MAG_LIM = 18. +SPIKE_SCALE_FACTOR = 0.3 +SPIKE_MAG_PIVOT = 13.8 +SPIKE_FLAG_VALUE = 128 +SPIKE_REG_FILE = spike.reg + + +## Messier parameters +[MESSIER_PARAMETERS] + +MESSIER_MAKE = True + +MESSIER_CAT_PATH = $SP_CONFIG/mask_default/Messier_catalog_updated.fits +MESSIER_PIXEL_SCALE = 0.187 +MESSIER_SIZE_PLUS = 0. +MESSIER_FLAG_VALUE = 16 + +## NGC mask +[NGC_PARAMETERS] + +NGC_MAKE = True + +NGC_CAT_PATH = $SP_CONFIG/mask_default/ngc_cat.fits +NGC_SIZE_PLUS = 0. +NGC_FLAG_VALUE = 32 + + +## External flag +[EXTERNAL_FLAG] + +EF_MAKE = False + + +## Missing data parameters +[MD_PARAMETERS] + +MD_MAKE = False + +MD_THRESH_FLAG = 0.3 +MD_THRESH_REMOVE = 0.75 +MD_REMOVE = False + + +## Other parameters +[OTHER] + +KEEP_REG_FILE = False +KEEP_INDIVIDUAL_MASK = False + +TEMP_DIRECTORY = .temp_tiles diff --git a/workflow/config/cfis/default.conv b/workflow/config/cfis/default.conv new file mode 100644 index 000000000..2590b9cba --- /dev/null +++ b/workflow/config/cfis/default.conv @@ -0,0 +1,5 @@ +CONV NORM +# 3x3 ``all-ground'' convolution mask with FWHM = 2 pixels. +1 2 1 +2 4 2 +1 2 1 diff --git a/workflow/config/cfis/default.param b/workflow/config/cfis/default.param new file mode 100644 index 000000000..09ad8405e --- /dev/null +++ b/workflow/config/cfis/default.param @@ -0,0 +1,68 @@ +NUMBER #Running object number +EXT_NUMBER #FITS extension number + +FLUX_AUTO #Flux within a Kron-like elliptical aperture [count] +FLUXERR_AUTO #RMS error for AUTO flux [count] +MAG_AUTO #Kron-like elliptical aperture magnitude [mag] +MAGERR_AUTO #RMS error for AUTO magnitude [mag] +FLUX_WIN #Gaussian-weighted flux [count] +FLUXERR_WIN #RMS error for WIN flux [count] +MAG_WIN #Gaussian-weighted magnitude [mag] +MAGERR_WIN #RMS error for MAG_WIN [mag] +FLUX_APER(1) +FLUXERR_APER(1) + +FLUX_RADIUS #Fraction-of-light radii [pixel] + +SNR_WIN #Gaussian-weighted SNR + +BACKGROUND #Background at centroid position [count] +THRESHOLD #Detection threshold above background [count] + +X_IMAGE #Object position along x [pixel] +Y_IMAGE #Object position along y [pixel] + +X_WORLD #Barycenter position along world x axis [deg] +Y_WORLD #Barycenter position along world y axis [deg] + +X2_IMAGE #Variance along x [pixel**2] +Y2_IMAGE #Variance along y [pixel**2] +XY_IMAGE #Covariance between x and y [pixel**2] +ERRX2_IMAGE #Variance of position along x [pixel**2] +ERRY2_IMAGE #Variance of position along y [pixel**2] +ERRXY_IMAGE #Covariance of position between x and y [pixel**2] + +XWIN_IMAGE #Windowed position estimate along x [pixel] +YWIN_IMAGE #Windowed position estimate along y [pixel] + +XWIN_WORLD #Windowed position along world x axis [deg] +YWIN_WORLD #Windowed position along world y axis [deg] + +X2WIN_IMAGE #Windowed variance along x [pixel**2] +Y2WIN_IMAGE #Windowed variance along y [pixel**2] +XYWIN_IMAGE #Windowed covariance between x and y [pixel**2] +ERRX2WIN_IMAGE #Variance of windowed pos along x [pixel**2] +ERRY2WIN_IMAGE #Variance of windowed pos along y [pixel**2] +ERRXYWIN_IMAGE #Covariance of windowed pos between x and y [pixel**2] + +MU_THRESHOLD #Analysis threshold above background [mag * arcsec**(-2)] +MU_MAX #Peak surface brightness above background [mag * arcsec**(-2)] + +FLAGS #Extraction flags +FLAGS_WIN #Flags for WINdowed parameters + +# The following flag requires a flag image +IMAFLAGS_ISO #FLAG-image flags OR'ed over the iso. profile !!! REQUIRE FLAG_IMAGE !!! + +FWHM_IMAGE #FWHM assuming a gaussian core [pixel] +FWHM_WORLD #FWHM assuming a gaussian core [deg] +ELONGATION #A_IMAGE/B_IMAGE +ELLIPTICITY #1 - B_IMAGE/A_IMAGE + +VIGNET(51,51) #Pixel data around detection [count] + +# For GaaP photometry +A_WORLD +B_WORLD +THETA_J2000 + diff --git a/workflow/config/cfis/default.psfex b/workflow/config/cfis/default.psfex new file mode 100644 index 000000000..a9d1a906c --- /dev/null +++ b/workflow/config/cfis/default.psfex @@ -0,0 +1,85 @@ +# Default configuration file for PSFEx 3.17.1 +# EB 2017-11-30 +# + +#-------------------------------- PSF model ---------------------------------- + +BASIS_TYPE PIXEL # NONE, PIXEL, GAUSS-LAGUERRE or FILE +BASIS_NUMBER 20 # Basis number or parameter +BASIS_NAME basis.fits # Basis filename (FITS data-cube) +BASIS_SCALE 1.0 # Gauss-Laguerre beta parameter +NEWBASIS_TYPE NONE # Create new basis: NONE, PCA_INDEPENDENT + # or PCA_COMMON +NEWBASIS_NUMBER 8 # Number of new basis vectors +PSF_SAMPLING 1. # Sampling step in pixel units (0.0 = auto) +PSF_PIXELSIZE 1.0 # Effective pixel size in pixel step units +PSF_ACCURACY 0.01 # Accuracy to expect from PSF "pixel" values +PSF_SIZE 51,51 # Image size of the PSF model +PSF_RECENTER N # Allow recentering of PSF-candidates Y/N ? +MEF_TYPE INDEPENDENT # INDEPENDENT or COMMON + +#------------------------- Point source measurements ------------------------- + +CENTER_KEYS XWIN_IMAGE,YWIN_IMAGE # Catalogue parameters for source pre-centering +PHOTFLUX_KEY FLUX_AUTO # Catalogue parameter for photometric norm. +PHOTFLUXERR_KEY FLUXERR_AUTO # Catalogue parameter for photometric error + +#----------------------------- PSF variability ------------------------------- + +PSFVAR_KEYS XWIN_IMAGE,YWIN_IMAGE # Catalogue or FITS (preceded by :) params +PSFVAR_GROUPS 1,1 # Group tag for each context key +PSFVAR_DEGREES 2 # Polynom degree for each group +PSFVAR_NSNAP 9 # Number of PSF snapshots per axis +HIDDENMEF_TYPE COMMON # INDEPENDENT or COMMON +STABILITY_TYPE EXPOSURE # EXPOSURE or SEQUENCE + +#----------------------------- Sample selection ------------------------------ + +SAMPLE_AUTOSELECT N # Automatically select the FWHM (Y/N) ? + +BADPIXEL_FILTER N # Filter bad-pixels in samples (Y/N) ? +BADPIXEL_NMAX 0 # Maximum number of bad pixels allowed + +#----------------------- PSF homogeneisation kernel -------------------------- + +HOMOBASIS_TYPE NONE # NONE or GAUSS-LAGUERRE +HOMOBASIS_NUMBER 10 # Kernel basis number or parameter +HOMOBASIS_SCALE 1.0 # GAUSS-LAGUERRE beta parameter +HOMOPSF_PARAMS 2.0, 3.0 # Moffat parameters of the idealised PSF +HOMOKERNEL_DIR # Where to write kernels (empty=same as input) +HOMOKERNEL_SUFFIX .homo.fits # Filename extension for homogenisation kernels + +#----------------------------- Output catalogs ------------------------------- + +OUTCAT_TYPE FITS_LDAC # NONE, ASCII_HEAD, ASCII, FITS_LDAC + +#------------------------------- Check-plots ---------------------------------- + +CHECKPLOT_DEV NULL # NULL, XWIN, TK, PS, PSC, XFIG, PNG, + # JPEG, AQT, PDF or SVG +CHECKPLOT_RES 0 # Check-plot resolution (0 = default) +CHECKPLOT_ANTIALIAS Y # Anti-aliasing using convert (Y/N) ? +CHECKPLOT_TYPE NONE # FWHM,ELLIPTICITY,COUNTS, COUNT_FRACTION, CHI2, RESIDUALS +CHECKPLOT_TYPE FWHM,ELLIPTICITY,COUNTS, COUNT_FRACTION, CHI2, RESIDUALS + # or NONE +CHECKPLOT_NAME fwhm, ellipticity, counts, countfrac, chi2, resi + +#------------------------------ Check-Images --------------------------------- + +# Note: Check-image types can be set the ShapePipe config file, psfex_runner section +####### +CHECKIMAGE_TYPE NONE # CHI,PROTOTYPES,SAMPLES,RESIDUALS,SNAPSHOTS + # or MOFFAT,-MOFFAT,-SYMMETRICAL +#CHECKIMAGE_NAME chi.fits,proto.fits,samp.fits,resi.fits,snap.fits + # Check-image filenames +#CHECKIMAGE_CUBE N # Save check-images as datacubes (Y/N) ? + +#----------------------------- Miscellaneous --------------------------------- + +PSF_SUFFIX .psf # Filename extension for output PSF filename +VERBOSE_TYPE NORMAL # can be QUIET,NORMAL,LOG or FULL +WRITE_XML N # Write XML file (Y/N)? + +NTHREADS 1 # Number of simultaneous threads for + # the SMP version of PSFEx + # 0 = automatic diff --git a/workflow/config/cfis/default_exp.sex b/workflow/config/cfis/default_exp.sex new file mode 100644 index 000000000..b87275ecb --- /dev/null +++ b/workflow/config/cfis/default_exp.sex @@ -0,0 +1,133 @@ +# Default configuration file for SExtractor 2.19.5 +# EB 2017-11-30 +# + +#-------------------------------- Catalog ------------------------------------ + +CATALOG_TYPE FITS_LDAC + +PARAMETERS_NAME default.param + +#------------------------------- Extraction ---------------------------------- + +DETECT_TYPE CCD # CCD (linear) or PHOTO (with gamma correction) +DETECT_MINAREA 5 # min. # of pixels above threshold +DETECT_MAXAREA 0 # max. # of pixels above threshold (0=unlimited) +THRESH_TYPE RELATIVE # threshold type: RELATIVE (in sigmas) + # or ABSOLUTE (in ADUs) +DETECT_THRESH 1.5 # or , in mag.arcsec-2 +ANALYSIS_THRESH 1.5 # or , in mag.arcsec-2 + +FILTER Y # apply filter for detection (Y or N)? +FILTER_NAME default.conv +FILTER_THRESH # Threshold[s] for retina filtering + +DEBLEND_NTHRESH 32 # Number of deblending sub-thresholds +DEBLEND_MINCONT 0.001 # Minimum contrast parameter for deblending + +CLEAN Y # Clean spurious detections? (Y or N)? +CLEAN_PARAM 1.0 # Cleaning efficiency + +MASK_TYPE CORRECT # type of detection MASKing: can be one of + # NONE, BLANK or CORRECT + +#-------------------------------- WEIGHTing ---------------------------------- + +WEIGHT_TYPE MAP_WEIGHT # type of WEIGHTing: NONE, BACKGROUND, + # MAP_RMS, MAP_VAR or MAP_WEIGHT +RESCALE_WEIGHTS Y # Rescale input weights/variances (Y/N)? +WEIGHT_IMAGE weight.fits # weight-map filename +WEIGHT_GAIN Y # modulate gain (E/ADU) with weights? (Y/N) +WEIGHT_THRESH # weight threshold[s] for bad pixels + +#-------------------------------- FLAGging ----------------------------------- + +FLAG_IMAGE flag.fits # filename for an input FLAG-image +FLAG_TYPE OR # flag pixel combination: OR, AND, MIN, MAX + # or MOST + +#------------------------------ Photometry ----------------------------------- + +PHOT_APERTURES 5 # MAG_APER aperture diameter(s) in pixels +PHOT_AUTOPARAMS 2.5, 3.5 # MAG_AUTO parameters: , +PHOT_PETROPARAMS 2.0, 3.5 # MAG_PETRO parameters: , + # +PHOT_AUTOAPERS 0.0,0.0 # , minimum apertures + # for MAG_AUTO and MAG_PETRO +PHOT_FLUXFRAC 0.5 # flux fraction[s] used for FLUX_RADIUS + +SATUR_KEY SATURATE # keyword for saturation level (in ADUs) + +MAG_ZEROPOINT 30.0 # magnitude zero-point +MAG_GAMMA 4.0 # gamma of emulsion (for photographic scans) + +GAIN_KEY GAIN # keyword for detector gain in e-/ADU +PIXEL_SCALE 0. # size of pixel in arcsec (0=use FITS WCS info) + +#------------------------- Star/Galaxy Separation ---------------------------- + +SEEING_FWHM 0.6 # stellar FWHM in arcsec +STARNNW_NAME default.nnw + +#------------------------------ Background ----------------------------------- + +BACK_TYPE AUTO # AUTO or MANUAL +BACK_VALUE 0.0 # Default background value in MANUAL mode +BACK_SIZE 64 # Background mesh: or , +BACK_FILTERSIZE 3 # Background filter: or , + +BACKPHOTO_TYPE GLOBAL # can be GLOBAL or LOCAL +BACKPHOTO_THICK 24 # thickness of the background LOCAL annulus +BACK_FILTTHRESH 0.0 # Threshold above which the background- + # map filter operates + +#------------------------------ Check Image ---------------------------------- + +####### +## AG : This parameter is set in pipeline config file. +####### +# CHECKIMAGE_TYPE NONE #BACKGROUND_RMS,BACKGROUND +# can be NONE, BACKGROUND, BACKGROUND_RMS, + # MINIBACKGROUND, MINIBACK_RMS, -BACKGROUND, + # FILTERED, OBJECTS, -OBJECTS, SEGMENTATION, + # or APERTURES +# CHECKIMAGE_NAME check.fits,back.fits +# Filename for the check-image + +#--------------------- Memory (change with caution!) ------------------------- + +MEMORY_OBJSTACK 3000 # number of objects in stack +MEMORY_PIXSTACK 300000 # number of pixels in stack +MEMORY_BUFSIZE 1024 # number of lines in buffer + +#------------------------------- ASSOCiation --------------------------------- + +ASSOC_NAME sky.list # name of the ASCII file to ASSOCiate +ASSOC_DATA 2,3,4 # columns of the data to replicate (0=all) +ASSOC_PARAMS 2,3,4 # columns of xpos,ypos[,mag] +ASSOCCOORD_TYPE PIXEL # ASSOC coordinates: PIXEL or WORLD +ASSOC_RADIUS 2.0 # cross-matching radius (pixels) +ASSOC_TYPE NEAREST # ASSOCiation method: FIRST, NEAREST, MEAN, + # MAG_MEAN, SUM, MAG_SUM, MIN or MAX +ASSOCSELEC_TYPE MATCHED # ASSOC selection type: ALL, MATCHED or -MATCHED + +#----------------------------- Miscellaneous --------------------------------- + +VERBOSE_TYPE NORMAL # can be QUIET, NORMAL or FULL +HEADER_SUFFIX .head # Filename extension for additional headers +WRITE_XML N # Write XML file (Y/N)? + +NTHREADS 1 # 1 single thread + +FITS_UNSIGNED N # Treat FITS integer values as unsigned (Y/N)? +INTERP_MAXXLAG 16 # Max. lag along X for 0-weight interpolation +INTERP_MAXYLAG 16 # Max. lag along Y for 0-weight interpolation +INTERP_TYPE ALL # Interpolation type: NONE, VAR_ONLY or ALL + +#--------------------------- Experimental Stuff ----------------------------- + +#PSF_NAME default.psf # File containing the PSF model +#PSF_NMAX 1 # Max.number of PSFs fitted simultaneously +#PATTERN_TYPE RINGS-HARMONIC # can RINGS-QUADPOLE, RINGS-OCTOPOLE, + # RINGS-HARMONICS or GAUSS-LAGUERRE +#SOM_NAME default.som # File containing Self-Organizing Map weights diff --git a/workflow/config/cfis/default_noimaflags.param b/workflow/config/cfis/default_noimaflags.param new file mode 100644 index 000000000..b251f5b74 --- /dev/null +++ b/workflow/config/cfis/default_noimaflags.param @@ -0,0 +1,65 @@ +NUMBER #Running object number +EXT_NUMBER #FITS extension number + +FLUX_AUTO #Flux within a Kron-like elliptical aperture [count] +FLUXERR_AUTO #RMS error for AUTO flux [count] +MAG_AUTO #Kron-like elliptical aperture magnitude [mag] +MAGERR_AUTO #RMS error for AUTO magnitude [mag] +FLUX_WIN #Gaussian-weighted flux [count] +FLUXERR_WIN #RMS error for WIN flux [count] +MAG_WIN #Gaussian-weighted magnitude [mag] +MAGERR_WIN #RMS error for MAG_WIN [mag] +FLUX_APER(1) +FLUXERR_APER(1) + +FLUX_RADIUS #Fraction-of-light radii [pixel] + +SNR_WIN #Gaussian-weighted SNR + +BACKGROUND #Background at centroid position [count] +THRESHOLD #Detection threshold above background [count] + +X_IMAGE #Object position along x [pixel] +Y_IMAGE #Object position along y [pixel] + +X_WORLD #Barycenter position along world x axis [deg] +Y_WORLD #Barycenter position along world y axis [deg] + +X2_IMAGE #Variance along x [pixel**2] +Y2_IMAGE #Variance along y [pixel**2] +XY_IMAGE #Covariance between x and y [pixel**2] +ERRX2_IMAGE #Variance of position along x [pixel**2] +ERRY2_IMAGE #Variance of position along y [pixel**2] +ERRXY_IMAGE #Covariance of position between x and y [pixel**2] + +XWIN_IMAGE #Windowed position estimate along x [pixel] +YWIN_IMAGE #Windowed position estimate along y [pixel] + +XWIN_WORLD #Windowed position along world x axis [deg] +YWIN_WORLD #Windowed position along world y axis [deg] + +X2WIN_IMAGE #Windowed variance along x [pixel**2] +Y2WIN_IMAGE #Windowed variance along y [pixel**2] +XYWIN_IMAGE #Windowed covariance between x and y [pixel**2] +ERRX2WIN_IMAGE #Variance of windowed pos along x [pixel**2] +ERRY2WIN_IMAGE #Variance of windowed pos along y [pixel**2] +ERRXYWIN_IMAGE #Covariance of windowed pos between x and y [pixel**2] + +MU_THRESHOLD #Analysis threshold above background [mag * arcsec**(-2)] +MU_MAX #Peak surface brightness above background [mag * arcsec**(-2)] + +FLAGS #Extraction flags +FLAGS_WIN #Flags for WINdowed parameters + +FWHM_IMAGE #FWHM assuming a gaussian core [pixel] +FWHM_WORLD #FWHM assuming a gaussian core [deg] +ELONGATION #A_IMAGE/B_IMAGE +ELLIPTICITY #1 - B_IMAGE/A_IMAGE + +VIGNET(51,51) #Pixel data around detection [count] + +# For GaaP photometry +A_WORLD +B_WORLD +THETA_J2000 + diff --git a/workflow/config/cfis/default_tile.sex b/workflow/config/cfis/default_tile.sex new file mode 100644 index 000000000..ff3b25213 --- /dev/null +++ b/workflow/config/cfis/default_tile.sex @@ -0,0 +1,133 @@ +# Default configuration file for SExtractor 2.19.5 +# EB 2017-11-30 +# + +#-------------------------------- Catalog ------------------------------------ + +CATALOG_TYPE FITS_LDAC + +PARAMETERS_NAME default.param + +#------------------------------- Extraction ---------------------------------- + +DETECT_TYPE CCD # CCD (linear) or PHOTO (with gamma correction) +DETECT_MINAREA 5 # min. # of pixels above threshold +DETECT_MAXAREA 0 # max. # of pixels above threshold (0=unlimited) +THRESH_TYPE RELATIVE # threshold type: RELATIVE (in sigmas) + # or ABSOLUTE (in ADUs) +DETECT_THRESH 1.5 # or , in mag.arcsec-2 +ANALYSIS_THRESH 1.5 # or , in mag.arcsec-2 + +FILTER Y # apply filter for detection (Y or N)? +FILTER_NAME default.conv +FILTER_THRESH # Threshold[s] for retina filtering + +DEBLEND_NTHRESH 32 # Number of deblending sub-thresholds +DEBLEND_MINCONT 0.0005 # Minimum contrast parameter for deblending + +CLEAN Y # Clean spurious detections? (Y or N)? +CLEAN_PARAM 1.0 # Cleaning efficiency + +MASK_TYPE CORRECT # type of detection MASKing: can be one of + # NONE, BLANK or CORRECT + +#-------------------------------- WEIGHTing ---------------------------------- + +WEIGHT_TYPE MAP_WEIGHT # type of WEIGHTing: NONE, BACKGROUND, + # MAP_RMS, MAP_VAR or MAP_WEIGHT +RESCALE_WEIGHTS Y # Rescale input weights/variances (Y/N)? +WEIGHT_IMAGE weight.fits # weight-map filename +WEIGHT_GAIN Y # modulate gain (E/ADU) with weights? (Y/N) +WEIGHT_THRESH # weight threshold[s] for bad pixels + +#-------------------------------- FLAGging ----------------------------------- + +FLAG_IMAGE flag.fits # filename for an input FLAG-image +FLAG_TYPE OR # flag pixel combination: OR, AND, MIN, MAX + # or MOST + +#------------------------------ Photometry ----------------------------------- + +PHOT_APERTURES 5 # MAG_APER aperture diameter(s) in pixels +PHOT_AUTOPARAMS 2.5, 3.5 # MAG_AUTO parameters: , +PHOT_PETROPARAMS 2.0, 3.5 # MAG_PETRO parameters: , + # +PHOT_AUTOAPERS 0.0,0.0 # , minimum apertures + # for MAG_AUTO and MAG_PETRO +PHOT_FLUXFRAC 0.5 # flux fraction[s] used for FLUX_RADIUS + +SATUR_KEY SATURATE # keyword for saturation level (in ADUs) + +MAG_ZEROPOINT 30.0 # magnitude zero-point +MAG_GAMMA 4.0 # gamma of emulsion (for photographic scans) + +GAIN_KEY GAIN # keyword for detector gain in e-/ADU +PIXEL_SCALE 0. # size of pixel in arcsec (0=use FITS WCS info) + +#------------------------- Star/Galaxy Separation ---------------------------- + +SEEING_FWHM 0.6 # stellar FWHM in arcsec +STARNNW_NAME default.nnw + +#------------------------------ Background ----------------------------------- + +BACK_TYPE MANUAL # AUTO or MANUAL +BACK_VALUE 0.0 # Default background value in MANUAL mode +BACK_SIZE 64 # Background mesh: or , +BACK_FILTERSIZE 3 # Background filter: or , + +BACKPHOTO_TYPE GLOBAL # can be GLOBAL or LOCAL +BACKPHOTO_THICK 24 # thickness of the background LOCAL annulus +BACK_FILTTHRESH 0.0 # Threshold above which the background- + # map filter operates + +#------------------------------ Check Image ---------------------------------- + +####### +## AG : This parameter is set in pipeline config file. +####### +# CHECKIMAGE_TYPE NONE #BACKGROUND_RMS,BACKGROUND +# can be NONE, BACKGROUND, BACKGROUND_RMS, + # MINIBACKGROUND, MINIBACK_RMS, -BACKGROUND, + # FILTERED, OBJECTS, -OBJECTS, SEGMENTATION, + # or APERTURES +# CHECKIMAGE_NAME check.fits,back.fits +# Filename for the check-image + +#--------------------- Memory (change with caution!) ------------------------- + +MEMORY_OBJSTACK 3000 # number of objects in stack +MEMORY_PIXSTACK 300000 # number of pixels in stack +MEMORY_BUFSIZE 1024 # number of lines in buffer + +#------------------------------- ASSOCiation --------------------------------- + +ASSOC_NAME sky.list # name of the ASCII file to ASSOCiate +ASSOC_DATA 2,3,4 # columns of the data to replicate (0=all) +ASSOC_PARAMS 2,3,4 # columns of xpos,ypos[,mag] +ASSOCCOORD_TYPE PIXEL # ASSOC coordinates: PIXEL or WORLD +ASSOC_RADIUS 2.0 # cross-matching radius (pixels) +ASSOC_TYPE NEAREST # ASSOCiation method: FIRST, NEAREST, MEAN, + # MAG_MEAN, SUM, MAG_SUM, MIN or MAX +ASSOCSELEC_TYPE MATCHED # ASSOC selection type: ALL, MATCHED or -MATCHED + +#----------------------------- Miscellaneous --------------------------------- + +VERBOSE_TYPE NORMAL # can be QUIET, NORMAL or FULL +HEADER_SUFFIX .head # Filename extension for additional headers +WRITE_XML N # Write XML file (Y/N)? + +NTHREADS 1 # 1 single thread + +FITS_UNSIGNED N # Treat FITS integer values as unsigned (Y/N)? +INTERP_MAXXLAG 16 # Max. lag along X for 0-weight interpolation +INTERP_MAXYLAG 16 # Max. lag along Y for 0-weight interpolation +INTERP_TYPE ALL # Interpolation type: NONE, VAR_ONLY or ALL + +#--------------------------- Experimental Stuff ----------------------------- + +#PSF_NAME default.psf # File containing the PSF model +#PSF_NMAX 1 # Max.number of PSFs fitted simultaneously +#PATTERN_TYPE RINGS-HARMONIC # can RINGS-QUADPOLE, RINGS-OCTOPOLE, + # RINGS-HARMONICS or GAUSS-LAGUERRE +#SOM_NAME default.som # File containing Self-Organizing Map weights diff --git a/workflow/config/cfis/final_cat.param b/workflow/config/cfis/final_cat.param new file mode 100644 index 000000000..8ddfebd37 --- /dev/null +++ b/workflow/config/cfis/final_cat.param @@ -0,0 +1,114 @@ +# coordinates +XWIN_WORLD +YWIN_WORLD + +# tile ID, for plot of tile-dependent additive bias. +# Can maybe be removed. +TILE_ID + +# flags +FLAGS +IMAFLAGS_ISO +NGMIX_MCAL_FLAGS + +# PSF ellipticity (original image PSF) +NGMIX_G1_PSF_ORIG_NOSHEAR +NGMIX_G2_PSF_ORIG_NOSHEAR + +# spread class +#SPREAD_CLASS + +# spread model flag and error +#SPREAD_MODEL +#SPREADERR_MODEL + +# Number of epochs (exposures) +N_EPOCH +NGMIX_N_EPOCH + +## Shape measurement outputs +## Ngmix: model fitting + +# galaxy ellipticity +NGMIX_G1_1M +NGMIX_G2_1M +NGMIX_G1_1P +NGMIX_G2_1P +NGMIX_G1_2M +NGMIX_G2_2M +NGMIX_G1_2P +NGMIX_G2_2P +NGMIX_G1_NOSHEAR +NGMIX_G2_NOSHEAR +#NGMIX_G1_ERR_1M +#NGMIX_G2_ERR_1M +#NGMIX_G1_ERR_1P +#NGMIX_G2_ERR_1P +#NGMIX_G1_ERR_2M +#NGMIX_G2_ERR_2M +#NGMIX_G1_ERR_2P +#NGMIX_G2_ERR_2P +NGMIX_G1_ERR_NOSHEAR +NGMIX_G2_ERR_NOSHEAR + +# flags +NGMIX_FLAGS_1M +NGMIX_FLAGS_1P +NGMIX_FLAGS_2M +NGMIX_FLAGS_2P +NGMIX_FLAGS_NOSHEAR + +# size and error +NGMIX_T_1M +NGMIX_T_1P +NGMIX_T_2M +NGMIX_T_2P +NGMIX_T_NOSHEAR +NGMIX_T_ERR_1M +NGMIX_T_ERR_1P +NGMIX_T_ERR_2M +NGMIX_T_ERR_2P +NGMIX_T_ERR_NOSHEAR +NGMIX_T_PSF_RECONV_1M +NGMIX_T_PSF_RECONV_1P +NGMIX_T_PSF_RECONV_2M +NGMIX_T_PSF_RECONV_2P +NGMIX_T_PSF_RECONV_NOSHEAR + +# flux and error +NGMIX_FLUX_1M +NGMIX_FLUX_1P +NGMIX_FLUX_2M +NGMIX_FLUX_2P +NGMIX_FLUX_NOSHEAR +NGMIX_FLUX_ERR_1M +NGMIX_FLUX_ERR_1P +NGMIX_FLUX_ERR_2M +NGMIX_FLUX_ERR_2P +NGMIX_FLUX_ERR_NOSHEAR + +# magnitudes +MAG_AUTO +MAGERR_AUTO +MAG_WIN +MAGERR_WIN +FLUX_AUTO +FLUXERR_AUTO +FLUX_APER +FLUXERR_APER +FLUX_RADIUS + +# SNR from SExtractor +SNR_WIN + +FWHM_IMAGE +FWHM_WORLD + +# PSF size measured on original image +NGMIX_T_PSF_ORIG_NOSHEAR + +# PSF size measured on reconvolved image +# NGMIX_T_PSF_RECONV_NOSHEAR + +# ngmix moment failure flag +NGMIX_MOM_FAIL diff --git a/workflow/config/cfis/mask_default/MEGAPRIME_star_i_13.8.reg b/workflow/config/cfis/mask_default/MEGAPRIME_star_i_13.8.reg new file mode 100644 index 000000000..4e4164aaf --- /dev/null +++ b/workflow/config/cfis/mask_default/MEGAPRIME_star_i_13.8.reg @@ -0,0 +1,24 @@ +-11.5 68 +-6 186.5 +7 188 +10 64.5 +31 55 +50 38.5 +56.5 11.5 +188 8 +192 -4 +59.5 -11.5 +45 -33 +13.5 -64 +5 -154 +-6 -155 +-11 -64.5 +-40 -44.5 +-51.5 -30.5 +-62.5 -22.5 +-68 -9.5 +-177 -2 +-176 3 +-78 12.5 +-67.5 14.5 +-38.5 50 diff --git a/workflow/config/cfis/mask_default/Messier_catalog.npy b/workflow/config/cfis/mask_default/Messier_catalog.npy new file mode 100644 index 000000000..ef07eb032 Binary files /dev/null and b/workflow/config/cfis/mask_default/Messier_catalog.npy differ diff --git a/workflow/config/cfis/mask_default/Messier_catalog_updated.fits b/workflow/config/cfis/mask_default/Messier_catalog_updated.fits new file mode 100644 index 000000000..6a9f00096 Binary files /dev/null and b/workflow/config/cfis/mask_default/Messier_catalog_updated.fits differ diff --git a/workflow/config/cfis/mask_default/default.ww b/workflow/config/cfis/mask_default/default.ww new file mode 100644 index 000000000..c2797f904 --- /dev/null +++ b/workflow/config/cfis/mask_default/default.ww @@ -0,0 +1,40 @@ +#--------------------------------- Weights ------------------------------------ + +WEIGHT_NAMES weightin.fits # Filename(s) of the input WEIGHT map(s) + +WEIGHT_MIN 0. # Pixel below those thresholds will be flagged +WEIGHT_MAX 1000. # Pixels above those thresholds will be flagged +WEIGHT_OUTFLAGS 1 # FLAG values for thresholded pixels + +#---------------------------------- Flags ------------------------------------- + +FLAG_NAMES flagin.fits # Filename(s) of the input FLAG map(s) + +FLAG_WMASKS 0xff # Bits which will nullify the WEIGHT-map pixels +FLAG_MASKS 0x01 # Bits which will be converted as output FLAGs +FLAG_OUTFLAGS 2 # Translation of the FLAG_MASKS bits + +#---------------------------------- Polygons ---------------------------------- + +POLY_NAMES "" # Filename(s) of input DS9 regions +POLY_OUTFLAGS # FLAG values for polygon masks +POLY_OUTWEIGHTS 0.0 # Weight values for polygon masks +POLY_INTERSECT Y # Use inclusive OR for polygon intersects (Y/N)? + +#---------------------------------- Output ------------------------------------ + +OUTWEIGHT_NAME "w.fits" # Output WEIGHT-map filename +OUTFLAG_NAME flag.fits # Output FLAG-map filename + +#----------------------------- Miscellaneous --------------------------------- + +GETAREA N # Compute area for flags and weights (Y/N)? +GETAREA_WEIGHT 0.0 # Weight threshold for area computation +GETAREA_FLAGS 1 # Bit mask for flag pixels not counted in area +MEMORY_BUFSIZE 256 # Buffer size in lines +VERBOSE_TYPE NORMAL # can be QUIET, NORMAL or FULL +WRITE_XML N # Write XML file (Y/N)? +XML_NAME ww.xml # Filename for XML output +XSL_URL file:///usr/local/share/weightwatcher/ww.xsl + # Filename for XSL style-sheet +NTHREADS 1 # 1 single thread \ No newline at end of file diff --git a/workflow/config/cfis/mask_default/halo_mask.reg b/workflow/config/cfis/mask_default/halo_mask.reg new file mode 100644 index 000000000..c44f25167 --- /dev/null +++ b/workflow/config/cfis/mask_default/halo_mask.reg @@ -0,0 +1,50 @@ + 274.66813 -1.25966 + 272.54579 32.47406 + 266.21222 65.67579 + 255.76731 97.82190 + 241.37579 128.40544 + 223.26462 156.94408 + 201.71942 182.98775 + 177.07997 206.12573 + 149.73486 225.99312 + 120.11532 242.27660 + 88.68848 254.71937 + 55.94996 263.12519 + 22.41606 267.36151 + -11.38436 267.36151 + -44.91826 263.12519 + -77.65678 254.71937 +-109.08362 242.27660 +-138.70315 225.99312 +-166.04827 206.12573 +-190.68772 182.98775 +-212.23292 156.94408 +-230.34409 128.40544 +-244.73561 97.82190 +-255.18052 65.67579 +-261.51409 32.47406 +-263.63643 -1.25966 +-261.51409 -34.99339 +-255.18052 -68.19511 +-244.73561 -100.34123 +-230.34409 -130.92476 +-212.23292 -159.46341 +-190.68772 -185.50708 +-166.04827 -208.64506 +-138.70315 -228.51245 +-109.08362 -244.79593 + -77.65678 -257.23870 + -44.91826 -265.64452 + -11.38436 -269.88084 + 22.41606 -269.88084 + 55.94996 -265.64452 + 88.68848 -257.23870 + 120.11532 -244.79593 + 149.73486 -228.51245 + 177.07997 -208.64506 + 201.71942 -185.50708 + 223.26462 -159.46341 + 241.37579 -130.92476 + 255.76731 -100.34123 + 266.21222 -68.19511 + 272.54579 -34.99339 diff --git a/workflow/config/cfis/mask_default/ngc_cat.fits b/workflow/config/cfis/mask_default/ngc_cat.fits new file mode 100644 index 000000000..f51546da7 Binary files /dev/null and b/workflow/config/cfis/mask_default/ngc_cat.fits differ diff --git a/workflow/config/cfis/star_selection.setools b/workflow/config/cfis/star_selection.setools new file mode 100644 index 000000000..8330a1eff --- /dev/null +++ b/workflow/config/cfis/star_selection.setools @@ -0,0 +1,103 @@ +## SETools configuration file for star/galaxy separation based on size/mag properties + +[MASK:preselect] +MAG_AUTO > 0 +MAG_AUTO < 21 +FWHM_IMAGE > 0.3 / 0.187 +FWHM_IMAGE < 1.5 / 0.187 +FLAGS == 0 +IMAFLAGS_ISO == 0 +NO_SAVE + +[MASK:flag] +FLAGS == 0 +IMAFLAGS_ISO == 0 +NO_SAVE + + +[MASK:star_selection] +# Star selection using the FWHM mode +MAG_AUTO > 18. +MAG_AUTO < 22. +FWHM_IMAGE <= mode(FWHM_IMAGE{preselect}) + 0.2 +FWHM_IMAGE >= mode(FWHM_IMAGE{preselect}) - 0.2 +FLAGS == 0 +IMAFLAGS_ISO == 0 + +[MASK:fwhm_mag_cut] +FWHM_IMAGE > 0 +FWHM_IMAGE < 40 +MAG_AUTO < 35 +FLAGS == 0 +IMAFLAGS_ISO == 0 +NO_SAVE + +# Split the 'star_selection' sample into +# two random sub-samples with ratio 80/20 +[RAND_SPLIT:star_split] +RATIO = 20 +MASK = star_selection + +# The following selection is only used for plotting + +[PLOT:size_mag] +TYPE = plot +FORMAT = png +X_1 = FWHM_IMAGE{fwhm_mag_cut} +Y_1 = MAG_AUTO{fwhm_mag_cut} +X_2 = FWHM_IMAGE{star_selection} +Y_2 = MAG_AUTO{star_selection} +MARKER_1 = + +MARKER_2 = . +MARKERSIZE_1 = 3 +MARKERSIZE_2 = 3 +LABEL_1 = All +LABEL_2 = "Stars, mean FWHM: @mean(FWHM_IMAGE{star_selection})*0.187@ arcsec" +TITLE = "Stellar locus" +XLABEL = "FWHM (pix)" +YLABEL = Mag + +[PLOT:hist_mag_stars] +TYPE = hist +FORMAT = png +Y = MAG_AUTO{star_selection} +BIN = 20 +LABEL = "stars" +XLABEL = "Magnitude" +YLABEL = "Number" +TITLE = "Magnitude of stars" + +[PLOT:fwhm_field] +TYPE = scatter +FORMAT = png +X = X_IMAGE{star_selection} +Y = Y_IMAGE{star_selection} +SCATTER = FWHM_IMAGE{star_selection}*0.186 +MARKER = . +LABEL = "FWHM (arcsec)" +TITLE = "FWHM of stars" +XLABEL = "X (pix)" +YLABEL = "Y (pix)" + +[PLOT:mag_star_field] +TYPE = scatter +FORMAT = png +X = X_IMAGE{star_selection} +Y = Y_IMAGE{star_selection} +SCATTER = MAG_AUTO{star_selection} +MARKER = . +LABEL = "Magnitude" +TITLE = "Magnitude of stars" +XLABEL = "X (pix)" +YLABEL = "Y (pix)" + +[STAT:star_stat] +"Nb objects full cat" = len(FWHM_IMAGE) +"Nb objects not masked" = len(FWHM_IMAGE{flag}) +"Nb stars" = len(FWHM_IMAGE{star_selection}) +"stars/deg^2" = len(FWHM_IMAGE{star_selection})/4612./0.187*3600.*1./2048./0.187*3600. +"Mean star fwhm selected (arcsec)" = mean(FWHM_IMAGE{star_selection})*0.187 +"Standard deviation fwhm star selected (arcsec)" = std(FWHM_IMAGE{star_selection})*0.187 +"Mode fwhm used (arcsec)" = mode(FWHM_IMAGE{preselect})*0.187 +"Min fwhm cut (arcesec)" = mode(FWHM_IMAGE{preselect})*0.187-0.1*0.187 +"Max fwhm cut (arcsec)" = mode(FWHM_IMAGE{preselect})*0.187+0.1*0.187 diff --git a/workflow/rules/exposure.smk b/workflow/rules/exposure.smk new file mode 100644 index 000000000..e41920137 --- /dev/null +++ b/workflow/rules/exposure.smk @@ -0,0 +1,366 @@ +"""Exposure chain — per exposure, keyed by exp base id (dedup is structural). + + exp_get_images -> exp_split -----> exp_mask -> exp_psf + -> exp_star_cat --/ + star_catalogue ---------------/ + +``star_catalogue`` is campaign-level, not per-exposure: one fetch of the whole +footprint's stars, which every exposure's ``exp_star_cat`` then cuts locally. + +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 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 +by the accumulating index. A temp() here would delete an exposure the moment +this invocation's readers finished and cascade destructive reruns across spatial +neighbours the next time a tile is appended. + +GROUPING (``group: "exp_short"``) covers exp_split and exp_mask, and only them — +one sbatch per exposure for two jobs whose medians are 1:28 and 1:54, well under +the 15-minute floor Alliance policy asks us to bundle away. The composition +rules are in prepare.smk's docstring; this chain is linear too, so the group +asks max(mem_mb) = 8000*attempt, max(threads) = 8, sum(runtime) = 240 min. + +The two rules NOT in it are structural, not taste: + * exp_psf is heavy (16 GB, 4 h) and never fuses with a short rule; + * exp_get_images cannot join, because ``exp_star_cat`` — a LOCALRULE, and so + ungroupable — sits between it and exp_mask. Pulling get_images in would make + the group both a dependency and a dependent of exp_star_cat, i.e. a cycle. + Starting the group at exp_split leaves star_cat's inputs entirely upstream + of it, so the group has one clean external edge. +Different exposures share no DAG edge, so this is one group job per exposure. + +NUMBER_LIST is set only for exp_split (its numbering scheme IS the exposure id); +never for get_images / exp_mask / exp_psf, whose per-CCD or download numbering +would make the #746 startup validation turn tolerated per-CCD attrition into a +whole-exposure hard failure. That is now a property of the committed configs +(config_exp_Sp.ini has NUMBER_LIST = $SP_UNIT_NUM; Gie/Ma/psfex have none). +""" + +rule exp_get_images: + output: + manifest = f"{EXP_DIR}/manifests/exp_get_images.json" + log: + f"{EXP_DIR}/logs/exp_get_images.json" + params: + pre = lambda wc: unit_pre("exp_get_images", "exp", wc.exp, + exp_name=EXP[wc.exp]), + script_hash = SCRIPT_HASH + threads: 1 + retries: 2 + resources: + mem_mb = lambda wc, attempt: 4000 * attempt, + runtime = 60 + shell: + sp_shell("exp_get_images", "config_exp_Gie.ini") + +# --- mask star catalogues --------------------------------------------------- +# Two rules, and the split between them is the design: the NETWORK is a function +# of the campaign's sky area, the per-exposure catalogue is a local cut. +# +# `star_catalogue` fetches the footprint's GSC 2.3 stars once, one Vizier query +# per HEALPix chunk, into a run-independent chunk store under config `star_cats`. +# `exp_star_cat` then reads the chunks covering an exposure's focal plane and +# cuts them to it — no network at all. workflow/scripts/star_cats.py holds both +# halves and the geometry they must agree on; its docstring is the reference for +# the chunking and the padding. +# +# The arithmetic: exposures overlap ~7-10 deep and a tile's exposures all look at +# the same square degree, so the old one-cone-per-exposure design re-fetched the +# same sky ~8 times over. Chunked by sky, a full-UNIONS footprint is ~1.5k +# queries where the exposure count would have been ~25k, and a campaign that +# grows within the fetched footprint issues none. + +# The container's certifi bundle. The host leaks SSL_CERT_FILE / CURL_CA_BUNDLE +# pointing at a path that does not exist inside the image, so requests is pointed +# at the bundle explicitly (proven in the p3-batch1 bash precedent). +STAR_CAT_CA = "/app/.venv/lib/python3.12/site-packages/certifi/cacert.pem" + +# The rules run star_cats.py inside the container (healpy, astroquery, astropy) +# but call apptainer THEMSELVES rather than letting the SDM wrap them +# (`container: None` on both): the CA bundle above and the exposure rule's +# host-side farm loop both need the explicit exec. bin/sp has loaded the +# apptainer module. PYTHONPATH pins this checkout's src/ for +# shapepipe.utilities.{vizier,cfis} — the mirror-retry query and the tile-ID +# grid convention are library code, not copies. +def in_container(cmd, *, network=False): + ca = ("," + ",".join(f"{k}={STAR_CAT_CA}" for k in + ("REQUESTS_CA_BUNDLE", "SSL_CERT_FILE", "CURL_CA_BUNDLE")) + if network else "") + return (f"apptainer exec --cleanenv --home {Path.home()}" + f" --bind /project --bind /scratch" + f" --env PYTHONPATH={REPO_DIR}/src{ca}" + f" '{config['container']}' {cmd}") + + +# The campaign's star catalogue: a first-class durable science product, keyed by +# sky rather than by run. Chunk-need is recomputed from the tile list on every +# run and only the missing chunks are fetched, so appending tiles costs exactly +# the chunks they add. +# +# A LOCALRULE (declared in the Snakefile): it is one job of network I/O, and the +# fetch loop is a 4-wide thread pool inside it — the same modest concurrency the +# per-exposure rule reached by accident through --local-cores, now an explicit +# number that does not scale with the head node's CPU count. +# +# `tile_list_hash` is what makes the incremental behaviour visible to the DAG. +# The tile list is parse-time config, not a rule input (and the profile drops the +# `input` rerun-trigger anyway), so appending tiles would otherwise leave this +# rule up to date against a footprint that has grown. Hashing the list into a +# param reruns it, and the rerun fetches only what is new. +STAR_CAT_MANIFEST = f"{RUN_DIR}/manifests/star_catalogue.json" + + +rule star_catalogue: + output: + manifest = STAR_CAT_MANIFEST + # No `log:` — see write_manifest() in star_cats.py: this rule runs no + # shapepipe_run and computes no completeness verdict, so there is nothing to + # split between a manifest and a log. Under `set -euo pipefail` it either + # completes or aborts, and snakemake's captured stderr is the evidence. + # `cmd` is a params value, so it is substituted AFTER the shell string is + # formatted: a `{output.manifest}` placeholder in here would survive + # literally (see unit_pre in the Snakefile). Hence the explicit path. + params: + cmd = in_container( + f"python {SCRIPTS}/star_cats.py fetch" + f" --tile-list '{config['tile_list']}' --store '{STAR_CATS}'" + f" --manifest '{STAR_CAT_MANIFEST}'", network=True), + tile_list_hash = hashlib.md5( + Path(config["tile_list"]).read_bytes()).hexdigest()[:12], + script_hash = STAR_CAT_HASH + container: + None + threads: 4 + retries: 2 + resources: + mem_mb = 4000, + runtime = 720 + shell: + "set -euo pipefail\n{params.cmd}" + + +# The per-exposure catalogue and the 40 per-CCD symlinks the mask module's +# numbering scheme needs. Local: one header read for the focal-plane footprint, +# a load of the chunks covering it, a radial cut. +# +# A LOCALRULE for the same reason as clean_exposure: seconds of work, and one +# sbatch per exposure would be ~20k submissions well under the 15-minute floor +# cluster policy asks us to bundle away. +# +# The per-unit farm is a REAL directory holding exactly this exposure's 40 +# numbers, and that is load-bearing: config_exp_Ma.ini reads it as an INPUT_DIR +# and the file handler INTERSECTS the numbers found across INPUT_DIRs, so a +# symlink to a shared whole-store pool contributes every other exposure's numbers +# and the intersection is empty ("numbers ... do not intersect", live). +# +# TWO declared outputs, and the second one is the point. +# +# The manifest keeps the "one rule, one manifest" currency of every other rule: +# written last, unique to this rule, a record of what the farm points at, and +# deleted by clean_exposure so a reclaimed exposure rebuilds its farm from the +# chunk store at no network cost. +# +# But a manifest attests FOREVER, and the two things it attests to both live +# outside the unit's manifests/ dir: the cut catalogue on /scratch (60-day purge) +# and the farm itself. Either can vanish under a manifest that still says +# "complete", and then exp_mask runs against nothing. So the ccd-0 farm link is +# declared too — one link stands for all 40, they are created by the same loop +# in the same instant, and declaring 40 buys nothing. Snakemake's existence test +# is os.path.exists, which FOLLOWS symlinks and is therefore False for a link +# whose target the purge removed. A purged cut or a deleted farm makes the rule +# out of date, it reruns, and it re-cuts or re-links as needed. +def star_cat_cmd(exp): + """The whole rule body, as bash — carried as a params value, never inlined + in ``shell:``: it contains literal ``{}`` (the manifest JSON) and snakemake + formats a shell string once, which would consume those braces.""" + cut_dir = f"{STAR_CATS}/exp" + cat = f"{cut_dir}/star_cat-{exp}.fits" + work = exp_dir(exp) + farm = f"{work}/star_cat_exp" + images = f"{work}/output/run_sp_exp_Gie/get_images_runner/output" + manifest = exp_manifest(exp, "exp_star_cat") + body = json.dumps({ + "stage": "exp_star_cat", "level": "exp", "unit": exp, + "status": "complete", "cat": cat, "link_dir": farm, "n_links": 40, + }, indent=2, sort_keys=True) + return "\n".join([ + "set -euo pipefail", + # LEGACY-SYMLINK HAZARD. Unit dirs built before this rule existed carry + # star_cat_exp as a SYMLINK into the old shared star-cat pool. `mkdir -p` + # is a no-op on an existing symlink-to-directory, so the 40-link loop + # below followed it and wrote this exposure's links INTO THE SHARED POOL + # (520 stray links found live). Replace the link — never `rm -rf` it, + # which would recurse into the pool, and never touch a real directory: + # a real farm is this rule's own output and `ln -sfn` refreshes it. + f"[ -L '{farm}' ] && rm -f '{farm}' || true", + f"mkdir -p '{cut_dir}' '{farm}' '{work}/manifests'", + in_container(f"python {SCRIPTS}/star_cats.py cut" + f" --images '{images}' --store '{STAR_CATS}'" + f" --out '{cat}'"), + f"test -s '{cat}'", + # The fan-out the file handler's NUMBERING_SCHEME wants: 40 links to the + # one focal-plane catalogue (pattern from the p3-batch1 precedent). + f"for ccd in $(seq 0 39); do ln -sfn '{cat}' " + f"'{farm}/star_cat-{exp}-'\"$ccd\"'.fits'; done", + # Byte-stable, and written only after the links exist: an unconditional + # write would move the mtime, which is a rerun-trigger. + f"tmp='{manifest}.tmp'", + "cat > \"$tmp\" <<'SP_STAR_CAT_JSON'", + body, + "SP_STAR_CAT_JSON", + f"cmp -s \"$tmp\" '{manifest}' && rm -f \"$tmp\" || mv -f \"$tmp\" '{manifest}'", + ]) + + +rule exp_star_cat: + input: + rules.exp_get_images.output.manifest, + # The chunks this cut reads. star_cats.py fails loudly on a chunk that is + # missing anyway, but the edge is what makes the fetch happen first. + rules.star_catalogue.output.manifest + output: + manifest = f"{EXP_DIR}/manifests/exp_star_cat.json", + # The sentinel: ccd-0 of the 40-link farm (see above). + link = f"{EXP_DIR}/star_cat_exp/star_cat-{{exp}}-0.fits" + # No `log:`, for the same reason as star_catalogue above. + params: + cmd = lambda wc: star_cat_cmd(wc.exp), + # star_cats.py is external to the shell string, so the `code` + # rerun-trigger does not see it — same reason SCRIPT_HASH exists. + script_hash = STAR_CAT_HASH + container: + None + threads: 1 + retries: 2 + resources: + mem_mb = 4000, + runtime = 10 + shell: + "{params.cmd}" + +# Split the multi-HDU exposure into single-CCD files (+ headers-*.npy, which the +# tiles' merge_headers reads). +rule exp_split: + group: "exp_short" + input: + rules.exp_get_images.output.manifest + output: + manifest = f"{EXP_DIR}/manifests/exp_split.json" + log: + f"{EXP_DIR}/logs/exp_split.json" + params: + pre = lambda wc: unit_pre("exp_split", "exp", wc.exp), + script_hash = SCRIPT_HASH + threads: 8 + resources: + mem_mb = lambda wc, attempt: 8000 * attempt, + runtime = 120 + shell: + sp_shell("exp_split", "config_exp_Sp.ini") + +rule exp_mask: + group: "exp_short" + input: + # Both inputs are real INPUT_DIRs of config_exp_Ma.ini: the split CCDs + # and this exposure's own star_cat_exp farm. + rules.exp_split.output.manifest, + rules.exp_star_cat.output.manifest + output: + manifest = f"{EXP_DIR}/manifests/exp_mask.json" + log: + f"{EXP_DIR}/logs/exp_mask.json" + params: + pre = lambda wc: unit_pre("exp_mask", "exp", wc.exp), + script_hash = SCRIPT_HASH + threads: 4 + resources: + mem_mb = lambda wc, attempt: 8000 * attempt, + runtime = 120 + shell: + sp_shell("exp_mask", "config_exp_Ma.ini") + +# SExtractor -> setools star selection -> PSFEx model -> psfex_interp, per CCD. +# setools may reject a sparse CCD (~0.2% attrition) — tolerated by the floor's +# :warn on psfex_interp_runner. +rule exp_psf: + input: + rules.exp_mask.output.manifest + output: + manifest = f"{EXP_DIR}/manifests/exp_psf.json" + log: + f"{EXP_DIR}/logs/exp_psf.json" + params: + pre = lambda wc: unit_pre("exp_psf", "exp", wc.exp), + script_hash = SCRIPT_HASH + threads: 8 + retries: 2 + benchmark: + # BESIDE manifests/, not inside it: clean_exposure deletes manifests/ + # wholesale, and this tsv is the measured-memory feed for mem_mb sizing + # (D4). Inside manifests/ it died with the first reclamation and took + # the campaign's only record of exp_psf's real footprint with it. + f"{EXP_DIR}/exp_psf.benchmark.tsv" + resources: + mem_mb = lambda wc, attempt: 16000 * attempt, + runtime = 240 + shell: + sp_shell("exp_psf", "config_exp_psfex.ini") + + +# --- 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 +# an exposure's consumer set closes over the CAMPAIGN. The index supplies that +# set (EXP_TILES, accumulated across invocations); the input is every consuming +# tile's tile_vignets manifest — vignets is the last stage that reads exposure +# products, everything after it reads tile-level files. +# +# Three properties make the late append behave (see clean_exposure.py): +# * the job deletes the exposure's manifests too, so a tile appended after the +# clean sees an unbuilt chain and regenerates it instead of running against +# an empty store. The tombstone deliberately does NOT stand in for those +# manifests — it is not an input to anything but itself. +# * a finished tile is not disturbed: Snakemake demands a missing intermediate +# only when something downstream of it must run. +# * params.consumers carries the consumer set, so growing it makes the +# tombstone stale under the default `params` rerun-trigger; the clean job +# reruns after the new tile's vignets, against the enlarged set. +# +# The tile side reads the exposure manifests through ancient() (see tile.smk), +# which is what keeps this deletion from rebuilding every neighbouring tile. +# This rule's OWN inputs are deliberately not ancient: a tile that really did +# rebuild its vignets must reschedule the cleans of the exposures it read. +# +# A localrule (declared in the Snakefile): it is an rmtree, not science, and one +# sbatch per exposure would be ~20k scheduler submissions at DR6 scale. Local +# execution serialises them under local-cores, which costs nothing at rmtree +# speed and never blocks the compute chains (this rule is in none of them). +rule clean_exposure: + input: + # ONLY the consumers this invocation may actually build. A consumer that + # is out of scope had its vignets manifest checked for existence at parse + # time (clean_targets' eligibility test) — declaring it here as well would + # pull that finished tile's whole chain into the DAG, where a rebuilt + # shared exposure then reruns it. That is how one damaged tile reached its + # 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] + output: + tombstone = f"{EXP_DIR}/cleaned.json" + params: + consumers = lambda wc: ",".join(clean_consumers(wc.exp)), + script_hash = CLEAN_HASH + threads: 1 + resources: + mem_mb = 2000, + runtime = 30 + shell: + f"python {SCRIPTS}/clean_exposure.py" + " --exp-dir $(dirname {output.tombstone}) --exp {wildcards.exp}" + " --tombstone {output.tombstone} --consumers '{params.consumers}'" diff --git a/workflow/rules/prepare.smk b/workflow/rules/prepare.smk new file mode 100644 index 000000000..96a0ae261 --- /dev/null +++ b/workflow/rules/prepare.smk @@ -0,0 +1,94 @@ +"""Invocation 1 — PREPARE: one static chain per tile. + + tile_get_images -> tile_uncompress -> tile_find_exposures + +Known from the tile list alone, cheap, wide, idempotent. find_exposures parses +the tile FITS HISTORY header into ``exp_numbers--.txt`` — the +data-derived tile->exposure edge that invocation 2's parse aggregates into the +index. Nibi compute nodes have internet, so downloads run in-DAG (no login-node +tier). + +All three rules carry ``group: "tile_prep"``, so one tile's whole chain is ONE +sbatch instead of three (medians 0:41 / ~0:40 / 0:15 — all far under the +15-minute floor Alliance policy asks us to bundle away, and at DR6 scale three +submissions per tile is a scheduler load out of all proportion to the work). +Group membership is per connected DAG component and distinct tiles share no +edge, so this is exactly one group job per tile, never a cross-tile bundle. +Rules stay group-compatible: shell only, no mid-chain localrules, no pipe outputs. + +Group resource composition (snakemake 9.23, ``GroupResources.basic_layered`` in +snakemake/resources.py): jobs are laid out per toposort level; within a level +non-additive resources (mem_mb, cpus) SUM — split into layers when a global +constraint is exceeded, the group's width being the widest layer — while the +additive resource ``runtime`` is maxed within a layer and SUMMED across layers. +This chain is strictly linear, one job per level, so the group asks for +max(mem_mb) = 8000*attempt, max(threads) = 4 and sum(runtime) = 150 min. +Attempt scaling survives grouping: ``GroupJob.attempt``'s setter clears the +cached group resources and re-sets ``attempt`` on every member (jobs.py), and +``GroupJob.restart_times`` is the max over members — so tile_get_images' +``retries: 2`` still governs. A retry re-runs the whole group, which is safe +because every rule ``rm -rf``s its own run dir at start. + +Star catalogues for masking are NOT a prepare-phase concern and not pre-run +input: the compute DAG fetches the campaign footprint's stars once +(``star_catalogue``) and cuts them per exposure (``exp_star_cat``), both in +exposure.smk, into a run-independent store. The tile side has no star-cat node +because it has no mask rule yet — see tile.smk. +""" + +# NUMBER_LIST is never set for get_images (download stage; nothing on disk to +# validate against — #746 would hard-fail the unit). The committed +# config_tile_Git.ini simply has no NUMBER_LIST, so this is now a property of the +# config, not of an injection step. +rule tile_get_images: + group: "tile_prep" + output: + manifest = f"{TILE_DIR}/manifests/tile_get_images.json" + log: + f"{TILE_DIR}/logs/tile_get_images.json" + params: + pre = lambda wc: unit_pre("tile_get_images", "tile", wc.tile), + script_hash = SCRIPT_HASH + threads: 1 + retries: 2 + resources: + mem_mb = lambda wc, attempt: 4000 * attempt, + runtime = 60 + shell: + sp_shell("tile_get_images", "config_tile_Git.ini") + +rule tile_uncompress: + group: "tile_prep" + input: + rules.tile_get_images.output.manifest + output: + manifest = f"{TILE_DIR}/manifests/tile_uncompress.json" + log: + f"{TILE_DIR}/logs/tile_uncompress.json" + params: + pre = lambda wc: unit_pre("tile_uncompress", "tile", wc.tile), + script_hash = SCRIPT_HASH + threads: 4 + resources: + mem_mb = lambda wc, attempt: 8000 * attempt, + runtime = 60 + shell: + sp_shell("tile_uncompress", "config_tile_Uz.ini") + +rule tile_find_exposures: + group: "tile_prep" + input: + rules.tile_uncompress.output.manifest + output: + manifest = f"{TILE_DIR}/manifests/tile_find_exposures.json" + log: + f"{TILE_DIR}/logs/tile_find_exposures.json" + params: + pre = lambda wc: unit_pre("tile_find_exposures", "tile", wc.tile), + script_hash = SCRIPT_HASH + threads: 1 + resources: + mem_mb = lambda wc, attempt: 2000 * attempt, + runtime = 30 + shell: + sp_shell("tile_find_exposures", "config_tile_Fe.ini") diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk new file mode 100644 index 000000000..2739a1932 --- /dev/null +++ b/workflow/rules/tile.smk @@ -0,0 +1,338 @@ +"""Tile post-chain — per tile: gather exposures, then detect / PSF / shape / catalogue. + + tile_exp_forest + tile_merge_headers -> tile_detect -> tile_vignets -> tile_ngmix x N + -> tile_merge_cats -> tile_make_cat + +The DAG edge to the exposures is always the exposures' MANIFESTS, looked up +through the index (TILE_EXP). The per-tile exposure "forest" — a symlink view of +exactly this tile's exposures' products — exists only so the ShapePipe configs +have one deterministic ``$SP_EXP`` to glob; it is NEVER the edge. Its 2-char +shard level is not cosmetic: ``exp_utils.get_exp_output_files`` hardwires +``///output/run_sp_*`` into its glob, so a flat forest +makes every tile gather stage fail "No split_exp_runner output found". + +All rules are group-compatible (shell only, no mid-chain localrules), and the two +short regions are grouped, per the composition rules in prepare.smk's docstring. +Distinct tiles share no edge, so each is one group job per tile: + +* ``group: "tile_gather"`` — tile_exp_forest (2 GB, 20 min) and + tile_merge_headers (median 0:38; 8 GB, 4 threads, 120 min). The group asks max + mem_mb = 8000*attempt, max threads = 4, sum runtime = 140. tile_detect also + consumes the forest, but a consumer OUTSIDE the group is just an ordinary DAG + edge on the group job — it does not pull tile_detect (16 GB) in. +* ``group: "tile_finish"`` — tile_merge_cats (median 0:15) and tile_make_cat + (1:34), two short jobs of identical shape (16 GB, 8 threads, 120 min each; + max mem_mb = 16000*attempt, max threads = 8, sum runtime = 240). + +Nothing else joins either. The heavy middle (tile_detect, tile_vignets, +tile_ngmix) is where the wall time is, and fusing a short rule onto one of those +would reserve its footprint for the short job too. + +Note there is no `tile_mask` rule: the committed config chain is the +"sx_nomask" tile_detect variant (config_tile_Sx.ini reads Git + Uz + Mh, no mask +run), and no tile-mask config was committed in the S2 sweep. Adding the masked +variant is a config + one rule, at the config selector the PRD describes. + +That rule also needs a tile-side analogue of ``exp_star_cat``: tile star cats key +on TILE id, so they are a separate cache namespace and a separate node, and the +earliest point it can run is after ``tile_uncompress`` (create_star_cat.py's +``-k tile`` mode reads the uncompressed tile image's primary header). The mask +config would then read a real per-unit ``$SP_RUN/star_cat_tiles`` directory, +built the same way and for the same reason (the file handler intersects numbers +across INPUT_DIRs, so a shared pool cannot be symlinked in wholesale). +""" + +def tile_exp(wc): + return TILE_EXP.get(wc.tile, []) + +# --- the tile->exposure edge, and why it is cut for finished tiles ---------- +# +# THE cascade fix. clean_exposure deletes the exposure's manifests on purpose: +# that is what makes a tile appended later rebuild the chain instead of running +# against an empty store. But those manifests are the tile side's inputs, and an +# exposure is read by ~7-10 tiles. So the moment ONE tile's chain rebuilt an +# exposure, every other tile reading it saw "input files updated by another job" +# and reran — and that rerun rebuilt ITS exposures, which reran THEIR other +# consumers, propagating across the whole exposure-overlap connected component. +# On fixture t4, asking for one damaged tile scheduled all four tiles' chains. +# +# Two mechanisms, and only the second one actually cuts it: +# +# 1. ancient() on every exposure-manifest edge. Correct on its own terms — a +# tile has no business rerunning because an exposure manifest is NEWER — and +# it is what keeps a pure-mtime disturbance (a re-touched manifest, a +# restored backup) from waking finished tiles. But ancient() governs +# TIMESTAMPS only. Snakemake propagates "my input is produced by a job that +# will run in this DAG" separately, and ancient does not suppress it +# (measured: t4 counts were identical with ancient alone). +# +# 2. Cutting the RECLAIMED edges of a FINISHED tile — the mechanism that works. +# A tile whose final_cat is on disk needs nothing further from its +# exposures: it has already extracted everything it will ever read. So for +# such a tile the input list drops the manifests that are GONE, and the +# propagation has nowhere to go. An UNfinished tile keeps its full edge set +# and therefore still drags in — and rebuilds — every exposure it needs, +# which is the accepted price of a late append, unchanged. +# +# Only the missing ones are dropped, never a manifest that still exists: a +# campaign that has cleaned nothing then declares exactly the edges it always +# did, and the cut cannot perturb it. +# +# The marker is final_cat, not the tile's own vignets manifest: on a tile whose +# catalogue was lost, the vignets manifest still exists while the vignette store +# (temp()) does not, so keying on vignets would cut the edge on exactly the tile +# that has to rerun, and run it against a deleted exposure store. +# +# THE CUT REQUIRES THE `input` RERUN-TRIGGER TO BE OFF (profiles/nibi sets the +# trigger list). Dropping an input is itself a change in the set of input files, +# which that trigger reads as a reason to rerun — reinstating the very cascade, +# now as "Set of input files has changed", and running finished tiles against a +# store that is gone. Measured on fixture t4, one damaged tile of four: 82 jobs +# with neither fix, 70 with the cut but the trigger on, 28 with both (= exactly +# the damaged tile's own chain, its two exposures, and the clean jobs). +# +# The cost, stated plainly: `--forcerun` on a tile whose final_cat exists will +# NOT rebuild its reclaimed exposures, because those edges are not in the DAG. +# Delete that tile's final_cat first and the whole chain comes back. +# +# What none of this weakens: clean_exposure's own inputs are neither ancient nor +# cut, so a tile that really did rebuild its vignets still reschedules the cleans +# of the exposures it read, and a grown consumer set still travels through +# params.consumers. +def tile_finished(tile): + return Path(final_cat(tile)).exists() + + +def exp_manifests(wc, stage): + paths = [exp_manifest(e, stage) for e in tile_exp(wc)] + if tile_finished(wc.tile): + paths = [p for p in paths if Path(p).exists()] + return [ancient(p) for p in paths] + +def tile_exp_split(wc): return exp_manifests(wc, "exp_split") +def tile_exp_mask(wc): return exp_manifests(wc, "exp_mask") +def tile_exp_psf(wc): return exp_manifests(wc, "exp_psf") +def tile_exp_all(wc): return tile_exp_split(wc) + tile_exp_mask(wc) + tile_exp_psf(wc) + + +# Build the per-tile symlink forest. Declaring the exposure manifests as input +# makes this wait on its exposures; the forest itself is only the $SP_EXP view. +# Its output stays a directory() (it has no ShapePipe run dir and no manifest — +# it is not a shapepipe_run at all). +rule tile_exp_forest: + group: "tile_gather" + input: + tile_exp_all + output: + forest = directory(f"{TILE_DIR}/exp_forest") + params: + cmd = lambda wc: (f"python {SCRIPTS}/build_forest.py --tile {wc.tile} " + f"--run-dir {RUN_DIR} --index {INDEX_DB}"), + # build_forest.py's own content hash rides here and nowhere else. + script_hash = FOREST_HASH + threads: 1 + resources: + mem_mb = 2000, + runtime = 20 + shell: + # --forest {output} lives in the shell string: snakemake formats shell + # ONCE, so an {output} placeholder inside params.cmd would survive + # literally and every forest job would race one './{output}'. + "{params.cmd} --forest {output.forest}" + +# Merge single-exposure WCS headers into the tile-level sqlite +# (log_exp_headers--.sqlite, which Sx / PiViVi / ngmix consume). +# Reads headers-*.npy through the forest -> the split manifests are the edge. +rule tile_merge_headers: + group: "tile_gather" + input: + forest = rules.tile_exp_forest.output.forest, + split = tile_exp_split, + # config_tile_Mh_exp.ini reads run_sp_tile_Fe output, which the PREPARE + # phase produced. Declaring the Fe manifest gives the COMPUTE DAG a + # regeneration path for it instead of a silent dependency on a + # previous invocation (prepare.smk is included in every parse, so the + # rule exists here too). Normally a satisfied no-op. + fe = f"{TILE_DIR}/manifests/tile_find_exposures.json", + output: + manifest = f"{TILE_DIR}/manifests/tile_merge_headers.json" + log: + f"{TILE_DIR}/logs/tile_merge_headers.json" + params: + pre = lambda wc: unit_pre("tile_merge_headers", "tile", wc.tile, + forest=forest_dir(wc.tile)), + script_hash = SCRIPT_HASH + threads: 4 + resources: + mem_mb = lambda wc, attempt: 8000 * attempt, + runtime = 120 + shell: + sp_shell("tile_merge_headers", "config_tile_Mh_exp.ini") + +# SExtractor object detection on the tile. +rule tile_detect: + input: + uz = f"{TILE_DIR}/manifests/tile_uncompress.json", + mh = rules.tile_merge_headers.output.manifest, + output: + manifest = f"{TILE_DIR}/manifests/tile_detect.json" + log: + f"{TILE_DIR}/logs/tile_detect.json" + params: + pre = lambda wc: unit_pre("tile_detect", "tile", wc.tile), + script_hash = SCRIPT_HASH + threads: 8 + resources: + mem_mb = lambda wc, attempt: 16000 * attempt, + runtime = 180 + shell: + sp_shell("tile_detect", "config_tile_Sx.ini") + +# PSFEx interpolation to galaxies + vignet postage stamps: the last stage that +# reads exposure products, and the bulk intra-tile intermediate. +# +# The vignette store is declared as a SECOND, temp(directory()) output alongside +# the manifest. This is the ONE scoped exception to "no directory() outputs" +# (D5): the store is ~tens of GB per tile and must be reclaimed when its last +# intra-tile reader finishes, but it must not become DAG currency — so the +# manifest stays the edge, and the directory rides along purely so native temp() +# fires at the right moment. Its readers (tile_ngmix, tile_make_cat) declare +# BOTH. `--notemp` keeps it for debugging. +rule tile_vignets: + input: + sx = rules.tile_detect.output.manifest, + forest = rules.tile_exp_forest.output.forest, + split = tile_exp_split, + psf = tile_exp_psf, + # config_tile_PiViVi.ini reads run_sp_tile_Fe output — same reason as + # tile_merge_headers above. + fe = f"{TILE_DIR}/manifests/tile_find_exposures.json", + output: + manifest = f"{TILE_DIR}/manifests/tile_vignets.json", + store = temp(directory(f"{TILE_DIR}/output/run_sp_tile_PiViVi")), + log: + f"{TILE_DIR}/logs/tile_vignets.json" + params: + pre = lambda wc: unit_pre("tile_vignets", "tile", wc.tile, + forest=forest_dir(wc.tile)), + script_hash = SCRIPT_HASH + threads: 16 + resources: + mem_mb = lambda wc, attempt: 32000 * attempt, + runtime = 240 + shell: + sp_shell("tile_vignets", "config_tile_PiViVi.ini") + +# ngmix shape measurement — N chunks per tile (D4). Each chunk computes its own +# CLOSED object-ID range at EXECUTION time from this tile's own sexcat: a params +# function cannot, because params evaluate before the sexcat exists. Closed, not +# open-ended: `ID_OBJ_MAX = -1` on the last chunk was the 13-hour straggler's +# root cause (ngmix treats id_obj_max <= 0 as unbounded). +# +# Chunks write nothing shared: each has its own run_sp_tile_ngmix_Ngu, and +# merge_sep_cats — DAG-serialised after all chunks — is the gather. +rule tile_ngmix: + input: + vignets = rules.tile_vignets.output.manifest, + store = rules.tile_vignets.output.store, + sx = rules.tile_detect.output.manifest, + output: + manifest = f"{TILE_DIR}/manifests/tile_ngmix_{{chunk}}.json", + # temp(directory()) for the same reason as the vignette store above. + chunkdir = temp(directory(f"{TILE_DIR}/output/run_sp_tile_ngmix_Ng{{chunk}}u")), + log: + f"{TILE_DIR}/logs/tile_ngmix_{{chunk}}.json" + params: + pre = lambda wc: unit_pre( + "tile_ngmix", "tile", wc.tile, + env={"SP_NGMIX_CHUNK": wc.chunk, "NGMIX_N_CHUNKS": NGMIX_CHUNKS}, + # Two steps, not `eval "$(...)"`: a command substitution inside eval + # discards the script's exit status, so a missing sexcat would fall + # through to shapepipe_run with an unset range and fail as something + # else. Capture, check, then eval — the range script fails as itself. + pre_run=[f'ngmix_range_out=$(python {SCRIPTS}/ngmix_range.py --run-dir ' + f'"$SP_RUN" --chunk {wc.chunk} --n-chunks {NGMIX_CHUNKS}) ' + f'|| exit 1', + 'eval "$ngmix_range_out"']), + script_hash = SCRIPT_HASH + threads: 4 + retries: 2 + benchmark: + f"{TILE_DIR}/manifests/tile_ngmix_{{chunk}}.benchmark.tsv" + resources: + mem_mb = lambda wc, attempt: 14000 * attempt, + runtime = 720 + shell: + sp_shell("tile_ngmix", "config_tile_Ng_template.ini") + + +def ngmix_manifests(wc): + return [f"{tile_dir(wc.tile)}/manifests/tile_ngmix_{k}.json" + for k in range(1, NGMIX_CHUNKS + 1)] + +def ngmix_chunkdirs(wc): + return [f"{tile_dir(wc.tile)}/output/run_sp_tile_ngmix_Ng{k}u" + for k in range(1, NGMIX_CHUNKS + 1)] + +# The gather: merge the N chunk catalogues. N_SPLIT_MAX comes from the workflow's +# own chunk count via $NGMIX_N_CHUNKS (env-expanded by the module). +rule tile_merge_cats: + group: "tile_finish" + input: + manifests = ngmix_manifests, + chunkdirs = ngmix_chunkdirs, + output: + manifest = f"{TILE_DIR}/manifests/tile_merge_cats.json" + log: + f"{TILE_DIR}/logs/tile_merge_cats.json" + params: + pre = lambda wc: unit_pre("tile_merge_cats", "tile", wc.tile, + env={"NGMIX_N_CHUNKS": NGMIX_CHUNKS}), + script_hash = SCRIPT_HASH + threads: 8 + resources: + mem_mb = lambda wc, attempt: 16000 * attempt, + runtime = 120 + shell: + sp_shell("tile_merge_cats", "config_merge_sep_cats.ini") + +# The run's science product. make_cat also reads the vignette store's +# psfex_interp output, so it — not ngmix — is the store's last reader. +# +# No protected(): the full default rerun-triggers govern, and protected() only +# ever forced people through a `--forcerun` detour. +rule tile_make_cat: + group: "tile_finish" + input: + ms = rules.tile_merge_cats.output.manifest, + store = rules.tile_vignets.output.store, + output: + manifest = f"{TILE_DIR}/manifests/tile_make_cat.json", + final_cat = f"{PROD_TILE_DIR}/final_cat-{{tile}}.fits", + log: + f"{TILE_DIR}/logs/tile_make_cat.json" + params: + pre = lambda wc: unit_pre("tile_make_cat", "tile", wc.tile), + script_hash = SCRIPT_HASH + threads: 8 + resources: + mem_mb = lambda wc, attempt: 16000 * attempt, + runtime = 120 + shell: + # sp_shell's body, plus the catalogue publish: a real file, so it is a + # real declared output (and it persists — never temp()). `--job-rc` is + # composed in for the same reason it is everywhere else (see sp_shell) — + # without it a job whose counts cleared their floors but whose + # shapepipe_run died writes a "complete" log for a manifest snakemake is + # about to delete. + "{params.pre}\n" + "rc=0\n" + 'shapepipe_run -c "$SP_CONFIG/config_tile_Mc.ini" -b {threads} || rc=$?\n' + f"python {SCRIPTS}/completeness.py check tile_make_cat {{output.manifest}}" + ' --log {log} --job-rc "$rc" || rc=1\n' + "if [ $rc -eq 0 ]; then\n" + ' cp -f "$(ls -1 "$SP_RUN"/output/run_sp_Mc/make_cat_runner/output/final_cat*.fits' + ' | head -1)" {output.final_cat}\n' + "fi\n" + "exit $rc\n" diff --git a/workflow/scripts/build_forest.py b/workflow/scripts/build_forest.py new file mode 100644 index 000000000..f780531e1 --- /dev/null +++ b/workflow/scripts/build_forest.py @@ -0,0 +1,53 @@ +#!/usr/bin/env python3 +"""Build a tile's exposure symlink forest (the SP_EXP view). + +A plain script (not a run: block) so the tile chain stays group-compatible. +Reads the tile's exposures from run_index.sqlite and symlinks each exposure's +``exp///output`` into ``///output`` by exact +name (no glob). The 2-char ```` shard level is NOT cosmetic: ShapePipe's +``exp_utils.get_exp_output_files`` hardwires the sharded v2.0 layout into its +$SP_EXP glob (``///output/run_sp_*/...``), so a flat +forest makes every tile gather stage fail "No split_exp_runner output found". +(The exposure STORE is sharded the same way, for the filesystem's sake.) +The forest is a convenience view; the DAG edge to the exposures is declared in +the rule's input (tile.smk), not here. +""" + +import argparse +import shutil +import sqlite3 +from pathlib import Path + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--tile", required=True) + p.add_argument("--run-dir", required=True, type=Path) + p.add_argument("--index", required=True, type=Path) + p.add_argument("--forest", required=True, type=Path) + args = p.parse_args() + + con = sqlite3.connect(args.index, timeout=60) + exps = [r[0] for r in con.execute( + "SELECT exp_id FROM tile_exposures WHERE tile_id=?", (args.tile,))] + con.close() + + args.forest.mkdir(parents=True, exist_ok=True) + for e in exps: + src = args.run_dir / "exp" / e[:2] / e / "output" + dst = args.forest / e[:2] / e / "output" # sharded: the module glob's shape + dst.parent.mkdir(parents=True, exist_ok=True) + # A symlink (the normal case) is unlinked; a REAL directory left behind + # by a hand-run or an older layout must be removed as a tree — unlink() + # raises IsADirectoryError on it and would kill the job. + if dst.is_symlink() or dst.exists(): + if dst.is_dir() and not dst.is_symlink(): + shutil.rmtree(dst) + else: + dst.unlink() + dst.symlink_to(src) + print(f"[build_forest] {args.tile}: {len(exps)} exposures -> {args.forest}") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/build_index.py b/workflow/scripts/build_index.py new file mode 100644 index 000000000..8dc9cb13a --- /dev/null +++ b/workflow/scripts/build_index.py @@ -0,0 +1,184 @@ +#!/usr/bin/env python3 +"""Build the run index (``run_index.sqlite``) that drives the compute DAG. + +The index is *parse-time data*, never a rule input: the Snakefile loads it once +at parse time into plain dicts, so extending a run's tile list changes which +jobs exist without touching the mtime chain of completed work. + +It ACCUMULATES ACROSS INVOCATIONS, but is AUTHORITATIVE FOR THE CURRENT TILE +LIST (D1). Precisely: + + * a tile in the current list that has find_exposures output has its ``tiles`` + row replaced and its ``tile_exposures`` edges DELETED AND REBUILT — the Fe + output is the truth, so a tile whose exposure list shrank or changed must + not keep stale edges to exposures it no longer reads; + * a tile NOT in the current list is left completely untouched, which is what + makes the index span the campaign and lets a later ``clean_exposure`` see + every consuming tile; + * ``exposures`` rows only ever accumulate (``INSERT OR IGNORE``). An exposure + no tile references any more is a harmless orphan: nothing reads that table + except by joining through ``tile_exposures``. + +It records: + + tiles(tile_id, ra_dir, n_exp, status) + exposures(exp_id) -- deduplicated union over tiles + tile_exposures(tile_id, exp_id) -- the tile->exposure edges + +The tile->exposure edges are *data-derived*: they are read from each tile's +``find_exposures`` output (``exp_numbers--.txt``), which +``find_exposures_runner`` produces by parsing the tile FITS ``HISTORY`` header. +So the build is not a DAG node: ``build()`` is imported and called at PARSE TIME +by the Snakefile of the COMPUTE invocation ONLY — ``SP_PHASE == "compute"``, +which ``bin/sp`` sets — after the PREPARE invocation has produced the tiles' +find_exposures output. No other parse builds anything: a prepare parse or a +passthrough invocation (``sp --unlock``, ``sp --dag``) just loads whatever is +already on disk. (There is no ``sp index`` verb; the CLI below stays for +hand-inspection.) It +iterates the *declared* tile list and checks each tile's Fe output at its +deterministic path (no globbing — O(tile-list) existence checks, respecting the +no-``ls``-at-scale ban). A bad tile costs only that tile: it is recorded in +``missing.json`` and the index is built over the rest. The build fails only if +the missing *fraction* exceeds ``missing_threshold`` (``None`` disables the check +entirely), so a keep-going download storm that lost a few tiles does not cost the +run. The threshold is checked BEFORE anything is written, so a failed build +leaves the previous index and ``missing.json`` intact. + +Exposure IDs are stored with their trailing single-char suffix stripped +(``2243881p`` -> ``2243881``); that ``exp_base`` is the dedup key and the +exposure-rule wildcard, matching the sharded ``exp///`` store. +""" + +import argparse +import json +import sqlite3 +import sys +from pathlib import Path + + +def read_exposure_list(exp_numbers_file: Path) -> list[tuple[str, str]]: + """Return ``(exp_id, name)`` pairs from one tile's find_exposures output. + + Each line is an exposure *name* like ``2243881p``; the bare base ID (suffix + stripped) is the dedup key everywhere in the DAG, but the original name is + kept in the index — the fabricated per-unit ``exp_numbers`` list must carry + it verbatim (``get_images`` matches ``.fits.fz`` in the store; the + bare ID matches nothing). + """ + pairs = [] + for line in exp_numbers_file.read_text().splitlines(): + name = line.strip() + if not name: + continue + pairs.append((name[:-1] if name[-1].isalpha() else name, name)) + return pairs + + +def exp_list_path(run_dir: Path, tile_id: str) -> Path: + """This tile's find_exposures output, at its deterministic path. + + Sharded store (D2), fixed run dir (RUN_DATETIME=False) — an existence check, + never a glob (no ``ls`` at scale). + """ + idra, iddec = tile_id.split(".") + return (run_dir / "tiles" / tile_id[:2] / tile_id / "output" / + "run_sp_tile_Fe" / "find_exposures_runner" / "output" / + f"exp_numbers-{idra}-{iddec}.txt") + + +def build(tile_ids: list[str], run_dir: Path, db_path: Path, + missing_threshold: float | None = 0.0) -> dict: + """Build the index over ``tile_ids``; return a summary dict. + + For each tile, check its ``exp_numbers--.txt`` at the tile's + deterministic ``find_exposures`` run dir (RUN_DATETIME=False, no glob). A + tile whose exposure list is missing is recorded in ``missing.json`` and the + index is built over the rest (a bad tile costs that tile, not the run). The + build is fatal only if the missing fraction exceeds ``missing_threshold``. + + ORDER MATTERS: the threshold is evaluated FIRST, from the missing set, and + the database + ``missing.json`` are written only if it passes. A build that + aborts must leave no trace — an aborted parse that had already mutated + durable state was the bug this ordering fixes. + + The write is idempotent, which is what makes it acceptable that the compute + parse runs it even under ``-n``: re-running over an unchanged tree produces + an identical database. + """ + missing = [t for t in tile_ids if not exp_list_path(run_dir, t).exists()] + frac = len(missing) / len(tile_ids) if tile_ids else 0.0 + if missing_threshold is not None and frac > missing_threshold: + raise SystemExit( + f"Missing exposure lists for {len(missing)}/{len(tile_ids)} tile(s) " + f"(fraction {frac:.3f} > threshold {missing_threshold}): {missing}. " + f"Re-run prepare_tiles for them, or raise --missing-threshold.") + + db_path.parent.mkdir(parents=True, exist_ok=True) + con = sqlite3.connect(db_path, timeout=60) + # No DROP: the index accumulates across invocations (D1). + con.executescript( + """ + CREATE TABLE IF NOT EXISTS tiles( + tile_id TEXT PRIMARY KEY, ra_dir TEXT, n_exp INTEGER); + CREATE TABLE IF NOT EXISTS exposures( + exp_id TEXT PRIMARY KEY, name TEXT NOT NULL); + CREATE TABLE IF NOT EXISTS tile_exposures( + tile_id TEXT, exp_id TEXT, + PRIMARY KEY (tile_id, exp_id)); + """ + ) + + missing_set = set(missing) + all_exposures: set[tuple[str, str]] = set() + for tile_id in tile_ids: + if tile_id in missing_set: + continue + ra_dir = tile_id.split(".")[0] + exp_pairs = read_exposure_list(exp_list_path(run_dir, tile_id)) + con.execute("INSERT OR REPLACE INTO tiles VALUES (?,?,?)", + (tile_id, ra_dir, len(exp_pairs))) + # Replace this tile's edge set wholesale. INSERT OR IGNORE alone only + # ever added, so a tile whose exposure list shrank kept edges to + # exposures it no longer reads — and those stale edges would block it in + # the report and pin those exposures against cleanup. + con.execute("DELETE FROM tile_exposures WHERE tile_id = ?", (tile_id,)) + con.executemany("INSERT INTO tile_exposures VALUES (?,?)", + [(tile_id, exp_id) for exp_id, _ in exp_pairs]) + all_exposures.update(exp_pairs) + + # Exposures accumulate: OR IGNORE, never REPLACE (the name never changes, + # and orphans left by a shrunken tile are harmless). + con.executemany("INSERT OR IGNORE INTO exposures VALUES (?,?)", + sorted(all_exposures)) + con.commit() + con.close() + + (db_path.parent / "missing.json").write_text(json.dumps(missing, indent=2)) + return {"n_tiles": len(tile_ids) - len(missing), + "n_exposures": len(all_exposures), + "n_missing": len(missing)} + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--tile-list", required=True, type=Path, + help="file of tile IDs, one per line") + p.add_argument("--run-dir", required=True, type=Path, + help="$SP_RUN: root of the tiles/ work-dir forest") + p.add_argument("--db", required=True, type=Path, + help="output run_index.sqlite path") + p.add_argument("--missing-threshold", type=float, default=0.0, + help="fatal if the missing-tile fraction exceeds this " + "(default 0.0: any missing tile is fatal)") + args = p.parse_args() + + tile_ids = [ln.strip() for ln in args.tile_list.read_text().splitlines() + if ln.strip()] + summary = build(tile_ids, args.run_dir, args.db, args.missing_threshold) + print(f"run_index: {summary['n_tiles']} tiles, " + f"{summary['n_exposures']} exposures, " + f"{summary['n_missing']} missing -> {args.db}", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/clean_exposure.py b/workflow/scripts/clean_exposure.py new file mode 100644 index 000000000..e65c0dbcf --- /dev/null +++ b/workflow/scripts/clean_exposure.py @@ -0,0 +1,135 @@ +#!/usr/bin/env python3 +"""Reclaim ONE exposure's store and leave a tombstone (PRD #848 D5, S5). + +Run as the shell of the in-DAG ``clean_exposure`` rule, never by hand: the rule's +``input:`` is every consuming tile's ``tile_vignets`` manifest, so by the time +this executes, every campaign tile that reads this exposure has already extracted +its postage stamps. Writer, then readers, then cleaner — DAG-ordered, race-free. + +What it deletes: the exposure's whole ``output/`` tree (the bulk store — +run_sp_exp_Gie/Sp/Ma/SxSePsfPi), its ``manifests/`` and its ``logs/``, and its +star-catalogue link farms (``star_cat_exp``, plus the legacy ``star_cat_tiles``). The farms are +reclaimed for consistency, not for bytes: ``exp_star_cat``'s manifest is deleted +here like every other, so the exposure's chain must read as unbuilt, and 40 +symlinks left behind are a farm no rule now owns. The catalogue itself lives in +the run-independent cache, so rebuilding the farm costs a relink and no query. + +Deletion is SYMLINK-SAFE: a target that is itself a symlink is ``unlink``ed, not +``rmtree``d. Legacy unit dirs carry ``star_cat_exp`` as a link into the old +shared pool, and an rmtree would recurse through it and delete the shared cache +for every other exposure in the campaign. + +Deleting the manifests is deliberate and load-bearing, not tidiness: + + * the manifests are the exposure rules' DECLARED outputs. If they survived, a + tile appended later would find the exposure chain "up to date" and run + tile_vignets against products that are no longer on disk. With them gone the + DAG sees the chain as unbuilt and regenerates it — the accepted cost of a + late append (D5), expressed as ordinary Snakemake bookkeeping rather than as + a special case. + * Snakemake only demands a missing intermediate when something downstream of it + needs to run, so tiles already finished are NOT rerun by their exposures' + manifests vanishing. + +``logs/`` goes with them, and for the same reason rather than for bytes. Each +log holds the completeness verdict of one stage, written on every run and kept by +snakemake through failures; a log left behind would attest "complete" for a store +that is no longer there, contradicting the unbuilt chain the DAG must now see. +Its content for a successful stage is byte-identical to the manifest beside it, +so absorbing the logs into the tombstone would duplicate what the manifests +already carry — they are deleted, not copied. + +Nothing is lost to the report: every ``manifests/*.json`` is copied verbatim into +the tombstone under ``manifests``, and ``run_report.py`` reads a cleaned +exposure's record out of the tombstone — it reports the unit as ``cleaned``, +warn counts and shortfalls intact, instead of "not run". + +The absorption is by GLOB, so it takes whatever is in ``manifests/``, keyed by +file stem; the report re-keys on each manifest's own ``stage`` field. A legacy +``.failed.json`` from the pre-``log:`` convention is therefore carried +through unremarkably — it should never be there (an exposure with a failed stage +has no complete vignets consumer and so is not eligible for cleaning), but it +costs nothing to be right about. + +Order matters, and it is the reverse of the obvious one: the tombstone is +written FIRST, complete, and only then is anything deleted. A crash between the +two leaves a tombstone beside a store that is still there — the next invocation +treats the exposure as cleaned and only the disk is lost. Deleting first would +put the crash window where the manifests are already gone and the record that +replaces them was never written, and the report would be blind to that exposure +forever. + +The exp_psf benchmark tsv lives beside ``manifests/`` and ``logs/``, not inside +either, so it survives this job — it is the measured-memory feed for resource sizing (D4). + +The tombstone records the consumer set it was cleaned against. The rule carries +that same set as a ``params`` value, so when the index grows a new consumer the +tombstone goes stale under the default ``params`` rerun-trigger and the clean job +is rescheduled after the new tile's vignets — the exposure is cleaned once per +consumer set, not once per campaign. +""" + +import argparse +import json +import shutil +import time +from pathlib import Path + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--exp-dir", required=True, type=Path) + p.add_argument("--exp", required=True) + p.add_argument("--tombstone", required=True, type=Path) + p.add_argument("--consumers", default="", + help="comma-separated tile ids this exposure was cleaned against") + args = p.parse_args() + + consumers = [t for t in args.consumers.split(",") if t] + + # Absorb the manifests before they go: the tombstone becomes the exposure's + # surviving record. + manifests = {} + mdir = args.exp_dir / "manifests" + if mdir.is_dir(): + for f in sorted(mdir.glob("*.json")): + try: + manifests[f.stem] = json.loads(f.read_text()) + except (OSError, json.JSONDecodeError) as exc: + manifests[f.stem] = {"unreadable": str(exc)} + + # is_symlink() first, and OR'd with exists(): exists() follows the link, so a + # dangling legacy star_cat_exp would otherwise be skipped and survive. + candidates = (args.exp_dir / "output", mdir, args.exp_dir / "logs", + args.exp_dir / "star_cat_exp", args.exp_dir / "star_cat_tiles") + targets = [t for t in candidates if t.is_symlink() or t.exists()] + + # Tombstone first, complete, fsync'd — then delete. See the module docstring: + # the crash window has to sit where the data still exists, not where the + # record does not. + args.tombstone.parent.mkdir(parents=True, exist_ok=True) + tmp = args.tombstone.with_suffix(".json.tmp") + tmp.write_text(json.dumps({ + "exp": args.exp, + "cleaned_at": time.strftime("%Y-%m-%dT%H:%M:%S"), + "consumers": consumers, + "removed": [str(t) for t in targets], + "manifests": manifests, + }, indent=2) + "\n") + tmp.replace(args.tombstone) # atomic: no half-written tombstone, ever + + removed = [] + for target in targets: + # NEVER rmtree a symlink: star_cat_exp is a link into the shared pool in + # legacy unit dirs, and rmtree would follow it and empty that pool. + if target.is_symlink(): + target.unlink() + else: + shutil.rmtree(target) + removed.append(str(target)) + print(f"[clean_exposure] {args.exp}: removed {len(removed)} tree(s) after " + f"{len(consumers)} consuming tile(s)") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/completeness.py b/workflow/scripts/completeness.py new file mode 100644 index 000000000..58c1ee9eb --- /dev/null +++ b/workflow/scripts/completeness.py @@ -0,0 +1,369 @@ +#!/usr/bin/env python3 +"""The count-floor completeness table — the single failure policy. + +This is the ported ``complete_check`` count table from the v2.0 bash layer +(``run_job_sp_canfar_v2.0.bash`` job dispatch, survey §4). It is the *only* +failure policy in the design: there is no 3-class taxonomy and no error-signature +whitelist. A stage is a real failure iff a mandatory runner produced fewer than +its ``floor`` files; per-CCD attrition (a sparse CCD setools rejects, ~0.2%) +sits between ``floor`` and ``expect`` and is tolerated. + +This file is also the ``check`` CLI — the second half of every rule's shell line +(PRD D2/D3). The rules capture ShapePipe's return code rather than ``&&``-ing +onto it, so the check runs — and the verdict is recorded — even when +``shapepipe_run`` failed:: + + rc=0 + shapepipe_run -c $SP_CONFIG/config_exp_Ma.ini -b {threads} || rc=$? + completeness.py check exp_mask {output} --log {log} --job-rc "$rc" || rc=1 + exit $rc + +It counts the unit's products under ``$SP_RUN`` and exits nonzero iff a mandatory +runner is below its floor OR ``--job-rc`` is nonzero — the verdict is COMPOSED of +the counts and shapepipe_run's own exit status, because a runner can raise after +the counted ones have written their files. + +WHERE it writes the verdict is the whole point, and it is two files with two +different jobs: + + * the LOG (``--log``, the rule's snakemake ``log:``) gets the full verdict on + EVERY run, success or failure — counts against floors, per-runner detail, + scraped failure reasons, ``job_rc`` when nonzero. Snakemake never deletes a + log file, so it survives the failed job that wrote it and is the post-mortem + evidence ``run_report.py`` reads. + * the MANIFEST (the rule's declared ``output:``) gets that same verdict ONLY + when it is a success. Snakemake deletes a failed job's declared output + natively, so nothing here has to unlink anything. + +``.json`` therefore means "this stage succeeded": a resume cannot schedule +a downstream stage on top of a failed one. The removal of a PREVIOUS success's +manifest is snakemake's to do, not this script's, and the one gap that leaves is +a head process SIGKILLed between the job's failure and that deletion — a +success-named manifest then outlives the failure it no longer describes. The +profile's ``rerun-incomplete`` covers the DAG side, and ``run_report.py`` takes +the WORST status across a stage's log and manifest, so the report is right even +in that window. + +Neither file carries wall-clock, and each is rewritten ONLY when its content +changes — identical on-disk state must leave a byte-identical manifest with an +UNMOVED mtime, or the mtime rerun-trigger churns the cone on every unrelated +``--forcerun``. + +Per-runner fields: + expect nominal file count for a fully complete unit (report yardstick) + floor the fail-loud minimum (below this the job exits nonzero) + warn if True the runner never fails the unit at all (bash ``:warn`` — + e.g. psfex_interp on tiles missing some epochs) + subpath count files in ``/output//`` instead of + ``/output/`` (bash ``:rand_split`` — setools split cats) + +Counts are file counts in the runner's output dir, matching the bash +``ls / | wc -l`` semantics (broken symlinks excluded by the caller). +""" + +import argparse +import json +import os +import re +import sys +from pathlib import Path + +# stage -> {runner_subdir: {expect, floor, [warn], [subpath]}} +COMPLETENESS = { + # --- tile prepare (phase A) --- + # get_images counts are CONFIG-FLAVOR-DEPENDENT: the v2.0 bash table said 4/6 + # for the canfar vos flavor; the nibi symlink configs produce one file per + # INPUT_FILE_PATTERN entry (tile: image+weight=2; exp: image+weight+flag=3), + # verified against the p3-batch1 baseline tree (100 files / 50 tiles). + "tile_get_images": {"get_images_runner": dict(expect=2, floor=2)}, + "tile_uncompress": {"uncompress_fits_runner": dict(expect=1, floor=1)}, + "tile_find_exposures": {"find_exposures_runner": dict(expect=1, floor=1)}, + + # --- exposure chain --- + "exp_get_images": {"get_images_runner": dict(expect=3, floor=3)}, + "exp_split": {"split_exp_runner": dict(expect=121, floor=41)}, + "exp_mask": {"mask_runner": dict(expect=40, floor=1)}, + # sextractor expect is nibi-flavor: 3 files/CCD (sexcat + background + + # background_rms; v2.0's 80 assumed 2/CCD), verified against the P0 tree + # AND the bash baseline (both 120/exposure). + "exp_psf": { + "sextractor_runner": dict(expect=120, floor=2), + "setools_runner": dict(expect=80, floor=2, subpath="rand_split"), + "psfex_runner": dict(expect=80, floor=2), + "psfex_interp_runner": dict(expect=40, floor=0, warn=True), + }, + + # --- tile post --- + "tile_merge_headers": {"merge_headers_runner": dict(expect=1, floor=1)}, + "tile_mask": {"mask_runner": dict(expect=1, floor=1)}, + "tile_detect": {"sextractor_runner": dict(expect=2, floor=2)}, + "tile_vignets": { + "psfex_interp_runner": dict(expect=1, floor=1), + "vignetmaker_runner_run_1": dict(expect=1, floor=1), + # 5 sqlites/tile on nibi (image/weight/flag/background/background_rms); + # v2.0's 4 was the canfar flavor. floor follows the tile-post pattern + # (expect=floor: all-or-nothing, every vignette feeds ngmix). + "vignetmaker_runner_run_2": dict(expect=5, floor=5), + }, + "tile_ngmix": {"ngmix_runner": dict(expect=1, floor=1)}, + "tile_merge_cats": {"merge_sep_cats_runner": dict(expect=1, floor=1)}, + "tile_make_cat": {"make_cat_runner": dict(expect=1, floor=1)}, +} + + +def count_products(run_dir, runner, spec): + """Count files in ``run_dir//output[/]/`` (live links only).""" + out = run_dir / runner / "output" + if "subpath" in spec: + out = out / spec["subpath"] + if not out.is_dir(): + return 0 + return sum(1 for p in out.iterdir() if p.exists()) # p.exists() drops dead links + + +def check_floor(stage, run_dir): + """Return (ok, details). ok is False iff a mandatory runner is below floor. + + ``details`` is a list of (runner, n_found, floor, expect, warn) tuples. + """ + table = COMPLETENESS[stage] + details, ok = [], True + for runner, spec in table.items(): + n = count_products(run_dir, runner, spec) + details.append((runner, n, spec["floor"], spec["expect"], + spec.get("warn", False))) + if not spec.get("warn", False) and n < spec["floor"]: + ok = False + return ok, details + + +# --- where a stage writes ------------------------------------------------- +# +# stage -> (level, run_sp_ dir under $SP_RUN/output/). These are the +# committed configs' RUN_NAMEs (RUN_DATETIME=False makes them fixed, PRD D2), so +# the check never resolves a run-log. The ngmix entry interpolates the same env +# var its config does, so chunk K's check looks at chunk K's dir. +STAGE_DIR = { + "tile_get_images": ("tile", "run_sp_tile_Git"), + "tile_uncompress": ("tile", "run_sp_tile_Uz"), + "tile_find_exposures": ("tile", "run_sp_tile_Fe"), + "exp_get_images": ("exp", "run_sp_exp_Gie"), + "exp_split": ("exp", "run_sp_exp_Sp"), + "exp_mask": ("exp", "run_sp_exp_Ma"), + "exp_psf": ("exp", "run_sp_exp_SxSePsfPi"), + "tile_merge_headers": ("tile", "run_sp_tile_Mh_exp"), + "tile_mask": ("tile", "run_sp_tile_Ma"), + "tile_detect": ("tile", "run_sp_tile_Sx"), + "tile_detect_uc": ("tile", "run_sp_tile_Uc"), + "tile_vignets": ("tile", "run_sp_tile_PiViVi"), + "tile_ngmix": ("tile", "run_sp_tile_ngmix_Ng${SP_NGMIX_CHUNK}u"), + "tile_merge_cats": ("tile", "run_sp_Ms"), + "tile_make_cat": ("tile", "run_sp_Mc"), +} + + +# --- failure reasons ------------------------------------------------------ + +# Lines worth showing a human who asks "why is this runner short?". Deliberately +# crude: the point is a pointer into the logs, not a taxonomy (there is no error +# whitelist in this design — the count floor is the policy). +_ERROR_RE = re.compile( + r"traceback|exception|\berror\b|\bfailed\b|no such file|not found|" + r"killed|out of memory|oom|segmentation fault|bad chi2", + re.IGNORECASE) +_TS_RE = re.compile(r"^\d{2}/\d{2}/\d{4} \d{2}:\d{2}:\d{2}\s*") +_NOISE_RE = re.compile(r"A total of 0 errors were recorded") + +MAX_LOG_FILES = 40 # logs are per-CCD; a handful is enough to characterise +MAX_TAIL_LINES = 120 # per file +MAX_REASONS = 3 # per runner + + +def _normalise(line: str) -> str: + """Collapse a log line to its shape, so 40 per-CCD copies dedupe to one.""" + line = _TS_RE.sub("", line.strip()) + line = re.sub(r"/\S+", "", line) # paths differ per CCD + line = re.sub(r"\d+", "N", line) + return line[:200] + + +def scrape_reasons(stage_dir, runner): + """Best-effort, bounded: distinct error-looking lines from a runner's logs. + + Two sources, in order of usefulness: the runner's per-process worker logs + (``/logs/process-*.log`` — where the module's own exception lands), + and the stage's ``logs/log_sp.log`` (where ShapePipe records its error + tally). Sorted, truncated, deduped by shape — a manifest must stay + byte-stable for a given tree. + """ + seen, reasons = {}, [] + candidates = [] + for d in (stage_dir / runner / "logs", stage_dir / "logs"): + if d.is_dir(): + candidates += sorted(p for p in d.iterdir() if p.is_file()) + for path in candidates[:MAX_LOG_FILES]: + try: + lines = path.read_text(errors="replace").splitlines()[-MAX_TAIL_LINES:] + except OSError: + continue + for raw in lines: + if not _ERROR_RE.search(raw) or _NOISE_RE.search(raw): + continue + shape = _normalise(raw) + if shape in seen: + seen[shape] += 1 + continue + seen[shape] = 1 + reasons.append([path.name, _TS_RE.sub("", raw.strip())[:300], shape]) + out = [] + for name, text, shape in reasons[:MAX_REASONS]: + n = seen[shape] + out.append(f"{name}: {text}" + (f" [x{n}]" if n > 1 else "")) + return out + + +# --- manifest ------------------------------------------------------------- + +def build_manifest(stage, run_dir, unit, stage_subdir=None): + """Count, classify and (on shortfall) scrape. Returns (manifest, ok). + + Stages absent from the table fall back to the zero-output floor: any product + anywhere under the stage dir passes, nothing at all fails. + """ + level, subdir = STAGE_DIR.get(stage, (None, None)) + subdir = stage_subdir or (os.path.expandvars(subdir) if subdir else None) + stage_dir = run_dir / "output" / subdir if subdir else run_dir + manifest = { + "stage": stage, + "level": level, + "unit": unit, + "run_dir": str(run_dir), + "stage_dir": str(stage_dir), + "runners": {}, + "failures": [], + } + + if stage not in COMPLETENESS: + produced = list(stage_dir.glob("**/output/*")) if stage_dir.is_dir() else [] + ok = bool(produced) + manifest["status"] = "complete" if ok else "failed" + manifest["n_products"] = len(produced) + if not ok: + manifest["failures"].append( + {"runner": None, "found": 0, "floor": 1, + "reasons": [f"zero output under {stage_dir}"]}) + return manifest, ok + + ok, details = check_floor(stage, stage_dir) + short = False + for runner, n, floor, expect, warn in details: + below = n < floor + if n < expect: + short = True + manifest["runners"][runner] = { + "found": n, "expect": expect, "floor": floor, "warn": warn, + "status": ("complete" if n >= expect else + "warn" if (warn or not below) else "below_floor"), + } + if below: + manifest["failures"].append({ + "runner": runner, "found": n, "floor": floor, "expect": expect, + "warn": warn, "reasons": scrape_reasons(stage_dir, runner), + }) + manifest["status"] = "failed" if not ok else ("warn" if short else "complete") + return manifest, ok + + +def write_if_changed(path: Path, text: str) -> None: + """Write only when the bytes differ — see the module docstring on mtime.""" + path.parent.mkdir(parents=True, exist_ok=True) + if not path.exists() or path.read_text() != text: + path.write_text(text) + + +def _unit_from_run_dir(run_dir): + """The human unit ID: the basename of ``$SP_RUN`` (``210.282``, ``2605805``). + + NOT ``SP_UNIT_NUM``, which carries ShapePipe's dashed numbering form + (``-210-282``) and would put ``210-282`` in the manifest — a key that joins + to nothing. ``run_report`` keys units on the store directory name, which is + exactly this basename, so the two now agree. + """ + return Path(str(run_dir)).name or "unknown" + + +def main(argv=None) -> int: + p = argparse.ArgumentParser(description="ShapePipe per-unit completeness check") + sub = p.add_subparsers(dest="cmd", required=True) + c = sub.add_parser("check", help="count products, write the manifest") + c.add_argument("stage") + c.add_argument("manifest", type=Path) + c.add_argument("--log", type=Path, required=True, + help="the rule's log: path — the verdict is written here every " + "run, success or failure") + c.add_argument("--run-dir", type=Path, default=None, + help="the unit's $SP_RUN (default: the env var)") + c.add_argument("--unit", default=None, + help="override the unit ID (default: basename of $SP_RUN)") + c.add_argument("--stage-dir", default=None, + help="override the run_sp_* subdir (default: the stage table)") + c.add_argument("--job-rc", type=int, default=0, + help="shapepipe_run's exit status, composed into the verdict") + args = p.parse_args(argv) + + run_dir = args.run_dir or Path(os.environ.get("SP_RUN", "")) + if not str(run_dir): + print("[completeness] FATAL: $SP_RUN unset and --run-dir not given", + file=sys.stderr) + return 2 + unit = args.unit or _unit_from_run_dir(run_dir) + + manifest, ok = build_manifest(args.stage, Path(run_dir), unit, args.stage_dir) + + # The verdict is COMPOSED of two independent statements: the count floors + # (above) and shapepipe_run's own exit status (here). Counts alone are not + # enough — a runner can raise AFTER the counted runners have written their + # files, so the floors are met while the job died. Without this, such a job + # publishes a SUCCESS manifest that snakemake then deletes as a failed job's + # output — and the log would claim "complete" for a stage with no manifest, + # which reads as a bookkeeping bug rather than as the failure it is. + # + # Recorded only when nonzero, which keeps every existing success manifest + # byte-identical (the mtime rerun-trigger reads those bytes). + if args.job_rc != 0: + ok = False + manifest["status"] = "failed" + manifest["job_rc"] = args.job_rc + manifest["failures"].append({ + "runner": "shapepipe_run", "found": 0, "floor": 1, "expect": 1, + "warn": False, + "reasons": [f"shapepipe_run exited {args.job_rc} " + f"(counts above floor)"], + }) + + text = json.dumps(manifest, indent=2, sort_keys=True) + "\n" + + # The log ALWAYS gets the verdict; the manifest gets it only on success. No + # unlink of anything: snakemake deletes a failed job's declared output, so + # the presence of .json IS the statement "this stage succeeded" and + # the DAG can never build on top of a failure. A failure->success transition + # publishes the manifest; a success->failure has the manifest removed for us, + # and the log is overwritten with the new verdict either way. + write_if_changed(args.log, text) + if ok: + write_if_changed(args.manifest, text) + + for runner, r in manifest["runners"].items(): + tag = {"complete": "OK", "warn": "warn", "below_floor": "<-- BELOW floor"} + print(f"[completeness] {runner}: {r['found']}/{r['expect']} " + f"(floor {r['floor']}) {tag[r['status']]}", file=sys.stderr) + print(f"[completeness] {args.stage} {unit}: {manifest['status']} " + f"-> {args.log}" + (f" + {args.manifest}" if ok else ""), file=sys.stderr) + for f in manifest["failures"]: + for reason in f["reasons"]: + print(f"[completeness] {f['runner']}: {reason}", file=sys.stderr) + return 0 if ok else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/workflow/scripts/container.py b/workflow/scripts/container.py new file mode 100644 index 000000000..bbd7e0bb9 --- /dev/null +++ b/workflow/scripts/container.py @@ -0,0 +1,500 @@ +#!/usr/bin/env python3 +"""Manage this user's copy of the ShapePipe container image. + +Two layers, and the second only exists when you ask for one: + +* the **SIF** (``~/.cache/shapepipe/shapepipe.sif``) -- a pristine, read-only + copy of the published image, pulled into your own cache. Per-user by + construction: one file, one owner, nobody else's refresh moves the ground + under a running job. +* an optional **sandbox** (``~/.cache/shapepipe/sandbox/``) -- the same image + unpacked into a writable directory, so a ``pip install`` inside it sticks. + The escape hatch for work that needs a package the image does not carry yet. + +Resolution order, shared by this CLI and by the workflow: **sandbox if it +exists, else the cached SIF if it exists, else the ``container:`` path in +workflow/config.yaml**. That last one is the current shared /project image, so +a checkout with an empty cache behaves exactly as it did before this verb +existed, and a package installed into your sandbox is there for your workflow +jobs too. + +Subcommands, exposed as ``sp container ``:: + + sp container pull # fetch the tag into the cache + sp container status # what is here, and how current it is + sp container sandbox # unpack the SIF into a writable dir + sp container exec # run something inside it + sp container exec --writable # ... with writes that persist + +``pull`` needs the network. Compute nodes on Alliance clusters generally have +none, so run it on a login node or inside an ``salloc`` allocation -- never +from a batch job. + +Deliberately **stdlib-only**: it runs on the bare host, outside the container, +where the science stack is not installed, and so must import without it. +""" + +import argparse +import os +import re +import shlex +import shutil +import subprocess +import sys +from pathlib import Path + + +class ContainerError(Exception): + """A misconfiguration the caller must fix (bad override, no image at all). + + Raised rather than ``sys.exit``ed so the Snakefile, which imports this + module at parse time, can turn it into a WorkflowError instead of a + SystemExit. ``main`` below turns it back into a one-line CLI error. + """ + +# The published image. CI pushes one tag per branch, sanitized; `-runtime` is +# the slim variant the workflow runs. +CONTAINER_URI = "docker://ghcr.io/cosmostat/shapepipe:develop-runtime" + +# The single source of truth for the fallback image: the workflow's own +# `container:` key, which is also what the Snakefile reads. Written down once, +# here, so the CLI and the workflow cannot disagree about the default. +CONFIG_FILE = Path(__file__).resolve().parents[1] / "config.yaml" +CONFIG_KEY = "container" + +# The profile whose `apptainer-args:` every workflow job runs under. `exec` reads +# it at runtime rather than restating it, so a one-off `sp container exec` and a +# job see the same environment (the PYTHONPATH pin above all: a divergence there +# means the one-off imports a different src/ than the workflow does). +PROFILE_FILE = Path(__file__).resolve().parents[2] / "profiles" / "nibi" / "config.yaml" + +# ~/.cache/shapepipe by default; SP_CACHE_DIR moves the whole cache (e.g. onto +# a filesystem with room), XDG_CACHE_HOME moves it with everything else. +CACHE_DIR = Path( + os.environ.get("SP_CACHE_DIR") + or Path(os.environ.get("XDG_CACHE_HOME", "~/.cache")) / "shapepipe" +).expanduser() + +# This user's pristine image. Override with ``SP_CONTAINER`` (absolute path). +DEFAULT_SIF = CACHE_DIR / "shapepipe.sif" + +# The optional writable unpacking of it. Override with ``SP_SANDBOX``. +DEFAULT_SANDBOX = CACHE_DIR / "sandbox" + +# Bind mounts for `exec`, matching the nibi profile's apptainer-args (the two +# cluster filesystems this workflow reads and writes, plus the home that holds +# ~/.ssl/cadcproxy.pem). Override wholesale with ``SP_APPTAINER_BINDS``. +DEFAULT_BINDS = "/project,/scratch,/home" + + +def configured_default(): + """Return the ``container:`` path from workflow/config.yaml, or ``None``. + + A deliberately minimal scalar read (the same one bin/sp does in sed): this + module is stdlib-only, so there is no yaml to import. + """ + try: + text = CONFIG_FILE.read_text() + except OSError: + return None + match = re.search(rf"^{CONFIG_KEY}:[ \t]*(\S+)", text, re.MULTILINE) + return match.group(1) if match else None + + +def profile_apptainer_args(): + """Return the profile's ``apptainer-args`` as a token list, or ``[]``. + + Minimal scalar read again (stdlib-only); a missing or unreadable profile + yields ``[]``, which callers replace with their own defaults. + """ + try: + text = PROFILE_FILE.read_text() + except OSError: + return [] + match = re.search(r'^apptainer-args:[ \t]*"(.*)"[ \t]*$', text, re.MULTILINE) + return shlex.split(match.group(1)) if match else [] + + +def local_sif(): + """Return this user's cached image path (may not exist yet).""" + override = os.environ.get("SP_CONTAINER") + return (Path(override) if override else DEFAULT_SIF).expanduser() + + +def local_sandbox(): + """Return this user's writable sandbox directory (may not exist).""" + override = os.environ.get("SP_SANDBOX") + return (Path(override) if override else DEFAULT_SANDBOX).expanduser() + + +def resolve_image(): + """Return ``(path, kind)`` for the image everything should run. + + ``kind`` is ``"sandbox"``, ``"sif"``, ``"configured"`` (the shared + /project image named in config.yaml -- the default when the cache is + empty) or ``"none"``. + """ + sandbox = local_sandbox() + if sandbox.is_dir(): + return str(sandbox), "sandbox" + sif = local_sif() + if sif.exists(): + return str(sif), "sif" + # An override that resolves to nothing is a typo, never an intention: falling + # through to the shared image would run the job against something other than + # what the user named. + if os.environ.get("SP_CONTAINER"): + raise ContainerError( + f"SP_CONTAINER={os.environ['SP_CONTAINER']} does not exist " + f"(resolved to {sif}). Unset it, or point it at an image that does." + ) + default = configured_default() + if default: + return default, "configured" + return "", "none" + + +def image_labels(image): + """Return the image's OCI labels as a dict, or ``{}`` if unreadable. + + Never raises: a missing file, a missing ``apptainer`` or a corrupt image + all mean "we don't know", which every caller treats as non-fatal. + """ + path = Path(image) + if not path.exists() or shutil.which("apptainer") is None: + return {} + try: + out = subprocess.run( + ["apptainer", "inspect", "--labels", str(path)], + capture_output=True, + text=True, + timeout=60, + ) + except (OSError, subprocess.SubprocessError): + return {} + if out.returncode != 0: + return {} + labels = {} + for line in out.stdout.splitlines(): + key, sep, value = line.partition(":") + if sep: + labels[key.strip()] = value.strip() + return labels + + +def image_revision(image): + """Return the commit the image was built from, or ``None``.""" + return image_labels(image).get("org.opencontainers.image.revision") + + +def _require_apptainer(): + """Exit unless ``apptainer`` is on PATH (bin/sp loads the module).""" + if shutil.which("apptainer") is None: + sys.exit("apptainer is not on PATH (bin/sp loads apptainer/1.4.5)") + + +def _git(*args, cwd=None): + """Run git, returning stripped stdout or ``None`` on any failure.""" + try: + out = subprocess.run( + ["git", *args], capture_output=True, text=True, cwd=cwd, timeout=30 + ) + except (OSError, subprocess.SubprocessError): + return None + return out.stdout.strip() if out.returncode == 0 else None + + +def compare_revision(revision, repo=None): + """Place an image revision relative to this checkout's HEAD. + + One of ``"in-sync"``, ``"behind"`` (the image predates HEAD), ``"ahead"``, + ``"diverged"``, or ``"unknown"`` (no label, no git, or a commit this clone + has never fetched). + """ + if not revision: + return "unknown" + repo = repo or Path(__file__).resolve().parents[2] + head = _git("rev-parse", "HEAD", cwd=repo) + if head is None: + return "unknown" + if head == revision: + return "in-sync" + if _git("cat-file", "-e", f"{revision}^{{commit}}", cwd=repo) is None: + return "unknown" + if _git("merge-base", "--is-ancestor", revision, head, cwd=repo) is not None: + return "behind" + if _git("merge-base", "--is-ancestor", head, revision, cwd=repo) is not None: + return "ahead" + return "diverged" + + +def cmd_pull(args): + """Pull ``--tag`` into the cache, atomically.""" + _require_apptainer() + sif = local_sif() + sif.parent.mkdir(parents=True, exist_ok=True) + # Pull to a sibling temp name and rename: an atomic rename within one + # directory, so an in-flight job sees either the whole old image or the + # whole new one. Pulling in place leaves the file half-written for the many + # minutes the pull takes. Jobs already running hold the old inode open. + tmp = sif.with_name(sif.name + f".pull.{os.getpid()}") + print(f"pulling {args.tag}\n -> {sif}") + try: + subprocess.run( + ["apptainer", "pull", "--force", "--name", str(tmp), args.tag], check=True + ) + os.replace(tmp, sif) + except subprocess.CalledProcessError as exc: + tmp.unlink(missing_ok=True) + sys.exit(f"pull failed ({exc.returncode}); {sif} is unchanged") + except KeyboardInterrupt: + tmp.unlink(missing_ok=True) + raise + labels = image_labels(sif) + print(f"revision: {labels.get('org.opencontainers.image.revision', 'unknown')}") + print(f"version: {labels.get('org.opencontainers.image.version', 'unknown')}") + return 0 + + +def cmd_sandbox(args): + """Unpack the image into a writable directory -- the opt-in escape hatch.""" + _require_apptainer() + sandbox = local_sandbox() + if sandbox.exists() and not args.force: + sys.exit( + f"sandbox already exists at {sandbox}\n" + "pass --force to discard it and rebuild from a clean image" + ) + source = args.source + if not source: + image, kind = resolve_image() + if kind == "sandbox": + # Rebuilding from the sandbox itself would just re-copy the drift. + image = str(local_sif()) if local_sif().exists() else ( + configured_default() or CONTAINER_URI + ) + source = image or CONTAINER_URI + sandbox.parent.mkdir(parents=True, exist_ok=True) + print(f"building sandbox from {source}\n -> {sandbox}") + # Build beside the target and swap it in, as `pull` does -- and for a + # sharper reason. A half-written .sif fails loudly, but a half-unpacked + # sandbox *directory* is still a directory, so resolve_image() would elect + # it and every job would silently run a broken tree. Staging also means a + # --force rebuild that fails leaves the sandbox you already had intact. + # + # `--fix-perms` so the tree can be deleted again later. No `--fakeroot`: an + # unprivileged build from an existing image goes through user namespaces, + # which is what the Alliance clusters provide. + staging = sandbox.with_name(f"{sandbox.name}.build.{os.getpid()}") + shutil.rmtree(staging, ignore_errors=True) + try: + subprocess.run( + ["apptainer", "build", "--sandbox", "--fix-perms", str(staging), source], + check=True, + ) + except subprocess.CalledProcessError as exc: + shutil.rmtree(staging, ignore_errors=True) + sys.exit(f"sandbox build failed ({exc.returncode}); {sandbox} is unchanged") + except (KeyboardInterrupt, OSError): + shutil.rmtree(staging, ignore_errors=True) + raise + + if sandbox.exists(): + print(f"replacing {sandbox}") + shutil.rmtree(sandbox, ignore_errors=True) + if sandbox.exists(): + shutil.rmtree(staging, ignore_errors=True) + sys.exit(f"could not remove {sandbox}; remove it by hand and retry") + os.replace(staging, sandbox) + print( + "\nthis sandbox now takes precedence over the SIF everywhere, including " + "workflow jobs.\ninstall into it with: sp container exec --writable pip " + "install \nreset to a clean image with: sp container pull && " + "sp container sandbox --force" + ) + return 0 + + +def cmd_status(args): + """Report which image layer is live, its revision, and how current it is.""" + sif = local_sif() + sandbox = local_sandbox() + # status is the verb you run WHEN something is wrong, so a broken override is + # reported here rather than raised. + try: + active, kind = resolve_image() + except ContainerError as exc: + print(f"active: NONE -- {exc}") + return 1 + + if sif.exists(): + print(f"SIF: {sif} ({sif.stat().st_size / 1e9:.1f} GB)") + else: + print(f"SIF: absent ({sif})") + if sandbox.is_dir(): + print(f"sandbox: {sandbox} (writable; may carry local modifications)") + else: + print("sandbox: none") + print(f"configured: {configured_default() or 'unset'} ({CONFIG_FILE})") + + if kind == "none": + print("\nactive: NONE -- no cached image and no container: in config.yaml") + print(f"run: sp container pull --tag {CONTAINER_URI}") + return 1 + + print(f"\nactive: {active} ({kind})") + if kind == "configured" and not Path(active).exists(): + print(" WARNING: that path does not exist on this host") + return 1 + + labels = image_labels(active) + revision = labels.get("org.opencontainers.image.revision") + source = "" + if revision is None and kind == "sandbox" and sif.exists(): + # Some sandbox trees do not carry the original labels through. The SIF + # beside it is the best remaining evidence of what it was built from -- + # a guess, so it is labelled as one rather than printed as fact. + revision = image_revision(sif) + if revision: + source = " (inferred from the SIF beside it, not read from the sandbox)" + print(f"revision: {revision or 'unknown'}{source}") + print(f"version: {labels.get('org.opencontainers.image.version', 'unknown')}") + if kind == "sandbox": + print( + " (that revision is what the sandbox was built from; " + "anything\n installed into it since is in no label)" + ) + verdict = compare_revision(revision) + explain = { + "in-sync": "matches this checkout's HEAD", + "behind": ( + "older than this checkout's HEAD -- `sp container pull` fetches " + f"{CONTAINER_URI}, which is not built from this branch unless you " + "pass --tag" + ), + "ahead": "newer than this checkout's HEAD", + "diverged": "on a different branch from this checkout", + "unknown": "cannot compare (no label, or a commit this clone lacks)", + }[verdict] + print(f"checkout: {verdict} ({explain})") + return 0 + + +def cmd_exec(args): + """Run a command inside the resolved image -- the one-off path.""" + _require_apptainer() + command = [a for a in args.command if a != "--"] + if not command: + sys.exit("nothing to run; pass a command after `exec`") + + # Same environment the workflow's jobs get: the profile's apptainer-args + # verbatim (--cleanenv, the PYTHONPATH pin, --home, its binds). Explicit + # binds still win, and a profile that cannot be read falls back to the old + # standalone defaults. + env_args = profile_apptainer_args() + explicit_binds = args.bind or os.environ.get("SP_APPTAINER_BINDS") + if not env_args: + print( + f"warning: could not read apptainer-args from {PROFILE_FILE}; " + "falling back to --cleanenv and the default binds, which may not " + "match what jobs run under", + file=sys.stderr, + ) + env_args = ["--cleanenv"] + explicit_binds = explicit_binds or DEFAULT_BINDS + if explicit_binds: + env_args += ["--bind", explicit_binds] + + if args.writable: + # A SIF is a read-only filesystem, so `--writable` against one fails + # obscurely; only a sandbox takes writes. + sandbox = local_sandbox() + if not sandbox.is_dir(): + sys.exit( + f"--writable needs a sandbox, and there is none at {sandbox}\n" + "build one with: sp container sandbox" + ) + image, extra = str(sandbox), ["--writable"] + else: + image, kind = resolve_image() + if kind == "none": + sys.exit("no image resolved; run: sp container pull") + extra = [] + + cmd = ["apptainer", "exec", *extra, *env_args, image, *command] + return subprocess.run(cmd).returncode + + +def cmd_resolve(args): + """Print just the resolved image path -- what the Snakefile consumes.""" + image, kind = resolve_image() + if kind == "none": + sys.exit("no image resolved; run: sp container pull") + print(image) + return 0 + + +def build_parser(): + parser = argparse.ArgumentParser( + prog="sp container", + description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + sub = parser.add_subparsers(dest="subcommand", required=True) + + p_pull = sub.add_parser( + "pull", + help="fetch the image into the cache (login node or salloc: compute " + "nodes have no network)", + ) + p_pull.add_argument( + "--tag", default=CONTAINER_URI, help=f"image to pull (default: {CONTAINER_URI})" + ) + p_pull.set_defaults(func=cmd_pull) + + p_status = sub.add_parser( + "status", help="report which image layer is live and how current it is" + ) + p_status.set_defaults(func=cmd_status) + + p_sandbox = sub.add_parser( + "sandbox", help="unpack the image into a writable directory (opt-in)" + ) + p_sandbox.add_argument( + "--source", help="image to unpack (default: the cached SIF, or the config path)" + ) + p_sandbox.add_argument( + "--force", + action="store_true", + help="discard an existing sandbox and rebuild from a clean image", + ) + p_sandbox.set_defaults(func=cmd_sandbox) + + p_exec = sub.add_parser("exec", help="run a command inside the resolved image") + p_exec.add_argument("--bind", help=f"bind mounts (default: {DEFAULT_BINDS})") + p_exec.add_argument( + "--writable", + action="store_true", + help="run against the sandbox so writes (e.g. pip install) persist", + ) + p_exec.add_argument("command", nargs=argparse.REMAINDER) + p_exec.set_defaults(func=cmd_exec) + + p_resolve = sub.add_parser( + "resolve", help="print the resolved image path (what the workflow runs)" + ) + p_resolve.set_defaults(func=cmd_resolve) + + return parser + + +def main(argv=None): + args = build_parser().parse_args(argv) + try: + return args.func(args) + except ContainerError as exc: + sys.exit(str(exc)) + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/workflow/scripts/ngmix_range.py b/workflow/scripts/ngmix_range.py new file mode 100644 index 000000000..9e3058052 --- /dev/null +++ b/workflow/scripts/ngmix_range.py @@ -0,0 +1,55 @@ +#!/usr/bin/env python3 +"""Print one ngmix chunk's closed object-ID range as bash exports. + +Run inside the chunk's shell, from the tile's own sexcat, because the range is +only knowable at execution time (PRD D4):: + + eval "$(ngmix_range.py --run-dir $SP_RUN --chunk 3 --n-chunks 8)" + # -> export NGMIX_ID_MIN=751; export NGMIX_ID_MAX=1125 + +SExtractor's NUMBER column (ngmix's obj_id) runs 1..N contiguous, so covering +[1, N] processes every object exactly once. The first N-1 chunks get equal +shares; the last takes the remainder, with a CLOSED upper bound at n_obj — +never ID_OBJ_MAX = -1, which ngmix reads as unbounded and would double-count. +Each tile is sized from its OWN count (the bash monolith used the campaign +average and left the last chunk open). +""" + +import argparse +from pathlib import Path + + +def id_ranges(n_obj: int, n_chunks: int) -> list[tuple[int, int]]: + base, rem = divmod(n_obj, n_chunks) + ranges, lo = [], 1 + for k in range(1, n_chunks + 1): + hi = lo + base + (rem if k == n_chunks else 0) - 1 + ranges.append((lo, hi)) + lo = hi + 1 + return ranges + + +def object_count(run_dir: Path) -> int: + """NAXIS2 of the last HDU of this tile's sexcat (get_number_objects.py).""" + from astropy.io import fits + + cats = sorted((run_dir / "output" / "run_sp_tile_Sx").glob( + "sextractor_runner/output/sexcat*.fits")) + if not cats: + raise SystemExit(f"[ngmix_range] FATAL: no sexcat under {run_dir}") + with fits.open(cats[0]) as hdul: + return int(hdul[-1].header["NAXIS2"]) + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--run-dir", required=True, type=Path) + p.add_argument("--chunk", required=True, type=int) + p.add_argument("--n-chunks", required=True, type=int) + a = p.parse_args() + lo, hi = id_ranges(object_count(a.run_dir), a.n_chunks)[a.chunk - 1] + print(f"export NGMIX_ID_MIN={lo}; export NGMIX_ID_MAX={hi}") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/run_report.py b/workflow/scripts/run_report.py new file mode 100644 index 000000000..22201f7e7 --- /dev/null +++ b/workflow/scripts/run_report.py @@ -0,0 +1,333 @@ +#!/usr/bin/env python3 +"""``sp report`` — the run's success/failure tables, read from the manifests. + +NOT a DAG node. A report rule that declared all tiles' outputs as inputs would be +a descendant of every job, so one hard failure under --keep-going would poison +its cone and the report would never run — the exact scenario it exists for. So it +is a plain script, runnable at any time including mid-run; the Snakefile's +onsuccess/onerror hooks call it so every invocation ends with one. + +It reads two things and nothing else (PRD D3): + + * the **index** (``run_index.sqlite``) — the units the run declared, and the + tile->exposure edges that let an exposure failure be blamed on the tiles it + blocks; + * the **verdicts** written by ``completeness.py check`` — per-runner + found/expect/floor and log-scraped failure reasons — which land in two places + with two different lifetimes. Every run writes the unit's ``logs/.json`` + (the rule's snakemake ``log:``, which snakemake never deletes); a run whose + verdict is a success ADDITIONALLY writes ``manifests/.json``, the + rule's declared output, which snakemake deletes when the job fails. + +Both dirs are read here and the status comes from the file BODY, so a failed +stage speaks through its log while a successful one is corroborated by two +identical files. A unit with neither ran nothing. Also read, for continuity with +stores written before this convention: ``manifests/.failed.json``, which +is where failure evidence used to go. + +A reclaimed exposure has neither dir: ``clean_exposure`` deleted both, after +copying the manifests into ``/cleaned.json``. That tombstone is read as the +unit's record and the unit is reported as **cleaned** — not "not run", and it +blocks no tile. + +No disk scanning: counting products is the *check's* job, done once at the moment +the products were fresh. A unit with no manifest for a stage is "not run" — which +is a real and distinct answer from "ran and produced nothing". + +Records are discovered by glob (``tiles/**/manifests/*.json`` and +``tiles/**/logs/*.json``), not by constructed path: the unit stores are sharded +(``tiles///``) and the sharding depth is not this script's business. +""" + +import argparse +import json +import sqlite3 +import sys +from collections import defaultdict +from pathlib import Path + +# Stage order per level — the report's column order, and the definition of +# "expected" (a declared unit with no manifest for one of these is not run). +TILE_STAGES = ["tile_get_images", "tile_uncompress", "tile_find_exposures", + "tile_merge_headers", "tile_detect", "tile_vignets", + "tile_ngmix", "tile_merge_cats", "tile_make_cat"] +EXP_STAGES = ["exp_get_images", "exp_star_cat", "exp_split", "exp_mask", "exp_psf"] + +STATUSES = ("complete", "warn", "failed", "not_run") + + +def load_manifests(run_dir: Path, sub: str) -> dict: + """``{unit: {stage: verdict}}`` for one store (``tiles`` or ``exp``). + + Reads BOTH of a unit's record dirs: ``manifests/`` (the rules' declared + outputs, success-only) and ``logs/`` (the rules' ``log:``, written every run + and never deleted by snakemake, so this is where a failure survives). + + The unit key is the record dir's *parent directory name* — shard-depth + agnostic, and the only form that joins to the index (the record's own + ``unit`` field carries ``SP_UNIT_NUM``'s dashed form, ``210-282``, which is + not the index's ``210.282``). The stage comes from the body, never the + filename: ngmix chunks share a stage under per-chunk filenames, and a log + names the same stage as the manifest beside it. + + Several files therefore map to one (unit, stage), and the WORST status wins. + That is what collapses the ngmix chunks to one entry, and it is why a + successful stage's two byte-identical records cost nothing while a failure + always speaks. A body with no ``stage`` field is skipped: it is not one of + ours, which is what keeps a stray JSON deeper in the tree inert. + """ + out: dict = defaultdict(dict) + paths = sorted((run_dir / sub).glob("**/manifests/*.json")) \ + + sorted((run_dir / sub).glob("**/logs/*.json")) + for path in paths: + try: + m = json.loads(path.read_text()) + except (OSError, json.JSONDecodeError) as exc: + print(f"[run_report] unreadable record {path}: {exc}", file=sys.stderr) + continue + if not isinstance(m, dict) or "stage" not in m: + continue + unit = path.parent.parent.name + stage = m["stage"] + prev = out[unit].get(stage) + # Worst status wins when several records share a stage (ngmix chunks; + # a stage's manifest and its identical log). + rank = lambda d: STATUSES.index(d.get("status")) if d.get("status") in STATUSES else len(STATUSES) # noqa: E731 + if prev is None or rank(m) > rank(prev): + out[unit][stage] = m + return out + + +def absorb_tombstones(run_dir: Path, sub: str, manifests: dict) -> set: + """Fill in reclaimed units from their ``cleaned.json``; return their ids. + + A cleaned exposure has neither ``manifests/`` nor ``logs/`` — + ``clean_exposure`` deleted both, after copying every manifest verbatim into + the tombstone (the logs duplicate them). Read them back, or + the report inverts the truth exactly when reclamation works: the exposure + shows as "not run" and blocks the very tiles whose completion authorised the + deletion. + + Manifests on disk win if both exist — that is a re-built chain, and the + tombstone is then a stale record of the previous generation. + """ + cleaned = set() + for path in sorted((run_dir / sub).glob("**/cleaned.json")): + unit = path.parent.name + try: + tomb = json.loads(path.read_text()) + except (OSError, json.JSONDecodeError) as exc: + print(f"[run_report] unreadable tombstone {path}: {exc}", file=sys.stderr) + continue + if manifests.get(unit): + continue + for key, m in (tomb.get("manifests") or {}).items(): + if not isinstance(m, dict): + continue + manifests[unit][m.get("stage", key)] = m + cleaned.add(unit) + return cleaned + + +def shortfalls(m: dict) -> dict: + """``{runner: (found, expect, floor)}`` for every runner under expect.""" + return {r: (d["found"], d["expect"], d["floor"]) + for r, d in m.get("runners", {}).items() if d["found"] < d["expect"]} + + +def reasons(m: dict) -> list: + """Flattened failure reasons, runner-tagged, for the report's why column.""" + out = [] + for f in m.get("failures", []): + head = f"{f['runner']} {f['found']}/{f.get('expect', '?')} (floor {f['floor']})" + out += [f"{head}: {r}" for r in f["reasons"]] or [head] + return out + + +def tally_level(units, stages, manifests, cleaned=frozenset()) -> dict: + """Per-stage counts + named unit lists, for one level. + + ``cleaned`` units are counted by the status their absorbed manifests carry + (complete or warn — attrition is preserved) and additionally listed under + ``cleaned``, so a reclaimed campaign reads as reclaimed rather than as a + campaign that never ran. + """ + per_stage = {} + for stage in stages: + t = {"complete": 0, "warn": [], "failed": [], "not_run": [], "cleaned": []} + agg = defaultdict(lambda: {"found": 0, "expect": 0, "by_unit": {}}) + for u in units: + m = manifests.get(u, {}).get(stage) + if m is None: + t["not_run"].append(u) + continue + if u in cleaned: + t["cleaned"].append(u) + status = m.get("status", "failed") + status = status if status in ("complete", "warn") else "failed" + if status == "complete": + t["complete"] += 1 + else: + t[status].append(u) + if status == "failed": + # Failed units are named above, never folded into the attrition + # aggregate: a whole-unit failure is not per-CCD attrition, and + # mixing them hides real deletion bugs behind a big denominator. + continue + for runner, d in m.get("runners", {}).items(): + a = agg[runner] + a["found"] += d["found"] + a["expect"] += d["expect"] + if d["found"] < d["expect"]: + a["by_unit"][u] = d["expect"] - d["found"] + for a in agg.values(): + if not a["by_unit"]: + del a["by_unit"] + t["products"] = dict(agg) + per_stage[stage] = t + return per_stage + + +def unit_rows(units, stages, manifests) -> list: + """One row per non-clean unit: its first bad stage, shortfalls, why.""" + rows = [] + for u in units: + got = manifests.get(u, {}) + bad = [s for s in stages + if got.get(s) is None or got[s].get("status") != "complete"] + if not bad: + continue + stage = bad[0] + m = got.get(stage) + rows.append({ + "unit": u, + "stage": stage, + "status": "not_run" if m is None else m.get("status", "failed"), + "shortfalls": shortfalls(m) if m else {}, + "reasons": reasons(m) if m else [], + "n_bad_stages": len(bad), + }) + return rows + + +def print_table(title, rows, limit=25): + print(f"\n{title} ({len(rows)} affected)") + if not rows: + print(" — none") + return + print(f" {'unit':<14} {'stage':<20} {'status':<8} why") + for r in rows[:limit]: + short = ", ".join(f"{k} {v[0]}/{v[1]}" for k, v in r["shortfalls"].items()) + why = (r["reasons"][0] if r["reasons"] else short) or "-" + print(f" {r['unit']:<14} {r['stage']:<20} {r['status']:<8} {why[:90]}") + if len(rows) > limit: + print(f" … and {len(rows) - limit} more (see the JSON report)") + + +def print_stage_table(title, per_stage, n_units): + print(f"\n{title} ({n_units} units declared)") + print(f" {'stage':<20} {'ok':>6} {'warn':>6} {'fail':>6} {'not run':>8} " + f"{'cleaned':>8} attrition") + for stage, t in per_stage.items(): + att = [f"{r} {a['found']}/{a['expect']}" + for r, a in t["products"].items() if a["found"] < a["expect"]] + print(f" {stage:<20} {t['complete']:>6} {len(t['warn']):>6} " + f"{len(t['failed']):>6} {len(t['not_run']):>8} " + f"{len(t.get('cleaned', [])):>8} {', '.join(att)[:60]}") + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--run-dir", required=True, type=Path) + p.add_argument("--index", required=True, type=Path) + p.add_argument("--status", default="manual") + p.add_argument("--out", type=Path, default=None) + p.add_argument("--limit", type=int, default=25, + help="rows per stdout table (the JSON report is complete)") + args = p.parse_args() + + tiles, exps, tile_exp = [], [], defaultdict(list) + if args.index.exists(): + con = sqlite3.connect(args.index, timeout=60) + tiles = [r[0] for r in con.execute("SELECT tile_id FROM tiles ORDER BY 1")] + exps = [r[0] for r in con.execute("SELECT exp_id FROM exposures ORDER BY 1")] + for tile_id, exp_id in con.execute("SELECT tile_id, exp_id FROM tile_exposures"): + tile_exp[tile_id].append(exp_id) + con.close() + else: + print(f"[run_report] no index at {args.index} — reporting manifests only", + file=sys.stderr) + + tile_m = load_manifests(args.run_dir, "tiles") + exp_m = load_manifests(args.run_dir, "exp") + # Reclaimed exposures speak through their tombstones (D5, S5). + cleaned_exp = absorb_tombstones(args.run_dir, "exp", exp_m) + tiles = tiles or sorted(tile_m) + exps = exps or sorted(exp_m) + + missing_json = args.index.parent / "missing.json" + missing = json.loads(missing_json.read_text()) if missing_json.exists() else [] + + report = { + "status": args.status, + "n_tiles": len(tiles), "n_exposures": len(exps), + "missing_tiles": missing, + "tile_stages": tally_level(tiles, TILE_STAGES, tile_m), + "exp_stages": tally_level(exps, EXP_STAGES, exp_m, cleaned_exp), + "cleaned_exposures": sorted(cleaned_exp), + "tiles": unit_rows(tiles, TILE_STAGES, tile_m), + "exposures": unit_rows(exps, EXP_STAGES, exp_m), + } + + # Blame propagation: a BLOCKING exposure blocks every tile that reads it. + # Without this, a tile stalled at tile_vignets looks like its own failure. + # + # Blocking means "failed" or "never ran" — NOT "warn". Warn is the expected + # per-CCD attrition (setools rejecting a sparse CCD, psfex_interp short an + # epoch); it is present in essentially every exposure at production scale, so + # counting it here made every exposure block every tile and the table said + # nothing. + # Judged over ALL the exposure's stages, not just the first bad one, so an + # exposure that warns early and fails late still blocks. + # A CLEANED exposure never blocks: its store is gone precisely because every + # consuming tile already had its vignets. Its absorbed manifests are read + # above, so a cleaned exposure that genuinely failed still shows in the + # tables — it just does not get to hold complete tiles hostage. + def _blocks(unit): + if unit in cleaned_exp: + return False + for stage in EXP_STAGES: + m = exp_m.get(unit, {}).get(stage) + if m is None or m.get("status", "failed") == "failed": + return True + return False + + bad_exp = {e for e in exps if _blocks(e)} + blocked = {t: sorted(set(tile_exp.get(t, [])) & bad_exp) for t in tiles} + report["tiles_blocked_by_exposures"] = {t: e for t, e in blocked.items() if e} + + done = report["tile_stages"]["tile_make_cat"]["complete"] + report["final_cats"] = {"present": done, "of": len(tiles)} + + out = args.out or (args.index.parent / "run_report.json") + out.parent.mkdir(parents=True, exist_ok=True) + out.write_text(json.dumps(report, indent=2, sort_keys=True) + "\n") + + print(f"[run_report] status={args.status} {done}/{len(tiles)} final cats" + + (f" ({len(missing)} tiles missing exposure lists)" if missing else "") + + (f" ({len(cleaned_exp)} exposures reclaimed)" if cleaned_exp else "")) + print_stage_table("EXPOSURES", report["exp_stages"], len(exps)) + print_stage_table("TILES", report["tile_stages"], len(tiles)) + print_table("exposures not complete", report["exposures"], args.limit) + print_table("tiles not complete", report["tiles"], args.limit) + nb = report["tiles_blocked_by_exposures"] + if nb: + print(f"\ntiles waiting on incomplete exposures ({len(nb)})") + for t, e in list(nb.items())[:args.limit]: + print(f" {t:<14} {', '.join(e[:6])}" + + (f" (+{len(e) - 6})" if len(e) > 6 else "")) + print(f"\n[run_report] -> {out}") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/star_cats.py b/workflow/scripts/star_cats.py new file mode 100644 index 000000000..56886cfa0 --- /dev/null +++ b/workflow/scripts/star_cats.py @@ -0,0 +1,354 @@ +#!/usr/bin/env python3 +"""The campaign's GSC 2.3 star catalogue, as a HEALPix-chunked sky store. + +Masking needs, for every exposure, the bright stars over its focal plane. The +sky does not change between exposures, so the network cost of that is a property +of the campaign's SKY AREA, not of its exposure count: exposures overlap each +other ~7-10 deep, and a tile's exposures all look at the same square degree. + +So the store is chunked by sky, not by exposure. One GSC 2.3 cone query per +HEALPix pixel of NSIDE=32, written run-independently under the ``star_cats`` +config root and never fetched twice. A campaign that grows past the fetched +footprint queries only the chunks its new tiles add; one that grows within it +queries nothing. + +Two numbers set the scale. A full-UNIONS footprint is ~1.5k chunks against ~25k +exposures, so the QUERY COUNT drops ~16x. The queried AREA drops ~4x: the old +design covered the footprint ~8 times over (that is just the exposure overlap +depth), the new one ~2 times, the 2x being the price of bounding a HEALPix +quadrilateral by the cone Vizier speaks (see ``pixel_cone``) — 5.6-9 deg^2 for a +3.36 deg^2 pixel, ~40-60k rows and ~3-4 MB per chunk. + +Two subcommands, one module, deliberately: ``fetch`` and ``cut`` must agree +EXACTLY on which pixel holds which star, and a shared NSIDE constant in one file +is the only version of that agreement which cannot drift. + + fetch --tile-list ... --store ... --manifest ... + The campaign side. Turns the tile list into the set of pixels its + exposures can possibly need, fetches the missing ones, writes a manifest. + + cut --images ... --store ... --out ... + The per-exposure side, purely local: read the focal-plane footprint from + the exposure's image headers, load the chunks covering it, deduplicate, + and cut to the focal-plane disc. Reproduces byte-for-byte the same sky + selection the old one-query-per-exposure cone did. + +Geometry, and why the fetch pad is what it is. Chunk-need is computed from the +TILE list rather than from exposure pointings, because tile IDs are the one thing +known before any download: a pointing center means reading a FITS header of an +image get_images has not fetched yet, and the DAG needs the chunk set at parse +time. Tiles sit on a fixed 0.5 deg grid (``cfis.get_tile_coord_from_nixy``), so +each tile is a disc of half-diagonal 0.354 deg; find_exposures gives a tile every +exposure whose footprint covers it, and the MegaCam focal plane is a disc of +radius 0.73 deg (measured on the cached catalogues). An exposure center is +therefore at most 0.354 + 0.73 deg from the tile center, and its stars 0.73 deg +beyond that: 1.81 deg, padded to ``PAD_DEG`` = 2.0. The pad is a perimeter cost +— negligible for a contiguous campaign, and paid once. + +The pad is a bound, not a promise: ``cut`` verifies that every chunk covering the +exposure it was handed is on disk, and fails loudly if one is not. A missing +chunk means the geometry above is wrong, and that must not degrade quietly into +an under-masked exposure. +""" + +import argparse +import json +import os +import sys +from concurrent.futures import ThreadPoolExecutor +from pathlib import Path + +import numpy as np +import healpy as hp +from astropy import units as u +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.table import Table, vstack +from astropy.wcs import WCS + +# GSC 2.3. The same catalogue the mask module's own CDS path uses +# (mask.py: _CDS_cat_ID), so the store is a drop-in for it. +CAT_ID = "I/305/out" + +# NSIDE=32 -> 3.36 deg^2 per pixel, 12288 pixels over the sky. Chosen so one +# chunk is a couple of MegaCam focal planes: small enough that a Vizier query +# stays within a small multiple of the per-exposure queries this replaces, large +# enough that a full-UNIONS footprint is ~1.5k chunks rather than ~25k. The +# cone-vs-quadrilateral overhead is scale-free, so NSIDE trades query count +# against query size and nothing else. NESTED, so a chunk id is a hierarchical +# sky address and a future NSIDE change is a subdivision. +NSIDE = 32 +NEST = True + +# Angular padding on the disc used to select chunks (see the module docstring). +PAD_DEG = 2.0 + +# The MegaCam focal-plane disc, and the margin added to a pixel's own bounding +# cone. Both in degrees. +MARGIN_DEG = 0.02 + +# GSC 2.3's object id: the deduplication key where chunk cones overlap. +ID_COL = "GSC2.3" + + +# --- the chunk store -------------------------------------------------------- + + +def store_dir(store: Path) -> Path: + """Chunks live under the catalogue and resolution that produced them, so a + later NSIDE or catalogue change is a new directory beside the old one rather + than a silent reinterpretation of files already on disk.""" + return Path(store) / CAT_ID.replace("/", "_") / f"nside{NSIDE}" + + +def chunk_path(store: Path, ipix: int) -> Path: + return store_dir(store) / f"star_chunk-{ipix:06d}.fits" + + +def chunks_for_disc(ra_deg: float, dec_deg: float, radius_deg: float) -> list[int]: + """Every pixel that touches the disc, as sorted ids. + + ``inclusive=True`` makes this a conservative superset — the guarantee ``cut`` + relies on is that no star inside the disc lives in a pixel this omits. + """ + vec = hp.ang2vec(ra_deg, dec_deg, lonlat=True) + return sorted(int(i) for i in hp.query_disc( + NSIDE, vec, np.radians(radius_deg), inclusive=True, fact=4, nest=NEST)) + + +def pixel_cone(ipix: int) -> tuple[float, float, float]: + """(ra, dec, radius_arcmin) of a cone that CONTAINS pixel ``ipix``. + + Vizier speaks cones, HEALPix speaks quadrilaterals, so the query is the + pixel's bounding cone: its center, and the largest center-to-boundary + distance plus a margin. The cone spills over the pixel edges, which costs a + little duplication between neighbours and buys the containment ``cut`` + depends on. The duplicates are removed on read, by ``ID_COL``. + """ + ra_c, dec_c = hp.pix2ang(NSIDE, ipix, nest=NEST, lonlat=True) + ra_b, dec_b = hp.vec2ang(hp.boundaries(NSIDE, ipix, step=8, nest=NEST).T, + lonlat=True) + center = SkyCoord(ra_c * u.deg, dec_c * u.deg) + radius = center.separation(SkyCoord(ra_b * u.deg, dec_b * u.deg)).deg.max() + return float(ra_c), float(dec_c), float((radius + MARGIN_DEG) * 60.0) + + +def write_atomic(table: Table, path: Path) -> None: + """Publish ``table`` at ``path`` all-or-nothing. + + The store's only cache test is ``Path.exists``, so a write killed part-way + (timeout, OOM, node failure) would otherwise leave a truncated FITS that + every later run trusts forever. The temp keeps the ``.fits`` suffix because + astropy picks its writer from the extension, and is dot-prefixed and + PID-tagged so it stays out of the ``star_chunk-*`` globs and two concurrent + writers cannot collide. + """ + tmp = path.parent / f".tmp-{os.getpid()}-{path.name}" + try: + table.write(tmp, overwrite=True) + os.replace(tmp, path) + finally: + if tmp.exists(): + tmp.unlink() + + +def read_chunks(store: Path, ipixels: list[int]) -> Table: + """Load and deduplicate the given chunks. + + A missing chunk is fatal (see the module docstring): it means the fetch + footprint did not cover this exposure, and an under-masked exposure is worse + than a failed job. + """ + missing = [i for i in ipixels if not chunk_path(store, i).exists()] + if missing: + raise SystemExit( + f"star chunk(s) {missing} not in {store_dir(store)}. The campaign's " + f"star_catalogue fetch did not cover this exposure — re-run it " + f"(and check that its tile list contains this exposure's tiles).") + + table = vstack([Table.read(chunk_path(store, i)) for i in ipixels], + metadata_conflicts="silent") + _, keep = np.unique(np.asarray(table[ID_COL]), return_index=True) + return table[np.sort(keep)] + + +# --- exposure footprint ----------------------------------------------------- + + +def _wcs(header) -> WCS: + """Build the WCS by hand, from the linear terms only. + + Same construction as ``scripts/python/create_star_cat.py``: it sidesteps + distortion-convention incompatibilities between headers and astropy, and a + focal-plane footprint needs nothing finer. + """ + w = WCS(naxis=2) + w.wcs.ctype = [header["CTYPE1"], header["CTYPE2"]] + try: + w.wcs.cunit = [header["CUNIT1"], header["CUNIT2"]] + except KeyError: + w.wcs.cunit = ["deg", "deg"] + w.wcs.crpix = [header["CRPIX1"], header["CRPIX2"]] + w.wcs.crval = [header["CRVAL1"], header["CRVAL2"]] + w.wcs.cd = [[header["CD1_1"], header["CD1_2"]], + [header["CD2_1"], header["CD2_2"]]] + return w + + +def focal_plane_disc(image: Path, n_ccd: int = 40) -> tuple[float, float, float]: + """(ra, dec, radius_deg) of the disc covering all CCDs of one exposure.""" + centers, radii = [], [] + for ext in range(1, n_ccd + 1): + h = fits.getheader(image, ext) + w = _wcs(h) + (ra_c, dec_c), (ra_0, dec_0) = w.all_pix2world( + [[h["NAXIS1"] / 2.0, h["NAXIS2"] / 2.0], [0, 0]], 1) + centers.append(SkyCoord(ra_c * u.deg, dec_c * u.deg)) + radii.append(centers[-1].separation( + SkyCoord(ra_0 * u.deg, dec_0 * u.deg)).deg) + + ras = np.array([c.ra.deg for c in centers]) + decs = np.array([c.dec.deg for c in centers]) + center = SkyCoord(ras.mean() * u.deg, decs.mean() * u.deg) + seps = center.separation(SkyCoord(ras * u.deg, decs * u.deg)).deg + return (float(center.ra.deg), float(center.dec.deg), + float(np.max(seps + np.array(radii)))) + + +def exposure_image(images_dir: Path) -> Path: + """The one multi-extension exposure image in a get_images output dir. + + That dir is a symlink farm holding ``image-.fitsfz`` plus its weight and + flag; only the image carries the 40 CCD WCSs. + """ + found = sorted(p for p in Path(images_dir).iterdir() if "image" in p.name) + if not found: + raise SystemExit(f"no image file in {images_dir}") + return found[0] + + +# --- fetch ------------------------------------------------------------------ + + +def campaign_chunks(tile_ids: list[str]) -> list[int]: + """Every chunk the campaign's exposures can need, from the tile list alone.""" + from shapepipe.utilities.cfis import get_tile_coord_from_nixy + + needed: set[int] = set() + for tile_id in tile_ids: + nix, niy = tile_id.split(".") + ra, dec = get_tile_coord_from_nixy(nix, niy) + needed.update(chunks_for_disc(ra.degree, dec.degree, PAD_DEG)) + return sorted(needed) + + +def fetch(args: argparse.Namespace) -> None: + from shapepipe.utilities.vizier import query_vizier + + tile_ids = [ln.strip() for ln in Path(args.tile_list).read_text().splitlines() + if ln.strip()] + needed = campaign_chunks(tile_ids) + out_dir = store_dir(args.store) + out_dir.mkdir(parents=True, exist_ok=True) + todo = [i for i in needed if not chunk_path(args.store, i).exists()] + print(f"star chunks: {len(needed)} needed for {len(tile_ids)} tiles, " + f"{len(todo)} to fetch -> {out_dir}", file=sys.stderr) + + def one(ipix: int) -> int: + ra, dec, radius_arcmin = pixel_cone(ipix) + table = query_vizier(ra, dec, radius_arcmin, CAT_ID) + write_atomic(table, chunk_path(args.store, ipix)) + print(f"chunk {ipix}: {len(table)} rows " + f"(ra={ra:.4f} dec={dec:.4f} r={radius_arcmin:.1f}')", + file=sys.stderr) + return len(table) + + # A handful of concurrent queries, never one per exposure: the same modest + # concurrency the per-exposure rule reached through --local-cores, now an + # explicit number instead of an accident of the head node's CPU count. + if todo: + with ThreadPoolExecutor(max_workers=args.workers) as pool: + list(pool.map(one, todo)) + + write_manifest(args, tile_ids, needed, len(todo)) + + +def write_manifest(args, tile_ids, needed, n_fetched) -> None: + """The rule's declared output. + + Not a ``completeness.py`` verdict: this rule runs no ``shapepipe_run`` and + has no per-runner count floors, so there is nothing to compose and no + separate ``log:`` — under ``set -euo pipefail`` the job either completes or + aborts at the failing query, and snakemake's captured stderr is the evidence. + The manifest keeps the workflow's "one rule, one manifest" currency: written + last, and only when the content changed, so an unchanged campaign leaves the + mtime where it was rather than churning the `mtime` rerun-trigger. + """ + body = json.dumps({ + "stage": "star_catalogue", + "level": "campaign", + "status": "complete", + "catalogue": CAT_ID, + "nside": NSIDE, + "nest": NEST, + "pad_deg": PAD_DEG, + "store": str(store_dir(args.store)), + "n_tiles": len(tile_ids), + "n_chunks": len(needed), + "n_fetched": n_fetched, + "chunks": needed, + }, indent=2, sort_keys=True) + path = Path(args.manifest) + path.parent.mkdir(parents=True, exist_ok=True) + if not path.exists() or path.read_text() != body: + path.write_text(body) + + +# --- cut -------------------------------------------------------------------- + + +def cut(args: argparse.Namespace) -> None: + image = exposure_image(args.images) + ra, dec, radius = focal_plane_disc(image) + ipixels = chunks_for_disc(ra, dec, radius) + table = read_chunks(args.store, ipixels) + + center = SkyCoord(ra * u.deg, dec * u.deg) + stars = SkyCoord(np.asarray(table["RAJ2000"]) * u.deg, + np.asarray(table["DEJ2000"]) * u.deg) + inside = table[center.separation(stars).deg <= radius] + + print(f"{image.name}: ra={ra:.4f} dec={dec:.4f} r={radius:.4f} deg, " + f"{len(ipixels)} chunks -> {len(inside)} stars", file=sys.stderr) + out = Path(args.out) + out.parent.mkdir(parents=True, exist_ok=True) + write_atomic(inside, out) + + +# --- CLI -------------------------------------------------------------------- + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + sub = p.add_subparsers(dest="cmd", required=True) + + f = sub.add_parser("fetch", help="fetch the campaign footprint's chunks") + f.add_argument("--tile-list", required=True, type=Path) + f.add_argument("--store", required=True, type=Path) + f.add_argument("--manifest", required=True, type=Path) + f.add_argument("--workers", type=int, default=4) + f.set_defaults(func=fetch) + + c = sub.add_parser("cut", help="cut one exposure's catalogue from the store") + c.add_argument("--images", required=True, type=Path, + help="a get_images output dir holding image-.fitsfz") + c.add_argument("--store", required=True, type=Path) + c.add_argument("--out", required=True, type=Path) + c.set_defaults(func=cut) + + args = p.parse_args() + args.func(args) + + +if __name__ == "__main__": + main()