diff --git a/astra.yaml b/astra.yaml index 4675e1eb7..b4c7def52 100644 --- a/astra.yaml +++ b/astra.yaml @@ -215,6 +215,9 @@ analyses: - id: sky_masks type: data source: UNIONS healsparse mask maps + - id: exposure_psf_ccds + type: data + source: per-exposure split WCS (headers-.npy) and the persisted validation_psf members (exp_persist) outputs: - id: masked_measurement_inputs type: data @@ -225,7 +228,97 @@ analyses: catalogue's MASK_n columns. inputs: [exposure_flags, sky_masks] decisions: [pixel_mask_source, psf_star_mask_veto, sky_mask_application, mask_default_cut] + - id: defect_map + type: data + format: hsp + description: >- + The campaign's instrument-flag defect map: a boolean healsparse + map, True = masked, the union of one fragment per exposure. A + product; nothing in the workflow reads it. Data runs only. + inputs: [exposure_flags] + decisions: [defect_map_from_flags] + - id: nexp_map + type: data + format: hsp + description: >- + The campaign's exposure-count map: a uint16 healsparse map at the + defect map's resolution holding, per pixel, the number of exposures + with a valid PSF model covering it. The count sp_validation's + npoint cut reads. A product; nothing in the workflow reads it. + Built for every fitted-PSF campaign; image simulations on the true + PSF record no footprints and skip it. + inputs: [exposure_psf_ccds] + decisions: [nexp_map_valid_psf_ccds] decisions: + nexp_map_valid_psf_ccds: + label: The exposure-count map stamps only CCDs with a valid PSF model + rationale: >- + exp_footprint records the four sky corners of each CCD whose + validation_psf catalogue exp_persist packed; psfex_interp writes that + catalogue on success and returns without it on NOT_ENOUGH_STARS, + BAD_CHI2 or FILE_NOT_FOUND, so the member list is the valid-PSF set. + Corners come from the split's WCS at pixel edges. nexp_map stamps + value 1 per CCD polygon into a uint16 map; the CCDs of one exposure + do not overlap, so a pixel's value is the number of exposures that + can contribute a PSF-corrected epoch there. It is built from every + record on the products root, reclaimed exposures included, and + written raw (no median smoothing). There is one count: of exposures, + not of pointings. The resolution is the mask ladder's, shared with + the defect map. + Values: + config.yaml#exposure_maps.nexp.enabled = true. + default: valid_psf_ccds + options: + valid_psf_ccds: + label: Count exposures whose CCD has a valid PSF model + all_ccds: + label: Count every exposure whose CCD covers the pixel + excluded: true + excluded_reason: >- + Counts epochs that cannot contribute a shape: a CCD without a + PSF model yields no PSF-corrected measurement, so the count + overstates depth exactly where the PSF fit failed. + defect_map_from_flags: + label: Instrument flags rasterized conservatively into a healsparse defect map + rationale: >- + The per-CCD flag images (bad columns, saturated + pixels, bleed trails; the exposure_flags input) otherwise never leave the + pixel domain, and the footprint built from CCD-corner WCS cannot + subtract them. exp_defect_map samples every nonzero-flag CCD pixel + on an oversample x oversample grid spanning its full extent, + corners included (exposure_maps.defect.oversample 3 in workflow/config.yaml), + maps the samples through the CCD's image-split WCS, and masks every + healpix pixel they touch, at the mask ladder's convention (nside + 131072 over coverage 128, True = masked). Any flag bit counts. + defect_map_merge ORs the campaign's fragments into + defect_map/defect_map_.hsp, with a sidecar recording which exposures + are in; a removed or changed fragment forces a rebuild, since a + union cannot be un-OR-ed. On 2079612p CCD 0, oversample 2, 3 and 5 + miss 169, 54 and 5 of 14229 reference healpix pixels. The map is a + product no workflow stage consumes (config_tile_Mc.ini's + MASK_EXT_PATHS does not name it); image simulations get none, their + flag images being all zero. + Values: + config.yaml#exposure_maps.nside = 131072; + config.yaml#exposure_maps.nside_coverage = 128; + config.yaml#exposure_maps.defect.oversample = 3. + default: rasterize_conservative_os3 + options: + rasterize_conservative_os3: + label: "Any-touch rasterization, 3x3 samples per CCD pixel" + centre_sampling: + label: Mask a healpix pixel only where its centre is flagged + excluded: true + excluded_reason: >- + Erases one-pixel bad columns, the thin geometry #878 exists to + keep. + none_footprint_from_wcs: + label: No defect map; footprint from CCD-corner WCS alone + excluded: true + excluded_reason: >- + The footprint silently includes defective pixels: percent-level + area, in exactly the thin small-scale geometry a window function + needs. pixel_mask_source: label: Only the instrument flag image reaches pixels rationale: >- diff --git a/docs/source/exposure_maps.md b/docs/source/exposure_maps.md new file mode 100644 index 000000000..86462833b --- /dev/null +++ b/docs/source/exposure_maps.md @@ -0,0 +1,33 @@ +# Exposure-level maps + +A workflow campaign builds two HealSparse maps from its exposures, beside the +merged catalogues on the products root. Both are per-exposure records combined +by one campaign job, at the mask ladder's resolution (`nside` 131072 over +`nside_coverage` 128, set once under `exposure_maps:` in +`workflow/config.yaml`), and both are produced automatically whenever the input +can support them — no configuration is needed. + +| map | chain | product | built for | +| --- | --- | --- | --- | +| defect | `exp_split` → `exp_defect_map` → `defect_map_merge` | `defect_map/defect_map_.hsp` | data runs | +| exposure count | `exp_persist` → `exp_footprint` → `nexp_map` | `nexp_map/nexp_map_.hsp` | fitted-PSF runs | + +**The defect map** is boolean, `True` = masked: every healpix pixel touched by +a flagged CCD pixel (bad columns, saturated pixels, bleed trails) of any +exposure. It carries the instrument flags, which otherwise never leave the pixel +domain, into the same form as the sky masks. Image simulations' flag images are +blank, so sims build none. + +**The exposure-count map** counts, per pixel, the exposures with a valid PSF +model covering it — the count behind sp_validation's `npoint >= 3` cut. +`exp_footprint` records the sky corners of each CCD whose PSF fit succeeded; +`nexp_map` stamps every record on the products root into the map, so it grows as +tiles are appended. `psf_model: fake` fits no PSF, so sims build none. +`exposure_maps.nexp.enabled: false` opts a campaign out; the map is rebuilt +whole, so a campaign appended in many small batches may prefer to build it once +at the end. + +Nothing in the workflow reads either map back. Plot the count map by hand with +`plot_coverage_map -i /nexp_map/nexp_map_.hsp ...`, using the +sky windows under `exposure_maps.nexp.plot`. The design and the measurements +are in `workflow/README.md`. diff --git a/docs/source/toc.rst b/docs/source/toc.rst index 01c5a1596..130905306 100644 --- a/docs/source/toc.rst +++ b/docs/source/toc.rst @@ -26,6 +26,7 @@ configuration testing workflow + exposure_maps pipeline_tutorial .. toctree:: diff --git a/pyproject.toml b/pyproject.toml index 53f0dae61..f034ef297 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -85,9 +85,6 @@ canfar_submit_job = "shapepipe.canfar_run:run_job" canfar_monitor = "shapepipe.canfar_run:run_log" canfar_monitor_log = "shapepipe.canfar_run:run_monitor_log" get_ccds_with_psf = "shapepipe.get_ccds_run:run_ccd_psf_handler" -download_headers = "shapepipe.coverage_run:run_download_headers" -extract_field_corners = "shapepipe.coverage_run:run_extract_corners" -build_coverage_map = "shapepipe.coverage_run:run_build_coverage" plot_coverage_map = "shapepipe.coverage_run:run_plot_coverage" [tool.uv] diff --git a/scripts/sh/build_and_plot_coverage_maps.sh b/scripts/sh/build_and_plot_coverage_maps.sh deleted file mode 100755 index b7ea9c8b2..000000000 --- a/scripts/sh/build_and_plot_coverage_maps.sh +++ /dev/null @@ -1,95 +0,0 @@ -#!/bin/bash - -set -euo pipefail - -# Build and plot per-CCD coverage (nexp) maps for a range of catalogue -# versions. Each map counts, per sky pixel, the number of exposures with a -# valid PSF model. The per-version flow is: -# 0. get_ccds_with_psf -> ccds_with_psf_.txt (valid CCD IDs) -# 1. download_headers -> per-exposure header text files -# 2. extract_field_corners-> exp_ra_dec_.txt (per-CCD corners) -# 3. build_coverage_map -> coverage_.x.hsp -# 4. plot_coverage_map -> SGC / NGC plots -# Steps 0-2 are expected to have been run already on canfar (they need -# VOSpace access); this script rebuilds from the extracted per-CCD corners and -# plots. Uncomment the marked lines to run the full chain. - -# Configuration variables - -## Catalogue versions -VERSIONS=("v1.3" "v1.4" "v1.5" "v1.6") - -## Output directory of plots -OUTPUT_DIR="coverage_map_plots" - -# healsparse resolution. nside=131072 gives ~0.1" per pixel, chosen to -# match the UNIONS bit-mask resolution so coverage and mask align pixel- -# wise. CoverageMapBuilder defaults to nside=2048 for lighter-weight -# offline use; override here for the production build. -BUILD_NSIDE=131072 -BUILD_CHANNELS=128 - -# Common parameters -VERBOSE="-v" - -## Colorbar -PLOT_COLORBAR="-C" -PLOT_MIN=1 -PLOT_MAX=5 - -# Plot parameters -## SGC region -SGC_RA_MIN=-20 -SGC_RA_MAX=45 -SGC_DEC_MIN=18 -SGC_DEC_MAX=40 - -## NGC region -NGC_RA_MIN=110 -NGC_RA_MAX=270 -NGC_DEC_MIN=28 -NGC_DEC_MAX=90 - -# Create output directory if it doesn't exist -mkdir -p "${OUTPUT_DIR}" - -# Loop over versions -for VERSION in "${VERSIONS[@]}"; do - echo "Processing ${VERSION}..." - - # Define file paths - CCD_LIST="ccds_with_psf_${VERSION}.txt" - HEADER_DIR="headers_${VERSION}" - INPUT_CORNERS="exp_ra_dec_${VERSION}.txt" - COVERAGE_MAP="coverage_${VERSION}.x.hsp" - - # Step 0-2: build the per-CCD corner file (canfar / VOSpace only). - # Uncomment to run the full chain from scratch. - # get_ccds_with_psf -V "${VERSION}" -o "${CCD_LIST}" ${VERBOSE} - # download_headers -i "${CCD_LIST}" -o "${HEADER_DIR}" ${VERBOSE} - # extract_field_corners -i "${HEADER_DIR}" -l "${CCD_LIST}" \ - # -o "${INPUT_CORNERS}" ${VERBOSE} - - # Build coverage map - echo " Building coverage map from ${INPUT_CORNERS}..." - CMD="build_coverage_map -i ${INPUT_CORNERS} -o ${COVERAGE_MAP} -c ${BUILD_CHANNELS} -n ${BUILD_NSIDE} ${VERBOSE}" - echo "$CMD" - $CMD - - # Plot SGC region - echo " Plotting SGC region..." - CMD="plot_coverage_map -i ${COVERAGE_MAP} -o ${OUTPUT_DIR}/coverage_${VERSION}_SGC.png ${VERBOSE} ${PLOT_COLORBAR} -R ${SGC_RA_MIN} -r ${SGC_RA_MAX} -D ${SGC_DEC_MIN} -d ${SGC_DEC_MAX} -m ${PLOT_MIN} -M ${PLOT_MAX}" - echo "$CMD" - $CMD - - # Plot NGC region - echo " Plotting NGC region..." - CMD="plot_coverage_map -i ${COVERAGE_MAP} -o ${OUTPUT_DIR}/coverage_${VERSION}_NGC.png ${VERBOSE} ${PLOT_COLORBAR} -R ${NGC_RA_MIN} -r ${NGC_RA_MAX} -D ${NGC_DEC_MIN} -d ${NGC_DEC_MAX} -m ${PLOT_MIN} -M ${PLOT_MAX}" - echo "$CMD" - $CMD - - echo " Done with ${VERSION}" - echo -done - -echo "All versions processed successfully!" diff --git a/src/shapepipe/coverage_run.py b/src/shapepipe/coverage_run.py index 15aaa1a37..8b6b14e10 100644 --- a/src/shapepipe/coverage_run.py +++ b/src/shapepipe/coverage_run.py @@ -1,79 +1,20 @@ """COVERAGE_RUN -Call coverage processing classes. +Console entry point for plotting a coverage map. + +Building the exposure-count map is the Snakemake ``nexp_map`` rule's job; +plotting a finished ``.hsp`` is a human act on a durable product, so it stays a +hand-run command, on the same argument that keeps ``run_report.py`` out of the +DAG. The plot windows for the UNIONS SGC and NGC fields are in +``workflow/config.yaml``'s ``exposure_maps.nexp.plot`` block. Author: Martin Kilbinger """ -import sys - -from shapepipe.utilities.header_downloader import HeaderDownloader -from shapepipe.utilities.field_corners_extractor import FieldCornersExtractor -from shapepipe.utilities.coverage_map_builder import CoverageMapBuilder from shapepipe.utilities.coverage_plotter import CoveragePlotter -def run_download_headers(args=None): - """Run Download Headers. - - Download FITS headers from VOSpace for exposures in a CCD list. - - Parameters - ---------- - args : list, optional - command line arguments - - Returns - ------- - int - exit code - - """ - obj = HeaderDownloader() - return obj.run(args=args) - - -def run_extract_corners(args=None): - """Run Extract Corners. - - Extract per-CCD sky-footprint corner coordinates from FITS headers. - - Parameters - ---------- - args : list, optional - command line arguments - - Returns - ------- - int - exit code - - """ - obj = FieldCornersExtractor() - return obj.run(args=args) - - -def run_build_coverage(args=None): - """Run Build Coverage. - - Build HealSparse coverage maps from per-CCD corner coordinates. - - Parameters - ---------- - args : list, optional - command line arguments - - Returns - ------- - int - exit code - - """ - obj = CoverageMapBuilder() - return obj.run(args=args) - - def run_plot_coverage(args=None): """Run Plot Coverage. @@ -92,13 +33,3 @@ def run_plot_coverage(args=None): """ obj = CoveragePlotter() return obj.run(args=args) - - -def main(argv=None): - """Main. - - Main program. - - """ - # Scripts to call coverage classes are created by pyproject.toml - return 0 diff --git a/src/shapepipe/utilities/ccd_footprint.py b/src/shapepipe/utilities/ccd_footprint.py new file mode 100644 index 000000000..3e176998b --- /dev/null +++ b/src/shapepipe/utilities/ccd_footprint.py @@ -0,0 +1,90 @@ +"""CCD_FOOTPRINT + +Project a single CCD onto the sky: image shape from a header, sky corners +from a WCS. + +The two functions here are the geometry the exposure-count map rests on. The +``exp_footprint`` rule (``workflow/scripts/exp_footprint.py``) calls them once +per CCD of an exposure, reading the headers out of ``headers-.npy``; the +``nexp_map`` rule then stamps the corners into the campaign's HealSparse map. + +Author: Mike Hudson, Martin Kilbinger + +""" + + +def _image_shape(header): + """Return the CCD image ``(nx, ny)`` pixel shape from a header. + + For fpack tile-compressed HDUs (``ZIMAGE = T``), ``NAXIS1``/``NAXIS2`` + describe the compressed *binary table* (row width in bytes, row count), + not the image, and astropy's WCS never maps the true dimensions into + ``pixel_shape``. The true image size is carried by ``ZNAXIS1``/``ZNAXIS2``, + which are preferred here; a plain image HDU falls back to + ``NAXIS1``/``NAXIS2``. + + Both branches are live. ``headers-.npy`` carries the *decompressed* + header astropy hands back for a tile-compressed HDU, so it takes the plain + ``NAXIS`` path; a header read straight off a fpacked file takes the ``Z*`` + one. + + Parameters + ---------- + header : astropy.io.fits.Header + single-HDU header + + Returns + ------- + tuple + ``(nx, ny)`` image pixel dimensions + + Raises + ------ + ValueError + if neither the ``Z*`` nor the plain ``NAXIS`` dimensions are present + """ + if header.get("ZIMAGE", False): + keys = ("ZNAXIS1", "ZNAXIS2") + else: + keys = ("NAXIS1", "NAXIS2") + + if keys[0] not in header or keys[1] not in header: + raise ValueError( + f"Header is missing image dimensions ({keys[0]}/{keys[1]})" + ) + + return int(header[keys[0]]), int(header[keys[1]]) + + +def _ccd_corners(w, shape): + """Return (ra, dec) of the 4 corners of a single CCD. + + The polygon is built from the CCD's pixel bounds ``shape = (nx, ny)`` using + pixel edges ``(-0.5, n - 0.5)`` so the quadrilateral covers the full CCD + area rather than pixel centres. Corners are returned in a consistent + counterclockwise pixel order (bottom-left, bottom-right, top-right, + top-left); since ``pixel_to_world`` maps this order to a simple, + non-self-intersecting sky quadrilateral for a CCD-sized field, and + HealSparse ``Polygon`` is orientation-agnostic, the resulting polygon is a + valid convex footprint. + + Parameters + ---------- + w : astropy.wcs.WCS + WCS of a single CCD + shape : tuple + ``(nx, ny)`` CCD image pixel dimensions + + Returns + ------- + tuple + ``(ra_list, dec_list)`` each of length 4, in degrees + """ + nx, ny = shape + + # Pixel-edge corners, counterclockwise: BL, BR, TR, TL + px = [-0.5, nx - 0.5, nx - 0.5, -0.5] + py = [-0.5, -0.5, ny - 0.5, ny - 0.5] + + sky = w.pixel_to_world(px, py) + return list(sky.ra.deg), list(sky.dec.deg) diff --git a/src/shapepipe/utilities/coverage_map_builder.py b/src/shapepipe/utilities/coverage_map_builder.py index 5031c8791..e62fc7ca0 100644 --- a/src/shapepipe/utilities/coverage_map_builder.py +++ b/src/shapepipe/utilities/coverage_map_builder.py @@ -1,25 +1,25 @@ """COVERAGE_MAP_BUILDER -Build HealSparse coverage maps from per-CCD corner coordinates. +Stamp per-CCD sky footprints into a HealSparse exposure-count (nexp) map. -Each input row is one CCD footprint. Since the CCDs of a single exposure do -not overlap, stamping value 1 per CCD polygon and accumulating makes the map -count, per sky pixel, the number of exposures with a valid PSF model covering -that pixel. +Each polygon is one CCD footprint. Since the CCDs of a single exposure do not +overlap, stamping value 1 per CCD polygon and accumulating makes the map count, +per sky pixel, the number of exposures with a valid PSF model covering that +pixel. + +Arrays in, map out. The only caller is the ``nexp_map`` rule +(``workflow/scripts/nexp_map.py``), which reads the corners out of the +per-exposure ``exp_footprint.json`` records. What is here is survey geometry — +the RA-seam guard, the pole guard, the nside validation and the median +smoothing — none of which depends on how the corners arrived. Author: Mike Hudson, Martin Kilbinger """ -import sys -from os.path import exists - -import numpy as np import healsparse as hsp import hpgeom as hpg - -from cs_util import args as cs_args -from cs_util import logging +import numpy as np # Declination beyond which a polygon is considered too close to a pole for # HealSparse's planar polygon fill to be reliable. CCD footprints are ~10 @@ -52,163 +52,134 @@ def unwrap_ra(ra): return ra +def check_nside(nside_coverage, nside): + """Raise unless both HealSparse resolutions are positive powers of two. + + Checked before any polygon is stamped: healsparse itself complains only + once a map is being made, by which time the caller has already paid for + reading every footprint record on the products root. -class CoverageMapBuilder(object): - """Coverage Map Builder Class. + Parameters + ---------- + nside_coverage : int + HealSparse coverage nside + nside : int + HealSparse map nside + + Raises + ------ + ValueError + if either value is not a positive power of two + + """ + for name, n in (("nside_coverage", nside_coverage), ("nside", nside)): + if n <= 0 or n & (n - 1): + raise ValueError(f"{name} must be a power of 2, got {n}") + + +def build_map( + ccd_ids, ra, dec, nside_coverage, nside, verbose=False +): + """Stamp one polygon per CCD footprint into a HealSparse nexp map. + + Since the CCDs of a single exposure do not overlap, accumulating value 1 + per CCD polygon makes the map count, per sky pixel, the number of + exposures with a valid PSF model covering it. + + Parameters + ---------- + ccd_ids : array_like + length-N CCD IDs, used only in the skipped-polygon warnings + ra : array_like + ``(N, 4)`` polygon RA corners, in degrees + dec : array_like + ``(N, 4)`` polygon Dec corners, in degrees + nside_coverage : int + HealSparse coverage nside + nside : int + HealSparse map nside + verbose : bool, optional + print progress + + Returns + ------- + healsparse.HealSparseMap + the nexp map - Builds HealSparse coverage maps from field corner coordinates. """ + check_nside(nside_coverage, nside) - def __init__(self): - """Initialize the builder.""" - self.params_default() - - def params_default(self): - """Set default parameters and command line options.""" - - self._params = { - "input_file": "exp_ra_dec.txt", - "output_file": "coverage.hsp", - "nside_coverage": 32, - "nside": 2048, - "apply_median_filter": False, - "n_median_iterations": 2, - "create_boolean": False, - "boolean_threshold": 3, - "boolean_output": None, - "create_plot": False, - "plot_output": None, - "plot_region": None, - "verbose": False, - } - - self._short_options = { - "input_file": "-i", - "output_file": "-o", - "nside_coverage": "-c", - "nside": "-n", - "apply_median_filter": "-m", - "n_median_iterations": "-N", - "create_boolean": "-b", - "boolean_threshold": "-t", - "boolean_output": "-B", - "create_plot": "-p", - "plot_output": "-P", - "plot_region": "-g", - } - - self._types = { - "nside_coverage": "int", - "nside": "int", - "apply_median_filter": "bool", - "n_median_iterations": "int", - "create_boolean": "bool", - "boolean_threshold": "int", - "create_plot": "bool", - "verbose": "bool", - } - - self._help_strings = { - "input_file": "input file with field corners; default is {}", - "output_file": "output HealSparse map file; default is {}", - "nside_coverage": "HealSparse coverage nside; default is {}", - "nside": "HealSparse map nside; default is {}", - "apply_median_filter": "apply median filter to smooth map; default is {}", - "n_median_iterations": "number of median filter iterations; default is {}", - "create_boolean": "create boolean coverage map; default is {}", - "boolean_threshold": "threshold value for boolean map; default is {}", - "boolean_output": "output file for boolean map; default is _bool.hsp", - "create_plot": "create plot of coverage map; default is {}", - "plot_output": "output file for plot; default is .png", - "plot_region": "predefined region for plot (NGC, SGC, fullsky); default is {}", - } - - def set_params_from_command_line(self, args): - """Set Params From Command line. - - Only use when calling using python from command line. - Does not work from ipython or jupyter. - - Parameters - ---------- - args : list - command line arguments - - """ - # Read command line options - options = cs_args.parse_options( - self._params, - self._short_options, - self._types, - self._help_strings, - args=args, + ra = np.atleast_2d(np.asarray(ra, dtype=float)) + dec = np.atleast_2d(np.asarray(dec, dtype=float)) + + if verbose: + print( + f"Creating HealSparse map (nside_coverage={nside_coverage}," + f" nside={nside})" ) - self._params = options - - # Save calling command to log_build_coverage_map; args excludes the program - # name, which log_command takes from argv[0] - logging.log_command(["build_coverage_map", *args]) - - def update_params(self): - """Update parameters. - - Set derived parameters based on input parameters. - """ - # Set boolean output filename if not specified - if self._params["boolean_output"] is None: - output_file = self._params["output_file"] - if output_file.endswith(".hsp"): - base = output_file[:-4] - else: - base = output_file - self._params["boolean_output"] = f"{base}_bool.hsp" - - # Set plot output filename if not specified - if self._params["plot_output"] is None: - output_file = self._params["output_file"] - if output_file.endswith(".hsp"): - base = output_file[:-4] - else: - base = output_file - self._params["plot_output"] = f"{base}.png" - - def check_params(self): - """Check parameters for validity.""" - if not exists(self._params["input_file"]): - raise FileNotFoundError( - f"Input file not found: {self._params['input_file']}" - ) + m = hsp.HealSparseMap.make_empty(nside_coverage, nside, np.uint16) + + if verbose: + print("Adding polygons to map") - # Check nside values are powers of 2 - nside_coverage = self._params["nside_coverage"] - nside = self._params["nside"] + n_added = 0 + n_skipped = 0 + for i in range(len(ccd_ids)): + dec_i = dec[i] - if not (nside_coverage & (nside_coverage - 1) == 0): - raise ValueError( - f"nside_coverage must be a power of 2, got {nside_coverage}" + # Pole guard: HealSparse's planar polygon fill degrades near the + # poles. CCD footprints never reach here, so warn and skip. + if np.any(np.abs(dec_i) >= _DEC_POLE_LIMIT): + print( + f"Warning: skipping CCD {ccd_ids[i]} with |dec| >= " + f"{_DEC_POLE_LIMIT} (too close to a pole)" ) + n_skipped += 1 + continue + + # RA-wrap guard: put corners on a common branch across the seam. + ra_i = unwrap_ra(ra[i]) - if not (nside & (nside - 1) == 0): - raise ValueError(f"nside must be a power of 2, got {nside}") + m += hsp.Polygon(ra=list(ra_i), dec=list(dec_i), value=1) + n_added += 1 - def median_filter(self, hsp_map): - """Median Filter. + if verbose and i % 1000 == 0: + print(f"{i:6d} / {len(ccd_ids):6d}") - Apply median filter to HealSparse map using neighbors. + print(f"Added {n_added} polygons to map") + if n_skipped > 0: + print(f"Skipped {n_skipped} polygons near a pole") - Parameters - ---------- - hsp_map : healsparse.HealSparseMap - input map + return m + + +def median_filter(hsp_map, n_iterations=1): + """Smooth an nexp map: each pixel becomes its neighbourhood median. + + The neighbourhood is a pixel and its eight HEALPix neighbours; out-of-map + neighbours read as the map's sentinel and so pull an edge pixel down, + which is the intended behaviour for a coverage map (a lone pixel is noise, + not depth). The ``nexp_map`` rule does NOT smooth — the map it writes + is the raw exposure count, which is what sp_validation's mask application + wants — so this is the offline step, for looking at a map rather than + applying one. + + Parameters + ---------- + hsp_map : healsparse.HealSparseMap + input nexp map + n_iterations : int, optional + number of smoothing passes - Returns - ------- - healsparse.HealSparseMap - filtered map + Returns + ------- + healsparse.HealSparseMap + filtered map - """ - nside = hsp_map.nside_sparse + """ + nside = hsp_map.nside_sparse + for _ in range(n_iterations): new_hsp = hsp_map.copy() pixs = hsp_map.valid_pixels n = hpg.neighbors(nside, pixs) @@ -218,156 +189,6 @@ def median_filter(self, hsp_map): mn = hsp_map[n] new = np.median(mn, axis=1).astype(np.uint16) new_hsp[pixs] = new + hsp_map = new_hsp - return new_hsp - - def run(self, args=None): - """Run. - - Main execution method. - - Parameters - ---------- - args : list, optional - command line arguments - - Returns - ------- - int - exit code (0 for success) - - """ - if args is None: - args = sys.argv[1:] - - # Set parameters from command line - self.set_params_from_command_line(args) - self.update_params() - self.check_params() - - # Get parameters - input_file = self._params["input_file"] - output_file = self._params["output_file"] - nside_coverage = self._params["nside_coverage"] - nside = self._params["nside"] - apply_median_filter = self._params["apply_median_filter"] - n_median_iterations = self._params["n_median_iterations"] - create_boolean = self._params["create_boolean"] - boolean_threshold = self._params["boolean_threshold"] - boolean_output = self._params["boolean_output"] - create_plot = self._params["create_plot"] - plot_output = self._params["plot_output"] - plot_region = self._params["plot_region"] - verbose = self._params["verbose"] - - if verbose: - print(f"Reading per-CCD corners from: {input_file}") - - # Load per-CCD corner data. Column 0 is a "-" string - # ID; the remaining 8 columns are the 4 RA then 4 Dec corners. - ccd_ids = np.atleast_1d(np.loadtxt(input_file, usecols=(0), dtype=str)) - corners = np.atleast_2d( - np.loadtxt(input_file, usecols=range(1, 9), dtype=float) - ) - ra_all = corners[:, :4] - dec_all = corners[:, 4:] - - print(f"Loaded {len(ccd_ids)} CCD footprints") - - if verbose: - print( - f"Creating HealSparse map (nside_coverage={nside_coverage}, nside={nside})" - ) - - # Create empty map - m = hsp.HealSparseMap.make_empty(nside_coverage, nside, np.uint16) - - # Add one polygon per CCD footprint - if verbose: - print("Adding polygons to map") - - n_added = 0 - n_skipped = 0 - for i in range(len(ccd_ids)): - dec = dec_all[i] - - # Pole guard: HealSparse's planar polygon fill degrades near the - # poles. CCD footprints never reach here, so warn and skip. - if np.any(np.abs(dec) >= _DEC_POLE_LIMIT): - print( - f"Warning: skipping CCD {ccd_ids[i]} with |dec| >= " - f"{_DEC_POLE_LIMIT} (too close to a pole)" - ) - n_skipped += 1 - continue - - # RA-wrap guard: put corners on a common branch across the seam. - ra = unwrap_ra(ra_all[i]) - - m += hsp.Polygon(ra=list(ra), dec=list(dec), value=1) - n_added += 1 - - if verbose and i % 1000 == 0: - print(f"{i:6d} / {len(ccd_ids):6d}") - - print(f"Added {n_added} polygons to map") - if n_skipped > 0: - print(f"Skipped {n_skipped} polygons near a pole") - - # Apply median filter if requested - if apply_median_filter: - if verbose: - print( - f"Applying median filter ({n_median_iterations} iterations)" - ) - - for i in range(n_median_iterations): - m = self.median_filter(m) - if verbose: - print(f" Iteration {i+1}/{n_median_iterations} complete") - - # Write main coverage map - if verbose: - print(f"Writing coverage map to: {output_file}") - - m.write(output_file, clobber=True) - print(f"Coverage map saved to {output_file}") - - # Create and write boolean map if requested - if create_boolean: - if verbose: - print( - f"Creating boolean map (threshold={boolean_threshold})" - ) - - c = hsp.HealSparseMap.make_empty(nside_coverage, nside, bool) - c[m.valid_pixels] = m[m.valid_pixels] >= boolean_threshold - - if verbose: - print(f"Writing boolean map to: {boolean_output}") - - c.write(boolean_output, clobber=True) - print(f"Boolean map saved to {boolean_output}") - - # Create plot if requested - if create_plot: - if verbose: - print(f"Creating plot of coverage map") - - try: - from shapepipe.utilities.coverage_plotter import CoveragePlotter - plotter = CoveragePlotter() - plotter.plot_coverage_map( - m, - plot_output, - region=plot_region, - vmax=boolean_threshold if create_boolean else 3, - colorbar=True, - colorbar_label="Coverage depth", - ) - print(f"Plot saved to {plot_output}") - except ImportError as e: - print(f"Warning: Could not create plot: {e}") - print("Install the cs_util package for plotting support.") - - return 0 + return hsp_map diff --git a/src/shapepipe/utilities/field_corners_extractor.py b/src/shapepipe/utilities/field_corners_extractor.py deleted file mode 100644 index a9b4e6c40..000000000 --- a/src/shapepipe/utilities/field_corners_extractor.py +++ /dev/null @@ -1,494 +0,0 @@ -"""FIELD_CORNERS_EXTRACTOR - -Extract per-CCD sky-footprint corner coordinates from FITS headers. - -For each CCD (image HDU) in an exposure's multi-HDU header, the four corners -of the CCD are projected from pixel bounds to the sky, producing one row per -CCD keyed on its ``-`` ID. - -Author: Mike Hudson, Martin Kilbinger - -""" - -import sys -import os -import glob -import re -from os.path import exists -from multiprocessing import Pool, cpu_count - -import numpy as np -from astropy import wcs -from astropy.io.fits import Header - -from cs_util import args as cs_args -from cs_util import logging - -# Filename suffix convention for header files: <6+ digit exposure number>.txt -_EXPNUM_RE = re.compile(r'(\d+)\.txt$') - - -def _expnum_from_path(path): - """Extract exposure number from a header filename like ``1234567.txt``.""" - match = _EXPNUM_RE.search(path) - if match is None: - raise ValueError(f"Could not extract exposure number from {path!r}") - return int(match.group(1)) - - -def _image_shape(header): - """Return the CCD image ``(nx, ny)`` pixel shape from a header. - - For fpack tile-compressed HDUs (``ZIMAGE = T``), ``NAXIS1``/``NAXIS2`` - describe the compressed *binary table* (row width in bytes, row count), - not the image, and astropy's WCS never maps the true dimensions into - ``pixel_shape``. The true image size is carried by ``ZNAXIS1``/``ZNAXIS2``, - which are preferred here; a plain image HDU falls back to - ``NAXIS1``/``NAXIS2``. - - Parameters - ---------- - header : astropy.io.fits.Header - single-HDU header - - Returns - ------- - tuple - ``(nx, ny)`` image pixel dimensions - - Raises - ------ - ValueError - if neither the ``Z*`` nor the plain ``NAXIS`` dimensions are present - """ - if header.get("ZIMAGE", False): - keys = ("ZNAXIS1", "ZNAXIS2") - else: - keys = ("NAXIS1", "NAXIS2") - - if keys[0] not in header or keys[1] not in header: - raise ValueError( - f"Header is missing image dimensions ({keys[0]}/{keys[1]})" - ) - - return int(header[keys[0]]), int(header[keys[1]]) - - -def _parse_header_to_wcs(path): - """Parse a multi-HDU header text file into ``(wcs, (nx, ny))`` per HDU. - - The primary HDU is skipped. Each remaining HDU yields its WCS together with - the CCD image pixel shape (``ZNAXIS1/2`` for compressed HDUs, else - ``NAXIS1/2``); the shape is read from the header because the WCS drops the - ``Z*`` keywords. - """ - with open(path, "r") as f: - string = f.read() - tokens = re.split(r"^(END\s+)", string, flags=re.MULTILINE) - result = [] - for i in range(2, len(tokens) - 1, 2): - header = Header.fromstring(tokens[i] + tokens[i + 1], sep="\n") - result.append((wcs.WCS(header), _image_shape(header))) - return result - - -def _ccd_corners(w, shape): - """Return (ra, dec) of the 4 corners of a single CCD. - - The polygon is built from the CCD's pixel bounds ``shape = (nx, ny)`` using - pixel edges ``(-0.5, n - 0.5)`` so the quadrilateral covers the full CCD - area rather than pixel centres. Corners are returned in a consistent - counterclockwise pixel order (bottom-left, bottom-right, top-right, - top-left); since ``pixel_to_world`` maps this order to a simple, - non-self-intersecting sky quadrilateral for a CCD-sized field, and - HealSparse ``Polygon`` is orientation-agnostic, the resulting polygon is a - valid convex footprint. - - Parameters - ---------- - w : astropy.wcs.WCS - WCS of a single CCD - shape : tuple - ``(nx, ny)`` CCD image pixel dimensions - - Returns - ------- - tuple - ``(ra_list, dec_list)`` each of length 4, in degrees - """ - nx, ny = shape - - # Pixel-edge corners, counterclockwise: BL, BR, TR, TL - px = [-0.5, nx - 0.5, nx - 0.5, -0.5] - py = [-0.5, -0.5, ny - 0.5, ny - 0.5] - - sky = w.pixel_to_world(px, py) - return list(sky.ra.deg), list(sky.dec.deg) - - -class FieldCornersExtractor(object): - """Field Corners Extractor Class. - - Extracts RA/Dec coordinates of field corners from FITS headers. - """ - - def __init__(self): - """Initialize the extractor.""" - self.params_default() - - def params_default(self): - """Set default parameters and command line options.""" - - self._params = { - "input_dir": "header", - "output_file": "exp_ra_dec.txt", - "ccd_list": None, - "resume": False, - "n_processes": 1, - "verbose": False, - } - - self._short_options = { - "input_dir": "-i", - "output_file": "-o", - "ccd_list": "-l", - "resume": "-r", - "n_processes": "-n", - } - - self._types = { - "resume": "bool", - "n_processes": "int", - "verbose": "bool", - } - - self._help_strings = { - "input_dir": "input directory containing header files; default is {}", - "output_file": "output file for per-CCD corners; default is {}", - "ccd_list": "file of valid CCD IDs (output of get_ccds_with_psf); when given, only listed CCDs are written; on --resume, CCDs new to an expanded list are added without duplicating existing rows; default is all CCDs", - "resume": "resume from existing output file; default is {}", - "n_processes": f"number of parallel processes (1=serial, 0=auto={cpu_count()}); default is {{}}", - } - - def set_params_from_command_line(self, args): - """Set Params From Command line. - - Only use when calling using python from command line. - Does not work from ipython or jupyter. - - Parameters - ---------- - args : list - command line arguments - - """ - # Read command line options - options = cs_args.parse_options( - self._params, - self._short_options, - self._types, - self._help_strings, - args=args, - ) - self._params = options - - # Save calling command to log_extract_field_corners; args excludes the program - # name, which log_command takes from argv[0] - logging.log_command(["extract_field_corners", *args]) - - def update_params(self): - """Update parameters. - - Set derived parameters based on input parameters. - """ - # Ensure input directory ends without trailing slash - if self._params["input_dir"].endswith("/"): - self._params["input_dir"] = self._params["input_dir"][:-1] - - def check_params(self): - """Check parameters for validity.""" - if not exists(self._params["input_dir"]): - raise FileNotFoundError( - f"Input directory not found: {self._params['input_dir']}" - ) - - if self._params["ccd_list"] is not None and not exists( - self._params["ccd_list"] - ): - raise FileNotFoundError( - f"CCD list file not found: {self._params['ccd_list']}" - ) - - # Set n_processes to cpu_count if 0 - if self._params["n_processes"] == 0: - self._params["n_processes"] = cpu_count() - - if self._params["n_processes"] < 0: - raise ValueError( - f"n_processes must be >= 0, got {self._params['n_processes']}" - ) - - @staticmethod - def load_ccd_list(path): - """Load CCD List. - - Read valid CCD IDs (one ``-`` per line) into a set. - - Parameters - ---------- - path : str - path to the CCD list file - - Returns - ------- - set - valid CCD IDs - - """ - with open(path) as f: - return {line.strip() for line in f if line.strip()} - - @staticmethod - def process_single_header(args): - """Process Single Header. - - Worker function to process a single header file into per-CCD corners. - Static method so it can be pickled for multiprocessing. - - Parameters - ---------- - args : tuple - ``(path, verbose, valid_ccds)`` where ``path`` is the header file - path, ``verbose`` is a bool, and ``valid_ccds`` is a set of CCD IDs - to keep (or ``None`` to keep all) - - Returns - ------- - list or None - list of ``(ccd_id, ra_list, dec_list)`` for the exposure's CCDs on - success, ``None`` on failure - - """ - path, verbose, valid_ccds = args - expnum = _expnum_from_path(path) - - try: - wcs_shapes = _parse_header_to_wcs(path) - except Exception as e: - if verbose: - print(f"Failed to process {expnum}: {e}") - return None - - rows = [] - for ccd_idx, (w, shape) in enumerate(wcs_shapes): - ccd_id = f"{expnum}-{ccd_idx}" - if valid_ccds is not None and ccd_id not in valid_ccds: - continue - try: - ra, dec = _ccd_corners(w, shape) - except Exception as e: - if verbose: - print(f"Failed to process CCD {ccd_id}: {e}") - continue - rows.append((ccd_id, ra, dec)) - - return rows - - def get_done_ccds(self): - """Get Done CCDs. - - Read the set of CCD IDs already present in the output file. Resume is - keyed on individual CCD IDs, not exposure numbers: a CCD counts as done - only if its own row is present. This keeps resume correct when a write - was interrupted mid-exposure (the missing CCDs are filled in) and when - a rerun uses an expanded ``--ccd_list`` (the newly requested CCDs are - added), and it never duplicates a row. - - Returns - ------- - set - CCD IDs already written - - """ - output_file = self._params["output_file"] - - if not exists(output_file): - return set() - - try: - ids = np.atleast_1d( - np.loadtxt(output_file, usecols=(0), dtype=str) - ) - return set(ids.tolist()) - except Exception as e: - if self._params["verbose"]: - print(f"Could not read existing output file: {e}") - return set() - - def run(self, args=None): - """Run. - - Main execution method. - - Parameters - ---------- - args : list, optional - command line arguments - - Returns - ------- - int - exit code (0 for success) - - """ - if args is None: - args = sys.argv[1:] - - # Set parameters from command line - self.set_params_from_command_line(args) - self.update_params() - self.check_params() - - # Get parameters - input_dir = self._params["input_dir"] - output_file = self._params["output_file"] - resume = self._params["resume"] - verbose = self._params["verbose"] - - # Load the set of valid CCD IDs to keep, if given - valid_ccds = None - if self._params["ccd_list"] is not None: - valid_ccds = self.load_ccd_list(self._params["ccd_list"]) - print(f"{len(valid_ccds)} valid CCDs in {self._params['ccd_list']}") - - # Find all header files - paths = glob.glob(f"{input_dir}/*.txt") - n = len(paths) - - if n == 0: - print(f"No header files found in {input_dir}/") - return 1 - - print(f"{n} header files found") - - # On resume, read the CCD IDs already written so we can skip them - # per-CCD (not per-exposure): a partially written exposure is completed - # rather than skipped, and no row is ever duplicated. - done_ccds = set() - if resume: - done_ccds = self.get_done_ccds() - print(f"{len(done_ccds)} CCDs already done") - - # Every header still needs parsing (a header may hold both done and - # not-yet-done CCDs); the per-CCD filter below decides what is written. - todo = paths - - # Get n_processes - n_processes = self._params["n_processes"] - - # Process headers - if n_processes == 1: - # Serial processing - results = self._process_serial(todo, verbose, valid_ccds) - else: - # Parallel processing - print(f"Using {n_processes} parallel processes") - results = self._process_parallel( - todo, n_processes, verbose, valid_ccds - ) - - # Flatten per-exposure CCD lists (dropping failed exposures), then drop - # CCDs already present in the output. - rows = [ - row - for res in results if res is not None - for row in res - if row[0] not in done_ccds - ] - - # Sort by exposure number, then CCD index - rows.sort(key=lambda r: (int(r[0].split("-")[0]), - int(r[0].split("-")[1]))) - - # Write results to file: " ra1 ra2 ra3 ra4 dec1 dec2 dec3 dec4" - mode = "a" if resume else "w" - with open(output_file, mode, buffering=1) as f: - for ccd_id, ra, dec in rows: - f.write(f"{ccd_id} ") - np.savetxt(f, ra, fmt="%9.5f", newline=" ") - np.savetxt(f, dec, fmt="%9.5f", newline=" ") - f.write("\n") - - n_exp_success = sum(1 for res in results if res is not None) - n_failed = len(todo) - n_exp_success - - print(f"Processed {n_exp_success} exposures, {len(rows)} new CCDs") - if n_failed > 0: - print(f"Failed to process {n_failed} exposures") - - print(f"Results written to {output_file}") - - return 0 - - def _process_serial(self, todo, verbose, valid_ccds): - """Process Serial. - - Process headers serially. - - Parameters - ---------- - todo : list - list of header file paths to process - verbose : bool - verbose output - valid_ccds : set or None - CCD IDs to keep, or ``None`` to keep all - - Returns - ------- - list - list of per-exposure CCD-row lists (or ``None`` for failures) - - """ - results = [] - n_todo = len(todo) - - for i, p in enumerate(todo): - result = self.process_single_header((p, verbose, valid_ccds)) - results.append(result) - - if verbose and i % 100 == 0: - print(f"{i:6d} / {n_todo:6d}") - - return results - - def _process_parallel(self, todo, n_processes, verbose, valid_ccds): - """Process Parallel. - - Process headers in parallel using multiprocessing. - - Parameters - ---------- - todo : list - list of header file paths to process - n_processes : int - number of parallel processes - verbose : bool - verbose output - valid_ccds : set or None - CCD IDs to keep, or ``None`` to keep all - - Returns - ------- - list - list of per-exposure CCD-row lists (or ``None`` for failures) - - """ - # Prepare arguments for worker function - args = [(p, verbose, valid_ccds) for p in todo] - - # Create pool and process - with Pool(processes=n_processes) as pool: - results = pool.map(self.process_single_header, args) - - return results diff --git a/src/shapepipe/utilities/header_downloader.py b/src/shapepipe/utilities/header_downloader.py deleted file mode 100644 index 21b9726b2..000000000 --- a/src/shapepipe/utilities/header_downloader.py +++ /dev/null @@ -1,292 +0,0 @@ -"""HEADER_DOWNLOADER - -Download FITS headers from VOSpace for exposures listed in a CCD file. - -Author: Mike Hudson, Martin Kilbinger - -""" - -import sys -import os -from os.path import exists - -import numpy as np -import vos - -from cs_util import args as cs_args -from cs_util import logging - - -class HeaderDownloader(object): - """Header Downloader Class. - - Downloads FITS headers from VOSpace for exposures in a CCD list. - """ - - def __init__(self): - """Initialize the downloader.""" - self.params_default() - - def params_default(self): - """Set default parameters and command line options.""" - - self._params = { - "input_file": None, - "output_dir": "header", - "vospace_path": "vos:cfis/pitcairn", # UNIONS/CFIS default; override for other surveys - "overwrite": False, - "dir_for_links": None, - "verbose": False, - } - - self._short_options = { - "input_file": "-i", - "output_dir": "-o", - "vospace_path": "-p", - "overwrite": "-O", - "dir_for_links": "-d", - } - - self._types = { - "overwrite": "bool", - "verbose": "bool", - } - - self._help_strings = { - "input_file": "input CCD list file (txt or csv); required", - "output_dir": "output directory for headers; default is {}", - "vospace_path": "VOSpace base path; default is {}", - "overwrite": "overwrite existing header files; default is {}", - "dir_for_links": "directory to check for existing headers to link instead of download; default is {}", - } - - def set_params_from_command_line(self, args): - """Set Params From Command line. - - Only use when calling using python from command line. - Does not work from ipython or jupyter. - - Parameters - ---------- - args : list - command line arguments - - """ - # Read command line options - options = cs_args.parse_options( - self._params, - self._short_options, - self._types, - self._help_strings, - args=args, - ) - self._params = options - - # Save calling command to log_download_headers; args excludes the program - # name, which log_command takes from argv[0] - logging.log_command(["download_headers", *args]) - - def update_params(self): - """Update parameters. - - Set derived parameters based on input parameters. - """ - # Ensure output directory ends without trailing slash - if self._params["output_dir"].endswith("/"): - self._params["output_dir"] = self._params["output_dir"][:-1] - - def check_params(self): - """Check parameters for validity.""" - if self._params["input_file"] is None: - raise ValueError("Input file is required (use -i or --input_file)") - - if not exists(self._params["input_file"]): - raise FileNotFoundError( - f"Input file not found: {self._params['input_file']}" - ) - - # Check if dir_for_links exists if specified - if self._params["dir_for_links"] is not None: - if not exists(self._params["dir_for_links"]): - raise FileNotFoundError( - f"Directory for links not found: {self._params['dir_for_links']}" - ) - if not os.path.isdir(self._params["dir_for_links"]): - raise ValueError( - f"Path is not a directory: {self._params['dir_for_links']}" - ) - - # Create output directory if it doesn't exist - if not exists(self._params["output_dir"]): - os.makedirs(self._params["output_dir"]) - if self._params["verbose"]: - print(f"Created output directory: {self._params['output_dir']}") - - def get_exposures(self, ccd_list_file): - """Get Exposures. - - Extract unique exposure numbers from CCD list file. - - Parameters - ---------- - ccd_list_file : str - path to CCD list file (txt or csv) - - Returns - ------- - np.array - unique exposure numbers - - """ - exps = [] - - # Check if CSV format - if ccd_list_file.endswith(".csv"): - import astropy.table - - t = astropy.table.Table.read(ccd_list_file) - expccd = t["CCD"].data - r = np.char.split(expccd, sep="-") - - for i, r1 in enumerate(r): - exp = int(r1[0]) - exps.append(exp) - else: - # Text format; atleast_1d so a single-line file (0-d array) is - # still iterable. - f = np.atleast_1d( - np.loadtxt(ccd_list_file, dtype="str", encoding="ascii") - ) - r = np.char.split(f, sep="-") - - for i, r1 in enumerate(r): - exp = int(r1[0]) - exps.append(exp) - - exps = np.array(exps) - uniq = np.unique(exps) - - return uniq - - def get_fits_header(self, expnum, client): - """Get FITS Header. - - Download FITS header from VOSpace, or create a symbolic link - if the header exists in dir_for_links. - - Parameters - ---------- - expnum : int - exposure number - client : vos.Client - VOSpace client - - Returns - ------- - bool - True if successful, False otherwise - - """ - vospace_path = self._params["vospace_path"] - output_dir = self._params["output_dir"] - overwrite = self._params["overwrite"] - dir_for_links = self._params["dir_for_links"] - - source = f"{vospace_path}/{expnum:d}p.fits.fz" - dest = f"{output_dir}/{expnum:d}.txt" - - if exists(dest) and not overwrite: - return True - - # Check if header exists in dir_for_links - if dir_for_links is not None: - link_source = os.path.abspath(f"{dir_for_links}/{expnum:d}.txt") - if exists(link_source): - try: - # Remove existing file/link if overwrite is True - if exists(dest): - os.remove(dest) - # Create symbolic link - os.symlink(link_source, dest) - return True - except Exception as e: - print(f"Could not create symlink from {link_source}: {e}") - # Fall through to download if symlink fails - - # Download from VOSpace atomically: copy to a temp file in the same - # directory, then rename on success. An interrupted transfer leaves - # only the temp file behind, so resume never treats a partial download - # as complete. - tmp_dest = f"{dest}.part" - try: - client.copy(source, tmp_dest, head=True) - os.rename(tmp_dest, dest) - return True - except Exception as e: - print(f"Could not copy {source}: {e}") - if exists(tmp_dest): - os.remove(tmp_dest) - return False - - def run(self, args=None): - """Run. - - Main execution method. - - Parameters - ---------- - args : list, optional - command line arguments - - Returns - ------- - int - exit code (0 for success) - - """ - if args is None: - args = sys.argv[1:] - - # Set parameters from command line - self.set_params_from_command_line(args) - self.update_params() - self.check_params() - - # Get parameters - input_file = self._params["input_file"] - output_dir = self._params["output_dir"] - verbose = self._params["verbose"] - - if verbose: - print(f"Reading CCD list from: {input_file}") - - # Extract unique exposures - exps = self.get_exposures(input_file) - - print(f"Found {len(exps)} unique exposures") - - if verbose: - print(f"Downloading headers to: {output_dir}/") - - # Initialize VOSpace client once - client = vos.Client() - - # Download headers - n_success = 0 - n_failed = 0 - - for i, exp in enumerate(exps): - success = self.get_fits_header(exp, client) - if success: - n_success += 1 - else: - n_failed += 1 - - if verbose and i % 100 == 0: - print(f"{i:6d} / {len(exps):6d}") - - print(f"Downloaded {n_success} headers") - if n_failed > 0: - print(f"Failed to download {n_failed} headers") - - return 0 diff --git a/tests/module/test_coverage.py b/tests/module/test_coverage.py index f7808a375..d16cd518f 100644 --- a/tests/module/test_coverage.py +++ b/tests/module/test_coverage.py @@ -1,56 +1,35 @@ -"""UNIT / PROPERTY TESTS FOR THE COVERAGE-MASK FEATURE. - -Covers the pure and lightly-fixtured logic behind the per-CCD coverage nexp -masks: exposure-number parsing, the CCD-list -> unique-exposure reduction, the -handler's missing-CCD subtraction (pinned against the real ID format), per-CCD -corner extraction from plain and fpack-compressed multi-HDU headers, -``--ccd_list`` filtering, the per-CCD resume path, the builder's per-CCD row -parsing, the accumulated exposure-count (nexp) contract, the RA-wrap and pole -guards, the power-of-two ``nside`` validation, the atomic-download rename, and -the ``-h``/argv handling of the console runners. The heavy paths (real VOSpace -download, production-resolution map building, plotting, multiprocessing) are -exercised end to end by the pipeline, not here. +"""UNIT / PROPERTY TESTS FOR THE COVERAGE-MASK GEOMETRY. + +Covers the pure logic the coverage nexp mask rests on, on both sides of the +``exp_footprint`` -> ``nexp_map`` handover: per-CCD image-shape resolution +from plain and fpack-compressed headers, the sky-corner projection, and then +the accumulated exposure-count (nexp) contract, the RA-wrap and pole guards, +the power-of-two ``nside`` validation and the median smoothing of a finished +map. Everything here takes arrays or headers in memory; the way a record +reaches these functions is ``tests/unit/test_exp_footprint.py``'s subject, and +the heavy paths (production-resolution map building, plotting) are exercised +end to end by the pipeline. """ -import sys -from types import SimpleNamespace - +import healsparse as hsp +import hpgeom as hpg import numpy as np import numpy.testing as npt import pytest from astropy import wcs -from hypothesis import given -from hypothesis import strategies as st +from astropy.io.fits import Header -from shapepipe.utilities import summary -from shapepipe.utilities.ccd_psf_handler import CcdPsfHandler +from shapepipe.utilities.ccd_footprint import _ccd_corners, _image_shape from shapepipe.utilities.coverage_map_builder import ( - CoverageMapBuilder, + build_map, + check_nside, + median_filter, unwrap_ra, ) -from shapepipe.utilities.coverage_plotter import CoveragePlotter -from shapepipe.utilities.field_corners_extractor import ( - FieldCornersExtractor, - _ccd_corners, - _expnum_from_path, - _image_shape, - _parse_header_to_wcs, -) -from shapepipe.utilities.header_downloader import HeaderDownloader - - -@pytest.fixture(autouse=True) -def _run_in_tmp_path(tmp_path, monkeypatch): - """Run each test in its own tmp dir: the runners log to ``log_`` in cwd.""" - monkeypatch.chdir(tmp_path) def _tan_wcs(crval_ra, crval_dec=0.0, nx=2080, ny=4612): - """Build a minimal 2-D TAN WCS with a populated pixel shape. - - ``nx``/``ny`` set ``pixel_shape`` so ``_header_text`` can emit matching - ``NAXIS1``/``NAXIS2`` keywords. - """ + """Build a minimal 2-D TAN WCS with a populated pixel shape.""" w = wcs.WCS(naxis=2) w.wcs.ctype = ["RA---TAN", "DEC--TAN"] w.wcs.crval = [crval_ra, crval_dec] @@ -60,244 +39,61 @@ def _tan_wcs(crval_ra, crval_dec=0.0, nx=2080, ny=4612): return w -def _header_text(wcs_list): - """Render a multi-HDU header text file: a primary HDU then one plain image - HDU per WCS, each carrying ``NAXIS1``/``NAXIS2``. - """ - blocks = ["SIMPLE = T"] - for w in wcs_list: - nx, ny = w.pixel_shape - header = w.to_header() - header["NAXIS"] = 2 - header["NAXIS1"] = nx - header["NAXIS2"] = ny - blocks.append("\n".join(str(card) for card in header.cards)) - return "".join(f"{block}\nEND \n" for block in blocks) +def _plain_header(w): + """A plain image-HDU header: ``NAXIS1``/``NAXIS2`` are the image.""" + nx, ny = w.pixel_shape + header = w.to_header() + header["NAXIS"] = 2 + header["NAXIS1"] = nx + header["NAXIS2"] = ny + return header -def _compressed_header_text(w): - """Render a one-CCD header mimicking an fpack tile-compressed HDU. +def _compressed_header(w): + """A header shaped like an fpack tile-compressed HDU. ``NAXIS1``/``NAXIS2`` describe the compressed binary table (byte width, row - count); the true image dimensions live in ``ZNAXIS1``/``ZNAXIS2``. This is - the shape of headers fetched from 'p.fits.fz' with ``head=True``. + count); the true image dimensions live in ``ZNAXIS1``/``ZNAXIS2``. """ nx, ny = w.pixel_shape - wcs_cards = "\n".join(str(card) for card in w.to_header().cards) - ccd = ( - "XTENSION= 'BINTABLE'\n" - "BITPIX = 8\n" - "NAXIS = 2\n" - "NAXIS1 = 8\n" - f"NAXIS2 = {ny}\n" - "PCOUNT = 1000000\n" - "GCOUNT = 1\n" - "TFIELDS = 1\n" - "ZIMAGE = T\n" - "ZBITPIX = -32\n" - "ZNAXIS = 2\n" - f"ZNAXIS1 = {nx}\n" - f"ZNAXIS2 = {ny}\n" - f"{wcs_cards}" - ) - return f"SIMPLE = T\nEND \n{ccd}\nEND \n" - - -# --------------------------------------------------------------------------- -# _expnum_from_path -# --------------------------------------------------------------------------- - -@pytest.mark.parametrize( - "path, expected", - [ - ("1234567.txt", 1234567), - ("/a/b/2143523.txt", 2143523), - ("headers/0000042.txt", 42), - ], -) -def test_expnum_from_path_extracts_trailing_number(path, expected): - """The trailing ``.txt`` is parsed as the exposure number.""" - assert _expnum_from_path(path) == expected - - -@pytest.mark.parametrize("path", ["no_number.txt", "1234567.fits", "abc.txt"]) -def test_expnum_from_path_raises_without_number(path): - """A filename without a trailing numeric stem raises ``ValueError``.""" - with pytest.raises(ValueError): - _expnum_from_path(path) - - -@given(st.integers(min_value=0, max_value=99999999)) -def test_expnum_from_path_roundtrips(expnum): - """Any exposure number round-trips through the filename convention.""" - assert _expnum_from_path(f"vos_headers/{expnum}.txt") == expnum - - -# --------------------------------------------------------------------------- -# HeaderDownloader.get_exposures -# --------------------------------------------------------------------------- - -def test_get_exposures_reduces_to_unique_exposures(tmp_path): - """A ``-`` CCD list collapses to its unique exposure numbers.""" - ccd_list = tmp_path / "ccds.txt" - ccd_list.write_text("2143523-0\n2143523-5\n2143524-3\n2143524-8\n") - - exps = HeaderDownloader().get_exposures(str(ccd_list)) - - npt.assert_array_equal(exps, np.array([2143523, 2143524])) - - -def test_get_exposures_single_line(tmp_path): - """A one-line CCD list (0-d loadtxt array) still yields one exposure.""" - ccd_list = tmp_path / "ccds.txt" - ccd_list.write_text("2143523-0\n") - - exps = HeaderDownloader().get_exposures(str(ccd_list)) - - npt.assert_array_equal(exps, np.array([2143523])) - - -def test_get_exposures_csv_matches_txt(tmp_path): - """The CSV and text code paths yield the same unique exposures.""" - txt = tmp_path / "ccds.txt" - txt.write_text("2143523-0\n2143523-5\n2143524-3\n") - csv = tmp_path / "ccds.csv" - csv.write_text("CCD\n2143523-0\n2143523-5\n2143524-3\n") - - dl = HeaderDownloader() - - npt.assert_array_equal( - dl.get_exposures(str(txt)), dl.get_exposures(str(csv)) - ) - - -# --------------------------------------------------------------------------- -# HeaderDownloader.get_fits_header — atomic rename -# --------------------------------------------------------------------------- - -def test_get_fits_header_writes_atomically(tmp_path): - """A successful download copies to ``.part`` then renames to the dest.""" - dl = HeaderDownloader() - dl._params["output_dir"] = str(tmp_path) - dl._params["overwrite"] = False - dl._params["dir_for_links"] = None - - dest = tmp_path / "42.txt" - tmp_dest = tmp_path / "42.txt.part" - - def fake_copy(source, target, head=True): - # The copy must land on the temp path, not the final destination. - assert target == str(tmp_dest) - with open(target, "w") as f: - f.write("HEADER") - - client = SimpleNamespace(copy=fake_copy) - - assert dl.get_fits_header(42, client) is True - assert dest.exists() - assert not tmp_dest.exists() - assert dest.read_text() == "HEADER" - - -def test_get_fits_header_failed_copy_leaves_no_dest(tmp_path): - """A failed download leaves no destination file (only, if any, ``.part``).""" - dl = HeaderDownloader() - dl._params["output_dir"] = str(tmp_path) - dl._params["overwrite"] = False - dl._params["dir_for_links"] = None - - def failing_copy(source, target, head=True): - raise RuntimeError("transfer interrupted") - - client = SimpleNamespace(copy=failing_copy) - - assert dl.get_fits_header(42, client) is False - assert not (tmp_path / "42.txt").exists() - assert not (tmp_path / "42.txt.part").exists() - - -# --------------------------------------------------------------------------- -# CcdPsfHandler.get_ccds_with_psf — missing-CCD subtraction -# --------------------------------------------------------------------------- - -def test_get_ccds_with_psf_subtracts_missing(monkeypatch): - """Valid CCDs are all exposure single-HDUs minus the missing set. - - The real ``summary.get_all_shdus`` is used so the cross-component - ``-`` ID format is pinned end to end. - """ - handler = CcdPsfHandler() - - # Two exposures, 3 CCDs each -> 6 candidate CCDs; two are missing. - monkeypatch.setattr(handler, "get_exp", lambda patches: {"100", "200"}) - monkeypatch.setattr( - handler, - "get_exp_shdu_missing", - lambda patches: {"100-1", "200-2"}, - ) - - result = handler.get_ccds_with_psf(["P1"], n_CCD=3) - - # get_all_shdus yields "-" for ccd in range(n_CCD). - assert result == {"100-0", "100-2", "200-0", "200-1"} - # Guard the assumption that the missing IDs share the produced format. - assert set(summary.get_all_shdus({"100"}, 3)) == {"100-0", "100-1", "100-2"} - - -@pytest.mark.parametrize( - ("version", "n_patch"), - [("v1.3", 7), ("v1.4", 7), ("v1.5", 8), ("v1.6", 9)], -) -def test_version_to_patch_count(version, n_patch): - """Each v1.x catalogue version maps to its patch count.""" - handler = CcdPsfHandler() - handler._params["version_cat"] = version - handler.update_params() - assert handler._params["n_patch"] == n_patch - assert len(handler._params["patches"]) == n_patch - - -def test_invalid_version_raises(): - """An unknown catalogue version fails loudly.""" - handler = CcdPsfHandler() - handler._params["version_cat"] = "v9.9" - with pytest.raises(ValueError, match="v9.9"): - handler.update_params() + header = w.to_header() + header["NAXIS"] = 2 + header["NAXIS1"] = 8 + header["NAXIS2"] = ny + header["ZIMAGE"] = True + header["ZNAXIS"] = 2 + header["ZNAXIS1"] = nx + header["ZNAXIS2"] = ny + return header # --------------------------------------------------------------------------- # image-shape resolution (fpack ZNAXIS vs plain NAXIS) # --------------------------------------------------------------------------- -def test_image_shape_prefers_znaxis_for_compressed_header(tmp_path): +def test_image_shape_prefers_znaxis_for_compressed_header(): """A compressed HDU reports ZNAXIS dims, not the binary-table NAXIS.""" - w = _tan_wcs(100.0, 20.0, nx=2080, ny=4612) - path = tmp_path / "1234567.txt" - path.write_text(_compressed_header_text(w)) - - (parsed_w, shape), = _parse_header_to_wcs(str(path)) + header = _compressed_header(_tan_wcs(100.0, 20.0, nx=2080, ny=4612)) # WCS pixel_shape would wrongly report the compressed byte width (8). - assert parsed_w.pixel_shape == (8, 4612) + assert wcs.WCS(header).pixel_shape == (8, 4612) # _image_shape recovers the true image dimensions. - assert shape == (2080, 4612) + assert _image_shape(header) == (2080, 4612) -def test_image_shape_falls_back_to_naxis(tmp_path): - """A plain image HDU (no ZIMAGE) uses NAXIS1/NAXIS2.""" - w = _tan_wcs(100.0, 20.0, nx=2080, ny=4612) - path = tmp_path / "1234567.txt" - path.write_text(_header_text([w])) +def test_image_shape_falls_back_to_naxis(): + """A plain image HDU (no ZIMAGE) uses NAXIS1/NAXIS2. - (_, shape), = _parse_header_to_wcs(str(path)) + This is the live path for the workflow: ``headers-.npy`` carries the + decompressed header astropy hands back for a tile-compressed HDU. + """ + header = _plain_header(_tan_wcs(100.0, 20.0, nx=2080, ny=4612)) - assert shape == (2080, 4612) + assert _image_shape(header) == (2080, 4612) def test_image_shape_raises_without_dimensions(): """A header with no image dimensions raises a clear ``ValueError``.""" - from astropy.io.fits import Header - header = Header() header["CTYPE1"] = "RA---TAN" @@ -305,18 +101,14 @@ def test_image_shape_raises_without_dimensions(): _image_shape(header) -def test_compressed_header_corners_are_full_width(tmp_path): +def test_compressed_header_corners_are_full_width(): """Corners from a compressed header span the true CCD width, not 8 px.""" - w = _tan_wcs(100.0, 20.0, nx=2080, ny=4612) - path = tmp_path / "1234567.txt" - path.write_text(_compressed_header_text(w)) + header = _compressed_header(_tan_wcs(100.0, 20.0, nx=2080, ny=4612)) - (parsed_w, shape), = _parse_header_to_wcs(str(path)) - ra, dec = _ccd_corners(parsed_w, shape) + ra, dec = _ccd_corners(wcs.WCS(header), _image_shape(header)) # RA extent must reflect ~2080 px * 1e-5 deg/px * cos(dec), not 8 px. - ra_extent = max(ra) - min(ra) - assert ra_extent > 0.01 + assert max(ra) - min(ra) > 0.01 # --------------------------------------------------------------------------- @@ -339,257 +131,56 @@ def test_ccd_corners_returns_four_corners_around_centre(): npt.assert_allclose(max(dec) - min(dec), ny * 1e-5, rtol=1e-3) -def test_parse_header_to_wcs_returns_one_wcs_per_extension(tmp_path): - """One (wcs, shape) pair is returned per CCD HDU; the primary is skipped.""" - ccd_wcs = [_tan_wcs(10.0), _tan_wcs(20.0), _tan_wcs(30.0)] - - path = tmp_path / "1234567.txt" - path.write_text(_header_text(ccd_wcs)) - - result = _parse_header_to_wcs(str(path)) - - assert len(result) == len(ccd_wcs) - assert all(shape == (2080, 4612) for _, shape in result) - - # --------------------------------------------------------------------------- -# FieldCornersExtractor.process_single_header — per-CCD rows and filtering +# nside validation # --------------------------------------------------------------------------- -def test_process_single_header_emits_one_row_per_ccd(tmp_path): - """Without a CCD list, every HDU yields a ``-`` row.""" - ccd_wcs = [_tan_wcs(10.0), _tan_wcs(20.0), _tan_wcs(30.0)] - path = tmp_path / "1234567.txt" - path.write_text(_header_text(ccd_wcs)) - - rows = FieldCornersExtractor.process_single_header( - (str(path), False, None) - ) - - assert [r[0] for r in rows] == [ - "1234567-0", - "1234567-1", - "1234567-2", - ] - assert all(len(r[1]) == 4 and len(r[2]) == 4 for r in rows) - - -def test_process_single_header_filters_to_ccd_list(tmp_path): - """A ``valid_ccds`` set keeps only the listed CCDs of the exposure.""" - ccd_wcs = [_tan_wcs(10.0), _tan_wcs(20.0), _tan_wcs(30.0)] - path = tmp_path / "1234567.txt" - path.write_text(_header_text(ccd_wcs)) - - rows = FieldCornersExtractor.process_single_header( - (str(path), False, {"1234567-0", "1234567-2"}) - ) - - assert [r[0] for r in rows] == ["1234567-0", "1234567-2"] - - -def test_load_ccd_list_reads_ids(tmp_path): - """The CCD list loader returns the stripped, non-blank IDs as a set.""" - path = tmp_path / "ccds.txt" - path.write_text("100-0\n100-1\n\n200-5\n") - - assert FieldCornersExtractor.load_ccd_list(str(path)) == { - "100-0", - "100-1", - "200-5", - } - - -def test_run_extract_writes_per_ccd_rows(tmp_path): - """End to end: run() writes one per-CCD row filtered by the CCD list.""" - header_dir = tmp_path / "headers" - header_dir.mkdir() - ccd_wcs = [_tan_wcs(10.0), _tan_wcs(20.0), _tan_wcs(30.0)] - (header_dir / "1234567.txt").write_text(_header_text(ccd_wcs)) - - ccd_list = tmp_path / "ccds.txt" - ccd_list.write_text("1234567-0\n1234567-2\n") - - out = tmp_path / "corners.txt" - - extractor = FieldCornersExtractor() - extractor.run( - args=[ - "-i", str(header_dir), - "-l", str(ccd_list), - "-o", str(out), - ] - ) - - lines = out.read_text().splitlines() - assert len(lines) == 2 - ids = [line.split()[0] for line in lines] - assert ids == ["1234567-0", "1234567-2"] - # Each row: 1 ID + 4 RA + 4 Dec = 9 columns. - assert all(len(line.split()) == 9 for line in lines) - - -# --------------------------------------------------------------------------- -# FieldCornersExtractor resume path (per-CCD done-set) -# --------------------------------------------------------------------------- - -def _write_headers(header_dir, expnums, n_ccd=3): - """Write one plain multi-HDU header per exposure into ``header_dir``.""" - header_dir.mkdir(exist_ok=True) - for j, expnum in enumerate(expnums): - ccd_wcs = [_tan_wcs(10.0 + j + i) for i in range(n_ccd)] - (header_dir / f"{expnum}.txt").write_text(_header_text(ccd_wcs)) - - -def test_get_done_ccds_reads_present_ids(tmp_path): - """The done-set is the exact set of CCD IDs already in the output.""" - out = tmp_path / "corners.txt" - out.write_text( - "1234567-0 10 10 10 10 20 20 20 20\n" - "1234567-2 10 10 10 10 20 20 20 20\n" - ) - - extractor = FieldCornersExtractor() - extractor._params["output_file"] = str(out) - - assert extractor.get_done_ccds() == {"1234567-0", "1234567-2"} - - -def test_resume_adds_new_exposure_without_duplicating(tmp_path): - """Resume appends a new exposure and leaves existing rows untouched.""" - header_dir = tmp_path / "headers" - _write_headers(header_dir, [1000001, 1000002]) - - out = tmp_path / "corners.txt" - # Pre-populate with exposure 1's three CCD rows. - pre = FieldCornersExtractor().process_single_header( - (str(header_dir / "1000001.txt"), False, None) - ) - with open(out, "w") as f: - for ccd_id, ra, dec in pre: - f.write(f"{ccd_id} " + " ".join(f"{v:.5f}" for v in ra + dec) + "\n") - - FieldCornersExtractor().run( - args=["-i", str(header_dir), "-o", str(out), "-r"] - ) - - ids = [line.split()[0] for line in out.read_text().splitlines()] - # No duplicate exposure-1 rows; exposure 2's three CCDs added. - assert ids.count("1000001-0") == 1 - assert sorted(ids) == [ - "1000001-0", "1000001-1", "1000001-2", - "1000002-0", "1000002-1", "1000002-2", - ] - - -def test_resume_completes_partial_exposure(tmp_path): - """An exposure interrupted mid-write is completed, not skipped.""" - header_dir = tmp_path / "headers" - _write_headers(header_dir, [1000001]) - - out = tmp_path / "corners.txt" - # Simulate an interrupt: only the first of exposure 1's CCDs was written. - pre = FieldCornersExtractor().process_single_header( - (str(header_dir / "1000001.txt"), False, None) - ) - ccd_id, ra, dec = pre[0] - with open(out, "w") as f: - f.write(f"{ccd_id} " + " ".join(f"{v:.5f}" for v in ra + dec) + "\n") - - FieldCornersExtractor().run( - args=["-i", str(header_dir), "-o", str(out), "-r"] - ) - - ids = [line.split()[0] for line in out.read_text().splitlines()] - # The missing CCDs are filled in; the present one is not duplicated. - assert sorted(ids) == ["1000001-0", "1000001-1", "1000001-2"] - assert ids.count("1000001-0") == 1 - - -# --------------------------------------------------------------------------- -# CoverageMapBuilder.check_params (nside validation) -# --------------------------------------------------------------------------- - -def test_check_params_accepts_power_of_two_nside(tmp_path): +def test_check_nside_accepts_powers_of_two(): """Powers of two for both nside values pass validation.""" - infile = tmp_path / "corners.txt" - infile.write_text("100-0 0.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0\n") - - builder = CoverageMapBuilder() - builder._params["input_file"] = str(infile) - builder._params["nside_coverage"] = 32 - builder._params["nside"] = 2048 - - builder.check_params() + check_nside(32, 2048) @pytest.mark.parametrize( - "nside_coverage, nside", [(33, 2048), (32, 100), (32, 3000)] + "nside_coverage, nside", [(33, 2048), (32, 100), (32, 3000), (0, 2048)] ) -def test_check_params_rejects_non_power_of_two_nside( - tmp_path, nside_coverage, nside -): +def test_check_nside_rejects_non_power_of_two(nside_coverage, nside): """A non-power-of-two nside raises ``ValueError``.""" - infile = tmp_path / "corners.txt" - infile.write_text("100-0 0.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0\n") - - builder = CoverageMapBuilder() - builder._params["input_file"] = str(infile) - builder._params["nside_coverage"] = nside_coverage - builder._params["nside"] = nside - with pytest.raises(ValueError): - builder.check_params() + check_nside(nside_coverage, nside) -def test_check_params_missing_input_file_raises(tmp_path): - """A missing input file raises ``FileNotFoundError``.""" - builder = CoverageMapBuilder() - builder._params["input_file"] = str(tmp_path / "does_not_exist.txt") - - with pytest.raises(FileNotFoundError): - builder.check_params() +def test_build_map_validates_nside_before_stamping(): + """``build_map`` refuses a bad nside before it reaches healsparse.""" + with pytest.raises(ValueError, match="nside"): + build_map(["100-0"], [0.0, 0.1, 0.1, 0.0], [0.0, 0.0, 0.1, 0.1], + 32, 3000) # --------------------------------------------------------------------------- -# CoverageMapBuilder — parsing, nexp contract, RA-wrap and pole guards +# build_map — nexp contract, RA-wrap and pole guards # --------------------------------------------------------------------------- -def _read_map(path): - import healsparse as hsp - - return hsp.HealSparseMap.read(str(path)) - - -def test_build_map_single_row(tmp_path): - """A one-row corners file exercises the atleast_1d/2d parsing guards.""" - infile = tmp_path / "corners.txt" - infile.write_text("100-0 9.9 10.1 10.1 9.9 0.0 0.0 0.2 0.2\n") - out = tmp_path / "cov.hsp" - - CoverageMapBuilder().run( - args=["-i", str(infile), "-o", str(out), "-c", "32", "-n", "1024"] +def test_build_map_single_footprint(): + """One footprint given as flat corner lists exercises the shape guards.""" + m = build_map( + ["100-0"], [9.9, 10.1, 10.1, 9.9], [0.0, 0.0, 0.2, 0.2], 32, 1024 ) - m = _read_map(out) assert 0 < len(m.valid_pixels) < 200 assert m[m.valid_pixels].max() == 1 -def test_build_map_nexp_counts_overlapping_exposures(tmp_path): - """Two overlapping CCDs from different exposures give value 2 in overlap.""" +def test_build_map_nexp_counts_overlapping_exposures(): + """Two overlapping CCDs from two exposures give value 2 in the overlap.""" # Two 0.4x0.4 deg CCDs offset by 0.2 deg in RA -> a central overlap strip. - infile = tmp_path / "corners.txt" - infile.write_text( - "100-0 9.8 10.2 10.2 9.8 19.8 19.8 20.2 20.2\n" - "200-0 10.0 10.4 10.4 10.0 19.8 19.8 20.2 20.2\n" - ) - out = tmp_path / "cov.hsp" - - CoverageMapBuilder().run( - args=["-i", str(infile), "-o", str(out), "-c", "32", "-n", "1024"] + m = build_map( + ["100-0", "200-0"], + [[9.8, 10.2, 10.2, 9.8], [10.0, 10.4, 10.4, 10.0]], + [[19.8, 19.8, 20.2, 20.2], [19.8, 19.8, 20.2, 20.2]], + 32, + 1024, ) - m = _read_map(out) values = m[m.valid_pixels] # The overlap is covered by both exposures (value 2); the union edges by # one (value 1). Both must be present; nothing exceeds 2. @@ -613,30 +204,19 @@ def test_unwrap_ra_leaves_normal_polygon_unchanged(): ) -def test_build_map_invariant_under_ra_shift(tmp_path): +def test_build_map_invariant_under_ra_shift(): """A seam CCD's footprint matches an identical CCD shifted +10 deg in RA. - Rotating both the seam CCD (via +360/unwrap) and a reference CCD onto the - same RA and comparing pixel counts fails if the seam polygon were filling + Comparing the seam polygon's pixel count against a reference polygon of the + same size at the same declination fails if the seam polygon were filling the ~360 deg complement. This is the real RA-wrap regression guard. """ - # Seam CCD straddling RA=0, and the same CCD translated to RA~10. - seam = tmp_path / "seam.txt" - seam.write_text("100-0 359.9 0.1 0.1 359.9 20.0 20.0 20.2 20.2\n") - ref = tmp_path / "ref.txt" - ref.write_text("200-0 9.9 10.1 10.1 9.9 20.0 20.0 20.2 20.2\n") - - m_seam = tmp_path / "seam.hsp" - m_ref = tmp_path / "ref.hsp" - CoverageMapBuilder().run( - args=["-i", str(seam), "-o", str(m_seam), "-c", "32", "-n", "1024"] - ) - CoverageMapBuilder().run( - args=["-i", str(ref), "-o", str(m_ref), "-c", "32", "-n", "1024"] - ) + dec = [20.0, 20.0, 20.2, 20.2] + m_seam = build_map(["100-0"], [359.9, 0.1, 0.1, 359.9], dec, 32, 1024) + m_ref = build_map(["200-0"], [9.9, 10.1, 10.1, 9.9], dec, 32, 1024) - n_seam = len(_read_map(m_seam).valid_pixels) - n_ref = len(_read_map(m_ref).valid_pixels) + n_seam = len(m_seam.valid_pixels) + n_ref = len(m_ref.valid_pixels) # Same-size footprints at the same declination: pixel counts agree to # within a few boundary pixels, and are nowhere near a hemisphere. @@ -644,17 +224,14 @@ def test_build_map_invariant_under_ra_shift(tmp_path): assert 0 < n_seam < 200 -def test_build_map_pole_guard_skips_polygon(tmp_path, capsys): +def test_build_map_pole_guard_skips_polygon(capsys): """A polygon with |dec| near 90 deg is skipped with a warning.""" - infile = tmp_path / "corners.txt" - infile.write_text( - "100-0 10.0 10.2 10.2 10.0 89.5 89.5 89.7 89.7\n" - "200-0 10.0 10.2 10.2 10.0 20.0 20.0 20.2 20.2\n" - ) - out = tmp_path / "cov.hsp" - - CoverageMapBuilder().run( - args=["-i", str(infile), "-o", str(out), "-c", "32", "-n", "1024"] + build_map( + ["100-0", "200-0"], + [[10.0, 10.2, 10.2, 10.0], [10.0, 10.2, 10.2, 10.0]], + [[89.5, 89.5, 89.7, 89.7], [20.0, 20.0, 20.2, 20.2]], + 32, + 1024, ) captured = capsys.readouterr() @@ -663,42 +240,17 @@ def test_build_map_pole_guard_skips_polygon(tmp_path, capsys): # --------------------------------------------------------------------------- -# console-runner argv handling (regression guard for the -h entry points) +# median smoothing (offline; the nexp_map rule writes the raw count) # --------------------------------------------------------------------------- -def test_run_help_flag_exits_cleanly(monkeypatch): - """``run()`` with no args reads ``sys.argv[1:]`` so ``-h`` exits 0. +def test_median_filter_fills_a_lone_hole(): + """A lone low pixel in a covered disc is pulled up to its neighbours.""" + nside = 1024 + m = hsp.HealSparseMap.make_empty(32, nside, np.uint16) + m[hpg.query_circle(nside, 10.0, 20.0, 0.1)] = 2 - Guards the regression where the runners parsed the full ``sys.argv`` - (including ``argv[0]``), which made ``-h`` collide with the integer - ``-c`` option and exit non-zero. - """ - monkeypatch.setattr(sys, "argv", ["extract_field_corners", "-h"]) - - with pytest.raises(SystemExit) as excinfo: - FieldCornersExtractor().run() - - assert excinfo.value.code == 0 - - -@pytest.mark.parametrize( - "runner, prog", - [ - (CcdPsfHandler, "get_ccds_with_psf"), - (HeaderDownloader, "download_headers"), - (FieldCornersExtractor, "extract_field_corners"), - (CoverageMapBuilder, "build_coverage_map"), - (CoveragePlotter, "plot_coverage_map"), - ], -) -def test_command_logged_under_program_name(tmp_path, runner, prog): - """Each runner logs its command line to ``log_``. - - ``args`` excludes the program name while ``log_command`` names the file - after ``argv[0]``; guards against the log being named after the first - flag (``log_-o``). - """ - runner().set_params_from_command_line(["-o", "out"]) + centre = hpg.angle_to_pixel(nside, 10.0, 20.0) + m[[centre]] = 1 + assert m[centre] == 1 - assert (tmp_path / f"log_{prog}").read_text() == f"{prog} -o out\n" - assert sorted(f.name for f in tmp_path.iterdir()) == [f"log_{prog}"] + assert median_filter(m)[centre] == 2 diff --git a/tests/module/test_nexp_map.py b/tests/module/test_nexp_map.py new file mode 100644 index 000000000..6ee419988 --- /dev/null +++ b/tests/module/test_nexp_map.py @@ -0,0 +1,207 @@ +"""The workflow's exposure-count map: exposure footprint records -> a HealSparse nexp map. + +``tests/module/test_coverage.py`` pins the nexp contract on ``build_map`` +directly. This module pins the SAME contract through the path the Snakemake +workflow takes — one ``exp_footprint.json`` per exposure on the products root, +globbed and fed to ``build_map`` as arrays — because that path has two +properties neither visible in the map nor testable on arrays: + + * it is CAMPAIGN-CUMULATIVE by construction. The script globs every record + under ``/exp/*/*/manifests/``, sharded, including exposures + whose scratch stores were reclaimed. A regression that fed it only the + declared rule inputs would build a map of one batch and look fine. + * the records are RAW SKY. Unwrapping across RA=0 happens once, inside + ``build_map``; a record written pre-unwrapped, or unwrapped twice, silently + moves a footprint by 360 degrees. + +Two exposures, one CCD each, offset so they overlap — the same fixture geometry +as the ``build_map`` test, so the two routes are comparable by eye. +""" + +import hashlib +import importlib.util +import json +import sys +from pathlib import Path + +import pytest + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPT = REPO_ROOT / "workflow" / "scripts" / "nexp_map.py" + + +def _load(): + """Import the script by path — ``workflow/scripts`` is not a package.""" + assert SCRIPT.exists(), f"{SCRIPT} not found; the rule calls it by path" + spec = importlib.util.spec_from_file_location("_nexp_map", SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +nexp_map = _load() + + +def write_footprint(products, exp, ccds): + """One exposure's record, at the sharded path exp_footprint writes to.""" + path = (products / "exp" / exp[:2] / exp / "manifests" + / "exp_footprint.json") + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text(json.dumps({ + "stage": "exp_footprint", "level": "exp", "unit": exp, + "status": "complete", + "n_ccd_headers": len(ccds), "n_valid_psf": len(ccds), + "ccds": [{"id": f"{exp}-{i}", "ra": ra, "dec": dec} + for i, (ra, dec) in enumerate(ccds)], + "ccds_no_psf": [], + }, indent=2, sort_keys=True)) + return path + + +def box(ra_min, ra_max, dec_min, dec_max): + """A rectangular CCD footprint, corners counterclockwise from bottom-left.""" + return ([ra_min, ra_max, ra_max, ra_min], + [dec_min, dec_min, dec_max, dec_max]) + + +def run(products, tmp_path, monkeypatch, nside="1024"): + out = tmp_path / "nexp_map" / "nexp_map_fixture-run.hsp" + manifest = tmp_path / "nexp_map" / "nexp_map_fixture-run.json" + monkeypatch.setattr(sys, "argv", [ + "nexp_map.py", + "--products-dir", str(products), + "--out", str(out), "--manifest", str(manifest), + "--nside-coverage", "32", "--nside", nside]) + nexp_map.main() + import healsparse as hsp + return hsp.HealSparseMap.read(str(out)), json.loads(manifest.read_text()) + + +def test_two_exposures_accumulate_to_nexp(tmp_path, monkeypatch): + """Overlapping CCDs from two exposures give 2 in the overlap, 1 outside. + + The exposures live in DIFFERENT shard directories, which is the campaign + layout: a glob that assumed one shard would find only one of them. + """ + products = tmp_path / "products" + write_footprint(products, "1000001", [box(9.8, 10.2, 19.8, 20.2)]) + write_footprint(products, "2000001", [box(10.0, 10.4, 19.8, 20.2)]) + + m, manifest = run(products, tmp_path, monkeypatch) + + values = m[m.valid_pixels] + assert values.max() == 2 + assert values.min() == 1 + assert (values == 2).sum() > 0 + assert (values == 1).sum() > 0 + + assert manifest["n_exposures"] == 2 + assert manifest["exposures"] == ["1000001", "2000001"] + assert manifest["n_ccds"] == 2 + assert manifest["nside_coverage"] == 32 + assert manifest["map"] == str( + tmp_path / "nexp_map" / "nexp_map_fixture-run.hsp") + assert manifest["n_coverage_pixels"] == int(m.coverage_mask.sum()) + + +def test_every_record_on_the_root_is_used(tmp_path, monkeypatch): + """The map is campaign-cumulative: nothing selects a subset of the records. + + A third exposure appears on the products root with no involvement from any + caller — the state a reclaimed or out-of-scope exposure is in — and it must + still be in the map. + """ + products = tmp_path / "products" + write_footprint(products, "1000001", [box(9.8, 10.2, 19.8, 20.2)]) + write_footprint(products, "2000001", [box(10.0, 10.4, 19.8, 20.2)]) + write_footprint(products, "3000001", [box(40.0, 40.4, 19.8, 20.2)]) + + m, manifest = run(products, tmp_path, monkeypatch) + + assert manifest["exposures"] == ["1000001", "2000001", "3000001"] + # The far-away exposure is real sky in the map, not just a row in the record. + assert m.get_values_pos(40.2, 20.0, lonlat=True) == 1 + + +def test_seam_record_is_unwrapped_by_the_builder(tmp_path, monkeypatch): + """A raw-sky record across RA=0 lands where it belongs, not 360 deg away.""" + products = tmp_path / "products" + write_footprint(products, "1000001", + [([359.8, 0.2, 0.2, 359.8], [19.8, 19.8, 20.2, 20.2])]) + + m, _ = run(products, tmp_path, monkeypatch) + + assert m.get_values_pos(0.0, 20.0, lonlat=True) == 1 + assert m.get_values_pos(359.9, 20.0, lonlat=True) == 1 + # Nothing was stamped on the far side of the sky. + assert m.get_values_pos(180.0, 20.0, lonlat=True) == 0 + + +def test_no_records_is_a_loud_failure(tmp_path, monkeypatch): + """An empty products root must not write an empty map and exit green.""" + products = tmp_path / "products" + products.mkdir() + with pytest.raises(SystemExit) as exc: + run(products, tmp_path, monkeypatch) + assert "nothing to build a map from" in str(exc.value) + + +def test_records_naming_no_ccd_is_a_loud_failure(tmp_path, monkeypatch): + """Records that name no CCD are the quieter empty map, and equally fatal. + + Every exposure on the root having lost every CCD is a broken PSF stage, not + a survey with no coverage — but the .hsp it would write is valid and + plausible, and its consumer would mask everything. + """ + products = tmp_path / "products" + write_footprint(products, "1000001", []) + write_footprint(products, "2000001", []) + + with pytest.raises(SystemExit) as exc: + run(products, tmp_path, monkeypatch) + assert "not one names a CCD with a PSF model" in str(exc.value) + + +def test_manifest_records_the_map_digest(tmp_path, monkeypatch): + """The manifest names the map generation it describes, by content.""" + products = tmp_path / "products" + write_footprint(products, "1000001", [box(9.8, 10.2, 19.8, 20.2)]) + + _, manifest = run(products, tmp_path, monkeypatch) + + out = Path(manifest["map"]) + assert manifest["map_sha256"] == hashlib.sha256(out.read_bytes()).hexdigest() + + +def test_interrupted_write_keeps_the_published_map(tmp_path, monkeypatch): + """A job killed after writing the new map, before publishing it, leaves the + previous map and its manifest exactly as they were.""" + products = tmp_path / "products" + write_footprint(products, "1000001", [box(9.8, 10.2, 19.8, 20.2)]) + _, manifest = run(products, tmp_path, monkeypatch) + out = Path(manifest["map"]) + manifest_path = tmp_path / "nexp_map" / "nexp_map_fixture-run.json" + old_map, old_manifest = out.read_bytes(), manifest_path.read_bytes() + + build_map = nexp_map.build_map + + def killed_after_write(*args, **kwargs): + hsp_map = build_map(*args, **kwargs) + write = hsp_map.write + + def write_then_die(path, **kw): + write(path, **kw) + raise RuntimeError("killed after the map was written") + + hsp_map.write = write_then_die + return hsp_map + + monkeypatch.setattr(nexp_map, "build_map", killed_after_write) + write_footprint(products, "2000001", [box(40.0, 40.4, 19.8, 20.2)]) + with pytest.raises(RuntimeError, match="killed"): + run(products, tmp_path, monkeypatch) + + assert out.read_bytes() == old_map + assert manifest_path.read_bytes() == old_manifest + assert sorted(p.name for p in out.parent.iterdir()) == [ + "nexp_map_fixture-run.hsp", "nexp_map_fixture-run.json"] diff --git a/tests/unit/test_campaign_lineage.py b/tests/unit/test_campaign_lineage.py index 492ec0ec8..cbf076e11 100644 --- a/tests/unit/test_campaign_lineage.py +++ b/tests/unit/test_campaign_lineage.py @@ -39,8 +39,18 @@ "prod_exp_dir": (("2605805",), False), "prod_exp_manifest": (("2605805", "exp_persist"), False), "prod_exp_tar": (("2605805",), False), + "defect_map": ((), True), + "defect_map_sidecar": ((), True), + "nexp_map": ((), True), + "nexp_map_manifest": ((), True), } PRODUCT_TEMPLATES = ("PROD_TILE_DIR", "PROD_EXP_DIR") +EXPOSURE_MAP_PATHS = { + "defect_map": f"{PRODUCTS}/defect_map/defect_map_{CAMPAIGN}.hsp", + "defect_map_sidecar": f"{PRODUCTS}/defect_map/defect_map_{CAMPAIGN}.json", + "nexp_map": f"{PRODUCTS}/nexp_map/nexp_map_{CAMPAIGN}.hsp", + "nexp_map_manifest": f"{PRODUCTS}/nexp_map/nexp_map_{CAMPAIGN}.json", +} def _code_lines(path): @@ -128,6 +138,12 @@ def test_product_templates_are_rooted_in_products_dir(helpers, name): assert str(helpers[name]).startswith(PRODUCTS + "/"), helpers[name] +@pytest.mark.parametrize("name,expected", EXPOSURE_MAP_PATHS.items()) +def test_exposure_map_paths_use_products_dir_and_run(helpers, name, expected): + """Each exposure-level map carries the run name; its record stays beside it.""" + assert str(helpers[name]()) == expected + + def _machine_outputs(): config = yaml.safe_load((WORKFLOW / "config.yaml").read_text()) for machine, entry in (config.get("machines") or {}).items(): diff --git a/tests/unit/test_defect_map_contains_flags.py b/tests/unit/test_defect_map_contains_flags.py new file mode 100644 index 000000000..ecc7c7c25 --- /dev/null +++ b/tests/unit/test_defect_map_contains_flags.py @@ -0,0 +1,161 @@ +"""An exposure's defect fragment contains its flag image. + +``rasterize_ccd`` (``workflow/scripts/defect_map_exp.py``) turns one CCD's flag +image and WCS into the healpix pixels its flags touch. The map it feeds is +CONSERVATIVE by contract (``defect-fragment-contains-flags``): a flagged CCD +pixel that lands in an unmasked healpix pixel is a hole in the footprint that +nothing downstream can see. So the invariant here is containment, checked +against the sky geometry itself and not against the function's own sampling: + + * every flagged pixel's centre AND its four corners — the vertices that bound + its footprint — lie in masked healpix pixels, with the sky positions taken + through astropy's 0-based ``pixel_to_world_values``, independent of the + 1-based ``all_pix2world`` path the rasterizer uses; + * a clean flag image yields no pixels at all. + +The fixture carries the three shapes a MegaCam flag image has: a ONE-PIXEL bad +column (the geometry that centre sampling erases), a saturated blob, and +isolated pixels, on the first and last row and column and on a lattice of hot +pixels spaced wider than a healpix pixel. The lattice is what makes the test +sharp: a lone pixel straddling a healpix boundary has no flagged neighbour +whose centre already masks the far side, so centre-only sampling misses it. The WCS is a rotated TAN +at 0.187 arcsec/pixel, MegaCam's scale, so a healpix pixel at the ladder's nside +covers ~74 CCD pixels, as on the sky. ``CHUNK`` is shrunk so the batching runs +several batches. The oversample is read from ``workflow/config.yaml``, the +value a campaign runs with. + +Needs healpy, healsparse and astropy, so it runs inside the container and +skips outside. +""" + +import importlib.util +import sys +from pathlib import Path + +import numpy as np +import pytest + + +pytestmark = [pytest.mark.unions, pytest.mark.decision("masking.defect_map_from_flags")] + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPTS = REPO_ROOT / "workflow" / "scripts" +SCRIPT = SCRIPTS / "defect_map_exp.py" +CONFIG = REPO_ROOT / "workflow" / "config.yaml" + +NY, NX = 320, 240 # a few hundred pixels, not a 2048 x 4612 chip +SCALE_DEG = 0.187 / 3600.0 # MegaCam +ROTATION_DEG = 23.0 # off-axis, so healpix boundaries cut obliquely +CRVAL = (150.3, 31.7) + + +def _load(): + """Import the rule's script by path — scripts/ is not a package.""" + assert SCRIPT.exists(), f"{SCRIPT} not found; the rule calls it by path" + spec = importlib.util.spec_from_file_location("_defect_map_exp", SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +@pytest.fixture(scope="module") +def raster(): + pytest.importorskip("healpy") + pytest.importorskip("healsparse") + pytest.importorskip("astropy") + return _load() + + +@pytest.fixture(scope="module") +def ladder(): + """``(nside, oversample)`` as the campaign config sets them.""" + yaml = pytest.importorskip("yaml") + block = yaml.safe_load(CONFIG.read_text())["exposure_maps"] + return int(block["nside"]), int(block["defect"]["oversample"]) + + +def _wcs(): + from astropy.wcs import WCS + + theta = np.deg2rad(ROTATION_DEG) + wcs = WCS(naxis=2) + wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] + wcs.wcs.crval = list(CRVAL) + wcs.wcs.crpix = [NX / 2 + 0.5, NY / 2 + 0.5] + wcs.wcs.cd = SCALE_DEG * np.array([[-np.cos(theta), np.sin(theta)], + [np.sin(theta), np.cos(theta)]]) + return wcs + + +def _flags(): + flags = np.zeros((NY, NX), dtype=np.int16) + flags[:, 101] = 1 # one-pixel bad column + yy, xx = np.mgrid[:NY, :NX] + flags[(yy - 210) ** 2 + (xx - 60) ** 2 <= 7 ** 2] = 2 # saturated blob + for row, col in [(0, 0), (0, NX - 1), (NY - 1, 0), (NY - 1, NX - 1), + (NY - 1, 150), (40, NX - 1), (77, 33)]: + flags[row, col] = 8 # isolated pixels + # Hot pixels spaced wider than a healpix pixel (11 x 0.187" > 1.61"), so + # one straddling a healpix boundary has no flagged neighbour whose centre + # already masks the far side. This is what makes centre sampling fail. + flags[5::11, 3::11] = 8 + return flags + + +def _write(tmp_path: Path, flags): + """The two files ``rasterize_ccd`` reads: the flag split and the image + split whose header alone carries the WCS.""" + from astropy.io import fits + + flag_path = tmp_path / "flag-2079612-0.fits" + image_path = tmp_path / "image-2079612-0.fits" + fits.PrimaryHDU(data=flags).writeto(flag_path) + fits.PrimaryHDU(header=_wcs().to_header()).writeto(image_path) + return flag_path, image_path + + +def _rasterize(raster, ladder, tmp_path, flags, monkeypatch): + nside, oversample = ladder + monkeypatch.setattr(raster, "CHUNK", 97) + off_x, off_y = raster.offsets(oversample) + flag_path, image_path = _write(tmp_path, flags) + return raster.rasterize_ccd(flag_path, image_path, nside, off_x, off_y) + + +def _vertex_pixels(nside, flags): + """``(healpix id per vertex, (row, col) per vertex)`` for the centre and the + four corners of every flagged pixel, through astropy's 0-based convention.""" + import healpy as hp + + rows, cols = np.nonzero(flags) + offsets = [(0.0, 0.0), (-0.5, -0.5), (-0.5, 0.5), (0.5, -0.5), (0.5, 0.5)] + dx = np.array([o[0] for o in offsets]) + dy = np.array([o[1] for o in offsets]) + x = (cols[:, None] + dx[None, :]).ravel() + y = (rows[:, None] + dy[None, :]).ravel() + ra, dec = _wcs().pixel_to_world_values(x, y) + ids = hp.ang2pix(nside, ra, dec, lonlat=True, nest=True) + where = np.repeat(np.stack([rows, cols], axis=1), len(offsets), axis=0) + return ids, where + + +def test_fragment_contains_every_flagged_pixel(raster, ladder, tmp_path, + monkeypatch): + """Contract defect-fragment-contains-flags: no flagged vertex lands unmasked.""" + nside, _ = ladder + flags = _flags() + masked = _rasterize(raster, ladder, tmp_path, flags, monkeypatch) + + ids, where = _vertex_pixels(nside, flags) + missed = ~np.isin(ids, masked) + assert not missed.any(), ( + f"{missed.sum()} of {missed.size} flagged-pixel vertices fall in " + f"unmasked healpix pixels, at (row, col) " + f"{sorted(set(map(tuple, where[missed].tolist())))[:10]}") + + +def test_clean_flag_image_masks_nothing(raster, ladder, tmp_path, + monkeypatch): + masked = _rasterize(raster, ladder, tmp_path, + np.zeros((NY, NX), dtype=np.int16), monkeypatch) + assert masked.size == 0 diff --git a/tests/unit/test_defect_map_reconcile.py b/tests/unit/test_defect_map_reconcile.py new file mode 100644 index 000000000..80cafc5a1 --- /dev/null +++ b/tests/unit/test_defect_map_reconcile.py @@ -0,0 +1,355 @@ +"""``merge_defect_map`` agrees the campaign's defect map, or rebuilds it. + +The union of the per-exposure defect fragments (CosmoStat/shapepipe#878) is the +one campaign product whose reconciliation is ASYMMETRIC, and that asymmetry is +the reason this file exists. A union cannot be un-OR-ed: two exposures both set +a pixel and nothing in the map records which. So + + * a NEW fragment is OR-ed in on the spot — the cheap, common path; + * a fragment that LEFT the campaign, or one that CHANGED on disk, forces a + REBUILD from every fragment; + * neither leaves the map UNTOUCHED, mtime included, because mtime is a + rerun trigger and an unconditional rewrite makes every invocation look + like a change. + +Those are the branches of ``reconcile_plan`` whose failure mode is silent: an +edit that stopped treating a removal as a rebuild leaves the map carrying bits +from exposures the campaign no longer has, and nothing downstream — nothing in +this workflow reads the map at all — would ever notice. Hence the pins here. + +The SIDECAR is pinned alongside, for the fifth case the plan cannot see: the +record carries two fields about the CAMPAIGN (how many exposures it has, which +of them have no fragment) that can move while the map cannot. The docstring +promises a short map says so on disk; that only holds if a no-op still refreshes +the record. + +EVERY CASE DRIVES ``main``, the rule's real entry point, over a one-tile index, +and asserts on what it prints and writes. A test that planned and applied by +itself would stay green with ``main`` broken. + +Fragments are made with ``HealSparseMap.make_empty`` and a handful of pixel ids, +each with the manifest ``exp_defect_map`` would write beside it. One end-to-end +case rasterizes a single flagged pixel through ``defect_map_exp`` itself. Needs +healsparse, healpy and astropy, so it runs inside the container and skips +outside. +""" + +import importlib.util +import json +import sys +from pathlib import Path + +import pytest + + +pytestmark = pytest.mark.unions + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPTS = REPO_ROOT / "workflow" / "scripts" +SCRIPT = SCRIPTS / "merge_defect_map.py" + +NSIDE = 4096 +NSIDE_COV = 32 + + +def _load(): + """Import the rule's script by path — scripts/ is not a package.""" + assert SCRIPT.exists(), f"{SCRIPT} not found; the rule calls it by path" + sys.path.insert(0, str(SCRIPTS)) + try: + spec = importlib.util.spec_from_file_location("_merge_defect_map", + SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + finally: + sys.path.remove(str(SCRIPTS)) + return module + + +@pytest.fixture(scope="module") +def merge(): + pytest.importorskip("healsparse") + pytest.importorskip("healpy") # defect_map_exp, where the digest lives + pytest.importorskip("astropy") + return _load() + + +def _fragment(merge, root: Path, exp: str, pixels, nside=NSIDE) -> Path: + """Write one exposure's fragment where ``fragment_path`` expects it, and + its manifest where ``manifest_path`` does, as ``exp_defect_map`` would.""" + import hashlib + import numpy as np + import healsparse as hsp + + path = merge.fragment_path(root, exp) + path.parent.mkdir(parents=True, exist_ok=True) + frag = hsp.HealSparseMap.make_empty(NSIDE_COV, nside, np.bool_, + bit_packed=True) + frag[np.asarray(pixels, dtype=np.int64)] = True + frag.write(str(path), clobber=True) + manifest = merge.manifest_path(root, exp) + manifest.parent.mkdir(parents=True, exist_ok=True) + manifest.write_text(json.dumps( + {"map": str(path), + "sha256": hashlib.sha256(path.read_bytes()).hexdigest()})) + return path + + +def _valid(path: Path): + import healsparse as hsp + return set(int(p) for p in hsp.HealSparseMap.read(str(path)).valid_pixels) + +def _main(merge, root: Path, exps, monkeypatch, missing=(), nside=NSIDE): + """Run the rule's real entry point over a one-tile campaign of ``exps`` + (plus ``missing``, indexed but with no fragment); returns its stdout.""" + import contextlib + import io + import sqlite3 + + tile_list, index_db = root / "tiles.txt", root / "index.sqlite" + tile_list.write_text("000.000\n") + index_db.unlink(missing_ok=True) + with sqlite3.connect(index_db) as con: + con.execute("CREATE TABLE tiles(tile_id TEXT PRIMARY KEY, ra_dir TEXT, " + "n_exp INTEGER)") + con.execute("INSERT INTO tiles VALUES ('000.000', '000', 1)") + con.execute("CREATE TABLE tile_exposures(tile_id TEXT, exp_id TEXT)") + con.executemany("INSERT INTO tile_exposures VALUES ('000.000', ?)", + [(e,) for e in [*exps, *missing]]) + monkeypatch.setattr(sys, "argv", [ + str(SCRIPT), "--products-dir", str(root), + "--tile-list", str(tile_list), "--index-db", str(index_db), + "--output", str(root / "defect_map_test.hsp"), + "--sidecar", str(root / "defect_map_test.json"), + "--nside", str(nside), "--nside-coverage", str(NSIDE_COV)]) + out = io.StringIO() + with contextlib.redirect_stdout(out): + merge.main() + return out.getvalue() + + +def _sidecar(root: Path) -> dict: + return json.loads((root / "defect_map_test.json").read_text()) + + +def _spy_reads(merge, monkeypatch) -> list: + """The exposures ``accumulate`` reads into the map, in order.""" + reads = [] + real = merge.accumulate + + def spy(target, paths, nside_coverage, nside): + paths = list(paths) + reads.extend(exp for exp, _ in paths) + return real(target, paths, nside_coverage, nside) + + monkeypatch.setattr(merge, "accumulate", spy) + return reads + + +BOTH = ["2079612p", "2079613p"] + + +@pytest.fixture +def campaign(merge, tmp_path, monkeypatch): + """Two exposures, disjoint pixels, merged once. The starting state.""" + _fragment(merge, tmp_path, "2079612p", [10, 11, 12]) + _fragment(merge, tmp_path, "2079613p", [20, 21]) + out = _main(merge, tmp_path, BOTH, monkeypatch) + assert "rebuilt from 2 fragment(s) (no map on disk)" in out, out + return tmp_path + + +def test_first_merge_is_the_union(merge, campaign): + assert _valid(campaign / "defect_map_test.hsp") == {10, 11, 12, 20, 21} + assert set(_sidecar(campaign)["exposures"]) == set(BOTH) + + +def test_append_reads_only_the_new_fragment(merge, campaign, monkeypatch): + """A grown campaign is an APPEND, not a rebuild — that is the cheap path.""" + _fragment(merge, campaign, "2079614p", [30]) + reads = _spy_reads(merge, monkeypatch) + out = _main(merge, campaign, [*BOTH, "2079614p"], monkeypatch) + assert "1 fragment(s) appended" in out, out + assert reads == ["2079614p"] + assert _valid(campaign / "defect_map_test.hsp") == {10, 11, 12, 20, 21, 30} + assert set(_sidecar(campaign)["exposures"]) == {*BOTH, "2079614p"} + + +def test_removal_forces_a_rebuild_and_drops_the_pixels(merge, campaign, + monkeypatch): + """The case a union cannot do incrementally, and the reason for rebuild.""" + out = _main(merge, campaign, ["2079612p"], monkeypatch) + assert ("rebuilt from 1 fragment(s) (1 exposure(s) left the campaign)" + in out), out + assert _valid(campaign / "defect_map_test.hsp") == {10, 11, 12} + assert set(_sidecar(campaign)["exposures"]) == {"2079612p"} + + +def test_changed_fragment_forces_a_rebuild(merge, campaign, monkeypatch): + """A pixel moved behind an unchanged size and mtime: the DIGEST is the + criterion, and the change drops pixel 21, so only a rebuild — not an + append over the old map — gives the right union. + """ + import os + path = merge.fragment_path(campaign, "2079613p") + st = path.stat() + _fragment(merge, campaign, "2079613p", [20, 22]) + os.utime(path, ns=(st.st_atime_ns, st.st_mtime_ns)) + assert path.stat().st_size == st.st_size, "fixture must preserve the size" + out = _main(merge, campaign, BOTH, monkeypatch) + assert ("rebuilt from 2 fragment(s) (1 fragment(s) changed on disk)" + in out), out + assert _valid(campaign / "defect_map_test.hsp") == {10, 11, 12, 20, 22} + + +def test_moved_defect_changes_the_manifest_and_the_merge_sees_it( + merge, tmp_path, monkeypatch): + """End to end through both scripts' entry points: a flag that moves while + the per-CCD healpix counts stay put changes the fragment's MANIFEST — the + DAG edge the merge waits on — and the merge folds the new fragment in.""" + import os + import numpy as np + from astropy.io import fits + from astropy.wcs import WCS + + spec = importlib.util.spec_from_file_location( + "_defect_map_exp", SCRIPTS / "defect_map_exp.py") + raster = importlib.util.module_from_spec(spec) + spec.loader.exec_module(raster) + + exp = "2079612p" + split = raster.split_dir(tmp_path / "scratch") + split.mkdir(parents=True) + wcs = WCS(naxis=2) + wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] + wcs.wcs.crval = [150.0, 30.0] + wcs.wcs.crpix = [150.0, 150.0] + wcs.wcs.cdelt = [-0.187 / 3600, 0.187 / 3600] + fits.PrimaryHDU(header=wcs.to_header()).writeto( + split / f"image-{exp}-0.fits") + flag = split / f"flag-{exp}-0.fits" + frag = merge.fragment_path(tmp_path, exp) + manifest = merge.manifest_path(tmp_path, exp) + + def rasterize(row_col): + flags = np.zeros((300, 300), dtype=np.int16) + flags[row_col] = 1 + fits.PrimaryHDU(flags).writeto(flag, overwrite=True) + monkeypatch.setattr(sys, "argv", [ + "defect_map_exp.py", "--exp-dir", str(tmp_path / "scratch"), + "--exp", exp, "--dest", str(frag.parent), + "--manifest", str(manifest), "--n-ccds", "1", + "--nside", str(NSIDE), "--nside-coverage", str(NSIDE_COV)]) + raster.main() + + rasterize((50, 50)) + _main(merge, tmp_path, [exp], monkeypatch) + before_manifest, before_pixels, st = (manifest.read_bytes(), _valid(frag), + frag.stat()) + + rasterize((250, 250)) + os.utime(frag, ns=(st.st_atime_ns, st.st_mtime_ns)) + assert _valid(frag) != before_pixels + assert len(_valid(frag)) == len(before_pixels), "fixture must keep counts" + assert frag.stat().st_size == st.st_size, "fixture must keep the size" + assert manifest.read_bytes() != before_manifest, ( + "the fragment moved but its manifest — the merge's DAG edge — did not") + + _main(merge, tmp_path, [exp], monkeypatch) + assert _valid(tmp_path / "defect_map_test.hsp") == _valid(frag) + + +def test_no_op_leaves_the_map_untouched(merge, campaign, monkeypatch): + """UNTOUCHED, not rewritten identically: mtime is a rerun trigger.""" + output = campaign / "defect_map_test.hsp" + sidecar = campaign / "defect_map_test.json" + before = output.stat().st_mtime_ns, sidecar.stat().st_mtime_ns + reads = _spy_reads(merge, monkeypatch) + out = _main(merge, campaign, BOTH, monkeypatch) + assert out.startswith("[merge_defect_map] unchanged"), out + assert "sidecar refreshed" not in out + assert reads == [] + assert (output.stat().st_mtime_ns, sidecar.stat().st_mtime_ns) == before + + +def test_no_op_still_refreshes_a_stale_sidecar(merge, campaign, monkeypatch): + """The campaign moved, the map could not: the RECORD must still say so. + + Tiles whose exposures were all reclaimed by a workflow predating this rule + add nothing to merge and nothing to remove — an empty plan — but they change + what the campaign asked for. A sidecar that kept reporting the old counts + would make a short map look complete on disk. + """ + output = campaign / "defect_map_test.hsp" + before_map = output.stat().st_mtime_ns + out = _main(merge, campaign, BOTH, monkeypatch, missing=["2079999p"]) + assert out.startswith("[merge_defect_map] unchanged"), out + assert "sidecar refreshed" in out + + after = _sidecar(campaign) + assert after["exposures_without_fragment"] == ["2079999p"] + assert after["campaign_exposures"] == 3 + assert after["generation"] == merge.map_generation(output), ( + "the refreshed sidecar must still describe the map beside it") + assert output.stat().st_mtime_ns == before_map, "map must not move" + + +def test_interrupted_publish_is_detected_and_rebuilt(merge, tmp_path, + monkeypatch): + """A kill between replacing the map and replacing the sidecar leaves a + NEW map under an OLD sidecar. Returning to the campaign that old sidecar + describes must rebuild, not read as a no-op over a map that still carries + the other campaign's pixels. + + The interruption is simulated by putting the old sidecar back after a + complete merge: the resulting pair is exactly what the kill leaves. + """ + _fragment(merge, tmp_path, "2079612p", [10, 11]) + _fragment(merge, tmp_path, "2079613p", [20]) + output = tmp_path / "defect_map_test.hsp" + sidecar = tmp_path / "defect_map_test.json" + + _main(merge, tmp_path, ["2079612p"], monkeypatch) + old_sidecar = sidecar.read_bytes() + _main(merge, tmp_path, ["2079612p", "2079613p"], monkeypatch) + assert _valid(output) == {10, 11, 20} + sidecar.write_bytes(old_sidecar) + + out = _main(merge, tmp_path, ["2079612p"], monkeypatch) + assert "rebuilt" in out, out + assert _valid(output) == {10, 11} + assert set(json.loads(sidecar.read_text())["exposures"]) == {"2079612p"} + + +def test_resolution_change_is_not_a_no_op(merge, tmp_path, monkeypatch): + """A new ``--nside`` over unchanged fragments must not pass as a no-op. + + The fragments did not change, so without a resolution check the plan is + empty, the map stays at the old nside and the sidecar is rewritten to + claim the new one. The plan has to see the map's resolution and rebuild, + and the rebuild refuses fragments at another resolution before either + file is touched. + """ + import healsparse as hsp + + _fragment(merge, tmp_path, "2079612p", [10, 11]) + _main(merge, tmp_path, ["2079612p"], monkeypatch) + output = tmp_path / "defect_map_test.hsp" + sidecar = tmp_path / "defect_map_test.json" + before = sidecar.read_bytes(), output.stat().st_mtime_ns + + with pytest.raises(SystemExit) as exc: + _main(merge, tmp_path, ["2079612p"], monkeypatch, nside=2 * NSIDE) + assert "re-rasterize 2079612p" in str(exc.value) + assert (sidecar.read_bytes(), output.stat().st_mtime_ns) == before + assert json.loads(sidecar.read_text())["nside"] == NSIDE + assert hsp.HealSparseCoverage.read(str(output)).nside_sparse == NSIDE + + +def test_nside_mismatch_is_an_error_not_an_upgrade(merge, campaign, + monkeypatch): + """The ladder's resolution is a campaign decision, not a per-fragment one.""" + _fragment(merge, campaign, "2079615p", [5], nside=NSIDE // 2) + with pytest.raises(SystemExit) as exc: + _main(merge, campaign, [*BOTH, "2079615p"], monkeypatch) + assert "re-rasterize 2079615p" in str(exc.value) diff --git a/tests/unit/test_exp_footprint.py b/tests/unit/test_exp_footprint.py new file mode 100644 index 000000000..5d114ee1f --- /dev/null +++ b/tests/unit/test_exp_footprint.py @@ -0,0 +1,275 @@ +"""The exp_footprint record, and the CCD-index contract it rests on. + +``workflow/scripts/exp_footprint.py`` joins two records that name CCDs in two +different ways and never cross-check each other: + + * ``headers-.npy`` names a CCD by its POSITION in the array — the same + ``idx-1`` split_exp used for ``image--.fits``; + * ``exp_persist.json`` names a CCD inside a FILENAME, + ``validation_psf--.fits``, written by psfex_interp. + +The whole design depends on those being the same integer, and no code asserts +it: split_exp writes the image and the array element in one loop, psfex_interp +inherits the numbering through the file handler, and the exposure-count map +would be silently WRONG — right pixels, wrong exposure count — if they ever +diverged by a permutation. This module is where that contract is pinned. + +The fixtures are real astropy WCSs, one per CCD with a DISTINCT centre, so a +permutation of the array shows up as corners on the wrong ``id`` rather than as +a passing test. One CCD straddles RA=0 deliberately: the record is raw sky, and +the seam is the map builder's business (``unwrap_ra``), not this script's. + +Deliberately not a module test: nothing here runs shapepipe, only its two +geometry helpers, so it belongs with the fast structural suite. +""" + +import importlib.util +import json +import sys +from pathlib import Path + +import numpy as np +import pytest +from astropy.io.fits import Header +from astropy.wcs import WCS + +REPO_ROOT = Path(__file__).resolve().parents[2] +SCRIPTS = REPO_ROOT / "workflow" / "scripts" +SCRIPT = SCRIPTS / "exp_footprint.py" + +EXP = "2605805" +# MegaCam-ish: 2048 x 4612 pixels at 0.187"/pixel, i.e. ~0.106 x 0.240 deg. +NX, NY = 2048, 4612 +PIXSCALE = 0.187 / 3600.0 + + +def _load(): + """Import the script by path — ``workflow/scripts`` is not a package. + + Its ``persist_exp`` import is a sibling 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("_exp_footprint", SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + finally: + sys.path.remove(str(SCRIPTS)) + return module + + +exp_footprint = _load() + + +def ccd_header(ra_centre, dec_centre): + """A single-CCD header + WCS, as split_exp stores them. + + The header carries plain ``NAXIS1/2`` and no ``ZIMAGE``, which is what + astropy hands back for a tile-compressed HDU once it is decompressed — the + case ``_image_shape`` takes its second branch for. + """ + w = WCS(naxis=2) + w.wcs.ctype = ["RA---TAN", "DEC--TAN"] + w.wcs.crpix = [NX / 2.0, NY / 2.0] + w.wcs.crval = [ra_centre, dec_centre] + w.wcs.cdelt = [-PIXSCALE, PIXSCALE] + header = w.to_header() + header["NAXIS"] = 2 + header["NAXIS1"] = NX + header["NAXIS2"] = NY + return {"WCS": w, "header": header.tostring()} + + +def headers_array(centres): + """``headers-.npy``'s content: object array, one entry per CCD.""" + arr = np.zeros(len(centres), dtype="O") + for i, (ra, dec) in enumerate(centres): + arr[i] = ccd_header(ra, dec) + return arr + + +# Six CCDs on a row, each 0.3 deg from the last so no two share a footprint, and +# CCD 3 sitting on the RA=0 seam. +CENTRES = [(10.0, 30.0), (10.3, 30.0), (10.6, 30.0), + (0.0, 30.0), (11.2, 30.0), (11.5, 30.0)] + + +@pytest.fixture +def store(tmp_path): + """A scratch exposure store with a headers npy, and a products root.""" + npy_dir = tmp_path / "exp" / EXP / exp_footprint.HEADERS_DIR + npy_dir.mkdir(parents=True) + np.save(npy_dir / f"headers-{EXP}.npy", headers_array(CENTRES)) + return tmp_path + + +def persist_manifest(path, ccds, patterns=("validation_psf-*.fits",), + label=True): + """An ``exp_persist`` manifest naming exactly ``ccds`` as PSF-bearing. + + Shaped as persist_exp.py writes it, including the decoy member: a keep list + of several patterns packs files this script must ignore, and reading a CCD + index out of one of them would be a real bug. + + ``label=False`` drops the ``product`` field, which is how a tar packed + before that field existed reads back — the case the name glob catches. + """ + files = [] + for c in ccds: + entry = {"name": f"validation_psf-{EXP}-{c}.fits", + "pattern": "validation_psf-*.fits", "bytes": 1} + if label: + entry["product"] = "psf_validation" + files.append(entry) + files.append({"name": f"{EXP}-0.psf", "product": "psf_model", + "pattern": "*.psf", "bytes": 1}) + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text(json.dumps( + {"stage": "exp_persist", "unit": EXP, "status": "complete", + "patterns": list(patterns), "n_files": len(files), "files": files}, + indent=2, sort_keys=True)) + return path + + +def run(store, persist, out, monkeypatch): + """Invoke the script's CLI, as the rule's shell does.""" + monkeypatch.setattr(sys, "argv", [ + "exp_footprint.py", + "--exp-dir", str(store / "exp" / EXP), + "--exp", EXP, + "--persist-manifest", str(persist), + "--manifest", str(out)]) + exp_footprint.main() + return json.loads(out.read_text()) + + +def test_ccd_index_alignment(store, monkeypatch): + """npy index ``i`` == ``-`` == ``validation_psf--.fits``. + + The three CCDs with a PSF are non-adjacent on purpose: an off-by-one or a + "renumber the survivors 0..n" bug passes a contiguous fixture and fails + this one. + """ + persist = persist_manifest(store / "prod" / "exp_persist.json", [0, 2, 5]) + record = run(store, persist, store / "prod" / "exp_footprint.json", + monkeypatch) + + assert record["n_ccd_headers"] == len(CENTRES) + assert record["n_valid_psf"] == 3 + assert [c["id"] for c in record["ccds"]] == [ + f"{EXP}-0", f"{EXP}-2", f"{EXP}-5"] + assert record["ccds_no_psf"] == [f"{EXP}-1", f"{EXP}-3", f"{EXP}-4"] + + # The corners on `-i` are the corners of ARRAY ELEMENT i, not of the + # i-th surviving CCD: each fixture CCD has its own centre, so this catches + # any permutation. Checked against the same two helpers the script imports, + # which is the point — the contract under test is the INDEXING, not the + # projection arithmetic (tests/module/test_coverage.py owns that). + arr = headers_array(CENTRES) + for entry in record["ccds"]: + i = int(entry["id"].rsplit("-", 1)[1]) + shape = exp_footprint._image_shape(Header.fromstring(arr[i]["header"])) + ra, dec = exp_footprint._ccd_corners(arr[i]["WCS"], shape) + assert entry["ra"] == pytest.approx(list(ra)) + assert entry["dec"] == pytest.approx(list(dec)) + # ... and that centre really is CCD i's, not its neighbour's. + assert np.mean(entry["dec"]) == pytest.approx(CENTRES[i][1], abs=1e-3) + + +def test_record_is_raw_sky_across_the_ra_seam(store, monkeypatch): + """CCD 3 straddles RA=0 and the record says so, unwrapped by nobody. + + The seam is handled once, in ``coverage_map_builder.unwrap_ra``, at the + moment a polygon is stamped. Unwrapping here as well would put negative RA + into a durable record that other consumers read, and double-unwrapping is + not idempotent. + """ + persist = persist_manifest(store / "prod" / "exp_persist.json", [3]) + record = run(store, persist, store / "prod" / "exp_footprint.json", + monkeypatch) + + ra = record["ccds"][0]["ra"] + assert record["ccds"][0]["id"] == f"{EXP}-3" + assert min(ra) >= 0.0 and max(ra) < 360.0, "raw sky, not a shifted branch" + assert max(ra) - min(ra) > 180.0, "the fixture must actually cross RA=0" + # The four corners land on both sides of the seam, ~0.05 deg out. + assert sorted(round(r) for r in ra) == [0, 0, 360, 360] + + +def test_byte_stable_rerun_keeps_the_mtime(store, monkeypatch): + """A rerun over an unchanged store must not move the manifest's mtime. + + mtime is a rerun trigger, and this manifest is an input of clean_exposure: + an unconditional rewrite would make every downstream reclamation look out + of date once per invocation. Same contract as persist_exp.py's. + """ + persist = persist_manifest(store / "prod" / "exp_persist.json", [0, 2, 5]) + out = store / "prod" / "exp_footprint.json" + + first = run(store, persist, out, monkeypatch) + raw = out.read_bytes() + before = out.stat().st_mtime_ns + + second = run(store, persist, out, monkeypatch) + assert second == first + assert out.read_bytes() == raw + assert out.stat().st_mtime_ns == before + # No stray tmp left behind on the persistent root. + assert not list(out.parent.glob("*.tmp")) + + +def test_unlabelled_members_are_found_by_name(store, monkeypatch): + """A tar packed before persist_exp labelled its members still reads. + + Membership is the product label OR the file name, merge_star_cat's test + exactly. The label is what persist_exp writes today; the name glob is what + an older manifest, or a keep list written as a raw glob, leaves behind. + Neither may make an exposure silently contribute no sky. + """ + persist = persist_manifest(store / "prod" / "exp_persist.json", [0, 2, 5], + label=False) + record = run(store, persist, store / "prod" / "exp_footprint.json", + monkeypatch) + assert [c["id"] for c in record["ccds"]] == [f"{EXP}-{i}" + for i in (0, 2, 5)] + + +def test_the_keep_list_cannot_turn_the_psf_set_off(store, monkeypatch): + """`persist_exp:` is optional retention and no precondition of this rule. + + exp_persist packs every CCD's psf_validation whatever the keep list says + (its ALWAYS), so a manifest whose `patterns` name only other products still + carries the valid-PSF set, and this rule must read it rather than refuse. + """ + persist = persist_manifest(store / "prod" / "exp_persist.json", [0, 2, 5], + patterns=("*.psf",)) + record = run(store, persist, store / "prod" / "exp_footprint.json", + monkeypatch) + assert [c["id"] for c in record["ccds"]] == [f"{EXP}-{i}" + for i in (0, 2, 5)] + + +def test_psf_for_a_ccd_the_split_never_wrote_is_fatal(store, monkeypatch): + """The one disagreement the index alignment cannot absorb, made loud.""" + persist = persist_manifest(store / "prod" / "exp_persist.json", [0, 99]) + with pytest.raises(SystemExit) as exc: + run(store, persist, store / "prod" / "exp_footprint.json", monkeypatch) + assert "99" in str(exc.value) + + +def test_missing_headers_array_is_fatal(store, monkeypatch): + """A purged or unbuilt split store fails here, never against VOS. + + This is the other half of the rule declaring only its DURABLE input + (exposure.smk): a persist manifest outliving its scratch store is a real + state, and it must cost one error line rather than an exposure rebuild. + """ + npy = (store / "exp" / EXP / exp_footprint.HEADERS_DIR + / f"headers-{EXP}.npy") + npy.unlink() + persist = persist_manifest(store / "prod" / "exp_persist.json", [0]) + with pytest.raises(SystemExit) as exc: + run(store, persist, store / "prod" / "exp_footprint.json", monkeypatch) + assert "no WCS array" in str(exc.value) diff --git a/tests/unit/test_footprint_edges.py b/tests/unit/test_footprint_edges.py new file mode 100644 index 000000000..80ced7128 --- /dev/null +++ b/tests/unit/test_footprint_edges.py @@ -0,0 +1,201 @@ +"""Which exposure-footprint records the workflow depends on, and how. + +The footprint record (``exp_footprint.json``) is exp_footprint's declared +output, hanging off exp_persist's manifest and through it off the exposure's +scratch chain. Two rules name it — ``rule all`` and ``clean_exposure`` — through +``footprint_edge()``, which must drop it once the exposure's store is reclaimed: +a named record keeps the producer chain in the DAG, and any upstream params +change then reschedules the exposure's download, split and PSF fit. + +nexp_map reads every record on the products root, edge or not, so what +reruns it is a fingerprint of that whole set on its params +(``nexp_map_exposures()``): it must move when a record arrives off the DAG, and +must NOT move when a record this invocation writes appears or when its store is +later reclaimed — either would rerun a many-hour job over the same records. + +The Snakefile helpers are lifted out by name and evaluated against a temporary +run root and products root, so these tests exercise what the helpers RETURN for +real files on disk. +""" + +import glob +import hashlib +import re +from pathlib import Path +from types import SimpleNamespace + +import pytest + +REPO_ROOT = Path(__file__).resolve().parents[2] +SNAKEFILE = REPO_ROOT / "workflow" / "Snakefile" + +EXP = "2243881" +HELPERS = ("exp_dir", "exp_manifest", "prod_exp_dir", "prod_exp_manifest", + "tombstone", "exp_store_reclaimed", "footprint_edge", + "unit_fingerprint", "nexp_map_exposures", "nexp_map", + "nexp_map_manifest", "nexp_map_targets") +# The ready tiles' exposures, as psf_exposures() would return them. +IN_SCOPE = ["2243881", "2243882"] + + +def _snakefile_def(name): + """The text of one top-level ``def`` in the Snakefile.""" + m = re.search(rf"^def {name}\(.*?(?=^\S)", SNAKEFILE.read_text(), + re.M | re.S) + assert m, f"Snakefile no longer defines {name}()" + return m.group(0) + + +@pytest.fixture +def roots(tmp_path): + """The lifted helpers, bound to a temporary scratch and products root.""" + warnings = [] + ns = {"Path": Path, "glob": glob, "hashlib": hashlib, + "RUN_DIR": tmp_path / "run", "PRODUCTS_DIR": tmp_path / "products", + "psf_exposures": lambda: list(IN_SCOPE), + "MAPS_NEXP": True, "CAMPAIGN": "campaign-sentinel", + "workflow": SimpleNamespace(is_main_process=True), + "logger": SimpleNamespace(warning=warnings.append), + "warnings": warnings} + for name in HELPERS: + exec(_snakefile_def(name), ns) + return ns + + +def touch(path): + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text("{}") + + +def test_live_store_depends_on_its_footprint(roots): + """A live exposure is asked for its record: the thing to build, and what + orders reclamation after the header read.""" + touch(roots["exp_manifest"](EXP, "exp_psf")) + assert roots["footprint_edge"](EXP) == [ + roots["prod_exp_manifest"](EXP, "exp_footprint")] + + +def test_unbuilt_store_depends_on_its_footprint(roots): + """A fresh campaign's exposure has no files yet and must still be asked.""" + assert roots["footprint_edge"](EXP) == [ + roots["prod_exp_manifest"](EXP, "exp_footprint")] + + +def test_cleaned_store_names_no_footprint(roots): + """Tombstoned: the record is on the products root, and naming it would + keep the exposure's producer chain in the DAG.""" + touch(roots["prod_exp_manifest"](EXP, "exp_persist")) + touch(roots["prod_exp_manifest"](EXP, "exp_footprint")) + touch(roots["tombstone"](EXP)) + assert roots["footprint_edge"](EXP) == [] + + +def test_purged_store_names_no_footprint(roots): + """Purged without a tombstone: persisted, and exp_psf's scratch manifest + gone. The tombstone alone cannot see this case.""" + touch(roots["prod_exp_manifest"](EXP, "exp_persist")) + touch(roots["prod_exp_manifest"](EXP, "exp_footprint")) + assert roots["footprint_edge"](EXP) == [] + + +def test_both_footprint_edges_use_footprint_edge(): + """``rule all`` and ``clean_exposure`` reach the record only through the + helper, so the reclaimed-store cut cannot hold on one edge and not the + other.""" + targets = _snakefile_def("footprint_targets") + assert "footprint_edge(" in targets + rules = (REPO_ROOT / "workflow" / "rules" / "exposure.smk").read_text() + clean = re.search(r"^rule clean_exposure:.*?(?=^rule |\Z)", rules, + re.M | re.S).group(0) + assert "footprint_edge(wc.exp)" in clean + code = [line for src in (targets, clean) for line in src.splitlines() + if not line.lstrip().startswith("#")] + assert not any('"exp_footprint"' in line for line in code), code + + +def fingerprint(roots): + return roots["unit_fingerprint"](roots["nexp_map_exposures"]()) + + +def reclaim(roots, exp): + """exp's store after clean_exposure: record kept, tombstone written.""" + touch(roots["prod_exp_manifest"](exp, "exp_persist")) + touch(roots["tombstone"](exp)) + + +def test_footprint_under_reclaimed_exposure_moves_the_fingerprint(roots): + """A record that reaches the root with no edge — written while the map + was off, then its store reclaimed — must still rerun the map.""" + for exp in IN_SCOPE: + reclaim(roots, exp) + touch(roots["prod_exp_manifest"]("2243881", "exp_footprint")) + before = fingerprint(roots) + touch(roots["prod_exp_manifest"]("2243882", "exp_footprint")) + assert fingerprint(roots) != before + + +def test_out_of_scope_footprint_moves_the_fingerprint(roots): + """A record of an exposure this tile list never names is in the glob, so + it is in the fingerprint.""" + before = fingerprint(roots) + touch(roots["prod_exp_manifest"]("9900001", "exp_footprint")) + assert fingerprint(roots) != before + + +def test_writing_a_declared_footprint_keeps_the_fingerprint(roots): + """The params are recorded at parse time, before this invocation writes + its records; the next parse must not see a change.""" + before = fingerprint(roots) + for exp in IN_SCOPE: + touch(roots["prod_exp_manifest"](exp, "exp_footprint")) + assert fingerprint(roots) == before + + +def test_reclaiming_keeps_the_fingerprint(roots): + """An exposure moving from edge to glob-only is the same record.""" + for exp in IN_SCOPE: + touch(roots["prod_exp_manifest"](exp, "exp_footprint")) + before = fingerprint(roots) + for exp in IN_SCOPE: + reclaim(roots, exp) + assert fingerprint(roots) == before + + +def test_nexp_map_params_carry_the_fingerprint(): + rule = (REPO_ROOT / "workflow" / "rules" / "exposure.smk").read_text() + params = re.search(r"^rule nexp_map:.*?^ params:\n(.*?)^ \w", + rule, re.M | re.S).group(1) + assert "unit_fingerprint(nexp_map_exposures())" in params, params + + +def test_no_record_anywhere_requests_no_map(roots): + """Every in-scope exposure reclaimed before exp_footprint existed: no record + on the root and none declared. Requesting the map would be a job that fails + on every invocation; it is not requested, and the loss is said.""" + for exp in IN_SCOPE: + reclaim(roots, exp) + assert roots["nexp_map_targets"]() == [] + assert len(roots["warnings"]) == 1 + assert "2 of this campaign's 2 exposure(s)" in roots["warnings"][0] + assert "no map is built" in roots["warnings"][0] + + +def test_partial_records_build_an_undercounting_map_and_say_so(roots): + """One exposure reclaimed with its record, one without: the map is built + from what there is, and the parse warns that it undercounts.""" + for exp in IN_SCOPE: + reclaim(roots, exp) + touch(roots["prod_exp_manifest"]("2243881", "exp_footprint")) + assert roots["nexp_map_targets"]() == [ + roots["nexp_map"](), roots["nexp_map_manifest"]()] + assert len(roots["warnings"]) == 1 + assert "1 of this campaign's 2 exposure(s)" in roots["warnings"][0] + assert "undercounts" in roots["warnings"][0] + + +def test_live_exposures_request_the_map_silently(roots): + """A fresh campaign: every exposure is declared, nothing is lost.""" + assert roots["nexp_map_targets"]() == [ + roots["nexp_map"](), roots["nexp_map_manifest"]()] + assert roots["warnings"] == [] diff --git a/tests/unit/test_nexp_map_gates.py b/tests/unit/test_nexp_map_gates.py new file mode 100644 index 000000000..654fb695e --- /dev/null +++ b/tests/unit/test_nexp_map_gates.py @@ -0,0 +1,139 @@ +"""Fitted-PSF gates on footprints, reclamation and the exposure-count map. + +Failure modes: fake PSFs request nonexistent persistence products; a fitted +PSF loses a durable edge before cleanup; a fake-PSF (sim) campaign requests +an exposure-count map it cannot build, or a fitted one silently loses it. Lift the +Snakefile's functions and the rule's input lambdas into sentinel roots, as +test_campaign_lineage does, without parsing a second workflow DAG. +""" + +import re +import textwrap +from pathlib import Path +from types import SimpleNamespace + +import pytest + +WORKFLOW = Path(__file__).resolve().parents[2] / "workflow" +SNAKEFILE = WORKFLOW / "Snakefile" +PRODUCTS = "/sentinel/products" +EXPOSURES = ["2605805", "2700001"] + + +class WorkflowError(Exception): + """Stand-in for snakemake.exceptions.WorkflowError, which the image lacks.""" + + +def _namespace(psf_model, tmp_path, *, main=True): + """Bind the real PSF gate and target helpers to a small indexed campaign.""" + tile_exp = { + "210.282": [EXPOSURES[1], EXPOSURES[0]], + "211.282": [EXPOSURES[0]], + "999.999": ["9900001"], + } + ns = { + "PSF_MODEL": psf_model, + "PRODUCTS_DIR": Path(PRODUCTS), + "RUN_DIR": tmp_path / "scratch", + "CAMPAIGN": "campaign-sentinel", + "TILES_READY": ["210.282", "211.282"], + "READY_SET": {"210.282", "211.282"}, + "tile_exposures": tile_exp.__getitem__, + "clean_consumers": lambda exp: ["210.282", "999.999"], + "workflow": SimpleNamespace(is_main_process=main), + "WorkflowError": WorkflowError, + "Path": Path, + "script_hash": lambda name: "hash-sentinel", + "glob": __import__("glob"), + "logger": SimpleNamespace(warning=lambda message: None), + } + text = SNAKEFILE.read_text() + assignment = re.search(r"^PERSISTS_PSF\s*=.*$", text, re.M) + assert assignment, "Snakefile must bind PERSISTS_PSF" + exec(assignment.group(0), ns) + # Fake PSFs are only legal with image sims, which rasterize no defects. + ns["INPUT_TYPE"] = "image_sims" if psf_model == "fake" else "data" + defects = re.search(r"^MAPS_DEFECTS\s*=.*$", text, re.M) + assert defects, "Snakefile must bind MAPS_DEFECTS" + exec(defects.group(0), ns) + for name in ("prod_exp_dir", "prod_exp_manifest", "exp_dir", "tombstone", + "exp_manifest", "exp_store_reclaimed", "footprint_edge", + "tile_dir", "tile_manifest", "psf_exposures", + "footprint_targets", "flag", "prod_exp_tar", + "prod_exp_fragment", "nexp_map", "nexp_map_manifest", + "nexp_map_exposures", "nexp_map_targets"): + definition = re.search( + rf"^def {name}\(.*?(?=^\S|\Z)", text, re.M | re.S) + assert definition, f"Snakefile must define {name}()" + exec(definition.group(0), ns) + return ns + + +def _clean_inputs(ns): + """Evaluate clean_exposure's actual input functions, not a copy of them.""" + text = (WORKFLOW / "rules" / "exposure.smk").read_text() + rule = re.search(r"^rule clean_exposure:\n(.*?)(?=^rule |\Z)", + text, re.M | re.S) + assert rule, "exposure.smk must define clean_exposure" + inputs = re.search(r"^ input:\n(.*?)(?=^ \w+:)", + rule.group(1), re.M | re.S) + assert inputs, "clean_exposure must declare its ordering edges" + return eval("[\n" + textwrap.dedent(inputs.group(1)) + "\n]", ns) + + +def _nexp_gate(ns, nexp): + """Execute the Snakefile's exposure-count-map gate against a `nexp:` block.""" + text = SNAKEFILE.read_text() + gate = re.search(r"^_NEXP\s*=.*\n^MAPS_NEXP\s*=.*\n", text, re.M) + assert gate, "Snakefile must gate exposure_maps.nexp on a fitted PSF" + ns["_MAPS"] = {} if nexp is None else {"nexp": nexp} + exec(gate.group(0), ns) + + +@pytest.mark.parametrize("psf_model", ["fake", "psfex"]) +@pytest.mark.parametrize("main", [True, False]) +def test_footprint_targets_require_fitted_psf(psf_model, main, tmp_path): + """Only a head-process fitted-PSF run requests the ready exposures.""" + ns = _namespace(psf_model, tmp_path, main=main) + expected = [] + if psf_model == "psfex" and main: + expected = [ + f"{PRODUCTS}/exp/{exp[:2]}/{exp}/manifests/exp_footprint.json" + for exp in EXPOSURES + ] + assert ns["footprint_targets"]() == expected + + +@pytest.mark.parametrize("psf_model", ["fake", "psfex"]) +def test_clean_exposure_waits_on_persist_iff_psf(psf_model, tmp_path): + """Reclamation waits for both durable reads iff a PSF is fitted.""" + ns = _namespace(psf_model, tmp_path) + wc = SimpleNamespace(exp=EXPOSURES[0]) + edges = [path for function in _clean_inputs(ns) for path in function(wc)] + expected = [str(tmp_path / "scratch" / "tiles" / "21" / "210.282" + / "manifests" / "tile_vignets.json")] + if psf_model == "psfex": + stages = ["exp_persist"] + if ns["MAPS_DEFECTS"]: + stages.append("exp_defect_map") + stages.append("exp_footprint") + expected += [ + f"{PRODUCTS}/exp/26/2605805/manifests/{stage}.json" + for stage in stages + ] + assert edges == expected + + +@pytest.mark.parametrize("psf_model", ["fake", "psfex"]) +@pytest.mark.parametrize("nexp", [None, {}, {"enabled": True}, + {"enabled": False}, {"enabled": "true"}, + {"enabled": "false"}]) +def test_nexp_map_follows_the_psf_model(psf_model, nexp, tmp_path): + """On by default for a fitted PSF, off only by explicit opt-out; a fake-PSF + (sim) campaign skips it silently, whatever the block says.""" + ns = _namespace(psf_model, tmp_path) + _nexp_gate(ns, nexp) + opted_out = (nexp or {}).get("enabled") in (False, "false") + built = psf_model == "psfex" and not opted_out + expected = [ns["nexp_map"](), ns["nexp_map_manifest"]()] if built else [] + assert ns["nexp_map_targets"]() == expected diff --git a/tests/unit/test_run_config.py b/tests/unit/test_run_config.py index 3290f9fc3..ffa6f16c6 100644 --- a/tests/unit/test_run_config.py +++ b/tests/unit/test_run_config.py @@ -9,6 +9,7 @@ through both, and an optional path left holding a ``$`` is reported. """ +import re import importlib.util from pathlib import Path @@ -96,3 +97,36 @@ def test_optional_path_with_an_unknown_variable_is_reported(): inputs={"masks": "/m/$nope"}, container="/c/$nope.sif")) assert set(run_config.unresolved(cfg)) == { "outputs.products_dir", "inputs.masks", "container"} + + +def test_shipped_config_carries_no_retired_key(): + assert run_config.retired(yaml.safe_load(CONFIG_YAML.read_text())) == [] + + +def test_retired_keys_are_found_wherever_they_sit(): + """`coverage:` / `defect_map:` at the top level, on a machine, or on a + machine's input_type block are each reported with their replacement.""" + config = { + "coverage": {"enabled": False}, + "exposure_maps": {"nexp": {"enabled": True}}, + "machines": { + "nibi": {"defect_map": {"oversample": 3}, + "data": {"coverage": {"nside": 131072}}}, + "candide": {"image_sims": {"retrieve": "symlink"}}, + }, + } + found = dict(run_config.retired(config)) + assert set(found) == {"coverage", "machines.nibi.defect_map", + "machines.nibi.data.coverage"} + assert "exposure_maps.nexp" in found["coverage"] + assert "exposure_maps.defect" in found["machines.nibi.defect_map"] + assert "exposure_maps.nside" in found["machines.nibi.data.coverage"] + + +def test_snakefile_refuses_retired_keys_at_parse_time(): + """The Snakefile raises on retired() before anything else reads config.""" + text = (REPO_ROOT / "workflow" / "Snakefile").read_text() + check = text.index("run_config.retired(config)") + defaults = re.search(r"run_config\.apply_(machine_)?defaults\(config\)", text) + assert defaults is not None and check < defaults.start() + assert "raise WorkflowError" in text[check:check + 400] diff --git a/tests/workflow/harness.py b/tests/workflow/harness.py index e3e487b28..e90c534a4 100644 --- a/tests/workflow/harness.py +++ b/tests/workflow/harness.py @@ -151,6 +151,11 @@ def persist_manifest(self, exp): return (self.products_dir / "exp" / exp[:2] / exp / "manifests" / "exp_persist.json") + def exp_manifest(self, exp, stage): + """Return one exposure's persistent manifest for a named stage.""" + return (self.products_dir / "exp" / exp[:2] / exp + / "manifests" / f"{stage}.json") + @dataclass class ResolvedDAG: diff --git a/tests/workflow/params_pin.json b/tests/workflow/params_pin.json index 3adfce597..4fa3df537 100644 --- a/tests/workflow/params_pin.json +++ b/tests/workflow/params_pin.json @@ -4,11 +4,15 @@ "all": "572122d8d1901e12ff591b30adf405f8920be8183129befececbb819f3392ed8", "clean_exposure": "22cb76b13a5205d20a02a9bd3b8c8bea5ea2801555e24dd2b84b11145f7e79d9", "clean_tile": "a5c07b0461526ed407df36a291deb866d4181524b4c8e0fbc3dd047fd9d28479", + "defect_map_merge": "397ea7c2cfdaaac5c5246ee2eaac71152636c488dee4d5e53ffe4477e26af7f6", + "exp_defect_map": "4a411b00a444a96b9bd133c7c8db5d36e3b0e187227d4acc08f790ba28bb5b9f", + "exp_footprint": "1b6a96dfaf637bd2c190d6ea4670043871716758d96967a5c37a220843af1ff4", "exp_get_images": "a9e9eca0a99a348e43e0cd44d0c0aae0f890773d03e35007fb77f886b7d21f5f", "exp_persist": "302e2837542bc1102430c27c81c600b7cda32e8bddcb5fd60d33950987609fff", "exp_psf": "87b9fd799a54130eb3dd9c78bb5124838854c23f365df7a991a53f4c6ac4c3ff", "exp_split": "8fc97dd48e76cd6c7fdf9ee388ba1c1a3dd63428cbe1961641fbebe76b05786e", "final_cat_merge": "e7f46859c4503a2220713d7bb2507555515d0a9632d780b20f14c59e32210023", + "nexp_map": "34340fa3ebca7cf3efc81fcf3da6c10627894c61170d717e3572bd60fa778bb6", "prepare_all_tiles": "b8f872a22adf014e25a7fa5198f49b71a6fe9e56042ed82b682bc8763970a844", "star_cat_merge": "6277450958474af5270982fa35360f2f237a29f7533c526ee9265dfd5acc07a0", "tile_detect": "602bee6fbbb7060e1169ff9bb8521a01194b29ab78e03b0f1df3a580832cd079", @@ -23,7 +27,7 @@ "tile_vignets": "b9853b2a1830840a21d8ec1f0877e780868a07da57a29f531df47f26fe8692ac" }, "schema": 1, - "sha256": "4235cff2328f83cc78ea73c7fd350dbdf4cd4b73c61b6b64b8db1e73e47c1aa3", + "sha256": "5c8f0b40ef0442a79a3801c66d41231f4ab92101b98d4fa1b5971c1d16c9050e", "unit_pre": { "exp_get_images": "dac6685a207dae3ea81d296ce53e9ba75d4636cf1f81e6399f2dda452680c992", "exp_psf": "87fc8ea1153a709ab7ba56542cfacc041f5f8b7f6cf2f38a9bc9d8004d1cecb4", diff --git a/tests/workflow/test_dag.py b/tests/workflow/test_dag.py index e582138bc..f7777ca8f 100644 --- a/tests/workflow/test_dag.py +++ b/tests/workflow/test_dag.py @@ -15,7 +15,8 @@ "tile_ngmix", "tile_merge_cats", "tile_make_cat", "clean_tile", "final_cat_merge", } -PSF_RULES = {"exp_persist", "star_cat_merge"} +PSF_RULES = {"exp_persist", "star_cat_merge", "exp_footprint", "nexp_map"} +DEFECT_RULES = {"exp_defect_map", "defect_map_merge"} # MAPS_DEFECTS: data only def test_rule_set_matches_input_mode(campaign, dag): @@ -23,6 +24,8 @@ def test_rule_set_matches_input_mode(campaign, dag): expected = BASE_RULES.copy() if campaign.psf_model != "fake": expected |= PSF_RULES + if campaign.input_type == "data": + expected |= DEFECT_RULES assert dag.rule_names == expected assert "merge_final_cats" not in dag.declared_rule_names @@ -39,6 +42,10 @@ def test_clean_exposure_waits_on_persist_iff_psf(campaign, dag): ] if campaign.psf_model != "fake": expected.append(campaign.persist_manifest(exp)) + if campaign.input_type == "data": + expected.append(campaign.exp_manifest(exp, "exp_defect_map")) + if campaign.psf_model != "fake": + expected.append(campaign.exp_manifest(exp, "exp_footprint")) assert Counter(map(str, job.input)) == Counter(map(str, expected)), ( "clean-exposure-waits-on-persist-iff-psf", exp, list(job.input) ) @@ -109,6 +116,17 @@ def test_missing_run_fails_during_parse(campaign, resolve_dag): pytest.fail("a campaign without run: must fail at parse time") +@pytest.mark.parametrize("key", ["coverage", "defect_map"]) +def test_retired_map_keys_fail_during_parse(campaign, resolve_dag, key): + """A run config still carrying a retired block is refused, naming its + replacement, instead of parsing with the block silently ignored.""" + campaign.config[key] = {"enabled": False} + campaign.write_config() + with pytest.raises(WorkflowError, match=rf"{key} -> exposure_maps\."): + with resolve_dag(campaign): + pytest.fail(f"a run config with {key}: must fail at parse time") + + def test_mccd_is_refused_during_parse(tmp_path, resolve_dag): """MCCD products are unreadable to persistence and the star merge.""" campaign = Campaign(tmp_path / "campaign", "data", "mccd") diff --git a/universes/committed.yaml b/universes/committed.yaml index 977f02878..de635a9be 100644 --- a/universes/committed.yaml +++ b/universes/committed.yaml @@ -11,6 +11,8 @@ analyses: psf_star_mask_veto: instrument_flags_only sky_mask_application: catalogue_columns mask_default_cut: r_mask_bits + defect_map_from_flags: rasterize_conservative_os3 + nexp_map_valid_psf_ccds: valid_psf_ccds detection: decisions: detection_threshold_policy: megapipe_tiles diff --git a/workflow/CONTRACTS b/workflow/CONTRACTS index 9538e78cf..0877ad091 100644 --- a/workflow/CONTRACTS +++ b/workflow/CONTRACTS @@ -9,8 +9,10 @@ prose names it. A campaign's name is the run config's `run:` and nothing else. The Snakefile binds it once, `CAMPAIGN = config["run"]`, and every product that carries a campaign name takes it from there: `final_cat_.hdf5` and its -`patches/` group, `full_starcat_.hdf5`, and the `$run` in the -machine defaults' `products_dir`/`index_db`. No rule or script reads a +`patches/` group, `full_starcat_.hdf5`, +`defect_map/defect_map_.hsp`, `nexp_map/nexp_map_.hsp`, and the +`$run` in the machine defaults' `products_dir`/`index_db`. No rule or script +reads a `campaign` key or derives a name from a directory (`PRODUCTS_DIR.name`); two sources that can disagree would file one campaign's merge under another's name. Every campaign product path is rooted in `PRODUCTS_DIR`. Enforced by @@ -53,3 +55,15 @@ a rule that is not in the DAG; naming the manifest for a reclaimed store lets a tests/workflow/test_dag.py::test_clean_exposure_waits_on_persist_iff_psf, which also checks that the consumer edges are exactly the in-scope vignets manifests. +`clean_exposure` also waits on the exposure's `exp_footprint` manifest while +the store is live (`footprint_edge`, Snakefile): the footprint is derived from +headers in the store this job deletes. Enforced by +tests/unit/test_nexp_map_gates.py and tests/unit/test_footprint_edges.py. + +@sc [label:custody] nexp-map-requires-fitted-psf +`footprint_targets()` requests the ready tiles' exposures only for a fitted +PSF model in the head process, and `nexp_map_targets()` requests the map only +under the same condition (`MAPS_NEXP`). Under `psf_model: fake` no exposure has +a valid-PSF CCD set to record, so a sim campaign requests neither — whatever +`exposure_maps.nexp.enabled` says — and parses without error. Enforced by +tests/unit/test_nexp_map_gates.py and tests/workflow/test_dag.py. diff --git a/workflow/README.md b/workflow/README.md index a5ced4172..d465b0b69 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -61,8 +61,10 @@ each overlay file to its `cfis/` original confined to input naming. `psf_model: fake` is the simulations' true PSF: the exposure stage runs only SExtractor (for the background maps the vignets read), and `tile_vignets` runs `fake_interp_runner`, which writes the `galaxy_psf` product from `psf_dict`. -With no PSF model there is nothing to persist per exposure, so `exp_persist` and -`star_cat_merge` do not run and `clean_exposure` does not wait on them. +With no PSF model there is nothing to persist per exposure and no valid-PSF CCD +set to record, so `exp_persist`, `exp_footprint` and `star_cat_merge` do not run +and `clean_exposure` does not wait on them; there is no exposure-count map to +build. Simulations that contain stars can run `psfex` or `mccd` exactly as the data do. One campaign per shear branch, each with its own run config: @@ -230,7 +232,7 @@ workflow/ bin/sp committed launcher (module load + /project venv + launch code snapshot + run/report/container/cancel) rules/ prepare.smk tile get_images/uncompress/find_exposures - exposure.smk per-exposure: get_images, split, psf, persist (no temp()); campaign star_cat_merge + exposure.smk per-exposure: get_images, split, psf, persist, footprint, defect_map (no temp()); campaign star_cat_merge, defect_map_merge, nexp_map tile.smk per-tile: exp forest, merge_headers, detect, vignets, ngmix, merge, make_cat; campaign final_cat_merge scripts/ build_index.py prepare-phase run_index.sqlite builder (plain script) @@ -243,6 +245,10 @@ workflow/ merge_star_cat.py ALL exposures' validation_psf, out of the tars -> full_starcat_.hdf5 merge_final_cat.py ALL tiles' final_cat -> final_cat_.hdf5 (the final_cat_merge rule) clean_exposure.py ONE exposure's store + manifests + logs -> tombstone (the clean_exposure rule) + defect_map_exp.py ONE exposure's per-CCD instrument flags -> a boolean healsparse fragment (the exp_defect_map rule) + merge_defect_map.py ALL exposures' fragments -> defect_map/defect_map_.hsp (the defect_map_merge rule) + exp_footprint.py ONE exposure's per-CCD sky corners, for the CCDs with a PSF (the exp_footprint rule) + nexp_map.py EVERY exposure footprint on products_dir -> nexp_map/nexp_map_.hsp (the nexp_map rule) profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; keep-going ``` @@ -293,6 +299,86 @@ profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; kee catalogue server, staged, or rasterized, which is why the old `star_catalogue` / `exp_star_cat` / `exp_mask` rules and their cache root are gone. +- **Two HealSparse maps are built from the exposures.** Both are per-exposure + records on the persistent root combined by one campaign job, both are + products nothing in the workflow reads back, and both sit at one shared + resolution, `config.yaml`'s `exposure_maps:` `nside`/`nside_coverage` — + the mask ladder's 131072 over 128, so they align pixel-wise with the UNIONS + bit masks and with each other. + + *The defect map: the one pixel-domain mask leaves the pixel domain.* The instrument + flag image is the exception to everything above: bad columns, saturated + pixels and bleed trails, split per CCD by `exp_split` and read by SExtractor + as `IMAFLAGS_ISO`, and never sky-fixed. The survey footprint is built from the + CCD corner WCS in the headers, so it cannot subtract them — the footprint + would silently include defective pixels, and the lost area, though only + percent-level, carries exactly the thin small-scale geometry an accurate + window function needs ([#878](https://github.com/CosmoStat/shapepipe/issues/878)). + So `exp_defect_map` rasterizes each exposure's flags into a boolean healsparse + fragment on the persistent root, and `defect_map_merge` unions the campaign's + fragments into `/defect_map/defect_map_.hsp`. Same form as every + other map here — nside 131072 over coverage 128, `True` = masked — so it drops + into the ladder unchanged. **It is a product, not an input:** nothing in the + workflow reads it back, and `config_tile_Mc.ini` deliberately does not name it + in `MASK_EXT_PATHS`, because that ladder names maps that exist before the run + and this one exists only after it; that file's header carries the recipe for + adding it once a campaign has produced one. The rule hangs off `exp_split`, + not off `exp_psf`, so re-rasterizing the campaign at a different fidelity + never touches the PSF chain, and `clean_exposure` takes its manifest as an + input for a LIVE exposure, so reclamation cannot overtake the copy — and only + for a live one: an exposure whose store went to the /scratch purge (no + tombstone, nothing left to rasterize) is asked for its existing fragment if it + has one and for nothing if it does not, the same split `defect_map_merge`'s + input makes, because requiring a manifest behind a vanished split dir would + rebuild the whole exposure chain from VOS to reclaim it. Its resolution and its + oversampling ride on `params`; `config.yaml`'s `exposure_maps:` block carries + both, the measured convergence table behind the default, and the measurement + on one real exposure (34 s, 0.62 GB, a 2.0 MB fragment; the rasterization is + batched, so a fully flagged chip — the worst case, and one a real exposure + carries whenever a chip is dead — is 0.74 GB rather than several). The union + RECONCILES like `final_cat_merge` — a new exposure is OR-ed in on the spot, an + exposure that left the campaign or a fragment that changed forces a rebuild + (a union cannot be un-OR-ed), and a no-op leaves the file untouched — against + a sidecar `defect_map/defect_map_.json` that records which exposures are + already in it. Memory is flat in the exposure count: fragments are read one at + a time and reduced to their pixel ids, so the job holds one accumulator (the + campaign's footprint, ~3 GB at DR6 scale) and one 2 MB fragment. What the map + and the sidecar say is a function of the input set; the map's BYTES are not, + because reaching a state by append rather than by rebuild round-trips it + through healsparse's writer (`merge_defect_map.py` measures the difference and + says what would have to change if anything ever consumed the map). + `tests/unit/test_defect_map_reconcile.py` pins the four reconcile branches and + the sidecar refresh. + + *The exposure-count map.* `exp_footprint` writes one JSON per exposure to + `/exp///manifests/exp_footprint.json`, giving the + four sky corners of every CCD that got a PSF model. It reads the valid-PSF CCD + set off `exp_persist.json`'s tar members — exact, because `psfex_interp` + returns *without* writing `validation_psf-*.fits` on NOT_ENOUGH_STARS, + BAD_CHI2 or FILE_NOT_FOUND — and the WCS off `headers-.npy`, written by + `exp_split`. It needs nothing of `persist_exp:`: those catalogues are the + `psf_validation` product `exp_persist` packs for every exposure whatever the + keep list says. Like `exp_persist` it runs whatever `clean:` and + `exposure_maps.nexp.enabled` say, because its input is on /scratch and the + purge takes it; `clean_exposure` takes its manifest as an input. + One further job, `nexp_map`, stamps every footprint into + `/nexp_map/nexp_map_.hsp` + — a uint16 map counting, per sky pixel, the exposures with a valid PSF there, + beside a record `nexp_map_.json` of which exposures it holds. There + is one count, of exposures: it is what sp_validation's `npoint >= 3` cut + reads (`notebooks/demo_apply_hsp_masks.py`). The job is + **campaign-cumulative**: its declared inputs are the in-scope footprints, but + the script reads *every* record on the products root, reclaimed exposures + included, so appending tiles grows the map instead of replacing it. Like the + defect map it is built whenever the campaign can support it — every + fitted-PSF run; `psf_model: fake` records no footprints and skips it — and + `exposure_maps.nexp.enabled: false` opts out: it is rebuilt whole rather + than reconciled, so a campaign appended in many small batches may prefer to + build it once at the end. Its memory is the map itself, sized on the + campaign's footprint like `defect_map_merge` (~48 GB of uint16 over the DR6 + footprint). Plotting stays out of the DAG: run `plot_coverage_map -i + /nexp_map/nexp_map_.hsp ...` by hand, with the sky + windows under `exposure_maps.nexp.plot` in `config.yaml`. - **External masks are wired, on the tile side only (data runs).** `inputs.masks` is a third input root beside tiles and exposures, set per machine in the `machines:` table, exported as `$SP_INPUT_MASKS` and @@ -378,6 +464,9 @@ profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; kee and an unknown *name* is a parse-time error listing the valid ones. The list is exposure-side only; tile-side retention is #844 follow-up. - **The campaign ends in two merged catalogues, and the workflow makes both.** + (Four campaign products, counting the two maps above — but those are maps + for the footprint, not catalogues, and nothing downstream of them lives + here.) Everything above is per unit; the two products downstream analysis actually opens are per *campaign*, and until these rules existed each was a manual pass after the run. diff --git a/workflow/Snakefile b/workflow/Snakefile index 42fdf4fbe..d1678f04e 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -31,8 +31,10 @@ the manifest says "this stage succeeded", the log says "here is what happened" (the contract is argued in completeness.py's docstring). """ +import configparser import fnmatch import functools +import glob import hashlib import json import os @@ -61,6 +63,16 @@ configfile: str(Path(workflow.snakefile).parent / "config.yaml") if os.environ.get("SP_RUN_CONFIG"): workflow.configfile(os.environ["SP_RUN_CONFIG"]) +# A retired block parses without complaint and does nothing, so a run config +# still carrying one would believe it had switched a map off (or resized it). +# Checked on the merged config, before machine defaults are applied. +_retired = run_config.retired(config) +if _retired: + raise WorkflowError( + "Retired config key(s), no longer read — move them: " + + "; ".join(f"{path} -> {repl}" for path, repl in _retired) + + f" ({_RUN_CONFIG}).") + # input_type picks where the pixels come from (and config/cfis vs # config/cfis_image_sims); everything after ingestion is the same chain. INPUT_TYPES = {"data", "image_sims"} @@ -459,6 +471,10 @@ MERGE_FINAL_HASH = ":".join(( path_hash(Path(workflow.basedir).parent / "scripts" / "python" / "create_final_cat.py"), path_hash(CONFIG_DIR / "final_cat.param"))) +DEFECT_HASH = script_hash("defect_map_exp.py") +MERGE_DEFECT_HASH = script_hash("merge_defect_map.py") +FOOTPRINT_HASH = script_hash("exp_footprint.py") +NEXP_MAP_HASH = script_hash("nexp_map.py") # ngmix_range.py earns a hash for a stronger reason than the others. What it # emits is not a stale RESULT but a stale BOUNDARY, and a tile's eight chunks are # a PARTITION of its object IDs: resume a tile across an edit to the split and @@ -630,6 +646,353 @@ for _entry in PERSIST_EXP: except KeyError as _exc: raise WorkflowError(f"config persist_exp: {_exc.args[0]}") +# --- exposure-level HealSparse maps ---------------------------------------- +# Two campaign products are derived exposure by exposure and combined once per +# campaign, both on the persistent root and both at ONE resolution: +# +# exp_defect_map -> defect_map_merge the instrument-flag defect map (#878): +# boolean, True = masked, OR of fragments +# exp_footprint -> nexp_map the exposure-count map: per sky pixel, +# how many exposures with a valid PSF +# model cover it +# +# THE RESOLUTION IS THE LADDER'S, not a choice made here: nside_sparse 131072 +# over nside_coverage 128, matching every map under `inputs.masks` and every +# entry in config_tile_Mc.ini's MASK_EXT_PATHS, so both maps align pixel-wise +# with the UNIONS bit masks and with each other. A coarser map looks entirely +# reasonable and does not align, and no consumer would notice. It is config +# (`exposure_maps:`) so that a campaign against a different ladder can say so, +# not so that one map can drift from the other. Checked here, at parse time: +# healsparse requires powers of two and says so only once a job has reached a +# node. +_MAPS = dict(config.get("exposure_maps") or {}) +MAP_NSIDE = int(_MAPS.get("nside", 131072)) +MAP_NSIDE_COVERAGE = int(_MAPS.get("nside_coverage", 128)) +for _key, _n in (("nside_coverage", MAP_NSIDE_COVERAGE), ("nside", MAP_NSIDE)): + if _n <= 0 or _n & (_n - 1): + raise WorkflowError( + f"config exposure_maps.{_key} must be a power of 2, got {_n}") +if MAP_NSIDE_COVERAGE >= MAP_NSIDE: + raise WorkflowError( + f"config exposure_maps: nside_coverage={MAP_NSIDE_COVERAGE} must be " + f"below nside={MAP_NSIDE}.") + +# THE DEFECT MAP. The exposures' flag images are the campaign's one masking +# input that never leaves the pixel domain, so the survey footprint — built from +# CCD corner WCS — cannot subtract them. exp_defect_map rasterizes them per +# exposure and defect_map_merge unions the fragments; both are argued in +# workflow/scripts/defect_map_exp.py and merge_defect_map.py. +# +# DATA RUNS ONLY. The image simulations ship a flag file per exposure +# (simu_flag-.fits, split and read by SExtractor like the survey's) but +# every pixel in it is zero — measured on 2079614 and 2079616 across three +# skills_out branches — so the union would be an empty map at the cost of a +# rasterization per exposure. defect_map_exposures() is where this gate acts; +# every defect-map target, input and sizing, and clean_exposure's edge, reads +# through it. +MAPS_DEFECTS = INPUT_TYPE == "data" + +DEFECT_OVERSAMPLE = int(dict(_MAPS.get("defect") or {}).get("oversample", 3)) +if DEFECT_OVERSAMPLE < 2: + raise WorkflowError( + f"config exposure_maps.defect.oversample={DEFECT_OVERSAMPLE}: 2 is the " + f"geometric floor (the four corners of a CCD pixel bound its " + f"footprint); 1 would sample pixel centres alone and lose thin defects " + f"entirely.") + +def _split_n_hdu(): + """How many CCDs exp_split writes, from ITS OWN config. + + Read rather than repeated, because it is the number defect_map_exp.py + checks the split dir against: a constant here that drifted from + config_exp_Sp.ini would either fail every exposure or restore the silent + half-exposure fragment the check exists to stop. It rides on the rule's + params, so changing N_HDU re-rasterizes. + """ + parser = configparser.ConfigParser() + parser.read(CONFIG_DIR / "config_exp_Sp.ini") + try: + return int(parser["SPLIT_EXP_RUNNER"]["N_HDU"]) + except (KeyError, ValueError) as exc: + raise WorkflowError( + f"config_exp_Sp.ini: no readable N_HDU ({exc}); exp_defect_map " + f"needs it to tell a complete split dir from a truncated one.") + + +DEFECT_N_CCDS = _split_n_hdu() + + +def defect_map(): + """The campaign's merged defect map. A PRODUCT, not an input: nothing in + this workflow opens it, and config_tile_Mc.ini deliberately does not name + it (that file's header carries the recipe for adding it after a campaign + has produced one).""" + return f"{PRODUCTS_DIR}/defect_map/defect_map_{CAMPAIGN}.hsp" + + +def defect_map_sidecar(): + """The record of which exposures are already in the map — what makes an + append cheap and a removal correct (merge_defect_map.py argues it).""" + return f"{PRODUCTS_DIR}/defect_map/defect_map_{CAMPAIGN}.json" + + +def prod_exp_fragment(exp): + """The healsparse fragment exp_defect_map writes. Not a declared output of + anything — the rule declares its manifest, exactly as exp_persist does.""" + return f"{prod_exp_dir(exp)}/defect/defect-{exp}.hsp" + + +@functools.lru_cache(maxsize=1) +def defect_map_exposures(): + """The exposures the union covers: every exposure of TILES_READY. + + SIMPLER THAN THE STAR CATALOGUE'S SET, and the difference is where the two + rules hang. exp_persist waits on exp_psf, so a reclaimed exposure's manifest + cannot be requested without rebuilding a four-hour chain from VOS — hence + star_cat_inputs()'s manifest-or-tar split. exp_defect_map waits on + exp_split, and its manifest lives on the persistent root, so an exposure + that already has one is a DAG leaf and an exposure that does not is one + split away. No special case is needed and none is invented. + + THE IDS, NOT THE PATHS: the fingerprint on `params` must move when the SET + moves and not when a path does. + + Empty for image simulations (MAPS_DEFECTS): their flag images are blank. + """ + if not MAPS_DEFECTS: + return [] + return sorted({e for t in TILES_READY for e in tile_exposures(t)}) + + +@functools.lru_cache(maxsize=1) +def defect_map_inputs(): + """What defect_map_merge waits for: the fragment manifests of the exposures + whose stores this invocation can still reach. + + THE SAME SPLIT AS star_cat_inputs(), AND FOR THE SAME REASON. A live + exposure is depended on through its MANIFEST, which is what orders the merge + after the rasterization. A reclaimed one is depended on through its + FRAGMENT — the .hsp, which is not a declared output of any rule, so a + fragment that exists is a DAG leaf snakemake requires and builds nothing + for. + + Naming the manifest for a reclaimed exposure would be the avalanche + persist_targets() drops such exposures to avoid: this rule's input is + exp_split's manifest, which went with the scratch store, so ANY rerun + trigger on exp_defect_map — and `params` is a live trigger, carrying the + script hash, the nside and the oversampling this file advertises as cheap to + change — would schedule exp_get_images and exp_split from VOS, a four-hour + chain per exposure, campaign-wide, on a one-character edit. The fragment is + already on the persistent root and cannot be rebuilt in place; requiring it + asks for nothing. + + An exposure reclaimed by a workflow PREDATING this rule has neither, and is + left out of the edge entirely; the job skips it and records it as missing. + """ + live, reclaimed = [], [] + for exp in defect_map_exposures(): + if not exp_store_reclaimed(exp): + live.append(prod_exp_manifest(exp, "exp_defect_map")) + elif Path(prod_exp_fragment(exp)).exists(): + reclaimed.append(prod_exp_fragment(exp)) + return live + reclaimed + + +def defect_map_targets(): + """The merged map, whenever this campaign has an exposure to rasterize. + + SILENCE IS THE WRONG ANSWER when there is none. Every exposure of a campaign + reclaimed by a workflow predating this rule is out of defect_map_inputs(), + so the DAG carries no defect job at all and an operator sees nothing — not + "no map is possible here". Say it once, at parse time, like the missing + index note above. + """ + if not workflow.is_main_process: + return [] + if not defect_map_inputs(): + if defect_map_exposures(): + logger.warning( + f"defect map: none of this campaign's " + f"{len(defect_map_exposures())} exposure(s) can contribute — " + f"each is reclaimed and has no fragment on the persistent root, " + f"so no fragment can be rasterized without rebuilding its chain " + f"from VOS. No defect map will be produced for this campaign " + f"(workflow/scripts/merge_defect_map.py argues why that is not " + f"an error).") + return [] + return [defect_map(), defect_map_sidecar()] + + +def defect_map_fragment_targets(): + """The per-exposure fragments `rule all` requests DIRECTLY. + + Same argument as persist_targets(): the fragment must leave /scratch whether + or not the campaign ever merges or reclaims, because the flag splits go with + the store at the purge. Reaching them only through defect_map_merge would + make a campaign that never merged lose them. + + It is defect_map_inputs() verbatim, so it inherits that function's split: + a live exposure's manifest, which is the thing to BUILD, and a reclaimed + one's existing fragment, which is already where it needs to be and asks for + nothing. + + HEAD PROCESS ONLY, for the same cost reason as persist_targets(). + """ + if not workflow.is_main_process: + return [] + return defect_map_inputs() + + +# THE EXPOSURE-COUNT MAP. exp_footprint records, per exposure, the four sky +# corners of every CCD with a valid PSF model; nexp_map stamps every record on +# the products root into one map. There is one count and it is of EXPOSURES: +# the CCDs of one exposure do not overlap, so stamping 1 per CCD polygon counts +# exposures per pixel, and it is what sp_validation's `npoint >= 3` cut reads. +# +# exp_footprint needs no precondition on `persist_exp:`: exp_persist packs every +# CCD's psf_validation catalogue whatever that list says (persist_exp.py's +# ALWAYS set), so its manifest always carries the valid-PSF set the record +# selects on. +# +# FITTED-PSF RUNS ONLY, by the same rule as MAPS_DEFECTS: the map is built +# whenever the campaign can support it. psf_model=fake (the image simulations' +# true PSF) fits no PSF model, so no exposure records a footprint +# (footprint_targets) and there is nothing to count; such a campaign skips the +# map rather than failing. `nexp.enabled: false` opts a fitted-PSF campaign +# out — the map is rebuilt whole, not reconciled, so a campaign appended in +# many small batches may want to build it once at the end; the per-exposure +# records are written either way. flag() because `--config` delivers booleans +# as strings. +_NEXP = dict(_MAPS.get("nexp") or {}) +MAPS_NEXP = PERSISTS_PSF and flag(_NEXP.get("enabled", True), default=True) + + +def nexp_map(): + """The campaign's exposure-count map; it carries run: so it identifies its + campaign outside the products root.""" + return f"{PRODUCTS_DIR}/nexp_map/nexp_map_{CAMPAIGN}.hsp" + + +def nexp_map_manifest(): + """The record of which exposures the map contains, beside it.""" + return f"{PRODUCTS_DIR}/nexp_map/nexp_map_{CAMPAIGN}.json" + + +def footprint_targets(): + """Which exposures this invocation must record a per-CCD sky footprint for. + + REQUESTED BY `rule all` DIRECTLY, AND NOT GATED ON `nexp.enabled`. The + record is derived from `headers-.npy`, which lives on /scratch and + which the 60-day purge takes whether or not this campaign ever builds a map + — so building it only when a map is asked for means that turning the map on + a month later finds the WCS gone and no way back short of re-downloading the + exposures. This is persist_targets' argument below, applied to the other + scratch-only input the map needs: milliseconds and a few KB now, against an + irrecoverable loss later. + + Scope, and the reclaimed-store exclusion, are persist_targets' exactly; + footprint_edge() below argues the exclusion. + + The exposures are psf_exposures(): the record is the valid-PSF CCD set, + read off exp_persist's manifest, so psf_model=fake (no PSF model, no + exp_persist) records none. + + HEAD PROCESS ONLY, for the same reason as clean_targets(). + """ + if not workflow.is_main_process: + return [] + return sorted(m for e in psf_exposures() for m in footprint_edge(e)) + + +def footprint_edge(exp): + """The footprint manifest to depend on for this exposure: itself while the + store is live, nothing once it is reclaimed. + + The one edge `rule all` and clean_exposure both take to an exposure's + footprint. A live exposure is asked for its manifest — the thing to build, + and what orders clean_exposure after the read of headers-.npy. + + A RECLAIMED exposure is asked for nothing. Its manifest is exp_footprint's + declared output, which hangs off exp_persist's manifest and through it off + exp_psf's, so naming it keeps the whole producer chain in the DAG: any + change upstream — a `persist_exp:` edit changes exp_persist's params — + reschedules exp_get_images, exp_split and exp_psf from VOS to rebuild a + record that is already on the persistent root. Unlike the defect fragment, + the record has no undeclared leaf to depend on instead; the manifest IS the + record. Nothing is lost by the omission: the store the edge ordered against + is already gone, and nexp_map reads every durable record whether or not it + is an edge. + + exp_store_reclaimed(), not the tombstone, is the test, for the purge case + its docstring names: a purged store has no tombstone, and a manifest asked + for there is a chain rebuilt, or, with no record written, an exp_footprint + job that fails reading headers that are gone. + """ + return ([] if exp_store_reclaimed(exp) + else [prod_exp_manifest(exp, "exp_footprint")]) + + +@functools.lru_cache(maxsize=1) +def nexp_map_exposures(): + """The exposure IDs nexp_map's map is built from — what its fingerprint is + taken over. + + nexp_map.py globs EVERY footprint record on the persistent root, but its + declared inputs are only footprint_targets(): a record of a reclaimed or + out-of-scope exposure is no edge, so without this set on `params` a record + added while the map was off, and then reclaimed, never reruns the map. The + set is the records already on the root plus the ones this invocation is + declared to write, which is exactly what the job will glob when it runs. + + IDS, NOT STAMPS, for unit_fingerprint's reason. A record written in this + invocation has no stamp at parse time and one at the next parse, and an + exposure moves from edge to glob-only when its store is reclaimed; a + fingerprint over mtimes or digests would move on both and rerun a + many-hour job over the same records. A declared record's content is + covered by its input edge instead, and a reclaimed exposure's record has + no producer left in this workflow to change it. + """ + # prod_exp_manifest("*", ...) is the record's path with the shard and the + # exposure id both "*": the same layout exp_footprint writes, as a glob. + durable = {Path(m).parent.parent.name + for m in glob.glob(prod_exp_manifest("*", "exp_footprint"))} + declared = {e for e in psf_exposures() if footprint_edge(e)} + return sorted(durable | declared) + + +def nexp_map_targets(): + """The campaign map, whenever the campaign fits a PSF (MAPS_NEXP) and has a + footprint record to stamp. + + NO RECORD, NO JOB. An exposure reclaimed by a workflow predating + exp_footprint has no record and cannot get one without rebuilding its chain + from VOS, so a resumed campaign can reach this point with nothing on the + root and nothing declared; requesting the map there is a job that fails + ("nothing to build a map from") on every invocation. + + SILENCE IS THE WRONG ANSWER for the in-scope exposures with no record: the + map is built without them and UNDERCOUNTS where they lie. Say how many, once, + at parse time, as defect_map_targets() does for its own lost exposures. + + HEAD PROCESS ONLY, for the reason clean_targets() gives: this feeds `rule + all` at module level, so it is evaluated on every per-job re-parse too, none + of which can schedule `all`. + """ + if not MAPS_NEXP or not workflow.is_main_process: + return [] + have = set(nexp_map_exposures()) + lost = [e for e in psf_exposures() if e not in have] + if lost: + logger.warning( + f"exposure-count map: {len(lost)} of this campaign's " + f"{len(psf_exposures())} exposure(s) have no footprint record — " + f"each was reclaimed before exp_footprint could read its headers — " + f"so the map undercounts wherever they lie" + + ("" if have else "; with no record at all, no map is built") + + ".") + if not have: + return [] + return [nexp_map(), nexp_map_manifest()] def persist_targets(): """Which exposures this invocation must pack PSF products off scratch for. @@ -942,6 +1305,27 @@ STAR_CAT_PRODUCT = _persist.ALWAYS STAR_CAT_PATTERN = _persist.resolve(_persist.ALWAYS) TILE_BYTES_DEFAULT = 46_000_000 +# THE TWO EXPOSURE-LEVEL MAPS. The campaign's footprint is what sizes each, +# not the exposure count: a map holds one block per coverage pixel the campaign +# touches, (nside/nside_coverage)^2 sparse pixels. The defect accumulator is one +# BIT per sparse pixel — 128 KiB per coverage pixel at the ladder's 131072/128 — +# and its loop holds one fragment (~2 MB) at a time; the exposure-count map is +# uint16, 2 MiB per coverage pixel, sixteen times the defect map's. Two ways to +# know the count, in order: the record the last build wrote beside its map, +# which is the campaign's own measured footprint; and before there is one, the +# exposures times the 13 coverage pixels ONE exposure touched (measured, +# 2079612p), capped at the survey FOOTPRINT — exposures overlap almost +# completely, so the per-exposure sum is only honest while the campaign is +# small, and unbounded it would ask for the whole sky (map_cov_pixels() carries +# the arithmetic). +DEFECT_MEM_BASE_MB = 800 # interpreter + astropy + healpy + healsparse +NEXP_MEM_BASE_MB = 800 +MAP_COV_PER_EXP = 13 +# The merge holds no partial state, so it must stay inside the walltime that +# Alliance policy lets a job run without checkpointing. A campaign that needs +# longer needs a resumable accumulator, not a bigger number here. +DEFECT_RUNTIME_CAP_MIN = 660 # 11 h + def _size(path, default): """Bytes on disk, or the documented per-unit default if it is not there.""" @@ -951,6 +1335,51 @@ def _size(path, default): return default +# The measured DR6 FOOTPRINT in nside-128 coverage pixels: ~23k, taken from the +# UNIONS ugriz maps staged under `inputs.masks`, which cover the same sky the +# exposures do. It is the ceiling on the pre-sidecar estimate below, and it is a +# FOOTPRINT rather than a whole sky: the survey is ~5000 deg^2, an eighth of the +# 12 * 128^2 = 196608 coverage pixels a full-sky cap would allow. +MAP_COV_FOOTPRINT = 23_000 + + +def map_cov_pixels(record, n_exposures): + """Coverage pixels a campaign map occupies: `n_coverage_pixels` from the + `record` the last build wrote beside it, else an estimate from the exposures. + + ONCE THERE IS A RECORD the count is the map's own, measured. Before that it + is a prior, and the prior must be a FOOTPRINT estimate, not a sum over + exposures. Exposures overlap almost completely — each tile is covered by 7-10 + of them and they tile the same sky — so MAP_COV_PER_EXP * n_exposures + passes the survey footprint at ~1800 exposures and, uncapped at the full sky, + asks for 25.8 GB of defect accumulator (and ~52 GB of rule, doubled again on + a retry) for a job the merge measures at ~3 GB. Capping at the measured DR6 + footprint keeps the first-ever build of a large campaign asking for what it + needs; a small campaign is still sized on its own exposures, where the + per-exposure figure is the honest one. + """ + try: + return int(json.loads(Path(record).read_text())["n_coverage_pixels"]) + except (OSError, ValueError, KeyError): + return min(MAP_COV_PER_EXP * n_exposures, MAP_COV_FOOTPRINT) + + +def defect_map_cov_bytes(): + """Bytes the union's accumulator occupies: coverage pixels x the bit-packed + block size.""" + block = (MAP_NSIDE // MAP_NSIDE_COVERAGE) ** 2 // 8 + return map_cov_pixels(defect_map_sidecar(), + len(defect_map_exposures())) * block + + +def nexp_map_cov_bytes(): + """Bytes the exposure-count map occupies: coverage pixels x the uint16 + block size — ~48 GB over the DR6 footprint.""" + block = (MAP_NSIDE // MAP_NSIDE_COVERAGE) ** 2 * 2 + return map_cov_pixels(nexp_map_manifest(), + len(nexp_map_exposures())) * block + + def star_cat_max_bytes(): """The LARGEST exposure's psf_validation members — what sizes the merge. @@ -1240,26 +1669,35 @@ include: "rules/tile.smk" # they would be ~20k (clean_exposure) or ~23k (clean_tile) sbatch submissions at # DR6 scale for work shorter than the scheduling latency. # -# Both are DAG LEAVES, so neither constrains a `group:` label. (A mid-chain -# 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.) +# Both clean rules are DAG LEAVES, so neither constrains a `group:` label. (A +# mid-chain 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.) # # 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 +# latency that would submit it (exposure.smk argues the placement in full). +# exp_footprint sits beside it and is smaller still: one pickle load and ~160 +# pixel_to_world calls. Both sit 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 they add no new grouping constraint. +# +# The two map merges are deliberately NOT local: defect_map_merge holds the +# campaign's footprint as a bit accumulator, and nexp_map stamps ~1M polygons at +# nside=131072 (exposure.smk). +localrules: all, prepare_all_tiles, clean_exposure, clean_tile, exp_persist, exp_footprint rule all: input: [final_cat(t) for t in TILES_READY], persist_targets(), + defect_map_fragment_targets(), + footprint_targets(), star_cat_targets(), final_cat_targets(), + defect_map_targets(), + nexp_map_targets(), clean_targets(), clean_tile_targets(), diff --git a/workflow/config.yaml b/workflow/config.yaml index 1d4ebd49c..02598259f 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -180,6 +180,107 @@ machines: persist_exp: - psf_model +# EXPOSURE-LEVEL HEALSPARSE MAPS. Two campaign products are built exposure by +# exposure and combined once per campaign, on products_dir: +# +# defect /defect_map/defect_map_.hsp boolean, True = masked +# nexp /nexp_map/nexp_map_.hsp uint16 exposure count +# +# nside / nside_coverage are SHARED, and they are the mask ladder's: 131072 over +# 128 (~1.6"/pixel) matches every map under `inputs.masks` and every entry in +# config_tile_Mc.ini's MASK_EXT_PATHS, so both maps drop into that ladder, and +# align with each other, without a resolution change. A coarser map looks +# entirely reasonable and does not align, and no consumer would notice. They are +# config so a campaign against a different ladder can say so, not so one map can +# drift from the other. Both must be powers of two, nside_coverage below nside; +# the Snakefile checks at parse time. +# +# Neither map is an input: nothing in the workflow reads either back. In +# particular config_tile_Mc.ini's MASK_EXT_PATHS names neither -- that ladder +# names maps that exist before the run starts, and these exist only after it. +# That file's header carries the recipe for adding the defect map once a +# campaign has produced one. +# +# DEFECT (CosmoStat/shapepipe#878). The exposures' flag images -- bad columns, +# saturated pixels, bleed trails -- are the one masking input this campaign has +# that never leaves the pixel domain: exp_split splits them per CCD and +# SExtractor reads them as IMAFLAGS_ISO, and that is the end of it. The survey +# footprint is built from the CCD corner WCS in the headers, so it cannot +# subtract them and would silently include defective pixels. So +# `exp_defect_map` rasterizes each exposure's flags into a boolean healsparse +# fragment on products_dir, at +# /exp///defect/defect-.hsp -- its own subdir +# beside the exposure's manifests/ and its PSF tar, not beside the tar itself -- +# and `defect_map_merge` unions the campaign's fragments. Data runs only: the +# image simulations' flag images are blank (MAPS_DEFECTS, Snakefile). +# +# Measured on one real exposure (2079612p, 40 CCDs, 4.1% of pixels flagged) on +# this login node inside the campaign container: 34 s at oversample 3, peak RSS +# 0.62 GB, 359858 healpix pixels, a 2.0 MB fragment over 13 coverage pixels +# (21 s / 0.41 GB / 357896 pixels at oversample 2; 95 s / 1.29 GB / 360840 at +# oversample 5 -- so 3 is within 0.3% of 5 over the whole exposure, at a third +# of the cost). The rasterization is BATCHED, so peak memory is the batch and +# not the exposure: a fully flagged MegaCam chip -- 9.4M pixels, the worst case +# there is, and one a real exposure carries whenever a chip is dead -- measures +# 0.74 GB and 18 s on its own. At DR6 scale (~20k exposures) the products are +# ~40 GB and ~60k inodes, +# against a ~1 M-inode group quota: THREE per exposure, not one -- the defect/ +# directory, the fragment in it, and the manifest under manifests/. +# +# `defect.oversample` is samples per CCD pixel per axis, endpoints included, so +# 2 means the four corners -- the geometric floor, since a healpix pixel at +# nside 131072 (1.61") is ~74 times a MegaCam pixel (0.187") and cannot hide +# inside one. Raising it only fills boundary and rounding gaps, at linear cost. +# Measured on 2079612p CCD 0 (430886 flagged pixels) against a 12x12 interior +# reference of 14229 healpix pixels: +# +# oversample samples/pixel healpix pixels missed vs reference +# 2 (corners) 4 14081 169 (1.2%) +# 3 9 14198 54 (0.4%) +# 5 25 14254 5 (0.04%) +# +# The rasterization is CONSERVATIVE by design -- a healpix pixel is masked if +# any part of the flagged region touches it -- because the centre-based +# alternative erases a one-pixel bad column, which is the geometry #878 exists +# to keep. +# +# NEXP, the exposure-count map: per sky pixel, the number of exposures with a +# VALID PSF MODEL covering it -- the count sp_validation's `npoint >= 3` cut +# reads. `exp_footprint` records, per exposure, the sky corners of every CCD +# with a PSF validation catalogue (psfex_interp writes validation_psf-*.fits on +# success and returns WITHOUT writing on NOT_ENOUGH_STARS / BAD_CHI2 / +# FILE_NOT_FOUND); exp_persist packs those catalogues whatever `persist_exp:` +# says, so there is nothing to configure for the records, and they are written +# all along whatever `enabled` says, because their input is on /scratch and the +# purge takes it. ONE campaign-level job, `nexp_map`, stamps every footprint on +# the products root into the map. Fitted-PSF runs only, like the defect map's +# data-only rule: psf_model: fake writes no PSF model and records no +# footprint, so an image-sims campaign on the true PSF skips the map. +# `nexp.enabled: false` opts a fitted-PSF campaign out; the map is rebuilt whole +# on every invocation that adds records, so a campaign appended in many small +# batches may prefer to build it once at the end. +# +# Plot windows are the two UNIONS sky regions, with the colourbar clipped to the +# 1-5 exposure range. Plotting stays OUT of the DAG (a human act on a durable +# product) and nothing reads `plot:`; it records the windows to hand +# `plot_coverage_map`. +# @sc [decision:masking.defect_map_from_flags] +exposure_maps: + nside: 131072 + nside_coverage: 128 + defect: + oversample: 3 + # @sc [decision:masking.nexp_map_valid_psf_ccds] + nexp: + enabled: true + plot: + colorbar: true + n_exp_min: 1 + n_exp_max: 5 + regions: + SGC: {ra_min: -20, ra_max: 45, dec_min: 18, dec_max: 40} + NGC: {ra_min: 110, ra_max: 270, dec_min: 28, dec_max: 90} + # Per-exposure store reclamation once every tile reading it has its vignets. clean: true diff --git a/workflow/config/cfis/config_tile_Mc.ini b/workflow/config/cfis/config_tile_Mc.ini index a92b099a6..c3c2f7b36 100644 --- a/workflow/config/cfis/config_tile_Mc.ini +++ b/workflow/config/cfis/config_tile_Mc.ini @@ -122,5 +122,33 @@ N_EPOCH_SLOTS = 12 # PSF-star selection deliberately does NOT consume these: MASK_PATHS in # config_exp_psfex.ini stays commented out, keeping the star diet narrow # (instrument flags only). See that file's header. +# +# THE CAMPAIGN'S OWN DEFECT MAP IS NOT IN THE LADDER BELOW, AND THAT IS +# DELIBERATE (CosmoStat/shapepipe#878). The workflow emits one more boolean +# healsparse map in exactly this form -- nside 131072, coverage 128, True = +# masked -- from the exposures' instrument flag images: bad columns, saturated +# pixels and bleed trails, the one masking input that otherwise never leaves the +# pixel domain. But it is a campaign PRODUCT, not an input. Every path above +# exists before the run starts; that one exists only after it, so wiring it here +# would name a file the first run of a fresh campaign cannot have, and every +# tile would fail on a missing map. +# +# TO ADD IT, AFTER A CAMPAIGN HAS PRODUCED ONE. It lands at +# /defect_map/defect_map_.hsp, being `run:` in the run +# config. Copy or symlink it into the `inputs.masks` root -- the same root +# everything above resolves through, so the ladder keeps one location -- and +# append one entry: +# +# MASK_EXT_PATHS = ..., defect:$SP_INPUT_MASKS/defect_map_.hsp +# +# then add the matching `MASK_defect` line to final_cat.param, which is what +# turns the query into a column. The label is free; `defect` reads better than +# an `n` name, since this map is not one bit of the UNIONS ladder. +# +# TWO THINGS TO KNOW BEFORE YOU DO. It is the campaign's OWN footprint, so a +# catalogue built from a different tile list reads False -- not clean, unknown -- +# wherever that campaign had no exposure; and the rasterization is conservative, +# widening a one-pixel bad column to the 1.61" healpix resolution, which is the +# price of keeping thin defects at all. Nothing here cuts on it either way. # @sc [decision:masking.sky_mask_application] MASK_EXT_PATHS = n1:$SP_INPUT_MASKS/mask_ugriz_nside131072_n1.hsp, n2:$SP_INPUT_MASKS/mask_ugriz_nside131072_n2.hsp, n4:$SP_INPUT_MASKS/mask_ugriz_nside131072_n4.hsp, n8:$SP_INPUT_MASKS/mask_ugriz_nside131072_n8.hsp, n16:$SP_INPUT_MASKS/mask_ugriz_nside131072_n16.hsp, n32:$SP_INPUT_MASKS/mask_ugriz_nside131072_n32.hsp, n64:$SP_INPUT_MASKS/mask_ugriz_nside131072_n64.hsp, n128:$SP_INPUT_MASKS/mask_ugriz_nside131072_n128.hsp, n256:$SP_INPUT_MASKS/mask_ugriz_nside131072_n256.hsp, n1024:$SP_INPUT_MASKS/mask_ugriz_nside131072_n1024.hsp, n2048:$SP_INPUT_MASKS/mask_ugriz_nside131072_n2048.hsp diff --git a/workflow/rules/exposure.smk b/workflow/rules/exposure.smk index 221d10928..6a4fb275d 100644 --- a/workflow/rules/exposure.smk +++ b/workflow/rules/exposure.smk @@ -1,16 +1,23 @@ """Exposure chain — per exposure, keyed by exp base id (dedup is structural). - exp_get_images -> exp_split -> exp_psf -> exp_persist + exp_get_images -> exp_split -> exp_psf -> exp_persist -> exp_footprint + `-> exp_defect_map Each in the exposure's own sharded work dir, chained by manifests; every config reads fixed ``$SP_RUN/output/run_sp_exp_*`` INPUT_DIRs, so nothing resolves a run log. There is no `prepare_exposures` aggregation target: these chains hang off the compute DAG (`all` <- final_cat <- tile chain <- exposure manifests). -NO MASK RULE, and that is the design (PR #847). ShapePipe generates no masks. -The only mask that reaches pixels is the instrument flag image delivered with -the exposure, which ``exp_split`` splits per CCD alongside image and weight and -SExtractor reads directly. Sky-fixed masks are healsparse maps, queried once per +NO MASK-GENERATION RULE, and that is the design (PR #847). ShapePipe generates +no masks. The only mask that reaches pixels is the instrument flag image +delivered with the exposure, which ``exp_split`` splits per CCD alongside image +and weight and SExtractor reads directly. ``exp_defect_map`` does not generate +that mask, it EXPORTS it: the flag image is the one masking input the campaign +has that never becomes sky-fixed, so the survey footprint cannot subtract it +(#878), and the rule rasterizes it into the same healsparse form as every other +mask here. Nothing in this workflow reads the result back. + +Sky-fixed masks are healsparse maps, queried once per object: ``mask_query`` (inside exp_psf's config chain) writes ``MASK_EXT`` onto each CCD's SExtractor catalogue, carried for transparency and measurement (selection's only mask cut is ``IMAFLAGS_ISO``; imposing ``MASK_EXT`` is @@ -19,12 +26,15 @@ opt-in, see ``star_selection.setools``), and ``make_cat`` writes the per-band 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. +Three per-exposure rules here write to the PERSISTENT root, before the purge +(or clean_exposure) can take their inputs off /scratch. ``exp_persist`` packs +the PSF products named by `persist_exp:` into one tar per exposure. 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. ``exp_footprint`` and ``exp_defect_map`` are +the per-exposure halves of the two exposure-level HealSparse maps — the +exposure-count map and the defect map — whose campaign halves, ``nexp_map`` and +``defect_map_merge``, close this file. 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, @@ -178,6 +188,144 @@ rule exp_persist: " {params.patterns}" +# --- per-CCD sky footprints (the exposure-count map's raw material) --------- +# One JSON per exposure recording, for each CCD that HAS a PSF model, the four +# sky corners of that CCD. nexp_map (below) stamps every such record into the +# campaign's exposure-count map. The record is a local read of two things this +# workflow already wrote; workflow/scripts/exp_footprint.py argues both inputs +# and the CCD-index invariant tests/unit/test_exp_footprint.py pins. +# +# A LOCALRULE (declared in the Snakefile), by exp_persist's arithmetic and then +# some: one pickle load and ~160 pixel_to_world calls, milliseconds, against a +# scheduling latency of seconds and ~20k exposures at DR6 scale. +# +# ONE DECLARED INPUT, AND IT IS THE PERSIST MANIFEST — deliberately NOT +# exp_psf's as well, even though the WCS array this rule reads is written by +# exp_split and lives in the same scratch store. exp_persist already orders this +# rule after the whole PSF chain, so the second edge would buy no ordering; what +# it WOULD buy is a scratch manifest (the one clean_exposure deletes) in the +# input list of a rule whose output is durable. A persistent-root manifest +# outliving a purged scratch store is a real state, and in it that edge would +# schedule a four-hour VOS rebuild of the exposure to satisfy a few KB of +# provenance. Declaring only the durable input makes the same state a loud +# one-line failure from the script instead. +# +# WRITES TO THE PERSISTENT ROOT, beside exp_persist's manifest and for the same +# reason: the record must outlive both reclamation and the purge — the map is +# campaign-cumulative, so a record written today is read by every map built +# after it. +# @sc [decision:masking.nexp_map_valid_psf_ccds] +rule exp_footprint: + input: + lambda wc: prod_exp_manifest(wc.exp, "exp_persist") + output: + manifest = f"{PROD_EXP_DIR}/manifests/exp_footprint.json" + # No `log:`, for exp_persist's reason: the failure modes are "no WCS array" + # and "the two stages disagree about the focal plane", both of which the + # script reports on stderr and neither of which has a per-CCD verdict. + params: + exp_dir = lambda wc: exp_dir(wc.exp), + persist = lambda wc: prod_exp_manifest(wc.exp, "exp_persist"), + script_hash = FOOTPRINT_HASH + threads: 1 + retries: 2 + resources: + mem_mb = 2000, + runtime = 10 + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/exp_footprint.py" + " --exp-dir '{params.exp_dir}' --exp {wildcards.exp}" + " --persist-manifest '{params.persist}'" + " --manifest {output.manifest}" + + +# --- the pixel-domain mask leaves the pixel domain (#878) ------------------- +# Another rule here that writes to the PERSISTENT root, and another one +# clean_exposure must wait for. It rasterizes this exposure's per-CCD instrument +# flag splits — the ONE masking input the campaign has that never becomes a +# sky-fixed map — into a boolean healsparse fragment at the mask ladder's own +# resolution, so the footprint can finally subtract bad columns, saturated +# pixels and bleed trails. The full argument, the WCS source and the measured +# oversampling table are in workflow/scripts/defect_map_exp.py. +# +# AFTER exp_split, NOT AFTER exp_psf, and it is deliberately not chained behind +# the PSF work: the flag splits exist the moment the split finishes, and hanging +# a two-minute rasterization off a four-hour rule would make re-rasterizing the +# campaign at a different oversampling cost the PSF chain. It runs in parallel +# with exp_psf, and both are ordered before clean_exposure. +# +# ONE DECLARED OUTPUT, AND IT IS A MANIFEST, exactly as exp_persist: the +# fragment is written beside it on the persistent root and the manifest records +# the per-CCD healpix counts. Byte-stable, so a no-op rerun does not move the +# mtime clean_exposure reads. +# +# NOT A LOCALRULE, and this is where it parts company with exp_persist. That +# rule is a tar of a few MB — seconds, far shorter than the scheduling latency. +# This one is 40 CCDs of WCS transforms and ang2pix over ~16 M flagged pixels: +# 34 s measured end to end on a real exposure at oversample 3 (2079612p, on this +# login node, inside the campaign container; 21 s at oversample 2). That is real +# work, it is CPU-bound, and running ~20k of them under local-cores in the head +# process would serialise the campaign behind them. +# +# THE RESOLUTION AND THE OVERSAMPLING RIDE ON params. Both are what the fragment +# IS, and both are decisions that may be revisited; on params they re-rasterize +# (a minute) and leave the split and the PSF chain alone. +# @sc [decision:masking.defect_map_from_flags] +rule exp_defect_map: + input: + rules.exp_split.output.manifest + output: + manifest = f"{PROD_EXP_DIR}/manifests/exp_defect_map.json" + # No `log:`, for exp_persist's reason: the failure modes are "no flag split + # under the store", "fewer flag splits than N_HDU" and "a flag split with no + # image beside it", all reported on stderr, none with a per-CCD verdict + # worth a completeness record. + params: + exp_dir = lambda wc: exp_dir(wc.exp), + dest = lambda wc: f"{prod_exp_dir(wc.exp)}/defect", + # config_exp_Sp.ini's own N_HDU, read at parse time. A split dir short + # of it is a hard error, not a smaller fragment: half an exposure's + # defects, written "complete", is a hole in the footprint nothing + # downstream can see (defect_map_exp.py's ccd_files argues it). + n_ccds = DEFECT_N_CCDS, + nside = MAP_NSIDE, + nside_cov = MAP_NSIDE_COVERAGE, + oversample = DEFECT_OVERSAMPLE, + script_hash = DEFECT_HASH + threads: 1 + retries: 2 + resources: + # Measured on 2079612p inside the container: peak RSS 0.62 GB at + # oversample 3 (0.41 GB at 2), dominated by one BATCH of a CCD's sample + # arrays plus the bit-packed fragment's 13 coverage pixels. + # + # FLAT IN BOTH DIRECTIONS THAT COULD BLOW IT: in the number of CCDs, + # which are rasterized one at a time, and in how badly any one of them + # is flagged, which is batched at defect_map_exp.CHUNK source pixels. + # The second is the one worth requesting for — a MegaCam exposure + # routinely carries a dead or saturated chip, 9.4M flagged pixels, and + # unbatched that is several GB and three OOMs, after which the exposure + # has no fragment AND cannot be reclaimed (clean_exposure waits on this + # manifest). Measured on exactly that case, a fully flagged chip against + # a real WCS: 0.74 GB. So 2000 covers the worst CCD at ~2.7x, not the + # measured average at ~3x. It scales with `oversample`, which is why + # that knob is not free. + mem_mb = lambda wc, attempt: 2000 * attempt, + # 34 s measured end to end on 2079612p at oversample 3; a factor of + # ~35 for a dirtier exposure (a fully flagged chip is 18 s on its own), + # a higher oversampling and a busy filesystem. + runtime = 20 + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/defect_map_exp.py" + " --exp-dir '{params.exp_dir}' --exp {wildcards.exp}" + " --dest '{params.dest}' --manifest {output.manifest}" + " --n-ccds {params.n_ccds}" + " --nside {params.nside} --nside-coverage {params.nside_cov}" + " --oversample {params.oversample}" + + # --- reclamation (D5) ------------------------------------------------------- # The one exception to "no reclamation in this file": clean_exposure OWNS # exposure-level deletion, and it is a real job, not temp() bookkeeping, because @@ -228,7 +376,40 @@ rule clean_exposure: else [prod_exp_manifest(wc.exp, "exp_persist")] if not exp_store_reclaimed(wc.exp) else [prod_exp_tar(wc.exp)] - if Path(prod_exp_tar(wc.exp)).exists() else []) + if Path(prod_exp_tar(wc.exp)).exists() else []), + # The fragment must be off /scratch before the store goes too — the flag + # splits go with it — but this edge is CONDITIONAL, and it is exactly + # defect_map_inputs()' split (Snakefile), for exactly its reason. + # + # Naming the exp_defect_map manifest unconditionally reopens the + # avalanche that function is written to avoid. An exposure whose store + # went to the 60-day /scratch purge, or to a `clean: false` run, has NO + # TOMBSTONE — exp_store_reclaimed()'s docstring names that case — so + # clean_targets() still asks for one, and the manifest it would then + # require sits behind exp_split's manifest, which went with the store: + # snakemake schedules exp_get_images and exp_split from VOS, four hours + # per exposure, campaign-wide, on the first run of this branch. + # + # So: a LIVE exposure is asked for its manifest (the thing to build, and + # the thing that orders this rule after the rasterization); a RECLAIMED + # one is asked for its FRAGMENT if it has one — already on the + # persistent root, no rule's declared output, hence a leaf that requires + # nothing — and for nothing at all if it has neither, which is an + # exposure reclaimed by a workflow predating this rule and whose flags + # are gone either way. Blocking its tombstone would pin its scratch + # store forever without recovering a single flag. Image simulations + # rasterize nothing (MAPS_DEFECTS, Snakefile), so they wait on nothing. + lambda wc: ([] if not MAPS_DEFECTS + else [prod_exp_manifest(wc.exp, "exp_defect_map")] + if not exp_store_reclaimed(wc.exp) + else [prod_exp_fragment(wc.exp)] + if Path(prod_exp_fragment(wc.exp)).exists() else []), + # And the footprint, for exp_persist's ordering reason one layer further + # out: it is derived from headers-.npy, which lives in the store + # this job deletes. Reclamation must not overtake the read; a reclaimed + # store has no read left to order against, and naming its footprint + # reopens the exposure chain: footprint_edge() (Snakefile). + lambda wc: (footprint_edge(wc.exp) if PERSISTS_PSF else []), output: tombstone = f"{EXP_DIR}/cleaned.json" params: @@ -325,3 +506,155 @@ rule star_cat_merge: " --output {output.star_cat}" " --campaign '{params.campaign}'" " --snapshot-json '{params.snapshot}'" + + +# --- the campaign's defect map (#878) --------------------------------------- +# ONE job per campaign: every exposure's fragment, OR-ed into +# `/defect_map/defect_map_.hsp`. It is the exposure side's third +# campaign product, and the only one nothing downstream in this workflow opens: +# it exists so the survey footprint can subtract the pixel-domain masking the +# CCD corner WCS cannot see. merge_defect_map.py argues the reconcile, the +# memory-flat accumulation and why the map is a campaign PRODUCT rather than an +# entry in config_tile_Mc.ini's MASK_EXT_PATHS. +# +# TWO DECLARED OUTPUTS, the map and the SIDECAR that records which exposures are +# already in it. They are written together and they are only meaningful +# together: a union cannot be un-OR-ed, so the record of what went in is what +# makes an append cheap and a removal correct. Declaring both means a failed job +# takes both, and the next invocation rebuilds rather than reconciling against a +# record that does not describe the map beside it. +# +# THE INPUT IS THE FRAGMENT MANIFESTS, not the fragments: the manifest is what +# exp_defect_map declares, so it is the edge that orders this after the +# rasterizations. Reclaimed exposures need no special case here — unlike the +# star catalogue's tars, the fragment manifest lives on the persistent root and +# survives reclamation, and the rule that writes it hangs off exp_split, not off +# a store that reclamation took. An exposure cleaned by a workflow that predates +# this rule simply has no fragment; the job skips it and says so. +# +# THE PATHS DO NOT REACH THE SHELL (~20k of them at DR6 scale, an order of +# magnitude over MAX_ARG_STRLEN): the job is handed the tile list and the index +# and derives the same set, with the set's FINGERPRINT on `params` as the rerun +# trigger. Same discipline as the two merges above. +# +# NOT A LOCALRULE: the accumulator is the campaign's footprint at nside 131072, +# ~3 GB resident at DR6 scale. +# @sc [decision:masking.defect_map_from_flags] +rule defect_map_merge: + input: + lambda wc: defect_map_inputs() + output: + defect_map = defect_map(), + sidecar = defect_map_sidecar() + params: + products_dir = str(PRODUCTS_DIR), + tile_list = str(config["tile_list"]), + index_db = str(INDEX_DB), + nside = MAP_NSIDE, + nside_cov = MAP_NSIDE_COVERAGE, + inputs = unit_fingerprint(defect_map_exposures()), + script_hash = MERGE_DEFECT_HASH + threads: 1 + # Declared so the attempt scaling above is not dead code. One retry, not the + # two the exposure rules take: a failed attempt here has already cost hours, + # both declared outputs go with it, and the retry starts from an empty + # accumulator — there is nothing to salvage and little to gain from a third. + retries: 1 + resources: + # Sized on the FOOTPRINT, not on the exposure count: the accumulator is + # one bit per sparse pixel of every touched coverage pixel, and the loop + # holds one fragment at a time (merge_defect_map.py's memory argument). + # defect_map_cov_bytes() is that arithmetic, from the sidecar's own + # recorded coverage count once there is one and, before that, from a + # per-exposure figure capped at the MEASURED DR6 footprint; the + # Snakefile's sizing block carries both and argues the cap. + # capped_mem() for the same reason star_cat_merge and final_cat_merge + # take it: defect_map_cov_bytes() scales as nside^2 and exposure_maps.nside + # is an advertised knob, so one ladder change turns this into a request + # no partition can schedule — a job that sits PENDING while the campaign + # looks alive, instead of a diagnosable OOM and a parse-time warning. + mem_mb = lambda wc, attempt: capped_mem( + attempt * (DEFECT_MEM_BASE_MB + + 2 * defect_map_cov_bytes() // 1_000_000), + "defect_map_merge"), + # Dominated by reading fragments (~2 MB each) and setting their pixels; + # ~1 s per exposure measured, over a floor that covers writing the map. + # + # AND IT IS THE EXPENSIVE CASE THAT SETS IT. An append reads exactly the + # appended exposures and finishes in minutes; a REBUILD — any exposure + # leaving the campaign, any fragment changed — reads all of them, ~40 + # GB at DR6 scale, and that is the ~6 h this formula sizes for at 20k + # exposures. Capped below 12 h because the job holds no partial state + # and Alliance policy asks anything longer to checkpoint: it cannot, so + # it must not ask. If a campaign ever needs longer than DEFECT_RUNTIME_ + # CAP_MIN, the accumulator has to become resumable (write map and + # sidecar every N fragments) rather than the cap being raised. + runtime = lambda wc, attempt: min( + attempt * (20 + len(defect_map_exposures()) // 60), + DEFECT_RUNTIME_CAP_MIN) + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/merge_defect_map.py" + " --products-dir '{params.products_dir}'" + " --tile-list '{params.tile_list}' --index-db '{params.index_db}'" + " --output {output.defect_map} --sidecar {output.sidecar}" + " --nside {params.nside} --nside-coverage {params.nside_cov}" + + +# --- the campaign's exposure-count map -------------------------------------- +# ONE job per campaign: every exposure footprint on the persistent root, stamped +# into `/nexp_map/nexp_map_.hsp` — per sky pixel, the +# number of exposures with a valid PSF model covering it. Fitted-PSF runs only +# (MAPS_NEXP, Snakefile, which argues the gate and the opt-out). +# +# CAMPAIGN-CUMULATIVE. The declared inputs are the IN-SCOPE footprint manifests +# of live stores — ordering, and reruns when one is rewritten, without dragging +# out-of-scope tiles or reclaimed exposures' chains into the DAG. The SCRIPT then +# reads every footprint record on the persistent root, reclaimed exposures +# included: their records outlive their scratch stores and are still valid sky. +# `params.footprints` fingerprints that whole set by exposure id +# (nexp_map_exposures(), Snakefile), so a record that reaches the root without +# being an edge still reruns the map. Appending tiles grows the map instead of +# replacing it, and rebuilding is one job rather than a campaign. +# +# REBUILT, NOT RECONCILED, unlike defect_map_merge: a count can be appended to +# but a changed footprint would have to be subtracted, and the records are small +# enough that restamping all of them is the simpler correct answer. +# +# NOT A LOCALRULE. At DR6 scale this stamps ~1M polygons at nside=131072 in a +# Python loop (coverage_map_builder.build_map). Memory is the map itself — +# uint16 over the campaign's footprint (nexp_map_cov_bytes(), Snakefile) — and +# the runtime is a first sizing, unmeasured above a few thousand CCDs. +# +# PLOTS STAY OUT OF THE DAG. `plot_coverage_map -i ...`, by hand, with +# the windows in config.yaml's `exposure_maps.nexp.plot` block — the same +# argument that keeps run_report.py a standalone script: it is a human act on a +# durable product. +# @sc [decision:masking.nexp_map_valid_psf_ccds] +rule nexp_map: + input: + footprint_targets() + output: + hsp = nexp_map(), + manifest = nexp_map_manifest() + # No `log:`: this is one job with one verdict, and its stderr is the job's. + params: + products = str(PRODUCTS_DIR), + nside_coverage = MAP_NSIDE_COVERAGE, + nside = MAP_NSIDE, + footprints = lambda wc: unit_fingerprint(nexp_map_exposures()), + script_hash = NEXP_MAP_HASH + threads: 1 + retries: 1 + resources: + mem_mb = lambda wc, attempt: capped_mem( + attempt * (NEXP_MEM_BASE_MB + + 2 * nexp_map_cov_bytes() // 1_000_000), + "nexp_map"), + runtime = 240 + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/nexp_map.py" + " --products-dir '{params.products}'" + " --out {output.hsp} --manifest {output.manifest}" + " --nside-coverage {params.nside_coverage} --nside {params.nside}" diff --git a/workflow/scripts/defect_map_exp.py b/workflow/scripts/defect_map_exp.py new file mode 100644 index 000000000..7784c300f --- /dev/null +++ b/workflow/scripts/defect_map_exp.py @@ -0,0 +1,367 @@ +#!/usr/bin/env python3 +"""Rasterize ONE exposure's instrument flags into a boolean healsparse fragment. + +Run as the shell of the in-DAG ``exp_defect_map`` rule, never by hand. + +WHY THIS EXISTS (CosmoStat/shapepipe#878). The survey footprint is built as a +positive coverage map minus the healsparse mask bits, and every masking input +the campaign has is already sky-fixed and queried per object — except one. The +instrument flag image delivered with each exposure (bad columns, saturated +pixels, bleed trails) never leaves the PIXEL domain: ``exp_split`` splits it per +CCD, SExtractor reads it as ``IMAFLAGS_ISO``, and that is the end of it. The +coverage map is built from the CCD corner WCS in the headers, so it cannot +subtract those pixels and the footprint silently includes them. The lost area is +percent-level, but it is exactly the thin, small-scale geometry an accurate +window function needs. Only ShapePipe ever opens these files, so the map has to +come from here. + +WHAT IT WRITES. ``/defect-.hsp``: a boolean ``HealSparseMap``, +``nside_sparse`` 131072 over ``nside_coverage`` 128, ``bit_packed``, ``True`` +where any flag bit is set. Those are not free choices — they are the +convention every other map in this campaign's ladder already follows (the +UNIONS ugriz bit maps under ``inputs.masks``, and ``config_tile_Mc.ini``'s +``MASK_EXT_PATHS``), so the fragments and the map they union into drop into that +ladder without a resolution change. ``True = masked`` likewise. + +ANY BIT, NOT A BIT TABLE. The flag image is a bitmask, but the CFIS flags are +"this pixel is not to be trusted" in several flavours (the campaign's own +exposures carry values 1, 2, 3, 8 and 11), and nothing downstream distinguishes +them: SExtractor's ``IMAFLAGS_ISO`` cut is nonzero-vs-zero. So the fragment is +``!= 0`` and the map is boolean. A per-bit ladder would be a different product +answering a question nobody has asked yet. + +WHERE THE WCS COMES FROM, AND WHY NOT FROM THE FLAG FILE. The flag mosaic's +per-CCD HDUs carry NO WCS at all — checked on a real exposure: ``CTYPE`` empty, +``CRVAL`` 0, the identity transform. Only the image mosaic is astrometric, so +``split_exp`` saves headers for the image suffix alone. The fragment therefore +takes each CCD's WCS from its ``image--.fits`` split, which sits +beside the flag split and carries the full SCAMP header (``RA---TAN`` with PV +distortion). The image PIXELS are never read: only ``fits.getheader``. + +NOT ``headers-.npy``, which is the other thing ``exp_split`` writes and +would be one small file instead of forty header reads. It is a pickled object +array of ``astropy.wcs.WCS`` INSTANCES, so reading it unpickles astropy objects +across whatever container rebuild happens next; the FITS headers are +self-describing text and cost milliseconds. ``merge_headers`` may live with the +pickle because it is the file's author's own consumer; a second consumer should +not inherit the coupling. + +HOW A PIXEL IS RASTERIZED, AND WHY IT IS SAMPLED RATHER THAN INTEGRATED. A +healpix pixel at nside 131072 is 1.61 arcsec across; a MegaCam pixel is 0.187 +arcsec. So one healpix pixel covers ~74 CCD pixels, and the question is never +"which healpix pixels does this CCD pixel cover" but "which healpix pixels does +the flagged REGION touch". The answer is taken by sampling each flagged CCD +pixel on an ``oversample`` x ``oversample`` grid spanning its full extent +(``linspace(-0.5, 0.5)``, so the four corners are always sampled), converting to +sky through that CCD's WCS and binning with ``ang2pix``. + +THE RASTERIZATION IS CONSERVATIVE, deliberately: a healpix pixel is masked if +any part of the flagged region falls in it. The alternative — mask when the +healpix pixel's CENTRE is flagged — erases a one-pixel bad column entirely, +which is precisely the geometry #878 exists to keep. The cost is that a thin +defect is widened to the healpix resolution; at 1.61 arcsec that is the price +of the ladder's nside. + +WHAT ``oversample`` BUYS, MEASURED (exposure 2079612p, CCD 0, 430886 flagged +pixels; the reference is a 12x12 interior grid, 14229 healpix pixels): + + oversample samples/pixel healpix pixels missed vs reference + 2 (corners) 4 14081 169 (1.2%) + 3 9 14198 54 (0.4%) + 5 25 14254 5 (0.04%) + +and the cost is linear in the sample count. 3 is the default in ``config.yaml`` +for that reason, and it rides on the rule's ``params`` so raising it re-rasterizes +without touching anything upstream. 2 is the geometric floor: the four corners +of a CCD pixel bound its footprint, and a healpix pixel 74 times its area cannot +sit inside it — the residual 1.2% is boundary and rounding, not a class of +missed defect. + +MEMORY IS BOUNDED BY THE BATCH, NOT BY THE EXPOSURE. Peak RSS is one CCD's +flag image plus one batch of ``CHUNK`` source pixels' worth of coordinates — +measured 0.62 GB on a real 4.1%-flagged exposure at oversample 3, and 0.74 GB on +a FULLY flagged CCD, which is the case the batching exists for (a dead or +saturated MegaCam chip is 9.4M flagged pixels and, unbatched, several GB; +``rasterize_ccd`` argues it). The bound is the batch, not the exposure. + +BYTE-STABLE, tmp-then-``cmp``-then-``mv``, the pattern ``persist_exp`` uses: +healsparse's FITS output carries no timestamp (checked), so re-rasterizing an +unchanged store produces an identical file and leaves its mtime alone. mtime is +a rerun trigger and ``clean_exposure`` waits on this rule's manifest, so an +unconditional rewrite would make every reclamation look out of date once per +invocation. + +THE MANIFEST IS THE ONLY DECLARED OUTPUT and it lives on the PERSISTENT root +beside the fragment (``/exp///manifests/``), not in +the exposure's scratch ``manifests/`` which ``clean_exposure`` deletes wholesale +— same placement, and same reason, as ``exp_persist``. It records the per-CCD +healpix counts, so a reader can see which CCD contributed what without opening +the map, and the SHA-256 of the fragment file. The digest is what makes the +manifest a faithful DAG edge: a re-rasterization that moves a defect while +keeping every count changes the fragment, and without the digest it would leave +the manifest byte-identical and untouched, so ``defect_map_merge`` would never +learn of it. The merge reads the same digest back as its record of what went +into the union. + +ORDERED BEFORE RECLAMATION. ``clean_exposure`` takes this manifest as an input, +exactly as it takes ``exp_persist``'s: the flag splits live on /scratch and go +with the store, so the fragment must be on /project before anything is deleted. +""" + +import argparse +import filecmp +import hashlib +import json +import re +import sys +import warnings +from pathlib import Path + +import numpy as np +from astropy.io import fits +from astropy.wcs import WCS + +import healpy as hp +import healsparse as hsp + +# The split stage's run dir (RUN_NAME in config_exp_Sp.ini) and its module. +# Hardcoded for the same reason persist_exp.py hardcodes its own: this rule +# rasterizes the SPLIT stage's flag images and nothing else, and a knob here +# would be a knob for "rasterize some other stage". +RUN_NAME = "run_sp_exp_Sp" +MODULE = "split_exp_runner" + +# `flag-2079612-13.fits` -> ccd 13. The number string is the exposure's +# ($SP_UNIT_NUM, `-2079612`), so the CCD is what follows the last dash. +_CCD = re.compile(r"-(\d+)\.fits$") + + +def split_dir(exp_dir: Path) -> Path: + return exp_dir / "output" / RUN_NAME / MODULE / "output" + + +def ccd_files(exp_dir: Path, n_ccds: int) -> list: + """``(ccd, flag path, image path)`` for every CCD this exposure split, in + CCD order — all ``n_ccds`` of them or none at all. + + COUNTED AGAINST THE EXPECTED CCD COUNT, not against whatever is on disk, and + that is the whole guard. ``split_exp`` writes image, weight and flag for + each of ``N_HDU`` CCDs in one pass, so the reachable failure is not "a flag + without its image" — it is an incompletely MATERIALISED split dir: an + age-based scratch purge deleting files one at a time, a store copied or + restored half way, a truncated rsync. Globbing for ``flag-*.fits`` and + rasterizing whatever comes back turns that into a fragment covering half the + exposure, written with ``"status": "complete"`` and with nothing downstream + able to notice — a hole in the footprint, which is precisely what this rule + exists to prevent. So the expected count comes in on ``params`` (the + Snakefile reads ``N_HDU`` from ``config_exp_Sp.ini``) and a short split is a + hard error. + + The image split is looked up beside each flag for its header alone, and a + missing one is the same hard error for the same reason. + """ + root = split_dir(exp_dir) + out = [] + for flag in sorted(root.glob("flag-*.fits")): + match = _CCD.search(flag.name) + if not match: + continue + image = flag.with_name(flag.name.replace("flag-", "image-", 1)) + if not image.exists(): + sys.exit(f"defect_map_exp: {flag.name} has no {image.name} beside " + f"it in {root}; the WCS lives on the image split (see the " + f"module docstring)") + out.append((int(match.group(1)), flag, image)) + if len(out) != n_ccds: + sys.exit(f"defect_map_exp: {root} holds {len(out)} flag split(s), not " + f"the {n_ccds} this exposure was split into; the split dir is " + f"incomplete and a fragment built from it would be a hole in " + f"the footprint marked complete (see ccd_files' docstring). " + f"Re-run exp_split for this exposure: delete its scratch " + f"manifests/exp_split.json and the workflow rebuilds the " + f"split. Until it does, clean_exposure waits on this rule " + f"and the store stays — a damaged store is not reclaimed " + f"silently.") + return sorted(out) + + +def offsets(oversample: int) -> tuple: + """Sample offsets within one CCD pixel, in pixel units. + + ``linspace`` with both endpoints, so the CORNERS are always sampled: they + are what bounds the pixel's footprint, and the interior samples only fill + boundary and rounding gaps (the docstring's table measures how many). + """ + if oversample < 2: + sys.exit(f"defect_map_exp: oversample={oversample} would sample the " + f"pixel centre alone and lose the pixel's extent; 2 is the " + f"geometric floor (its four corners)") + step = np.linspace(-0.5, 0.5, oversample) + grid_x, grid_y = np.meshgrid(step, step) + return grid_x.ravel(), grid_y.ravel() + + +# Flagged CCD pixels converted per batch. Peak RSS is set by THIS, not by how +# bad the CCD is: one batch at oversample 3 is 500k x 9 samples x two float64 +# coordinate arrays in and two out, plus astropy's PV/SIP temporaries, and the +# accumulated result is a deduplicated int64 pixel list bounded by the CCD's +# healpix footprint (124601 ids for a WHOLE CCD at nside 131072), not by the +# sample count. Measured on a fully flagged MegaCam chip (2048 x 4612 = 9.4M +# pixels, the worst case there is) against 2079612p CCD 1's real WCS: 0.74 GB +# peak and 18.3 s at 500k, 1.13 GB and 18.9 s at 1M. Time is flat in the batch +# size and memory is linear in it, so the smaller batch is free. +CHUNK = 500_000 + + +def rasterize_ccd(flag_path: Path, image_path: Path, nside: int, + off_x, off_y) -> np.ndarray: + """The healpix pixel ids (NEST, ``nside``) this CCD's flags touch. + + @sc [decision:masking.defect_map_from_flags,label:convention] defect-fragment-contains-flags + Conservative by construction: every CCD pixel with a nonzero flag, centre + and four corners, lands in a returned healpix pixel. The corners are what + bound a pixel's footprint, so ``offsets`` always samples them; sampling + centres alone erases one-pixel bad columns and lone hot pixels that + straddle a healpix boundary. Pixel coordinates are 1-based into + ``all_pix2world(..., 1)``. Checked by + ``tests/unit/test_defect_map_contains_flags.py``. + + ONE CCD AT A TIME AND, WITHIN IT, ONE BATCH AT A TIME. The first is why the + loop in ``main`` is a loop; the second is why this one is. A MegaCam + exposure routinely carries a dead or saturated chip, and a FULLY flagged CCD + is 2048 x 4612 = 9.4M nonzero pixels — 85M samples at oversample 3. Held in + one shot that is ~680 MB per coordinate array in and the same again out of + ``all_pix2world``, plus astropy's own PV/SIP temporaries: several GB, well + over the rule's request, on all three attempts. Batched it is 0.74 GB + (measured), inside the request with room to spare. The exposure would then + never get a fragment AND, because ``clean_exposure`` waits on this rule's + manifest, never be reclaimed either. The 0.62 GB measured on a 4.4%-flagged + exposure says nothing about that case; ``CHUNK`` does. + + The batches are ``np.unique``-reduced as they go, so what survives across + them is the CCD's healpix footprint and not its samples. + """ + with warnings.catch_warnings(): + # SCAMP headers carry a deprecated RADECSYS and a redundant SIP block + # beside the PV distortion astropy actually uses; both are FITSFixedWarning + # noise on every one of 40 CCDs and neither changes the transform. + warnings.simplefilter("ignore") + wcs = WCS(fits.getheader(image_path)) + data = fits.getdata(flag_path) + rows, cols = np.nonzero(data) + del data + if rows.size == 0: + return np.empty(0, dtype=np.int64) + found = np.empty(0, dtype=np.int64) + for start in range(0, rows.size, CHUNK): + stop = start + CHUNK + # 1-based FITS pixel coordinates, sampled across each flagged pixel's + # extent. + x = (cols[start:stop, None] + 1.0 + off_x[None, :]).ravel() + y = (rows[start:stop, None] + 1.0 + off_y[None, :]).ravel() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + ra, dec = wcs.all_pix2world(x, y, 1) + del x, y + pixels = hp.ang2pix(nside, ra, dec, lonlat=True, nest=True) + del ra, dec + found = np.union1d(found, pixels) + del pixels + return found + + +def file_digest(path: Path) -> str: + """SHA-256 of a file's bytes, hex: the fragment's identity in its manifest + and in the merge's sidecar.""" + with open(path, "rb") as fh: + return hashlib.file_digest(fh, "sha256").hexdigest() + + +def write_stable(tmp: Path, dest: Path) -> None: + """Move ``tmp`` onto ``dest``, or drop it when the bytes already match.""" + if dest.exists() and filecmp.cmp(tmp, dest, shallow=False): + tmp.unlink() + else: + tmp.replace(dest) + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--exp-dir", required=True, type=Path, + help="the exposure's scratch store") + parser.add_argument("--exp", required=True) + parser.add_argument("--dest", required=True, type=Path, + help="/exp///defect; the " + "fragment is /defect-.hsp") + parser.add_argument("--manifest", required=True, type=Path) + parser.add_argument("--n-ccds", required=True, type=int, + help="how many CCDs exp_split wrote (N_HDU in " + "config_exp_Sp.ini); a split dir short of this is " + "an error, not a smaller fragment") + parser.add_argument("--nside", type=int, default=131072, + help="nside_sparse; the mask ladder's resolution") + parser.add_argument("--nside-coverage", type=int, default=128) + parser.add_argument("--oversample", type=int, default=3, + help="samples per CCD pixel per axis (see the module " + "docstring's measured table)") + args = parser.parse_args() + + ccds = ccd_files(args.exp_dir, args.n_ccds) + if not ccds: + sys.exit(f"defect_map_exp: {args.exp}: no flag split under " + f"{split_dir(args.exp_dir)}") + + off_x, off_y = offsets(args.oversample) + fragment = hsp.HealSparseMap.make_empty( + args.nside_coverage, args.nside, np.bool_, bit_packed=True) + per_ccd = {} + for ccd, flag_path, image_path in ccds: + pixels = rasterize_ccd(flag_path, image_path, args.nside, off_x, off_y) + per_ccd[str(ccd)] = int(pixels.size) + if pixels.size: + fragment[pixels] = True + + args.dest.mkdir(parents=True, exist_ok=True) + frag_path = args.dest / f"defect-{args.exp}.hsp" + tmp = frag_path.with_name(frag_path.name + ".tmp") + try: + fragment.write(str(tmp), clobber=True) + write_stable(tmp, frag_path) + finally: + tmp.unlink(missing_ok=True) + + body = { + "stage": "exp_defect_map", "level": "exp", "unit": args.exp, + "status": "complete", + "map": str(frag_path), + "nside": args.nside, + "nside_coverage": args.nside_coverage, + "oversample": args.oversample, + "n_ccds": len(ccds), + # Per CCD, so a reader can see WHICH CCD contributed what without + # opening the map — a CCD at zero is a real thing (a clean chip) and a + # whole exposure at zero is not. + "pixels_per_ccd": per_ccd, + "n_pixels": int(fragment.n_valid), + "n_coverage_pixels": int(fragment.coverage_mask.sum()), + "bytes": frag_path.stat().st_size, + # The fragment's CONTENT, so the manifest changes whenever the fragment + # does (see the module docstring). + "sha256": file_digest(frag_path), + } + args.manifest.parent.mkdir(parents=True, exist_ok=True) + tmp = args.manifest.with_name(args.manifest.name + ".tmp") + try: + tmp.write_text(json.dumps(body, indent=2, sort_keys=True) + "\n") + write_stable(tmp, args.manifest) + finally: + tmp.unlink(missing_ok=True) + + print(f"[defect_map_exp] {args.exp}: {len(ccds)} CCD(s), " + f"{body['n_pixels']} healpix pixel(s) over " + f"{body['n_coverage_pixels']} coverage pixel(s), " + f"{body['bytes'] / 1e6:.1f} MB -> {frag_path}") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/exp_footprint.py b/workflow/scripts/exp_footprint.py new file mode 100644 index 000000000..7b3b5257f --- /dev/null +++ b/workflow/scripts/exp_footprint.py @@ -0,0 +1,211 @@ +#!/usr/bin/env python3 +"""Record ONE exposure's per-CCD sky footprint, for the CCDs that have a PSF model. + +Run as the shell of the in-DAG ``exp_footprint`` rule, never by hand. + +The footprint records the four sky corners of each CCD with a valid PSF model. +It joins the persisted PSF member list to the split stage's local WCS array and +writes ``/exp///manifests/exp_footprint.json``. +``nexp_map`` reads these durable records to count exposures per sky pixel; +``clean_exposure`` waits for the record before deleting its scratch inputs. + +THE TWO INPUTS, AND WHY EACH IS THE RIGHT ONE. + +1. WHICH CCDS HAVE A PSF -> ``exp_persist.json``'s ``files[].name``, the + ``validation_psf--.fits`` members of the persisted tar. This is + EXACT, not inferred: ``psfex_interp`` returns WITHOUT writing that file on + NOT_ENOUGH_STARS, BAD_CHI2 and FILE_NOT_FOUND, and writes it on success, so + the member list IS the valid-PSF set. It is also DURABLE (persistent root, + survives ``clean_exposure``) and a DECLARED RULE OUTPUT, so the DAG orders + this rule after it for free. ``exp_psf.json``'s psfex_interp entry is a + COUNT, not a list of names, and lives inside the directory + ``clean_exposure`` deletes wholesale. + +2. WCS -> ``headers-.npy``, written by ``split_exp`` beside the per-CCD + images. A length-N ``dtype=object`` array; element ``i`` is + ``{"WCS": WCS(h), "header": h.tostring()}``, loaded with ``allow_pickle=True`` + exactly as ``merge_headers.py`` does. + +THE LOAD-BEARING INVARIANT, PINNED BY tests/unit/test_exp_footprint.py: +array index ``i`` == the CCD index in ``image--.fits`` == ``-`` +== ``validation_psf--.fits``. ``split_exp`` writes both the image and the +array element from the same ``idx-1`` in one loop, so the alignment is +structural. The tests enforce this cross-component contract between +split_exp's numbering, psfex_interp's filenames and this record's ids. +The footprint assigns each CCD's id from its position in the array; the +persisted filenames select which of those indices have a PSF model. + +WHY THE SHAPE COMES FROM THE STORED HEADER AND NOT FROM THE WCS. astropy hands +back the DECOMPRESSED header for a tile-compressed HDU, so this npy carries true +``NAXIS1/2`` and no ``ZIMAGE``. ``_image_shape`` reads the dimensions from the +header, including ``ZNAXIS1/2`` when given a compressed binary-table header. +``WCS`` drops the ``Z*`` keywords, so the stored header is the shape source. + +DATA AND MANIFEST IN ONE FILE. 40 rows of 8 floats is a few KB, and inodes are +what bind on /project (persist_exp.py argues the quota arithmetic). Per-exposure +JSON keeps parallel exposure jobs independent; appending to one shared text +file would make them contend on the same record. +``ccds_no_psf`` is carried explicitly so the record is self-describing about +attrition — a reader can tell "this CCD is not in the map" from "this CCD was +never looked at". + +Written byte-stable, no timestamp, tmp-then-``cmp``-then-``mv``, exactly as +``persist_exp.py`` does and for the same reason: mtime is a rerun trigger, and an +unconditional rewrite would make ``clean_exposure`` look out of date once per +invocation. +""" + +import argparse +import filecmp +import json +import sys +from fnmatch import fnmatch +from pathlib import Path + +import numpy as np +from astropy.io.fits import Header + +from shapepipe.utilities.ccd_footprint import ( + _ccd_corners, + _image_shape, +) + +# Same directory; the rule invokes this file by path, so it is sys.path[0]. +import persist_exp + +# The split stage's run dir (RUN_NAME in config_exp_Sp.ini) and its module +# output dir. Hardcoded rather than passed, for persist_exp.py's reason: this +# rule reads the split stage's headers and nothing else, and a knob here would +# be a knob for "read some other stage". +HEADERS_DIR = "output/run_sp_exp_Sp/split_exp_runner/output" + +# The members this rule selects, named and resolved through the catalogue +# exp_persist packs by, so the glob has one definition. They are always there to +# find: persist_exp packs this product for every exposure whatever `persist_exp:` +# says (its ALWAYS), and fails the pack rather than writing a manifest without it. +MEMBER_PRODUCT = persist_exp.ALWAYS +MEMBER_PATTERN = persist_exp.resolve(MEMBER_PRODUCT) + +# Only the PREFIX of a member name is ours to know; the rest is the unit id and +# the CCD index, which is what makes the set a set of indices. Derived from the +# pattern so it cannot drift from it. +PSF_PREFIX = MEMBER_PATTERN.split("-", 1)[0] + + +def valid_ccds(manifest, exp): + """The CCD indices with a PSF model, read off an ``exp_persist`` manifest. + + BY PRODUCT NAME, OR FAILING THAT BY FILE NAME — merge_star_cat.is_member()'s + test, and for its reasons: persist_exp labels every member with the product + it came from, but a tar packed before that field existed carries no label, + and a keep list written as a raw glob labels its members with the glob. + + No precondition on the keep list. exp_persist packs this product for every + exposure whatever `persist_exp:` says, so an empty answer here means the + exposure genuinely lost every CCD, not that the config unpacked them. + """ + body = json.loads(Path(manifest).read_text()) + ccds = set() + for entry in body.get("files") or []: + name = entry.get("name", "") + if not (entry.get("product") == MEMBER_PRODUCT + or fnmatch(name, MEMBER_PATTERN)): + continue + if not name.startswith(f"{PSF_PREFIX}-{exp}-"): + continue + stem = name[len(f"{PSF_PREFIX}-{exp}-"):].removesuffix(".fits") + if not stem.isdigit(): + sys.exit(f"exp_footprint: {exp}: cannot read a CCD index out of " + f"tar member {name!r}") + ccds.add(int(stem)) + return ccds + + +def footprint(headers, exp, with_psf): + """One record's ``ccds`` and ``ccds_no_psf``, in array order. + + ``headers`` is the ``headers-.npy`` array; index ``i`` IS the CCD + index (see the module docstring). Corners are computed only for the CCDs in + ``with_psf`` — a CCD with no PSF model contributes nothing to a map that + counts exposures with a valid PSF, and computing its corners anyway would + invite a later reader to use them. + """ + ccds, no_psf = [], [] + for i, entry in enumerate(headers): + ccd_id = f"{exp}-{i}" + if i not in with_psf: + no_psf.append(ccd_id) + continue + # The WCS is the pickled object split_exp built; the SHAPE must come + # from the stored header text (module docstring). + shape = _image_shape(Header.fromstring(entry["header"])) + ra, dec = _ccd_corners(entry["WCS"], shape) + # float(), because _ccd_corners hands back numpy scalars and this record + # has to be byte-stable: a plain double's repr is, a numpy type's + # serialisation is json's business rather than ours. + ccds.append({"id": ccd_id, + "ra": [float(x) for x in ra], + "dec": [float(x) for x in dec]}) + return ccds, no_psf + + +def write_stable(path, body): + """Write JSON tmp-then-``cmp``-then-``mv``; an unchanged body keeps its mtime.""" + path.parent.mkdir(parents=True, exist_ok=True) + tmp = path.with_name(path.name + ".tmp") + try: + tmp.write_text(json.dumps(body, indent=2, sort_keys=True) + "\n") + if path.exists() and filecmp.cmp(tmp, path, shallow=False): + tmp.unlink() # unchanged: leave the mtime alone + else: + tmp.replace(path) + finally: + tmp.unlink(missing_ok=True) + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--exp-dir", required=True, type=Path, + help="the exposure's scratch store") + p.add_argument("--exp", required=True) + p.add_argument("--persist-manifest", required=True, type=Path, + help="/exp///manifests/" + "exp_persist.json; names the valid-PSF CCDs") + p.add_argument("--manifest", required=True, type=Path) + args = p.parse_args() + + npy = args.exp_dir / HEADERS_DIR / f"headers-{args.exp}.npy" + if not npy.exists(): + sys.exit(f"exp_footprint: {args.exp}: no WCS array at {npy} — the " + f"split stage's store is gone or never wrote one " + f"(config_exp_Sp.ini's OUTPUT_SUFFIX must include `image`).") + headers = np.load(npy, allow_pickle=True) + + with_psf = valid_ccds(args.persist_manifest, args.exp) + # A PSF file for a CCD the split never produced means the two stages + # disagree about the focal plane — the one failure the id alignment cannot + # absorb, so it is loud rather than silently dropped. + stray = sorted(i for i in with_psf if i >= len(headers)) + if stray: + sys.exit(f"exp_footprint: {args.exp}: {args.persist_manifest} names " + f"CCD(s) {stray} but {npy} holds only {len(headers)}") + + ccds, no_psf = footprint(headers, args.exp, with_psf) + + write_stable(args.manifest, { + "stage": "exp_footprint", "level": "exp", "unit": args.exp, + "status": "complete", + "source_manifest": str(args.persist_manifest), + "n_ccd_headers": len(headers), + "n_valid_psf": len(ccds), + "ccds": ccds, + "ccds_no_psf": no_psf, + }) + + warn = f" ({len(no_psf)} without a PSF model)" if no_psf else "" + print(f"[exp_footprint] {args.exp}: {len(ccds)}/{len(headers)} CCD " + f"footprint(s){warn} -> {args.manifest}") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/merge_defect_map.py b/workflow/scripts/merge_defect_map.py new file mode 100644 index 000000000..973796402 --- /dev/null +++ b/workflow/scripts/merge_defect_map.py @@ -0,0 +1,472 @@ +#!/usr/bin/env python3 +"""Union the campaign's per-exposure defect fragments into ONE healsparse map. + +Run as the shell of the campaign-level ``defect_map_merge`` rule, never by hand. + +WHAT IT PRODUCES, AND FOR WHOM. +``/defect_map/defect_map_.hsp``: +a boolean ``HealSparseMap``, ``True`` where any exposure of the campaign flagged +the sky, at the mask ladder's own resolution (nside_sparse 131072 over +nside_coverage 128, ``bit_packed``). It is the pixel-domain half of the masking +the footprint could not otherwise see (CosmoStat/shapepipe#878; +``defect_map_exp.py`` argues why the map has to come from here at all). + +IT IS A CAMPAIGN PRODUCT, NOT AN INPUT. Nothing in this workflow consumes it: +``config_tile_Mc.ini``'s ``MASK_EXT_PATHS`` names maps that exist before the run +starts, and this one exists only after it. Adding it to that ladder is a +deliberate, later config edit against a path a campaign has actually produced — +that file's header carries the recipe, and deliberately not the entry. + +TRUE MEANS MASKED, matching every other map in the ladder. A position outside +the campaign's coverage reads ``False``, indistinguishable from clean, exactly +as the UNIONS bit maps behave outside theirs; the union's coverage is the +campaign's exposures and nothing else says where that is. + +IT RECONCILES, AND THE ASYMMETRY IS THE WHOLE DESIGN. The output must be a +function of the input set — that is what makes the rule's fingerprint mean +anything — but a union is not invertible: + + * an exposure whose fragment is NEW is OR-ed into the map on the spot. This is + the common case (a campaign grows by appending tiles) and it reads exactly + the appended exposures; + * an exposure that LEFT the campaign, or whose fragment CHANGED, forces a + REBUILD from every fragment, because nothing can un-OR a pixel that two + exposures both set. Rebuilding is honest about that rather than leaving a + stale bit nothing would ever notice; + * a plan with neither leaves the file UNTOUCHED — not rewritten identically, + untouched, so its mtime cannot move. mtime is a rerun trigger. + +WHAT IS AND IS NOT A FUNCTION OF THE INPUT SET, in the same terms +``hdf5_reconcile.py`` sets them. The map's CONTENT is: the same fragments give +the same valid pixels, the same counts and the same sidecar, whether they +arrived at once or one append at a time. Its BYTES are not — reaching a state by +append rather than by rebuild round-trips the map through healsparse's reader +and writer, which can lay the same pixels out in a different number of 2880-byte +FITS blocks (measured: 1,586,880 B rebuilt vs 1,589,760 B appended, identical +``valid_pixels``). That is the trade for not re-reading the campaign. It means +``write_stable`` below can move the map's mtime on a rebuild that changed +nothing — cheap while nothing consumes the map, and the thing to fix (by +rebuilding whenever the append path would rewrite anyway) if something ever +does. The no-op case is unaffected: it compares the PLAN, not the bytes, and +never reads the map's pixels, only its coverage table. + +Which exposures are already in the map is recorded in a SIDECAR beside it +(``defect_map_.json``), each with its fragment's SHA-256 as its +``exp_defect_map`` manifest records it. The CONTENT digest, not a size/mtime +stamp: a re-rasterization can move a defect while keeping the file's size, and +the manifest carries the digest precisely so that such a change reaches this +rule as a changed input. Reading it from the manifest costs one small JSON read +per exposure, where hashing every fragment would read the whole campaign on +every invocation. The map itself cannot carry that +record: a healsparse FITS header is no place for twenty thousand exposures. The +sidecar is therefore a DECLARED OUTPUT of the rule alongside the map; losing one +without the other would be a map nobody can reconcile, and snakemake removing +both on a failure is the correct recovery (the next run rebuilds). + +MAP AND SIDECAR ARE BOUND BY A GENERATION ID, because two files cannot be +replaced atomically together. Both are written to temporaries first, so a kill +before publication leaves the previous pair; a kill BETWEEN the two renames +still leaves a new map under the old sidecar, and snakemake cannot be relied on +to clean up after its own head process dies. So each merge stamps +``generation`` — a digest of the resolution and every exposure's fragment +digest, the whole of what determines the map's content — into both the map's +healsparse metadata (``DMAPGEN`` in its primary FITS header) and the sidecar. +``reconcile_plan`` reads the header and rebuilds on any disagreement. The id is +a function of the input set, not a fresh token, so an unchanged campaign +rewrites nothing. + +MEMORY IS FLAT IN THE NUMBER OF EXPOSURES, which is the reason for the loop +below and not for a list comprehension over ``HealSparseMap.read``. Fragments +are never held together and are never OR-ed as maps: each is read, reduced to +its ``valid_pixels`` (~360k int64, ~3 MB on a real exposure), set into the +accumulator, and dropped. The job's footprint is therefore ONE accumulator plus +ONE fragment, whether the campaign is 127 exposures or 20k. The accumulator is +a function of the campaign's FOOTPRINT, not of its exposure count: a coverage +pixel costs (nside/nside_coverage)^2 / 8 bytes bit-packed — 128 KiB at the +ladder's resolution — so a DR6-scale footprint (~23k coverage pixels, measured +on the UNIONS ugriz maps) is ~3 GB resident and the rule is sized on exactly +that count. + +WHICH EXPOSURES — AND WHY THE JOB DERIVES THE SET. The campaign's: every +exposure of every tile both declared in ``tile_list`` and present in the index, +through ``build_index.campaign_exposures`` so there is one definition and not +two that can drift. It is derived rather than passed because at DR6 scale the +set is ~20k paths and a shell command reaches ``execve`` as a single argv entry +capped at 128 KiB; the rule's ``input`` is the DAG edge and its ``params`` +carries a fingerprint of the same ids. + +A CAMPAIGN EXPOSURE WITH NO FRAGMENT IS SKIPPED, NOT AN ERROR, and that is the +one place this differs from ``merge_final_cat``. Fragments accumulate on the +persistent root across campaigns and survive reclamation, but an exposure +reclaimed by a workflow PREDATING this rule has none and never will without a +rebuild from VOS. Failing would make the map unbuildable for exactly the +campaigns that most want it; the count of skipped exposures is reported and +recorded in the sidecar instead. +""" + +import argparse +import filecmp +import hashlib +import json +import sys +from pathlib import Path + +import numpy as np + +import healsparse as hsp +from astropy.io import fits + +# Same directory; the rule invokes this file by path, so it is sys.path[0]. +# The vocabulary below (`Plan`, `empty`, `describe`, plan-then-apply, +# untouched-on-a-no-op) is hdf5_reconcile.py's, deliberately. +import build_index +# The fragment's writer defines its digest; one definition, so the digest a +# manifest records and the one computed here for a manifest without it agree. +from defect_map_exp import file_digest + + +def fragment_path(products_dir: Path, exp: str) -> Path: + """Where ``exp_defect_map`` wrote this exposure's fragment.""" + return (products_dir / "exp" / exp[:2] / exp / "defect" + / f"defect-{exp}.hsp") + + +def manifest_path(products_dir: Path, exp: str) -> Path: + """``exp_defect_map``'s declared output: the Snakefile's + ``prod_exp_manifest(exp, "exp_defect_map")``.""" + return (products_dir / "exp" / exp[:2] / exp / "manifests" + / "exp_defect_map.json") + + +def fragment_digest(products_dir: Path, exp: str) -> str: + """The fragment's SHA-256, as its manifest records it. + + A manifest that lacks the digest, or cannot be read, is answered by hashing + the fragment itself: an exposure whose store is reclaimed cannot be + re-rasterized to supply one. + """ + try: + return json.loads(manifest_path(products_dir, exp).read_text())["sha256"] + except (OSError, ValueError, KeyError): + return file_digest(fragment_path(products_dir, exp)) + + +def fragments(products_dir: Path, tile_list: Path, index_db: Path) -> tuple: + """``({exp: fragment path}, [exposures with no fragment])``, in exposure order.""" + have, missing = {}, [] + for exp in build_index.campaign_exposures(tile_list, index_db): + path = fragment_path(products_dir, exp) + if path.exists(): + have[exp] = path + else: + missing.append(exp) + return have, missing + + +class Plan: + """What reconciling this campaign into this map requires. + + ``append`` is the cheap path — OR these fragments into the map on disk. + ``rebuild`` is the honest one: a union cannot drop a pixel, so a removal or + a changed fragment means reading every fragment again. + """ + + def __init__(self, append, rebuild, reason): + self.append, self.rebuild, self.reason = append, rebuild, reason + + def empty(self): + return not (self.append or self.rebuild) + + def describe(self): + if self.rebuild: + return f"rebuilt from {len(self.rebuild)} fragment(s) ({self.reason})" + return f"{len(self.append)} fragment(s) appended" + + +# The map's metadata key for the generation id; a FITS keyword, so <= 8 chars. +GENERATION_KEY = "DMAPGEN" + + +def generation(digests: dict, nside: int, nside_coverage: int) -> str: + """The id binding a map to its sidecar: a digest of everything that + determines the map's content, so the same inputs give the same id.""" + blob = json.dumps({"nside": nside, "nside_coverage": nside_coverage, + "exposures": digests}, sort_keys=True) + return hashlib.sha256(blob.encode()).hexdigest() + + +def map_generation(output: Path): + """The generation id stamped in the map's primary header, or ``None``.""" + try: + return fits.getheader(str(output), 0).get(GENERATION_KEY) + except OSError: + return None + + +def read_sidecar(path: Path) -> dict: + try: + return json.loads(path.read_text()) + except (OSError, ValueError): + return {} + + +def reconcile_plan(output: Path, sidecar: Path, digests: dict, nside: int, + nside_coverage: int) -> Plan: + """Compare what is on disk with the campaign, WITHOUT writing anything. + + @sc [decision:masking.defect_map_from_flags,label:convention] defect-map-is-the-sidecars-union + The map is the OR of exactly the fragments its sidecar records, and never + un-OR-ed: an exposure that left the campaign, or a fragment that changed + on disk, is a rebuild from every fragment, because a union cannot tell + which exposure set a pixel. Only a pure addition may append in place. + Checked by ``tests/unit/test_defect_map_reconcile.py``. + + ``digests`` is ``{exp: fragment_digest}`` over the campaign's fragments; a + fragment has changed exactly when its digest differs from the recorded one. + + A missing map, or a sidecar that does not describe it, is a rebuild: the two + are written together and either one alone is not evidence about the other. + "Describes it" is checked, not assumed: the map's ``DMAPGEN`` header must + equal the sidecar's ``generation`` (the module docstring argues why). + + So is a map at another resolution than the one requested. Unchanged + fragments would otherwise make that an empty plan, leaving the map at its + old resolution under a sidecar rewritten to claim the new one. The rebuild + then reads the fragments, and ``accumulate`` refuses any at the wrong + resolution before either file is written. The resolution is read from the + map itself — its coverage table, a few kilobytes — not from the sidecar. + """ + record = read_sidecar(sidecar) + known = record.get("exposures") or {} + if not output.exists() or not known: + return Plan([], sorted(digests), "no map on disk") + + if map_generation(output) != record.get("generation"): + return Plan([], sorted(digests), + "map and sidecar are from different merges") + + cov = hsp.HealSparseCoverage.read(str(output)) + if (cov.nside_sparse, cov.nside_coverage) != (nside, nside_coverage): + return Plan([], sorted(digests), + f"map is nside {cov.nside_sparse} / nside_coverage " + f"{cov.nside_coverage}, not the requested {nside} / " + f"{nside_coverage}") + + gone = sorted(set(known) - set(digests)) + if gone: + return Plan([], sorted(digests), + f"{len(gone)} exposure(s) left the campaign") + changed = sorted(exp for exp, digest in digests.items() + if exp in known and known[exp] != digest) + if changed: + return Plan([], sorted(digests), + f"{len(changed)} fragment(s) changed on disk") + return Plan(sorted(set(digests) - set(known)), [], "") + + +def union_coverage(paths, nside_coverage) -> np.ndarray: + """The coverage pixels every fragment in ``paths`` touches, together. + + A PRE-PASS, so the accumulator is allocated ONCE. ``target[pixels] = True`` + into a coverage pixel the map has not seen yet makes healsparse GROW its + sparse array, which copies it; at DR6 scale that array is gigabytes and a + rebuild discovers coverage pixels all the way through the campaign, so the + copies dominate everything the "reading fragments dominates" comment below + models. Seeding the coverage up front turns O(fragments) reallocations of a + growing array into one allocation of the final one. + + It costs a second read of each fragment's COVERAGE TABLE only — + ``HealSparseCoverage.read`` never touches the sparse array — which is + kilobytes against the megabytes the accumulation itself reads. + + MEASURED, on synthetic fragments at the campaign's own resolution (nside + 131072 / coverage 128), 400 fragments discovering 5131 coverage pixels — a + 656 MB accumulator: accumulation 8.4 s unseeded, 7.1 s seeded, with a 0.7 s + coverage pre-pass. So the reallocations are ~16% of the accumulation here, + not the dominant term a naive "copy the array once per new coverage pixel" + reading predicts (healsparse grows the sparse array in blocks). The win + grows with the accumulator; the pre-pass does not. Both paths produced + identical maps. + + Only the rebuild path uses it: an append starts from the map on disk, whose + coverage is already most of the footprint, and reads a handful of fragments. + """ + mask = None + for path in paths: + cov = hsp.HealSparseCoverage.read(str(path)) + if cov.nside_coverage != nside_coverage: + # accumulate() is the one place that reports a fragment built at the + # wrong resolution, with the exposure id and what to do about it. + # Here it is only a seed: give up on it and let that error stand. + return None + mask = (cov.coverage_mask.copy() if mask is None + else mask | cov.coverage_mask) + return None if mask is None else np.where(mask)[0] + + +def accumulate(target, paths, nside_coverage, nside) -> None: + """OR each fragment into ``target``, ONE AT A TIME (see the docstring). + + A fragment at the wrong resolution is a hard error rather than a silent + upgrade: the ladder's nside is a campaign-wide decision, and a fragment that + disagrees with it was rasterized by a differently-configured run. + """ + for exp, path in paths: + frag = hsp.HealSparseMap.read(str(path)) + if (frag.nside_sparse, frag.nside_coverage) != (nside, nside_coverage): + sys.exit(f"merge_defect_map: {path} is nside_sparse " + f"{frag.nside_sparse} / nside_coverage " + f"{frag.nside_coverage}, not {nside} / {nside_coverage}; " + f"re-rasterize {exp} before merging") + pixels = frag.valid_pixels + del frag + if pixels.size: + target[pixels] = True + del pixels + + +def write_stable(tmp: Path, dest: Path) -> None: + if dest.exists() and filecmp.cmp(tmp, dest, shallow=False): + tmp.unlink() + else: + tmp.replace(dest) + + +def build_record(digests: dict, missing: list, nside_coverage: int, nside: int, + n_pixels: int, n_coverage: int) -> dict: + """The sidecar: what is in the map, and what the campaign wanted but lacked. + + Built apart from writing the map because it is not only the map's record. + Two of its fields describe the CAMPAIGN — how many exposures it has and + which of them have no fragment — and those can move while the map itself + cannot: add tiles whose exposures were all reclaimed by a workflow + predating this rule and there is nothing to append, nothing to rebuild, and + a sidecar still reporting the previous campaign's counts. The docstring says + a short map should say so ON DISK; that means the record has to be rewritten + even when the map is untouched. + """ + return { + "campaign_exposures": len(digests) + len(missing), + "nside": nside, + "nside_coverage": nside_coverage, + "n_pixels": n_pixels, + "n_coverage_pixels": n_coverage, + # Must equal the DMAPGEN header of the map beside it. + "generation": generation(digests, nside, nside_coverage), + # What the NEXT invocation reconciles against; sorted so the sidecar is + # byte-stable for a given campaign state. + "exposures": dict(sorted(digests.items())), + # Recorded rather than merely printed: a map short of exposures should + # say so on disk, not only in a job log nobody keeps. + "exposures_without_fragment": sorted(missing), + } + + +def dump_record(record: dict) -> str: + return json.dumps(record, indent=2, sort_keys=True) + "\n" + + +def write_sidecar(sidecar: Path, record: dict) -> None: + tmp = sidecar.with_name(sidecar.name + ".tmp") + try: + tmp.write_text(dump_record(record)) + write_stable(tmp, sidecar) + finally: + tmp.unlink(missing_ok=True) + + +def apply_plan(output: Path, sidecar: Path, plan: Plan, have: dict, + digests: dict, missing: list, nside_coverage: int, + nside: int) -> dict: + """Carry the plan out on tmp copies, then move both files into place. + + Both temporaries are complete before either rename, so a crash before + publication leaves the previous PAIR intact. A crash between the two + renames leaves a new map under the old sidecar; their generation ids then + disagree and the next ``reconcile_plan`` rebuilds. + """ + if plan.rebuild: + todo = plan.rebuild + target = hsp.HealSparseMap.make_empty( + nside_coverage, nside, np.bool_, bit_packed=True, + cov_pixels=union_coverage( + [have[exp] for exp in todo], nside_coverage)) + else: + target = hsp.HealSparseMap.read(str(output)) + todo = plan.append + accumulate(target, [(exp, have[exp]) for exp in todo], + nside_coverage, nside) + + record = build_record(digests, missing, nside_coverage, nside, + int(target.n_valid), + int(target.coverage_mask.sum())) + + target.metadata = {GENERATION_KEY: record["generation"]} + map_tmp = output.with_name(output.name + ".tmp") + sidecar_tmp = sidecar.with_name(sidecar.name + ".tmp") + try: + target.write(str(map_tmp), clobber=True) + sidecar_tmp.write_text(dump_record(record)) + write_stable(map_tmp, output) + write_stable(sidecar_tmp, sidecar) + finally: + map_tmp.unlink(missing_ok=True) + sidecar_tmp.unlink(missing_ok=True) + return record + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--products-dir", required=True, type=Path, + help="the persistent root; fragments are found beneath it") + parser.add_argument("--tile-list", required=True, type=Path, + help="the campaign's tile list (config tile_list)") + parser.add_argument("--index-db", required=True, type=Path, + help="the campaign's run index (config outputs.index_db)") + parser.add_argument("--output", required=True, type=Path) + parser.add_argument("--sidecar", required=True, type=Path, + help="the reconciliation record written beside the map") + parser.add_argument("--nside", type=int, default=131072) + parser.add_argument("--nside-coverage", type=int, default=128) + args = parser.parse_args() + + have, missing = fragments(args.products_dir, args.tile_list, args.index_db) + if not have: + # An empty map would satisfy every downstream existence check and mask + # nothing anywhere. + sys.exit(f"merge_defect_map: no campaign exposure in {args.tile_list} " + f"has a defect fragment under {args.products_dir}") + + digests = {exp: fragment_digest(args.products_dir, exp) for exp in have} + args.output.parent.mkdir(parents=True, exist_ok=True) + plan = reconcile_plan(args.output, args.sidecar, digests, args.nside, + args.nside_coverage) + if plan.empty(): + # The MAP is untouched — that is what an empty plan means, and its mtime + # must not move. The RECORD still can be stale: campaign_exposures and + # exposures_without_fragment describe the campaign, not the map, so + # tiles whose exposures all lack fragments change them without changing + # a single bit of the union. Rewrite it alone when it differs; + # write_stable drops the tmp when it does not. + old_record = read_sidecar(args.sidecar) + record = build_record(digests, missing, args.nside_coverage, args.nside, + int(old_record.get("n_pixels", 0)), + int(old_record.get("n_coverage_pixels", 0))) + stale = record != old_record + if stale: + write_sidecar(args.sidecar, record) + print(f"[merge_defect_map] unchanged: {args.output} " + f"({len(have)} exposure(s)" + f"{'; sidecar refreshed' if stale else ''})") + return + record = apply_plan(args.output, args.sidecar, plan, have, digests, + missing, args.nside_coverage, args.nside) + warn = (f"; {len(missing)} campaign exposure(s) have no fragment" + if missing else "") + print(f"[merge_defect_map] {plan.describe()} -> {args.output} " + f"({record['n_pixels']} healpix pixel(s) over " + f"{record['n_coverage_pixels']} coverage pixel(s)){warn}") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/nexp_map.py b/workflow/scripts/nexp_map.py new file mode 100644 index 000000000..bf996e429 --- /dev/null +++ b/workflow/scripts/nexp_map.py @@ -0,0 +1,173 @@ +#!/usr/bin/env python3 +"""Build the campaign's HealSparse exposure-count (nexp) map from the exposure footprints. + +Run as the shell of the in-DAG ``nexp_map`` rule, never by hand. + +WHAT THE MAP IS. Per sky pixel, the number of exposures with a VALID PSF MODEL +covering it. The CCDs of one exposure do not overlap, so stamping value 1 per CCD +polygon and accumulating counts exposures. It is the one coverage count there +is: sp_validation's ``npoint >= 3`` cut reads exactly this, and applies it as a +structural mask (``notebooks/demo_apply_hsp_masks.py``). + +CAMPAIGN-CUMULATIVE, AND THAT IS THE POINT. This script GLOBS every +``/exp/*/*/manifests/exp_footprint.json`` — not just the ones the +rule declared as inputs, and INCLUDING exposures whose scratch stores have been +reclaimed. The declared inputs are the in-scope footprints of live stores, which +buys ordering without dragging out-of-scope tiles or reclaimed chains into the +DAG, and the rule's params fingerprint the ids of every record this glob will +find, which is what reruns the map when one arrives off the DAG; the records +themselves live on the persistent root and stay valid sky forever. So +appending tiles GROWS the map rather than replacing it, which is what a survey +coverage mask should do. The rule's comment says the same thing +where a reader of the DAG will meet it. + +THE STAMPING IS NOT HERE. It is ``shapepipe.utilities.coverage_map_builder. +build_map``, which is where the RA-seam guard, the pole guard and the polygon +accumulation live. This script is the JSON-to-arrays half, and the records are +raw sky: unwrapping across RA=0 happens once, inside build_map, at the moment a +polygon is stamped. + +NSIDE IS NOT DEFAULTED HERE. Both values are required arguments, carried from +`exposure_maps:` in config.yaml — the one resolution the defect map shares — +and the difference they make is invisible in the output: nside=131072 is +~1.6"/pixel, chosen to match the UNIONS bit-mask resolution so coverage and mask +align pixel-wise. A map built coarser would look entirely reasonable and would +not align, and the consumer would not notice. +""" + +import argparse +import filecmp +import hashlib +import json +import sys +from pathlib import Path + +import numpy as np + +from shapepipe.utilities.coverage_map_builder import build_map + +# Where exp_footprint writes, relative to the products root. The shard and the +# exposure id are both globbed: this map is the whole campaign's, so it is +# deliberately not built from a list of units. +FOOTPRINT_GLOB = "exp/*/*/manifests/exp_footprint.json" + + +def read_footprints(products_dir): + """Every exposure footprint on the persistent root, as flat arrays. + + Returns ``(ccd_ids, ra, dec, units)``: length-M id and ``(M, 4)`` corner + arrays over all CCDs of all exposures, plus the sorted exposure ids that + contributed. Records are read in sorted path order so the polygon order — + and hence the map — does not depend on readdir order. + """ + paths = sorted(Path(products_dir).glob(FOOTPRINT_GLOB)) + if not paths: + sys.exit(f"nexp_map: no {FOOTPRINT_GLOB} under {products_dir}; " + f"there is nothing to build a map from") + + ccd_ids, ra, dec, units = [], [], [], [] + for path in paths: + body = json.loads(path.read_text()) + units.append(body["unit"]) + for ccd in body["ccds"]: + ccd_ids.append(ccd["id"]) + ra.append(ccd["ra"]) + dec.append(ccd["dec"]) + + # Records but no CCDs is the other empty map, and it is the quieter one: it + # means every exposure on the root lost every CCD, which is a broken PSF + # stage rather than a survey with no coverage. Written out it would be a + # valid, plausible-looking .hsp full of nothing, and its consumer masks + # everything. + if not ccd_ids: + sys.exit(f"nexp_map: {len(paths)} footprint record(s) under " + f"{products_dir}, and not one names a CCD with a PSF model; " + f"there is nothing to stamp") + + return (np.array(ccd_ids, dtype=str), + np.array(ra, dtype=float).reshape(-1, 4), + np.array(dec, dtype=float).reshape(-1, 4), + sorted(units)) + + +def write_stable(path, body): + """Write JSON tmp-then-``cmp``-then-``mv``; an unchanged body keeps its mtime.""" + path.parent.mkdir(parents=True, exist_ok=True) + tmp = path.with_name(path.name + ".tmp") + try: + tmp.write_text(json.dumps(body, indent=2, sort_keys=True) + "\n") + if path.exists() and filecmp.cmp(tmp, path, shallow=False): + tmp.unlink() + else: + tmp.replace(path) + finally: + tmp.unlink(missing_ok=True) + + +def publish_map(hsp_map, out): + """Write the map to a sibling temporary, then rename it over ``out``. + + The rename is the publication: a job killed before it leaves the previous + map in place, with the manifest that describes it. Returns the published + file's sha256, which the manifest records to bind itself to this map. + """ + out.parent.mkdir(parents=True, exist_ok=True) + tmp = out.with_name(f".{out.name}.tmp") + try: + hsp_map.write(str(tmp), clobber=True) + digest = hashlib.sha256() + with open(tmp, "rb") as f: + for block in iter(lambda: f.read(1 << 24), b""): + digest.update(block) + tmp.replace(out) + finally: + tmp.unlink(missing_ok=True) + return digest.hexdigest() + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--products-dir", required=True, type=Path, + help="the persistent root; every exposure footprint under " + "it goes into the map") + p.add_argument("--out", required=True, type=Path, + help="the HealSparse map, " + "/nexp_map/nexp_map_.hsp") + p.add_argument("--manifest", required=True, type=Path) + p.add_argument("--nside-coverage", required=True, type=int) + p.add_argument("--nside", required=True, type=int) + p.add_argument("--verbose", action="store_true") + args = p.parse_args() + + ccd_ids, ra, dec, units = read_footprints(args.products_dir) + print(f"[nexp_map] {len(units)} exposure(s), {len(ccd_ids)} CCD " + f"footprint(s) from {args.products_dir}") + + hsp_map = build_map(ccd_ids, ra, dec, args.nside_coverage, args.nside, + verbose=args.verbose) + + map_sha256 = publish_map(hsp_map, args.out) + + # The manifest is written AFTER the map, and it is the rule's record of + # which exposures the map contains — the question a mask's consumer asks + # months later, and one nothing else on disk can answer once the campaign + # has grown past it. Its digest names the map it describes, so a map and + # manifest from different generations disagree visibly. + write_stable(args.manifest, { + "stage": "nexp_map", "level": "campaign", "status": "complete", + "map": str(args.out), + "map_sha256": map_sha256, + "nside_coverage": args.nside_coverage, + "nside": args.nside, + # What the rule's memory request is sized from on the next build + # (map_cov_pixels(), Snakefile). + "n_coverage_pixels": int(hsp_map.coverage_mask.sum()), + "n_exposures": len(units), + "n_ccds": len(ccd_ids), + "exposures": units, + }) + print(f"[nexp_map] -> {args.out}") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/run_config.py b/workflow/scripts/run_config.py index 831147986..940cd21fa 100644 --- a/workflow/scripts/run_config.py +++ b/workflow/scripts/run_config.py @@ -127,6 +127,36 @@ def _dollar_keys(value, prefix): return [prefix] if isinstance(value, str) and "$" in value else [] +# Config blocks the workflow no longer reads, and what replaced each. A key +# nothing reads is a switch that silently does nothing (`coverage.enabled: +# false` would leave the exposure-count map on), so retired() makes it an error. +RETIRED = { + "coverage": "exposure_maps.nexp (on/off: exposure_maps.nexp.enabled; " + "resolution: the shared exposure_maps.nside / " + "exposure_maps.nside_coverage)", + "defect_map": "exposure_maps.defect (oversample: " + "exposure_maps.defect.oversample; resolution: the shared " + "exposure_maps.nside / exposure_maps.nside_coverage)", +} + + +def retired(config): + """Every place a RETIRED key appears — top level, a `machines:` entry, or + one of its per-input_type blocks — as `(path, replacement)` pairs.""" + found = [(key, RETIRED[key]) for key in RETIRED if key in config] + for machine, entry in (config.get("machines") or {}).items(): + if not isinstance(entry, dict): + continue + found += [(f"machines.{machine}.{key}", RETIRED[key]) + for key in RETIRED if key in entry] + for input_type, block in entry.items(): + if isinstance(block, dict): + found += [(f"machines.{machine}.{input_type}.{key}", + RETIRED[key]) + for key in RETIRED if key in block] + return found + + def unresolved(config): """REQUIRED keys that are unset, the placeholder, or hold an unexpanded $variable, then every MACHINE_KEYS value (recursively through diff --git a/workflow/scripts/run_report.py b/workflow/scripts/run_report.py index 11cdcffe4..2feb2d09d 100644 --- a/workflow/scripts/run_report.py +++ b/workflow/scripts/run_report.py @@ -61,12 +61,12 @@ "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. +# exp_persist, exp_footprint and exp_defect_map are DELIBERATELY NOT in that +# list. This report disk-scans the scratch run_dir, and those are the exposure +# manifests that live on products_dir instead — that placement is what makes +# them survive clean_exposure. Listed here they 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