Snakemake orchestration for the real-data pipeline

.gitignore

@@ -139,3 +139,4 @@ code
 # .felt here is a machine-local symlink into it. Never track it in this repo.
 /.felt/
 /.felt
+.snakemake/

profiles/nibi/config.yaml

@@ -0,0 +1,87 @@
+# 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; the container is set in the Snakefile (local .sif).
+
+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 <workdir>, 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 <cwd>; restores cadcproxy.pem for vos/vcp
+#  PYTHONPATH          SETTLED CALL 3 — P0 pins THIS branch's src/ (identical to
+#                      develop@97e16d50: orchestration commits never touch src/).
+#                      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). 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
+
+# D3 DEPENDS ON THIS. Snakemake's default is to DELETE the declared outputs of a
+# failed job — and this workflow's declared output IS the manifest, the only
+# record of *why* a unit failed and the only thing `sp report` reads. Without
+# keep-incomplete, every failed unit reads back as "not_run" and the post-mortem
+# debris is gone. Safe here because each rule `rm -rf`s its own run dir at start,
+# so a rerun never sees stale products.
+keep-incomplete: true
+
+# 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 (uncompress, merges) into their chunky
+# neighbours (PRD D-profile) is NOT set here — it requires labels in the rule
+# files themselves, which is out of scope for this profile-only pass. Deferred
+# to whichever slice next touches workflow/rules/*.smk.

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")

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"]

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 = (

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,
     )
 

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")