diff --git a/scripts/python/create_final_cat.py b/scripts/python/create_final_cat.py index 2b583b857..ea8a4d723 100755 --- a/scripts/python/create_final_cat.py +++ b/scripts/python/create_final_cat.py @@ -144,13 +144,17 @@ def read_param_file(path, verbose=False): print("No parameters read", end="") print(" into merged catalogue") - param_list_unique = list(set(param_list)) - + # Ordered dedup. list(set(...)) reordered the columns by the process's + # string hash seed, so two runs of this tool over the same inputs produced + # files whose datasets differed in column ORDER — which is part of a + # structured dtype, and therefore part of the file. + param_list_unique = list(dict.fromkeys(param_list)) + if verbose: n = len(param_list) - len(param_list_unique) - if n > 1: - print("Removed {n} duplicate entries") - + if n > 0: + print(f"Removed {n} duplicate entries") + return param_list_unique @@ -312,16 +316,20 @@ def read_data(fits_file, params): if params["param_list"] is None: params["param_list"] = [col for col in data.keys()] - try: - extracted_data = {col: data[col] for col in params["param_list"]} - dtype = data.dtype - except: - print(f"Error for ID {id}, path {fits_file}") - for col in params["param_list"]: - if col not in data: - print(col, end=" ") - print() - continue + # RAISE, do not print and fall through. The bare `except:` this replaces + # left extracted_data and dtype unbound, so the caller's own error was an + # UnboundLocalError from the return statement below, naming neither the + # file nor the column that was actually missing. + present = set(data.dtype.names or ()) + missing = [col for col in params["param_list"] if col not in present] + if missing: + raise KeyError( + f"{fits_file}: missing {len(missing)} of the " + f"{len(params['param_list'])} requested column(s): " + f"{' '.join(missing)}" + ) + extracted_data = {col: data[col] for col in params["param_list"]} + dtype = data.dtype return extracted_data, dtype @@ -330,16 +338,29 @@ def copy_data(param_list, extracted_data, dtype): """Copy Data. """ + # THE REQUESTED COLUMNS ONLY, IN THE PARAMETER FILE'S ORDER. Two things + # are being fixed here and they are easy to conflate. Allocating with the + # source's full dtype and filling only the requested columns left every + # other column as uninitialised memory — meaningless values, and different + # bytes on every run over the same inputs. And ordering the result by the + # SOURCE catalogue's columns made the output dtype a property of the + # catalogue rather than of the parameter file: two tiles written by + # different ShapePipe versions, whose catalogues order or extend their + # columns differently, then landed in one merged file with two different + # structured dtypes, which np.concatenate refuses. The parameter file is + # the schema; it says which columns AND in what order. + wanted = set(dtype.names or ()) + columns = [col for col in param_list if col in wanted] + subset = np.dtype([(col, dtype[col]) for col in columns]) + # Initialize new data structure structured_data = np.empty( len(extracted_data[param_list[0]]), - dtype=dtype, + dtype=subset, ) # Loop over parameters - for col in param_list: - if not col in extracted_data: - print(f"Column {col} not in file with ID {id}") + for col in columns: structured_data[col] = extracted_data[col] #if isinstance(extracted_data[col][0], (np.ndarray, tuple, list)): @@ -467,12 +488,14 @@ def process(params): structured_data = copy_data(params["param_list"], extracted_data, dtype) - # Create a new dataset + # Create a new dataset. dtype comes from the array copy_data + # built, not from the source catalogue: they differ now that + # copy_data allocates the requested columns alone. try: patch_group.create_dataset( str(id), data=structured_data, - dtype=dtype, + dtype=structured_data.dtype, ) except: print(f"Error for {id}: Could not create dataset in group {patch}") diff --git a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py index 7bc0cb76b..1e0cc7256 100644 --- a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py +++ b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py @@ -17,6 +17,36 @@ from shapepipe.pipeline import file_io +def _stack(chunks, dtype=None): + """Concatenate one column's per-catalogue arrays into a single array. + + THE COLUMN ACCUMULATORS ARE LISTS OF ARRAYS, ONE PER INPUT CATALOGUE, and + not lists of values, because these classes are the last step of a whole + campaign. ``x += list(data["X"])`` turns 4 bytes of float32 payload into a + 32-byte python object plus an 8-byte pointer in a list that overallocates — + measured at ~10x the input bytes end to end, which put a full-survey merge + (~20k exposures x 40 CCDs) at ~400 GB of RAM and made it unrunnable on any + node. One array per catalogue plus one concatenate at the end holds ~1x, and + produces the identical output: np.array() over a list of numpy scalars and + np.concatenate() over the arrays they came from agree on dtype and on order. + + IT EMPTIES THE LIST IT IS GIVEN, and that is not a side effect to tidy away + later — it is half the saving. np.concatenate holds the chunks and the + result at once, so a caller that stacks sixteen columns while all sixteen + chunk lists are still alive peaks at twice the campaign. Released column by + column, the peak is one campaign plus one column. Callers stack once, at the + end, and do not touch the accumulators afterwards. + + An empty input list is a merge over no catalogues, which the callers guard + against; it returns an empty array so the output column still exists. + """ + if not chunks: + return np.array([], dtype=dtype or np.float64) + out = np.concatenate(chunks) + del chunks[:] + return out + + class MergeStarCatMCCD(object): """Merge Star Catalogue MCCD. @@ -315,66 +345,37 @@ def process(self): model_var.append(model_var_val) model_var_size.append(model_var_val.size) + # ONE ARRAY PER CATALOGUE PER COLUMN (see _stack): the per-value + # python lists this replaces cost ~10x the input bytes. # positions - x += list( - starcat_j[self._hdu_table].data["GLOB_POSITION_IMG_LIST"][:, 0] - ) - y += list( - starcat_j[self._hdu_table].data["GLOB_POSITION_IMG_LIST"][:, 1] - ) + pos = starcat_j[self._hdu_table].data["GLOB_POSITION_IMG_LIST"] + x.append(np.asarray(pos[:, 0])) + y.append(np.asarray(pos[:, 1])) # RA and DEC positions try: - ra += list(starcat_j[self._hdu_table].data["RA_LIST"][:]) - dec += list(starcat_j[self._hdu_table].data["DEC_LIST"][:]) + ra.append(np.asarray(starcat_j[self._hdu_table].data["RA_LIST"][:])) + dec.append(np.asarray(starcat_j[self._hdu_table].data["DEC_LIST"][:])) except Exception: - ra += list( - np.zeros( - starcat_j[self._hdu_table] - .data["GLOB_POSITION_IMG_LIST"][:, 0] - .shape, - dtype=int, - ) - ) - dec += list( - np.zeros( - starcat_j[self._hdu_table] - .data["GLOB_POSITION_IMG_LIST"][:, 0] - .shape, - dtype=int, - ) - ) + ra.append(np.zeros(pos[:, 0].shape, dtype=int)) + dec.append(np.zeros(pos[:, 0].shape, dtype=int)) # shapes (convert sigmas to T = 2 sigma^2) - g1_psf += list( - starcat_j[self._hdu_table].data["PSF_MOM_LIST"][:, 0] - ) - g2_psf += list( - starcat_j[self._hdu_table].data["PSF_MOM_LIST"][:, 1] - ) - size_psf += list( - cs_size.sigma_to_T( - starcat_j[self._hdu_table].data["PSF_MOM_LIST"][:, 2] - ) - ) - g1 += list(starcat_j[self._hdu_table].data["STAR_MOM_LIST"][:, 0]) - g2 += list(starcat_j[self._hdu_table].data["STAR_MOM_LIST"][:, 1]) - size += list( - cs_size.sigma_to_T( - starcat_j[self._hdu_table].data["STAR_MOM_LIST"][:, 2] - ) - ) + psf_mom = starcat_j[self._hdu_table].data["PSF_MOM_LIST"] + star_mom = starcat_j[self._hdu_table].data["STAR_MOM_LIST"] + g1_psf.append(np.asarray(psf_mom[:, 0])) + g2_psf.append(np.asarray(psf_mom[:, 1])) + size_psf.append(np.asarray(cs_size.sigma_to_T(psf_mom[:, 2]))) + g1.append(np.asarray(star_mom[:, 0])) + g2.append(np.asarray(star_mom[:, 1])) + size.append(np.asarray(cs_size.sigma_to_T(star_mom[:, 2]))) # flags - flag_psf += list( - starcat_j[self._hdu_table].data["PSF_MOM_LIST"][:, 3] - ) - flag_star += list( - starcat_j[self._hdu_table].data["STAR_MOM_LIST"][:, 3] - ) + flag_psf.append(np.asarray(psf_mom[:, 3])) + flag_star.append(np.asarray(star_mom[:, 3])) # ccd id list - ccd_nb += list(starcat_j[self._hdu_table].data["CCD_ID_LIST"]) + ccd_nb.append(np.asarray(starcat_j[self._hdu_table].data["CCD_ID_LIST"])) starcat_j.close() @@ -447,15 +448,21 @@ def process(self): ) # Mask and transform to numpy arrays - flagmask = np.abs(np.array(flag_star) - 1) * np.abs( - np.array(flag_psf) - 1 - ) - psf_e1 = np.array(g1_psf)[flagmask.astype(bool)] - psf_e2 = np.array(g2_psf)[flagmask.astype(bool)] - psf_r2 = np.array(size_psf)[flagmask.astype(bool)] - star_e1 = np.array(g1)[flagmask.astype(bool)] - star_e2 = np.array(g2)[flagmask.astype(bool)] - star_r2 = np.array(size)[flagmask.astype(bool)] + # Concatenate once, here: everything below already wanted arrays and + # was calling np.array() on python lists to get them (see _stack). + x, y, ra, dec = _stack(x), _stack(y), _stack(ra), _stack(dec) + g1_psf, g2_psf, size_psf = _stack(g1_psf), _stack(g2_psf), _stack(size_psf) + g1, g2, size = _stack(g1), _stack(g2), _stack(size) + flag_psf, flag_star = _stack(flag_psf), _stack(flag_star) + ccd_nb = _stack(ccd_nb) + + flagmask = np.abs(flag_star - 1) * np.abs(flag_psf - 1) + psf_e1 = g1_psf[flagmask.astype(bool)] + psf_e2 = g2_psf[flagmask.astype(bool)] + psf_r2 = size_psf[flagmask.astype(bool)] + star_e1 = g1[flagmask.astype(bool)] + star_e2 = g2[flagmask.astype(bool)] + star_r2 = size[flagmask.astype(bool)] rmse, mean, std_dev = MSC.stats_calculator(star_e1, psf_e1) self._w_log.info( @@ -551,72 +558,134 @@ def __init__( self._hdu_table = hdu_table self._input_cat_type = input_cat_type + # The columns this class writes, and where each comes from. Kept as data + # rather than as sixteen repeated lines, because a two-pass merge would + # otherwise state every column three times: to size it, to allocate it and + # to fill it. + _COLUMNS = ( + ("X", "X"), ("Y", "Y"), ("RA", "RA"), ("DEC", "DEC"), + ("HSM_G1_PSF", "HSM_G1_PSF"), ("HSM_G2_PSF", "HSM_G2_PSF"), + ("HSM_T_PSF", "HSM_T_PSF"), ("HSM_G1_STAR", "HSM_G1_STAR"), + ("HSM_G2_STAR", "HSM_G2_STAR"), ("HSM_T_STAR", "HSM_T_STAR"), + ("HSM_FLAG_PSF", "HSM_FLAG_PSF"), ("HSM_FLAG_STAR", "HSM_FLAG_STAR"), + ) + # Present in psfex_interp output, absent from pix2wcs-converted files + # (MKDEBUG); zero-filled when missing rather than failing the merge. + _OPTIONAL = (("MAG", "MAG"), ("SNR", "SNR"), ("ACCEPTED", "ACCEPTED")) + + def _ccd_nb(self, path): + """The CCD number this catalogue's rows carry, parsed from its name.""" + return re.split(r"\-([0-9]*)\-([0-9]+)\.", path)[-2] + def process(self): """Process. Process merging. + TWO PASSES, AND NEITHER HOLDS THE CAMPAIGN TWICE. The first reads only + the FITS HEADER of every input — NAXIS2, the row count — and never + touches a data block; the second allocates the output columns once, at + their exact final length, and fills them slice by slice. Peak memory is + therefore ONE output plus ONE input catalogue. + + What this replaces, in two steps, is instructive about the cost of the + obvious code. Accumulating each column into a python LIST OF VALUES — + ``x += list(data["X"])`` — turned 4 bytes of float32 payload into a + 32-byte object plus an 8-byte pointer, measured at ~10x the input bytes + end to end and putting a full-survey merge (~20k exposures x 40 CCDs) at + ~400 GB. Accumulating one ARRAY PER CATALOGUE and concatenating once + brought that to ~5.5x. This pass structure removes what was left of the + accumulation: there are no chunks, and no concatenate that must hold its + inputs and its result at the same time. + + ``self._input_file_list`` MUST BE ITERABLE TWICE, which the module + runner's list is. A one-shot generator is not, and would silently merge + nothing on the second pass — hence the explicit row-count check below. """ - x, y, ra, dec = [], [], [], [] - g1_psf, g2_psf, size_psf = [], [], [] - g1, g2, size = [], [], [] - flag_psf, flag_star = [], [] - mag, snr, psfex_acc = [], [], [] - ccd_nb = [] - self._w_log.info( f"Merging {len(self._input_file_list)} star catalogues" ) + # --- pass 1: row counts and dtypes, from headers alone -------------- + # THE OPTIONAL COLUMNS ARE A PER-FILE QUESTION, NOT A PER-MERGE ONE. + # A pix2wcs-converted catalogue has no MAG/SNR/ACCEPTED while an + # ordinary one does, and a merge can be handed both. Deciding from the + # first file alone got it wrong in both directions: converted-first + # zero-filled the real values of every ordinary file behind it, and + # ordinary-first raised KeyError on the first converted one. So the + # dtype comes from ANY file that carries the column, and pass 2 asks + # each file for itself. + names, dtypes, opt_dtypes, n_total = [], None, {}, 0 for name in self._input_file_list: try: - starcat_j = fits.open(name[0], memmap=False, ignore_missing_simple=True) - except OSError as e: + with fits.open(name[0], memmap=False, + ignore_missing_simple=True) as starcat_j: + hdu = starcat_j[self._hdu_table] + n_rows = hdu.header["NAXIS2"] + # ColDefs.dtype describes the table without reading it. + # NOTE: it is the RAW storage dtype and ignores TSCAL/TZERO, + # so a scaled column would be allocated narrower than the + # values .data returns. Latent, not live: no validation_psf + # column is scaled. Read the dtype off .data if one ever is. + cols = hdu.columns.dtype + if dtypes is None: + dtypes = cols + for _, col in self._OPTIONAL: + if col not in opt_dtypes and col in (cols.names or ()): + opt_dtypes[col] = cols[col] + except OSError: print(f"Error while opening file '{name[0]}'") #raise continue - - data_j = starcat_j[self._hdu_table].data - - # positions - x += list(data_j["X"]) - y += list(data_j["Y"]) - ra += list(data_j["RA"]) - dec += list(data_j["DEC"]) - - # shapes (size column already holds T = 2 sigma^2) - g1_psf += list(data_j["HSM_G1_PSF"]) - g2_psf += list(data_j["HSM_G2_PSF"]) - size_psf += list(data_j["HSM_T_PSF"]) - g1 += list(data_j["HSM_G1_STAR"]) - g2 += list(data_j["HSM_G2_STAR"]) - size += list(data_j["HSM_T_STAR"]) - - # flags - flag_psf += list(data_j["HSM_FLAG_PSF"]) - flag_star += list(data_j["HSM_FLAG_STAR"]) - - # misc - - # MKDEBUG: The following columns do not exist (yet) - # for psf converted (pix2wcs) files. - try: - mag += list(data_j["MAG"]) - except: - mag += list(np.zeros_like(data_j["X"])) - try: - snr += list(data_j["SNR"]) - except: - snr += list(np.zeros_like(data_j["X"])) + names.append(name[0]) + n_total += n_rows + + if dtypes is None: + raise ValueError("merge_starcat: no readable input catalogue") + + # --- allocate once, at the exact final length ----------------------- + data = {out: np.empty(n_total, dtype=dtypes[col]) + for out, col in self._COLUMNS} + for out, col in self._OPTIONAL: + # A column no file carries still gets a column, zero-filled, in the + # positional dtype the old code used for it. + data[out] = np.empty(n_total, dtype=opt_dtypes.get(col, dtypes["X"])) + # CCD_NB is one string per catalogue, repeated over its rows; its width + # is the widest CCD number in the campaign, which pass 1 already knows. + width = max((len(self._ccd_nb(n)) for n in names), default=1) + data["CCD_NB"] = np.empty(n_total, dtype=f"U{width}") + + # --- pass 2: fill --------------------------------------------------- + at = 0 + for name in self._input_file_list: try: - psfex_acc += list(data_j["ACCEPTED"]) - except: - psfex_acc += list(np.zeros_like(data_j["X"])) + starcat_j = fits.open(name[0], memmap=False, + ignore_missing_simple=True) + except OSError: + continue + data_j = starcat_j[self._hdu_table].data + n_rows = len(data_j) + sl = slice(at, at + n_rows) + + have = set(data_j.dtype.names or ()) + for out, col in self._COLUMNS: + data[out][sl] = data_j[col] + for out, col in self._OPTIONAL: + # THIS file's schema, not the merge's: zero-fill only the files + # that actually lack the column. + data[out][sl] = data_j[col] if col in have else 0 + data["CCD_NB"][sl] = self._ccd_nb(name[0]) + + at += n_rows + starcat_j.close() - # CCD number - ccd_nb += [re.split(r"\-([0-9]*)\-([0-9]+)\.", name[0])[-2]] * len( - data_j["RA"] - ) + if at != n_total: + # The two passes disagreed: an input changed under us, or the list + # was a one-shot iterable. Either way the output would be padded + # with uninitialised memory, so say so rather than write it. + raise ValueError( + f"merge_starcat: pass 1 counted {n_total} rows, pass 2 filled " + f"{at} — is the input list iterable more than once?") # Prepare output FITS catalogue # MKDEBUG: SEx_cat=True -> False @@ -627,25 +696,8 @@ def process(self): SEx_catalogue=False, ) - # Collect columns (size stored as T = 2 sigma^2) - data = { - "X": x, - "Y": y, - "RA": ra, - "DEC": dec, - "HSM_G1_PSF": g1_psf, - "HSM_G2_PSF": g2_psf, - "HSM_T_PSF": size_psf, - "HSM_G1_STAR": g1, - "HSM_G2_STAR": g2, - "HSM_T_STAR": size, - "HSM_FLAG_PSF": flag_psf, - "HSM_FLAG_STAR": flag_star, - "MAG": mag, - "SNR": snr, - "ACCEPTED": psfex_acc, - "CCD_NB": ccd_nb, - } + # `data` was built by the two passes above (size stored as T = 2 + # sigma^2); every column is already an array of its final length. # Write file # MKDEBUG for psf conv (pix2WCS) files do not write as SExtractorCat; @@ -791,29 +843,36 @@ def process(self): data_j = starcat_j[self._hdu_table].data # positions - x += list(data_j["XWIN_IMAGE"]) - y += list(data_j["YWIN_IMAGE"]) - ra += list(data_j["XWIN_WORLD"]) - dec += list(data_j["YWIN_WORLD"]) - + x.append(np.asarray(data_j["XWIN_IMAGE"])) + y.append(np.asarray(data_j["YWIN_IMAGE"])) + ra.append(np.asarray(data_j["XWIN_WORLD"])) + dec.append(np.asarray(data_j["YWIN_WORLD"])) + + # PRE-EXISTING BUG, LEFT ALONE DELIBERATELY: these four REBIND the + # accumulators initialised above rather than appending to them, so + # only the LAST input file's ellipticities reach the output while + # every other column carries the whole merge. Setools is not wired + # to any workflow path today; fixing it is its own change with its + # own verification, and doing it silently inside a memory rewrite + # would bury it. m11, m20, m02 = self.get_moments(data_j) eps1, eps2 = self.get_ellipticity(m11, m20, m02, "epsilon") chi1, chi2 = self.get_ellipticity(m11, m20, m02, "chi") - size += list(data_j["FLUX_RADIUS"]) + size.append(np.asarray(data_j["FLUX_RADIUS"])) # flags - flags += list(data_j["FLAGS_WIN"]) - flags_ext += list(data_j["IMAFLAGS_ISO"]) + flags.append(np.asarray(data_j["FLAGS_WIN"])) + flags_ext.append(np.asarray(data_j["IMAFLAGS_ISO"])) # misc - mag += list(data_j["MAG_WIN"]) - snr += list(data_j["SNR_WIN"]) + mag.append(np.asarray(data_j["MAG_WIN"])) + snr.append(np.asarray(data_j["SNR_WIN"])) # CCD number - ccd_nb += [re.split(r"\-([0-9]*)\-([0-9]+)\.", name[0])[-2]] * len( - data_j["XWIN_IMAGE"] - ) + ccd_nb.append(np.full( + len(data_j["XWIN_IMAGE"]), + re.split(r"\-([0-9]*)\-([0-9]+)\.", name[0])[-2])) # Prepare output FITS catalogue output = file_io.FITSCatalogue( @@ -825,20 +884,20 @@ def process(self): # Collect columns # convert back to sigma for consistency data = { - "X": x, - "Y": y, - "RA": ra, - "DEC": dec, + "X": _stack(x), + "Y": _stack(y), + "RA": _stack(ra), + "DEC": _stack(dec), "EPS1": eps1, "EPS2": eps2, "CHI1": chi1, "CHI2": chi2, - "SIZE": size, - "FLAGS": flags, - "FLAGS_EXT": flags_ext, - "MAG": mag, - "SNR": snr, - "CCD_NB": ccd_nb, + "SIZE": _stack(size), + "FLAGS": _stack(flags), + "FLAGS_EXT": _stack(flags_ext), + "MAG": _stack(mag), + "SNR": _stack(snr), + "CCD_NB": _stack(ccd_nb, dtype="U1"), } # Write file diff --git a/tests/unit/test_hdf5_reconcile_props.py b/tests/unit/test_hdf5_reconcile_props.py new file mode 100644 index 000000000..a990c3f89 --- /dev/null +++ b/tests/unit/test_hdf5_reconcile_props.py @@ -0,0 +1,336 @@ +"""Property-based state machine over ``workflow/scripts/hdf5_reconcile.py``. + +The module's contract is that an hdf5 catalogue reconciled against a campaign +is a FUNCTION OF ITS INPUT SET — the same units with the same sources give the +same datasets, the same dtypes and the same count attribute, however they got +there. That is a claim about every reachable sequence of appends, refreshes and +removals, not about the three the unit tests happen to walk, so it is tested +here against a model: a random sequence of campaign edits, each followed by a +real plan/apply against a real file on disk, with the model asserted after +every step. + +The operations are the four things a campaign can do between invocations — +add a unit, change a unit's source, drop a unit, change the column set — plus +a no-op, which is the one that must leave the file's mtime alone. + +Source mtimes are set EXPLICITLY with ``os.utime`` rather than left to the +clock. ``stamp()`` is (size, mtime_ns), so a test that rewrote a file with the +same length inside one filesystem tick would silently exercise "nothing +changed" while believing it exercised a refresh. +""" + +import importlib.util +import os +import sys +from pathlib import Path + +import numpy as np +import pytest +from hypothesis import HealthCheck, settings +from hypothesis import strategies as st +from hypothesis.stateful import ( + RuleBasedStateMachine, + initialize, + invariant, + precondition, + rule, +) + +h5py = pytest.importorskip("h5py") + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPTS = REPO_ROOT / "workflow" / "scripts" + + +def _load(name): + path = SCRIPTS / f"{name}.py" + assert path.exists(), f"{path} not found; the rules call it by path" + sys.path.insert(0, str(SCRIPTS)) + try: + spec = importlib.util.spec_from_file_location(f"_{name}", path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + finally: + sys.path.remove(str(SCRIPTS)) + return module + + +reconcile = _load("hdf5_reconcile") + +GROUP = "cat/campaign_a" +COUNT_ATTR = "n_units" +UNITS = ["u0", "u1", "u2", "u3"] +# Two column sets, so a schema change is a real change of dtype and width. +COLUMN_SETS = [("RA", "DEC", "E1"), ("RA", "DEC", "E1", "FWHM")] +# What a rebuild may cost over a from-scratch build of the same campaign: the +# metadata h5py's group copy writes for a moved dataset. Measured at ~1.4 kB +# and constant in the number of rebuilds; the allowance is generous because the +# property being defended is "does not grow with history", not an exact size. +COPY_SLACK = 8192 + + +def _array(columns, rows, seed): + rng = np.random.default_rng(seed) + dtype = [(c, " None: + np.save(path, array, allow_pickle=False) + os.utime(path, ns=(mtime_ns, mtime_ns)) + + +def _read(unit, source): + return np.load(source, allow_pickle=False) + + +def _build(output: Path, sources: dict, columns): + """Plan and apply once, exactly as the two merge rules do; return the plan.""" + units = sorted(sources.items()) + digest = reconcile.schema_digest(columns) + todo = reconcile.plan(output, GROUP, units, digest) + if todo.empty(): + return todo + reconcile.apply(output, GROUP, todo, units, _read, digest, COUNT_ATTR) + return todo + + +class ReconcileMachine(RuleBasedStateMachine): + """A campaign that changes under a catalogue that must keep up with it.""" + + @initialize() + def setup(self): + self.dir = Path( + __import__("tempfile").mkdtemp(prefix="reconcile-props-") + ) + self.output = self.dir / "cat.h5" + self.columns = COLUMN_SETS[0] + self.sources = {} # unit -> source path + self.expected = {} # unit -> array as last written + self.clock = 1_000_000_000_000_000_000 + # Compaction is only claimed of the rebuild path (a plan that removes + # or refreshes). An add-only plan copies the file and appends, so its + # layout carries whatever the previous writes left behind. + self.rebuilt = False + + def teardown(self): + __import__("shutil").rmtree(self.dir, ignore_errors=True) + + # --- the campaign's moves ------------------------------------------- + def _tick(self): + self.clock += 1_000_000_000 + return self.clock + + def _step(self, changed): + before = (self.output.stat().st_mtime_ns + if self.output.exists() else None) + todo = _build(self.output, self.sources, self.columns) + self.rebuilt = bool(todo.remove or todo.refresh) + if not changed and before is not None: + assert self.output.stat().st_mtime_ns == before, ( + "a no-op reconcile rewrote the file; mtime is a rerun trigger" + ) + + @rule(pick=st.integers(0, 2**16), rows=st.integers(1, 5), + seed=st.integers(0, 2**16)) + @precondition(lambda self: len(self.sources) < len(UNITS)) + def add_unit(self, pick, rows, seed): + free = sorted(set(UNITS) - set(self.sources)) + unit = free[pick % len(free)] + path = self.dir / f"{unit}.npy" + array = _array(self.columns, rows, seed) + _write_source(path, array, self._tick()) + self.sources[unit] = path + self.expected[unit] = array + self._step(changed=True) + + @rule(pick=st.integers(0, 2**16), rows=st.integers(1, 5), + seed=st.integers(0, 2**16), resize=st.booleans()) + @precondition(lambda self: bool(self.sources)) + def modify_source(self, pick, rows, seed, resize): + unit = sorted(self.sources)[pick % len(self.sources)] + old = self.expected[unit] + rows = rows if resize else len(old) + array = _array(self.columns, rows, seed) + _write_source(self.sources[unit], array, self._tick()) + self.expected[unit] = array + self._step(changed=True) + + @rule(pick=st.integers(0, 2**16)) + @precondition(lambda self: bool(self.sources)) + def remove_unit(self, pick): + unit = sorted(self.sources)[pick % len(self.sources)] + self.sources.pop(unit).unlink() + self.expected.pop(unit) + self._step(changed=True) + + @rule() + def change_columns(self): + """Flip to the other column set — a digest change, so every unit + refreshes.""" + columns = next(c for c in COLUMN_SETS if c != self.columns) + self.columns = columns + # A schema change is a change to how the SOURCES are read, so the + # sources are rewritten under the new column set as the campaign would. + for i, (unit, path) in enumerate(sorted(self.sources.items())): + array = _array(columns, len(self.expected[unit]), 4242 + i) + _write_source(path, array, self._tick()) + self.expected[unit] = array + self._step(changed=True) + + @rule() + def no_op(self): + self._step(changed=False) + + # --- what must be true after every step ------------------------------ + @invariant() + def file_matches_campaign(self): + if not self.expected: + return + assert self.output.exists() + with h5py.File(self.output, "r") as f: + assert set(f[GROUP]) == set(self.expected), ( + "datasets and campaign units disagree") + assert f.attrs[COUNT_ATTR] == len(self.expected) + assert (f.attrs["param_digest"] + == reconcile.schema_digest(self.columns)) + dtypes = set() + for unit, want in self.expected.items(): + got = f[GROUP][unit][...] + assert got.dtype.names == want.dtype.names + np.testing.assert_array_equal(got, want) + dtypes.add(got.dtype) + stamp = reconcile.stamp(self.sources[unit]) + assert (int(f[GROUP][unit].attrs["src_bytes"]), + int(f[GROUP][unit].attrs["src_mtime_ns"])) == stamp + assert len(dtypes) == 1, ( + "sources share a column list; datasets must share a dtype") + + @invariant() + def compact(self): + """A rebuild does not carry the old file's dead space forward. + + HDF5 never reclaims a deleted dataset's space, which is why ``apply`` + builds the tmp FRESH whenever a plan removes or refreshes anything + instead of copying and editing in place. If that path stopped firing, + a long-lived campaign would grow by one unit per refresh forever. + + The bound is a from-scratch build of the same campaign plus a fixed + allowance: moving a dataset across with h5py's group copy costs a + little more metadata than creating it from an array does, measured at + ~1.4 kB here and — see the cycle test below — independent of how many + times the file has been rebuilt. What must never hold is growth that + tracks the history. + """ + if not self.rebuilt or not self.output.exists(): + return + fresh = self.dir / "fresh.h5" + fresh.unlink(missing_ok=True) + try: + _build(fresh, self.sources, self.columns) + if not fresh.exists(): + return + assert (self.output.stat().st_size + <= fresh.stat().st_size + COPY_SLACK), ( + "a rebuilt file is carrying dead space: " + f"{self.output.stat().st_size} bytes against " + f"{fresh.stat().st_size} from scratch") + finally: + fresh.unlink(missing_ok=True) + + +ReconcileMachine.TestCase.settings = settings( + max_examples=150, + stateful_step_count=14, + deadline=None, + suppress_health_check=[HealthCheck.too_slow, HealthCheck.data_too_large], +) +TestReconcileMachine = ReconcileMachine.TestCase + + +def test_crash_between_tmp_and_replace_leaves_the_file_untouched(): + """A failed rename must leave the previous catalogue byte-identical. + + Not a hypothesis case: the interesting axis is the crash point, and there + is one. ``os.replace`` is made to raise where the tmp is moved into place. + """ + import shutil + import tempfile + + work = Path(tempfile.mkdtemp(prefix="reconcile-crash-")) + try: + output = work / "cat.h5" + columns = COLUMN_SETS[0] + sources = {} + for i, unit in enumerate(UNITS[:2]): + path = work / f"{unit}.npy" + _write_source(path, _array(columns, 3, i), + 1_000_000_000_000_000_000 + i) + sources[unit] = path + _build(output, sources, columns) + before = output.read_bytes() + before_mtime = output.stat().st_mtime_ns + + # A third unit arrives, and the rename fails. + path = work / "u2.npy" + _write_source(path, _array(columns, 3, 99), 1_000_000_000_000_000_099) + sources["u2"] = path + units = sorted(sources.items()) + digest = reconcile.schema_digest(columns) + todo = reconcile.plan(output, GROUP, units, digest) + assert todo.add == ["u2"] + + real_replace = Path.replace + + def boom(self, target): + raise OSError("simulated crash between write and rename") + + Path.replace = boom + try: + with pytest.raises(OSError): + reconcile.apply(output, GROUP, todo, units, _read, digest, + COUNT_ATTR) + finally: + Path.replace = real_replace + + assert output.read_bytes() == before, "the old catalogue was modified" + assert output.stat().st_mtime_ns == before_mtime + assert not (work / "cat.h5.tmp").exists(), "tmp outlived the failure" + finally: + shutil.rmtree(work, ignore_errors=True) + + +def test_repeated_refresh_does_not_grow_the_file(): + """The leak the rebuild path exists to prevent, asserted directly. + + Twelve refreshes of one unit in a two-unit campaign. If ``apply`` ever + copied the file and edited it in place, each would strand the previous + dataset's bytes and the size would climb monotonically. + """ + import shutil + import tempfile + + work = Path(tempfile.mkdtemp(prefix="reconcile-growth-")) + try: + output = work / "cat.h5" + columns = COLUMN_SETS[0] + sources = {} + for i, unit in enumerate(("u0", "u1")): + path = work / f"{unit}.npy" + _write_source(path, _array(columns, 4, i), 10**18 + i) + sources[unit] = path + _build(output, sources, columns) + + sizes = [] + for k in range(12): + _write_source(sources["u0"], _array(columns, 4, 100 + k), + 10**18 + 100 + k) + todo = _build(output, sources, columns) + assert todo.refresh == ["u0"], todo.describe() + sizes.append(output.stat().st_size) + assert len(set(sizes)) == 1, f"file size drifted across refreshes: {sizes}" + finally: + shutil.rmtree(work, ignore_errors=True) diff --git a/tests/unit/test_persist_exp_props.py b/tests/unit/test_persist_exp_props.py new file mode 100644 index 000000000..4fc0bc039 --- /dev/null +++ b/tests/unit/test_persist_exp_props.py @@ -0,0 +1,360 @@ +"""Property-based state machine over ``workflow/scripts/persist_exp.py``. + +``exp_persist`` packs one exposure's keepable PSF products into a tar on +/project and writes a manifest describing it, and its central promise is that +RETENTION IS ADDITIVE: an existing tar is a floor, so shrinking the campaign's +keep list can never delete a product from the backed-up filesystem. That is a +claim about every sequence of keep lists and store states the campaign can +walk through, so it is tested here against a model — random keep lists over a +random set of present products, packed repeatedly, with the tar and the +manifest asserted after every pack. + +The keep lists mix product NAMES (``psf_model``), RAW GLOBS (``*.fits``) and +overlapping combinations of the two, because overlap is the case that once +failed every exposure in a campaign: two patterns matching one file is one +file, not a name collision. A genuine collision — two DIFFERENT source paths +landing on one flat member name — must still be fatal, and has its own test. + +The script is driven through ``main()`` with a patched ``sys.argv`` rather than +a subprocess: the rule invokes it as a script, but a subprocess per hypothesis +step would put this file out of reach of a login node's time budget. +""" + +import fnmatch +import hashlib +import importlib.util +import json +import shutil +import sys +import tarfile +import tempfile +from pathlib import Path + +import pytest +from hypothesis import HealthCheck, given, settings +from hypothesis import strategies as st +from hypothesis.stateful import ( + RuleBasedStateMachine, + initialize, + precondition, + rule, +) + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPTS = REPO_ROOT / "workflow" / "scripts" + + +def _load(name): + path = SCRIPTS / f"{name}.py" + assert path.exists(), f"{path} not found; the rule calls it by path" + sys.path.insert(0, str(SCRIPTS)) + try: + spec = importlib.util.spec_from_file_location(f"_{name}", path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + finally: + sys.path.remove(str(SCRIPTS)) + return module + + +persist = _load("persist_exp") +ALWAYS = persist.ALWAYS + +# One concrete file name per catalogued product, in the module output dir the +# real chain writes it to. The names are shaped like the campaign's (module +# tag, exposure, CCD) so the catalogue's globs match them for the same reason +# they match the real thing. +LAYOUT = { + "star_selection": ("setools", "mask", "star_selection-2079614-5.fits"), + "star_train": ("setools", "rand_split", + "star_split_ratio_80-2079614-5.fits"), + "star_test": ("setools", "rand_split", + "star_split_ratio_20-2079614-5.fits"), + "star_stats": ("setools", "stat", "star_stat-2079614-5.txt"), + "psf_model": ("psfex", "", "star_split_ratio_80-2079614-5.psf"), + "psfex_cat": ("psfex", "", "psfex_cat-2079614-5.cat"), + "psf_validation": ("psfex_interp", "", "validation_psf-2079614-5.fits"), +} +OPTIONAL = sorted(set(LAYOUT) - {ALWAYS}) +# What a campaign can write in `persist_exp:` — names, raw globs, and one name +# the catalogue does not know, which must be refused before any work happens. +ENTRIES = OPTIONAL + ["*.fits", "*.psf", "star_*", "validation_psf-*.fits"] +UNKNOWN = "psf_residuals" + +EXP = "2079614" + + +def _md5(path: Path) -> str: + return hashlib.md5(path.read_bytes()).hexdigest() + + +def _members(tar: Path) -> list: + with tarfile.open(tar) as tf: + return [ti.name for ti in tf.getmembers() if ti.isfile()] + + +class Store: + """One exposure's scratch store, its destination, and how to pack it.""" + + def __init__(self): + self.root = Path(tempfile.mkdtemp(prefix="persist-exp-props-")) + self.exp_dir = self.root / "exp" / EXP + self.dest = self.root / "products" / "psf" + self.manifest = self.root / "products" / "manifests" / f"{EXP}.json" + self.tar = self.dest / f"{EXP}.tar" + + def close(self): + shutil.rmtree(self.root, ignore_errors=True) + + def path_of(self, product: str) -> Path: + module, sub, name = LAYOUT[product] + base = (self.exp_dir / "output" / persist.RUN_NAME + / f"run_sp_{module}" / "output") + return (base / sub / name) if sub else (base / name) + + def write(self, product: str, payload: bytes) -> None: + path = self.path_of(product) + path.parent.mkdir(parents=True, exist_ok=True) + path.write_bytes(payload) + + def drop(self, product: str) -> None: + self.path_of(product).unlink(missing_ok=True) + + def pack(self, keep: list) -> int: + """Run the script's ``main`` as the rule does. 0 on success.""" + argv = ["persist_exp.py", "--exp-dir", str(self.exp_dir), + "--exp", EXP, "--dest", str(self.dest), + "--manifest", str(self.manifest)] + for entry in keep: + argv += ["--pattern", entry] + old = sys.argv + sys.argv = argv + try: + persist.main() + return 0 + except SystemExit as exc: + return 1 if exc.code not in (0, None) else 0 + finally: + sys.argv = old + + +class PersistExpMachine(RuleBasedStateMachine): + """A store that gains and loses products under a keep list that changes.""" + + @initialize() + def setup(self): + self.store = Store() + self.present = set() + self.prior_members = [] # members of the last tar written + self.prior_labels = {} # member -> product recorded for it + + def teardown(self): + self.store.close() + + # --- the store and the config move ---------------------------------- + @rule(product=st.sampled_from(sorted(LAYOUT)), size=st.integers(1, 64)) + def add_product(self, product, size): + self.store.write(product, bytes([len(product) % 251]) * size) + self.present.add(product) + + @rule(product=st.sampled_from(sorted(LAYOUT))) + def drop_product(self, product): + self.store.drop(product) + self.present.discard(product) + + @rule(keep=st.lists(st.sampled_from(ENTRIES), max_size=4, unique=True)) + def pack(self, keep): + self._pack_and_check(keep) + + @rule(keep=st.lists(st.sampled_from(ENTRIES), max_size=3, unique=True)) + def pack_with_unknown_product(self, keep): + """An unknown product name is refused before anything is written.""" + before = (_md5(self.store.tar) if self.store.tar.exists() else None) + code = self.store.pack(keep + [UNKNOWN]) + assert code != 0, "an unknown product name was accepted" + after = (_md5(self.store.tar) if self.store.tar.exists() else None) + assert after == before, "a refused keep list still touched the tar" + + @rule() + @precondition(lambda self: ALWAYS in self.present) + def pack_twice_unchanged(self): + """A rerun over an unchanged store must not move a single byte.""" + keep = sorted(OPTIONAL)[:2] + self._pack_and_check(keep) + tar_md5, man_md5 = _md5(self.store.tar), _md5(self.store.manifest) + tar_mtime = self.store.tar.stat().st_mtime_ns + man_mtime = self.store.manifest.stat().st_mtime_ns + assert self.store.pack(keep) == 0 + assert _md5(self.store.tar) == tar_md5, "the tar is not byte-stable" + assert _md5(self.store.manifest) == man_md5, "the manifest is not byte-stable" + assert self.store.tar.stat().st_mtime_ns == tar_mtime, ( + "an unchanged rerun rewrote the tar; mtime is a rerun trigger") + assert self.store.manifest.stat().st_mtime_ns == man_mtime, ( + "an unchanged rerun rewrote the manifest") + + # --- what a pack must leave behind ----------------------------------- + def _pack_and_check(self, keep): + had_tar = self.store.tar.exists() + tar_before = _md5(self.store.tar) if had_tar else None + man_before = (_md5(self.store.manifest) + if self.store.manifest.exists() else None) + code = self.store.pack(keep) + + if ALWAYS not in self.present: + # The star catalogue's input is not optional: the job fails and + # nothing downstream may be told the store is safe to reclaim. + assert code != 0, ( + f"{ALWAYS} is missing and the pack still succeeded") + assert (_md5(self.store.tar) if self.store.tar.exists() + else None) == tar_before, "a failed pack touched the tar" + assert (_md5(self.store.manifest) + if self.store.manifest.exists() + else None) == man_before, ( + "a failed pack wrote a manifest; clean_exposure would take " + "that as permission to delete the store") + return + + assert code == 0, f"pack failed with {ALWAYS} present and keep={keep}" + assert self.store.tar.exists() and self.store.manifest.exists() + members = _members(self.store.tar) + assert len(members) == len(set(members)), ( + f"duplicate member names in the tar: {members}") + + # ADDITIVE: an existing tar is a floor. + assert set(members) >= set(self.prior_members), ( + "members vanished from the tar: " + f"{sorted(set(self.prior_members) - set(members))}") + + body = json.loads(self.store.manifest.read_text()) + listed = {f["name"] for f in body["files"]} + assert listed == set(members), ( + "manifest and tar disagree about what was packed: " + f"{sorted(listed ^ set(members))}") + assert body["n_files"] == len(members) + assert body["unit"] == EXP and body["status"] == "complete" + + entries = [ALWAYS] + [e for e in keep if e != ALWAYS] + assert body["products"] == entries + for f in body["files"]: + if f["src"] is None: # carried from the previous tar + assert f["name"] in self.prior_members + assert f["product"] == self.prior_labels.get(f["name"], "?") + continue + assert f["product"] in entries, ( + f"{f['name']} labelled {f['product']!r}, not in the keep list") + assert fnmatch.fnmatch(f["name"], persist.resolve(f["product"])), ( + f"{f['name']} does not match {f['product']!r}'s glob") + assert Path(f["src"]).exists() + assert f["bytes"] == Path(f["src"]).stat().st_size + + # Every present product the keep list asks for is in there. + for entry in entries: + glob = persist.resolve(entry) + for product in self.present: + if fnmatch.fnmatch(LAYOUT[product][2], glob): + assert LAYOUT[product][2] in listed, ( + f"{product} matched {entry!r} but was not packed") + + self.prior_members = members + self.prior_labels = {f["name"]: f["product"] for f in body["files"]} + + +PersistExpMachine.TestCase.settings = settings( + max_examples=120, + stateful_step_count=12, + deadline=None, + suppress_health_check=[HealthCheck.too_slow, HealthCheck.data_too_large], +) +TestPersistExpMachine = PersistExpMachine.TestCase + + +@pytest.fixture() +def store(): + s = Store() + yield s + s.close() + + +def _seed(store, products=(ALWAYS,)): + for i, product in enumerate(products): + store.write(product, bytes([i + 1]) * (16 + i)) + + +def test_corrupt_existing_tar_is_refused_and_left_alone(store): + """A tar that cannot be read may still hold the only copy of something.""" + _seed(store, (ALWAYS, "psf_model")) + assert store.pack(["psf_model"]) == 0 + store.tar.write_bytes(b"not a tar at all, not even close" * 8) + corrupt = store.tar.read_bytes() + man_before = _md5(store.manifest) + + assert store.pack(["psf_model"]) != 0, "a corrupt tar was overwritten" + assert store.tar.read_bytes() == corrupt, "the corrupt tar was modified" + assert _md5(store.manifest) == man_before, ( + "a manifest was written over a tar that could not be read") + assert not store.tar.with_name(store.tar.name + ".tmp").exists() + + +def test_two_sources_with_one_member_name_is_fatal(store): + """Members are flat, so a real name clash would silently overwrite.""" + _seed(store, (ALWAYS,)) + # The same file name under a second module output dir. + clash = (store.exp_dir / "output" / persist.RUN_NAME / "run_sp_setools" + / "output" / "new_cat" / LAYOUT[ALWAYS][2]) + clash.parent.mkdir(parents=True, exist_ok=True) + clash.write_bytes(b"a different file with the same name") + + assert store.pack([]) != 0, "two different sources shared a member name" + assert not store.manifest.exists() + assert not store.tar.exists() + + +def test_shrinking_the_keep_list_cannot_delete_a_product(store): + """The property the additive rule exists for, stated end to end.""" + _seed(store, (ALWAYS, "psf_model", "star_train")) + assert store.pack(["psf_model", "star_train"]) == 0 + wide = set(_members(store.tar)) + assert LAYOUT["psf_model"][2] in wide + + # The campaign changes its mind, and the scratch store is gone. + for product in ("psf_model", "star_train"): + store.drop(product) + assert store.pack([]) == 0 + assert set(_members(store.tar)) == wide, ( + "shrinking persist_exp: deleted products from the backed-up tar") + body = json.loads(store.manifest.read_text()) + carried = {f["name"] for f in body["files"] if f["src"] is None} + assert LAYOUT["psf_model"][2] in carried + assert {f["name"]: f["product"] for f in body["files"]}[ + LAYOUT["psf_model"][2]] == "psf_model", ( + "a carried member lost the product label the old manifest had") + + +@settings(max_examples=80, deadline=None, + suppress_health_check=[HealthCheck.too_slow]) +@given(st.lists(st.sampled_from( + [ALWAYS, "*.fits", "validation_psf-*.fits", "psf_validation", + "star_*", "*.psf", "psf_model", "star_train"]), + min_size=1, max_size=5)) +def test_overlapping_patterns_never_fail(keep): + """Two patterns matching one file is one file, not a name collision. + + Overlap is ordinary — ``validation_psf-*.fits`` beside ``*.fits`` is a + perfectly reasonable way to say "the validation catalogues, and everything + else FITS while we are here" — and treating the second match as a clash + once failed every exposure in a campaign. + + A fresh store per example, because the additive rule makes packing + stateful and this property is about ONE pack. + """ + s = Store() + try: + _seed(s, tuple(LAYOUT)) + assert s.pack(keep) == 0, f"overlapping keep list failed: {keep}" + members = _members(s.tar) + assert len(members) == len(set(members)), members + # Every product present matched something, so all seven are packed. + assert set(members) == {name for _, _, name in LAYOUT.values()} & set( + members) + finally: + s.close() diff --git a/tests/unit/test_star_cat_columns.py b/tests/unit/test_star_cat_columns.py new file mode 100644 index 000000000..20975815d --- /dev/null +++ b/tests/unit/test_star_cat_columns.py @@ -0,0 +1,95 @@ +"""The star catalogue's 16 columns are defined twice, and must not drift. + +Two writers emit a full_starcat, for two consumers that have to agree about it: + + * ``MergeStarCatPSFEX`` (``src/shapepipe/modules/merge_starcat_package``), + which the ``merge_starcat`` MODULE RUNNER calls, writing the flat FITS table + sp_validation opens today; + * ``workflow/scripts/merge_star_cat.py``, the Snakemake workflow's + ``star_cat_merge`` rule, writing the per-exposure hdf5 that replaces it + (CosmoStat/sp_validation#340 moves the readers). + +They were one definition until the workflow stopped calling the module class: +the rule reads validation_psf members out of the per-exposure tars, keeps their +native dtypes and reconciles its output, none of which the class does or should +do. Two implementations is the right answer for the behaviour; two COLUMN LISTS +is not, and nothing else would notice them diverging — a column added to one +writer would simply be absent from the other's product, discovered by whoever +next tried to compute rho statistics from the wrong one. + +Hence this module, which asserts the one thing they must share. It does NOT +assert the dtypes: the whole point of the hdf5 writer is that they differ (the +FITS one widens every float to 1D). Only the names, and their order. + +Deliberately import-light on the workflow side: merge_star_cat.py pulls in h5py +and astropy, which the class does too, so a container-free run is not on offer +here and is not worth contorting for. +""" + +import importlib.util +import sys +from pathlib import Path + +import pytest + + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPTS = REPO_ROOT / "workflow" / "scripts" +SCRIPT = SCRIPTS / "merge_star_cat.py" + + +def _load_workflow_merge(): + """Import the rule's script by path — ``workflow/scripts`` is not a package. + + Its own imports (build_index, hdf5_reconcile, persist_exp) are siblings it + reaches through ``sys.path[0]``, which is how the rule invokes it, so the + directory goes on the path here too. + """ + assert SCRIPT.exists(), f"{SCRIPT} not found; the rule calls it by path" + sys.path.insert(0, str(SCRIPTS)) + try: + spec = importlib.util.spec_from_file_location("_merge_star_cat", SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + finally: + sys.path.remove(str(SCRIPTS)) + return module + + +@pytest.fixture(scope="module") +def writers(): + """The two column lists, each in the order its writer emits them.""" + h5py = pytest.importorskip("h5py") # noqa: F841 - workflow dep + pytest.importorskip("astropy") + workflow = _load_workflow_merge() + from shapepipe.modules.merge_starcat_package.merge_starcat import ( + MergeStarCatPSFEX, + ) + # The class carries (output name, source column) pairs plus its optional + # set and appends CCD_NB last; the script carries output names throughout. + module_columns = ( + tuple(out for out, _ in MergeStarCatPSFEX._COLUMNS) + + tuple(out for out, _ in MergeStarCatPSFEX._OPTIONAL) + + ("CCD_NB",) + ) + return module_columns, tuple(workflow.ALL_COLUMNS) + + +def test_column_names_and_order_agree(writers): + """Same names, same order — the schema both products promise.""" + module_columns, workflow_columns = writers + assert workflow_columns == module_columns + + +def test_sixteen_columns(writers): + """The count is itself the documented contract (README, config.yaml).""" + module_columns, workflow_columns = writers + assert len(module_columns) == 16 + assert len(workflow_columns) == 16 + + +def test_ccd_nb_is_last(writers): + """CCD_NB is appended per input file rather than read from one, in both.""" + module_columns, workflow_columns = writers + assert module_columns[-1] == "CCD_NB" + assert workflow_columns[-1] == "CCD_NB" diff --git a/workflow/README.md b/workflow/README.md index 55611cedc..c00d66c18 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -107,7 +107,9 @@ and the run fails if either phase failed. ## The launch code snapshot `sp run` copies the code it is about to launch — `workflow/` (config symlinks -dereferenced), `src/` and the profile — into `/code`, records HEAD +dereferenced), `src/`, the repo's `scripts/` (`final_cat_merge` loads +`scripts/python/create_final_cat.py` by path) and the profile — into +`/code`, records HEAD plus a dirty flag in `/code/snapshot.json`, and runs the campaign entirely out of that copy. It matters because a campaign is not one process: the SLURM executor re-invokes snakemake on every job's node, so jobs re-parse the @@ -156,15 +158,18 @@ workflow/ bin/sp committed launcher (module load + /project venv + launch code snapshot + run/report/container/cancel) rules/ prepare.smk tile get_images/uncompress/find_exposures - exposure.smk per-exposure: get_images, split, psf (no temp()) - tile.smk per-tile: exp forest, merge_headers, detect, vignets, ngmix, merge, make_cat + exposure.smk per-exposure: get_images, split, psf, persist (no temp()); campaign star_cat_merge + tile.smk per-tile: exp forest, merge_headers, detect, vignets, ngmix, merge, make_cat; campaign final_cat_merge scripts/ - sp_rule.py the thin per-unit wrapper (isolation furniture, config copy, log-sync, count check) build_index.py prepare-phase run_index.sqlite builder (plain script) build_forest.py per-tile exposure symlink forest (group-compatible shell) completeness.py the ported count table (shared by sp_rule + run_report) run_report.py standalone report (NOT a DAG node; run_report hooks call it) container.py image layers + the resolution order behind `sp container` (stdlib-only) + persist_exp.py ONE exposure's keepable PSF products -> one tar on products_dir (the exp_persist rule) + hdf5_reconcile.py bring an hdf5 catalogue into agreement with a campaign (shared by both merges) + merge_star_cat.py ALL exposures' validation_psf, out of the tars -> full_starcat_.hdf5 + merge_final_cat.py ALL tiles' final_cat -> final_cat_.hdf5 (the final_cat_merge rule) clean_exposure.py ONE exposure's store + manifests + logs -> tombstone (the clean_exposure rule) profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; keep-going ``` @@ -240,6 +245,99 @@ profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; kee trigger reads that cut as a reason to rerun the very tiles it protects. Know the consequence — `--forcerun` on a tile whose `final_cat` exists will not rebuild its reclaimed exposures. Delete the `final_cat` first. +- **PSF products leave scratch before the purge does.** `exp_persist` packs + the products named by `persist_exp:` in `config.yaml` from the exposure's + scratch store into ONE uncompressed tar, + `/exp///psf/.tar` (inodes, not bytes, bind + on /project), and writes ONE manifest beside it recording the patterns, the + members and their sizes. The + threat it answers is the /scratch purge, not `clean_exposure` — the store goes + in 60 days whether or not the workflow reclaimed it — so it runs even with + `clean: false`, requested directly by `rule all`. `clean_exposure` takes its + manifest as an input, so reclamation can never overtake the copy. It is a + rule of its own rather than a `cp` on the end of `exp_psf` because the keep + list rides on `params`: adding a pattern reruns seconds of packing, not four + hours of PSF fitting per exposure. A pattern that matches nothing is a + recorded warning (setools rejects sparse CCDs); matching nothing at all is a + failure. A `localrule`, by the same arithmetic as `clean_exposure`. +- **The star catalogue's inputs are always kept; `persist_exp:` is what you + keep on top.** `exp_persist` packs `psf_validation` — the psfex_interp + validation catalogue, one per CCD — for every exposure whatever the config + says, because `star_cat_merge` stacks exactly those into the campaign's + `full_starcat`. They are that catalogue's provenance, and they are what keeps + appending a tile next month cheap rather than a rebuild from VOS. About 2 MB + per exposure: ~40 GB and ~40k inodes at DR6 scale, against a ~1 M-inode group + quota. `persist_exp:` is purely additive, and an empty list is legal — the tar + then holds the merge's inputs and nothing else. +- **The keep list names products, not globs.** Entries are names from a + catalogue in `workflow/scripts/persist_exp.py`, which is the single source of + truth for what each one means and what keeping it buys + ([#844](https://github.com/CosmoStat/shapepipe/issues/844)); `config.yaml`'s + block is that catalogue rendered, and `persist_exp.py --list-products` prints + it. Sizes are per exposure, 40 CCDs, measured on smk-m2. + + | product | glob | per exposure | what it buys | + |---|---|---|---| + | `psf_model` | `*.psf` | 2.8 MB | re-interpolate the PSF anywhere later, no rebuild | + | `psfex_cat` | `psfex_cat-*.cat` | unmeasured | which stars PSFEx's outlier rejection clipped | + | `star_selection` | `star_selection-*.fits` | 24.5 MB | which stars the selection cuts rejected, and why | + | `star_train` | `star_split_ratio_80-*.fits` | 19.9 MB | the 80% sample PSFEx fitted | + | `star_test` | `star_split_ratio_20-*.fits` | 7.1 MB | the 20% sample `psf_validation` corresponds to | + | `star_stats` | `star_stat-*.txt` | unmeasured | setools' per-CCD counts, density and FWHM cuts | + + The default is `psf_model`. `psf_validation` is in the catalogue too but needs + no naming; naming it anyway is harmless. **Retention is additive**: an + existing tar is a floor, so shrinking the list adds nothing and removes + nothing. Dropping a product is a deliberate act on `products_dir`, not a + config edit — otherwise editing a config would delete products from the + backed-up filesystem whose scratch originals are long gone. A raw glob is still accepted as an + escape hatch — anything with a glob metacharacter or a dot is read as one — + and an unknown *name* is a parse-time error listing the valid ones. The list + is exposure-side only; tile-side retention is #844 follow-up. +- **The campaign ends in two merged catalogues, and the workflow makes both.** + Everything above is per unit; the two products downstream analysis actually + opens are per *campaign*, and until these rules existed each was a manual pass + after the run. + `star_cat_merge` collects every exposure's every CCD's `psf_validation` into + `/full_starcat_.hdf5`, one dataset per exposure at + `exposures/` — the rho/tau statistics input. It reads the members + straight out of the per-exposure tars (`tarfile`; unpacking ~800k files to + merge them would defeat the tar's whole purpose), keeps their native dtypes, + and stores `CCD_NB` as an int. sp_validation still opens the old flat FITS + name, `full_starcat-0000000.fits`; its readers move to this file under + [sp_validation#340](https://github.com/CosmoStat/sp_validation/issues/340), + the same migration that retires the `patches/` key on the tile side. The rule + exists whenever the campaign has a persisted exposure. + **Two writers, one schema.** The module runner still emits the flat FITS + table through `MergeStarCatPSFEX`, and this rule emits the hdf5; they are + separate implementations on purpose, because only one of them reads tars, + keeps native dtypes and reconciles. Their 16 COLUMN NAMES must not drift + apart, and nothing else would notice if they did — a column added to one + writer would just be missing from the other's product. `tests/unit/` + `test_star_cat_columns.py` is what holds them together. + `final_cat_merge` collects every ready tile's `final_cat-.fits` into + `/final_cat_.hdf5`: one dataset per tile under a group + named for the campaign, the `final_cat.param` columns, an `n_tiles` attribute. + That schema is what sp_validation's reader opens, so it is fixed; the column + extraction reuses `scripts/python/create_final_cat.py` while the file is + written here, because that script's own discovery walks a directory layout + this workflow does not have. `campaign:` in `config.yaml` names the group and + defaults to the persistent root's basename. + BOTH RECONCILE, through one shared module (`hdf5_reconcile.py`) so the + campaign's two products cannot disagree about what an output owes its inputs. + Each adds the units that have no dataset, drops datasets whose unit left the + campaign, re-reads one whose source changed (every dataset records its + source's size and mtime) or whose column set moved (a digest on the file's + root), and leaves the rest unread — because re-reading a campaign to add one + unit is ~800 GB of IO at DR6 scale. The *content* is still a function of the + input set; the byte layout is not, and a no-op leaves the file untouched + rather than rewritten. + Both rerun when the set changes: the unit ids' fingerprint rides on `params`. + Neither is a `localrule` — one job over ~20k units is real work — and neither + puts its input paths in its shell, which is not fastidiousness: ~20k paths is + an order of magnitude over Linux's 128 KiB `MAX_ARG_STRLEN` for a single argv + entry, so each job is handed the tile list and the run index and derives the + same set from them. - **A dead tile can be told to stop pinning exposures.** An exposure is cleanable only once every consuming tile has its vignets, so one permanently-failed tile holds its ~80 exposures for the life of the diff --git a/workflow/Snakefile b/workflow/Snakefile index bbdfeefb6..1a45598ed 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -31,6 +31,7 @@ the manifest says "this stage succeeded", the log says "here is what happened" (the contract is argued in completeness.py's docstring). """ +import fnmatch import functools import hashlib import json @@ -40,6 +41,10 @@ import sys from pathlib import Path from snakemake.exceptions import WorkflowError +# Explicit rather than relying on the name snakemake injects into this +# namespace: the one place we log at parse time is a branch that only a +# non-default keep list reaches, and a NameError there would be found by a user. +from snakemake.logging import logger # Resolved relative to THIS file, not the working directory: snakemake runs with # --directory on /scratch (bin/sp) so .snakemake/ state never lands on /project @@ -92,11 +97,24 @@ RUN_DIR = Path(OUTPUTS["run_dir"]) # second path: one root, exactly the pre-D5 layout. PRODUCTS_DIR = Path(OUTPUTS.get("products_dir") or RUN_DIR) INDEX_DB = Path(OUTPUTS["index_db"]) +# The campaign's NAME — what the two campaign-level merges label their output +# with (`final_cat_.hdf5`, and the group inside it that holds the +# campaign's per-tile datasets). It +# defaults to the persistent root's basename, which is already how every +# campaign here is named (smk-g4, smk-g5, smk-g6: run_dir, products_dir and +# index all end in it), so the common case needs no key at all. Set `campaign:` +# in config.yaml when the two must differ. +CAMPAIGN = config.get("campaign") or PRODUCTS_DIR.name SCRIPTS = Path(workflow.basedir) / "scripts" # The config chain is the repo's committed directory (D2). The configs and # rules that set their environment variables must be versioned together. There # is no `config_src` knob. CONFIG_DIR = Path(workflow.basedir) / "config" / "cfis" +# `sp run`'s code snapshot (bin/sp), a sibling of workflow.basedir inside +# $STATE_DIR/code. Read by the campaign-level merges so the products they write +# carry the code that produced them; absent when this Snakefile is driven +# outside `sp run`, which the merges tolerate (hdf5_reconcile.code_provenance). +SNAPSHOT_JSON = Path(workflow.basedir).parent / "snapshot.json" sys.path.insert(0, str(SCRIPTS)) import build_index # noqa: E402 @@ -123,8 +141,14 @@ if PHASE not in ("prepare", "compute", "passthrough"): raise WorkflowError( f"SP_PHASE={PHASE!r} is not one of prepare, compute, passthrough.") +# DEDUPED, order preserved. The tile list is appended to by hand across a +# campaign, so a tile can appear twice — harmless for a per-tile target, which +# is the same path requested twice, but not for the campaign-level merges: their +# fingerprint counts what is in this list and the job derives a deduped set from +# the same file (build_index.campaign_tiles), so a duplicate line would make the +# two disagree about a set they must name identically. with open(config["tile_list"]) as f: - TILES = [ln.strip() for ln in f if ln.strip()] + TILES = list(dict.fromkeys(ln.strip() for ln in f if ln.strip())) # --- ngmix scatter (D4) ---------------------------------------------------- # Native directive: `--set-scatter ngmix=N` overrides it, N=1 degenerates to one @@ -263,6 +287,7 @@ EXP_DIR = str(RUN_DIR / "exp" / "{shard}" / "{exp}") # The persistent root mirrors the scratch one, shard for shard, so the two trees # read as the same campaign seen from two filesystems. PROD_TILE_DIR = str(PRODUCTS_DIR / "tiles" / "{shard}" / "{tile}") +PROD_EXP_DIR = str(PRODUCTS_DIR / "exp" / "{shard}" / "{exp}") def tile_dir(tile): return f"{RUN_DIR}/tiles/{tile[:2]}/{tile}" @@ -276,6 +301,18 @@ def tile_manifest(tile, stage): def exp_manifest(exp, stage): return f"{exp_dir(exp)}/manifests/{stage}.json" +def prod_exp_dir(exp): + """The exposure's dir on the PERSISTENT root — where exp_persist writes. + + Sharded identically to the scratch one, so the two trees read as the same + campaign seen from two filesystems, exposure side as well as tile side.""" + return f"{PRODUCTS_DIR}/exp/{exp[:2]}/{exp}" + +def prod_exp_manifest(exp, stage): + """A manifest that must SURVIVE reclamation, so it is not in the exposure's + scratch manifests/ dir (clean_exposure deletes that wholesale).""" + return f"{prod_exp_dir(exp)}/manifests/{stage}.json" + def forest_dir(tile): return f"{tile_dir(tile)}/exp_forest" @@ -308,12 +345,52 @@ def unit_num(unit): # gets its own hash, on the forest rule only. def script_hash(name): """The 12-hex fingerprint of one script under workflow/scripts/.""" - return hashlib.md5((SCRIPTS / name).read_bytes()).hexdigest()[:12] + return path_hash(SCRIPTS / name) + + +def path_hash(path): + """The 12-hex fingerprint of any file the rules depend on but do not own. + + A rule whose behaviour comes from more than its own script needs all of it + in one trigger. final_cat_merge is the case: what it writes is decided by + scripts/python/create_final_cat.py (the column extraction) and by + config/cfis/final_cat.param (which columns), and NEITHER is under + workflow/scripts/ nor a declared input. Without them in the hash, this PR's + own edits to both would have left every finished campaign's hdf5 untouched + and nothing would have said so. + + A MISSING FILE IS NOT A PARSE ERROR. This runs at module level, so raising + here kills EVERY invocation of the workflow — `sp --unlock`, `sp report`, + a dry run — over a file only one rule needs. Worse, it kills them with a + bare FileNotFoundError, which is exactly the diagnosis merge_final_cat.py + already carries and would print if the job were allowed to reach it. So the + hash degrades to a sentinel and says so once; the parse survives, the rule + still exists, and the job fails with the message written for it. + """ + path = Path(path) + try: + return hashlib.md5(path.read_bytes()).hexdigest()[:12] + except OSError: + if workflow.is_main_process: + logger.warning( + f"missing: {path} — it is part of a rule's rerun trigger, so " + f"that rule cannot tell whether it is out of date. The job " + f"that needs the file will say so when it runs.") + return "missing" SCRIPT_HASH = script_hash("completeness.py") FOREST_HASH = script_hash("build_forest.py") CLEAN_HASH = script_hash("clean_exposure.py") CLEAN_TILE_HASH = script_hash("clean_tile.py") +PERSIST_HASH = script_hash("persist_exp.py") +MERGE_STAR_HASH = script_hash("merge_star_cat.py") +# Three files, one trigger: the rule's script, the column extraction it calls, +# and the parameter file that says which columns (path_hash argues why). +MERGE_FINAL_HASH = ":".join(( + script_hash("merge_final_cat.py"), + path_hash(Path(workflow.basedir).parent / "scripts" / "python" + / "create_final_cat.py"), + path_hash(CONFIG_DIR / "final_cat.param"))) # ngmix_range.py earns a hash for a stronger reason than the others. What it # emits is not a stale RESULT but a stale BOUNDARY, and a tile's eight chunks are # a PARTITION of its object IDs: resume a tile across an edit to the split and @@ -431,6 +508,426 @@ def clean_targets(): out.append(tombstone(exp)) return sorted(out) +# --- persisted exposure products (D5) -------------------------------------- +# The keep list is config, not a rule input, and it is READ HERE so that exactly +# one place converts it into the form the rule carries. +# +# OPTIONAL RETENTION, and only that. What star_cat_merge needs — every CCD's +# psf_validation — is packed by exp_persist whatever this list says +# (persist_exp.py's ALWAYS argues why: provenance for the merged catalogue, and +# a cheap tile append later). So an EMPTY list is a coherent instruction and not +# a switch that turns persistence off: the tar then holds the star catalogue's +# inputs and nothing else, and exp_persist still runs for every exposure. +PERSIST_EXP = list(config.get("persist_exp") or []) + +# The keep list names PRODUCTS (`psf_model`), not globs (`*.psf`); the +# catalogue that maps one to the other lives in persist_exp.py, which is also +# what the rule runs, so there is one definition and not a copy here. +# UNKNOWN NAMES DIE AT PARSE TIME, listing the valid ones — a typo in a keep +# list would otherwise be a silently-empty keep or a per-exposure failure an +# hour into a campaign. +import persist_exp as _persist # noqa: E402 + +for _entry in PERSIST_EXP: + try: + _persist.resolve(_entry) + except KeyError as _exc: + raise WorkflowError(f"config persist_exp: {_exc.args[0]}") + + +def persist_targets(): + """Which exposures this invocation must pack PSF products off scratch for. + + `rule all` requests these DIRECTLY rather than reaching them only through + clean_exposure. Persistence and reclamation are different concerns — the + /scratch purge takes the store whether or not `clean:` is on — and hanging + the copy off the clean rule alone would mean a campaign run with clean:false + persists nothing and loses everything at the purge. + + Scope is the ready tiles' exposures, which `all` already builds through the + tile chain, so nothing new is pulled into the DAG by asking. + + EXCEPT AN EXPOSURE WHOSE STORE IS GONE. Its exp_psf manifest is not there, + so requesting its persist manifest would make the DAG rebuild the whole + exposure chain from VOS — the avalanche tile.smk's reclaimed-edge cut exists + to prevent, arriving through a new target instead. exp_store_reclaimed() + below is that test, and what it does NOT test is the tombstone. + + HEAD PROCESS ONLY, for the same reason as clean_targets() above. + """ + if not workflow.is_main_process: + return [] + return persist_manifests() + + +@functools.lru_cache(maxsize=1) +def persist_manifests(): + """persist_targets() without the head-process guard, memoised. + + The guard on persist_targets() is a cost decision, not a correctness one: + `rule all` reads it at MODULE level, so a job parse would pay a whole + campaign's index walk for a target it can never schedule. star_cat_merge + reads the same list through an INPUT FUNCTION, which snakemake evaluates + only for parses that actually build that job — the head process, and the one + merge job's own re-parse under the slurm executor, which genuinely needs it. + So this half carries no guard and the memo keeps either parse to one walk. + """ + exps = {e for t in TILES_READY for e in tile_exposures(t)} + return sorted(prod_exp_manifest(e, "exp_persist") for e in exps + if not exp_store_reclaimed(e)) + + +def exp_store_reclaimed(exp): + """True when this exposure's PSF products exist ONLY on the persistent root. + + The one condition both the persist target list and star_cat_merge's edge + choice turn on, and it deliberately does NOT read the tombstone. + + THE TOMBSTONE ALONE IS NOT THE EVIDENCE. It says clean_exposure ran, and + clean_exposure is only one of the two ways a scratch store disappears. + + WHAT THE TEST HAS TO SEPARATE is a store that is GONE from one that has not + been BUILT yet, and no single file says that. This runs at parse time, before + any exp_psf job of a fresh campaign has run, so "the exp_psf manifest is + missing" alone would skip every exposure of a new campaign and persist + nothing at all. The question is therefore: is there evidence this exposure + once had a store? Two files carry it, and either will do: + + * the TOMBSTONE — clean_exposure ran, so the store was built and reclaimed, + and whatever was going to be packed was packed before it went; + * a PERSISTED MANIFEST with no exp_psf manifest beside it — exp_persist + ran, so the store existed, and it is not there now. This is the purge + case, and it is the one the tombstone cannot see: /scratch is purged on + a 60-day window whether or not this workflow reclaimed anything, and it + leaves nothing behind. Keying on the tombstone alone meant that after a + purge, or on any campaign run with `clean: false`, every exposure looked + live, its persist manifest was requested, its exp_psf manifest was not + there, and snakemake rebuilt the entire exposure chain from VOS. + + An exposure with a LIVE store and a manifest is not reclaimed and is still + asked for, which is what lets an edit to `persist_exp:` re-pack in seconds + rather than be silently ignored — the whole reason exp_persist is a rule of + its own. An exposure with neither file is asked for too: either it has not + run yet, or it was purged having saved nothing, and only the DAG can tell + those apart by trying. + """ + return (Path(tombstone(exp)).exists() + or (Path(prod_exp_manifest(exp, "exp_persist")).exists() + and not Path(exp_manifest(exp, "exp_psf")).exists())) + +# --- the campaign-level merges --------------------------------------------- +# Two rules, one job each per campaign, both writing to the persistent root, and +# both the LAST link of a chain whose per-unit half the workflow already had: +# the exposure side ends in one `full_starcat_.hdf5` (every CCD's PSF +# validation catalogue — the rho/tau statistics input) and the tile side in one +# `final_cat_.hdf5` (every tile's final catalogue — the shear +# catalogue sp_validation reads). Until they existed the workflow's product set +# was two files short of what the old `combine_runs.bash` + `create_final_cat.py` +# chain delivered, and every campaign ended with a manual merge. +# +# NEITHER IS A LOCALRULE, and the arithmetic runs the opposite way from +# exp_persist's. Those rules are ~20k jobs of seconds each, so submitting them +# costs more in scheduling latency than the work; these are ONE job each over the +# whole campaign — ~800k catalogues stacked in memory, or ~20k catalogues read +# end to end at DR6 scale. That is a compute job, and it belongs on a node. +# +# NEITHER PUTS ITS INPUT PATHS IN ITS SHELL. `{input}` at DR6 scale is ~20k paths +# in a single argv entry, an order of magnitude over Linux's 128 KiB +# MAX_ARG_STRLEN, and the job would die on exec. So each rule's `input` is the +# DAG EDGE (what must exist first) and each script rediscovers the same set from +# the tile list and the index; what travels is a FINGERPRINT of that set's unit +# ids, on `params`, which is what makes the merge rerun when the set changes and +# not otherwise. +# The scripts' docstrings argue the rediscovery — it is also what lets the merges +# cover exposures whose scratch stores reclamation has since taken. + + +def unit_fingerprint(units): + """A short digest of a set of UNIT IDS, for a merge rule's `params`. + + `params` is a rerun trigger and a set of ids is not: appending a tile grows + the set, moves the digest and reruns the merge, while a rerun over the same + set leaves it alone. Sorted before hashing because the ORDER is not part of + what changed. + + IDS RATHER THAN THE RULE'S `input` PATHS. A path can change while the set + does not — star_cat_merge's edge for one exposure flips from its manifest to + its tar when the store is reclaimed — and a merge that reruns over identical + content on every reclamation pass is a rerun trigger firing on bookkeeping. + The ids are also exactly what the job derives on its own side, so both + halves agree on the set and on how it is named. + """ + joined = "\n".join(sorted(str(u) for u in units)) + return f"{len(units)}:{hashlib.md5(joined.encode()).hexdigest()[:12]}" + + +# NO GATE ON THE KEEP LIST. star_cat_merge used to exist only when +# `persist_exp:` named something validation_psf-shaped, which made the +# campaign's star catalogue an opt-in and a typo away from silently absent. +# exp_persist now packs psf_validation unconditionally, so the merge is +# requested whenever the campaign has a persisted exposure at all, and +# star_cat_targets() below is the only condition left. + + +def full_starcat(): + """The campaign's merged star catalogue — the rho/tau statistics input. + + hdf5, one dataset per exposure, named for the campaign exactly as the shear + catalogue beside it is. The old flat FITS table it replaces was called + `full_starcat-0000000.fits` and sp_validation still opens that name; + CosmoStat/sp_validation#340 moves its readers to this file, the same + migration that retires the `patches/` key on the galaxy side.""" + return f"{PRODUCTS_DIR}/full_starcat_{CAMPAIGN}.hdf5" + + +def final_cat_hdf5(): + """The campaign's merged shear catalogue — sp_validation's galaxy_cat_path.""" + return f"{PRODUCTS_DIR}/final_cat_{CAMPAIGN}.hdf5" + + +def prod_exp_tar(exp): + """The tar exp_persist writes. Not a declared output of anything — see + star_cat_inputs().""" + return f"{prod_exp_dir(exp)}/psf/{exp}.tar" + + +@functools.lru_cache(maxsize=1) +def star_cat_inputs(): + """What star_cat_merge waits for: every exposure of TILES_READY whose PSF + products are on the persistent root, live and reclaimed alike. + + THE SAME SET merge_star_cat.py derives at job time, and that equality is + load-bearing — the fingerprint on `params` is taken over THIS list, so + anything the job stacked that was not in it would be rows no rerun trigger + could see. The job states the rule from its own side: same tile list, same + index, exp_persist manifest present on the persistent root. By the time it + runs, every exposure below has one. + + RECLAIMED EXPOSURES BELONG IN THE STAR CATALOGUE. Carrying their PSF + products off scratch is exactly what exp_persist is for, and a merge that + dropped them would shrink the campaign's star catalogue every time + reclamation ran. But their exp_psf manifest is gone, so REQUESTING their + exp_persist manifest rebuilds the whole exposure chain from VOS — the + avalanche persist_targets() drops them to avoid. + ancient() DOES NOT HELP: it suppresses the timestamp comparison, not the + missing input, and snakemake schedules the chain anyway. Measured on smk-g6 + with one reclaimed exposure given a manifest by hand: the dry run grew + exp_get_images, exp_split, exp_psf and exp_persist jobs. + + So a reclaimed exposure is depended on through its TAR instead. The tar is + not a declared output of any rule (exp_persist declares only its manifest, + deliberately — persist_exp.py says why), so a tar that exists is a DAG leaf: + snakemake requires it and builds nothing. A live exposure keeps its manifest + edge, which is what orders the merge after the packing; its tar does not + exist yet, so it could not serve as the edge. + + An exposure reclaimed by a workflow PREDATING exp_persist has neither tar nor + manifest and is in no set at all. Nothing short of rebuilding its chain from + VOS recovers it; the merge reports how many exposures it found. + """ + live, reclaimed = [], [] + for exp in sorted({e for t in TILES_READY for e in tile_exposures(t)}): + if not exp_store_reclaimed(exp): + live.append(prod_exp_manifest(exp, "exp_persist")) + elif Path(prod_exp_tar(exp)).exists(): + reclaimed.append(prod_exp_tar(exp)) + return live + reclaimed + + +@functools.lru_cache(maxsize=1) +def star_cat_exposures(): + """The exposure IDs star_cat_merge stacks — what its fingerprint is taken + over. + + THE IDS, NOT THE PATHS, and the difference is a rerun. An exposure's edge + FLIPS from its manifest to its tar the moment its store is reclaimed, so a + fingerprint over paths moves on every reclamation pass and reruns the merge + over content that did not change. The ids move only when the set does, which + is what the trigger is for. It is also what merge_star_cat.py derives on the + job side, so the two agree on the set AND on how it is named. + """ + return sorted(e for e in {e for t in TILES_READY for e in tile_exposures(t)} + if Path(prod_exp_manifest(e, "exp_persist")).exists() + or not exp_store_reclaimed(e)) + + +# --- sizing the two merges (D4) --------------------------------------------- +# MEASURED, not guessed, and measured as a SLOPE rather than a single number: +# these are the only two rules whose one job's footprint grows with the whole +# campaign, so a constant is wrong by however much the campaign is not the one +# it was tuned on. +# +# Both slopes were measured on this login node, inside the campaign container, +# against synthetic tars for the star side and against smk-g6's real +# catalogues for the tile side. Peak RSS is getrusage(RUSAGE_CHILDREN). +# +# STAR SIDE, AND IT IS FLAT IN THE CAMPAIGN. The merge writes one hdf5 dataset +# per exposure and reads one exposure at a time, so it is sized on the LARGEST +# exposure's members — ~2 MB — not on the campaign's. What follows is the +# history of how that came to be true, because the numbers are the argument. +# +# The FITS full_starcat this replaced was one flat table, so the job held the +# whole campaign. Two fixture points, 20 and 80 exposures of 40 CCDs x 400 stars +# (1.6 MB of members per exposure, against the 2.0 MB measured on smk-m2), +# across the rewrites this PR made to MergeStarCatPSFEX: +# +# input members python lists arrays+concat two passes +# 32.3 MB 383 MB 238 MB 221 MB +# 129.0 MB 1313 MB 740 MB 661 MB +# slope 10.1x 5.5x 4.8x +# +# The tenfold was one python float object (32 bytes) plus a list pointer (8) per +# 4 bytes of float32 payload. Arrays per catalogue removed that; counting rows +# from the headers and filling a preallocated array removed the rest. What +# remained at 4.8x was the OUTPUT: file_io writes every float column as FITS 1D, +# so float32 became a float64 table astropy then buffered. +# +# Per-exposure hdf5 removes the term entirely rather than shrinking it — and +# with it the ~240 GB a DR6-scale flat table would have wanted. Those +# improvements stay upstream regardless: the module runner still merges to one +# FITS table, and they are its fix. +# +# TILE SIDE, and it is the reassuring one. Two points against real smk-g6 +# catalogues, 2 tiles (73.9 MB in, largest 39.6 MB) and 6 tiles (235.5 MB in, +# largest 47.7 MB): peak RSS 129 MB and 139 MB. FLAT IN THE NUMBER OF TILES — +# the merge holds one catalogue at a time — so it is sized on the LARGEST tile, +# not the total, at ~3x it plus the interpreter. +# THE CEILING ON ANY REQUEST, and it is not a formatting nicety: a mem_mb above +# the partition maximum is a job SLURM will never schedule and snakemake will +# never diagnose — it sits PENDING with a reason nobody reads while the campaign +# looks alive. The two merge formulas grow with the campaign, so at some size +# they WILL cross it; capping turns "silently never runs" into "runs on the +# biggest node there is, and possibly dies with a diagnosable OOM". +# +# Nibi's standard compute node is 766 GB (192 cores, 4 GB/core); 750000 leaves +# room for the OS and the slurm accounting overhead. Override with `max_mem_mb:` +# for a cluster with smaller nodes, or to reserve headroom. +MAX_MEM_MB = int(config.get("max_mem_mb", 750_000)) +_capped_warned = set() + + +def capped_mem(mb, rule): + """min(mb, MAX_MEM_MB), and say so ONCE at parse time when it bites.""" + mb = int(mb) + if mb > MAX_MEM_MB: + if rule not in _capped_warned and workflow.is_main_process: + _capped_warned.add(rule) + logger.warning( + f"{rule}: sized at {mb} MB, capped to max_mem_mb={MAX_MEM_MB} " + f"(Nibi's standard node is 766 GB). The job will run with less " + f"memory than the measurement says it wants — expect an OOM, " + f"and split the campaign or fix the merge rather than raising " + f"this number past what a node has.") + return MAX_MEM_MB + return mb + + +STAR_MEM_BASE_MB = 500 # interpreter + astropy + h5py, rounded up +STAR_MEM_FACTOR = 6 # x the LARGEST exposure's members +FINAL_MEM_BASE_MB = 800 +FINAL_MEM_FACTOR = 4 # x the LARGEST tile; ~3 measured +# What one unit costs when its product is not on disk yet to be measured — a +# fresh campaign sizes its merge before anything has been packed or made. The +# exposure figure is psf_validation's alone (the only members the star merge +# reads), not a whole tar's; both are the measured medians in config.yaml's +# persist_exp block and the D5 notes. +EXP_BYTES_DEFAULT = 2_000_000 +# The product whose members star_cat_merge stacks — named once, here and in +# merge_star_cat.py, and resolved through persist_exp.py's catalogue. +STAR_CAT_PRODUCT = _persist.ALWAYS +STAR_CAT_PATTERN = _persist.resolve(_persist.ALWAYS) +TILE_BYTES_DEFAULT = 46_000_000 + + +def _size(path, default): + """Bytes on disk, or the documented per-unit default if it is not there.""" + try: + return Path(path).stat().st_size + except OSError: + return default + + +def star_cat_max_bytes(): + """The LARGEST exposure's psf_validation members — what sizes the merge. + + The star merge holds ONE exposure at a time now that its output is hdf5 + with a dataset per exposure, so its memory is flat in the campaign exactly + as the tile side's is. Sizing on the total would ask a node for a campaign's + worth of memory to hold ~2 MB. + """ + return max(_star_cat_exposure_bytes() or [EXP_BYTES_DEFAULT]) + + +def star_cat_bytes(): + """Total bytes of the members the star merge will actually read. + + THE TAR'S SIZE IS THE WRONG NUMBER, and increasingly wrong as the keep list + grows: the merge reads the psf_validation members and nothing else, while + the tar also holds whatever `persist_exp:` retains. With the default + retention that is 2.4x too much, and with the star_* products on it is ~36x + — a memory request that misses by more than an order of magnitude, and one + that would jump the moment an exposure got packed, since an unpacked one + contributed the per-exposure default instead. So the MANIFEST is read and + only the psf_validation members are counted; persist_exp records the product + each member came from, exactly so this is answerable without opening a tar. + + One json parse per exposure at DAG build, and only for the parse that builds + this job. An exposure not yet packed has no manifest and contributes the + measured default, which is the psf_validation figure and not the tar's. + """ + return sum(_star_cat_exposure_bytes()) + + +@functools.lru_cache(maxsize=1) +def _star_cat_exposure_bytes(): + """Per exposure, the bytes of the members the star merge will read.""" + out = [] + for exp in star_cat_exposures(): + manifest = Path(prod_exp_manifest(exp, "exp_persist")) + if not manifest.exists(): + out.append(EXP_BYTES_DEFAULT) + continue + try: + body = json.loads(manifest.read_text()) + # By product name or, for a manifest written before that field + # existed or by a raw-glob keep list, by file name — the same test + # merge_star_cat.is_member() applies, so the sizing counts exactly + # the members the job will read. + out.append(sum(f["bytes"] for f in body["files"] + if f.get("product") == STAR_CAT_PRODUCT + or fnmatch.fnmatch(f["name"], STAR_CAT_PATTERN))) + except (OSError, ValueError, KeyError): + out.append(EXP_BYTES_DEFAULT) + return out + + +def final_cat_max_bytes(): + """The LARGEST tile catalogue the hdf5 merge will read — what sizes it.""" + return max([_size(final_cat(t), TILE_BYTES_DEFAULT) for t in TILES_READY] + or [TILE_BYTES_DEFAULT]) + + +def star_cat_targets(): + """`full_starcat` when there is anything to stack into it, else nothing. + + One way to get nothing, and it is a state rather than an error: every + exposure in scope is already tombstoned — a + campaign resumed after reclamation, whose exposures were cleaned by a + workflow that predates exp_persist and therefore left neither tar nor + manifest to read. A rule with an empty input list would still be a JOB, and + it would write an empty star catalogue over a good one. + """ + if not workflow.is_main_process: + return [] + return [full_starcat()] if star_cat_inputs() else [] + + +def final_cat_targets(): + """The merged hdf5, whenever this campaign has a tile to put in it.""" + if not workflow.is_main_process or not TILES_READY: + return [] + return [final_cat_hdf5()] + # --- tile reclamation (D5) -------------------------------------------------- # A separate flag from `clean:` (config.yaml carries the full # argument): exposure reclamation costs nothing but a rebuild if a tile is @@ -601,11 +1098,21 @@ include: "rules/tile.smk" # localrule would: a local job cannot be fused into a submitted group. The old # star-catalogue rules were exactly that, and they are gone with the internal # mask generation.) -localrules: all, prepare_all_tiles, clean_exposure, clean_tile +# +# exp_persist joins them for the same arithmetic — one tar of a few MB per +# exposure, ~20k of them at DR6 scale, each far shorter than the scheduling +# latency that would submit it (exposure.smk argues the placement in full). It +# sits mid-chain between exp_psf and clean_exposure, but both of those are +# outside every group already (exp_psf is heavy, clean_exposure is local), so it +# adds no new grouping constraint. +localrules: all, prepare_all_tiles, clean_exposure, clean_tile, exp_persist rule all: input: [final_cat(t) for t in TILES_READY], + persist_targets(), + star_cat_targets(), + final_cat_targets(), clean_targets(), clean_tile_targets(), diff --git a/workflow/bin/sp b/workflow/bin/sp index fdc3d4f59..e3aec91a4 100755 --- a/workflow/bin/sp +++ b/workflow/bin/sp @@ -72,7 +72,10 @@ STATE_DIR="${SP_STATE_DIR:-${RUN_DIR}-state}"; mkdir -p "$STATE_DIR" # WHAT. `sp run` copies the code it is about to launch into $STATE_DIR/code and # runs the campaign entirely out of that copy: the Snakefile, the rules, the # scripts, the ini chain (symlinks DEREFERENCED -- workflow/config/cfis points -# into example/, and the copy must be self-contained), src/, and the profile. +# into example/, and the copy must be self-contained), src/, the repo's own +# scripts/ (final_cat_merge loads scripts/python/create_final_cat.py by path -- +# it is a script, not an installed module, and the hdf5 layout it defines must +# be pinned to the campaign like everything else here), and the profile. # Every workflow-internal path hangs off `workflow.basedir`, which IS the # snapshot, so they all follow it for free; the profile's PYTHONPATH pin is the # one that cannot (YAML splices nothing) and is rewritten below. @@ -91,10 +94,10 @@ snapshot_code() { mkdir -p "$SNAPSHOT" if command -v rsync >/dev/null 2>&1; then rsync -a --delete --copy-links --exclude '__pycache__' --exclude '*.egg-info' \ - "$HERE" "$REPO/src" "$REPO/profiles" "$SNAPSHOT/" + "$HERE" "$REPO/src" "$REPO/scripts" "$REPO/profiles" "$SNAPSHOT/" else rm -rf "$SNAPSHOT"; mkdir -p "$SNAPSHOT" - cp -rL "$HERE" "$REPO/src" "$REPO/profiles" "$SNAPSHOT/" + cp -rL "$HERE" "$REPO/src" "$REPO/scripts" "$REPO/profiles" "$SNAPSHOT/" find "$SNAPSHOT" -name __pycache__ -type d -prune -exec rm -rf {} + fi diff --git a/workflow/config.yaml b/workflow/config.yaml index 63b29413c..38979b2d0 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -57,6 +57,19 @@ outputs: # (-state; bin/sp explains why). products_dir: /project/def-mjhudson/cdaley/sp-products/smk-g6 +# The campaign's NAME. It labels the two campaign-level merges' output — +# /final_cat_.hdf5 and the group inside it holding that +# campaign's per-tile datasets — +# and nothing else; the per-unit stores are named by their own IDs. UNSET means +# the persistent root's basename, which is already how every campaign here is +# named (run_dir, products_dir and index_db all end in smk-g6), so this key only +# earns its place when the two must differ. +# +# It is not a rule input, so renaming a campaign mid-flight changes the merged +# catalogue's PATH and therefore builds a new one; the per-tile catalogues it +# reads are untouched. +# campaign: smk-g6 + # There is no config_src knob: the config chain is workflow/config/cfis, resolved # relative to the Snakefile. The configs interpolate $SP_RUN / $SP_UNIT_NUM / # $SP_CONFIG / $SP_EXP / $NGMIX_* and the rules export them -- configs and rules @@ -68,6 +81,100 @@ outputs: # would otherwise have to rebuild from tile headers. index_db: /project/def-mjhudson/cdaley/sp-products/smk-g6/index/run_index.sqlite +# OPTIONAL per-exposure retention: what to carry onto the persistent root ON TOP +# OF the star catalogue's own inputs, before the scratch store goes +# (`exp_persist`, exposure.smk). +# +# WHAT IS ALWAYS KEPT, AND IS NOT A CHOICE HERE: psf_validation, the psfex_interp +# validation catalogue, one per CCD. `star_cat_merge` stacks every one of them +# into /full_starcat_.hdf5, so they are that +# catalogue's PROVENANCE — a merged star catalogue with no per-exposure inputs +# beside it cannot be audited, re-cut, or recomputed after a purge — and they are +# what keeps APPENDING TILES CHEAP, since a tile added next month brings +# exposures whose catalogues must join the existing stack. ~2 MB per exposure: +# ~40 GB and ~40k inodes at DR6 scale, against a ~1 M-inode group quota. That is +# the price of being able to say where the number came from, and it is paid. +# +# SO THIS LIST IS PURELY ADDITIVE, and an empty one is a coherent instruction: +# the tar then holds the star catalogue's inputs and nothing else. +# +# Entries are PRODUCT NAMES — not globs. The catalogue below is the rendering of +# workflow/scripts/persist_exp.py's PRODUCTS table, which is the single source of +# truth for what each name means and what keeping it buys +# (CosmoStat/shapepipe#844); print it any time with +# +# workflow/bin/sp container exec python workflow/scripts/persist_exp.py --list-products +# +# product glob size/exposure +# -------------- -------------------------- ------------- +# star_selection star_selection-*.fits 24.5 MB +# setools' PRE-SPLIT selection. The only file that answers which +# stars the selection cuts rejected and why; the split samples have +# already lost the rejects. +# star_train star_split_ratio_80-*.fits 19.9 MB +# the 80% TRAINING sample, the stars PSFEx actually fitted. Rows +# duplicate star_selection. +# star_test star_split_ratio_20-*.fits 7.1 MB +# the 20% VALIDATION sample — the positions psf_validation's rows +# correspond to. Rows duplicate star_selection. +# star_stats star_stat-*.txt unmeasured +# setools' per-CCD STAT block: star counts, stars/deg^2, FWHM mode +# and cuts. The selection's summary without its catalogue. +# psf_model *.psf 2.8 MB +# the PSFEx model itself. Keeping it means the PSF can be +# re-interpolated at ANY position later without rebuilding the +# exposure chain — the single most capability-adding entry here. +# psfex_cat psfex_cat-*.cat unmeasured +# PSFEx's own output catalogue (FITS_LDAC): the per-star FLAGS_PSF +# and CHI2_PSF, i.e. WHICH stars outlier rejection clipped. Not +# recoverable from anything else — the .psf header keeps only the +# LOADED/ACCEPTED counts. +# psf_validation validation_psf-*.fits 2.0 MB +# the psfex_interp validation catalogue, one per CCD: the input to +# the rho/tau statistics, and to the star_cat_merge rule that stacks +# them into the campaign's full_starcat. +# +# A RAW GLOB IS STILL ACCEPTED, as an escape hatch for a file the catalogue does +# not name yet: anything carrying a glob metacharacter or a dot is taken as a +# glob rather than a name (`*.psf` is a glob, `psf_model` is the name for it). +# An unknown NAME is a parse-time error listing the valid ones. +# +# Matches are packed, flat, into ONE uncompressed tar per exposure: +# /exp///psf/.tar, with a manifest listing the +# members (and the product each came from) beside it. One tar rather than loose +# copies because inodes, not bytes, bind on /project. FITS members read straight +# from the tar: fits.open(io.BytesIO(tarfile.open(t).extractfile(m).read())). +# +# WHY COPY RATHER THAN EXEMPT THESE FROM CLEANUP. Reclamation is not the threat. +# run_dir is /scratch and is PURGED on a 60-day window whether or not +# clean_exposure ever ran; products_dir is /project, backed up and not purged. +# The only way a per-exposure product outlives its campaign is to leave the +# filesystem. (Ordering is free: clean_exposure takes the exp_persist manifest +# as an input, so a store is never reclaimed before its keepers are written.) +# +# EDITING THIS LIST IS CHEAP. It rides on exp_persist's `params`, so a change +# reruns the packing (seconds) and NOT exp_psf (four hours per exposure). That +# separation is the whole reason exp_persist is a rule of its own. +# +# THE DEFAULT is psf_model (2.8 MB per exposure on top of psf_validation's 2.0): +# the model that lets the PSF be re-interpolated at any position later without +# rebuilding the exposure chain from VOS, which is the single most +# capability-adding thing an exposure can keep. Add psfex_cat for a production +# run if you want to know which stars PSFEx clipped; the star_* products are for +# selection studies and cost an order of magnitude more. +# +# THIS LIST IS EXPOSURE-SIDE ONLY. Tile-side retention is not configurable: the +# only tile product that persists today is final_cat, written by tile_make_cat +# straight to products_dir. A tile keep list is #844 follow-up. +# +# NOTE ON products_dir DEFAULTING TO run_dir (a fixture or smoke test): the tar +# then lands beside the store on the same filesystem and buys nothing, and the +# manifest sits in the exposure's own manifests/ dir, which clean_exposure +# deletes wholesale — so a one-root run re-persists after every reclamation. +# Harmless, and exactly the pre-D5 behaviour a one-root run asks for. +persist_exp: + - psf_model + # 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, @@ -132,6 +239,15 @@ clean_tiles: true # it is dead, not while you are still debugging it. clean_ignore_tiles: [] +# The ceiling on any rule's mem_mb. A request above the partition maximum is a +# job SLURM never schedules and snakemake never diagnoses: it sits PENDING while +# the campaign looks alive. Nibi's standard compute node is 766 GB (192 cores at +# 4 GB/core), so 750000 leaves room for the OS and slurm's own overhead. The two +# campaign-level merges size themselves from the campaign's bytes and will cross +# this at survey scale — the cap turns "never runs" into "runs on the biggest +# node there is", with a parse-time warning saying which rule was capped. +max_mem_mb: 750000 + # ngmix within-tile chunking: static N chunks (closed ID ranges computed # per-tile, in-job, from the tile's own sexcat). ngmix_chunks: 8 diff --git a/workflow/config/cfis/final_cat.param b/workflow/config/cfis/final_cat.param index 00bcb3f73..f3fa39677 100644 --- a/workflow/config/cfis/final_cat.param +++ b/workflow/config/cfis/final_cat.param @@ -8,7 +8,26 @@ TILE_ID # flags FLAGS -IMAFLAGS_ISO +# NO IMAFLAGS_ISO, AND NO MASK COLUMN AT ALL — READ THIS BEFORE ADDING ONE. +# The tile-side SExtractor runs with FLAG_IMAGE = False and DOT_PARAM_FILE = +# default_noimaflags.param (config_tile_Sx.ini), so IMAFLAGS_ISO is never +# written into a tile catalogue and asking for it here only made the merge +# fail. Instrument flags reach the pipeline on the EXPOSURE side, where +# exp_split delivers the flag image and SExtractor reads it. +# +# Its intended replacement is make_cat's per-band MASK_ columns, queried +# from the sky-fixed healsparse maps named by MASK_EXT_PATHS. THE WORKFLOW SETS +# NO SUCH PATHS: config_tile_Mc.ini has no MASK_EXT_PATHS entry, so +# save_mask_ext_data is never called, no MASK_ column exists in any tile +# catalogue this workflow has produced, and smk-g6's carry none (checked). +# Naming one here would fail every merge on every campaign. +# +# So the merged catalogue carries NO mask information today, and that is a +# CONFIG gap and not a gap in this file: turning it on is setting +# MASK_EXT_PATHS in config_tile_Mc.ini (`band:path` pairs, the same grammar as +# the commented MASK_PATHS in config_exp_psfex.ini) and adding the matching +# MASK_ names here, in that order. No healsparse map is staged under +# /project/def-mjhudson yet. NGMIX_MCAL_FLAGS # PSF ellipticity (original image PSF) @@ -113,5 +132,6 @@ NGMIX_T_PSF_ORIG_NOSHEAR # PSF size measured on reconvolved image # NGMIX_T_PSF_RECONV_NOSHEAR -# ngmix moment failure flag -NGMIX_MOM_FAIL +# ngmix metacalibration type failure flag (renamed from NGMIX_MOM_FAIL in +# f0fca23e; catalogues written before that commit carry the old name) +NGMIX_MCAL_TYPES_FAIL diff --git a/workflow/rules/exposure.smk b/workflow/rules/exposure.smk index 4e1ef45a4..708840fdf 100644 --- a/workflow/rules/exposure.smk +++ b/workflow/rules/exposure.smk @@ -1,6 +1,6 @@ """Exposure chain — per exposure, keyed by exp base id (dedup is structural). - exp_get_images -> exp_split -> exp_psf + exp_get_images -> exp_split -> exp_psf -> exp_persist Each in the exposure's own sharded work dir, chained by manifests; every config reads fixed ``$SP_RUN/output/run_sp_exp_*`` INPUT_DIRs, so nothing resolves a @@ -17,6 +17,13 @@ the per-band ``MASK_`` columns on the tile side. Neither needs a rule, a star catalogue, or a network fetch — hence no ``star_catalogue`` / ``exp_star_cat`` here, and no ``exp_mask``. +``exp_persist`` is the one rule here that writes to the PERSISTENT root: it +packs the PSF products named by `persist_exp:` into one tar per exposure off +/scratch before the purge (or clean_exposure) can take them. It is a separate +rule from exp_psf precisely so that editing that list costs a re-pack and not a +four-hour refit; the full +argument is in workflow/scripts/persist_exp.py. + NO temp() anywhere in this file, ever (D5). Exposures overlap tiles by construction (~7-10 tiles each), so their consumer set closes over the CAMPAIGN, not over one invocation — reclamation here is clean_exposure's job (S5), driven @@ -106,6 +113,62 @@ rule exp_psf: sp_shell("exp_psf", f"config_exp_{PSF_MODEL}.ini") +# --- persistence (D5) ------------------------------------------------------- +# The counterpart of reclamation, and it must come first in the DAG: this packs +# the exposure's keepable PSF products into one tar on the persistent root, and +# clean_exposure below takes its manifest as an input so the store is never +# reclaimed before the keepers have left /scratch. The purge would take them +# anyway — that, not clean_exposure, is what this rule exists for +# (persist_exp.py's docstring argues both halves, and config.yaml's +# `persist_exp:` block carries the keep list and its candidates). +# +# A LOCALRULE (declared in the Snakefile), by exactly the arithmetic that made +# clean_exposure one: the body is a `tar` of a few MB from one shared filesystem +# to another, seconds of work, and one sbatch per exposure would be ~20k +# submissions at DR6 scale for jobs shorter than the scheduling latency. The +# grouping constraint that binds mid-chain localrules (this file's docstring) +# does not bite here: exp_persist's only neighbours are exp_psf, which is too +# heavy to ever fuse, and clean_exposure, which is local itself. +# +# ONE DECLARED OUTPUT, AND IT IS A MANIFEST, NOT THE TAR OR A directory(). The +# tar is not declared: a directory output would attest that a directory exists, +# where what we want written down is WHICH files were packed and how big each was — +# the provenance a rho-statistics run months from now needs in order to know +# what it is reading. The manifest is byte-stable, so a no-op rerun does not +# move its mtime and does not make clean_exposure look out of date. +# +# THE KEEP LIST RIDES ON params. That is the entire reason this is not three +# lines of tar appended to exp_psf's shell: `params` is a rerun trigger, so +# adding a pattern reruns the packing and leaves the PSF chain alone. +rule exp_persist: + input: + rules.exp_psf.output.manifest + output: + manifest = f"{PROD_EXP_DIR}/manifests/exp_persist.json" + # No `log:`: the script's only failure modes are "nothing matched" and a + # name collision, both of which it reports on stderr and neither of which + # has a per-CCD verdict worth a completeness record. + params: + # Only the OPTIONAL retention list travels: psf_validation is packed + # by persist_exp.py whatever this says. It still rides on params, so + # adding a product re-packs (seconds) rather than re-fitting the PSF. + patterns = " ".join(f"--pattern '{p}'" for p in PERSIST_EXP), + exp_dir = lambda wc: exp_dir(wc.exp), + dest = lambda wc: f"{prod_exp_dir(wc.exp)}/psf", + script_hash = PERSIST_HASH + threads: 1 + retries: 2 + resources: + mem_mb = 2000, + runtime = 10 + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/persist_exp.py" + " --exp-dir '{params.exp_dir}' --exp {wildcards.exp}" + " --dest '{params.dest}' --manifest {output.manifest}" + " {params.patterns}" + + # --- reclamation (D5) ------------------------------------------------------- # The one exception to "no reclamation in this file": clean_exposure OWNS # exposure-level deletion, and it is a real job, not temp() bookkeeping, because @@ -137,7 +200,14 @@ rule clean_exposure: # spatial neighbours. In-scope consumers keep their edge: they may run in # this DAG, so the clean must be ordered after them. lambda wc: [tile_manifest(t, "tile_vignets") - for t in clean_consumers(wc.exp) if t in READY_SET] + for t in clean_consumers(wc.exp) if t in READY_SET], + # The keepers must be off /scratch before the store goes. Unlike the + # consumer edges above, this edge does not depend on scope: it is the + # same exposure's own rule, so it drags nothing into the DAG that this + # exposure's chain did not already put there. It is UNCONDITIONAL now: + # exp_persist always packs the star catalogue's inputs, so there is no + # keep list under which this rule has nothing to wait for. + lambda wc: [prod_exp_manifest(wc.exp, "exp_persist")] output: tombstone = f"{EXP_DIR}/cleaned.json" params: @@ -151,3 +221,86 @@ rule clean_exposure: f"python {SCRIPTS}/clean_exposure.py" " --exp-dir $(dirname {output.tombstone}) --exp {wildcards.exp}" " --tombstone {output.tombstone} --consumers '{params.consumers}'" + + +# --- the campaign's star catalogue ------------------------------------------ +# ONE job per campaign: every exposure's every CCD's `validation_psf--.fits`, +# collected into `/full_starcat_.hdf5`, one dataset per +# exposure. That file is the rho/tau statistics input; the old bash chain built +# a flat FITS table with `combine_runs.bash psf` + a `merge_starcat_runner` +# pass, and the workflow emitted neither. sp_validation still opens the FITS +# name today — CosmoStat/sp_validation#340 moves its readers to this file, the +# same migration that retires the `patches/` key on the tile side. +# +# ONE DATASET PER EXPOSURE, NOT ONE TABLE, and it is the same decision as the +# tile side's: it makes the file RECONCILABLE. A flat table had to be restacked +# from every exposure the campaign had ever seen to add one — ~40 GB of members +# at DR6 scale to add ~2 MB — and held the whole campaign in memory while it did +# so. Reconciled, an append reads the appended exposures and nothing else, and +# the job holds one exposure at a time. hdf5_reconcile.py is the shared +# machinery; merge_star_cat.py argues the format and the tar reading. +# +# THE INPUT IS star_cat_inputs() (Snakefile): every exposure of TILES_READY whose +# PSF products are on the persistent root — the live ones through the exp_persist +# manifest edge `rule all` already requests, the RECLAIMED ones through their TAR, +# which no rule declares and which therefore requires nothing to be built. That +# asymmetry is not a flourish; requesting a reclaimed exposure's manifest +# rebuilds its whole chain from VOS, and ancient() does not prevent it (measured +# — the Snakefile carries the numbers). Nothing new enters the DAG either way. It +# is read through an INPUT FUNCTION rather than at module level so that only a +# parse which actually builds this job pays for the walk. +# +# THE PATHS DO NOT REACH THE SHELL, and that is not a style choice: ~20k manifest +# paths is an order of magnitude over Linux's 128 KiB MAX_ARG_STRLEN for a single +# argv entry, so `{input}` here would be a job that dies on exec at DR6 scale. +# The job is handed the two small files the Snakefile itself started from — the +# tile list and the index — and derives THE SAME SET from them; `params.inputs` +# carries that set's FINGERPRINT, which is the rerun trigger. The equality is +# the point: a job that stacked anything the fingerprint did not see would be +# rows no rerun trigger could notice, which is what a glob over products_dir +# would have given on a root shared with an earlier, larger tile list. +# +# NOT A LOCALRULE. exp_persist is local because it is 20k jobs of seconds; this +# is one job that reads the campaign's tars end to end. Its MEMORY is flat in +# the campaign (one exposure at a time) and sized on the largest exposure; its +# RUNTIME is the total. +# +# NO JOB AT ALL when every exposure in scope is tombstoned with no tar left +# behind: star_cat_targets() (Snakefile) simply does not request the output. +rule star_cat_merge: + input: + lambda wc: star_cat_inputs() + output: + star_cat = full_starcat() + params: + products_dir = str(PRODUCTS_DIR), + tile_list = str(config["tile_list"]), + index_db = str(INDEX_DB), + campaign = CAMPAIGN, + snapshot = str(SNAPSHOT_JSON), + inputs = unit_fingerprint(star_cat_exposures()), + script_hash = MERGE_STAR_HASH + threads: 1 + resources: + # Sized on the campaign's own member bytes, slope and intercept + # measured (the Snakefile's sizing block carries both points, and the + # ceiling this rule runs into at DR6 scale). Still * attempt, because a + # measured slope on synthetic tars is not a guarantee about real ones. + # Sized on the LARGEST exposure, not the total: the merge holds one + # exposure at a time (the Snakefile's sizing block carries the history). + mem_mb = lambda wc, attempt: capped_mem(attempt * ( + STAR_MEM_BASE_MB + + STAR_MEM_FACTOR * star_cat_max_bytes() // 1_000_000), + "star_cat_merge"), + # ~2 min per GB of members on the measurement above, doubled, over a + # floor that covers the fixed cost of opening ~40 members per exposure. + runtime = lambda wc, attempt: attempt * ( + 30 + 4 * star_cat_bytes() // 1_000_000_000) + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/merge_star_cat.py" + " --products-dir '{params.products_dir}'" + " --tile-list '{params.tile_list}' --index-db '{params.index_db}'" + " --output {output.star_cat}" + " --campaign '{params.campaign}'" + " --snapshot-json '{params.snapshot}'" diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index 68d6da9ba..51e1ff3f5 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -898,3 +898,77 @@ rule clean_tile: f"python {SCRIPTS}/clean_tile.py" " --tile-dir $(dirname {output.tombstone}) --tile {wildcards.tile}" " --tombstone {output.tombstone}" + + +# --- the campaign's shear catalogue ----------------------------------------- +# ONE job per campaign, the tile-side twin of exposure.smk's star_cat_merge, and +# the same three design calls hold: the input is the list `rule all` already +# requests (every ready tile's final_cat), the paths never reach the shell +# (MAX_ARG_STRLEN), and a fingerprint on `params` is what makes it rerun when a +# tile is appended. The job derives the same set the fingerprint was taken over +# from the tile list and the index rather than globbing products_dir — on a +# products root shared with an earlier, larger tile list a glob would merge tiles +# no rerun trigger ever saw. +# +# THE OUTPUT SCHEMA IS AN INTERFACE, NOT A CHOICE. sp_validation opens this file +# as its `galaxy_cat_path`: one dataset per tile under a named group, the +# columns of workflow/config/cfis/final_cat.param, an `n_tiles` attribute on the +# root. The group is named for the CAMPAIGN, which is the only unit this +# workflow has above the tile. So the rule reuses +# scripts/python/create_final_cat.py's column extraction rather than restating +# it, and writes the file itself — merge_final_cat.py argues that split, the one +# legacy literal in the schema, and the two places where the reference +# implementation had to be pinned down to be reproducible. +# +# THE INPUT IS final_cat, NOT the tile_make_cat manifest, for the same reason +# clean_tile's is: final_cat on the persistent root IS the campaign's +# tile-finished marker (see final_cat() in the Snakefile), and it is the file +# this rule actually reads. +# +# NOT A LOCALRULE, and here the reason is IO rather than memory: a first build +# reads every tile's catalogue end to end — ~32-46 MB per tile, so ~2 GB for a +# 64-tile campaign and ~800 GB at DR6's 23k tiles. It RECONCILES rather than +# rebuilds or appends: a tile with no dataset is added, a dataset whose tile +# left the campaign is deleted, a dataset whose source catalogue changed is +# re-read, and one that agrees with its source is left alone. So an append +# reads the appended tiles and nothing else, while the file still cannot drift +# from its inputs the way an append-only tool does (merge_final_cat.py argues +# what is and is not a function of the input set here). Memory is one tile's +# catalogue at a time plus the hdf5 write buffer, which is why mem_mb is modest +# where star_cat_merge's is not — and why runtime, which is sized on the whole +# campaign, is the pessimistic first-build case. +rule final_cat_merge: + input: + lambda wc: [final_cat(t) for t in TILES_READY] + output: + merged = final_cat_hdf5() + params: + products_dir = str(PRODUCTS_DIR), + tile_list = str(config["tile_list"]), + index_db = str(INDEX_DB), + param_file = str(CONFIG_DIR / "final_cat.param"), + campaign = CAMPAIGN, + snapshot = str(SNAPSHOT_JSON), + inputs = unit_fingerprint(TILES_READY), + script_hash = MERGE_FINAL_HASH + threads: 1 + resources: + # Sized on the LARGEST tile, not the total: the merge holds one + # catalogue at a time, and the measurement is flat in the tile count + # (the Snakefile's sizing block carries both points). + mem_mb = lambda wc, attempt: capped_mem(attempt * ( + FINAL_MEM_BASE_MB + + FINAL_MEM_FACTOR * final_cat_max_bytes() // 1_000_000), + "final_cat_merge"), + # Runtime, unlike memory, is the TOTAL: every tile is read end to end. + # ~1 min per 10 tiles on the measurement, triply generous, over a floor. + runtime = lambda wc, attempt: attempt * (30 + len(TILES_READY) // 3) + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/merge_final_cat.py" + " --products-dir '{params.products_dir}'" + " --tile-list '{params.tile_list}' --index-db '{params.index_db}'" + " --output {output.merged}" + " --campaign '{params.campaign}'" + " --param-file '{params.param_file}'" + " --snapshot-json '{params.snapshot}'" diff --git a/workflow/scripts/build_index.py b/workflow/scripts/build_index.py index 2dd191306..24290b561 100644 --- a/workflow/scripts/build_index.py +++ b/workflow/scripts/build_index.py @@ -152,6 +152,54 @@ def build(tile_ids: list[str], run_dir: Path, db_path: Path, "n_missing": len(missing)} +# --- reading it back, for the campaign-level merges ------------------------- +# The Snakefile loads this index into dicts at parse time and derives the +# campaign's unit sets from them (TILES_READY, and the exposures those tiles +# read). A merge JOB has to derive the same two sets, and cannot be handed them +# on its command line — ~20k paths is an order of magnitude over Linux's 128 KiB +# MAX_ARG_STRLEN for a single argv entry. So it is given the two things the +# Snakefile itself started from, the tile list and this database, and rebuilds +# the sets here. Both halves therefore read the schema through one module rather +# than two hand-written queries that could drift apart. + + +def campaign_tiles(tile_list: Path, db_path: Path) -> list[str]: + """The campaign's ready tiles: declared in the list AND indexed. + + Exactly the Snakefile's TILES_READY, computed the same way from the same two + files — a declared tile with no indexed exposure list cannot have been + computed, so it has no catalogue to merge. + """ + # DEDUPED, order preserved. The tile list is appended to by hand across a + # campaign, so a tile can appear twice; a merge would then try to write that + # tile's dataset twice and die on the second. Deduping here rather than at + # the call sites keeps the answer the same for every reader of the index. + seen, declared = set(), [] + with open(tile_list) as f: + for line in f: + tile = line.strip() + if tile and tile not in seen: + seen.add(tile) + declared.append(tile) + con = sqlite3.connect(db_path, timeout=60) + indexed = {r[0] for r in con.execute("SELECT DISTINCT tile_id FROM tile_exposures")} + con.close() + return [t for t in declared if t in indexed] + + +def campaign_exposures(tile_list: Path, db_path: Path) -> list[str]: + """Every exposure the campaign's ready tiles read, sorted. + + Exactly the set the Snakefile's persist_manifests() builds its manifest + paths from. + """ + tiles = set(campaign_tiles(tile_list, db_path)) + con = sqlite3.connect(db_path, timeout=60) + rows = con.execute("SELECT tile_id, exp_id FROM tile_exposures").fetchall() + con.close() + return sorted({e for t, e in rows if t in tiles}) + + def main() -> None: p = argparse.ArgumentParser(description=__doc__) p.add_argument("--tile-list", required=True, type=Path, diff --git a/workflow/scripts/hdf5_reconcile.py b/workflow/scripts/hdf5_reconcile.py new file mode 100644 index 000000000..c99101860 --- /dev/null +++ b/workflow/scripts/hdf5_reconcile.py @@ -0,0 +1,265 @@ +#!/usr/bin/env python3 +"""Bring an hdf5 catalogue into agreement with a campaign, one dataset per unit. + +Shared by the two campaign-level merges — ``merge_final_cat.py`` (one dataset +per tile) and ``merge_star_cat.py`` (one per exposure) — because they want +exactly the same thing of their output and disagreeing about it would be a bug +waiting to happen rather than a difference worth having. + +WHY RECONCILE RATHER THAN REBUILD. The output must be a function of the input +set — that is what makes the rules' fingerprints mean anything — but reading +every unit to add one is ~800 GB of IO at DR6 scale for a few tens of MB of new +data. So the file is brought INTO AGREEMENT with the campaign instead: + + * a unit with no dataset is read and added; + * a dataset whose unit has left the campaign is deleted; + * a dataset whose SOURCE has changed is re-read. Each records its source's + size and mtime as attributes, and a mismatch is what changed means. This is + the only reason a finished unit is read twice, and it is why the file + cannot drift from its inputs the way an append-only tool does; + * a dataset whose column set was written under a DIFFERENT SCHEMA is re-read. + The column set is the one input nothing else can see: it is not a source + file, so no stamp moves when it changes. It travels as a digest on the + file's root. + * a dataset that agrees with its source and its schema is left alone, unread. + +An append therefore READS exactly the appended units. It still WRITES the whole +file: the existing one is copied so the result can be moved into place +atomically, which costs one pass over it and, briefly, twice its size on disk. +That is the cheap half by orders of magnitude — copying a 1 GB hdf5 against +re-reading 800 GB of catalogues — but it is not free, and `apply` refuses rather +than filling the filesystem when the free space is not there. + +WHAT IS AND IS NOT A FUNCTION OF THE INPUT SET. The file's CONTENT is: the same +units with the same sources give the same datasets, the same columns and the +same count attribute, whether they arrived at once or one batch at a time. Its +BYTE LAYOUT is not, because hdf5 lays a group out in the order things were +added. That is the trade for not re-reading the campaign, and it is why the +no-op case compares ACTIONS rather than bytes. + +UNTOUCHED ON A NO-OP, which is stronger than byte-stable and cheaper to +establish. Reconciling is PLANNED against a read-only open; an empty plan never +opens the file for writing, so its mtime cannot move — and mtime is a rerun +trigger, so an unconditional rewrite would make every invocation look like a +change. +""" + +import hashlib +import json +import shutil +import sys +from pathlib import Path + +import h5py + + +def schema_digest(columns) -> str: + """A fingerprint of the COLUMN SET the datasets were written with.""" + return hashlib.md5("\n".join(columns).encode()).hexdigest()[:16] + + +def code_provenance(snapshot_json) -> dict: + """The launch code's identity, to stamp onto the merged file's root. + + ``snapshot_json`` is ``sp run``'s code snapshot (``bin/sp``'s + ``$STATE_DIR/code/snapshot.json``), passed through by the calling rule. A + workflow driven outside ``sp run`` has no such file — the merge still + succeeds, and the caller writes ``code_head = "unknown"`` rather than + failing an otherwise-good build. + """ + if not snapshot_json or not Path(snapshot_json).exists(): + return {"head": "unknown"} + data = json.loads(Path(snapshot_json).read_text()) + out = {k: data[k] for k in ("head", "branch", "dirty", "taken_at") + if k in data} + if data.get("dirty") and data.get("dirty_files"): + out["dirty_files"] = data["dirty_files"] + return out + + +def stamp(path: Path) -> tuple: + """A source's identity, as recorded on the dataset built from it. + + Size and mtime, not a checksum: the question is "did this change since we + read it", which mtime answers for a pipeline that writes a file once. A + campaign that rewrote a source in place with identical size and mtime would + defeat it, and nothing does. + """ + st = Path(path).stat() + return st.st_size, st.st_mtime_ns + + +class Plan: + """What reconciling requires: three unit lists. + + ``add`` and ``refresh`` are both "read the source and write the dataset"; + they are separate only so the log can say which happened, because a refresh + means a finished unit's source moved under us and that is worth seeing. + """ + + def __init__(self, add, refresh, remove): + self.add, self.refresh, self.remove = add, refresh, remove + + def empty(self): + return not (self.add or self.refresh or self.remove) + + def describe(self): + return (f"{len(self.add)} added, {len(self.refresh)} refreshed, " + f"{len(self.remove)} removed") + + +def plan(output: Path, group_path: str, units: list, digest: str) -> Plan: + """Compare the file on disk with the campaign, WITHOUT writing anything.""" + if not output.exists(): + return Plan([u for u, _ in units], [], []) + + want = {unit for unit, _ in units} + add, refresh = [], [] + with h5py.File(output, "r") as f: + stale_schema = f.attrs.get("param_digest") != digest + have = dict(f[group_path].items()) if group_path in f else {} + present = set(have) + for unit, source in units: + if unit not in present: + add.append(unit) + elif stale_schema: + refresh.append(unit) + else: + attrs = have[unit].attrs + if (int(attrs.get("src_bytes", -1)), + int(attrs.get("src_mtime_ns", -1))) != stamp(source): + refresh.append(unit) + return Plan(add, refresh, sorted(present - want)) + + +# Twice the file, plus a tenth of it again: the copy and the original coexist, +# and hdf5 is not a format to run to the last byte of a filesystem on. +FREE_SPACE_MARGIN = 2.1 + + +def check_free_space(output: Path) -> None: + """Refuse to start a rewrite the filesystem cannot hold. + + A merge that fills /project does not just fail: it fails everything else + writing there at the same time, and it can leave a truncated tmp beside a + catalogue people trust. Cheaper to say so first. + """ + if not output.exists(): + return + size = output.stat().st_size + free = shutil.disk_usage(output.parent).free + if free < size * FREE_SPACE_MARGIN: + sys.exit( + f"hdf5_reconcile: {output.parent} has {free / 1e9:.1f} GB free and " + f"this merge needs about {size * FREE_SPACE_MARGIN / 1e9:.1f} GB — " + f"it rewrites {output.name} ({size / 1e9:.1f} GB) through a tmp " + f"copy beside it. Free space or move products_dir; the existing " + f"catalogue is untouched.") + + +def check_sole_group(output: Path, group_path: str) -> None: + """One file, one campaign — refuse to half-update a file holding two. + + Renaming `campaign:` mid-flight points the rule at a NEW group inside the + SAME file (the path carries the campaign only on the tile side, where the + group does). Reconciling would then add a second group beside the first, + leave the first frozen and stale, and set a count attribute describing only + one of them. Nothing downstream reads such a file correctly, and no rule + here means to produce one. Say what is there and stop. + """ + if not output.exists() or "/" not in group_path: + return + parent, leaf = group_path.rsplit("/", 1) + with h5py.File(output, "r") as f: + if parent not in f: + return + others = sorted(k for k in f[parent] if k != leaf) + if others: + sys.exit( + f"hdf5_reconcile: {output} already holds {parent}/" + f"{', '.join(others)} beside {group_path}. One file is one " + f"campaign: reconciling would freeze the other group and count " + f"only this one. Point `campaign:` back, or write to a new path.") + + +def apply(output: Path, group_path: str, todo: Plan, units: list, read, + digest: str, count_attr: str, provenance: dict | None = None) -> None: + """Carry the plan out on a tmp file, then move it into place. + + ``read(unit, source)`` returns the structured array for one unit; it is + called only for the units the plan names, which is what makes an append + cheap. + + ``provenance`` (``code_provenance()``'s return) is stamped onto the file's + root as ``code_head``/``code_branch``/``code_dirty``/``code_snapshot_at``, + plus ``code_dirty_files`` (newline-joined) when the snapshot was dirty. It + is written here, alongside ``count_attr`` and ``param_digest``, rather than + on every no-op invocation: reconciling is planned against a read-only open, + and an empty plan must leave the file's mtime alone (see the module + docstring), so a run that changes no data never touches the file even if + the code that would have produced it has moved on. + + TWO WAYS TO BUILD THE TMP, and which one is used is about SPACE, not speed. + HDF5 never reclaims the space a deleted dataset occupied, so a file that is + copied and then edited in place grows for the life of the campaign — every + refresh of a unit leaks that unit. So: + + * a plan that only ADDS copies the existing file and appends to it. There + is nothing to reclaim, and copying beats rewriting. It is still a pass + over the whole file — an append is cheap in READS, not in writes. + * a plan that removes or refreshes anything builds the tmp FRESH, moving + the datasets it keeps across with h5py's own group copy — a + dataset-level copy inside the library that never reads a row into numpy + — and writing only the units that actually changed. The result is + compact. + + Either way the tmp is moved into place at the end, so a crash mid-merge + leaves the old catalogue intact rather than a half-written one. A SIGKILL + between writing the tmp and renaming it leaves the tmp behind — one file, + beside the catalogue, overwritten by the next run; the rename itself is + atomic, which is the property that matters. + """ + sources = dict(units) + rewrite = bool(todo.remove or todo.refresh) + check_free_space(output) + check_sole_group(output, group_path) + written = set(todo.add) | set(todo.refresh) + keep = [u for u, _ in units if u not in written] + tmp = output.with_name(output.name + ".tmp") + try: + tmp.unlink(missing_ok=True) + if output.exists() and not rewrite: + shutil.copy2(output, tmp) + with h5py.File(tmp, "a") as f: + group = (f[group_path] if group_path in f + else f.create_group(group_path)) + if rewrite and output.exists(): + with h5py.File(output, "r") as src: + for unit in keep: + # File.copy, not Dataset.copy — the latter does not + # exist, and the difference only shows when a plan both + # rewrites and keeps something. + src.copy(f"{group_path}/{unit}", group, name=unit) + for unit in todo.add + todo.refresh: + source = sources[unit] + data = read(unit, source) + dset = group.create_dataset(unit, data=data, dtype=data.dtype) + # The dataset's own record of what it was read from; this is + # what lets a later invocation leave it alone. + dset.attrs["src_bytes"], dset.attrs["src_mtime_ns"] = \ + stamp(source) + f.attrs[count_attr] = len(group) + f.attrs["param_digest"] = digest + if provenance: + f.attrs["code_head"] = provenance.get("head", "unknown") + for key, attr in (("branch", "code_branch"), + ("dirty", "code_dirty"), + ("taken_at", "code_snapshot_at")): + if key in provenance: + f.attrs[attr] = provenance[key] + if provenance.get("dirty_files"): + f.attrs["code_dirty_files"] = \ + "\n".join(provenance["dirty_files"]) + tmp.replace(output) # atomic: same filesystem + finally: + tmp.unlink(missing_ok=True) diff --git a/workflow/scripts/merge_final_cat.py b/workflow/scripts/merge_final_cat.py new file mode 100644 index 000000000..aaca2298c --- /dev/null +++ b/workflow/scripts/merge_final_cat.py @@ -0,0 +1,189 @@ +#!/usr/bin/env python3 +"""Collect the campaign's per-tile final catalogues into ONE hdf5 file. + +Run as the shell of the campaign-level ``final_cat_merge`` rule, never by hand. + +WHAT IT PRODUCES, AND FOR WHOM. ``/final_cat_.hdf5``: +one dataset per tile, carrying the columns named by +``workflow/config/cfis/final_cat.param``, plus an ``n_tiles`` attribute on the +file root. sp_validation opens that file as its ``galaxy_cat_path`` +(``sp_validation/catalog.py``), so its SCHEMA is an interface and not a choice — +see ``SPVAL_GROUP`` below for the one legacy literal in it. + +(sp_validation's own ``merge_catalogues`` is a different layer entirely: it +works over already-calibrated ``shape_catalog_comprehensive_*.fits``. It does +not do this merge, and this does not do that one.) + +WHAT IT REUSES, AND WHAT IT DOES NOT. The column extraction is +``create_final_cat.py``'s — ``read_param_file`` for the parameter list, +``read_data`` and ``copy_data`` for pulling those columns out of one catalogue +with their FITS dtypes — so the column grammar keeps exactly one definition. +Those three are REPRODUCIBLE FUNCTIONS, and this PR is what made them so: the +parameter list comes back ordered rather than through a set, ``copy_data`` +allocates the requested columns alone rather than leaving every other column of +the source as uninitialised memory, and a missing column raises with its own +name instead of falling out of a bare ``except:`` as an UnboundLocalError. The +fixes are upstream, in that script, because a hand-run of it deserves them as +much as this rule does. +Its ``process()`` is NOT used and neither is any of its discovery: that function +walks a directory tree the workflow does not have and never will, and it groups +by a unit ShapePipe v2 no longer has. This script walks the workflow's own +products tree instead (``tiles/<2-char prefix>//final_cat-.fits``) and +writes the hdf5 itself. + +WHERE ``create_final_cat.py`` IS FOUND. Beside this workflow, at +``/scripts/python/create_final_cat.py`` — resolved relative to THIS file, +so it follows the launch code snapshot (``bin/sp``) exactly as +``workflow/scripts/*`` does, and a campaign never reads a mid-run edit. It is +loaded by path rather than imported: it is a script, not an installed module, +and the container's ``shapepipe`` install does not carry it. + +IT RECONCILES, IT NEITHER REBUILDS NOR BLINDLY APPENDS, and the machinery for +that is ``hdf5_reconcile.py``, shared with the star side so the campaign's two +products cannot disagree about what an output owes its inputs. That module +carries the argument in full: an append reads the appended tiles, a source that +changed is re-read, a tile that left the campaign is deleted, a column-set +change refreshes everything, and a no-op leaves the file untouched. +``create_final_cat.py``'s own ``process()`` implements only the append-only half +— it skips a tile already in the file, whatever the file on disk now says — +which is right for a hand-driven update and wrong for a DAG output. (Its ``-s`` +single-ID mode implements ``check`` and ``remove``; ``add`` is accepted by the +argument validator and then falls through to the ordinary walk, so it is not a +way to add one tile by hand.) + +WHICH TILES — AND WHY THE JOB DERIVES THE SET RATHER THAN BEING TOLD IT. The set +is the CAMPAIGN's: every tile both declared in ``tile_list`` and present in the +index, which is exactly the Snakefile's TILES_READY, rebuilt here from the same +two files the Snakefile started from (``--tile-list`` and ``--index-db``, read +through ``build_index.campaign_tiles`` so there is one definition and not two +that can drift). It is derived rather than passed because at DR6 scale the set +is ~20k paths and a shell command reaches ``execve`` as a SINGLE argv entry +capped at 128 KiB by ``MAX_ARG_STRLEN``; the rule's ``input`` is the DAG edge +and its ``params`` carries a fingerprint of that same list, which is the rerun +trigger. + +THE TWO SETS ARE THE SAME SET, which is the point of deriving it this way rather +than globbing ``/tiles``: a products root shared with an earlier, +larger tile list would hand the job tiles the fingerprint never saw and no rerun +trigger would notice. A tile in the derived set whose catalogue is missing is a +hard error here, not a skip — under the DAG it cannot happen, since every one of +them is a declared input of this job. +""" + +import argparse +import importlib.util +import sys +from pathlib import Path + +# Same directory; the rule invokes this file by path, so it is sys.path[0]. +import build_index +import hdf5_reconcile + +# /scripts/python/create_final_cat.py, from /workflow/scripts/this. +CFC_PATH = (Path(__file__).resolve().parents[2] + / "scripts" / "python" / "create_final_cat.py") + + +def spval_group(campaign: str) -> str: + """The hdf5 group the campaign's per-tile datasets live under. + + ``patches/`` is a LEGACY KEY IN sp_validation's FILE SCHEMA, kept verbatim + only so its reader works unchanged (CosmoStat/sp_validation#340 tracks + removing it); it names nothing in this workflow, which has campaigns and + tiles and no other unit. This is the one place the literal appears — + everything else here says campaign. + """ + return f"patches/{campaign}" + + +def load_create_final_cat(): + """The hdf5 layout's definition, loaded by path (see the module docstring).""" + if not CFC_PATH.exists(): + sys.exit(f"merge_final_cat: {CFC_PATH} is not there — the launch code " + f"snapshot must carry scripts/python/ (see bin/sp).") + spec = importlib.util.spec_from_file_location("create_final_cat", CFC_PATH) + mod = importlib.util.module_from_spec(spec) + spec.loader.exec_module(mod) + return mod + + +def catalogues(products_dir: Path, tile_list: Path, index_db: Path) -> list: + """``(tile ID, path)`` for the campaign's tiles, in ID order. + + Not a glob over the products root: see the module docstring on why the set + is the campaign's and not the filesystem's. + """ + out, missing = [], [] + for tile in sorted(build_index.campaign_tiles(tile_list, index_db)): + path = (products_dir / "tiles" / tile[:2] / tile + / f"final_cat-{tile}.fits") + if path.exists(): + out.append((tile, path)) + else: + missing.append(tile) + if missing: + sys.exit(f"merge_final_cat: {len(missing)} campaign tile(s) have no " + f"final catalogue: {' '.join(missing[:5])}" + f"{' ...' if len(missing) > 5 else ''}") + return out + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--products-dir", required=True, type=Path, + help="the persistent root; per-tile catalogues are found " + "beneath it") + p.add_argument("--tile-list", required=True, type=Path, + help="the campaign's tile list (config tile_list)") + p.add_argument("--index-db", required=True, type=Path, + help="the campaign's run index (config outputs.index_db)") + p.add_argument("--output", required=True, type=Path) + p.add_argument("--campaign", required=True, + help="names the campaign's group in the output file") + p.add_argument("--param-file", required=True, type=Path, + help="workflow/config/cfis/final_cat.param — the column list") + p.add_argument("--hdu", type=int, default=1) + p.add_argument("--snapshot-json", type=Path, default=None, + help="sp run's code snapshot (bin/sp's " + "$STATE_DIR/code/snapshot.json); absent outside sp run") + args = p.parse_args() + + cfc = load_create_final_cat() + param_list = cfc.read_param_file(str(args.param_file), verbose=False) + if not param_list: + sys.exit(f"merge_final_cat: no columns read from {args.param_file}") + # read_data/copy_data read their knobs out of this dict, exactly as + # create_final_cat.py's own main() builds it. + params = {"hdu_num": args.hdu, "param_list": param_list, "verbose": False} + + tiles = catalogues(args.products_dir, args.tile_list, args.index_db) + if not tiles: + # An empty hdf5 would satisfy every downstream existence check and + # produce an empty shear catalogue. + sys.exit(f"merge_final_cat: no tile in {args.tile_list} is indexed in " + f"{args.index_db}, so there is nothing to merge") + + args.output.parent.mkdir(parents=True, exist_ok=True) + group_path = spval_group(args.campaign) + digest = hdf5_reconcile.schema_digest(param_list) + + def read_tile(tile, path): + """One tile's requested columns, via create_final_cat.py's own reader.""" + extracted, dtype = cfc.read_data(str(path), params) + return cfc.copy_data(params["param_list"], extracted, dtype) + + todo = hdf5_reconcile.plan(args.output, group_path, tiles, digest) + if todo.empty(): + print(f"[merge_final_cat] unchanged: {args.output} " + f"({len(tiles)} tile(s))") + return + hdf5_reconcile.apply(args.output, group_path, todo, tiles, read_tile, + digest, "n_tiles", + hdf5_reconcile.code_provenance(args.snapshot_json)) + print(f"[merge_final_cat] {todo.describe()} -> {args.output} " + f"({len(tiles)} tile(s), {len(param_list)} column(s), " + f"group {group_path})") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/merge_star_cat.py b/workflow/scripts/merge_star_cat.py new file mode 100644 index 000000000..c9b4e45f1 --- /dev/null +++ b/workflow/scripts/merge_star_cat.py @@ -0,0 +1,329 @@ +#!/usr/bin/env python3 +"""Collect the campaign's per-CCD PSF validation catalogues into ONE hdf5 file. + +Run as the shell of the campaign-level ``star_cat_merge`` rule, never by hand. + +WHAT IT PRODUCES, AND FOR WHOM. ``/full_starcat_.hdf5``: +one dataset per exposure at ``exposures/``, holding that exposure's every +CCD's ``validation_psf--.fits`` rows stacked, with a ``CCD_NB`` column +recording which CCD each row came from. It is the input to the rho/tau +statistics. Historically this was ``combine_runs.bash psf`` plus a +``merge_starcat_runner`` pass producing one flat FITS table, +``full_starcat-0000000.fits``, and sp_validation still opens that name today; +its readers move to this hdf5 under CosmoStat/sp_validation#340, the same +migration that retires the ``patches/`` key on the galaxy side. + +WHY HDF5, AND WHY ONE DATASET PER EXPOSURE. The campaign's two products should +behave the same way, and one flat table cannot: appending a tile meant +restacking every exposure the campaign had ever seen — ~40 GB of members at DR6 +scale to add ~2 MB. Per-exposure datasets make the file RECONCILABLE +(hdf5_reconcile.py carries that argument, and merge_final_cat.py is the same +machinery on the tile side), so an append reads the appended exposures and +nothing else while the file still cannot drift from its inputs. Memory follows: +one exposure at a time, not one campaign. + +NATIVE DTYPES. Columns are written as the validation_psf files store them — +float32 stays float32. The FITS writer this replaces widened every float column +to ``1D``, doubling both the file and the peak memory of the job that wrote it, +for no information. + +CCD_NB IS AN INTEGER. It is parsed out of the member name +(``validation_psf--.fits``), where it is always digits, so a string +buys nothing — and an int column costs 4 bytes a row against the 8 a +two-character fixed-width string does. + +IT READS THE TARS, IT DOES NOT UNPACK THEM. ``exp_persist`` packs each +exposure's keepers into one uncompressed tar on the persistent root +(``/exp///psf/.tar``) precisely because inodes, +not bytes, bind on /project. Unpacking ~20k tars x ~40 members to merge them +would materialise ~800k files on the filesystem that design exists to protect. +Members are read through the archive's own file object — seekable, the tar +being uncompressed by design — so the counting pass costs a header rather than +a member. + +THE OPTIONAL COLUMNS ARE A PER-FILE QUESTION. A pix2wcs-converted catalogue has +no MAG/SNR/ACCEPTED where an ordinary one does, and a campaign can hold both. +Deciding once for the merge is wrong in both directions: it either fails on the +first converted file or silently zeroes the real values of every ordinary one. +Each file is asked for its own schema, and only the files that lack a column are +zero-filled. + +WHICH EXPOSURES — AND WHY THE JOB DERIVES THE SET RATHER THAN BEING TOLD IT. +The set is the CAMPAIGN's: every exposure read by a tile that is both declared +in ``tile_list`` and present in the index, which is the Snakefile's TILES_READY +walked one edge further. This script rebuilds it from the same two files the +Snakefile started from (``--tile-list`` and ``--index-db``, both small, both on +the persistent root, both read through ``build_index.campaign_exposures`` so +there is one query and not two that can drift), and then takes the exposures +whose ``exp_persist`` manifest is on the persistent root. + +It is derived rather than passed because at DR6 scale the set is ~20k paths and +a shell command reaches ``execve`` as a SINGLE argv entry capped at 128 KiB by +``MAX_ARG_STRLEN``. So the rule's ``input`` is the DAG EDGE — what must exist +before this runs — and its ``params`` carries a FINGERPRINT of the same set, +which is what makes the merge rerun when the set changes. A glob over +``/exp`` would NOT be the same set: it would sweep in exposures of +an earlier, larger tile list sharing the products root, stacking rows the +fingerprint never saw and no rerun trigger would notice. + +THE MANIFEST, NOT THE TAR, IS WHAT IT READS FIRST: the manifest records what was +actually packed, member by member, with sizes and the product each came from, so +this script never guesses at tar contents. +""" + +import argparse +import json +import sys +import tarfile +from fnmatch import fnmatch +from pathlib import Path + +import numpy as np +from astropy.io import fits + +# Same directory; the rule invokes this file by path, so it is sys.path[0]. +import build_index +import hdf5_reconcile +import persist_exp + +# The members this merge consumes, named as the keep list names them and +# resolved through the same catalogue persist_exp packs by — so the glob has one +# definition and adding a product cannot leave the two disagreeing. They are +# always there to find: persist_exp packs this product for every exposure +# whatever `persist_exp:` says, and fails the pack rather than writing a +# manifest without it. +MEMBER_PRODUCT = persist_exp.ALWAYS +MEMBER_PATTERN = persist_exp.resolve(MEMBER_PRODUCT) + +# The group holding the per-exposure datasets. Unlike the galaxy side's +# `patches/`, this name is ours and says what it holds. +GROUP = "exposures" + +# The validation_psf table's HDU: what MergeStarCatPSFEX defaulted to and what +# psfex_interp writes — a SExtractor-style file, empty primary, header-carrying +# image extension, then the table. +HDU = 2 + +# The columns, in the order the FITS full_starcat carried them, which is the +# order every consumer has seen. The optional three are zero-filled per file. +COLUMNS = ("X", "Y", "RA", "DEC", + "HSM_G1_PSF", "HSM_G2_PSF", "HSM_T_PSF", + "HSM_G1_STAR", "HSM_G2_STAR", "HSM_T_STAR", + "HSM_FLAG_PSF", "HSM_FLAG_STAR") +# CANONICAL DTYPES, not whatever the first file that carries the column happens +# to use. These three are absent from pix2wcs-converted catalogues, so an +# exposure whose files all lack them would otherwise be allocated a fallback +# dtype while its neighbours got the real one — and datasets under exposures/* +# would then differ in dtype, which np.concatenate refuses and no digest can +# repair, since nothing about the schema CHANGED. Pinning the dtype here is what +# makes every exposure's dataset the same shape whatever its files carry. +OPTIONAL = {"MAG": np.float32, "SNR": np.float32, "ACCEPTED": np.int32} +CCD_COLUMN = "CCD_NB" +ALL_COLUMNS = COLUMNS + tuple(OPTIONAL) + (CCD_COLUMN,) + + +def ccd_number(member_name: str) -> int: + """The CCD this member's rows belong to: ``validation_psf--.fits``. + + Always digits, which is why the column is an int; a member name that does + not carry one is a tar we do not understand, and saying so beats writing a + sentinel into the catalogue. + """ + ccd = member_name.rsplit(".", 1)[0].rsplit("-", 1)[-1] + if not ccd.isdigit(): + sys.exit(f"merge_star_cat: cannot read a CCD number out of member " + f"name {member_name!r}") + return int(ccd) + + +def is_member(entry: dict) -> bool: + """Is this manifest entry one of the members this merge reads? + + BY PRODUCT NAME, OR FAILING THAT BY FILE NAME. persist_exp records the + product every member came from and always packs psf_validation, so the name + is the answer for anything it writes today. The glob is the fallback, and it + earns its place twice over: a tar packed before the product field existed + has no label at all, and a keep list written as a raw glob + (`validation_psf-*.fits` rather than `psf_validation`) labels its members + with the glob. Neither should make the campaign's star catalogue silently + empty. + """ + return (entry.get("product") == MEMBER_PRODUCT + or fnmatch(entry["name"], MEMBER_PATTERN)) + + +def manifests(products_dir: Path, tile_list: Path, index_db: Path) -> list: + """``(exposure, manifest path)`` for the campaign's packed exposures. + + Not a glob over the products root: see the module docstring on why the set + is the campaign's and not the filesystem's. + """ + out = [] + for exp in build_index.campaign_exposures(tile_list, index_db): + path = (products_dir / "exp" / exp[:2] / exp / "manifests" + / "exp_persist.json") + if path.exists(): + out.append((exp, path)) + return out + + +def tars(manifest_paths: list) -> tuple: + """``[(exposure, tar path)]`` for the merge, and the exposures with nothing. + + Every tar is checked for existence HERE, so a products root missing a file + fails before a single row is read rather than an hour in. The tar is also + the unit's SOURCE for reconciling: its size and mtime are what a later + invocation compares against to decide whether this exposure changed. + """ + chosen, empty = [], [] + for exp, man_path in manifest_paths: + man = json.loads(man_path.read_text()) + # MEMBERSHIP IS THE MEMBER NAME, and only the member name. is_member() + # will also accept a manifest's own product LABEL, which is the right + # test for "did this exposure keep the product" — but a label is not + # what read_exposure() selects on, and a mislabeled entry whose name + # does not match would put this exposure in the merge and then abort + # the whole campaign when the tar turned out to hold nothing selectable. + # So the two agree by construction: both ask the name. + if not any(fnmatch(f["name"], MEMBER_PATTERN) for f in man["files"]): + if any(is_member(f) for f in man["files"]): + # Labelled as the product, named as something else. Worth one + # line — it means a manifest we did not write, or a keep list + # whose glob does not match the member it matched. + print(f"[merge_star_cat] {exp}: manifest labels a " + f"{MEMBER_PRODUCT} member whose name does not match " + f"{MEMBER_PATTERN}; not merging it") + empty.append(exp) + continue + tar_path = Path(man["tar"]) + if not tar_path.exists(): + sys.exit(f"merge_star_cat: {man_path} names a tar that is not " + f"there: {tar_path}") + chosen.append((exp, tar_path)) + return chosen, empty + + +def read_exposure(exp: str, tar_path: Path) -> np.ndarray: + """One exposure's every CCD, stacked, as a structured array. + + TWO PASSES over the tar's members, and neither holds the exposure twice: + the first reads only each member's FITS HEADER — NAXIS2, the row count — + and the second allocates the columns once at their exact final length and + fills them slice by slice. Members are visited in sorted name order, so the + row order is a function of the tar's contents alone. + + NOTE ON WHEN THIS IS CALLED AGAIN. The unit's source is the TAR, so adding a + retention product re-packs it, moves its mtime, and refreshes this exposure + even though its validation members are byte-for-byte what they were. Reading + one exposure is seconds and the alternative — stamping the members rather + than the archive — buys a rarely-taken shortcut for a per-member bookkeeping + cost on every exposure. Not worth it. + """ + try: + tf = tarfile.open(tar_path) + except tarfile.TarError as exc: + sys.exit(f"merge_star_cat: cannot read {tar_path}: {exc}. That tar is " + f"this exposure's only copy of its PSF products — do not " + f"delete it; re-pack the exposure if its scratch store is " + f"still there, and treat the exposure as lost if it is not.") + with tf: + names = sorted(n for n in tf.getnames() + if fnmatch(n, MEMBER_PATTERN)) + if not names: + sys.exit(f"merge_star_cat: {tar_path} holds no {MEMBER_PATTERN}") + + # --- pass 1: row counts and dtypes, from headers alone -------------- + counts, dtypes, n_total = [], None, 0 + for name in names: + with fits.open(tf.extractfile(name), memmap=False, + ignore_missing_simple=True) as hdul: + hdu = hdul[HDU] + counts.append(hdu.header["NAXIS2"]) + # ColDefs.dtype describes the table without reading it. NOTE: + # it is the RAW storage dtype and ignores TSCAL/TZERO, so a + # scaled column would be allocated narrower than the values + # .data returns. Latent, not live: no validation_psf column is + # scaled. Read the dtype off .data if one ever is. + if dtypes is None: + dtypes = hdu.columns.dtype + n_total += counts[-1] + + fields = [(c, dtypes[c]) for c in COLUMNS] + # The optional three take their CANONICAL dtype, not one file's (see + # OPTIONAL): every exposure's dataset must have the same dtype whether + # or not its files carry the column. + fields += list(OPTIONAL.items()) + fields += [(CCD_COLUMN, np.int32)] + data = np.empty(n_total, dtype=np.dtype(fields)) + + # --- pass 2: fill --------------------------------------------------- + at = 0 + for name, n_rows in zip(names, counts): + with fits.open(tf.extractfile(name), memmap=False, + ignore_missing_simple=True) as hdul: + rows = hdul[HDU].data + have = set(rows.dtype.names or ()) + sl = slice(at, at + n_rows) + for col in COLUMNS: + data[col][sl] = rows[col] + for col in OPTIONAL: + # THIS file's schema, not the exposure's. + data[col][sl] = rows[col] if col in have else 0 + data[CCD_COLUMN][sl] = ccd_number(name) + at += n_rows + + if at != n_total: + raise ValueError(f"merge_star_cat: {tar_path}: pass 1 counted " + f"{n_total} rows, pass 2 filled {at}") + return data + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--products-dir", required=True, type=Path, + help="the persistent root; exp_persist manifests and tars " + "are found beneath it") + p.add_argument("--tile-list", required=True, type=Path, + help="the campaign's tile list (config tile_list)") + p.add_argument("--index-db", required=True, type=Path, + help="the campaign's run index (config outputs.index_db)") + p.add_argument("--output", required=True, type=Path) + p.add_argument("--campaign", required=True, + help="named in the log; the group name is fixed") + p.add_argument("--snapshot-json", type=Path, default=None, + help="sp run's code snapshot (bin/sp's " + "$STATE_DIR/code/snapshot.json); absent outside sp run") + args = p.parse_args() + + manifest_paths = manifests(args.products_dir, args.tile_list, args.index_db) + chosen, empty = tars(manifest_paths) + if not chosen: + # Not a no-op: an empty star catalogue would pass every downstream + # existence check and produce meaningless rho statistics. + sys.exit(f"merge_star_cat: no {MEMBER_PRODUCT} member in any of " + f"{len(manifest_paths)} exp_persist manifest(s) for this " + f"campaign. persist_exp packs {MEMBER_PRODUCT} for every " + f"exposure, so this means the manifests are not what we think " + f"they are.") + if empty: + print(f"[merge_star_cat] {len(empty)} exposure(s) persisted no " + f"{MEMBER_PRODUCT}: {', '.join(sorted(empty)[:5])}" + f"{' ...' if len(empty) > 5 else ''}") + + args.output.parent.mkdir(parents=True, exist_ok=True) + digest = hdf5_reconcile.schema_digest(ALL_COLUMNS) + todo = hdf5_reconcile.plan(args.output, GROUP, chosen, digest) + if todo.empty(): + print(f"[merge_star_cat] unchanged: {args.output} " + f"({len(chosen)} exposure(s))") + return + hdf5_reconcile.apply(args.output, GROUP, todo, chosen, read_exposure, + digest, "n_exposures", + hdf5_reconcile.code_provenance(args.snapshot_json)) + print(f"[merge_star_cat] {todo.describe()} -> {args.output} " + f"({len(chosen)} exposure(s), {len(ALL_COLUMNS)} column(s), " + f"campaign {args.campaign})") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/persist_exp.py b/workflow/scripts/persist_exp.py new file mode 100644 index 000000000..c3e14633d --- /dev/null +++ b/workflow/scripts/persist_exp.py @@ -0,0 +1,445 @@ +#!/usr/bin/env python3 +"""Pack ONE exposure's keepable PSF products into a tar off scratch, and record what went. + +Run as the shell of the in-DAG ``exp_persist`` rule, never by hand. + +WHY A COPY AND NOT AN EXEMPTION FROM CLEANUP. The obvious alternative — teach +``clean_exposure`` to spare these files — does not work, because reclamation is +not what threatens them. The exposure store lives on ``run_dir``, which is +/scratch: a 60-day purge takes everything there whether or not this workflow +ever cleaned it. ``products_dir`` is /project, backed up and not purged. So the +only way a per-exposure product outlives its campaign is to LEAVE THE +FILESYSTEM, and that is a copy. Reclamation ordering then falls out for free: +``clean_exposure`` takes this rule's manifest as an input, so the store is never +deleted before its keepers have been written elsewhere. + +WHY A SEPARATE RULE AND NOT A ``cp`` APPENDED TO ``exp_psf``. The list of what +to keep is a decision that will be revisited — rho statistics want one file +today, a residual study may want three tomorrow — and ``exp_psf`` is four hours +per exposure. The list rides on this rule's ``params``, so editing it makes +snakemake rerun THIS rule (seconds of cp) and leaves the PSF chain alone. Folded +into ``exp_psf``, the same edit would re-derive every PSF model in the campaign. + +WHAT IT SEARCHES. ``/output/run_sp_exp_SxSePsf/*/output/`` — the four +module output dirs of the PSF config (sextractor, setools, psfex, psfex_interp) +— RECURSIVELY. The recursion is not laziness: setools does not write flat, it +writes into ``mask/``, ``rand_split/``, ``new_cat/``, ``plot/`` and ``stat/`` +beneath its own output dir, so a caller who wrote ``star_split_ratio_80-*.fits`` +meaning "the training star sample" would match nothing under a non-recursive +glob. Patterns are therefore plain FILE names and the layout is ours to know, +not the config author's. + +RETENTION IS ADDITIVE, AND THAT IS A SAFETY PROPERTY. The keep list rides on +the rule's ``params``, so SHRINKING it reruns this script — and a naive rerun +would rewrite the tar without the products that were dropped, deleting them +from the backed-up filesystem because someone edited a config, with the scratch +store they came from usually long gone. An existing tar is therefore a FLOOR: +its members are carried into the new one whatever the current list says, and a +config change can only ever add. Removing a product is a deliberate act on +products_dir, not a config edit. + +THE KEEP LIST IS WHAT THE CAMPAIGN KEEPS ON TOP OF THE MERGE'S INPUTS. +``psf_validation`` is packed unconditionally (see ALWAYS below); ``persist_exp:`` +is purely optional retention, and an EMPTY one is a coherent instruction — the +tar then holds the star catalogue's inputs and nothing else. + +ZERO MATCHES FOR ONE PATTERN IS A WARNING, NOT A FAILURE. setools rejects sparse +CCDs (~0.2% attrition, tolerated by exp_psf's own count floor), so per-CCD +counts are not fixed, and a pattern naming an optional diagnostic may legitimately +find nothing. ZERO FILES IN TOTAL IS A FAILURE: it means the store was not what +we think it is, and writing a green manifest over that would let +``clean_exposure`` delete an exposure whose products were never saved. + +The manifest lists every member (name, pattern, source path, bytes), so a reader +knows what the tar holds without opening it. + +ONE UNCOMPRESSED TAR PER EXPOSURE, ``/.tar``, NOT LOOSE COPIES. +Inodes, not bytes, are what bind on /project: the group quota is ~1 M files, +and loose per-CCD copies are ~200 per exposure with all candidates on — ~25k for +a 64-tile campaign, ~2 M at DR6 scale, against ~7 GB of bytes. A tar collapses +that to one inode per exposure and costs nothing to read: FITS members go +``tarfile.open(t).extractfile(m).read()`` -> ``fits.open(io.BytesIO(...))``, +which is why a tar rather than a multi-HDU FITS bundle (the keep list mixes +FITS, ``.psf`` and ``.txt``; a FITS container could not hold the last two). +Uncompressed because FITS barely compresses and a plain tar is seekable. + +Members are FLAT — file name only, no module subtree — because the module a +file came from is already in its name and the consumer globs member names. A +name collision between two modules is therefore a hard error rather than a +silent overwrite; nothing in the current config can produce one, and if a +future one can we want to hear about it. + +The tar is written DETERMINISTICALLY (ownership zeroed, members in sorted +order, source mtimes kept), tmp-then-``cmp``-then-``mv``: a rerun over an +unchanged store produces a byte-identical tar and leaves the existing one's +mtime alone. + +The manifest is the rule's ONLY declared output, and it lives on the persistent +root beside the tar (``/exp///manifests/``, beside the tar's ``psf/``), NOT in +the exposure's scratch ``manifests/`` dir which ``clean_exposure`` deletes +wholesale. It is deliberately NOT a ``directory()`` output: what was copied, and +how big each file was, is provenance we want written down, and a directory +output attests only that some directory exists. + +It carries no timestamp and is written tmp-then-``cmp``-then-``mv`` (the pattern +``clean_exposure`` uses), so a rerun that packs the same files leaves the mtime +alone — mtime is a rerun trigger, and an unconditional rewrite would make every +downstream ``clean_exposure`` look out of date once per invocation. +""" + +import argparse +import filecmp +import json +import sys +import tarfile +from pathlib import Path + +# The PSF chain's run dir: RUN_NAME in config_exp_psfex.ini AND in +# config_exp_mccd.ini, which carry the same name on purpose so nothing +# downstream of exp_psf branches on the PSF model. Hardcoded rather than passed: +# this rule persists the PSF stage's products and nothing else, and a knob here +# would be a knob for "persist some other stage", which is a different rule. +# tests/unit/test_workflow_run_names.py holds this equal to the configs. +RUN_NAME = "run_sp_exp_SxSePsf" + +# --- the product catalogue (CosmoStat/shapepipe#844) ------------------------ +# THE SINGLE SOURCE OF TRUTH for what an exposure can keep. `persist_exp:` in +# config.yaml names PRODUCTS, not globs: `psf_model`, not `*.psf`. The glob is +# an implementation detail of the module that writes the file, and a keep list +# written in globs is a keep list nobody can read — the argument that produced +# #844 and the 2026-09-08 call's request to keep the PSF model, which had to be +# spelled `*.psf` to be said at all. +# +# Each entry is (glob, per-exposure size, what keeping it buys). Sizes are for +# 40 CCDs, measured on smk-m2 (127 exposures, 64 tiles); "?" means not yet +# measured. `persist_exp.py --list-products` renders this table, and +# config.yaml's block is that rendering rather than a second copy of it. +# +# ORDER IS THE ORDER OF THE CHAIN — sextractor, setools, psfex, psfex_interp — +# so the table reads as the pipeline runs. +# THE STAR CATALOGUE'S INPUTS ARE NOT A USER CHOICE. star_cat_merge stacks +# every CCD's psf_validation into the campaign's full_starcat, so exp_persist +# ALWAYS packs it, whatever `persist_exp:` says. Two reasons, and neither is +# about taste. It is the merged catalogue's PROVENANCE: a full_starcat with no +# per-exposure inputs beside it cannot be audited, re-cut or recomputed after a +# purge. And it is what keeps APPENDING TILES CHEAP: a tile added next month +# brings exposures whose validation catalogues must join the existing stack, and +# if the earlier ones are gone the merge either shrinks or rebuilds their chains +# from VOS. ~2 MB per exposure, so ~40 GB and ~40k inodes at DR6 scale, against +# a group quota of ~1 M inodes — the cost of being able to say where the number +# came from. +ALWAYS = "psf_validation" + +PRODUCTS = { + "star_selection": ( + "star_selection-*.fits", 24_500_000, + "setools' PRE-SPLIT selection. The only file that answers which stars " + "the selection cuts rejected and why; the split samples have already " + "lost the rejects."), + "star_train": ( + "star_split_ratio_80-*.fits", 19_900_000, + "the 80% TRAINING sample, the stars PSFEx actually fitted. Rows " + "duplicate star_selection."), + "star_test": ( + "star_split_ratio_20-*.fits", 7_100_000, + "the 20% VALIDATION sample — the positions psf_validation's rows " + "correspond to. Rows duplicate star_selection."), + "star_stats": ( + "star_stat-*.txt", None, + "setools' per-CCD STAT block: star counts, stars/deg^2, FWHM mode and " + "cuts. The selection's summary without its catalogue."), + "psf_model": ( + "*.psf", 2_800_000, + "the PSFEx model itself. Keeping it means the PSF can be " + "re-interpolated at ANY position later without rebuilding the exposure " + "chain — the single most capability-adding entry here."), + "psfex_cat": ( + "psfex_cat-*.cat", None, + "PSFEx's own output catalogue (FITS_LDAC): the per-star FLAGS_PSF and " + "CHI2_PSF, i.e. WHICH stars outlier rejection clipped. Not recoverable " + "from anything else — the .psf header keeps only the LOADED/ACCEPTED " + "counts."), + "psf_validation": ( + "validation_psf-*.fits", 2_000_000, + "the psfex_interp validation catalogue, one per CCD: the input to the " + "rho/tau statistics, and to the star_cat_merge rule that stacks them " + "into the campaign's full_starcat."), +} + +# PSFEx residual/check images and its XML diagnostics are deliberately absent: +# the committed default.psfex sets CHECKIMAGE_TYPE NONE and WRITE_XML N, so +# nothing is emitted to match. They are a config change first, a catalogue +# entry second. + +# A raw glob is still accepted, as an escape hatch for a file the catalogue does +# not name yet. The test is syntactic and deliberately cheap: a product name is +# a bare identifier, so anything carrying a glob metacharacter or a dot is a +# glob. That makes `*.psf`, `star_stat-*.txt` and `default.psfex` globs, and +# `psf_model` a name, with no ambiguity a user could stumble into. +_GLOBBY = set("*?[]. ") + + +def is_glob(entry: str) -> bool: + """True when this keep-list entry is a raw glob rather than a product name.""" + return any(ch in _GLOBBY for ch in entry) + + +def resolve(entry: str) -> str: + """The file-name glob for one keep-list entry, name or raw glob.""" + if is_glob(entry): + return entry + try: + return PRODUCTS[entry][0] + except KeyError: + raise KeyError( + f"unknown persist_exp product {entry!r}; the products are " + f"{', '.join(PRODUCTS)} (or write a raw glob such as '*.psf')" + ) from None + + +def product_of(entry: str) -> str: + """The NAME to record for an entry — the entry itself for a raw glob.""" + return entry + + +def render_products() -> str: + """The catalogue as a table, for --list-products and for config.yaml.""" + width = max(len(n) for n in PRODUCTS) + lines = [f"{'product'.ljust(width)} {'glob'.ljust(26)} size/exposure", + f"{'-' * width} {'-' * 26} -------------"] + for name, (glob, size, why) in PRODUCTS.items(): + size_s = "unmeasured" if size is None else f"{size / 1e6:.1f} MB" + lines.append(f"{name.ljust(width)} {glob.ljust(26)} {size_s}") + for i, chunk in enumerate(_wrap(why, 66)): + lines.append(f"{' ' * width} {chunk}") + return "\n".join(lines) + + +def _wrap(text: str, width: int) -> list: + out, line = [], "" + for word in text.split(): + if line and len(line) + 1 + len(word) > width: + out.append(line) + line = word + else: + line = f"{line} {word}".strip() + if line: + out.append(line) + return out + + +def collect(exp_dir: Path, patterns: list) -> tuple: + """Matched files per ENTRY, in a stable order, plus the entries that matched + nothing. Entries are product names or raw globs; resolve() takes either.""" + root = exp_dir / "output" / RUN_NAME + found, empty = {}, [] + for entry in patterns: + pat = resolve(entry) + # One glob per module output dir, recursive beneath it (see the module + # docstring on setools' subdirectories). sorted() over the union keeps + # the manifest byte-stable across filesystem readdir order. + hits = sorted({p for mod in sorted(root.glob("*/output")) + for p in mod.rglob(pat) if p.is_file()}) + if hits: + found[entry] = hits + else: + empty.append(entry) + return found, empty + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--exp-dir", type=Path, + help="the exposure's scratch store") + p.add_argument("--exp") + p.add_argument("--dest", type=Path, + help="/exp///psf; the tar is " + "/.tar") + p.add_argument("--manifest", type=Path) + p.add_argument("--pattern", action="append", default=[], + help=f"repeatable; a product name (see --list-products) or " + f"a raw file-name glob. {ALWAYS} is packed whether or " + f"not it is named — star_cat_merge needs it") + p.add_argument("--list-products", action="store_true", + help="print the product catalogue and exit") + args = p.parse_args() + + # --list-products is a QUERY, not a run: it answers "what can I keep?" and + # needs no exposure, so the run arguments are optional at the parser and + # required here instead. + if args.list_products: + print(render_products()) + return + missing = [f"--{n.replace('_', '-')}" for n in + ("exp_dir", "exp", "dest", "manifest") + if getattr(args, n) is None] + if missing: + p.error(f"the following arguments are required: {', '.join(missing)}") + + # The merge's input first and always, then whatever the campaign chose to + # keep on top of it (see ALWAYS). Deduped, so naming it explicitly in + # persist_exp: is harmless rather than a repeated pattern. + entries = [ALWAYS] + [e for e in args.pattern if e != ALWAYS] + + for entry in entries: # loud, and before any work + try: + resolve(entry) + except KeyError as exc: + sys.exit(f"persist_exp: {exc.args[0]}") + + found, empty = collect(args.exp_dir, entries) + # THE MERGE'S INPUTS ARE NOT ALLOWED TO BE MISSING, and this is a harder + # rule than "something matched". An exposure whose psfex_interp failed but + # whose PSFEx model landed has a non-empty match set under the default + # retention list, so it used to get a green manifest — and clean_exposure + # takes that manifest as its go-ahead and deletes the store, taking the + # stars with it. There is no recovering them afterwards short of rebuilding + # the chain from VOS, so a missing psf_validation fails the job here, while + # the store is still on disk. Retention products that match nothing stay + # warnings: they are optional by construction. + if ALWAYS not in found: + sys.exit(f"persist_exp: {args.exp}: nothing matched {ALWAYS} " + f"({resolve(ALWAYS)}) under {args.exp_dir}/output/{RUN_NAME}. " + f"That is the star catalogue's input and it is not optional — " + f"refusing to write a manifest that would let clean_exposure " + f"reclaim this store.") + if not found: + sys.exit(f"persist_exp: {args.exp}: no file matched any of " + f"{entries} under {args.exp_dir}/output/{RUN_NAME}") + + args.dest.mkdir(parents=True, exist_ok=True) + tar_path = args.dest / f"{args.exp}.tar" + # A file matched by TWO patterns is one file, not a collision. Keep lists + # overlap on purpose — `validation_psf-*.fits` alongside `*.fits` is a + # perfectly ordinary way to say "the validation catalogues, and everything + # else FITS while we are here" — and treating the second match as a name + # clash failed every exposure in the campaign. What must still be fatal is + # two DIFFERENT paths landing on one flat member name, which would silently + # overwrite; that is a same-name/different-source test, and the first + # pattern to match a file is the one recorded for it. + seen, files = {}, [] + for pat, hits in found.items(): + for src in hits: + if src.name in seen: + if seen[src.name][0] == src: + continue # same file, a second matching pattern + sys.exit(f"persist_exp: {args.exp}: two source files are both " + f"named {src.name} ({seen[src.name][0]} and {src}); tar " + f"members are flat, so this would silently overwrite") + seen[src.name] = (src, pat) + files.append({"name": src.name, "product": pat, + "pattern": resolve(pat), + "src": str(src), "bytes": src.stat().st_size}) + + # --- RETENTION IS ADDITIVE: an existing tar is a FLOOR, never a draft ---- + # Shrinking `persist_exp:` used to rerun this rule (the list rides on + # params, which is the whole point of the rule) and overwrite the tar with + # a smaller one — deleting products from the BACKED-UP filesystem because + # someone edited a config. The scratch store they came from is usually gone + # by then, so nothing could put them back. Whatever is already in the tar + # therefore stays in it: a config change can only ever ADD. + # + # Removing a product is consequently not a config edit. It is a deliberate + # act on products_dir, and it should look like one. + carried, prior_products = [], {} + if tar_path.exists(): + prior = args.manifest + if prior.exists(): + try: + prior_products = {f["name"]: f.get("product", "?") + for f in json.loads(prior.read_text())["files"]} + except (OSError, ValueError, KeyError): + pass # a damaged manifest loses only labels + try: + old_read = tarfile.open(tar_path) + except tarfile.TarError as exc: + sys.exit(f"persist_exp: {args.exp}: cannot read the existing " + f"{tar_path}: {exc}. Refusing to write a new one — the " + f"old tar is left exactly as it is, and it may still hold " + f"products nothing else has. Move it aside deliberately " + f"if you have decided it is lost.") + with old_read as tf: + for ti in tf.getmembers(): + if ti.name in seen or not ti.isfile(): + continue # a live source supersedes it + carried.append(ti.name) + files.append({"name": ti.name, + "product": prior_products.get(ti.name, "?"), + "pattern": None, "src": None, "bytes": ti.size}) + + files.sort(key=lambda f: f["name"]) + + def anonymous(ti: tarfile.TarInfo) -> tarfile.TarInfo: + # Ownership is the one thing that would differ between two writes of + # the same files from different accounts/nodes; drop it. mtime stays: + # it is the product's, and it is stable while the store is. + ti.uid = ti.gid = 0 + ti.uname = ti.gname = "" + return ti + + # tmp-then-cmp-then-mv, and the tmp NEVER outlives a failure: an orphaned + # .tmp on /project is an inode nothing revisits — the leak this whole tar + # design exists to avoid, one per failed attempt at DR6 scale. + tmp = tar_path.with_name(tar_path.name + ".tmp") + try: + # One pass in sorted member order, taking each member from whichever + # side has it: a live source on disk, or the existing tar. Members are + # copied across with their own TarInfo, so a carried member is + # byte-for-byte what it was and a rerun that changes nothing still + # produces an identical archive. + with tarfile.open(tmp, "w", format=tarfile.PAX_FORMAT) as tf: + # Already proven readable above, where the members were listed. + old_tar = (tarfile.open(tar_path) if carried else None) + try: + for f in files: + if f["name"] in seen: + tf.add(seen[f["name"]][0], arcname=f["name"], + filter=anonymous) + else: + ti = anonymous(old_tar.getmember(f["name"])) + tf.addfile(ti, old_tar.extractfile(f["name"])) + finally: + if old_tar is not None: + old_tar.close() + if tar_path.exists() and filecmp.cmp(tmp, tar_path, shallow=False): + tmp.unlink() # unchanged: leave the mtime alone + else: + tmp.replace(tar_path) # atomic: no half-written archive + finally: + tmp.unlink(missing_ok=True) + + body = { + "stage": "exp_persist", "level": "exp", "unit": args.exp, + "status": "complete", + "tar": str(tar_path), + "products": entries, + "patterns": [resolve(e) for e in entries], + # The warning the docstring argues for: named patterns that matched + # nothing. Present as a key even when empty, so a reader never has to + # wonder whether an old manifest predates the field. + "patterns_unmatched": empty, + "n_files": len(files), + "bytes": sum(f["bytes"] for f in files), + "files": files, + } + args.manifest.parent.mkdir(parents=True, exist_ok=True) + tmp = args.manifest.with_name(args.manifest.name + ".tmp") + try: + tmp.write_text(json.dumps(body, indent=2, sort_keys=True) + "\n") + if args.manifest.exists() and filecmp.cmp(tmp, args.manifest, shallow=False): + tmp.unlink() # unchanged: leave the mtime alone + else: + tmp.replace(args.manifest) + finally: + tmp.unlink(missing_ok=True) + + warn = (f" ({len(empty)} retention product(s) matched nothing: {empty})" + if empty else "") + if carried: + warn += f" ({len(carried)} member(s) carried from the existing tar)" + print(f"[persist_exp] {args.exp}: {len(files)} file(s), " + f"{body['bytes'] / 1e6:.1f} MB -> {tar_path}{warn}") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/run_report.py b/workflow/scripts/run_report.py index 24c07670b..11cdcffe4 100644 --- a/workflow/scripts/run_report.py +++ b/workflow/scripts/run_report.py @@ -61,6 +61,13 @@ "tile_ngmix", "tile_merge_cats", "tile_make_cat"] EXP_STAGES = ["exp_get_images", "exp_split", "exp_psf"] +# exp_persist is DELIBERATELY NOT in that list. This report disk-scans the +# scratch run_dir, and exp_persist's manifest is the one exposure manifest that +# lives on products_dir instead — that placement is what makes it survive +# clean_exposure. Listed here it would read as "not run" for every exposure in +# the campaign. Reporting on the persisted products means scanning the second +# root, which is a report this one does not yet do. + # The manifests clean_tile leaves on disk (workflow/scripts/clean_tile.py names # the mechanism that owns each). Their presence is therefore NOT evidence that a # tile's chain was rebuilt, which absorb_tombstones needs to know