diff --git a/astra.yaml b/astra.yaml index 0e9f4867e..46617f908 100644 --- a/astra.yaml +++ b/astra.yaml @@ -1421,22 +1421,25 @@ analyses: defect_fill: label: Image content of defect pixels before metacal rationale: >- - A defect pixel (nonzero instrument flag, zero exposure weight or - invalid background RMS) gets weight 0 in prepare_ngmix_weights. Its + A defect pixel (nonzero instrument flag, zero exposure weight, + invalid background RMS, or off-tile) gets weight 0 in prepare_ngmix_weights. Its image value still matters: ngmix's metacal deconvolves, shears and reconvolves an InterpolatedImage of the whole image and copies the weights through, so a zero-weight pixel's content spreads into the weighted pixels within about a PSF width. DES's ngmixer fills for that reason: "it may be important for codes that take moments or use - FFTs". In the committed default (`BLEND_HANDLING = noisefill`), - masked pixels are replaced with independent noise at the per-pixel - background RMS when supplied, or the stamp's robust noise scale - otherwise. With `BLEND_HANDLING = uberseg`, the image is left - untouched and weights are zeroed on defects and neighbour-side - pixels. [LINT] the prepare_ngmix_weights docstring says noisefill - keeps the weight of filled pixels (the code zeroes it), and the - ngmix_runner comment says noisefill fills neighbour pixels (it fills - flagged pixels and leaves neighbours untouched). The committed fill + FFTs". Under every BLEND_HANDLING, the committed uberseg included, + and the default DEFECT_FILL = + noise, defect pixels are replaced with independent noise at the + per-pixel background RMS when supplied, or the stamp's robust noise + scale otherwise. The tile VIGNET holds -1e30 on other detections' + footprints and beyond the tile's edge (written by SExtractor). The runs of entirely -1e30 stamp + rows and columns that start at a stamp border (the off-image part of + a rectangle clip) are off-tile: the epoch holds the object's own + light there, so they are defects, flagged 2**10. The other markers, + including a neighbour footprint that completes an interior row + beside an off-tile band, are neighbour pixels, not defects; + blend_handling decides their treatment. The committed fill uses the unsymmetrized defect set: DES symmetrized its masks, but four-fold symmetrization quadruples m and still leaves an additive c1 (symmetrized_4fold_noise). The cost of not symmetrizing, a hole in @@ -1458,7 +1461,11 @@ analyses: interpolate: label: Interpolate short bounded runs; noise-fill the rest description: >- - Not implemented. PR #916 measures c1 = -1.3e-3 for a column + DEFECT_FILL = interpolate: row or column runs of at most + MAX_INTERPOLATED_RUN = 3 defect pixels that stop short of the + stamp border are Clough-Tocher interpolated from the clean pixels + within SUPPORT_RADIUS = 4 px; other defects are noise-filled. + PR #916 measures c1 = -1.3e-3 for a column 8 px from a 0.5 arcsec galaxy and +2e-6 when the weight is also zeroed on its quarter-turn orbit. For a 3-px bleed 6 px from that galaxy, PR #916 measures m11 = +0.89% when the fill is @@ -1476,11 +1483,11 @@ analyses: with symmetrization against 0.7e-4 without it. insights: [mask_bad_column_symmetrize, mask_des_defect_practice, mask_fixed_orientation] raw: - label: "No fill: raw defect values (BLEND_HANDLING = uberseg)" + label: "No fill: raw defect values" description: >- - With BLEND_HANDLING = uberseg, defect pixels keep weight 0 but - their image values stay untouched in the image that metacal - deconvolves, shears and reconvolves. + Defect pixels keep weight 0 but their image values stay + untouched in the image that metacal deconvolves, shears and + reconvolves. excluded: true excluded_reason: >- Metacal acts on every pixel regardless of weight, so raw defects @@ -1490,12 +1497,15 @@ analyses: blend_handling: label: Neighbour treatment before metacal rationale: >- - Covers only pixels shared with a neighbour; defect_fill is coupled to - it through BLEND_HANDLING. noisefill (the default; the committed - config sets no key) leaves neighbours fully weighted and untouched. - uberseg zeroes the weight of pixels nearer a neighbour's coadd - segmentation footprint than the target's (DILATE_NEIGHBOUR, default 1) - and leaves the image untouched, as official uberseg does: + Covers only pixels shared with a neighbour; defect_fill does not + depend on it: defects are filled the same way under every + BLEND_HANDLING, and the epoch cuts never count neighbour pixels. + The committed choice is uberseg (the workflow's default + blend_handling), which ignores the -1e30 neighbour markers in the + tile VIGNET and zeroes the weight of pixels nearer a neighbour's + coadd segmentation footprint than the target's (DILATE_NEIGHBOUR, + default 1), leaving the neighbour side's image untouched, as + official uberseg does: esheldon/meds get_uberseg returns a weight map (a nearest-segment-pixel Voronoi split). DES Y1's fiducial metacal and the last Y3 config ran on uberseg-weighted stamps with raw neighbour @@ -1512,24 +1522,31 @@ analyses: with simulations (mask_des_y1_uberseg_only, mask_blend_bias_detection). Noise-filling the neighbour side would instead cut the target's own light along an unsheared boundary, a - sharp edge that rings in the FFTs. The recommended comparison arm is - uberseg (weight-only), with defect_fill held equal across arms. - default: none + sharp edge that rings in the FFTs. The alternative, noisefill + (ngmix's own default, and the workflow's blend_handling: + noisefill), gives weight 0 to the pixels marked -1e30 in the tile + VIGNET on other detections' footprints (about 93% of + neighbour-footprint pixels on a simulated tile), excluding the + off-tile rows and columns, which are defects, and replaces them with + noise; unmarked neighbour pixels keep their weight and light. It is + the comparison arm, with defect_fill held equal across arms. + default: uberseg options: - none: - label: No neighbour treatment (BLEND_HANDLING = noisefill) + noisefill: + label: Noise-fill the marked neighbour pixels (BLEND_HANDLING = noisefill) description: >- - Neighbour pixels keep their full weight and image values, so the - fit sees all neighbour light, which biases shapes toward - neighbours (Jarvis et al. 2016); this masks less than even their - plain segmentation map. A candidate cause of the FLAGS=2 B-modes + Pixels marked -1e30 in the tile VIGNET get weight 0 and noise. + The fill stops at the marked footprint, so unmarked neighbour + pixels and the neighbour's wings keep their weight and light, + which biases shapes toward neighbours (Jarvis et al. 2016). A + candidate cause of the FLAGS=2 B-modes investigated in #814. insights: [mask_uberseg_neighbour_bias] uberseg: label: UberSeg, weight-only description: >- - As in DES Y1/Y3. Needs the coadd segmentation stamp - (SEG_VIGNET_PATH). DILATE_NEIGHBOUR absorbs the coadd-vs-epoch + As in DES Y1/Y3. Needs the coadd segmentation stamp (the tile + catalogue's SEG_VIGNET). DILATE_NEIGHBOUR absorbs the coadd-vs-epoch overlay offset: ShapePipe reuses one coadd seg stamp for every epoch where MEDS reprojects it. insights: [mask_uberseg_weight_only, mask_uberseg_neighbour_bias, mask_des_y1_uberseg_only, mask_blend_bias_detection] @@ -1553,11 +1570,13 @@ analyses: central_defect_veto: label: Per-epoch veto on a defect near the stamp centre rationale: >- - Not on develop; implemented on feat/defect-fill-veto (7777181b). - There an epoch is dropped when a defect pixel lies strictly closer - to the stamp centre than its fill's radius, beside the - masked-fraction cut in the epoch loop. The veto reads only the - defect mask, so it selects on nothing shear-responsive; for the same + An epoch is dropped when a defect pixel lies strictly closer to the + stamp centre than its fill's radius, beside the masked-fraction cut + in the epoch loop; EPOCH_CENTRAL_DEFECT_RADIUS = 0 disables it. + The tile VIGNET's -1e30 neighbour markers are not defects: all epochs + share the tile VIGNET, so a neighbour inside the radius would drop + every epoch. Off-tile pixels are defects, so an object near the + tile edge is vetoed. The veto reads only the defect mask, so it selects on nothing shear-responsive; for the same reason the radii are fixed rather than scaled by galaxy size: 10 px for noise-filled pixels (EPOCH_CENTRAL_DEFECT_RADIUS) and 7 px for interpolated ones (EPOCH_INTERPOLATED_DEFECT_RADIUS, @@ -1568,16 +1587,18 @@ analyses: wide defects through the elliptical PSF sit at the 1% bound (m11 = -0.98%; -0.24% at 11 px), and noise fill needs 14 px for 0.7 and 0.9 arcsec galaxies. Measured on feat/defect-fill-veto and - feat/defect-interpolation. - default: disabled + feat/defect-interpolation. The 7 px radius applies only under + DEFECT_FILL = interpolate; with the committed noise fill every defect + takes the 10 px radius. + Values: + EPOCH_CENTRAL_DEFECT_RADIUS = 10; + EPOCH_INTERPOLATED_DEFECT_RADIUS = 7. + default: fixed_radii options: disabled: - label: No central veto (committed code) + label: No central veto (EPOCH_CENTRAL_DEFECT_RADIUS = 0) fixed_radii: label: Fixed radii, 10 px for noise-filled and 7 px for interpolated defects - description: >- - Implemented on feat/defect-fill-veto and feat/defect-interpolation, - not on develop. size_scaled_radius: label: Veto radius scaled by galaxy size excluded: true @@ -1587,14 +1608,16 @@ analyses: epoch_masked_fraction_cut: label: Per-epoch masked-fraction cut rationale: >- - [HARDCODED] on develop, an epoch whose stamp has more than 1/3 of its - pixels flagged (any nonzero flag bit, including the tile-coverage bit - 2**10 set where the tile vignet is off-image) is dropped from the - multi-epoch fit; an object with no surviving epoch has no shape. - Zero-weight and invalid-RMS pixels are not counted. On - feat/defect-fill-veto the cut counts the raw, unsymmetrized defect - set (flagged, zero-weight and invalid-RMS pixels) against - EPOCH_MASKED_FRACTION_CUT, default 1/3. Before the cut, an epoch is + An epoch whose stamp has more than EPOCH_MASKED_FRACTION_CUT + (default 1/3) of its pixels in the raw, unsymmetrized defect set is + dropped from the multi-epoch fit: flagged pixels (any nonzero + exposure flag bit), zero-weight pixels, invalid-RMS pixels and + off-tile pixels (whole -1e30 rows and columns at the tile VIGNET's + border, flag 2**10), the set defect_fill fills. On a 51-px stamp an object + within about 8.5 px of the tile edge fails the 1/3 cut. The other + -1e30 markers, neighbour footprints, are not counted: every epoch + shares the tile VIGNET, so a large neighbour would drop them all. An object with no surviving epoch has no + shape. Before the cut, an epoch is dropped silently if its galaxy stamp is all zeros or its background-subtracted sigma_mad is not positive. DES was stricter: Y1 rejected any epoch with a masked or zero-weight pixel, or with @@ -1610,12 +1633,12 @@ analyses: label: 1/3 of the stamp in the defect set description: >- Drop an epoch only if more than 1/3 of the stamp pixels are - defects (on develop, flagged pixels). + defects. ten_percent: label: 10% (DES Y3 / Y6) description: >- - Not a default; EPOCH_MASKED_FRACTION_CUT = 0.1 on - feat/defect-fill-veto. Matches DES Y3 max_zero_weight_frac and + Not a default; EPOCH_MASKED_FRACTION_CUT = 0.1. Matches DES Y3 + max_zero_weight_frac and Y6 max_masked_fraction. insights: [mask_multi_epoch_drop] any_masked: diff --git a/src/shapepipe/modules/make_cat_package/make_cat.py b/src/shapepipe/modules/make_cat_package/make_cat.py index 8a6f9548a..0c3a84dd6 100644 --- a/src/shapepipe/modules/make_cat_package/make_cat.py +++ b/src/shapepipe/modules/make_cat_package/make_cat.py @@ -113,7 +113,8 @@ def save_sextractor_data(final_cat_file, sexcat_path, remove_vignet=True): sexcat_path : str Path to SExtractor catalogue to save remove_vignet : bool - If ``True`` will not save the ``VIGNET`` field into the final catalogue + If ``True`` will not save the ``VIGNET`` and ``SEG_VIGNET`` stamp + fields into the final catalogue Returns ------- @@ -127,6 +128,7 @@ def save_sextractor_data(final_cat_file, sexcat_path, remove_vignet=True): data = np.copy(sexcat_file.get_data()) if remove_vignet: data = remove_field_name(data, "VIGNET") + data = remove_field_name(data, "SEG_VIGNET") cat_size = len(data) tile_name = os.path.basename(sexcat_path) diff --git a/src/shapepipe/modules/ngmix_package/defect_interpolation.py b/src/shapepipe/modules/ngmix_package/defect_interpolation.py new file mode 100644 index 000000000..ae0d53c41 --- /dev/null +++ b/src/shapepipe/modules/ngmix_package/defect_interpolation.py @@ -0,0 +1,168 @@ +"""DEFECT INTERPOLATION. + +Clough-Tocher interpolation of short defect runs before metacal, selected +with ``DEFECT_FILL = interpolate`` (see +:func:`shapepipe.modules.ngmix_package.ngmix.prepare_ngmix_weights`). + +""" + +import numpy as np +from scipy.interpolate import CloughTocher2DInterpolator +from scipy.ndimage import binary_dilation, label +from scipy.spatial import QhullError + +# Longest row or column run of defect pixels that is interpolated. +MAX_INTERPOLATED_RUN = 3 + +# Clean pixels within this Chebyshev distance (pixels) of the interpolated +# pixels support the interpolant. +SUPPORT_RADIUS = 4 + +_ROW_RUNS = np.array([[0, 0, 0], [1, 1, 1], [0, 0, 0]]) + + +def _short_row_runs(defect, max_run): + """Pixels in row runs of at most ``max_run`` defects that stop short of + both stamp borders.""" + labels, n_runs = label(defect, structure=_ROW_RUNS) + if n_runs == 0: + return np.zeros_like(defect) + short = np.bincount(labels.ravel(), minlength=n_runs + 1) <= max_run + short[0] = False + short[labels[:, 0]] = False + short[labels[:, -1]] = False + return short[labels] + + +def interpolable_defects(defect, max_run=MAX_INTERPOLATED_RUN): + """Defect pixels that ``DEFECT_FILL = interpolate`` interpolates. + + @sc [decision:shape_measurement.defect_fill] interpolable-defects + A defect pixel is interpolated when its row or its column run of defect + pixels is at most ``max_run`` (3) long and has clean pixels at both + ends. That covers columns, 3-px bleeds and isolated pixels, the widths + whose shear recovery is calibrated. Wider holes, and runs that reach the + stamp border (edge bands, corners), have clean light on one side only; + they are noise-filled and vetoed at the noise-fill radius (see + :func:`~shapepipe.modules.ngmix_package.ngmix.central_defect_vetoes`). + The rule reads only the mask and commutes with quarter turns of the + stamp. + + Parameters + ---------- + defect : numpy.ndarray of bool + Defect mask of one epoch stamp. + max_run : int, optional + Longest interpolated run; the default is ``MAX_INTERPOLATED_RUN``. + + Returns + ------- + numpy.ndarray of bool + ``True`` on the defect pixels to interpolate. + """ + defect = np.asarray(defect, dtype=bool) + return ( + _short_row_runs(defect, max_run) + | _short_row_runs(defect.T, max_run).T + ) + + +def fourfold(mask): + """Union of a square stamp mask with its quarter turns. + + Parameters + ---------- + mask : numpy.ndarray of bool + Square mask. + + Returns + ------- + numpy.ndarray of bool + ``mask`` ORed with its rotations by 90, 180 and 270 degrees about + the stamp centre. + + Raises + ------ + ValueError + If ``mask`` is not square. + """ + mask = np.asarray(mask, dtype=bool) + if mask.ndim != 2 or mask.shape[0] != mask.shape[1]: + raise ValueError( + f"A quarter-turn orbit needs a square stamp, not {mask.shape}" + ) + return mask | np.rot90(mask) | np.rot90(mask, 2) | np.rot90(mask, 3) + + +def _interpolate_once(planes, excluded, target): + """One Clough-Tocher interpolant of every plane at ``target``; NaN + elsewhere and where the support cannot reach.""" + out = np.full(planes.shape, np.nan) + support = binary_dilation( + target, structure=np.ones((3, 3), dtype=bool), + iterations=SUPPORT_RADIUS, + ) & ~excluded + points = np.argwhere(support).astype(float) + if len(points) < 3: + return out + query = np.argwhere(target) + try: + interpolant = CloughTocher2DInterpolator( + points, planes[:, support].T, fill_value=np.nan, + ) + except QhullError: + return out + out[:, query[:, 0], query[:, 1]] = interpolant(query.astype(float)).T + return out + + +def interpolate_defects(planes, excluded, target): + """Replace the ``target`` pixels of every plane by a Clough-Tocher + interpolant of the kept pixels around them. + + @sc [decision:shape_measurement.defect_fill] shared-rotation-averaged-interpolant + The support is the pixels within ``SUPPORT_RADIUS`` (4 px) of the target + outside ``excluded``; no excluded pixel enters it, so their values are + never read. For each quarter turn of the stamp, one Delaunay triangulation of the support + serves every plane, so the science image and the metacal noise image + see the same linear operator and fixnoise mirrors the science image's + interpolated noise. A regular grid's triangulation has degenerate + diagonals, so one orientation has a preferred direction; averaging the + four quarter-turned operators makes the fill commute with rotations of + the stamp. That removes up to 6e-5 of c2 for a column or 3-px bleed + 6 px from the object. + + Parameters + ---------- + planes : array_like + Stamp planes, shape ``(n, ny, nx)``. + excluded : numpy.ndarray of bool + Pixels that never support the interpolant, shape ``(ny, nx)``: every + defect, and any pixel whose light the image does not keep. + target : numpy.ndarray of bool + Defect pixels to interpolate (:func:`interpolable_defects`). + + Returns + ------- + numpy.ndarray + A copy of ``planes`` with ``target`` pixels interpolated; NaN at a + target pixel whose support is degenerate in some orientation. + """ + planes = np.asarray(planes, dtype=float) + excluded = np.asarray(excluded, dtype=bool) + target = np.asarray(target, dtype=bool) + out = planes.copy() + if not target.any(): + return out + turns = [ + np.rot90( + _interpolate_once( + np.rot90(planes, k, axes=(1, 2)), np.rot90(excluded, k), + np.rot90(target, k), + ), + -k, axes=(1, 2), + )[:, target] + for k in range(4) + ] + out[:, target] = np.mean(turns, axis=0) + return out diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 12b192fd8..937007be3 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -8,6 +8,7 @@ import os import re +from collections import Counter from typing import NamedTuple import ngmix @@ -21,11 +22,40 @@ from scipy.spatial import cKDTree from sqlitedict import SqliteDict +from shapepipe.modules.ngmix_package.defect_interpolation import ( + fourfold, + interpolable_defects, + interpolate_defects, +) from shapepipe.pipeline import file_io +# Flag bit set on an epoch's off-tile pixels (see :func:`split_tile_markers`). +OFF_TILE_FLAG = 2**10 + # Neighbour treatments selectable with the BLEND_HANDLING option. BLEND_HANDLINGS = ("noisefill", "uberseg") +# Defect fills selectable with the DEFECT_FILL option (see +# :func:`prepare_ngmix_weights`). +DEFECT_FILLS = ("noise", "interpolate") + +# Default of the EPOCH_MASKED_FRACTION_CUT option: an epoch is dropped when +# more than this fraction of its stamp is in :func:`defect_mask` (see +# :func:`prepare_postage_stamps`). +EPOCH_MASKED_FRACTION_CUT = 1 / 3 + +# Default of the EPOCH_CENTRAL_DEFECT_RADIUS option (pixels): an epoch is +# dropped when a pixel of :func:`defect_mask` lies closer than this to the +# stamp centre (see :func:`has_central_defect`). +# @sc [decision:shape_measurement.central_defect_veto] +EPOCH_CENTRAL_DEFECT_RADIUS = 10 + +# Default of the EPOCH_INTERPOLATED_DEFECT_RADIUS option (pixels): under +# DEFECT_FILL = interpolate, the veto radius for interpolated defect pixels +# (see :func:`central_defect_vetoes`). +# @sc [decision:shape_measurement.central_defect_veto] +EPOCH_INTERPOLATED_DEFECT_RADIUS = 7 + # @sc [decision:shape_measurement.metacal_scheme] METACAL_TYPES = ('noshear', '1p', '1m', '2p', '2m') @@ -508,22 +538,17 @@ class Tile_cat(): Parameters ---------- cat_path : str - Path to the tile SExtractor catalogue. - seg_cat_path : str, optional - Path to the coadd-frame segmentation VIGNET catalogue (a CLASSIC-mode - vignetmaker output cut from the tile ``SEGMENTATION`` check image), - row-aligned to ``cat_path``. When given, ``self.seg`` holds one integer - seg stamp per object for the ``"uberseg"`` blend handling; ``None`` - leaves ``self.seg`` unset (the noise-fill path is unaffected). + Path to the tile SExtractor catalogue. Its optional ``SEG_VIGNET`` + column, one integer coadd segmentation stamp per object on the grid + of its ``VIGNET``, becomes ``self.seg`` for the ``"uberseg"`` blend + handling; without it ``self.seg`` is ``None``. """ def __init__( self, cat_path, - seg_cat_path=None, ): self.cat_path = cat_path - self.seg_cat_path = seg_cat_path if cat_path: self.get_data(cat_path) @@ -540,44 +565,20 @@ def get_data(self, cat_path): self.ra = np.copy(data['XWIN_WORLD']) self.dec = np.copy(data['YWIN_WORLD']) - # Optional columns — may be absent in external (non-SExtractor) catalogs + # Optional columns — may be absent in external (non-SExtractor) catalogs. + # The stamp columns are the bulk of the table and are views into it, + # not copies, so it is held once (prepare_postage_stamps copies each + # object's stamp before changing it). self.flux = np.copy(data['FLUX_AUTO']) if 'FLUX_AUTO' in cols else None - self.vign = np.copy(data['VIGNET']) if 'VIGNET' in cols else None + self.vign = data['VIGNET'] if 'VIGNET' in cols else None - tile_cat.close() + # Coadd-frame segmentation stamp (integer labels, the catalogue's + # NUMBER), one per object on the grid of its VIGNET, overlaid + # unchanged on every epoch for uberseg neighbour masking + # (shapepipe#776). + self.seg = data['SEG_VIGNET'] if 'SEG_VIGNET' in cols else None - # Coadd-frame SExtractor segmentation stamp (integer labels), one per - # object and row-aligned to the tile catalogue, overlaid unchanged on - # every epoch for uberseg neighbour masking (shapepipe#776). None -> - # uberseg unavailable; the noise-fill path is unaffected. - self.seg = None - if self.seg_cat_path: - seg_cat = file_io.FITSCatalogue( - self.seg_cat_path, - SEx_catalogue=True, - ) - seg_cat.open() - seg_data = seg_cat.get_data() - # The seg VIGNETs are indexed by tile-catalogue row (self.seg[i]), - # so the two catalogues MUST be row-aligned. Fail loud at load if - # they are not — a silent length/order mismatch would hand every - # object the wrong footprint and quietly corrupt every mask. - if len(seg_data) != len(self.obj_id): - raise ValueError( - f"SEG_VIGNET_PATH '{self.seg_cat_path}' has" - + f" {len(seg_data)} rows but the tile catalogue has" - + f" {len(self.obj_id)}; the segmentation vignets must be" - + " row-aligned to the tile catalogue." - ) - if 'NUMBER' in seg_data.dtype.names: - if not np.array_equal(seg_data['NUMBER'], self.obj_id): - raise ValueError( - f"SEG_VIGNET_PATH '{self.seg_cat_path}' NUMBER column" - + " does not match the tile catalogue NUMBER; the" - + " segmentation vignets are misaligned or reordered." - ) - self.seg = np.copy(seg_data['VIGNET']) - seg_cat.close() + tile_cat.close() class Postage_stamp(): """Galaxy Postage Stamp. @@ -602,9 +603,15 @@ def __init__( self.psfs = [] self.weights = [] self.flags = [] + # Neighbour masks, one per epoch: the pixels marked -1e30 in the tile + # VIGNET on other detections' footprints (off-tile markers + # are flagged as defects instead; see split_tile_markers), + # MegaCam-flipped to the epoch. noisefill zero-weights and noise-fills them; uberseg and + # the epoch cuts do not read them (see prepare_ngmix_weights). + self.neighbours = [] self.bkg_rms = [] # Segmentation stamps, one per epoch, used only by the "uberseg" blend - # handling; empty for the default noise-fill path. All epochs carry the + # handling; empty under the default "noisefill". All epochs carry the # SAME coadd-frame seg stamp (shapepipe#776: one coadd seg per object, # no per-epoch reprojection), each MegaCam-flipped to match its galaxy # stamp so the overlay stays registered. @@ -621,6 +628,9 @@ def __init__( # CCD number of the first epoch, used only to build the per-object # position seed (see :func:`position_seed`). self.ccd = None + self.epoch_cuts = Counter( + considered=0, masked_fraction=0, central_veto=0 + ) self.bkg_sub = bkg_sub self.megacam_flip = megacam_flip @@ -708,15 +718,32 @@ class Ngmix(object): adaptive-moment centroid measured from the stamp pixels. See :func:`make_ngmix_observation`. blend_handling : {"noisefill", "uberseg"}, optional - Neighbour treatment; ``"noisefill"`` (default) is the historical - noise-fill, ``"uberseg"`` hard-masks neighbour-side pixels from the - coadd segmentation map and requires ``seg_cat_path``. - seg_cat_path : str, optional - Path to the coadd-frame segmentation VIGNET catalogue (see - :class:`Tile_cat`). Required when ``blend_handling="uberseg"``. + Neighbour treatment. ``"noisefill"`` (default) zero-weights and + noise-fills the pixels marked -1e30 in the tile VIGNET on other + detections' footprints; ``"uberseg"`` ignores those markers and + zeroes the weight of neighbour-side pixels from the coadd + segmentation stamps, the tile catalogue's ``SEG_VIGNET`` column (see + :class:`Tile_cat`), which it requires. Defect pixels are filled + under both (see :func:`prepare_ngmix_weights`). dilate_neighbour : int, optional Neighbour-mask dilation iterations for ``"uberseg"`` (see :func:`uberseg_weight`); the default is ``1``. + epoch_central_defect_radius : float, optional + Drop an epoch when a defect pixel lies closer than this many pixels + to the stamp centre (see :func:`has_central_defect`); the default is + ``EPOCH_CENTRAL_DEFECT_RADIUS``. + epoch_masked_fraction_cut : float, optional + Drop an epoch when more than this fraction of its stamp is defects + (see :func:`prepare_postage_stamps`); the default is + ``EPOCH_MASKED_FRACTION_CUT``. + defect_fill : {"noise", "interpolate"}, optional + How defect pixels are filled before metacal (see + :func:`prepare_ngmix_weights`); the default is ``"noise"``. + epoch_interpolated_defect_radius : float, optional + Under ``defect_fill="interpolate"``, drop an epoch when an + interpolated defect pixel lies closer than this many pixels to the + stamp centre (see :func:`central_defect_vetoes`); the default is + ``EPOCH_INTERPOLATED_DEFECT_RADIUS``. Notes ----- @@ -728,8 +755,7 @@ class Ngmix(object): IndexError If the length of the input file list is incorrect ValueError - If ``blend_handling`` is unknown, or ``"uberseg"`` is selected without - ``seg_cat_path``. + If ``blend_handling`` or ``defect_fill`` is unknown. """ @@ -747,9 +773,12 @@ def __init__( bkg_sub=True, centroid_source="wcs", blend_handling="noisefill", - seg_cat_path=None, dilate_neighbour=1, metacal_psf="fitgauss", + epoch_central_defect_radius=EPOCH_CENTRAL_DEFECT_RADIUS, + epoch_masked_fraction_cut=EPOCH_MASKED_FRACTION_CUT, + defect_fill="noise", + epoch_interpolated_defect_radius=EPOCH_INTERPOLATED_DEFECT_RADIUS, ): # Base count = catalogue + vignets, excluding the f_wcs headers (passed @@ -768,13 +797,10 @@ def __init__( f"Unknown BLEND_HANDLING '{blend_handling}'; expected one of" + f" {BLEND_HANDLINGS}" ) - - # Fail fast at construction (not deep in the per-epoch loop) when - # uberseg is requested without its segmentation input (shapepipe#776). - if blend_handling == "uberseg" and seg_cat_path is None: + if defect_fill not in DEFECT_FILLS: raise ValueError( - "blend_handling='uberseg' requires SEG_VIGNET_PATH (the coadd" - + " SExtractor segmentation vignets); none configured." + f"Unknown DEFECT_FILL '{defect_fill}'; expected one of" + + f" {DEFECT_FILLS}" ) self._tile_cat_path = input_file_list[0] @@ -816,9 +842,14 @@ def __init__( self._bkg_sub = bkg_sub self._centroid_source = centroid_source self._blend_handling = blend_handling - self._seg_cat_path = seg_cat_path self._dilate_neighbour = dilate_neighbour self._metacal_psf = metacal_psf + self._epoch_central_defect_radius = epoch_central_defect_radius + self._epoch_masked_fraction_cut = epoch_masked_fraction_cut + self._defect_fill = defect_fill + self._epoch_interpolated_defect_radius = ( + epoch_interpolated_defect_radius + ) self._w_log = w_log @@ -1200,9 +1231,18 @@ def process(self): @sc [decision:shape_measurement.fit_initialisation,decision:shape_measurement.ngmix_seed_mode] """ - tile_cat = Tile_cat(self._tile_cat_path, self._seg_cat_path) + tile_cat = Tile_cat(self._tile_cat_path) vignet_cat = self._vignet_cat + # Fail before the per-object loop, whose try/except would otherwise + # drop every object one by one (shapepipe#776). + if self._blend_handling == "uberseg" and tile_cat.seg is None: + raise ValueError( + "BLEND_HANDLING = uberseg needs the tile catalogue's" + + f" SEG_VIGNET column, which {self._tile_cat_path} lacks;" + + " write it at tile detection (SEG_VIGNET = True)." + ) + check_wcs_centroid_offset( self._centroid_source, tile_cat, vignet_cat.gal_vign_cat ) @@ -1218,6 +1258,8 @@ def process(self): id_first = -1 id_last = -1 count_batch = 0 + epoch_cuts = Counter(considered=0, masked_fraction=0, central_veto=0) + n_emptied = 0 saved_batch_cumul = 0 rows = chunk_rows( @@ -1251,10 +1293,18 @@ def process(self): self._bkg_sub, psf_obj, gal_obj, + epoch_central_defect_radius=self._epoch_central_defect_radius, + epoch_masked_fraction_cut=self._epoch_masked_fraction_cut, + defect_fill=self._defect_fill, + epoch_interpolated_defect_radius=( + self._epoch_interpolated_defect_radius + ), ) + epoch_cuts.update(stamp.epoch_cuts) if len(stamp.gals) == 0: n_no_epoch += 1 + n_emptied += stamp.epoch_cuts["considered"] > 0 continue # Per-object RNG, seeded from (ra, dec, ccd) — see @@ -1296,6 +1346,7 @@ def process(self): object_number=obj_id, dilate_neighbour=self._dilate_neighbour, metacal_psf=self._metacal_psf, + defect_fill=self._defect_fill, ) except Exception as ee: self._w_log.info( @@ -1364,6 +1415,13 @@ def process(self): + f" {n_fitted} fitted" ) + self._w_log.info( + "epoch cuts:" + + f" considered={epoch_cuts['considered']}" + + f" masked_fraction={epoch_cuts['masked_fraction']}" + + f" central_veto={epoch_cuts['central_veto']}" + + f" objects_emptied={n_emptied}" + ) log_run_health(self._w_log, count, n_fitted, n_flagged) vignet_cat.close() @@ -1385,11 +1443,75 @@ def prepare_postage_stamps( bkg_sub=True, psf_obj=None, gal_obj=None, + epoch_central_defect_radius=EPOCH_CENTRAL_DEFECT_RADIUS, + epoch_masked_fraction_cut=EPOCH_MASKED_FRACTION_CUT, + defect_fill="noise", + epoch_interpolated_defect_radius=EPOCH_INTERPOLATED_DEFECT_RADIUS, ): - """Prepare the per-object lists of exposures passed to ngmix. + """Gather one object's epoch stamps, dropping epochs its defects spoil. + + @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.epoch_masked_fraction_cut,decision:shape_measurement.defect_fill] epoch-cut-on-defect-mask + An epoch is dropped when more than ``epoch_masked_fraction_cut`` of its + stamp lies in :func:`defect_mask`: flagged, zero-weight and invalid-RMS + pixels, the set that :func:`prepare_ngmix_weights` zero-weights and + fills under every ``blend_handling``. Counting flags alone would keep + epochs whose filled area exceeds the cut. The zero-weight orbit of + interpolated pixels keeps its light and is not counted. The default is + 1/3; 10%, the DES Y3 and Y6 value, is the alternative to test. + + @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.epoch_masked_fraction_cut,decision:shape_measurement.blend_handling] neighbour-markers-are-not-defects + The tile VIGNET holds -1e30 on the footprints of other detections and + beyond the tile's edge (:func:`split_tile_markers`). The + footprint markers form the epoch's neighbour mask (``stamp.neighbours``, + MegaCam-flipped like the epoch), kept apart from its flag stamp, so + neither the masked-fraction cut nor the central veto counts them. Every + epoch shares the tile VIGNET: counting a neighbour within the veto + radius as a defect would drop every epoch of the object. + + @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.epoch_masked_fraction_cut,decision:shape_measurement.defect_fill] off-tile-pixels-are-defects + Beyond the tile's edge the epoch holds real data, the object's own light + cut off by the tile. Those pixels are flagged ``OFF_TILE_FLAG`` (2**10) + in the epoch's flag stamp and so join its defect set under every + ``blend_handling``: zero weight, the defect fill, and both epoch cuts. + On a 51-px stamp an object within about 8.5 px of the tile edge fails + the 1/3 cut. + + Parameters + ---------- + vignet : Vignet + Per-object vignet stores. + obj_id : int + Object ID (SExtractor ``NUMBER``). + i_tile : int + Row of the object in ``tile_cat``. + tile_cat : Tile_cat + Tile catalogue. + bkg_sub : bool, optional + Subtract the background vignet; the default is ``True``. + psf_obj, gal_obj : dict, optional + The object's PSF and galaxy vignet dicts, if already read. + epoch_central_defect_radius : float, optional + Drop an epoch with a defect pixel closer than this many pixels to the + stamp centre (:func:`has_central_defect`); 0 disables the veto. The + default is ``EPOCH_CENTRAL_DEFECT_RADIUS``. + epoch_masked_fraction_cut : float, optional + Drop an epoch with more than this fraction of its stamp in + :func:`defect_mask`. The default is ``EPOCH_MASKED_FRACTION_CUT``. + defect_fill : {"noise", "interpolate"}, optional + The fill :func:`prepare_ngmix_weights` will apply; it sets the veto + radius of each defect pixel (:func:`central_defect_vetoes`). The + default is ``"noise"``. + epoch_interpolated_defect_radius : float, optional + Veto radius for interpolated defect pixels under + ``defect_fill="interpolate"``. The default is + ``EPOCH_INTERPOLATED_DEFECT_RADIUS``. - @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.epoch_masked_fraction_cut] + Returns + ------- + Postage_stamp + The surviving epochs' stamps. """ + # define per-object lists of individual exposures to go into ngmix stamp = Postage_stamp(bkg_sub=bkg_sub) # Read each store's per-object dict ONCE: every sqlitedict access # unpickles the object's whole all-epoch dict, so keeping these out of @@ -1456,19 +1578,30 @@ def prepare_postage_stamps( tile_seg = Ngmix.MegaCamFlip(tile_seg, int(ccd_n)) flag_vign = flag_obj[expccd_name]['VIGNET'] - if tile_vign is not None: - flag_vign[np.where(tile_vign == -1e30)] = 2**10 - v_flag_tmp = flag_vign.ravel() - # remove objects that are more than 1/3 masked - if len(np.where(v_flag_tmp != 0)[0]) / v_flag_tmp.size > 1 / 3.0: - continue - + # Off-tile pixels are defects (off-tile-pixels-are-defects); the + # other -1e30 markers are neighbours (neighbour-markers-are-not-defects). + neighbour, off_tile = split_tile_markers(tile_vign, np.shape(gal_vign)) + flag_vign[off_tile] = OFF_TILE_FLAG weight_vign = weight_obj[expccd_name]['VIGNET'] bkg_rms_vign = ( bkg_rms_obj[expccd_name]['VIGNET'] if bkg_rms_obj is not None else None ) + # Drop the epoch when too much of it would be zero-weighted and + # filled (epoch-cut-on-defect-mask), or when a filled pixel would sit + # near the object (veto-radius-follows-the-fill). + masked = defect_mask(weight_vign, flag_vign, bkg_rms_vign) + stamp.epoch_cuts["considered"] += 1 + if masked.mean() > epoch_masked_fraction_cut: + stamp.epoch_cuts["masked_fraction"] += 1 + continue + if central_defect_vetoes( + masked, epoch_central_defect_radius, defect_fill, + epoch_interpolated_defect_radius, + ): + stamp.epoch_cuts["central_veto"] += 1 + continue # One unpickle per exposure (all CCDs), reused across this object's # epochs; the cache is per-call, so bounded by the object's exposures @@ -1502,6 +1635,7 @@ def prepare_postage_stamps( stamp.psfs.append(psf_obj[expccd_name]['VIGNET']) stamp.weights.append(weight_vign_scaled) stamp.flags.append(flag_vign) + stamp.neighbours.append(neighbour) stamp.bkg_rms.append(bkg_rms_vign_scaled) if tile_seg is not None: stamp.segs.append(tile_seg) @@ -1511,9 +1645,7 @@ def prepare_postage_stamps( # make_ngmix_observation), which raises if it is missing; the "hsm" # path ignores it, so read it leniently rather than coupling hsm to a # field it never uses. - stamp.offsets.append( - vignet.gal_vign_cat[str(obj_id)][expccd_name].get('OFFSET') - ) + stamp.offsets.append(gal_obj[expccd_name].get('OFFSET')) stamp.ra.append(tile_cat.ra[i_tile]) stamp.dec.append(tile_cat.dec[i_tile]) # CCD of the first surviving epoch — Fabian's coord_list[0] convention @@ -1524,6 +1656,51 @@ def prepare_postage_stamps( return stamp +def split_tile_markers(tile_vign, shape): + """Split the tile VIGNET's -1e30 markers into neighbour and off-tile. + + @sc [decision:shape_measurement.blend_handling,decision:shape_measurement.defect_fill] off-tile-is-marked-border-rows-and-columns + The tile VIGNET holds -1e30 on the footprints of other detections and on + stamp pixels beyond the tile's edge, as SExtractor writes it. A stamp + clipped by the tile's + rectangle loses whole rows and whole columns from its border, so the + off-tile pixels are the union of the runs of entirely -1e30 rows and + columns that start at a stamp border. The remaining markers are + neighbour pixels: a footprint touching the stamp border, and a footprint + that completes an interior row or column beside an off-tile band, stay + neighbours. + + Parameters + ---------- + tile_vign : numpy.ndarray or None + Tile VIGNET stamp, oriented like the epoch; ``None`` marks nothing. + shape : tuple of int + Stamp shape, used when ``tile_vign`` is ``None``. + + Returns + ------- + numpy.ndarray of bool + Neighbour pixels. + numpy.ndarray of bool + Off-tile pixels. + """ + if tile_vign is None: + return np.zeros(shape, dtype=bool), np.zeros(shape, dtype=bool) + marker = tile_vign == -1e30 + + def border_runs(full): + # Lines in the unbroken run of marked lines from either border. + lead = np.logical_and.accumulate(full) + trail = np.logical_and.accumulate(full[::-1])[::-1] + return lead | trail + + off_tile = ( + border_runs(marker.all(axis=1))[:, None] + | border_runs(marker.all(axis=0))[None, :] + ) + return marker & ~off_tile, off_tile + + def background_subtract(gal,bkg): """background subtraction. @@ -1767,9 +1944,9 @@ def uberseg_weight(weight, seg, object_number, dilate_neighbour=0): Because the partition is by distance to the nearest footprint, the pixels surviving around a compact central object form a single connected, roughly circular core; the "circularisation" is emergent geometry, not a - separate aperture. Unlike the noise-fill treatment, masked pixels are - handed to ngmix as a hard mask (weight = 0), never replaced by a noise - realisation. + separate aperture. Only the weight changes: neighbour-side pixels are + handed to ngmix as a hard mask (weight = 0) and keep their image values, + so metacal shears the neighbour's light along with the target's. Parameters ---------- @@ -1831,33 +2008,260 @@ def uberseg_weight(weight, seg, object_number, dilate_neighbour=0): return weight + +def defect_mask(weight, flag, bkg_rms=None): + """Defect pixels of one epoch stamp. + + @sc [decision:shape_measurement.defect_fill,label:physics] defect-set-unsymmetrized + A defect is a pixel with zero exposure weight, a nonzero exposure flag, + or (when a background RMS map is given) a non-finite or non-positive RMS. + This one set is zero-weighted and filled by :func:`prepare_ngmix_weights` + under every ``blend_handling`` and counted by the epoch cuts in + :func:`prepare_postage_stamps`. The tile VIGNET's neighbour markers are not + in it (neighbour-markers-are-not-defects); off-tile pixels are, as flag + ``OFF_TILE_FLAG`` (off-tile-pixels-are-defects). It is not ORed with its + rotations. For defects the central-defect veto keeps + (:func:`has_central_defect`), the unsymmetrized fill leaves + |c| <= 3e-4 per affected epoch on galaxies with half-light radius 0.3" + and 0.5" through a 0.7" PSF, round or elliptical. Symmetrizing would + quadruple the filled area near the object, and with it m: a 3-px bleed at + 10 px on the 0.5" galaxy gives m11 = -2.7% four-fold against -0.64% + unsymmetrized. Through a PSF with ellipticity (0.05, 0.02) it would not + cancel c either: c1 = 7.3e-4 four-fold against 0.7e-4. It would also + turn a 5-px edge band (9.8% of the stamp), which biases nothing, into a + 35.4% frame that fails the masked-fraction cut. + + Parameters + ---------- + weight : numpy.ndarray + Exposure weight stamp. + flag : numpy.ndarray + Exposure flag stamp. + bkg_rms : numpy.ndarray, optional + Background RMS stamp. + + Returns + ------- + numpy.ndarray of bool + ``True`` on defect pixels. + """ + defect = (weight == 0) | (flag != 0) + if bkg_rms is not None: + defect |= ~(np.isfinite(bkg_rms) & (bkg_rms > 0)) + return defect + + +def has_central_defect(defect, radius): + """Whether a defect pixel lies closer than ``radius`` to the stamp centre. + + @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.defect_fill] epoch-central-defect-veto + An epoch is dropped when any pixel of :func:`defect_mask` lies closer + than ``radius`` pixels to the stamp centre. A noise-filled hole in the + object's light is sheared by metacal but not by the sky, so the metacal + response is wrong for that epoch, and a one-sided hole adds an additive + term. The default radius (``EPOCH_CENTRAL_DEFECT_RADIUS``, 10 px or + 1.9") is the smallest at which columns, 3-px bleeds, single pixels and + edge bands give |m11|, |m22| < 1% and |c1|, |c2| < 5e-4 (full response + matrix) on galaxies with half-light radius 0.3" and 0.5" through a 0.7" + PSF. The bias is anisotropic: a column at 10 px on the 0.5" galaxy gives + m11 = -0.60%, m22 = +0.01% and c1 = -1.4e-4; at 9 px, m11 = -1.75% and + c1 = -1.1e-3. Through a PSF with ellipticity (0.05, 0.02), wide defects + at 10 px on the 0.5" galaxy sit at the bound: a 3-px bleed gives + m11 = -0.98% +/- 0.04% and the widest edge band the veto keeps (16 px) + -0.94% +/- 0.06%; at 11 px the bleed gives -0.24%. Larger galaxies need + more: at hlr 0.7" through a 0.9" PSF a 3-px bleed at 10 px gives + m11 = -6.8% and c1 = -4.9e-3, and passes from 14 px. The distance is + measured from the stamp centre, where the extractor places the object to + within half a pixel. The veto reads only the defect mask, never the + object's pixels, so it selects on nothing that responds to shear. + Radius 0 disables it. Guarded by ``tests/science/test_defect_veto.py``. + + Parameters + ---------- + defect : numpy.ndarray of bool + Defect mask of one epoch stamp (:func:`defect_mask`). + radius : float + Veto radius in pixels. + + Returns + ------- + bool + ``True`` if the epoch should be dropped. + """ + rows, cols = np.nonzero(defect) + centre_row = (defect.shape[0] - 1) / 2 + centre_col = (defect.shape[1] - 1) / 2 + distance = np.hypot(rows - centre_row, cols - centre_col) + return bool(np.any(distance < radius)) + + +def central_defect_vetoes( + defect, radius, defect_fill="noise", + interpolated_radius=EPOCH_INTERPOLATED_DEFECT_RADIUS, +): + """Whether the central-defect veto drops an epoch under ``defect_fill``. + + @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.defect_fill] veto-radius-follows-the-fill + Each defect pixel is vetoed (:func:`has_central_defect`) at the radius + calibrated for the operator that fills it. Under ``"noise"`` every + defect keeps ``radius`` (``EPOCH_CENTRAL_DEFECT_RADIUS``). Under + ``"interpolate"`` the pixels :func:`interpolable_defects` selects are + interpolated and vetoed at ``interpolated_radius`` + (``EPOCH_INTERPOLATED_DEFECT_RADIUS``); the others are noise-filled and + keep ``radius``. 7 px is the smallest interpolated radius at which + columns, full and finite 3-px bleeds and single pixels give + |m11|, |m22| < 1% and |c1|, |c2| < 5e-4 (full response matrix) on + galaxies with half-light radius 0.3" and 0.5" through a 0.7" PSF, round + or with ellipticity (0.05, 0.02), and on a 0.7" galaxy through a 0.9" + PSF. There the worst case is c1 = 3.6e-4 (3-px bleed, 0.7" galaxy, + elliptical PSF); on the 0.3" and 0.5" galaxies |c| <= 5e-5 and + |m| <= 0.21%. At 6 px those two still pass (|c| <= 3.1e-4), but a 3-px + bleed on the 0.7" galaxy gives c1 = 7.5e-4; at 5 px a 3-px bleed on the + 0.5" galaxy gives 9.5e-4. The radius does not scale with galaxy size: + interpolation needs one pixel more for the 0.7" galaxy (noise fill needs + four, 10 to 14 px), and a size-dependent veto would select on measured + size, which responds to shear. Guarded by + ``tests/science/test_defect_interpolation.py``. + + Parameters + ---------- + defect : numpy.ndarray of bool + Defect mask of one epoch stamp (:func:`defect_mask`). + radius : float + Veto radius in pixels for noise-filled defect pixels. + defect_fill : {"noise", "interpolate"}, optional + The fill the epoch will get; the default is ``"noise"``. + interpolated_radius : float, optional + Veto radius in pixels for interpolated defect pixels; the default is + ``EPOCH_INTERPOLATED_DEFECT_RADIUS``. + + Returns + ------- + bool + ``True`` if the epoch should be dropped. + """ + if defect_fill == "interpolate": + interpolated = interpolable_defects(defect) + return has_central_defect( + interpolated, interpolated_radius + ) or has_central_defect(defect & ~interpolated, radius) + return has_central_defect(defect, radius) + + +def fill_defects(image, defect, noise): + """Replace an image's defect pixels by a noise realisation. + + @sc [decision:shape_measurement.defect_fill,label:physics] defect-noise-fill + Defect pixels are replaced by ``noise``, an independent realisation at + the per-pixel background RMS. Metacal deconvolves, shears and reconvolves + the whole image without reading the weights, so a raw defect value (bad + column, bleed, cosmic ray) would leak into the fit. The noise fill leaves + a hole in the object's light that metacal shears and the sky does not, + which biases m for an epoch with a defect near the object; the + central-defect veto drops those epochs (epoch-central-defect-veto). + + Parameters + ---------- + image : numpy.ndarray + Stamp image. + defect : numpy.ndarray of bool + Pixels to replace: the defects (:func:`defect_mask`), and under + noisefill the marked neighbour pixels. + noise : numpy.ndarray + Noise realisation on the stamp grid. + + Returns + ------- + numpy.ndarray + ``image`` with ``defect`` pixels taken from ``noise``, in the dtype + of ``image``. + """ + return np.where(defect, noise, image).astype(image.dtype, copy=False) + + def prepare_ngmix_weights( gal, weight, flag, rng, bkg_rms=None, blend_handling="noisefill", seg=None, object_number=None, - dilate_neighbour=0, + dilate_neighbour=0, defect_fill="noise", neighbour=None, ): - """bookkeeping for ngmix weights. runs on a single galaxy and epoch - pixel scale and galaxy guess - TO DO: decide if we want galaxy guess stuff + """Build one epoch's image, weight map and noise image for ngmix. + + Defect pixels (:func:`defect_mask`: flagged, zero-weight or invalid-RMS + pixels) get weight 0 and are filled under either ``blend_handling``: + by an independent noise realisation at their background RMS + (:func:`fill_defects`), or under ``defect_fill="interpolate"``, for + short defect runs, by an interpolant of the clean pixels around them. + ``blend_handling`` decides only the neighbour pixels. + + @sc [decision:shape_measurement.defect_fill,decision:shape_measurement.blend_handling] defects-filled-whatever-the-blend-handling + Every pixel of ``defect_mask`` is zero-weighted and filled whatever + ``blend_handling`` is, by the same operator. ``blend_handling`` acts on + the neighbour pixels only, and the noise image covers the whole stamp. + + @sc [decision:shape_measurement.defect_fill,decision:shape_measurement.blend_handling] interpolation-support-is-the-kept-image + Under ``defect_fill="interpolate"`` the interpolant is supported on the + pixels whose light the image keeps: never a defect, and under + ``"noisefill"`` never a marked neighbour pixel, whose light is replaced + by noise; a Clough-Tocher interpolant next to a bright neighbour would + otherwise build the defect's value from light the image no longer + contains. Under ``"uberseg"`` the neighbour light stays in the image and + supports the interpolant. + + @sc [decision:shape_measurement.blend_handling] noisefill-fills-markers + Under ``"noisefill"`` the pixels of ``neighbour``, the tile VIGNET's + -1e30 neighbour markers, get weight 0 and are replaced by the same noise + realisation as the noise-filled defects, so no marked neighbour light + reaches metacal. Under the default noise fill, the image, weight map and + noise image are those the marked pixels would get as flagged defects; + only the epoch cuts treat them differently, by not counting them + (neighbour-markers-are-not-defects). Unmarked neighbour light stays + raw and weighted. + + @sc [decision:shape_measurement.blend_handling] uberseg-ignores-markers + Under ``"uberseg"`` the neighbour markers are ignored: :func:`uberseg_weight` + zeroes the weight of pixels nearer a neighbour's segmentation footprint + than the target's and leaves their image values raw, because that light + is real sky that metacal shears along with the target; filling it with + noise would cut the target along an unsheared edge. + + @sc [decision:shape_measurement.defect_fill,label:physics] interpolated-defect-weight-orbit + Under ``defect_fill="interpolate"`` the short defect runs + (:func:`interpolable_defects`) take the interpolant of the clean pixels + around them (:func:`interpolate_defects`) in the image and in the noise + image alike, and the other defects are noise-filled. Interpolation + restores the object's light, which leaves the hole in the likelihood: + ngmix fits a Gaussian to a non-Gaussian profile, and a one-sided + zero-weight hole pulls that fit. Once interpolated, a column 8 px from a + galaxy with half-light radius 0.5" through a 0.7" PSF still gives + c1 = -1.3e-3. The weight is therefore also zero on the quarter-turn + orbit of the interpolated pixels, whose light stays, which cancels that + spin-2 term: c1 = +2e-6 for the same column. Only the weights are + symmetrized, not the fill: interpolating the orbit as well would trade + true light for interpolated light over four times the area, raising m11 + for a 3-px bleed 6 px from the same galaxy from +0.19% to +0.89%. Parameters ---------- gal : numpy.ndarray + Background-subtracted galaxy stamp. weight : numpy.ndarray + Exposure weight stamp; zero marks a defect. flag : numpy.ndarray + Exposure flag stamp; nonzero marks a defect. rng : numpy.random.RandomState Random state for the noise realisations (seeded per object; see :func:`position_seed`). bkg_rms : numpy.ndarray, optional - Per-pixel background RMS map. If supplied, unmasked pixels use - ``1 / bkg_rms**2`` as the ngmix inverse variance. + Per-pixel background RMS map. If supplied, clean pixels use + ``1 / bkg_rms**2`` as the ngmix inverse variance, and non-finite or + non-positive values mark defects. Otherwise every clean pixel gets + ``1 / sigma_mad(gal)**2``. blend_handling : {"noisefill", "uberseg"}, optional - How to treat pixels shared with a neighbour. ``"noisefill"`` (default) - replaces flagged pixels with a noise realisation and keeps their - inverse-variance weight — the historical behaviour. ``"uberseg"`` - instead hard-masks (weight = 0) every pixel closer to a neighbour's - segmentation footprint than to the central object's, leaving the - image untouched (see :func:`uberseg_weight`). + Neighbour treatment. ``"noisefill"`` (default) zero-weights and + noise-fills the ``neighbour`` pixels. ``"uberseg"`` ignores + ``neighbour``, zeroes the weight of every pixel closer to a + neighbour's segmentation footprint than to the central object's and + keeps its raw image value (see :func:`uberseg_weight`). seg : numpy.ndarray, optional Segmentation map on the stamp grid (object NUMBERs). Required for ``blend_handling="uberseg"``; ignored otherwise. @@ -1867,49 +2271,90 @@ def prepare_ngmix_weights( dilate_neighbour : int, optional Neighbour-mask dilation iterations, passed to :func:`uberseg_weight` under ``blend_handling="uberseg"``; ignored otherwise. + defect_fill : {"noise", "interpolate"}, optional + ``"noise"`` (default) noise-fills every defect; ``"interpolate"`` + interpolates the short defect runs and noise-fills the rest. + neighbour : numpy.ndarray of bool, optional + The epoch's neighbour mask (``Postage_stamp.neighbours``); read only + under ``blend_handling="noisefill"``. ``None`` marks no pixel. Returns ------- numpy.ndarray - Galaxy image. For ``"noisefill"`` masked pixels are replaced by noise; - for ``"uberseg"`` the image is returned untouched. + Galaxy image with defect pixels, and under noisefill marked + neighbour pixels, filled. numpy.ndarray - Variance map for NGMIX. + Inverse-variance weight map for ngmix. numpy.ndarray - Noise image. + Noise image: an independent realisation over the whole stamp, for + metacal's ``fixnoise``, interpolated where the galaxy image is. + + Raises + ------ + ValueError + If ``blend_handling`` or ``defect_fill`` is unknown, or + ``"uberseg"`` lacks ``seg`` or ``object_number``. @sc [decision:masking.pixel_mask_source,decision:shape_measurement.blend_handling,decision:shape_measurement.defect_fill,decision:shape_measurement.galaxy_pixel_weights] """ + if defect_fill not in DEFECT_FILLS: + raise ValueError( + f"Unknown DEFECT_FILL '{defect_fill}'; expected one of" + + f" {DEFECT_FILLS}" + ) if blend_handling not in BLEND_HANDLINGS: raise ValueError( f"Unknown blend_handling '{blend_handling}'; expected one of" + f" {BLEND_HANDLINGS}" ) + if blend_handling == "uberseg" and (seg is None or object_number is None): + raise ValueError( + "blend_handling='uberseg' requires a segmentation map and the" + + " central object_number; none reached prepare_ngmix_weights." + + " The tile catalogue's SEG_VIGNET column carries the map (see" + + " CosmoStat/shapepipe#776)." + ) - mask = np.copy(weight) != 0 - mask[flag != 0] = False + defect = defect_mask(weight, flag, bkg_rms) + # Marked neighbour pixels that noisefill removes (noisefill-fills-markers). + removed_neighbour = ( + np.asarray(neighbour, dtype=bool) + if blend_handling == "noisefill" and neighbour is not None + else np.zeros_like(defect) + ) + clean = ~(defect | removed_neighbour) + interpolated = ( + interpolable_defects(defect) + if defect_fill == "interpolate" + else np.zeros_like(defect) + ) + # The quarter-turn orbit of the interpolated pixels carries no weight + # (interpolated-defect-weight-orbit). + weighted = ( + clean & ~fourfold(interpolated) if interpolated.any() else clean + ) if bkg_rms is None: sig_noise = sigma_mad(gal) # Guard the degenerate constant stamp (sigma_mad == 0): 0 * inf # would otherwise put NaN in a fully-masked weight map. weight_map = ( - mask.astype(float) / sig_noise ** 2 + weighted.astype(float) / sig_noise ** 2 if sig_noise > 0 else np.zeros_like(gal, dtype=float) ) else: - valid_rms = np.isfinite(bkg_rms) & (bkg_rms > 0) - mask &= valid_rms weight_map = np.zeros_like(gal, dtype=float) - weight_map[mask] = 1.0 / bkg_rms[mask] ** 2 + weight_map[weighted] = 1.0 / bkg_rms[weighted] ** 2 # Per-pixel noise sigma for the realisations below: metacal's # fixnoise bookkeeping (1/w + 1/w_noise) assumes the noise image # is a faithful realisation of the per-pixel variance the weights # claim; a scalar sigma there mis-reports errors and erodes the # inverse-variance advantage whenever the RMS map actually varies. + # Pixels without a valid RMS take the median over clean pixels. + valid_rms = np.isfinite(bkg_rms) & (bkg_rms > 0) sig_noise = ( - np.where(valid_rms, bkg_rms, np.median(bkg_rms[mask])) - if mask.any() + np.where(valid_rms, bkg_rms, np.median(bkg_rms[clean])) + if clean.any() else sigma_mad(gal) ) @@ -1922,34 +2367,34 @@ def prepare_ngmix_weights( noise_img = rng.standard_normal(gal.shape) * sig_noise noise_img_gal = rng.standard_normal(gal.shape) * sig_noise + gal_filled = fill_defects(gal, ~clean, noise_img_gal) + if interpolated.any(): + # One operator for the image and the noise image, supported on the + # pixels the image keeps (interpolation-support-is-the-kept-image); a + # pixel whose support is degenerate, or that noisefill removes as a + # neighbour, keeps its noise fill. + filled = interpolate_defects([gal, noise_img], ~clean, interpolated) + done = ( + interpolated + & ~removed_neighbour + & np.all(np.isfinite(filled), axis=0) + ) + gal_filled[done] = filled[0][done] + noise_img = np.where(done, filled[1], noise_img) - gal_masked = np.copy(gal) if blend_handling == "uberseg": - # Hard-mask neighbour-side pixels (weight -> 0) from the segmentation - # geometry; the image is left untouched (the masked pixels carry no - # weight, so ngmix ignores them in the likelihood). Bad/flagged - # pixels already sit at weight 0 from the mask above. - if seg is None or object_number is None: - raise ValueError( - "blend_handling='uberseg' requires a segmentation map and the" - + " central object_number; none reached prepare_ngmix_weights." - + " Set SEG_VIGNET_PATH on the ngmix run (see" - + " CosmoStat/shapepipe#776)." - ) weight_map = uberseg_weight( weight_map, seg, object_number, dilate_neighbour=dilate_neighbour ) - elif (~mask).any(): - # noisefill (default): replace masked pixels with a noise realisation. - gal_masked[~mask] = noise_img_gal[~mask] - return gal_masked, weight_map, noise_img + return gal_filled, weight_map, noise_img + def make_ngmix_observation( gal, weight, flag, psf, wcs, rng, bkg_rms=None, centroid_source="wcs", offset=None, blend_handling="noisefill", seg=None, object_number=None, - dilate_neighbour=0, + dilate_neighbour=0, defect_fill="noise", neighbour=None, ): """Build an ngmix Observation for a single galaxy epoch. @@ -1993,7 +2438,8 @@ def make_ngmix_observation( ``"hsm"``). blend_handling : {"noisefill", "uberseg"}, optional Neighbour treatment passed through to :func:`prepare_ngmix_weights`; - the default ``"noisefill"`` is the historical behaviour. + the default ``"noisefill"`` zero-weights and noise-fills the + ``neighbour`` pixels. seg : numpy.ndarray, optional Segmentation map on the stamp grid. Required for ``blend_handling="uberseg"`` (ignored otherwise). @@ -2003,6 +2449,11 @@ def make_ngmix_observation( dilate_neighbour : int, optional Neighbour-mask dilation iterations passed through to :func:`prepare_ngmix_weights` under ``blend_handling="uberseg"``. + defect_fill : {"noise", "interpolate"}, optional + Defect fill passed through to :func:`prepare_ngmix_weights`; the + default is ``"noise"``. + neighbour : numpy.ndarray of bool, optional + Neighbour mask passed through to :func:`prepare_ngmix_weights`. Returns ------- @@ -2031,16 +2482,18 @@ def make_ngmix_observation( gal_masked, weight_map, noise_img = prepare_ngmix_weights( gal, weight, flag, rng, bkg_rms=bkg_rms, blend_handling=blend_handling, seg=seg, object_number=object_number, - dilate_neighbour=dilate_neighbour, + dilate_neighbour=dilate_neighbour, defect_fill=defect_fill, + neighbour=neighbour, ) if centroid_source == "hsm": # Re-center the Jacobian on the HSM adaptive-moment centroid (pixel - # offset from the stamp center); fall back to the stamp center if - # HSM fails. + # offset from the stamp center), measured on the filled image so no + # raw defect value pulls it; fall back to the stamp center if HSM + # fails. try: _hsm = galsim.hsm.FindAdaptiveMom( - galsim.Image(gal, scale=1.0), strict=False + galsim.Image(gal_masked, scale=1.0), strict=False ) if _hsm.error_message != "": raise galsim.hsm.GalSimHSMError(_hsm.error_message) @@ -2261,7 +2714,7 @@ def make_runners(prior, flux_guess, rng): def do_ngmix_metacal( stamp, prior, flux_guess, rng, centroid_source="wcs", blend_handling="noisefill", object_number=None, dilate_neighbour=0, - metacal_psf="fitgauss", + metacal_psf="fitgauss", defect_fill="noise", ): """Do Ngmix Metacal. @@ -2286,9 +2739,9 @@ def do_ngmix_metacal( stamp pixels — see that function. blend_handling : {"noisefill", "uberseg"}, optional Neighbour treatment passed through to - :func:`make_ngmix_observation`; the default ``"noisefill"`` is the - historical behaviour. ``"uberseg"`` consumes ``stamp.segs`` and - ``object_number``. + :func:`make_ngmix_observation`; the default ``"noisefill"`` + zero-weights and noise-fills the pixels of ``stamp.neighbours``. + ``"uberseg"`` consumes ``stamp.segs`` and ``object_number``. object_number : int, optional Central object's segmentation label — its SExtractor ``NUMBER`` (``obj_id``), authoritative because seg labels are the NUMBERs of the @@ -2308,6 +2761,9 @@ def do_ngmix_metacal( round PSF that metacal reconvolves with after shearing, so it moves the metacal *response* (and therefore the recovered shear) but never reaches the deconvolution, which is by the PSF image. + defect_fill : {"noise", "interpolate"}, optional + Defect fill passed through to :func:`make_ngmix_observation`; the + default is ``"noise"``. Returns ------- @@ -2341,6 +2797,10 @@ def do_ngmix_metacal( seg=stamp.segs[n_e] if n_e < len(stamp.segs) else None, object_number=object_number, dilate_neighbour=dilate_neighbour, + defect_fill=defect_fill, + neighbour=( + stamp.neighbours[n_e] if n_e < len(stamp.neighbours) else None + ), ) gal_obs_list.append(gal_obs) diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 6d15905fa..86ec796de 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -11,7 +11,13 @@ from sqlitedict import SqliteDict from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.ngmix_package.ngmix import Ngmix, write_empty_tile_output +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_INTERPOLATED_DEFECT_RADIUS, + EPOCH_MASKED_FRACTION_CUT, + Ngmix, + write_empty_tile_output, +) @module_runner( @@ -86,24 +92,6 @@ def ngmix_runner( else: input_file_list = input_file_list[:wcs_idx] - # SEG_VIGNET_PATH (optional): coadd-frame SExtractor segmentation vignets - # (a CLASSIC-mode vignetmaker output), row-aligned to the tile catalogue. - # Required for BLEND_HANDLING = uberseg; when set, the file must exist for - # every tile (missing file -> error). Read on Tile_cat, not via Vignet, so - # it is threaded to Ngmix as its own argument rather than into - # input_file_list. - if config.has_option(module_config_sec, "SEG_VIGNET_PATH"): - seg_vignet_path = config.getexpanded( - module_config_sec, - "SEG_VIGNET_PATH", - ).format(file_number_string=file_number_string) - if not os.path.exists(seg_vignet_path): - raise FileNotFoundError( - f"Segmentation vignet file not found: {seg_vignet_path}" - ) - else: - seg_vignet_path = None - # Batch save option if config.has_option(module_config_sec, "SAVE_BATCH"): save_batch = config.getint( @@ -130,12 +118,15 @@ def ngmix_runner( else: centroid_source = "wcs" - # Neighbour treatment: "noisefill" (default, historical) replaces a - # neighbour's pixels with a noise realisation; "uberseg" hard-masks - # (weight -> 0) every pixel closer to a neighbour than to the central - # object, from the segmentation map. See the ngmix module docstrings. + # Neighbour treatment: "noisefill" (default) zero-weights and noise-fills + # the pixels marked -1e30 in the tile VIGNET (other detections' + # footprints); "uberseg" ignores those markers, zeroes the weight of every + # pixel closer to a neighbour than to the central object, from the + # segmentation map, and leaves its image raw. Defect pixels (flagged, + # zero-weight, invalid-RMS or off-tile) are zero-weighted and filled + # under both (DEFECT_FILL below); see prepare_ngmix_weights. if config.has_option(module_config_sec, "BLEND_HANDLING"): - blend_handling = config.get(module_config_sec, "BLEND_HANDLING") + blend_handling = config.getexpanded(module_config_sec, "BLEND_HANDLING") else: blend_handling = "noisefill" @@ -147,6 +138,47 @@ def ngmix_runner( else: dilate_neighbour = 1 + # EPOCH_CENTRAL_DEFECT_RADIUS (optional, pixels): drop an epoch when a + # defect pixel (flagged, zero-weight, invalid-RMS or off-tile; not a + # neighbour marker) lies closer than this to the stamp centre; 0 disables. + if config.has_option(module_config_sec, "EPOCH_CENTRAL_DEFECT_RADIUS"): + epoch_central_defect_radius = config.getfloat( + module_config_sec, "EPOCH_CENTRAL_DEFECT_RADIUS" + ) + else: + epoch_central_defect_radius = EPOCH_CENTRAL_DEFECT_RADIUS + + # EPOCH_MASKED_FRACTION_CUT (optional): drop an epoch when more than this + # fraction of its stamp is defects (flagged, zero-weight, invalid-RMS or + # off-tile; neighbour markers are not counted). + if config.has_option(module_config_sec, "EPOCH_MASKED_FRACTION_CUT"): + epoch_masked_fraction_cut = config.getfloat( + module_config_sec, "EPOCH_MASKED_FRACTION_CUT" + ) + else: + epoch_masked_fraction_cut = EPOCH_MASKED_FRACTION_CUT + + # DEFECT_FILL (optional): "noise" (default) noise-fills every defect; + # "interpolate" interpolates short defect runs (at most 3 px along a row + # or column) from the clean pixels around them and noise-fills the rest. + if config.has_option(module_config_sec, "DEFECT_FILL"): + defect_fill = config.get(module_config_sec, "DEFECT_FILL") + else: + defect_fill = "noise" + + # EPOCH_INTERPOLATED_DEFECT_RADIUS (optional, pixels): under + # DEFECT_FILL = interpolate, drop an epoch when an interpolated defect + # pixel lies closer than this to the stamp centre; noise-filled defect + # pixels keep EPOCH_CENTRAL_DEFECT_RADIUS. + if config.has_option( + module_config_sec, "EPOCH_INTERPOLATED_DEFECT_RADIUS" + ): + epoch_interpolated_defect_radius = config.getfloat( + module_config_sec, "EPOCH_INTERPOLATED_DEFECT_RADIUS" + ) + else: + epoch_interpolated_defect_radius = EPOCH_INTERPOLATED_DEFECT_RADIUS + # Check PSF vignets first: if all are empty dicts {}, the exposures for this # tile are absent from the PSF dictionary and no shape measurement is possible. # This check must come before reading image vignets to avoid a C-level malloc @@ -204,9 +236,12 @@ def ngmix_runner( bkg_sub=bkg_sub, centroid_source=centroid_source, blend_handling=blend_handling, - seg_cat_path=seg_vignet_path, dilate_neighbour=dilate_neighbour, metacal_psf=metacal_psf, + epoch_central_defect_radius=epoch_central_defect_radius, + epoch_masked_fraction_cut=epoch_masked_fraction_cut, + defect_fill=defect_fill, + epoch_interpolated_defect_radius=epoch_interpolated_defect_radius, ) # Process ngmix shape measurement and metacalibration diff --git a/src/shapepipe/modules/sextractor_package/match_catalogue.py b/src/shapepipe/modules/sextractor_package/match_catalogue.py index 68bca78f7..6c59d31e7 100644 --- a/src/shapepipe/modules/sextractor_package/match_catalogue.py +++ b/src/shapepipe/modules/sextractor_package/match_catalogue.py @@ -23,9 +23,13 @@ # The SEG_VIGNET label of a footprint whose SExtractor row has no partner in # the external catalogue, and so leaves the catalogue. Negative, so it never -# collides with a NUMBER; UberSeg only asks "self or not self". +# collides with a NUMBER (match_catalogue requires positive ones); UberSeg only +# asks "self or not self". UNMATCHED_LABEL = -1 +# The largest NUMBER the int32 NUMBER and SEG_VIGNET columns hold. +MAX_NUMBER = np.iinfo(np.int32).max + def mutual_nearest(x_a, y_a, x_b, y_b, radius): """One-to-one pairs of mutual nearest neighbours closer than ``radius``. @@ -100,9 +104,12 @@ def match_catalogue(cat_path, ext_cat_path, radius=1.0, min_fraction=0.98, external catalogue (``X_IMAGE``, ``Y_IMAGE``, same image grid) within ``radius`` pixels. Paired rows take the external ``NUMBER``; unpaired rows leave the catalogue. A ``SEG_VIGNET`` column, when present, is - relabelled to the new numbering, with the footprints of rows that left - marked ``UNMATCHED_LABEL``. The catalogue is rewritten in place; every - other HDU and column is kept. + relabelled through the whole map from old to new numbers, with the + footprints of rows that left marked ``UNMATCHED_LABEL``: as the external + numbers are unique integers in ``[1, MAX_NUMBER]``, each row's own + footprint carries its new ``NUMBER`` and no other footprint can, + whatever the two numberings share. + The catalogue is rewritten in place; every other HDU and column is kept. Both catalogues come from the same pixels, so pairs agree to ~1e-4 pixel, and on eight DR6 tiles at least 99.2% of each side pairs; the @@ -144,13 +151,25 @@ def match_catalogue(cat_path, ext_cat_path, radius=1.0, min_fraction=0.98, Raises ------ ValueError - If more than ``tolerated_unpaired`` rows of either side, and more - than ``1 - min_fraction`` of it, have no partner + If the external ``NUMBER`` repeats or is not an integer in + ``[1, MAX_NUMBER]``, or if more + than ``tolerated_unpaired`` rows of either side, and more than + ``1 - min_fraction`` of it, have no partner @sc [decision:detection.tile_detection] """ ext = asc.read(ext_cat_path, format="sextractor", include_names=["NUMBER", "X_IMAGE", "Y_IMAGE"]) + ext_number = np.asarray(ext["NUMBER"]) + if (not np.issubdtype(ext_number.dtype, np.integer) + or (ext_number <= 0).any() or (ext_number > MAX_NUMBER).any() + or len(np.unique(ext_number)) < len(ext)): + raise ValueError( + f"{ext_cat_path} has a NUMBER that is not an integer in" + + f" [1, {MAX_NUMBER}] or that repeats; the join needs unique" + + " numbers that fit the int32 NUMBER and SEG_VIGNET columns, so" + + " that no relabelled footprint takes another object's number." + ) with fits.open(cat_path) as hdul: hdus = [hdu.copy() for hdu in hdul] objects = next(h for h in hdus if h.name == "LDAC_OBJECTS") @@ -194,7 +213,7 @@ def too_many_unpaired(n_total): old_number = np.asarray(data["NUMBER"]) new_number = np.full(len(data), UNMATCHED_LABEL, np.int64) - new_number[i_sex] = np.asarray(ext["NUMBER"])[i_ext] + new_number[i_sex] = ext_number[i_ext] columns = [] for col in objects.columns: diff --git a/src/shapepipe/modules/sextractor_package/sextractor_script.py b/src/shapepipe/modules/sextractor_package/sextractor_script.py index 002bf14f7..8d05d7edb 100644 --- a/src/shapepipe/modules/sextractor_package/sextractor_script.py +++ b/src/shapepipe/modules/sextractor_package/sextractor_script.py @@ -16,6 +16,190 @@ from shapepipe.pipeline.sqlite_store import read_sqlitedict +def cut_stamps(array, col, row, stamp_size, fill): + """Cut one square stamp per object from a 2-D array. + + Each stamp is ``stamp_size`` pixels on a side with the 0-based pixel + (``row``, ``col``) at index ``stamp_size // 2`` on both axes; pixels off + the array take ``fill``. This is the window SExtractor cuts ``VIGNET`` + with, so stamps cut at the same centres from arrays on one pixel grid are + registered pixel for pixel. + + Parameters + ---------- + array : numpy.ndarray + 2-D array, shape ``(ny, nx)`` + col, row : array_like of int + 0-based column and row of each stamp's centre pixel + stamp_size : int + Side length of the stamps + fill : scalar + Value of the stamp pixels off the array + + Returns + ------- + numpy.ndarray + Stamps, shape ``(n_obj, stamp_size, stamp_size)``, dtype of ``array`` + + """ + ny, nx = array.shape + half = stamp_size // 2 + stamps = np.full((len(col), stamp_size, stamp_size), fill, array.dtype) + for i, (xi, yi) in enumerate(zip(col, row)): + x0, y0 = xi - half, yi - half + xc0, xc1 = max(0, x0), min(nx, x0 + stamp_size) + yc0, yc1 = max(0, y0), min(ny, y0 + stamp_size) + if xc0 >= xc1 or yc0 >= yc1: + continue + stamps[i, yc0 - y0:yc1 - y0, xc0 - x0:xc1 - x0] = ( + array[yc0:yc1, xc0:xc1] + ) + return stamps + + +def seg_vignet_column(seg_vignets): + """The ``SEG_VIGNET`` LDAC column: one int32 seg stamp per object. + + Parameters + ---------- + seg_vignets : numpy.ndarray + Segmentation stamps, shape ``(n_obj, stamp_size, stamp_size)``, cut + on ``VIGNET``'s grid (:func:`cut_stamps`) + + Returns + ------- + astropy.io.fits.Column + The column, laid out as ``VIGNET`` is + + """ + n_obj, ny, nx = seg_vignets.shape + return fits.Column( + name="SEG_VIGNET", + format=f"{ny * nx}J", + array=seg_vignets.astype(np.int32, copy=False).reshape(n_obj, -1), + dim=f"({nx},{ny})", + ) + + +# Double-precision positions add_seg_vignet centres SEG_VIGNET on. +DOUBLE_POSITIONS = ("X_IMAGE_DBL", "Y_IMAGE_DBL") + + +def vignet_centre(pos): + """The 0-based pixel SExtractor centres VIGNET on, along one axis. + + SExtractor 2.25.0 (``src/analyse.c``, ``ix=(int)(obj->mx+0.49999)``, + unchanged in Debian's 2.25.0+ds-3) truncates the 0-based barycentre plus + 0.49999. That is not round-to-nearest: a fractional part in + [0.5, 0.50001) goes down, where ``np.rint`` would go up (or to even). + + Parameters + ---------- + pos : array_like of float + 1-based double-precision position (``X_IMAGE_DBL`` or + ``Y_IMAGE_DBL``; SExtractor writes ``mx + 1``) + + Returns + ------- + numpy.ndarray + 0-based pixel index, int64 + + """ + mx = np.asarray(pos, np.float64) - 1.0 + return np.trunc(mx + 0.49999).astype(np.int64) + + +def seg_vignet_param_file(dot_param, output_path): + """Write a SExtractor parameter file that also asks for the double + positions. + + SExtractor centres VIGNET from its double-precision barycentre + (:func:`vignet_centre`), which the float32 ``X_IMAGE`` / ``Y_IMAGE`` + cannot resolve near half pixels; ``X_IMAGE_DBL`` / ``Y_IMAGE_DBL`` + carry it. + :func:`add_seg_vignet` reads them and drops them again. + + Parameters + ---------- + dot_param : str + Path to the configured parameter file + output_path : str + Path to write the extended parameter file to + + Returns + ------- + str + ``output_path`` + + """ + with open(dot_param) as f: + text = f.read() + if text and not text.endswith("\n"): + text += "\n" + with open(output_path, "w") as f: + f.write(text + "".join(f"{name}\n" for name in DOUBLE_POSITIONS)) + return output_path + + +def add_seg_vignet(cat_path, seg_path, w_log=None): + """Add the ``SEG_VIGNET`` column to a SExtractor catalogue. + + The SEGMENTATION check image, whose labels are the catalogue's + ``NUMBER``, is cut on the grid of each object's VIGNET, centred where + SExtractor centres VIGNET (:func:`vignet_centre` of ``X_IMAGE_DBL`` / + ``Y_IMAGE_DBL``, :func:`cut_stamps`), and written, int32 and 0 off the image, as + ``SEG_VIGNET`` in ``LDAC_OBJECTS``, which ngmix's UberSeg blend handling + reads. The double positions (:func:`seg_vignet_param_file`) are dropped; + every other HDU and column is kept. + + Parameters + ---------- + cat_path : str + Path to the SExtractor FITS-LDAC catalogue, rewritten in place + seg_path : str + Path to the SEGMENTATION check image + w_log : logging.Logger, optional + Pipeline logger + + Raises + ------ + ValueError + If the catalogue lacks the double positions + + """ + with fits.open(cat_path) as hdul: + hdus = [hdu.copy() for hdu in hdul] + objects = next(h for h in hdus if h.name == "LDAC_OBJECTS") + data = objects.data + missing = [name for name in DOUBLE_POSITIONS if name not in data.names] + if missing: + raise ValueError( + f"{cat_path} lacks {', '.join(missing)}, which SEG_VIGNET is" + + " centred on; run SExtractor with seg_vignet_param_file." + ) + col = vignet_centre(data["X_IMAGE_DBL"]) + row = vignet_centre(data["Y_IMAGE_DBL"]) + size = data["VIGNET"].shape[1] + seg_vignets = cut_stamps(fits.getdata(seg_path), col, row, size, 0) + kept = [c for c in objects.columns.columns + if c.name not in DOUBLE_POSITIONS] + new = fits.BinTableHDU.from_columns( + kept + [seg_vignet_column(seg_vignets)], + header=objects.header, + name="LDAC_OBJECTS", + ) + fits.HDUList( + [new if h is objects else h for h in hdus] + ).writeto(cat_path, overwrite=True) + if w_log: + centre = size // 2 + own = np.mean(seg_vignets[:, centre, centre] == data["NUMBER"]) + w_log.info( + f"SEG_VIGNET cut from {seg_path} for {len(col)} objects;" + + f" centre label = NUMBER for {own:.4f}" + ) + + def get_header_value(image_path, key): """Get Header Value. @@ -536,6 +720,7 @@ def get_check_image(self, check_image): if (len(check_image) == 1) & (check_image[0] == ""): check_type = ["NONE"] check_name = ["none"] + self.check_paths = {} else: check_type = [] check_name = [] @@ -549,6 +734,7 @@ def get_check_image(self, check_image): + self._num_str + ".fits" ) + self.check_paths = dict(zip(check_type, check_name)) self._cmd_line_extra += ( f' -CHECKIMAGE_TYPE {",".join(check_type)} ' diff --git a/src/shapepipe/modules/sextractor_runner.py b/src/shapepipe/modules/sextractor_runner.py index 65181712f..493bc5c96 100644 --- a/src/shapepipe/modules/sextractor_runner.py +++ b/src/shapepipe/modules/sextractor_runner.py @@ -79,6 +79,24 @@ def sextractor_runner( f_wcs_path = input_file_list[-1] input_file_list = list(input_file_list[:-1]) + # SEG_VIGNET (optional, environment-expanded boolean): add the + # SEGMENTATION check image's stamps, on each VIGNET's grid, as the + # SEG_VIGNET column ngmix's UberSeg blend handling reads. SExtractor then + # also writes the double-precision positions VIGNET is centred on, which + # add_seg_vignet reads and drops. + seg_vignet = config.has_option( + module_config_sec, "SEG_VIGNET" + ) and config.getexpandedboolean(module_config_sec, "SEG_VIGNET") + if seg_vignet: + if "SEGMENTATION" not in [key.upper() for key in check_image]: + raise ValueError( + "SEG_VIGNET needs the SEGMENTATION check image in CHECKIMAGE." + ) + dot_param = ss.seg_vignet_param_file( + dot_param, + f"{run_dirs['tmp']}/seg_vignet{file_number_string}.param", + ) + # Create sextractor caller class instance ss_inst = ss.SExtractorCaller( input_file_list, @@ -110,6 +128,15 @@ def sextractor_runner( # Parse SExtractor errors stdout, stderr = ss_inst.parse_errors(stderr, stdout) + # SEG_VIGNET is cut before the join, which relabels its stamps with the + # rows' new NUMBERs. + if seg_vignet: + ss.add_seg_vignet( + ss_inst.path_output_file, + ss_inst.check_paths["SEGMENTATION"], + w_log=w_log, + ) + # MATCH_CATALOGUE (optional, environment-expanded path; empty for none): # take membership and NUMBER from that external catalogue of the same # image, before the post-processing keys the epoch HDUs on NUMBER. diff --git a/src/shapepipe/pipeline/config.py b/src/shapepipe/pipeline/config.py index 2c294ae65..1f49fe043 100644 --- a/src/shapepipe/pipeline/config.py +++ b/src/shapepipe/pipeline/config.py @@ -97,6 +97,29 @@ def getexpanded(self, section, option, **kwargs): """ return self._get(section, _expandvars_strict, option, **kwargs) + def getexpandedboolean(self, section, option, **kwargs): + """Get Expanded Boolean. + + Expand enviroment variables in the value, then read it as a boolean + the way ``getboolean`` does. + + Parameters + ---------- + section : str + Configuration file section + option : str + Configuration file option + + Returns + ------- + bool + The expanded value as a boolean + + """ + return self._convert_to_boolean( + self.getexpanded(section, option, **kwargs) + ) + def getlist(self, section, option, delimiter=",", **kwargs): """Get List. diff --git a/tests/helpers/defect_response.py b/tests/helpers/defect_response.py new file mode 100644 index 000000000..c3e3ac157 --- /dev/null +++ b/tests/helpers/defect_response.py @@ -0,0 +1,89 @@ +"""Full-matrix metacal recovery for fixed detector defects.""" + +import numpy as np + +from shapepipe.modules.ngmix_package import ngmix as ngm +from shapepipe.testing.simulate import make_data +from tests.helpers.metacal_sim import build_stamp + + +def defect_response(bad, hlr=0.5, psf=0.7, seeds=range(4), + psf_shear=(0.0, 0.0), options=None, known_rms=True, + defect_value=1e3): + """Recover c and M = inverse(mean R) A - I with paired seed errors. + + Null pairs rotate the pixels by 90 degrees while pre-rotating the PSF + ellipticity oppositely, so the final PSF stays fixed in detector + coordinates. + This does NOT average away elliptical-PSF leakage. + + The flagged pixels ``bad`` hold ``defect_value`` (far above the galaxy's + peak) rather than sky, as a hot column or bleed would, so any defect value + that reaches metacal shows up as a bias. + """ + gamma = 0.02 + options = {} if options is None else options + samples = [] + for seed in seeds: + def arm(axis, sign, rotation=0): + shear = [0.0, 0.0] + if axis >= 0: + shear[axis] = sign * gamma + ps = tuple(v * (-1 if rotation else 1) for v in psf_shear) + data = list(make_data( + rng=np.random.RandomState(seed + 100), shear=shear, + psf_shear=ps, noise=1e-4, n_epochs=1, img_size=51, + gal_hlr=hlr, psf_fwhm=psf, return_centers=True, + )) + centre = data.pop()[0] + offset = np.array([centre.y - 26, centre.x - 26]) + if rotation: + data[0] = [np.rot90(a).copy() for a in data[0]] + data[1] = [np.rot90(a).copy() for a in data[1]] + offset = np.array([-offset[1], offset[0]]) + data[0] = [np.where(bad, defect_value, a) for a in data[0]] + data[4] = [bad.astype(np.int32)] + stamp = build_stamp(data) + stamp.offsets = [offset] + if known_rms: + stamp.bkg_rms = [np.full((51, 51), 1e-4)] + rng = np.random.RandomState(seed) + result, _, _ = ngm.do_ngmix_metacal( + stamp, ngm.get_prior(0.1857, rng), 1.0, rng, + centroid_source="wcs", **options, + ) + assert all(result[t]["flags"] == 0 for t in ngm.METACAL_TYPES) + e = np.asarray(result["noshear"]["g"]) + response = np.column_stack([ + (np.asarray(result[p]["g"]) - result[m]["g"]) / 0.02 + for p, m in (("1p", "1m"), ("2p", "2m")) + ]) + assert np.all(np.isfinite(e)) and np.all(np.isfinite(response)) + return e, response + null = [arm(-1, 0, k) for k in (0, 1)] + e0 = np.mean([a[0] for a in null], axis=0) + r0 = np.mean([a[1] for a in null], axis=0) + derivatives, responses = [], [] + for axis in (0, 1): + plus, minus = arm(axis, 1), arm(axis, -1) + derivatives.append((plus[0] - minus[0]) / (2 * gamma)) + responses.extend([plus[1], minus[1]]) + samples.append((e0, r0, np.column_stack(derivatives), + np.mean(responses, axis=0))) + e, r, a, rm = [np.array([s[i] for s in samples]) for i in range(4)] + assert np.linalg.svd(rm.mean(axis=0), compute_uv=False).min() > 0.1 + c = np.linalg.solve(r.mean(axis=0), e.mean(axis=0)) + matrix = np.linalg.solve(rm.mean(axis=0), a.mean(axis=0)) - np.eye(2) + draw = np.random.RandomState(91).randint(len(e), size=(1000, len(e))) + cb = np.linalg.solve(r[draw].mean(axis=1), + e[draw].mean(axis=1)[..., None])[..., 0] + mb = np.linalg.solve(rm[draw].mean(axis=1), a[draw].mean(axis=1)) + mb -= np.eye(2) + return dict(c=c.tolist(), m=np.diag(matrix).tolist(), + matrix=matrix.tolist(), + c_err=cb.std(axis=0).tolist(), + m_err=np.diag(mb.std(axis=0)).tolist(), + response=r.mean(axis=0).tolist(), + seed_groups=[dict(e=s[0].tolist(), R=s[1].tolist(), + A=s[2].tolist(), Rm=s[3].tolist()) + for s in samples]) diff --git a/tests/module/test_defect_interpolation.py b/tests/module/test_defect_interpolation.py new file mode 100644 index 000000000..1368e78cc --- /dev/null +++ b/tests/module/test_defect_interpolation.py @@ -0,0 +1,237 @@ +"""Defect interpolation (``DEFECT_FILL = interpolate``): which pixels are +interpolated, and the properties of the interpolant. + +:func:`interpolable_defects` picks the defect pixels that lie in a short row +or column run with clean pixels at both ends; :func:`interpolate_defects` +fills them from the nearby clean pixels with a Clough-Tocher interpolant, +averaged over the four quarter turns of the stamp and shared by every plane +(science image and metacal noise image). +""" + +import numpy as np +import numpy.testing as npt +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from shapepipe.modules.ngmix_package.defect_interpolation import ( + MAX_INTERPOLATED_RUN, + fourfold, + interpolable_defects, + interpolate_defects, +) + + +# --- interpolable_defects --------------------------------------------------- + +def _oracle(defect, max_run=MAX_INTERPOLATED_RUN): + """Brute force: walk each defect pixel's row and column run.""" + n0, n1 = defect.shape + out = np.zeros_like(defect) + for i, j in zip(*np.nonzero(defect)): + for di, dj in ((0, 1), (1, 0)): + lo, hi = (i, j), (i, j) + while (0 <= lo[0] - di and 0 <= lo[1] - dj + and defect[lo[0] - di, lo[1] - dj]): + lo = (lo[0] - di, lo[1] - dj) + while (hi[0] + di < n0 and hi[1] + dj < n1 + and defect[hi[0] + di, hi[1] + dj]): + hi = (hi[0] + di, hi[1] + dj) + length = hi[0] - lo[0] + hi[1] - lo[1] + 1 + bounded = (lo[0] - di >= 0 and lo[1] - dj >= 0 + and hi[0] + di < n0 and hi[1] + dj < n1) + if bounded and length <= max_run: + out[i, j] = True + return out + + +@given( + n=st.integers(5, 21), + density=st.floats(0.0, 0.6), + seed=st.integers(0, 2**31 - 1), +) +@settings(max_examples=60, deadline=None) +def test_interpolable_defects_are_the_short_bounded_runs(n, density, seed): + """A defect pixel is interpolated exactly when its row or column run is + at most MAX_INTERPOLATED_RUN long and has clean pixels at both ends; the + rule commutes with quarter turns. + + Failure modes: a run touching the stamp border (edge band, corner) or a + wide hole is interpolated from one side; a 3-px bleed is left to noise; + a clean pixel is selected; one axis is ignored, so the rule has a + preferred direction. + """ + defect = np.random.RandomState(seed).uniform(size=(n, n)) < density + out = interpolable_defects(defect) + npt.assert_array_equal(out, _oracle(defect)) + assert not out[~defect].any() + for k in range(1, 4): + npt.assert_array_equal( + interpolable_defects(np.rot90(defect, k)), np.rot90(out, k) + ) + + +@pytest.mark.parametrize("kind,expected", [ + ("column", True), ("bleed3", True), ("finite_bleed3", True), + ("pixel", True), ("bleed4", False), ("blob5", False), + ("edge_band", False), ("corner", False), +]) +def test_calibrated_widths_are_interpolated_and_wider_holes_are_not( + kind, expected, +): + """Columns, 3-px bleeds and single pixels are interpolated; 4-px bleeds, + blobs, edge bands and corners are not.""" + n, c = 51, 25 + defect = np.zeros((n, n), dtype=bool) + region = { + "column": np.s_[:, c + 8], + "bleed3": np.s_[:, c + 8:c + 11], + "finite_bleed3": np.s_[c - 5:c + 6, c + 8:c + 11], + "pixel": np.s_[c, c + 8], + "bleed4": np.s_[:, c + 8:c + 12], + "blob5": np.s_[c - 2:c + 3, c + 8:c + 13], + "edge_band": np.s_[:, -3:], + "corner": np.s_[:2, :2], + }[kind] + defect[region] = True + out = interpolable_defects(defect) + assert out[defect].all() if expected else not out.any() + + +# --- fourfold --------------------------------------------------------------- + +@given(n=st.integers(2, 20), seed=st.integers(0, 2**31 - 1)) +@settings(max_examples=30, deadline=None) +def test_fourfold_is_the_quarter_turn_orbit(n, seed): + """The union contains the mask, is invariant under quarter turns, and is + the smallest such set (idempotent).""" + mask = np.random.RandomState(seed).uniform(size=(n, n)) < 0.1 + out = fourfold(mask) + assert out[mask].all() + for k in range(1, 4): + npt.assert_array_equal(np.rot90(out, k), out) + npt.assert_array_equal(fourfold(out), out) + npt.assert_array_equal( + out, mask | np.rot90(mask) | np.rot90(mask, 2) | np.rot90(mask, 3) + ) + + +def test_fourfold_rejects_rectangular_stamps(): + with pytest.raises(ValueError, match="square"): + fourfold(np.zeros((11, 12), dtype=bool)) + + +# --- interpolate_defects ---------------------------------------------------- + +def _mask(n=31): + """A column, a finite 3-px bleed and a single pixel, all bounded.""" + defect = np.zeros((n, n), dtype=bool) + defect[:, 19] = True + defect[4:12, 7:10] = True + defect[22, 11] = True + return defect + + +@given(seed=st.integers(0, 2**31 - 1)) +@settings(max_examples=20, deadline=None) +def test_interpolation_reproduces_planes_without_reading_defects(seed): + """Linear planes are reproduced at the interpolated pixels; clean pixels, + and defect pixels not selected, are returned untouched; the inputs are + not modified; and no defect value (NaN or a sentinel) is ever read. + + Failure modes: defect pixels enter the support; coordinates are + transposed or mis-rotated; the fill smooths clean pixels. + """ + n = 31 + rng = np.random.RandomState(seed) + rows, cols = np.indices((n, n)) + coeff = rng.uniform(-2, 2, (2, 3)) + planes = np.array([a + b * rows + c * cols for a, b, c in coeff]) + defect = _mask(n) + defect[:, -2:] = True # an edge band and a blob next to the column: + defect[13:18, 21:26] = True # defects that are not targets + target = interpolable_defects(defect) + assert target[15, 19] + assert not target[13:18, 21:26].any() and not target[:, -2:].any() + contaminated = planes.copy() + contaminated[:, defect] = np.nan + saved = contaminated.copy() + + out = interpolate_defects(contaminated, defect, target) + + npt.assert_allclose(out[:, target], planes[:, target], atol=1e-5) + npt.assert_array_equal(out[:, ~target], contaminated[:, ~target]) + npt.assert_array_equal(contaminated, saved) + contaminated[:, defect] = 1e30 + npt.assert_array_equal( + interpolate_defects(contaminated, defect, target)[:, target], + out[:, target], + ) + + +def test_interpolation_commutes_with_quarter_turns(): + """Rotating the stamp and its mask rotates the fill. A regular grid's + Delaunay triangulation has degenerate diagonals, so a single-orientation + interpolant does not commute; the four-orientation average does. + + Failure mode: the average over orientations is skipped, so the fill has a + preferred direction. + """ + n = 31 + rows, cols = np.indices((n, n)) + image = np.exp(-((rows - 15.3) ** 2 + (cols - 14.6) ** 2) / 10.0) + image += 0.05 * np.random.RandomState(3).normal(size=(n, n)) + planes = image[None] + defect = _mask(n) + target = interpolable_defects(defect) + out = interpolate_defects(planes, defect, target) + for k in range(1, 4): + rotated = interpolate_defects( + np.rot90(planes, k, axes=(1, 2)), np.rot90(defect, k), + np.rot90(target, k), + ) + npt.assert_allclose( + rotated, np.rot90(out, k, axes=(1, 2)), atol=1e-12, rtol=0 + ) + + +def test_every_plane_sees_the_same_operator(): + """The fill is one linear operator applied to every plane: filling + a * image + b * noise gives a * fill(image) + b * fill(noise), up to the + Clough-Tocher gradient solver's tolerance. + + Failure mode: the noise image is filled differently from the science + image (another support, triangulation or orientation set), so metacal's + fixnoise no longer mirrors the science image's correlated noise. + """ + n = 31 + rng = np.random.RandomState(11) + image, noise = rng.normal(size=(2, n, n)) + defect = _mask(n) + target = interpolable_defects(defect) + both = interpolate_defects(np.array([image, noise]), defect, target) + mixed = interpolate_defects( + np.array([2.0 * image - 3.0 * noise]), defect, target + ) + npt.assert_allclose(mixed[0], 2.0 * both[0] - 3.0 * both[1], atol=1e-5) + alone = interpolate_defects(noise[None], defect, target) + npt.assert_allclose(alone[0], both[1], atol=1e-5) + + +def test_no_target_is_a_no_op(): + planes = np.random.RandomState(5).normal(size=(2, 15, 15)) + defect = np.zeros((15, 15), dtype=bool) + defect[:, -3:] = True + out = interpolate_defects(planes, defect, np.zeros_like(defect)) + npt.assert_array_equal(out, planes) + + +def test_degenerate_support_gives_nan_not_an_error(): + """With fewer than three non-collinear clean pixels the target is NaN, + which the caller replaces by noise.""" + defect = np.ones((5, 5), dtype=bool) + defect[2, 1] = defect[2, 3] = False + target = np.zeros_like(defect) + target[2, 2] = True + out = interpolate_defects(np.ones((1, 5, 5)), defect, target) + assert np.isnan(out[0, 2, 2]) diff --git a/tests/module/test_ngmix.py b/tests/module/test_ngmix.py index 959d007f1..f0844011f 100644 --- a/tests/module/test_ngmix.py +++ b/tests/module/test_ngmix.py @@ -1,5 +1,7 @@ """UNIT TESTS FOR MODULE PACKAGE: NGMIX.""" +from collections import Counter + from astropy.io import fits from astropy.wcs import WCS import galsim @@ -629,6 +631,7 @@ def test_process_counts_flagged_fits_across_batches(tmp_path, monkeypatch, flags galaxies = {str(i): {"exp-1": {"OFFSET": [0., 0.]}} for i in tile.obj_id} stamp = SimpleNamespace( gals=[np.ones((5, 5))], ra=[42.], dec=[30.], ccd=20, + epoch_cuts=Counter(considered=1), jacobs=[galsim.JacobianWCS(.186, 0., 0., .186)], ) psf = dict( @@ -642,13 +645,14 @@ def test_process_counts_flagged_fits_across_batches(tmp_path, monkeypatch, flags results.append((result, psf, psf)) fits_to_return = iter(results) monkeypatch.setattr(module, "Tile_cat", lambda *args: tile) - monkeypatch.setattr(module, "prepare_postage_stamps", lambda *args: stamp) + monkeypatch.setattr( + module, "prepare_postage_stamps", lambda *args, **kwargs: stamp, + ) monkeypatch.setattr( module, "do_ngmix_metacal", lambda *args, **kwargs: next(fits_to_return), ) inst = object.__new__(Ngmix) inst._tile_cat_path = "in-memory-tile" - inst._seg_cat_path = None inst._vignet_cat = SimpleNamespace( gal_vign_cat=galaxies, psf_vign_cat=galaxies, close=lambda: None, ) @@ -658,6 +662,12 @@ def test_process_counts_flagged_fits_across_batches(tmp_path, monkeypatch, flags inst._blend_handling = "noisefill" inst._dilate_neighbour = 1 inst._metacal_psf = "fitgauss" + inst._epoch_central_defect_radius = module.EPOCH_CENTRAL_DEFECT_RADIUS + inst._epoch_masked_fraction_cut = module.EPOCH_MASKED_FRACTION_CUT + inst._defect_fill = "noise" + inst._epoch_interpolated_defect_radius = ( + module.EPOCH_INTERPOLATED_DEFECT_RADIUS + ) inst._save_batch = 1 inst._zero_point = 30. inst._output_dir = str(tmp_path) @@ -725,6 +735,7 @@ def test_process_centroid_prior_is_each_objects_own_pixel_scale(monkeypatch): obj_id: SimpleNamespace( gals=[np.ones((5, 5))] * len(jacobs), jacobs=jacobs, ra=[10. * obj_id], dec=[30.], ccd=obj_id, + epoch_cuts=Counter(considered=len(jacobs)), ) for obj_id, jacobs in epochs.items() } @@ -737,14 +748,13 @@ def capture(stamp, prior, flux_guess, rng, **kwargs): monkeypatch.setattr(module, "Tile_cat", lambda *args: tile) monkeypatch.setattr( module, "prepare_postage_stamps", - lambda vignet, obj_id, *args: stamps[obj_id], + lambda vignet, obj_id, *args, **kwargs: stamps[obj_id], ) monkeypatch.setattr(module, "do_ngmix_metacal", capture) monkeypatch.setattr(Ngmix, "save_results", lambda self, res: None) monkeypatch.setattr(Ngmix, "log_mean_ellipticity", lambda self: None) inst = object.__new__(Ngmix) inst._tile_cat_path = "in-memory-tile" - inst._seg_cat_path = None inst._vignet_cat = SimpleNamespace( gal_vign_cat=galaxies, psf_vign_cat=galaxies, close=lambda: None, ) @@ -754,6 +764,12 @@ def capture(stamp, prior, flux_guess, rng, **kwargs): inst._blend_handling = "noisefill" inst._dilate_neighbour = 1 inst._metacal_psf = "fitgauss" + inst._epoch_central_defect_radius = module.EPOCH_CENTRAL_DEFECT_RADIUS + inst._epoch_masked_fraction_cut = module.EPOCH_MASKED_FRACTION_CUT + inst._defect_fill = "noise" + inst._epoch_interpolated_defect_radius = ( + module.EPOCH_INTERPOLATED_DEFECT_RADIUS + ) inst._save_batch = -1 inst._w_log = _RecordingLogger() diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py new file mode 100644 index 000000000..e8e39ffa6 --- /dev/null +++ b/tests/module/test_ngmix_defect_fill.py @@ -0,0 +1,1303 @@ +"""Defect fill and the epoch cuts (ngmix module). + +A defect is a stamp pixel with a nonzero flag, zero exposure weight or an +invalid background RMS. Before metacal, :func:`prepare_ngmix_weights` gives +every defect weight 0 and fills it, whatever ``BLEND_HANDLING`` is: with noise +at its background RMS (``DEFECT_FILL = noise``, the default) or, for short +defect runs, with an interpolant of the clean pixels around them +(``DEFECT_FILL = interpolate``). Under uberseg, pixels on the neighbour side +only lose their weight, and their image values stay raw. The per-epoch cuts +in :func:`prepare_postage_stamps` act on the same defect set: the +masked-fraction cut counts it, and the central-defect veto drops an epoch +with a defect near the stamp centre, at a radius set by the defect's fill. +The tile VIGNET's -1e30 neighbour markers are not defects: +noisefill zero-weights and noise-fills them, uberseg ignores them, and the +epoch cuts never count them. +""" + +import re +from pathlib import Path +from types import SimpleNamespace + +import numpy as np +import numpy.testing as npt +import pytest +from astropy.io import fits +from astropy.wcs import WCS +from hypothesis import given +from hypothesis import strategies as st +from modopt.math.stats import sigma_mad +from sqlitedict import SqliteDict +from shapepipe.modules.ngmix_package import ngmix as ngmix_module + +from shapepipe.modules.ngmix_package.defect_interpolation import ( + fourfold, + interpolable_defects, + interpolate_defects, +) +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_INTERPOLATED_DEFECT_RADIUS, + Ngmix, + make_ngmix_observation, + prepare_ngmix_weights, + prepare_postage_stamps, + split_tile_markers, + uberseg_weight, +) +from shapepipe.modules.sextractor_package.sextractor_script import cut_stamps + + +# --- prepare_ngmix_weights: the filled set is the defect set --------------- + +@st.composite +def defect_stamps(draw): + """Square stamp with random flagged, zero-weight and bad-RMS pixels.""" + n = draw(st.integers(min_value=5, max_value=21)) + pixels = st.lists( + st.tuples(st.integers(0, n - 1), st.integers(0, n - 1)), max_size=n + ) + weight = np.ones((n, n)) + flag = np.zeros((n, n), dtype=np.int32) + bkg_rms = np.ones((n, n)) + for i, j in draw(pixels): + weight[i, j] = 0.0 + for i, j in draw(pixels): + flag[i, j] = draw(st.sampled_from([1, 2, 2**10])) + for i, j in draw(pixels): + bkg_rms[i, j] = draw(st.sampled_from([0.0, -1.0, np.nan, np.inf])) + # Raw values far outside the unit-RMS noise, so a filled pixel is + # unambiguous. + gal = 1.0e3 + np.arange(n * n, dtype=float).reshape(n, n) + return gal, weight, flag, bkg_rms + + +def _uberseg_seg(n): + """Central object on the centre pixel, neighbour footprint in a corner.""" + seg = np.zeros((n, n), dtype=np.int32) + seg[n // 2, n // 2] = 1 + seg[:2, :2] = 2 + return seg + + +@given( + stamp=defect_stamps(), + blend_handling=st.sampled_from(["noisefill", "uberseg"]), + seed=st.integers(0, 2**31 - 1), +) +def test_filled_set_is_the_defect_set(stamp, blend_handling, seed): + """Filled pixels are exactly the defects (flag, zero weight, bad RMS). + They carry zero weight and look like noise. Every other pixel keeps its + raw value, under either BLEND_HANDLING. + + Failure modes: + * a defect source (flag, zero weight, bad RMS) is left out of the fill; + * the filled set grows beyond the defects (for example symmetrized); + * the filled set and the zero-weight defect set differ; + * the fill is skipped under uberseg; + * neighbour-side pixels are filled. + """ + gal, weight, flag, bkg_rms = stamp + n = gal.shape[0] + kwargs = ( + dict(seg=_uberseg_seg(n), object_number=1, dilate_neighbour=1) + if blend_handling == "uberseg" + else {} + ) + + gal_out, w_out, _ = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(seed), bkg_rms=bkg_rms, + blend_handling=blend_handling, **kwargs, + ) + + defect = ( + (weight == 0) | (flag != 0) | ~(np.isfinite(bkg_rms) & (bkg_rms > 0)) + ) + if defect.all(): + return # fully masked: no clean pixel left to set the noise level + filled = gal_out != gal + npt.assert_array_equal(filled, defect, "filled set is not the defect set") + assert np.all(np.abs(gal_out[filled]) < 10.0), "fill is not unit noise" + neighbour = ( + uberseg_weight(np.ones((n, n)), kwargs["seg"], 1, dilate_neighbour=1) + == 0.0 + if blend_handling == "uberseg" + else np.zeros((n, n), dtype=bool) + ) + # Zero weight on defects and neighbour side; neighbour pixels that are + # not defects keep their raw values (filled-set equality above). + npt.assert_array_equal(w_out == 0.0, defect | neighbour) + npt.assert_array_equal(w_out[~(defect | neighbour)], 1.0) + + +def test_uberseg_defect_in_neighbour_region_is_filled(): + """A defect pixel that also lies on the neighbour side is filled. + + Failure mode: the fill is restricted to pixels uberseg keeps, so a raw bad + pixel on the neighbour side still reaches metacal. + """ + n = 21 + gal = 1.0e3 + np.arange(n * n, dtype=float).reshape(n, n) + weight = np.ones((n, n)) + flag = np.zeros((n, n), dtype=np.int32) + flag[1, 1] = 1 # inside the neighbour footprint + gal_out, w_out, _ = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(0), bkg_rms=np.ones((n, n)), + blend_handling="uberseg", seg=_uberseg_seg(n), object_number=1, + ) + assert w_out[1, 1] == 0.0 and gal_out[1, 1] != gal[1, 1] + # A neighbour-side pixel that is not a defect: zero weight, raw value. + assert w_out[0, 3] == 0.0 and gal_out[0, 3] == gal[0, 3] + + +def test_committed_blend_handling_fills_defects(): + """The committed universe's blend_handling x defect_fill pair is what + prepare_ngmix_weights does: the workflow's default blend_handling is the + committed one, and under it every defect is zero-weighted and filled. + + Failure mode: the committed blend handling skips the defect fill, so the + record claims a fill the default campaign does not run (raw defects, + which astra excludes, reaching metacal). + """ + import yaml + + repo = Path(__file__).resolve().parents[2] + universe = yaml.safe_load((repo / "universes" / "committed.yaml").read_text()) + decisions = universe["analyses"]["shape_measurement"]["decisions"] + blend_handling = decisions["blend_handling"] + defect_fill = decisions["defect_fill"] + assert blend_handling in ngmix_module.BLEND_HANDLINGS + assert defect_fill in ngmix_module.DEFECT_FILLS + workflow = yaml.safe_load((repo / "workflow" / "config.yaml").read_text()) + assert workflow["blend_handling"] == blend_handling + + n = 21 + gal = 1.0e3 + np.arange(n * n, dtype=float).reshape(n, n) + weight = np.ones((n, n)) + flag = np.zeros((n, n), dtype=np.int32) + flag[3, 15] = 1 + flag[10, 2] = 2**10 + weight[17, 9] = 0.0 + defect = (weight == 0) | (flag != 0) + kwargs = ( + dict(seg=_uberseg_seg(n), object_number=1) + if blend_handling == "uberseg" + else {} + ) + gal_out, w_out, _ = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(0), bkg_rms=np.ones((n, n)), + blend_handling=blend_handling, defect_fill=defect_fill, **kwargs, + ) + assert np.all(w_out[defect] == 0.0) + assert np.all(gal_out[defect] != gal[defect]), "defects are left raw" + assert np.all(np.abs(gal_out[defect]) < 10.0), "fill is not unit noise" + + +# --- prepare_postage_stamps: the fraction cut counts the defect set -------- + +N_STAMP = 51 +RA, DEC = 150.0, 2.0 + + +def _fake_inputs(epochs): + """Minimal vignet / tile-catalogue stand-ins for prepare_postage_stamps. + + ``epochs`` maps ``"-"`` to ``(flag, weight)`` or + ``(flag, weight, bkg_rms)``. Returns ``(vignet, tile_cat, psf_obj, + gal_obj)``. + """ + rng = np.random.default_rng(1) + wcs = WCS(naxis=2) + wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] + wcs.wcs.crval = [RA, DEC] + wcs.wcs.crpix = [N_STAMP / 2, N_STAMP / 2] + wcs.wcs.cdelt = [-0.187 / 3600, 0.187 / 3600] + header = fits.Header({"FSCALE": 1.0}).tostring() + + def per_epoch(make): + return {k: {"VIGNET": make(k)} for k in epochs} + + with_rms = any(len(v) == 3 for v in epochs.values()) + psf_obj = per_epoch(lambda k: np.ones((N_STAMP, N_STAMP))) + gal_obj = per_epoch(lambda k: rng.normal(0.0, 1.0, (N_STAMP, N_STAMP))) + for epoch in gal_obj.values(): + epoch["OFFSET"] = np.zeros(2) + vignet = SimpleNamespace( + gal_vign_cat={"1": gal_obj}, + bkg_vign_cat=None, + bkg_rms_vign_cat=( + {"1": per_epoch( + lambda k: epochs[k][2] if len(epochs[k]) == 3 + else np.ones((N_STAMP, N_STAMP)) + )} + if with_rms + else None + ), + flag_vign_cat={"1": per_epoch(lambda k: epochs[k][0])}, + weight_vign_cat={"1": per_epoch(lambda k: epochs[k][1])}, + f_wcs_file={ + k.split("-")[0]: { + int(k.split("-")[1]): {"WCS": wcs, "header": header} + } + for k in epochs + }, + ) + tile_cat = SimpleNamespace( + vign=None, seg=None, ra=np.array([RA]), dec=np.array([DEC]) + ) + return vignet, tile_cat, psf_obj, gal_obj + + +def _two_sided_band(width): + """Mask of ``width`` columns on each side of the stamp, far from the + centre (at least 17 px for width 9).""" + band = np.zeros((N_STAMP, N_STAMP), dtype=bool) + band[:, :width] = True + band[:, -width:] = True + return band + + +def _epochs(): + """A clean epoch and four masked ones, all defects far from the centre. + + * band10: 10 flagged columns on one side, 19.6% (39% if symmetrized); + * flag18: 2 x 9 flagged columns, 35.3%; + * dead18: the same columns at zero weight, no flags; + * rms18: the same columns with an invalid background RMS. + """ + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + band10 = clean.copy() + band10[:, -10:] = 1 + two = _two_sided_band(9) + dead = ones.copy() + dead[two] = 0.0 + rms = ones.copy() + rms[two] = np.nan + return { + "2100001-10": (clean, ones), + "2100002-11": (band10, ones), + "2100003-12": (two.astype(np.int32), ones), + "2100004-13": (clean.copy(), dead), + "2100005-14": (clean.copy(), ones, rms), + } + + +def _surviving(epochs, **kwargs): + vignet, tile_cat, psf_obj, gal_obj = _fake_inputs(epochs) + stamp = prepare_postage_stamps( + vignet, 1, 0, tile_cat, bkg_sub=False, + psf_obj=psf_obj, gal_obj=gal_obj, **kwargs, + ) + names = {id(v[0]): name for name, v in epochs.items()} + return sorted(names[id(flag)] for flag in stamp.flags) + + +def test_epoch_cut_counts_the_defect_set(): + """At the default 1/3 cut, the three 35% epochs are dropped, whichever + defect source masks them, and the 20% band survives. + + Failure modes: the cut counts flags only (keeps dead18 and rms18), omits + one defect source, or counts a symmetrized set (drops band10); in each + case the cut disagrees with the set that is zero-weighted and filled + (epoch-cut-on-defect-mask). + """ + assert _surviving(_epochs()) == ["2100001-10", "2100002-11"] + + +def test_epoch_cut_threshold_is_the_configured_fraction(): + """At a 10% cut (the DES Y3 and Y6 value), the 19.6% band is dropped as + well, and only the clean epoch survives. + + Failure mode: the configured threshold is ignored. + """ + assert _surviving(_epochs(), epoch_masked_fraction_cut=0.1) == [ + "2100001-10" + ] + + +# --- prepare_postage_stamps: the central-defect veto ----------------------- + +def _veto_epochs(radius): + """A clean epoch, one with a single flagged pixel just inside + ``radius`` of the stamp centre, and one with a flagged column exactly + ``radius`` away. Both masked epochs are far below the fraction cut. + """ + centre = N_STAMP // 2 + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + near, far = clean.copy(), clean.copy() + near[centre, centre + radius - 1] = 1 + far[:, centre + radius] = 1 + return { + "2100001-10": (clean, ones), + "2100002-11": (near, ones), + "2100003-12": (far, ones), + } + + +def test_central_defect_vetoes_the_epoch(): + """At the default radius, a single defect pixel inside it drops the + epoch, and a column at the radius does not. + + Failure mode: the veto is skipped, so an epoch whose filled hole overlaps + the object's light enters the fit; or the boundary is inclusive, dropping + the column at the radius (epoch-central-defect-veto). + """ + epochs = _veto_epochs(EPOCH_CENTRAL_DEFECT_RADIUS) + assert _surviving(epochs) == ["2100001-10", "2100003-12"] + + +def test_central_defect_radius_is_the_configured_value(): + """Radius 0 disables the veto; a radius one pixel larger than the far + column's distance drops that epoch too. + + Failure mode: the configured radius is ignored. + """ + epochs = _veto_epochs(EPOCH_CENTRAL_DEFECT_RADIUS) + assert _surviving(epochs, epoch_central_defect_radius=0) == [ + "2100001-10", "2100002-11", "2100003-12" + ] + assert _surviving( + epochs, epoch_central_defect_radius=EPOCH_CENTRAL_DEFECT_RADIUS + 1 + ) == ["2100001-10"] + + +# --- prepare_postage_stamps: per-epoch OFFSET ------------------------------ + +def test_each_surviving_epoch_carries_its_own_offset(): + """``stamp.offsets`` holds each surviving epoch's vignette OFFSET, in the + order of ``stamp.flags``. + + Failure mode: the offset is dropped or read from another epoch, so the + default "wcs" centroid raises or puts the Jacobian origin off the object. + """ + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + epochs = {f"210000{i}-1{i}": (clean.copy(), ones) for i in range(3)} + vignet, tile_cat, psf_obj, gal_obj = _fake_inputs(epochs) + del vignet.gal_vign_cat # OFFSET must come from the supplied gal_obj. + for i, name in enumerate(epochs): + gal_obj[name]["OFFSET"] = np.array([0.1 * i, -0.1 * i]) + stamp = prepare_postage_stamps( + vignet, 1, 0, tile_cat, bkg_sub=False, + psf_obj=psf_obj, gal_obj=gal_obj, + ) + names = {id(flag): name for name, (flag, _) in epochs.items()} + assert len(stamp.offsets) == len(stamp.flags) == 3 + for flag, offset in zip(stamp.flags, stamp.offsets): + npt.assert_array_equal(offset, gal_obj[names[id(flag)]]["OFFSET"]) + + +# --- Ngmix.process: the per-tile epoch-cut tally --------------------------- + +class _RecordingLogger: + def __init__(self): + self.messages = [] + + def info(self, msg, *_args, **_kwargs): + self.messages.append(msg) + + warning = error = info + + +def test_process_logs_the_epoch_cut_tally(tmp_path, monkeypatch): + """One tile, four objects; the end-of-tile line counts each cut's drops. + + * object 1: clean, 18-column edge band, defect inside the radius -> one + epoch each for considered, masked_fraction, central_veto; survives. + * object 2: edge band and central defect -> both epochs dropped; emptied. + * object 3: one all-zero stamp, skipped before the cuts -> not considered, + and not emptied by the cuts. + * object 4: no PSF ('empty') -> never reaches the cuts. + + Failure modes: a cut's drops are not counted or land in the wrong + counter; epochs skipped before the cuts are counted as considered; an + object with no epoch at all is reported as emptied by the cuts; counts + from one object overwrite another's. + """ + radius = EPOCH_CENTRAL_DEFECT_RADIUS + centre = N_STAMP // 2 + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + wide18, near = clean.copy(), clean.copy() + wide18[:, -18:] = 1 + near[centre, centre + radius - 1] = 1 + objects = { + 1: {"2100001-10": (clean, ones), "2100002-11": (wide18, ones), + "2100003-12": (near, ones)}, + 2: {"2100004-13": (wide18.copy(), ones), + "2100005-14": (near.copy(), ones)}, + 3: {"2100006-15": (clean.copy(), ones)}, + } + stores = {} + for obj_id, epochs in objects.items(): + vignet, _, psf_obj, gal_obj = _fake_inputs(epochs) + if obj_id == 3: + gal_obj["2100006-15"]["VIGNET"] = np.zeros((N_STAMP, N_STAMP)) + stores[obj_id] = (vignet, psf_obj, gal_obj) + vignet = SimpleNamespace( + bkg_vign_cat=None, + bkg_rms_vign_cat=None, + psf_vign_cat={ + "4": "empty", **{str(i): s[1] for i, s in stores.items()} + }, + gal_vign_cat={ + "4": "empty", **{str(i): s[2] for i, s in stores.items()} + }, + flag_vign_cat={ + str(i): s[0].flag_vign_cat["1"] for i, s in stores.items() + }, + weight_vign_cat={ + str(i): s[0].weight_vign_cat["1"] for i, s in stores.items() + }, + f_wcs_file={ + k: v for s in stores.values() for k, v in s[0].f_wcs_file.items() + }, + close=lambda: None, + ) + tile_cat = SimpleNamespace( + obj_id=np.array([1, 2, 3, 4]), ra=np.full(4, RA), dec=np.full(4, DEC), + vign=None, seg=None, flux=None, + ) + + paths = [tmp_path / f"{name}.sqlite" for name in + ("gal", "psf", "weight", "flag", "headers")] + for path in paths: + SqliteDict(str(path)).close() + log = _RecordingLogger() + ngmix = Ngmix( + ["tile_cat.fits"] + [str(p) for p in paths[:4]], + str(tmp_path), "-001-001", 30.0, str(paths[4]), log, + bkg_sub=False, + ) + ngmix._vignet_cat.close() + ngmix._vignet_cat = vignet + monkeypatch.setattr(ngmix_module, "Tile_cat", lambda *a, **k: tile_cat) + + def no_fit(*_args, **_kwargs): + raise RuntimeError("metacal is not under test") + + monkeypatch.setattr(ngmix_module, "do_ngmix_metacal", no_fit) + for method in ("compile_results", "save_results", "log_mean_ellipticity"): + monkeypatch.setattr(Ngmix, method, lambda *_a, **_k: None) + + ngmix.process() + + lines = [m for m in log.messages if m.startswith("epoch cuts:")] + assert len(lines) == 1, log.messages + tally = dict( + (k, int(v)) for k, v in re.findall(r"(\w+)=(\d+)", lines[0]) + ) + assert tally == dict( + considered=5, masked_fraction=2, central_veto=2, objects_emptied=1 + ), lines[0] + + +# --- ngmix_runner: the epoch-cut options reach Ngmix ------------------------ + +class _OptionConfig: + """Config stub answering from a dict; absent options take the caller's + fallback.""" + + def __init__(self, options): + self._options = {"MAG_ZP": "30.0", "ID_OBJ_MIN": "-1", + "ID_OBJ_MAX": "-1", **options} + + def has_option(self, _sec, key): + return key in self._options + + def get(self, _sec, key): + return self._options[key] + + def getexpanded(self, _sec, key): + return self._options[key] + + def getfloat(self, _sec, key): + return float(self._options[key]) + + def getint(self, _sec, key): + return int(self._options[key]) + + def getboolean(self, _sec, key, fallback=False): + return fallback + + +def test_runner_threads_the_epoch_cut_options(tmp_path, monkeypatch): + """EPOCH_CENTRAL_DEFECT_RADIUS, EPOCH_MASKED_FRACTION_CUT, DEFECT_FILL and + EPOCH_INTERPOLATED_DEFECT_RADIUS reach Ngmix as configured, and default + to the module constants (and the noise fill) when absent. + + Failure mode: the runner drops or ignores an option, so a configured A/B + arm silently runs the default cuts. + """ + from shapepipe.modules import ngmix_runner as runner_module + + captured = [] + + class _Capture: + def __init__(self, *_args, **kwargs): + captured.append(kwargs) + + def process(self): + pass + + monkeypatch.setattr(runner_module, "Ngmix", _Capture) + inputs = [str(tmp_path / f"in{i}.sqlite") for i in range(7)] + for path in inputs: + SqliteDict(path).close() + + for options, radius, fraction, fill, interpolated_radius in ( + ({"EPOCH_CENTRAL_DEFECT_RADIUS": "7.5", + "EPOCH_MASKED_FRACTION_CUT": "0.1", + "DEFECT_FILL": "interpolate", + "EPOCH_INTERPOLATED_DEFECT_RADIUS": "5.5"}, + 7.5, 0.1, "interpolate", 5.5), + ({}, EPOCH_CENTRAL_DEFECT_RADIUS, + ngmix_module.EPOCH_MASKED_FRACTION_CUT, "noise", + EPOCH_INTERPOLATED_DEFECT_RADIUS), + ): + runner_module.ngmix_runner( + inputs, {"output": str(tmp_path)}, "-001-001", + _OptionConfig(options), "NGMIX_RUNNER", _RecordingLogger(), + ) + assert captured[-1]["epoch_central_defect_radius"] == radius + assert captured[-1]["epoch_masked_fraction_cut"] == fraction + assert captured[-1]["defect_fill"] == fill + assert (captured[-1]["epoch_interpolated_defect_radius"] + == interpolated_radius) + + +# --- make_ngmix_observation: the HSM centroid reads the filled image ------- + +def test_hsm_centroid_ignores_raw_defect_values(): + """With ``centroid_source="hsm"``, the Jacobian origin lands on the + object even when a flagged column a few pixels away holds a raw value + ten times the object's peak. + + Failure mode: HSM measures the raw stamp, so the defect drags the centroid + off the object (or HSM fails and falls back to the stamp centre). + """ + import galsim + + n = 51 + centre = (n - 1) / 2 + d_row, d_col = 1.3, -0.7 + rows, cols = np.mgrid[:n, :n] + gal = np.exp( + -((rows - centre - d_row) ** 2 + (cols - centre - d_col) ** 2) + / (2 * 2.0 ** 2) + ) + flag = np.zeros((n, n), dtype=np.int32) + flag[:, n // 2 + 6] = 1 + gal[flag != 0] = 10.0 + psf = np.exp(-((rows - centre) ** 2 + (cols - centre) ** 2) / 2.0) + obs = make_ngmix_observation( + gal, np.ones((n, n)), flag, psf / psf.sum(), + galsim.PixelScale(0.1857).jacobian(), np.random.RandomState(0), + bkg_rms=np.full((n, n), 1e-3), centroid_source="hsm", + ) + row, col = obs.jacobian.get_cen() + npt.assert_allclose( + [row - centre, col - centre], [d_row, d_col], atol=0.05 + ) + + +# --- DEFECT_FILL = interpolate ---------------------------------------------- + +def _hot_stamp(): + """A bright object with hot defects: a column and a finite 3-px bleed + (interpolated), an edge band and a 5x5 blob (noise-filled).""" + centre = N_STAMP // 2 + rows, cols = np.mgrid[:N_STAMP, :N_STAMP] + gal = 1e3 * np.exp( + -((rows - centre) ** 2 + (cols - centre) ** 2) / (2 * 3.0 ** 2) + ) + flag = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + flag[:, centre + 8] = 1 + flag[centre - 5:centre + 6, centre - 11:centre - 8] = 1 + flag[:, :4] = 1 + flag[centre + 9:centre + 14, centre + 12:centre + 17] = 1 + gal[flag != 0] = 5e4 + return gal, np.ones((N_STAMP, N_STAMP)), flag + + +@pytest.mark.parametrize("blend_handling", ["noisefill", "uberseg"]) +def test_interpolated_fill_and_its_weights(blend_handling): + """Short defect runs take the interpolant of the clean image; other + defects take noise; the weight is zero on the defects and on the + quarter-turn orbit of the interpolated pixels, whose light stays; the + metacal noise image is interpolated with the same operator. + + Failure modes: the orbit is not zero-weighted (a one-sided hole in the + likelihood biases c), or its light is replaced; the fill mask is + symmetrized; wide defects or edge bands are extrapolated; raw defect + values leak; the noise image keeps independent noise where the science + image is smooth. + """ + gal, weight, flag = _hot_stamp() + defect = flag != 0 + target = interpolable_defects(defect) + assert target.any() and (defect & ~target).any() + kwargs = ( + dict(seg=_uberseg_seg(N_STAMP), object_number=1, dilate_neighbour=1) + if blend_handling == "uberseg" + else {} + ) + gal_out, w_out, noise_out = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(4), + bkg_rms=np.ones((N_STAMP, N_STAMP)), blend_handling=blend_handling, + defect_fill="interpolate", **kwargs, + ) + + neighbour = ( + uberseg_weight(np.ones_like(gal), kwargs["seg"], 1, dilate_neighbour=1) + == 0.0 + if kwargs else np.zeros_like(defect) + ) + npt.assert_array_equal(w_out == 0.0, defect | fourfold(target) | neighbour) + npt.assert_array_equal(gal_out[~defect], gal[~defect]) + expected = interpolate_defects(gal[None], defect, target)[0] + npt.assert_allclose(gal_out[target], expected[target], rtol=1e-5) + assert np.all(np.abs(gal_out[defect & ~target]) < 10.0) + refilled = interpolate_defects(noise_out[None], defect, target)[0] + npt.assert_allclose(noise_out[target], refilled[target], atol=1e-5) + + +def test_noise_is_the_default_fill(): + """Leaving DEFECT_FILL unset is bit-identical to the noise fill.""" + gal, weight, flag = _hot_stamp() + rms = np.ones_like(gal) + default = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(9), bkg_rms=rms + ) + noise = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(9), bkg_rms=rms, + defect_fill="noise", + ) + for a, b in zip(default, noise): + npt.assert_array_equal(a, b) + + +def test_unknown_defect_fill_is_rejected(): + gal, weight, flag = _hot_stamp() + with pytest.raises(ValueError, match="DEFECT_FILL"): + prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(0), + defect_fill="interp", + ) + + +def _fill_veto_epochs(): + """A clean epoch; a column at the interpolated-defect radius; a column + one pixel inside it; a 5-px bleed (noise-filled) at the same radius.""" + centre = N_STAMP // 2 + radius = int(EPOCH_INTERPOLATED_DEFECT_RADIUS) + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + at, inside, wide = clean.copy(), clean.copy(), clean.copy() + at[:, centre + radius] = 1 + inside[:, centre + radius - 1] = 1 + wide[:, centre + radius:centre + radius + 5] = 1 + return { + "2100001-10": (clean, ones), + "2100002-11": (at, ones), + "2100003-12": (inside, ones), + "2100004-13": (wide, ones), + } + + +def test_the_veto_radius_follows_the_fill(): + """Under interpolation, an interpolated defect is vetoed inside + EPOCH_INTERPOLATED_DEFECT_RADIUS and a noise-filled one inside + EPOCH_CENTRAL_DEFECT_RADIUS. Under the noise fill every defect keeps the + noise radius. The masked-fraction cut counts the raw defect set in + both modes. + + Failure modes: the interpolated radius is applied to noise-filled pixels + (a wide hole next to the object survives), or not applied at all; the + configured radius is ignored; the boundary is inclusive; the fraction + cut counts the zero-weight orbit. + """ + epochs = _fill_veto_epochs() + assert _surviving(epochs, defect_fill="interpolate") == [ + "2100001-10", "2100002-11" + ] + assert _surviving(epochs) == ["2100001-10"] + assert _surviving( + epochs, defect_fill="interpolate", + epoch_interpolated_defect_radius=EPOCH_INTERPOLATED_DEFECT_RADIUS + 1, + ) == ["2100001-10"] + assert _surviving(_epochs(), defect_fill="interpolate") == [ + "2100001-10", "2100002-11" + ] + + +def test_process_threads_the_defect_fill(tmp_path, monkeypatch): + """Ngmix hands DEFECT_FILL to both the epoch cuts and metacal: under + interpolation the column at the interpolated-defect radius survives and + metacal is asked to interpolate; under noise it is vetoed. + + Failure mode: the option is dropped on the way to either, so the tile + silently runs the noise fill or the noise cuts. + """ + epochs = {k: v for k, v in _fill_veto_epochs().items() + if k in ("2100001-10", "2100002-11")} + vignet, tile_cat, psf_obj, _ = _fake_inputs(epochs) + vignet.psf_vign_cat = {"1": psf_obj} + vignet.close = lambda: None + tile_cat.obj_id = np.array([1]) + tile_cat.flux = None + monkeypatch.setattr(ngmix_module, "Tile_cat", lambda *a, **k: tile_cat) + for method in ("compile_results", "save_results", "log_mean_ellipticity"): + monkeypatch.setattr(Ngmix, method, lambda *_a, **_k: None) + calls = [] + + def record(stamp, *_args, **kwargs): + calls.append((len(stamp.gals), kwargs.get("defect_fill"))) + raise RuntimeError("metacal is not under test") + + monkeypatch.setattr(ngmix_module, "do_ngmix_metacal", record) + paths = [tmp_path / f"{name}.sqlite" for name in + ("gal", "psf", "weight", "flag", "headers")] + for path in paths: + SqliteDict(str(path)).close() + for fill in ("interpolate", "noise"): + ngmix = Ngmix( + ["tile_cat.fits"] + [str(p) for p in paths[:4]], + str(tmp_path), "-001-001", 30.0, str(paths[4]), + _RecordingLogger(), bkg_sub=False, defect_fill=fill, + ) + ngmix._vignet_cat.close() + ngmix._vignet_cat = vignet + ngmix.process() + assert calls == [(2, "interpolate"), (1, "noise")] + + +def test_ngmix_rejects_an_unknown_defect_fill(tmp_path): + paths = [tmp_path / f"{name}.sqlite" for name in + ("gal", "psf", "weight", "flag", "headers")] + for path in paths: + SqliteDict(str(path)).close() + with pytest.raises(ValueError, match="DEFECT_FILL"): + Ngmix( + ["tile_cat.fits"] + [str(p) for p in paths[:4]], + str(tmp_path), "-001-001", 30.0, str(paths[4]), + _RecordingLogger(), bkg_sub=False, defect_fill="interp", + ) + + +# --- The tile VIGNET's -1e30 neighbour markers are not defects ------------- +# +# The tile VIGNET carries -1e30 on the footprints of other detections. Every +# epoch shares that tile stamp, so a marker counted as a defect would drop +# every epoch of an object with a neighbour inside the veto radius. The +# markers are their own per-epoch mask (``stamp.neighbours``): noisefill +# zero-weights and noise-fills them, uberseg ignores them, and the epoch cuts +# never read them. + +_MARKER = -1.0e30 +_CENTRE = N_STAMP // 2 +# Flipped (ccd < 18) and unflipped (ccd >= 18) MegaCam CCDs. +_MARKER_EPOCH_NAMES = ["2100001-10", "2100002-20", "2100003-11"] + + +def _tile_with_neighbour(columns_from=_CENTRE + 3, rows=(_CENTRE - 2, + _CENTRE + 3)): + """Tile VIGNET with a -1e30 neighbour footprint whose nearest pixel is + 3 px from the stamp centre. The footprint is off-centre, so the MegaCam + flip moves it.""" + tile = np.random.default_rng(5).normal(0.0, 1.0, (N_STAMP, N_STAMP)) + tile[rows[0]:rows[1], columns_from:columns_from + 6] = _MARKER + return tile + + +def _marker_stamp(tile, epochs=None, **kwargs): + """Run prepare_postage_stamps on defect-free epochs under ``tile``.""" + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + if epochs is None: + epochs = {name: (clean.copy(), ones) for name in _MARKER_EPOCH_NAMES} + vignet, tile_cat, psf_obj, gal_obj = _fake_inputs(epochs) + tile_cat.vign = tile[np.newaxis] + stamp = prepare_postage_stamps( + vignet, 1, 0, tile_cat, bkg_sub=False, + psf_obj=psf_obj, gal_obj=gal_obj, **kwargs, + ) + return stamp, epochs, gal_obj + + +def _expected_neighbours(tile, name): + return Ngmix.MegaCamFlip(tile, int(name.split("-")[1])) == _MARKER + + +def test_neighbour_markers_near_the_centre_keep_every_epoch(): + """A neighbour footprint 3 px from the centre, and one covering 41% of + the stamp, drop no epoch: the masked-fraction cut and the central veto + count no marker. + + Failure mode: the markers are written into the flag stamp and counted as + defects, so every epoch (they all share the tile VIGNET) is dropped by + the veto or the fraction cut and the object loses its shape + (neighbour-markers-are-not-defects). + """ + small = _tile_with_neighbour() + # Large, but short of the stamp border: no row or column is entirely + # marked, so it is a neighbour, not off-tile. + large = _tile_with_neighbour() + large[1:-1, _CENTRE + 3:-1] = _MARKER + assert (large == _MARKER).mean() > 1 / 3 + for tile in (small, large): + stamp, _, _ = _marker_stamp(tile) + assert len(stamp.gals) == len(_MARKER_EPOCH_NAMES) + assert stamp.epoch_cuts["considered"] == len(_MARKER_EPOCH_NAMES) + assert stamp.epoch_cuts["masked_fraction"] == 0 + assert stamp.epoch_cuts["central_veto"] == 0 + + +def test_neighbour_markers_are_their_own_per_epoch_mask(): + """``stamp.neighbours`` holds the MegaCam-flipped marker mask of each + surviving epoch, and the flag stamps stay the exposure's own. + + Failure mode: the markers are merged into the flags, or the neighbour + mask is not flipped with its epoch and lands on the wrong pixels. + """ + tile = _tile_with_neighbour() + stamp, epochs, _ = _marker_stamp(tile) + names = {id(v[0]): name for name, v in epochs.items()} + assert len(stamp.neighbours) == len(stamp.flags) + for flag, neighbour in zip(stamp.flags, stamp.neighbours): + name = names[id(flag)] + npt.assert_array_equal(flag, 0) + npt.assert_array_equal(neighbour, _expected_neighbours(tile, name)) + assert not np.array_equal(stamp.neighbours[0], stamp.neighbours[1]) + + +def test_a_flagged_column_near_the_centre_is_still_vetoed(): + """With a neighbour footprint present, an epoch with a genuinely flagged + column 3 px from the centre is still dropped, and only that epoch. + + Failure mode: handling the markers apart also exempts real defects from + the central veto. + """ + clean = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + ones = np.ones((N_STAMP, N_STAMP)) + column = clean.copy() + column[:, _CENTRE - 3] = 1 + epochs = { + "2100001-10": (clean, ones), + "2100002-20": (column, ones), + "2100003-11": (clean.copy(), ones), + } + stamp, _, _ = _marker_stamp(_tile_with_neighbour(), epochs) + names = {id(v[0]): name for name, v in epochs.items()} + assert sorted(names[id(f)] for f in stamp.flags) == [ + "2100001-10", "2100003-11" + ] + assert stamp.epoch_cuts["central_veto"] == 1 + assert stamp.epoch_cuts["masked_fraction"] == 0 + + +def _stamp_epoch_weights(stamp, i, seed, **kwargs): + return prepare_ngmix_weights( + 1.0e3 + stamp.gals[i], stamp.weights[i], stamp.flags[i], + np.random.RandomState(seed), bkg_rms=stamp.bkg_rms[i], + neighbour=stamp.neighbours[i], **kwargs, + ) + + +def test_noisefill_fills_exactly_the_marked_pixels(): + """Under noisefill, a defect-free epoch has zero weight and noise + exactly on the marked neighbour pixels; every other pixel keeps its raw + value and its weight. + + Failure mode: the neighbour markers are dropped with the flags, so + noisefill no longer removes neighbour light (noisefill-fills-markers). + """ + tile = _tile_with_neighbour() + stamp, epochs, _ = _marker_stamp(tile) + assert len(stamp.gals) == len(_MARKER_EPOCH_NAMES) + for i in range(len(stamp.gals)): + gal = 1.0e3 + stamp.gals[i] + gal_out, w_out, _ = _stamp_epoch_weights( + stamp, i, seed=i, blend_handling="noisefill", + ) + neighbour = stamp.neighbours[i] + assert neighbour.any() + npt.assert_array_equal(gal_out != gal, neighbour) + npt.assert_array_equal(w_out == 0.0, neighbour) + assert np.all(np.abs(gal_out[neighbour]) < 10.0) + + +def test_uberseg_leaves_the_marked_pixels_raw(): + """Under uberseg, the markers mask nothing: with a seg map holding only + the central object, a defect-free epoch keeps every pixel raw and + weighted, marked or not. + + Failure mode: the markers reach the defect set or the fill, so uberseg + noise-fills the neighbour's light instead of leaving it to the seg-based + weight (uberseg-ignores-markers). + """ + tile = _tile_with_neighbour() + stamp, _, _ = _marker_stamp(tile) + seg = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + seg[_CENTRE - 1:_CENTRE + 2, _CENTRE - 1:_CENTRE + 2] = 1 + assert len(stamp.gals) == len(_MARKER_EPOCH_NAMES) + for i in range(len(stamp.gals)): + gal = 1.0e3 + stamp.gals[i] + gal_out, w_out, _ = _stamp_epoch_weights( + stamp, i, seed=i, blend_handling="uberseg", seg=seg, + object_number=1, + ) + npt.assert_array_equal(gal_out, gal) + assert np.all(w_out > 0.0) + + +def _develop_noisefill(gal, weight, flag, rng, bkg_rms=None): + """prepare_ngmix_weights under BLEND_HANDLING = noisefill on develop + (240b37e4), where the markers arrive as flag 2**10: the reference the + noisefill output is pinned to.""" + mask = np.copy(weight) != 0 + mask[flag != 0] = False + if bkg_rms is None: + sig_noise = sigma_mad(gal) + weight_map = mask.astype(float) / sig_noise ** 2 + else: + valid_rms = np.isfinite(bkg_rms) & (bkg_rms > 0) + mask &= valid_rms + weight_map = np.zeros_like(gal, dtype=float) + weight_map[mask] = 1.0 / bkg_rms[mask] ** 2 + sig_noise = np.where(valid_rms, bkg_rms, np.median(bkg_rms[mask])) + noise_img = rng.standard_normal(gal.shape) * sig_noise + noise_img_gal = rng.standard_normal(gal.shape) * sig_noise + gal_masked = np.copy(gal) + gal_masked[~mask] = noise_img_gal[~mask] + return gal_masked, weight_map, noise_img + + +@given( + seed=st.integers(0, 2**31 - 1), + with_rms=st.booleans(), + rms_scale=st.floats(0.5, 2.0), + with_defects=st.booleans(), + dtype=st.sampled_from([np.float64, np.float32]), +) +def test_noisefill_matches_develop_on_marked_neighbours( + seed, with_rms, rms_scale, with_defects, dtype, +): + """For an epoch with a marked neighbour, noisefill returns the image, + weight and noise image develop returned, bit for bit: with no defect, + and with a flagged column and a dead pixel under the default noise fill, + for float64 and float32 stamps (the filled image keeps the stamp's + dtype). + The stamps carry no off-tile pixels, and the equivalence is claimed for + such stamps only. Off-tile pixels reach this function as flag 2**10 + defects (off-tile-pixels-are-defects), as every marker did on develop; + which epochs survive the cuts differs from develop wherever there are + neighbour markers. + + Failure mode: carrying the markers apart from the flags changes what + noisefill does to neighbour pixels, their weights, the noise level or + the RNG stream (noisefill-matches-develop). + """ + rng = np.random.default_rng(seed) + n = 31 + gal = rng.normal(0.0, 1.0, (n, n)) + gal[n // 2 - 2:n // 2 + 3, n // 2 - 2:n // 2 + 3] += 50.0 + weight = np.ones((n, n)) + flag = np.zeros((n, n), dtype=np.int32) + neighbour = np.zeros((n, n), dtype=bool) + neighbour[n // 2 - 3:n // 2 + 4, n // 2 + 3:n // 2 + 9] = True + gal[neighbour] += 30.0 + if with_defects: + flag[:, 2] = 1 + weight[n - 3, n // 2] = 0.0 + gal = gal.astype(dtype) + bkg_rms = ( + rms_scale * (1.0 + 0.1 * rng.random((n, n))) if with_rms else None + ) + + new = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(seed), bkg_rms=bkg_rms, + blend_handling="noisefill", neighbour=neighbour, + ) + old = _develop_noisefill( + gal, weight, np.where(neighbour, 2**10, flag), + np.random.RandomState(seed), bkg_rms=bkg_rms, + ) + for a, b in zip(new, old): + assert a.dtype == b.dtype + npt.assert_array_equal(a, b) + + +def test_do_ngmix_metacal_threads_each_epochs_neighbour_mask(monkeypatch): + """Each epoch's neighbour mask reaches make_ngmix_observation. + + Failure mode: the mask is built but never used, so noisefill silently + stops filling neighbours. + """ + tile = _tile_with_neighbour() + stamp, _, _ = _marker_stamp(tile) + assert len(stamp.gals) == len(_MARKER_EPOCH_NAMES) + seen = [] + + class _Stop(Exception): + pass + + def fake_observation(*args, **kwargs): + seen.append(kwargs["neighbour"]) + if len(seen) == len(stamp.gals): + raise _Stop + return None + + monkeypatch.setattr( + ngmix_module, "make_ngmix_observation", fake_observation, + ) + monkeypatch.setattr(ngmix_module, "ObsList", list) + with pytest.raises(_Stop): + ngmix_module.do_ngmix_metacal( + stamp, None, 1.0, np.random.RandomState(0), + ) + for got, want in zip(seen, stamp.neighbours): + assert got is want + + +# --- Off-tile pixels are defects ------------------------------------------- +# +# The tile VIGNET also holds -1e30 beyond the tile's edge, where the epoch +# holds the object's own light, cut off. Those pixels are the runs of +# entirely -1e30 stamp rows and columns that start at a stamp border (the +# off-image part of a rectangle clip); +# they join the epoch's defect set as flag 2**10. The other markers are the +# neighbour mask. + +_OFF_TILE = 2**10 + + +def _off_tile_expected(tile, name): + flipped = Ngmix.MegaCamFlip(tile, int(name.split("-")[1])) == _MARKER + return flipped.all(axis=1)[:, None] | flipped.all(axis=0)[None, :] + + +def test_object_three_px_from_the_tile_edge_is_dropped(): + """With the tile edge 3 px from the object, every epoch is dropped: the + off-tile band counts toward the epoch cuts. + + Failure mode: off-tile pixels are treated as neighbour markers, so an + edge object is measured with a noise-filled band through its own light + (off-tile-pixels-are-defects). + """ + tile = np.random.default_rng(5).normal(0.0, 1.0, (N_STAMP, N_STAMP)) + tile[:, :_CENTRE - 2] = _MARKER + stamp, _, _ = _marker_stamp(tile) + assert len(stamp.gals) == 0 + assert stamp.epoch_cuts["considered"] == len(_MARKER_EPOCH_NAMES) + assert ( + stamp.epoch_cuts["masked_fraction"] + stamp.epoch_cuts["central_veto"] + == len(_MARKER_EPOCH_NAMES) + ) + + +def test_the_central_veto_sees_off_tile_pixels(): + """An off-tile band 12 px from the object (27% of the stamp) passes the + default cuts and is vetoed at radius 13. + + Failure mode: the central veto does not read the off-tile set. + """ + tile = np.random.default_rng(5).normal(0.0, 1.0, (N_STAMP, N_STAMP)) + tile[:, :_CENTRE - 11] = _MARKER + kept, _, _ = _marker_stamp(tile) + assert len(kept.gals) == len(_MARKER_EPOCH_NAMES) + vetoed, _, _ = _marker_stamp(tile, epoch_central_defect_radius=13) + assert len(vetoed.gals) == 0 + assert vetoed.epoch_cuts["central_veto"] == len(_MARKER_EPOCH_NAMES) + + +def test_corner_off_tile_region_and_border_neighbour_are_classified(): + """At a tile corner, the L-shaped off-tile region is flagged 2**10 + exactly, and a neighbour footprint touching the stamp border without + filling a row or column stays in the neighbour mask. + + Failure modes: off-tile pixels are classified by something other than + whole marked rows and columns (the L is missed or a border-touching + neighbour is swallowed); the classification ignores the MegaCam flip. + """ + tile = np.random.default_rng(5).normal(0.0, 1.0, (N_STAMP, N_STAMP)) + tile[:5, :] = _MARKER + tile[:, -5:] = _MARKER + tile[40:, :6] = _MARKER # neighbour on the bottom-left border + stamp, epochs, _ = _marker_stamp(tile) + names = {id(v[0]): name for name, v in epochs.items()} + assert len(stamp.gals) == len(_MARKER_EPOCH_NAMES) + for flag, neighbour in zip(stamp.flags, stamp.neighbours): + name = names[id(flag)] + off_tile = _off_tile_expected(tile, name) + assert off_tile.sum() == 5 * N_STAMP * 2 - 25 + npt.assert_array_equal(flag == _OFF_TILE, off_tile) + npt.assert_array_equal( + neighbour, _expected_neighbours(tile, name) & ~off_tile + ) + assert neighbour.sum() == 11 * 6 + + +@pytest.mark.parametrize("blend_handling", ["noisefill", "uberseg"]) +def test_off_tile_pixels_are_zero_weighted_and_filled(blend_handling): + """Off-tile pixels are zero-weighted and noise-filled under either + BLEND_HANDLING; under uberseg the neighbour markers stay raw. + + Failure mode: under uberseg the off-tile band keeps its weight, or the + fill differs between blend handlings. + """ + tile = np.random.default_rng(5).normal(0.0, 1.0, (N_STAMP, N_STAMP)) + tile[:5, :] = _MARKER + tile[20:24, 35:40] = _MARKER + stamp, _, _ = _marker_stamp(tile) + seg = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + seg[_CENTRE - 1:_CENTRE + 2, _CENTRE - 1:_CENTRE + 2] = 1 + kwargs = ( + dict(seg=seg, object_number=1) if blend_handling == "uberseg" else {} + ) + for i in range(len(stamp.gals)): + gal = 1.0e3 + stamp.gals[i] + gal_out, w_out, _ = _stamp_epoch_weights( + stamp, i, seed=i, blend_handling=blend_handling, **kwargs, + ) + off_tile = stamp.flags[i] == _OFF_TILE + assert off_tile.sum() == 5 * N_STAMP + removed = off_tile | ( + stamp.neighbours[i] if blend_handling == "noisefill" else False + ) + npt.assert_array_equal(gal_out != gal, removed) + npt.assert_array_equal(w_out == 0.0, removed) + + +# --- Tile VIGNETs marked as SExtractor marks them --------------------------- +# +# SExtractor writes -1e30 into the tile VIGNET off the image and on the +# segmentation footprints of other detections. The off-image pixels are the +# off-tile defects and the footprint pixels are the neighbour mask, exactly. + +DR6_PATCH = Path(__file__).parent / "data" / "dr6_202.301_seg_patch.fits" + + +def _marked_stamps(seg, x, y): + """Tile VIGNETs on a unit image, marked -1e30 off the image and on every + footprint but the one under the stamp's centre pixel, and each stamp's + off-image mask.""" + col = np.rint(np.asarray(x)).astype(np.int64) - 1 + row = np.rint(np.asarray(y)).astype(np.int64) - 1 + seg_stamps = cut_stamps(seg.astype(np.int64), col, row, N_STAMP, -1) + off_image = [s == -1 for s in seg_stamps] + own = seg[row, col] + vignets = np.ones(seg_stamps.shape, np.float32) + for vign, s, o in zip(vignets, seg_stamps, own): + vign[(s == -1) | ((s != 0) & (s != o))] = _MARKER + return vignets, off_image + + +def test_dr6_marked_stamps_split_into_off_image_and_neighbours(): + """On the real 202.301 segmentation patch, every marked stamp splits into + its off-image pixels (off-tile) and its other -1e30 pixels (neighbours). + + Failure modes: a float32 -1e30 is not recognised as a + marker; off-image pixels of an edge stamp land in the neighbour mask; + neighbour-footprint pixels become off-tile defects. + """ + with fits.open(DR6_PATCH) as hdul: + seg = hdul["SEG"].data + objects = hdul["OBJECTS"].data + x, y = np.array(objects["X_IMAGE"]), np.array(objects["Y_IMAGE"]) + vignets, off_image = _marked_stamps(seg, x, y) + assert vignets.dtype == np.float32 + n_edge = 0 + for vign, off in zip(vignets, off_image): + neighbour, off_tile = split_tile_markers(vign, vign.shape) + npt.assert_array_equal(off_tile, off) + npt.assert_array_equal(neighbour, (vign == _MARKER) & ~off) + n_edge += off.any() + assert n_edge >= 5 + assert sum( + split_tile_markers(v, v.shape)[0].sum() for v in vignets + ) > 1000 + + +def test_a_neighbour_completing_rows_beside_the_tile_edge_stays_a_neighbour(): + """An object 14 px from the tile's left edge, with a wide neighbour + footprint that runs from the tile edge across the stamp: in the rows of + that footprint every stamp pixel is -1e30, off the image or on the + neighbour. Only the off-image columns are off-tile; the footprint is the + neighbour mask, through prepare_postage_stamps. + + Failure mode: every entirely -1e30 row counts as off-tile, so the + neighbour's rows become defects that the epoch cuts count and the defect + fill interpolates (off-tile-is-marked-border-rows-and-columns). + """ + seg = np.zeros((80, 80), np.int32) + seg[20:24, 0:45] = 5 + seg[27:32, 13:18] = 1 + x, y = np.array([15.0, 30.0]), np.array([30.0, 22.0]) + vignets, off_image = _marked_stamps(seg, x, y) + tile, off = vignets[0], off_image[0] + footprint = (tile == _MARKER) & ~off + assert footprint.sum() == 4 * (N_STAMP - 11) + assert ((tile == _MARKER).all(axis=1) & ~off.all(axis=1)).sum() == 4 + + stamp, epochs, _ = _marker_stamp(tile) + names = {id(v[0]): name for name, v in epochs.items()} + assert len(stamp.flags) == len(_MARKER_EPOCH_NAMES) + for flag, neighbour in zip(stamp.flags, stamp.neighbours): + ccd = int(names[id(flag)].split("-")[1]) + npt.assert_array_equal(flag == _OFF_TILE, Ngmix.MegaCamFlip(off, ccd)) + npt.assert_array_equal(neighbour, Ngmix.MegaCamFlip(footprint, ccd)) + + +# --- DEFECT_FILL = interpolate beside a removed neighbour ------------------- + +def _defect_beside_neighbour(): + """A single flagged pixel 6 px right of the centre, with a bright marked + neighbour footprint starting on the next column.""" + n = 31 + c = n // 2 + gal = np.random.default_rng(3).normal(0.0, 1.0, (n, n)) + flag = np.zeros((n, n), dtype=np.int32) + flag[c, c + 6] = 1 + neighbour = np.zeros((n, n), dtype=bool) + neighbour[c - 2:c + 3, c + 7:c + 10] = True + gal[neighbour] += 1.0e4 + seg = np.zeros((n, n), dtype=np.int32) + seg[c - 1:c + 2, c - 1:c + 2] = 1 + return gal, flag, neighbour, seg, (c, c + 6) + + +def test_noisefill_interpolation_does_not_read_removed_neighbour_light(): + """Under noisefill, the interpolant of a defect beside a marked neighbour + is built from the pixels the image keeps: the neighbour's light, which + noisefill removes, is not in its support. Under uberseg the neighbour + light is raw in the image and supports the interpolant. + + Failure mode: the support includes the removed neighbour pixels, so the + defect is filled with light the image no longer contains + (interpolation-support-is-the-kept-image). + """ + gal, flag, neighbour, seg, pix = _defect_beside_neighbour() + weight = np.ones_like(gal) + target = np.zeros_like(neighbour) + target[pix] = True + + out, _, _ = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(0), + blend_handling="noisefill", neighbour=neighbour, + defect_fill="interpolate", + ) + expected = interpolate_defects(gal[None], (flag != 0) | neighbour, target) + assert out[pix] == pytest.approx(expected[0][pix]) + assert abs(out[pix]) < 100.0 + + out, _, _ = prepare_ngmix_weights( + gal, weight, flag, np.random.RandomState(0), + blend_handling="uberseg", seg=seg, object_number=1, + neighbour=neighbour, defect_fill="interpolate", + ) + expected = interpolate_defects(gal[None], flag != 0, target) + assert out[pix] == pytest.approx(expected[0][pix]) + assert out[pix] > 1000.0 diff --git a/tests/module/test_ngmix_uberseg.py b/tests/module/test_ngmix_uberseg.py index 7e7790c55..f846eb28a 100644 --- a/tests/module/test_ngmix_uberseg.py +++ b/tests/module/test_ngmix_uberseg.py @@ -6,10 +6,10 @@ Geometry assertions on a synthetic two-object stamp: neighbour-side pixels are zeroed, the surviving central core is a *single connected* region (the emergent "circularisation"), and the neighbour footprint is fully removed. -* :func:`prepare_ngmix_weights` — the ``noisefill`` default is byte-for-byte - unchanged (asserted against an independent recomputation of the legacy - three-line noise-fill on a shared RNG), while ``uberseg`` hard-masks the - weight (weight -> 0) and leaves the image untouched. +* :func:`prepare_ngmix_weights` under ``uberseg`` — neighbour-side pixels + lose their weight and keep their raw image values, while defect pixels are + noise-filled as under any blend handling (the defect fill itself is covered + in ``test_ngmix_defect_fill.py``). * The error contract when ``uberseg`` is selected without a segmentation map (the seg-map source is plumbing-gated upstream). """ @@ -17,11 +17,13 @@ import numpy as np import numpy.testing as npt import pytest +from astropy.io import fits from scipy import ndimage from sqlitedict import SqliteDict from shapepipe.modules.ngmix_package.ngmix import ( Ngmix, + Tile_cat, central_seg_label, prepare_ngmix_weights, seg_has_neighbour, @@ -226,7 +228,7 @@ def test_uberseg_matches_bruteforce_nearest_segment(): npt.assert_array_equal(out, brute) -# --- prepare_ngmix_weights: default unchanged, uberseg hard-masks ---------- +# --- prepare_ngmix_weights: uberseg zeroes neighbour weights only --------- def _gal_flag_weight(npix=41, seed=7): rng = np.random.default_rng(seed) @@ -239,35 +241,9 @@ def _gal_flag_weight(npix=41, seed=7): return gal, flag, weight -def test_noisefill_default_is_byte_identical_to_legacy(): - """The default path reproduces the legacy three-line noise-fill exactly - (same RNG stream): masked pixels replaced by noise, weight 1/sigma^2.""" - gal, flag, weight = _gal_flag_weight() - - gal_out, w_out, noise_out = prepare_ngmix_weights( - gal, weight, flag, np.random.RandomState(123), - ) - - # Independent recomputation of the legacy algorithm on the same stream. - from modopt.math.stats import sigma_mad - rng = np.random.RandomState(123) - mask = np.copy(weight) != 0 - mask[flag != 0] = False - sig = sigma_mad(gal) - w_exp = mask.astype(float) / sig ** 2 - noise_exp = rng.standard_normal(gal.shape) * sig - noise_gal = rng.standard_normal(gal.shape) * sig - gal_exp = np.copy(gal) - gal_exp[~mask] = noise_gal[~mask] - - npt.assert_array_equal(gal_out, gal_exp) - npt.assert_array_equal(w_out, w_exp) - npt.assert_array_equal(noise_out, noise_exp) - - def test_noisefill_ignores_seg_and_dilate_kwargs(): - """Under noisefill, passing seg / dilate_neighbour changes nothing: the - result matches the plain default call on the same RNG stream.""" + """Under BLEND_HANDLING = noisefill, passing seg / dilate_neighbour changes + nothing: the result matches the plain default call on the same RNG stream.""" gal, flag, weight = _gal_flag_weight() seg, _, _ = two_object_seg(npix=gal.shape[0], sep=12) @@ -282,9 +258,16 @@ def test_noisefill_ignores_seg_and_dilate_kwargs(): npt.assert_array_equal(a, b) -def test_uberseg_hard_masks_weight_and_leaves_image_untouched(): - """uberseg: image returned untouched, weight zeroed on neighbour-side and - flagged pixels, positive on the central core.""" +def test_uberseg_fills_defects_and_leaves_neighbour_pixels_raw(): + """uberseg: flagged pixels are noise-filled at weight 0; neighbour-side + pixels get weight 0 and keep their raw image values; the central core + keeps weight and image. + + Failure modes: the defect fill is skipped under uberseg, so raw bad + pixels reach metacal; or the neighbour side is noise-filled + (defects-filled-whatever-the-blend-handling, + uberseg-ignores-markers). + """ npix = 41 gal, flag, weight = _gal_flag_weight(npix=npix) seg, centre, neigh = two_object_seg(npix=npix, sep=12) @@ -294,13 +277,16 @@ def test_uberseg_hard_masks_weight_and_leaves_image_untouched(): blend_handling="uberseg", seg=seg, object_number=1, ) - # Image untouched under uberseg (no noise fill). - npt.assert_array_equal(gal_out, gal) - # Neighbour footprint hard-masked; central centre kept. + # Flagged pixels: zero weight, image replaced by noise. + for pix in [(5, 5), (30, 12)]: + assert w_out[pix] == 0.0 + assert gal_out[pix] != gal[pix] + # Neighbour footprint: zero weight, raw image (never noise-filled). assert np.all(w_out[seg == 2] == 0.0) + npt.assert_array_equal(gal_out[seg == 2], gal[seg == 2]) + # Central core: weight and image untouched. assert w_out[centre] > 0.0 - # Flagged bad pixels remain at weight 0 (folded into the base mask). - assert w_out[5, 5] == 0.0 + assert gal_out[centre] == gal[centre] def test_uberseg_requires_seg_and_object_number(): @@ -366,59 +352,108 @@ def test_check_central_seg_label_missing_number_raises(tmp_path): ngmix._check_central_seg_label(seg, obj_id=99) -# --- FIX 3a: uberseg without seg_cat_path fails at construction ------------- - -def test_ngmix_init_uberseg_without_seg_cat_raises(tmp_path): - """blend_handling='uberseg' with seg_cat_path=None raises at __init__.""" +# --- the seg stamps ride the tile catalogue as SEG_VIGNET ------------------- + +def _write_tile_cat(path, seg_vignets=None): + """A three-object LDAC tile catalogue, with SEG_VIGNET when given.""" + n_obj, size = 3, 5 + cols = [ + fits.Column(name="NUMBER", format="J", array=np.array([4, 5, 9])), + fits.Column(name="XWIN_WORLD", format="D", array=np.zeros(n_obj)), + fits.Column(name="YWIN_WORLD", format="D", array=np.zeros(n_obj)), + fits.Column(name="VIGNET", format=f"{size * size}E", + array=np.ones((n_obj, size * size), np.float32), + dim=f"({size},{size})"), + ] + if seg_vignets is not None: + cols.append(fits.Column(name="SEG_VIGNET", format=f"{size * size}J", + array=seg_vignets.reshape(n_obj, -1), + dim=f"({size},{size})")) + imhead = fits.BinTableHDU.from_columns( + [fits.Column(name="Field Header Card", format="1A", array=["x"])], + name="LDAC_IMHEAD", + ) + fits.HDUList([ + fits.PrimaryHDU(), imhead, + fits.BinTableHDU.from_columns(cols, name="LDAC_OBJECTS"), + ]).writeto(path) + return path + + +def test_tile_cat_reads_seg_vignet(tmp_path): + """Tile_cat's seg stamps are the catalogue's SEG_VIGNET, row for row; + a catalogue without the column has none.""" + seg = np.arange(3 * 25, dtype=np.int32).reshape(3, 5, 5) + with_seg = Tile_cat(str(_write_tile_cat(tmp_path / "a.fits", seg))) + npt.assert_array_equal(with_seg.seg, seg) + assert with_seg.seg.dtype.kind == "i" + assert Tile_cat(str(_write_tile_cat(tmp_path / "b.fits"))).seg is None + + +def test_tile_cat_holds_the_stamp_columns_once(tmp_path): + """VIGNET and SEG_VIGNET are views into the one loaded table: copying + either would double the bulk of every ngmix chunk's catalogue memory.""" + seg = np.zeros((3, 5, 5), np.int32) + tile = Tile_cat(str(_write_tile_cat(tmp_path / "a.fits", seg))) + # Interleaved fields of one record buffer: their extents overlap. + assert np.may_share_memory(tile.vign, tile.seg) + + +def test_uberseg_without_seg_vignet_fails_loudly(tmp_path): + """BLEND_HANDLING = uberseg on a catalogue without SEG_VIGNET raises + before any object is fitted, rather than dropping every object.""" + cat = _write_tile_cat(tmp_path / "tile_cat.fits") names = ("gal", "bkg", "psf", "weight", "flag", "headers") paths = [tmp_path / f"{name}.sqlite" for name in names] for path in paths: SqliteDict(str(path)).close() - with pytest.raises(ValueError, match="requires SEG_VIGNET_PATH"): - Ngmix( - ["tile_cat.fits"] + [str(p) for p in paths[:5]], - str(tmp_path), "-001-001", 30.0, str(paths[5]), - _RecordingLogger(), - blend_handling="uberseg", - seg_cat_path=None, - ) - - -# --- FIX 3b: runner raises when SEG_VIGNET_PATH is set but missing ---------- - -class _FakeConfig: - """Config stub: SEG_VIGNET_PATH is the only present option, and it resolves - to ``seg_path`` (a path the test leaves nonexistent).""" - - def __init__(self, seg_path): - self._seg_path = seg_path - - def getfloat(self, _sec, _key): - return 30.0 - - def getboolean(self, _sec, _key, fallback=False): - return fallback - - def has_option(self, _sec, key): - return key == "SEG_VIGNET_PATH" - - def getexpanded(self, _sec, _key): - return self._seg_path - - -def test_runner_missing_seg_vignet_file_raises(tmp_path): - """SEG_VIGNET_PATH configured but the resolved file is absent -> the runner - fails fast with FileNotFoundError, before constructing Ngmix.""" - from shapepipe.modules.ngmix_runner import ngmix_runner - - input_file_list = ["tile_cat.fits"] + [f"in{i}.sqlite" for i in range(6)] - seg_path = str(tmp_path / "does_not_exist.fits") - with pytest.raises(FileNotFoundError, match="Segmentation vignet file"): - ngmix_runner( - input_file_list, - {"output": str(tmp_path)}, - "-001-001", - _FakeConfig(seg_path), - "NGMIX_RUNNER", - _RecordingLogger(), - ) + ngmix = Ngmix( + [str(cat)] + [str(p) for p in paths[:5]], + str(tmp_path), "-001-001", 30.0, str(paths[5]), + _RecordingLogger(), blend_handling="uberseg", + ) + with pytest.raises(ValueError, match="SEG_VIGNET"): + ngmix.process() + + +def test_runner_reads_blend_handling_from_the_environment( + tmp_path, monkeypatch, +): + """The committed ini's ``${SP_BLEND_HANDLING:-noisefill}`` reaches Ngmix + as noisefill when the variable is unset and as its value when set.""" + from shapepipe.modules import ngmix_runner as runner_module + from shapepipe.pipeline.config import CustomParser + + names = ("cat", "gal", "bkg", "psf", "weight", "flag", "headers") + paths = [str(tmp_path / f"{name}.sqlite") for name in names] + for path in paths[1:]: + SqliteDict(path).close() + config = CustomParser() + config.add_section("NGMIX_RUNNER") + for key, value in {"MAG_ZP": "30.0", "ID_OBJ_MIN": "-1", + "ID_OBJ_MAX": "-1", + "BLEND_HANDLING": "${SP_BLEND_HANDLING:-noisefill}", + }.items(): + config.set("NGMIX_RUNNER", key, value) + + seen = [] + + class _Stop(Exception): + pass + + def fake_ngmix(*_args, **kwargs): + seen.append(kwargs["blend_handling"]) + raise _Stop + + monkeypatch.setattr(runner_module, "Ngmix", fake_ngmix) + runner = getattr(runner_module.ngmix_runner, "__wrapped__", + runner_module.ngmix_runner) + for env in (None, "uberseg"): + if env is None: + monkeypatch.delenv("SP_BLEND_HANDLING", raising=False) + else: + monkeypatch.setenv("SP_BLEND_HANDLING", env) + with pytest.raises(_Stop): + runner(paths, {"output": str(tmp_path)}, "-001-001", config, + "NGMIX_RUNNER", _RecordingLogger()) + assert seen == ["noisefill", "uberseg"] diff --git a/tests/module/test_pipeline.py b/tests/module/test_pipeline.py index 400ca2ec3..125f24dc6 100644 --- a/tests/module/test_pipeline.py +++ b/tests/module/test_pipeline.py @@ -134,6 +134,46 @@ def test_custom_parser_getexpanded_expands_fallback(monkeypatch): ) +def test_custom_parser_expands_a_default_only_when_the_variable_is_unset( + monkeypatch, +): + """``${VAR:-default}`` takes ``default`` when VAR is unset or empty, and + VAR's value otherwise; a bare unset ``$VAR`` still raises.""" + parser = config.CustomParser() + parser.add_section("S") + parser.set("S", "MODE", "${SP_TEST_MODE:-noisefill}") + parser.set("S", "FLAG", "${SP_TEST_FLAG:-False}") + parser.set("S", "BARE", "$SP_TEST_MODE") + + monkeypatch.delenv("SP_TEST_MODE", raising=False) + monkeypatch.delenv("SP_TEST_FLAG", raising=False) + assert parser.getexpanded("S", "MODE") == "noisefill" + assert parser.getexpandedboolean("S", "FLAG") is False + with pytest.raises(ValueError, match="SP_TEST_MODE"): + parser.getexpanded("S", "BARE") + + monkeypatch.setenv("SP_TEST_MODE", "") + assert parser.getexpanded("S", "MODE") == "noisefill" + + monkeypatch.setenv("SP_TEST_MODE", "uberseg") + monkeypatch.setenv("SP_TEST_FLAG", "True") + assert parser.getexpanded("S", "MODE") == "uberseg" + assert parser.getexpanded("S", "BARE") == "uberseg" + assert parser.getexpandedboolean("S", "FLAG") is True + + +def test_custom_parser_expands_a_variable_value_once(monkeypatch): + """A variable's value is inserted as is, never expanded again.""" + monkeypatch.setenv("SP_TEST_ROOT", "/data/$SP_TEST_OTHER") + monkeypatch.setenv("SP_TEST_OTHER", "expanded-again") + parser = config.CustomParser() + parser.add_section("S") + parser.set("S", "PATH", "${SP_TEST_ROOT:-/x}/a $SP_TEST_ROOT/b") + assert parser.getexpanded("S", "PATH") == ( + "/data/$SP_TEST_OTHER/a /data/$SP_TEST_OTHER/b" + ) + + def test_custom_parser_getlist_honours_custom_delimiter(): parser = config.CustomParser() diff --git a/tests/module/test_sextractor_seg_vignet.py b/tests/module/test_sextractor_seg_vignet.py new file mode 100644 index 000000000..f66faa99a --- /dev/null +++ b/tests/module/test_sextractor_seg_vignet.py @@ -0,0 +1,453 @@ +"""UNIT TESTS FOR SEXTRACTOR MODE'S SEG_VIGNET. + +``sextractor_script.add_seg_vignet`` cuts the SEGMENTATION check image on the +grid SExtractor cut each object's VIGNET on and adds it to the sexcat as the +int32 ``SEG_VIGNET`` column. SExtractor 2.25.0 centres VIGNET on +``(int)(mx + 0.49999)`` of its 0-based double-precision barycentre +(``src/analyse.c``), which the catalogue's float32 X_IMAGE / Y_IMAGE cannot +resolve near half pixels, so the runner asks SExtractor for X_IMAGE_DBL / +Y_IMAGE_DBL too (``seg_vignet_param_file``), and ``add_seg_vignet`` centres +on them and then drops them. + +The expected centres below are written out by hand from that rule, and the +expected stamps are cut by plain slicing, not by the module's helpers. +""" + +import numpy as np +import numpy.testing as npt +import pytest +from astropy.io import fits + +from shapepipe.modules.sextractor_package import sextractor_script as ss + +NX, NY = 60, 50 +STAMP = 7 +BIG = np.float32(-1e30) + +# (NUMBER, X_IMAGE_DBL, Y_IMAGE_DBL, VIGNET centre col, row (0-based)). +# Object 3's x is an exact half with an odd floor (np.rint goes to 32 - 1), +# object 4's y has a fractional part in [0.5, 0.50001) (np.rint goes up), and +# object 6's x is just past 0.50001, where float32 falls back below it. +OBJECTS = [ + (1, 20.2, 15.7, 19, 15), + (2, 1.3, 48.9, 0, 48), + (3, 32.5, 30.1, 31, 29), + (4, 44.8, 22.5000099, 44, 21), + (5, 59.9, 1.2, 59, 0), + (6, 31.5000101, 10.2, 31, 9), +] +DBL = ["X_IMAGE_DBL", "Y_IMAGE_DBL"] + + +@pytest.mark.parametrize("pos, centre", [ + (32.5, 31), # exact half, odd floor: down, not to even + (31.5, 30), # exact half, even floor + (31.49999, 30), + (31.5000099, 30), # fraction 0.5000099 < 0.50001: down + (31.5000101, 31), # fraction 0.5000101 > 0.50001: up + (31.4999, 30), + (20.2, 19), + (20.8, 20), + (0.6, 0), # mx = -0.4 truncates to 0 +]) +def test_vignet_centre_is_sextractors_rule(pos, centre): + assert ss.vignet_centre([pos]).tolist() == [centre] + + +def _cut(array, col, row, fill): + """The STAMP stamp of ``array`` with 0-based (row, col) at its centre.""" + half = STAMP // 2 + padded = np.pad(array, half, constant_values=fill) + return padded[row:row + STAMP, col:col + STAMP] + + +def _scene(tmp_path, flat=False): + """A sexcat whose VIGNETs are cut on the hand-written centres. + + ``flat`` makes the image uniform, so no pixel comparison can tell a tie + object's two candidate centres apart. + """ + rng = np.random.default_rng(3) + image = (np.full((NY, NX), 100, np.float32) if flat + else rng.normal(100, 5, (NY, NX)).astype(np.float32)) + seg = np.zeros((NY, NX), np.int32) + for number, _, _, col, row in OBJECTS: + seg[max(row - 2, 0):row + 3, max(col - 2, 0):col + 3] = number + number = np.array([o[0] for o in OBJECTS]) + x = np.array([o[1] for o in OBJECTS]) + y = np.array([o[2] for o in OBJECTS]) + + vignets = np.array([_cut(image, o[3], o[4], BIG) for o in OBJECTS]) + seg_true = np.array([_cut(seg, o[3], o[4], 0) for o in OBJECTS]) + vignets[(seg_true != 0) & (seg_true != number[:, None, None])] = BIG + + paths = {name: str(tmp_path / f"{name}-001-001.fits") + for name in ("segmentation", "sexcat")} + fits.PrimaryHDU(seg).writeto(paths["segmentation"]) + objects = fits.BinTableHDU.from_columns([ + fits.Column(name="NUMBER", format="J", array=number), + fits.Column(name="X_IMAGE", format="E", array=x.astype(np.float32)), + fits.Column(name="Y_IMAGE", format="E", array=y.astype(np.float32)), + fits.Column(name="VIGNET", format=f"{STAMP * STAMP}E", + array=vignets.reshape(len(number), -1), + dim=f"({STAMP},{STAMP})"), + fits.Column(name="THETA_J2000", format="E", + array=np.zeros(len(number))), + fits.Column(name="X_IMAGE_DBL", format="D", array=x), + fits.Column(name="Y_IMAGE_DBL", format="D", array=y), + ], name="LDAC_OBJECTS") + imhead = fits.BinTableHDU.from_columns( + [fits.Column(name="Field Header Card", format="80A", + array=np.array(["HISTORY x"]))], name="LDAC_IMHEAD") + fits.HDUList([fits.PrimaryHDU(), imhead, objects]).writeto( + paths["sexcat"]) + return paths, seg_true + + +def _add(paths): + ss.add_seg_vignet(paths["sexcat"], paths["segmentation"]) + with fits.open(paths["sexcat"]) as hdul: + return [h.name for h in hdul], hdul["LDAC_OBJECTS"].data.copy() + + +def test_neither_rint_nor_float32_positions_give_the_centres(): + """The premise: np.rint of the double positions misses objects 3 and 4, + and SExtractor's rule on the float32 positions misses object 6.""" + x = np.array([o[1] for o in OBJECTS]) + y = np.array([o[2] for o in OBJECTS]) + col = np.array([o[3] for o in OBJECTS]) + row = np.array([o[4] for o in OBJECTS]) + rint_miss = (np.rint(x - 1) != col) | (np.rint(y - 1) != row) + assert np.flatnonzero(rint_miss).tolist() == [2, 3] + x32 = x.astype(np.float32).astype(np.float64) + y32 = y.astype(np.float32).astype(np.float64) + f32_miss = ((ss.vignet_centre(x32) != col) + | (ss.vignet_centre(y32) != row)) + assert np.flatnonzero(f32_miss).tolist() == [5] + + +@pytest.mark.parametrize("flat", [False, True], ids=["noisy", "flat"]) +def test_seg_vignet_is_registered_with_vignet(tmp_path, flat): + """SEG_VIGNET is the check image on VIGNET's grid for every object, the + half-pixel ties included, even where the image cannot tell the two + candidate centres apart: int32, VIGNET's shape, 0 off the image.""" + paths, seg_true = _scene(tmp_path, flat=flat) + names, data = _add(paths) + assert names == ["PRIMARY", "LDAC_IMHEAD", "LDAC_OBJECTS"] + seg_vignets = data["SEG_VIGNET"] + assert seg_vignets.dtype.kind == "i" and seg_vignets.dtype.itemsize == 4 + assert seg_vignets.shape == data["VIGNET"].shape + npt.assert_array_equal(seg_vignets, seg_true) + centre = STAMP // 2 + npt.assert_array_equal(seg_vignets[:, centre, centre], data["NUMBER"]) + + +def test_seg_vignet_replaces_the_double_positions(tmp_path): + """The written catalogue has the columns SExtractor writes without the + double positions, plus SEG_VIGNET, and every other HDU and column as + they were.""" + paths, _ = _scene(tmp_path) + with fits.open(paths["sexcat"]) as hdul: + before = hdul["LDAC_OBJECTS"].data.copy() + imhead = hdul["LDAC_IMHEAD"].data.copy() + _, after = _add(paths) + kept = [name for name in before.names if name not in DBL] + assert after.names == kept + ["SEG_VIGNET"] + for name in kept: + npt.assert_array_equal(after[name], before[name]) + with fits.open(paths["sexcat"]) as hdul: + npt.assert_array_equal(hdul["LDAC_IMHEAD"].data, imhead) + + +def test_seg_vignet_needs_the_double_positions(tmp_path): + paths, _ = _scene(tmp_path) + with fits.open(paths["sexcat"]) as hdul: + objects = hdul["LDAC_OBJECTS"] + cols = [c for c in objects.columns if c.name not in DBL] + hdus = [hdul[0].copy(), hdul[1].copy(), + fits.BinTableHDU.from_columns(cols, name="LDAC_OBJECTS")] + fits.HDUList(hdus).writeto(paths["sexcat"], overwrite=True) + with pytest.raises(ValueError, match="X_IMAGE_DBL"): + _add(paths) + + +def test_param_file_adds_the_double_positions(tmp_path): + """The runner's parameter file is the configured one plus the two + double-precision positions, appended after its last parameter.""" + dot_param = tmp_path / "default.param" + dot_param.write_text("NUMBER # running number\nX_IMAGE\nVIGNET(51,51)") + out = ss.seg_vignet_param_file(str(dot_param), str(tmp_path / "out.param")) + lines = [line.split("#")[0].strip() + for line in open(out).read().splitlines()] + assert [line for line in lines if line] == [ + "NUMBER", "X_IMAGE", "VIGNET(51,51)"] + DBL + assert dot_param.read_text().count("DBL") == 0 + + +def test_sextractor_caller_names_its_check_images(tmp_path): + """The runner finds the SEGMENTATION check image by type.""" + caller = ss.SExtractorCaller( + [str(tmp_path / "image-001-001.fits"), + str(tmp_path / "weight-001-001.fits")], + str(tmp_path), "-001-001", "d.sex", "d.param", "d.conv", + True, False, False, False, False, False, False, + check_image=["BACKGROUND", "SEGMENTATION"], + ) + assert caller.check_paths == { + "BACKGROUND": f"{tmp_path}/background-001-001.fits", + "SEGMENTATION": f"{tmp_path}/segmentation-001-001.fits", + } + + +# --- SEG_VIGNET through the join to the UNIONS catalogue ------------------- +# +# UberSeg (uberseg_weight, seg_has_neighbour, Ngmix._check_central_seg_label) +# compares each stamp's labels with the row's NUMBER. The check image is +# labelled with SExtractor's NUMBERs; the join replaces NUMBER with the +# catalogue's, and relabels every stamp through the whole SExtractor -> UNIONS +# map, the footprints of rows it drops included (UNMATCHED_LABEL). The scene +# makes old and new numbers collide: each paired row's new NUMBER is another +# object's SExtractor NUMBER, one of them the dropped detection's. + +CROWD_STAMP = 9 +# (SExtractor NUMBER, x, y, UNIONS NUMBER or None for no partner) +CROWD = [ + (1, 10.0, 10.0, 2), + (2, 13.0, 10.0, None), + (3, 10.0, 13.0, 1), + (4, 22.0, 22.0, 3), +] +EXT_ONLY = (7, 5.0, 25.0) +EXT_HEADER = ( + "# 1 NUMBER Running object number\n" + "# 2 X_IMAGE Object position along x [pixel]\n" + "# 3 Y_IMAGE Object position along y [pixel]\n" +) + + +def _crowd_seg(): + seg = np.zeros((30, 30), np.int32) + for number, x, y, _ in CROWD: + col, row = int(x) - 1, int(y) - 1 + seg[row - 1:row + 2, col - 1:col + 2] = number + return seg + + +def _write_crowd(sexcat_path, seg_path): + """SExtractor's outputs for CROWD: the sexcat (with the double positions + seg_vignet_param_file asks for) and the SEGMENTATION check image.""" + number = np.array([o[0] for o in CROWD], np.int32) + x = np.array([o[1] for o in CROWD]) + y = np.array([o[2] for o in CROWD]) + n = len(number) + size = CROWD_STAMP * CROWD_STAMP + objects = fits.BinTableHDU.from_columns([ + fits.Column(name="NUMBER", format="J", array=number), + fits.Column(name="X_IMAGE", format="E", array=x.astype(np.float32)), + fits.Column(name="Y_IMAGE", format="E", array=y.astype(np.float32)), + fits.Column(name="VIGNET", format=f"{size}E", + array=np.ones((n, size), np.float32), + dim=f"({CROWD_STAMP},{CROWD_STAMP})"), + fits.Column(name="X_IMAGE_DBL", format="D", array=x), + fits.Column(name="Y_IMAGE_DBL", format="D", array=y), + ], name="LDAC_OBJECTS") + imhead = fits.BinTableHDU.from_columns( + [fits.Column(name="Field Header Card", format="80A", + array=np.array(["HISTORY x"]))], name="LDAC_IMHEAD") + fits.HDUList([fits.PrimaryHDU(), imhead, objects]).writeto(sexcat_path) + fits.PrimaryHDU(_crowd_seg()).writeto(seg_path) + + +def _write_crowd_external(path): + rows = [(new, x, y) for _, x, y, new in CROWD if new is not None] + rows.append(EXT_ONLY) + path.write_text(EXT_HEADER + "".join( + f"{n:10d} {x:11.4f} {y:11.4f}\n" for n, x, y in rows)) + return str(path) + + +def _old_stamps(): + """The check image's stamps in SExtractor's numbering, row by row.""" + col = np.array([int(o[1]) - 1 for o in CROWD]) + row = np.array([int(o[2]) - 1 for o in CROWD]) + return ss.cut_stamps(_crowd_seg(), col, row, CROWD_STAMP, 0) + + +def _run_runner(tmp_path, monkeypatch, match): + """sextractor_runner as the workflow runs it under blend_handling: + uberseg, with SExtractor replaced by the CROWD outputs.""" + import logging + + from shapepipe.modules import sextractor_runner as runner_module + from shapepipe.pipeline.config import CustomParser + + out = tmp_path / "output" + tmp = tmp_path / "tmp" + out.mkdir() + tmp.mkdir() + num = "-001-001" + + def fake_execute(command_line): + _write_crowd(out / f"sexcat{num}.fits", out / f"segmentation{num}.fits") + return "", "All done" + + monkeypatch.setattr(runner_module, "execute", fake_execute) + dot_param = tmp_path / "default.param" + dot_param.write_text("NUMBER\nX_IMAGE\nY_IMAGE\nVIGNET(9,9)\n") + ext = _write_crowd_external(tmp_path / "CFIS_cat-001-001.cat") + config = CustomParser() + config["SEXTRACTOR_RUNNER"] = { + "EXEC_PATH": "source-extractor", + "DOT_SEX_FILE": "d.sex", "DOT_PARAM_FILE": str(dot_param), + "DOT_CONV_FILE": "d.conv", "WEIGHT_IMAGE": "True", + "FLAG_IMAGE": "False", "PSF_FILE": "False", + "DETECTION_IMAGE": "False", "DETECTION_WEIGHT": "False", + "ZP_FROM_HEADER": "False", "BKG_FROM_HEADER": "False", + "CHECKIMAGE": "BACKGROUND, SEGMENTATION", + "SEG_VIGNET": "${SP_SEG_VIGNET:-False}", + "MATCH_CATALOGUE": "${SP_MATCH_CATALOGUE:-}", + "MATCH_RADIUS": "1.0", "MATCH_MIN_FRACTION": "0.98", + "MATCH_TOLERATED_UNPAIRED": "20", "MAKE_POST_PROCESS": "False", + } + monkeypatch.setenv("SP_SEG_VIGNET", "True") + monkeypatch.setenv("SP_MATCH_CATALOGUE", ext if match else "") + runner_module.sextractor_runner( + [str(tmp_path / f"image{num}.fits"), + str(tmp_path / f"weight{num}.fits")], + {"output": str(out), "tmp": str(tmp)}, num, config, + "SEXTRACTOR_RUNNER", logging.getLogger("test"), + ) + with fits.open(out / f"sexcat{num}.fits") as hdul: + return hdul["LDAC_OBJECTS"].data.copy() + + +def test_runner_seg_vignet_carries_the_final_number(tmp_path, monkeypatch): + """Through the runner, each row's own footprint carries its final + (UNIONS) NUMBER and nothing else does, though every new NUMBER is some + other object's SExtractor NUMBER; the dropped detection's footprint is + UNMATCHED_LABEL. Fails if the join runs before add_seg_vignet, or relabels + only the paired labels.""" + from shapepipe.modules.sextractor_package import match_catalogue as mc + + data = _run_runner(tmp_path, monkeypatch, match=True) + paired = [i for i, o in enumerate(CROWD) if o[3] is not None] + old = _old_stamps()[paired] + old_number = np.array([CROWD[i][0] for i in paired]) + new_number = np.array([CROWD[i][3] for i in paired]) + npt.assert_array_equal(data["NUMBER"], new_number) + seg = data["SEG_VIGNET"] + assert "X_IMAGE_DBL" not in data.names + centre = CROWD_STAMP // 2 + npt.assert_array_equal(seg[:, centre, centre], data["NUMBER"]) + dropped = CROWD[1][0] + for i in range(len(paired)): + npt.assert_array_equal(seg[i] == new_number[i], + old[i] == old_number[i]) + npt.assert_array_equal(seg[i] == mc.UNMATCHED_LABEL, + old[i] == dropped) + npt.assert_array_equal(seg[i] == 0, old[i] == 0) + # Rows 0 and 1 overlap the dropped detection, whose SExtractor label (2) + # is row 0's new NUMBER. + assert (seg[0] == mc.UNMATCHED_LABEL).any() + assert (old[0] == new_number[0]).any() + + +def test_runner_seg_vignet_without_a_join_keeps_sextractor_numbers( + tmp_path, monkeypatch, +): + """Image simulations (empty MATCH_CATALOGUE): every row stays, and the + stamps are the check image's, labelled with SExtractor's NUMBER.""" + data = _run_runner(tmp_path, monkeypatch, match=False) + npt.assert_array_equal(data["NUMBER"], [o[0] for o in CROWD]) + npt.assert_array_equal(data["SEG_VIGNET"], _old_stamps()) + + +def test_uberseg_sees_the_same_neighbours_after_the_join(tmp_path, + monkeypatch): + """UberSeg only asks own-versus-other, so the relabelled stamps give the + weights and the neighbour flag the SExtractor-labelled ones give.""" + from shapepipe.modules.ngmix_package.ngmix import ( + seg_has_neighbour, + uberseg_weight, + ) + + data = _run_runner(tmp_path, monkeypatch, match=True) + paired = [i for i, o in enumerate(CROWD) if o[3] is not None] + old = _old_stamps()[paired] + weight = np.ones((CROWD_STAMP, CROWD_STAMP)) + for i, j in enumerate(paired): + new_seg, number = data["SEG_VIGNET"][i], data["NUMBER"][i] + npt.assert_array_equal( + uberseg_weight(weight, new_seg, number), + uberseg_weight(weight, old[i], CROWD[j][0])) + assert (seg_has_neighbour(new_seg, number) + == seg_has_neighbour(old[i], CROWD[j][0])) + assert seg_has_neighbour(data["SEG_VIGNET"][0], data["NUMBER"][0]) + assert not seg_has_neighbour(data["SEG_VIGNET"][2], data["NUMBER"][2]) + + +def test_a_partial_relabel_would_collide(): + """The negative control: relabelling only the paired labels leaves the + dropped detection's SExtractor label 2 in row 0's stamp, equal to row 0's + new NUMBER, so UberSeg would take the neighbour for the object.""" + old = _old_stamps()[0] + partial = old.copy() + for number, _, _, new in CROWD: + if new is not None: + partial[old == number] = new + assert ((partial == 2) & (old != 1)).any() + + +def _join_with_numbers(case, numbers): + """Join CROWD to a catalogue whose three objects carry ``numbers`` + (written as given), and return the joined NUMBER column.""" + from shapepipe.modules.sextractor_package import match_catalogue as mc + + case.mkdir() + _write_crowd(case / "sexcat.fits", case / "seg.fits") + rows = zip(numbers, [10.0, 10.0, 22.0], [10.0, 13.0, 22.0]) + ext = case / "ext.cat" + ext.write_text(EXT_HEADER + "".join( + f"{n:>10} {x:11.4f} {y:11.4f}\n" for n, x, y in rows)) + mc.match_catalogue(str(case / "sexcat.fits"), str(ext)) + with fits.open(case / "sexcat.fits") as hdul: + return hdul["LDAC_OBJECTS"].data["NUMBER"].tolist() + + +@pytest.mark.parametrize("numbers", [ + [5, 5, 6], # repeated + [0, 5, 6], # not positive + [-3, 5, 6], + ["5.5", 7, 6], # not an integer + ["5.0", 7, 6], # read as a float + [2**31, 5, 6], # past int32 + [2**32 - 1, 5, 6], # would narrow to -1, UNMATCHED_LABEL +]) +def test_join_rejects_numbers_that_could_collide(tmp_path, numbers): + """Unique catalogue NUMBERs that the int32 columns hold exactly are what + keep a relabelled stamp's own footprint apart from every other one.""" + with pytest.raises(ValueError, match="NUMBER"): + _join_with_numbers(tmp_path / "case", numbers) + + +def test_join_accepts_the_largest_int32_number(tmp_path): + from shapepipe.modules.sextractor_package import match_catalogue as mc + + assert _join_with_numbers(tmp_path / "case", [1, mc.MAX_NUMBER, 6]) == [ + 1, mc.MAX_NUMBER, 6] + + +def test_seg_vignet_stays_out_of_the_final_catalogue(tmp_path): + """make_cat drops SEG_VIGNET with VIGNET: the stamps are ngmix inputs, + not catalogue columns.""" + from shapepipe.modules.make_cat_package import make_cat + + _write_crowd(tmp_path / "sexcat-001-001.fits", tmp_path / "seg.fits") + ss.add_seg_vignet(str(tmp_path / "sexcat-001-001.fits"), + str(tmp_path / "seg.fits")) + final = make_cat.prepare_final_cat_file(str(tmp_path), "-001-001") + make_cat.save_sextractor_data(final, str(tmp_path / "sexcat-001-001.fits")) + with fits.open(tmp_path / "final_cat-001-001.fits") as hdul: + names = hdul["RESULTS"].data.names + assert "VIGNET" not in names and "SEG_VIGNET" not in names + assert "NUMBER" in names diff --git a/tests/science/test_defect_interpolation.py b/tests/science/test_defect_interpolation.py new file mode 100644 index 000000000..d36be39b0 --- /dev/null +++ b/tests/science/test_defect_interpolation.py @@ -0,0 +1,119 @@ +"""Shear recovery for the defects kept under ``DEFECT_FILL = interpolate``. + +Physics invariant: every defect the veto keeps leaves both additive terms +|c1|, |c2| < 5e-4 and both diagonal multiplicative terms |m11|, |m22| < 1%, +from the full 2x2 response matrix. Interpolated defects (columns, full and +finite 3-px bleeds, single pixels) are checked at +``EPOCH_INTERPOLATED_DEFECT_RADIUS`` and two pixels beyond it, on galaxies +with half-light radius 0.3" and 0.5" through a 0.7" PSF, round and with +ellipticity (0.05, 0.02), and on a 0.7" galaxy through a 0.9" PSF. Defects +too wide to interpolate are noise-filled, and are checked at +``EPOCH_CENTRAL_DEFECT_RADIUS``. The cases sit at the radii themselves, so +lowering either below its calibrated value turns this red. + +Positive control: a 3-px bleed three pixels inside the interpolated-defect +radius breaks the bound. +""" + +import json + +import numpy as np +import pytest + +from shapepipe.modules.ngmix_package.defect_interpolation import ( + interpolable_defects, +) +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_INTERPOLATED_DEFECT_RADIUS, + EPOCH_MASKED_FRACTION_CUT, + central_defect_vetoes, + defect_mask, +) +from tests.helpers.defect_response import defect_response + +N = 51 +CENTRE = N // 2 +RI = int(np.ceil(EPOCH_INTERPOLATED_DEFECT_RADIUS)) +RN = int(np.ceil(EPOCH_CENTRAL_DEFECT_RADIUS)) +ROUND = (0.0, 0.0) +ELLIPTICAL = (0.05, 0.02) +SEEDS = range(6) +INTERPOLATE = {"defect_fill": "interpolate"} + + +def geometry(kind, distance): + """A detector defect whose nearest pixel is ``distance`` px from the + stamp centre.""" + bad = np.zeros((N, N), dtype=bool) + near = CENTRE + distance + if kind == "pixel": + bad[CENTRE, near] = True + elif kind == "column": + bad[:, near] = True + elif kind == "bleed": + bad[:, near:near + 3] = True + elif kind == "finite_bleed": + bad[CENTRE - 5:CENTRE + 6, near:near + 3] = True + elif kind == "wide_bleed": + bad[:, near:near + 5] = True + elif kind == "edge": + bad[:, near:] = True + return bad + + +INTERPOLATED = ("column", "bleed", "finite_bleed", "pixel") +CASES = ( + [(k, d, h, 0.7, ROUND) for k in INTERPOLATED + for d in (RI, RI + 2) for h in (0.3, 0.5)] + + [(k, RI, 0.5, 0.7, ELLIPTICAL) for k in INTERPOLATED] + + [(k, RI, 0.7, 0.9, psf) for k in ("bleed", "finite_bleed") + for psf in (ROUND, ELLIPTICAL)] + + [(k, RN, h, 0.7, ROUND) for k in ("wide_bleed", "edge") + for h in (0.3, 0.5)] +) + + +def case_id(case): + kind, distance, hlr, fwhm, psf_shear = case + psf = "elliptical" if any(psf_shear) else "round" + return f"{kind}-{distance}px-hlr{hlr}-psf{fwhm}-{psf}" + + +def recover(bad, hlr, fwhm, psf_shear, tmp_path): + result = defect_response(bad, hlr=hlr, psf=fwhm, seeds=SEEDS, + psf_shear=psf_shear, options=INTERPOLATE) + (tmp_path / "recovery.json").write_text(json.dumps(result, indent=2)) + return np.abs(result["m"]).max(), np.abs(result["c"]).max(), result + + +@pytest.mark.parametrize("kind,distance,hlr,fwhm,psf_shear", CASES, + ids=[case_id(c) for c in CASES]) +def test_kept_defects_recover_shear_on_both_axes(kind, distance, hlr, fwhm, + psf_shear, tmp_path): + """Failure modes: a veto radius is below its calibrated value; the + interpolated pixels' quarter-turn orbit keeps its weight (a one-sided + hole in the likelihood); the fill mask is symmetrized; a wide hole is + interpolated; raw defect values leak into metacal.""" + bad = geometry(kind, distance) + masked = defect_mask(np.ones((N, N)), bad.astype(np.int32)) + interpolated = interpolable_defects(masked) + assert interpolated.any() == (kind in INTERPOLATED) + assert masked.mean() <= EPOCH_MASKED_FRACTION_CUT + assert not central_defect_vetoes(masked, EPOCH_CENTRAL_DEFECT_RADIUS, + "interpolate") + m, c, result = recover(bad, hlr, fwhm, psf_shear, tmp_path) + assert m < 0.01, result + assert c < 5e-4, result + + +def test_vetoed_bleed_breaks_the_bound(tmp_path): + """Positive control: a 3-px bleed three pixels inside the + interpolated-defect radius, on the 0.3" galaxy, gives |c| > 1e-3 and + |m| > 1.5%. The veto drops it.""" + bad = geometry("bleed", RI - 3) + assert central_defect_vetoes(bad, EPOCH_CENTRAL_DEFECT_RADIUS, + "interpolate") + m, c, result = recover(bad, 0.3, 0.7, ROUND, tmp_path) + assert c > 1e-3, result + assert m > 0.015, result diff --git a/tests/science/test_defect_veto.py b/tests/science/test_defect_veto.py new file mode 100644 index 000000000..3384263e4 --- /dev/null +++ b/tests/science/test_defect_veto.py @@ -0,0 +1,104 @@ +"""Shear recovery for the defects the central-defect veto keeps (noise fill). + +Physics invariant: a defect the veto keeps -- a column, a 3-px bleed or a +single pixel at the veto radius or beyond, or an edge band -- leaves both +additive terms |c1|, |c2| < 5e-4 and both diagonal multiplicative terms +|m11|, |m22| < 1%, from the full 2x2 response matrix. The grid spans galaxies +with half-light radius 0.3" and 0.5" through a 0.7" PSF, round and with +ellipticity (0.05, 0.02). A noise-filled defect biases m anisotropically +(m11 and m22 differ by up to a factor of ten), so a scalar m would hide it. +The cases sit at ``EPOCH_CENTRAL_DEFECT_RADIUS`` itself, so lowering the +radius below the calibrated value turns this red. + +Positive control: the same column two pixels inside the radius breaks the +bound, so the recovery check can see the bias the veto removes. +""" + +import json + +import numpy as np +import pytest + +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_MASKED_FRACTION_CUT, + defect_mask, + has_central_defect, +) +from tests.helpers.defect_response import defect_response + +N = 51 +CENTRE = N // 2 +R = int(np.ceil(EPOCH_CENTRAL_DEFECT_RADIUS)) +ELLIPTICAL = (0.05, 0.02) +SEEDS = range(6) + + +def geometry(kind, distance): + """A detector defect whose nearest pixel is ``distance`` px from the + stamp centre; an edge band reaches in to that distance.""" + bad = np.zeros((N, N), dtype=bool) + if kind == "pixel": + bad[CENTRE, CENTRE + distance] = True + elif kind == "column": + bad[:, CENTRE + distance] = True + elif kind == "bleed": + bad[:, CENTRE + distance:CENTRE + distance + 3] = True + elif kind == "edge": + bad[:, CENTRE + distance:] = True + return bad + + +# Wide defects (a 3-px bleed, the widest edge band the veto keeps) at the +# radius itself on the 0.5" galaxy through the elliptical PSF sit at the +# bound (m11 = -0.98% +/- 0.04% and -0.94% +/- 0.06%; see +# :func:`has_central_defect`), so those two are checked one pixel out. +CASES = ( + [(k, d, h, (0.0, 0.0)) for k in ("column", "bleed", "pixel") + for d in (R, R + 2) for h in (0.3, 0.5)] + + [(k, R, 0.5, ELLIPTICAL) for k in ("column", "pixel")] + + [(k, R + 1, 0.5, ELLIPTICAL) for k in ("bleed", "edge")] + + [("edge", R, 0.5, (0.0, 0.0))] + + [("edge", N - CENTRE - 5, 0.5, psf) for psf in ((0.0, 0.0), ELLIPTICAL)] +) + + +def recover(bad, hlr, psf_shear, tmp_path): + result = defect_response(bad, hlr=hlr, psf=0.7, seeds=SEEDS, + psf_shear=psf_shear) + (tmp_path / "recovery.json").write_text(json.dumps(result, indent=2)) + return np.abs(result["m"]).max(), np.abs(result["c"]).max(), result + + +def case_id(case): + kind, distance, hlr, psf_shear = case + psf = "elliptical" if any(psf_shear) else "round" + return f"{kind}-{distance}px-hlr{hlr}-{psf}" + + +@pytest.mark.parametrize("kind,distance,hlr,psf_shear", CASES, + ids=[case_id(c) for c in CASES]) +def test_kept_defects_recover_shear_on_both_axes(kind, distance, hlr, + psf_shear, tmp_path): + """Failure modes: the veto radius is below the calibrated value; the + fill leaks raw defect values or leaves more of the object's light + missing; the check reads m11 alone and misses m22, or c1 alone. + """ + bad = geometry(kind, distance) + masked = defect_mask(np.ones((N, N)), bad.astype(np.int32)) + np.testing.assert_array_equal(masked, bad) + assert masked.mean() <= EPOCH_MASKED_FRACTION_CUT + assert not has_central_defect(masked, EPOCH_CENTRAL_DEFECT_RADIUS) + m, c, result = recover(bad, hlr, psf_shear, tmp_path) + assert m < 0.01, result + assert c < 5e-4, result + + +def test_vetoed_column_breaks_the_bound(tmp_path): + """Positive control: a column two pixels inside the radius, on the + 0.5" galaxy, gives |m| > 2% and |c| > 1e-3. The veto drops it.""" + bad = geometry("column", R - 2) + assert has_central_defect(bad, EPOCH_CENTRAL_DEFECT_RADIUS) + m, c, result = recover(bad, 0.5, (0.0, 0.0), tmp_path) + assert m > 0.02, result + assert c > 1e-3, result diff --git a/tests/unit/test_workflow_tile_detection.py b/tests/unit/test_workflow_tile_detection.py index 498a015de..ae75aac5c 100644 --- a/tests/unit/test_workflow_tile_detection.py +++ b/tests/unit/test_workflow_tile_detection.py @@ -153,3 +153,28 @@ def test_catalogue_retrieve_follows_the_prefix(source, retrieve): config = {"tile_detection": "unions_catalogue", "inputs": {"catalogues": source}} assert run_config.catalogue_source(config) == (source, retrieve) + + +# --- blend handling -------------------------------------------------------- + + +def test_blend_handlings_mirror_the_ngmix_module(): + """run_config.BLEND_HANDLINGS is a copy; the copy must stay true. + + The Snakefile validates `blend_handling:` against it in the launcher venv, + outside the container where shapepipe is importable, so the tuple is + mirrored rather than imported, and read out of the source text here for + the same reason: this file is container-free. + """ + import re + + src = (REPO_ROOT / "src" / "shapepipe" / "modules" / "ngmix_package" + / "ngmix.py").read_text() + match = re.search(r"^BLEND_HANDLINGS = \(([^)]*)\)", src, re.M) + assert match, "ngmix.py no longer defines BLEND_HANDLINGS" + assert _load("run_config").BLEND_HANDLINGS == tuple( + part.strip().strip('"\'') for part in match.group(1).split(",") + if part.strip()) + committed = _load("run_config").load( + str(REPO_ROOT / "workflow" / "config.yaml")) + assert committed["blend_handling"] == "uberseg" diff --git a/tests/workflow/README.md b/tests/workflow/README.md index c3576f08c..9748e9126 100644 --- a/tests/workflow/README.md +++ b/tests/workflow/README.md @@ -44,7 +44,7 @@ The API context remains open while tests inspect jobs and closes before the fixt ## Campaign-boundary pin -`params_pin.json` pins SHA-256 digests of: +`params_pin.json` pins the default campaign (data, psfex, `blend_handling: uberseg`), and `params_pin_noisefill.json` the same campaign with `blend_handling: noisefill`, so a campaign that ran as noisefill stays resumable. Each pins SHA-256 digests of: - `unit_pre()` rendered for every stage; - every rule's shell template; @@ -76,7 +76,10 @@ Apply mutations only to a disposable checkout, run the named test without `--upd | `test_clean_exposure_waits_on_persist_iff_psf` | Drop the persist edge; make it unconditional under fake PSFs; drop vignets consumers; remove the in-scope consumer filter. | | `test_final_cat_merge_reads_every_ready_tile` | Drop one ready tile; append an out-of-scope tile. | | `test_products_use_products_dir_and_run_name` | Rename either merged catalogue or the persist manifest; route products to scratch; derive `CAMPAIGN` from the products directory's basename. | +| `test_blend_handling_reaches_detection_and_ngmix_only_under_uberseg` | Make the Snakefile or config.yaml default `noisefill`; export `BLEND_ENV` under noisefill; drop `blend_env(...)` from either tile_detect or tile_ngmix; export the wrong value; set `NGMIX_SEG_MEM_MB` or `DETECT_SEG_MEM_MB` to 0; give a committed ini's `SEG_VIGNET` / `BLEND_HANDLING` a literal or the wrong default. | +| `test_unknown_blend_handling_fails_during_parse` | Remove the Snakefile's `blend_handling` check. | | `test_missing_run_fails_during_parse` | Remove `run` from `run_config.REQUIRED`; literal paths must still receive the required-key diagnostic, not a later `KeyError`. | | `test_unit_pre_changes_at_campaign_boundary` | Append a line to `unit_pre`; change one rule's `params.pre`; change one shell; change a rendered thread count. | +| `test_noisefill_plan_changes_at_campaign_boundary` | Export anything under noisefill; change one rule's `params.pre` or shell. | | `test_params_pin_ignores_fixture_root` | Remove fixture-root normalization. | | `test_params_is_a_rerun_trigger_under_both_profiles` | Remove `params` from candide or nibi's `rerun-triggers`. | diff --git a/tests/workflow/conftest.py b/tests/workflow/conftest.py index 134a9c62a..f08beaaa8 100644 --- a/tests/workflow/conftest.py +++ b/tests/workflow/conftest.py @@ -22,7 +22,7 @@ def campaign(request, tmp_path): @pytest.fixture def resolve_dag(monkeypatch): """Expose the resolver so tests can also assert parse-time failures.""" - return lambda campaign: resolve(campaign, monkeypatch) + return lambda campaign, **kwargs: resolve(campaign, monkeypatch, **kwargs) @pytest.fixture diff --git a/tests/workflow/harness.py b/tests/workflow/harness.py index 1c8420f26..abe0783ad 100644 --- a/tests/workflow/harness.py +++ b/tests/workflow/harness.py @@ -206,8 +206,12 @@ def load_profile(name): @contextmanager -def resolve(campaign, monkeypatch): - """Resolve ``all`` in an isolated state directory, without running jobs.""" +def resolve(campaign, monkeypatch, launch_env=None): + """Resolve ``all`` in an isolated state directory, without running jobs. + + ``launch_env`` adds variables to the scrubbed environment, as if the + launching shell had exported them. + """ from snakemake import workflow as sm_workflow scripts = REPO / "workflow" / "scripts" @@ -234,6 +238,8 @@ def resolve(campaign, monkeypatch): "XDG_CACHE_HOME": campaign.root / "cache", }.items(): patch.setenv(name, str(value)) + for name, value in (launch_env or {}).items(): + patch.setenv(name, str(value)) profile = load_profile("candide") try: with SnakemakeApi() as api: diff --git a/tests/workflow/params_pin.json b/tests/workflow/params_pin.json index 6de487277..1744bb06d 100644 --- a/tests/workflow/params_pin.json +++ b/tests/workflow/params_pin.json @@ -11,7 +11,7 @@ "final_cat_merge": "e7f46859c4503a2220713d7bb2507555515d0a9632d780b20f14c59e32210023", "prepare_all_tiles": "b8f872a22adf014e25a7fa5198f49b71a6fe9e56042ed82b682bc8763970a844", "star_cat_merge": "6277450958474af5270982fa35360f2f237a29f7533c526ee9265dfd5acc07a0", - "tile_detect": "8336b148769e9d43b64f0d945c05e7d163c64dba4e6287e8c1964d1502e3ffe1", + "tile_detect": "24d015063fb07ba4a48b961284fa1345675f13e167acc62a2c7b397795113b82", "tile_exp_forest": "7447ab4a1049de5f0b5c81e5f9ed2a8644c7bdab85cb0a060bfde81261a89b28", "tile_find_exposures": "8704317871744996c44351c2836fcb222d7986a602e9e046a90d684ba7b3c184", "tile_get_catalogue": "7e3f889a955a14b0b917a015c2a8bc90e433a89b1cc5e00b8ca94695c4c93d3c", @@ -19,12 +19,12 @@ "tile_make_cat": "a57518b04c11f70bf41b320532fddffdd1115fdac604d7dd3928eb61cca28f23", "tile_merge_cats": "ff21216ea804dccc2d2c290d2b2499d5d05f0c34c0a56993c233f43fe3c06bdb", "tile_merge_headers": "7a344849d62936e2f5598dc8731a2c4947eff2c4c7218b1a731dd2a9577e7111", - "tile_ngmix": "6ed10a3a4d3658ba0303fb0c23ad5100ec8cbdc638ca36c3b25875087b5c3415", + "tile_ngmix": "5b86694f3bac7b8c3f62d6c1df72b38a8b9173fb69ab53c75d7d70e38f5dc9ff", "tile_uncompress": "1e2b01acbf9708e0371070fb01c5b9568d7efb5f89d1835fc6bf91e2c8b60cb3", "tile_vignets": "9d4ae0d99c18217f2185f245281312454c8a219ec1628176e08c71a5efc4dc91" }, "schema": 1, - "sha256": "d67ee5da728e8fee513e1f01a522b386eb6a86d21c9a36250afd6013a06301d5", + "sha256": "7d03fd3acb7608ed1eda69b162f46b2c88eddb8bfa441cbf66245a65510d809b", "unit_pre": { "exp_get_images": "8dec850af212879f225fcf27a5f1281e1a075264c7b97d38c2214395d360168c", "exp_psf": "f2358ddf7385918dc5033d10b37f6dc97a15d02b071a3ea0a4619a5f7e6f5bec", diff --git a/tests/workflow/params_pin_noisefill.json b/tests/workflow/params_pin_noisefill.json new file mode 100644 index 000000000..6de487277 --- /dev/null +++ b/tests/workflow/params_pin_noisefill.json @@ -0,0 +1,43 @@ +{ + "algorithm": "sha256", + "rules": { + "all": "572122d8d1901e12ff591b30adf405f8920be8183129befececbb819f3392ed8", + "clean_exposure": "22cb76b13a5205d20a02a9bd3b8c8bea5ea2801555e24dd2b84b11145f7e79d9", + "clean_tile": "a5c07b0461526ed407df36a291deb866d4181524b4c8e0fbc3dd047fd9d28479", + "exp_get_images": "71e76ff7f96af1c5d2b85697cc5819e2271911f177a2075253f3b8cfa268c1a9", + "exp_persist": "302e2837542bc1102430c27c81c600b7cda32e8bddcb5fd60d33950987609fff", + "exp_psf": "2c4f6d00f1939ccbf05b4982a0202a0ff92a4373a727f4e00f3aaab4eba03352", + "exp_split": "6e954f8f3d06f44d3f164675912ce27d6216648855d04968f9168bd7d0f2c4fa", + "final_cat_merge": "e7f46859c4503a2220713d7bb2507555515d0a9632d780b20f14c59e32210023", + "prepare_all_tiles": "b8f872a22adf014e25a7fa5198f49b71a6fe9e56042ed82b682bc8763970a844", + "star_cat_merge": "6277450958474af5270982fa35360f2f237a29f7533c526ee9265dfd5acc07a0", + "tile_detect": "8336b148769e9d43b64f0d945c05e7d163c64dba4e6287e8c1964d1502e3ffe1", + "tile_exp_forest": "7447ab4a1049de5f0b5c81e5f9ed2a8644c7bdab85cb0a060bfde81261a89b28", + "tile_find_exposures": "8704317871744996c44351c2836fcb222d7986a602e9e046a90d684ba7b3c184", + "tile_get_catalogue": "7e3f889a955a14b0b917a015c2a8bc90e433a89b1cc5e00b8ca94695c4c93d3c", + "tile_get_images": "331a67e747f211ebf4c14b946a7af7f9fc9f55243d69d2791d74aecc3ca228c3", + "tile_make_cat": "a57518b04c11f70bf41b320532fddffdd1115fdac604d7dd3928eb61cca28f23", + "tile_merge_cats": "ff21216ea804dccc2d2c290d2b2499d5d05f0c34c0a56993c233f43fe3c06bdb", + "tile_merge_headers": "7a344849d62936e2f5598dc8731a2c4947eff2c4c7218b1a731dd2a9577e7111", + "tile_ngmix": "6ed10a3a4d3658ba0303fb0c23ad5100ec8cbdc638ca36c3b25875087b5c3415", + "tile_uncompress": "1e2b01acbf9708e0371070fb01c5b9568d7efb5f89d1835fc6bf91e2c8b60cb3", + "tile_vignets": "9d4ae0d99c18217f2185f245281312454c8a219ec1628176e08c71a5efc4dc91" + }, + "schema": 1, + "sha256": "d67ee5da728e8fee513e1f01a522b386eb6a86d21c9a36250afd6013a06301d5", + "unit_pre": { + "exp_get_images": "8dec850af212879f225fcf27a5f1281e1a075264c7b97d38c2214395d360168c", + "exp_psf": "f2358ddf7385918dc5033d10b37f6dc97a15d02b071a3ea0a4619a5f7e6f5bec", + "exp_split": "358fa8bbe59680d4f9839e007b343dd25fed733cf3156dccd30d305d33ac9480", + "tile_detect": "adad5d671fa65dd04433e1c82a725b845635e9d2e368a1e70833ed990df2d55a", + "tile_find_exposures": "c5922fb507f6fd040a179b53fba0818697661016c9688984b6ce249186dc6986", + "tile_get_catalogue": "b124d12252617abb8df8a98d6234ddfd6cc503e19fa72c9a0c48c21559d6f48c", + "tile_get_images": "45b47c44bfeb34ac973b89d4e752c8c82e78028f9f0da379c44627005b45e279", + "tile_make_cat": "7579a52e32e76523c0b77d75c468a47d7f2e2ca0c55865f40e9628a7481cd79d", + "tile_merge_cats": "69cfe94ba941d2ba7b1ce24883961b47b191b80aff381da1865c6bc38de0c354", + "tile_merge_headers": "8590cfa3281c88d43eb8c4fe5760176a13e928a419cbea9bea3a8f88b0634b42", + "tile_ngmix": "b6b75b553a62df54aea6c337dcde1e0cdcd0fe30000e86880c1a4689e42c4a72", + "tile_uncompress": "91a15538491e53ee2b0d52b0472909374270918e4b4c80880c4f5a8329e40161", + "tile_vignets": "a6111cef708aa3c9fb0144b49de9f781eef84d2096ba4f8a3de9bb45b1a3774d" + } +} diff --git a/tests/workflow/test_dag.py b/tests/workflow/test_dag.py index dd5fc7f70..6a8ccafd3 100644 --- a/tests/workflow/test_dag.py +++ b/tests/workflow/test_dag.py @@ -1,5 +1,6 @@ """Resolved-job checks for campaign scope, product paths, and PSF custody.""" +import os from collections import Counter from pathlib import Path @@ -137,3 +138,155 @@ def test_mccd_is_refused_during_parse(tmp_path, resolve_dag): with pytest.raises(WorkflowError, match=r"psf_model=mccd: PSF persistence"): with resolve_dag(campaign): pytest.fail("psf_model=mccd must be refused at parse time") + + +# --- blend_handling --------------------------------------------------------- + +# The memory each stage adds under uberseg for SEG_VIGNET. +SEG_MEM_MB = {"tile_detect": 3000, "tile_ngmix": 500} +BLEND_EXPORTS = {"tile_detect": {"SP_SEG_VIGNET": "True"}, + "tile_ngmix": {"SP_BLEND_HANDLING": "uberseg"}} +# The option each stage's committed ini reads the export through, and the +# value it must resolve to with and without it. +BLEND_OPTIONS = {"tile_detect": ("SEG_VIGNET", "False", "True"), + "tile_ngmix": ("BLEND_HANDLING", "noisefill", "uberseg")} + + +def _surface(campaign, resolve_dag): + """Every job's prologue, shell and memory, keyed by rule and wildcards.""" + with resolve_dag(campaign) as dag: + return dag.rule_names, { + (job.rule.name, tuple(sorted(job.wildcards_dict.items()))): ( + getattr(job.params, "pre", None), job.shellcmd, + job.resources.get("mem_mb")) + for job in dag.jobs + } + + +def _committed_value(shell, config_dir, option, env, monkeypatch): + """``option`` of the module section of the ini ``shell`` runs, expanded + under ``env`` as ShapePipe expands it.""" + from shapepipe.pipeline.config import CustomParser + + name = shell.split('shapepipe_run -c "$SP_CONFIG/')[1].split('"')[0] + parser = CustomParser() + parser.optionxform = str + assert parser.read(config_dir / name) + section = next(s for s in parser.sections() if s.endswith("_RUNNER")) + for key in ("SP_SEG_VIGNET", "SP_BLEND_HANDLING"): + monkeypatch.delenv(key, raising=False) + for key, value in env.items(): + monkeypatch.setenv(key, value) + return parser.getexpanded(section, option) + + +@pytest.mark.parametrize("detection", ["sextractor", "unions_catalogue"]) +def test_blend_handling_reaches_detection_and_ngmix_only_under_uberseg( + tmp_path, resolve_dag, monkeypatch, detection): + """A campaign without the knob plans exactly the uberseg campaign. Against + explicit noisefill, uberseg only adds its two exports to tile_detect's and + tile_ngmix's prologues and the seg stamps' memory to both, and each + export turns its committed ini's option from the noise-fill default to + the uberseg value.""" + surfaces = {} + for blend in (None, "noisefill", "uberseg"): + campaign = Campaign(tmp_path / str(blend), "data", "psfex") + campaign.config["tile_detection"] = detection + if blend is not None: + campaign.config["blend_handling"] = blend + campaign.write_config() + surfaces[blend] = _surface(campaign, resolve_dag) + + def normalized(blend): + rules, jobs = surfaces[blend] + root = str(tmp_path / str(blend)) + return rules, {key: tuple(v.replace(root, "") + if isinstance(v, str) else v for v in value) + for key, value in jobs.items()} + + assert normalized("uberseg") == normalized(None) + rules, noisefill = normalized("noisefill") + uberseg_rules, uberseg = normalized("uberseg") + assert uberseg_rules == rules + assert uberseg.keys() == noisefill.keys() + config_dir = Path(__file__).parents[2] / "workflow" / "config" / "cfis" + for key, (pre, shell, mem) in uberseg.items(): + rule = key[0] + nf_pre, nf_shell, nf_mem = noisefill[key] + exports = BLEND_EXPORTS.get(rule, {}) + lines = {f"export {name}='{value}'" for name, value in exports.items()} + + def without_exports(text): + if text is None: + return None + return "\n".join(line for line in text.split("\n") + if line not in lines) + + # The rendered shell carries the prologue, so it moves with it and + # nowhere else. + assert without_exports(shell) == nf_shell, key + if pre is None: + assert nf_pre is None and not exports, key + continue + assert lines <= set(pre.split("\n")), key + assert without_exports(pre) == nf_pre, key + if exports: + option, default, value = BLEND_OPTIONS[rule] + assert _committed_value(shell, config_dir, option, {}, + monkeypatch) == default + assert _committed_value(shell, config_dir, option, exports, + monkeypatch) == value + assert mem == nf_mem + SEG_MEM_MB.get(rule, 0), key + + +def test_unknown_blend_handling_fails_during_parse(tmp_path, resolve_dag): + campaign = Campaign(tmp_path / "campaign", "data", "psfex") + campaign.config["blend_handling"] = "mof" + campaign.write_config() + with pytest.raises(WorkflowError, match=r"Invalid blend_handling='mof'"): + with resolve_dag(campaign): + pytest.fail("an unknown blend_handling must fail at parse time") + + +INHERITED_BLEND_ENV = { + f"{prefix}{name}": value + for prefix in ("", "APPTAINERENV_", "SINGULARITYENV_") + for name, value in (("SP_SEG_VIGNET", "True"), + ("SP_BLEND_HANDLING", "uberseg")) +} + + +def test_noisefill_ignores_blend_variables_in_the_launch_shell( + tmp_path, resolve_dag): + """A noisefill campaign launched from a shell that still exports the + uberseg variables plans exactly the campaign launched from a clean one, + and the parse leaves none of them in the environment jobs inherit (the + slurm executor submits with --export=ALL from this process).""" + campaigns = {} + for name in ("clean", "dirty"): + campaigns[name] = Campaign(tmp_path / name, "data", "psfex") + campaigns[name].config["blend_handling"] = "noisefill" + campaigns[name].write_config() + clean, dirty = campaigns["clean"], campaigns["dirty"] + _, clean_jobs = _surface(clean, resolve_dag) + + with resolve_dag(dirty, launch_env=INHERITED_BLEND_ENV) as dag: + leaked = sorted(set(INHERITED_BLEND_ENV) & set(os.environ)) + dirty_jobs = { + (job.rule.name, tuple(sorted(job.wildcards_dict.items()))): ( + getattr(job.params, "pre", None), job.shellcmd, + job.resources.get("mem_mb")) + for job in dag.jobs + } + assert leaked == [] + + def normalized(jobs, root): + return {key: tuple(v.replace(str(root), "") + if isinstance(v, str) else v for v in value) + for key, value in jobs.items()} + + assert (normalized(dirty_jobs, tmp_path / "dirty") + == normalized(clean_jobs, tmp_path / "clean")) + for _, shell, _ in dirty_jobs.values(): + assert "SP_BLEND_HANDLING" not in (shell or "") + assert "SP_SEG_VIGNET" not in (shell or "") diff --git a/tests/workflow/test_params_pin.py b/tests/workflow/test_params_pin.py index faeb43158..c97a8177b 100644 --- a/tests/workflow/test_params_pin.py +++ b/tests/workflow/test_params_pin.py @@ -7,6 +7,9 @@ from tests.workflow.params import params_pin PIN = Path(__file__).with_name("params_pin.json") +# blend_handling: noisefill's plan, pinned on its own so a campaign that ran as +# noisefill stays resumable with the knob set explicitly. +NOISEFILL_PIN = Path(__file__).with_name("params_pin_noisefill.json") BOUNDARY_MESSAGE = ( "params.pre is a rerun trigger under both profiles; this change reruns " "every finished unit of a resumed campaign — land it at a campaign " @@ -23,6 +26,21 @@ def test_unit_pre_changes_at_campaign_boundary(psfex_dag, pytestconfig): assert actual == json.loads(PIN.read_text()), BOUNDARY_MESSAGE +def test_noisefill_plan_changes_at_campaign_boundary(tmp_path, resolve_dag, + pytestconfig): + """The same pin for a campaign that sets blend_handling: noisefill.""" + campaign = Campaign(tmp_path / "campaign", "data", "psfex") + campaign.config["blend_handling"] = "noisefill" + campaign.write_config() + with resolve_dag(campaign) as dag: + actual = params_pin(dag) + if pytestconfig.getoption("--update-params-pin"): + NOISEFILL_PIN.write_text( + json.dumps(actual, indent=2, sort_keys=True) + "\n") + assert NOISEFILL_PIN.is_file(), BOUNDARY_MESSAGE + assert actual == json.loads(NOISEFILL_PIN.read_text()), BOUNDARY_MESSAGE + + def test_params_pin_ignores_fixture_root(tmp_path, resolve_dag): """A different temporary campaign directory cannot require a new pin.""" pins = [] diff --git a/universes/committed.yaml b/universes/committed.yaml index 05d88e097..fa4893a84 100644 --- a/universes/committed.yaml +++ b/universes/committed.yaml @@ -53,9 +53,9 @@ analyses: psf_likelihood_noise: psf_noise_1em5 megacam_ccd_flip: megapipe_flip defect_fill: noise - blend_handling: none + blend_handling: uberseg epoch_masked_fraction_cut: one_third - central_defect_veto: disabled + central_defect_veto: fixed_radii catalogue_assembly: decisions: star_galaxy_classification: deferred_downstream diff --git a/workflow/README.md b/workflow/README.md index 4f0b135c3..19e7ae683 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -35,6 +35,13 @@ uv pip install 'snakemake>=9,<10' 'snakemake-executor-plugin-slurm>=2.7,<3' # (DR6), or the join fails. `sextractor` (the default for image sims) keeps # SExtractor's own NUMBER. Either way make_cat writes # TILE_UNIQUE_ID = tile_id * 10**6 + NUMBER. +# `blend_handling` is ngmix's neighbour treatment: `uberseg` (the default), +# under which tile_detect adds the coadd segmentation stamps to the sexcat as +# SEG_VIGNET and ngmix masks neighbours from them, or `noisefill`. It is chosen +# per campaign: flipping it on an existing run dir reruns tile_detect and the +# shape chain for every finished tile, and reclaimed exposures make that +# destructive, so a comparison needs its own `run:`; a campaign that ran as +# noisefill resumes unchanged with `blend_handling: noisefill`. # The committed launcher loads apptainer/1.4.5 + the /project venv, so a # fresh shell always has the right state. diff --git a/workflow/Snakefile b/workflow/Snakefile index ac9d0760c..f488ed5d5 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -188,6 +188,25 @@ try: except ValueError as err: raise WorkflowError(str(err)) from None +# How ngmix treats a stamp's neighbours; tile.smk's blend handling section +# carries it to tile_detect and tile_ngmix. +BLEND_HANDLING = config.get("blend_handling", "uberseg") +if BLEND_HANDLING not in run_config.BLEND_HANDLINGS: + raise WorkflowError( + f"Invalid blend_handling={BLEND_HANDLING!r}; expected one of " + f"{sorted(run_config.BLEND_HANDLINGS)}.") + +# The committed inis read SP_SEG_VIGNET and SP_BLEND_HANDLING with the modules' +# own noise-fill defaults, and only tile.smk's blend_env exports them (under +# uberseg, the workflow's default). Drop any copy the +# launching shell carries, bare or as the APPTAINERENV_/SINGULARITYENV_ forms +# that survive --cleanenv, so a job inherits none: the slurm executor submits +# with --export=ALL from this process's environment, and each job re-parses +# this file before its shell runs. +for _prefix in ("", "APPTAINERENV_", "SINGULARITYENV_"): + for _name in ("SP_SEG_VIGNET", "SP_BLEND_HANDLING"): + os.environ.pop(_prefix + _name, None) + # SP_PHASE is set by bin/sp and NOWHERE else: `prepare`/`compute` on the two # invocations of `sp run`, `passthrough` on direct commands (`sp --unlock`, `sp # --dag`, `sp exp_psf ...`). It gates the two parse-time side effects — the index diff --git a/workflow/config.yaml b/workflow/config.yaml index f144a8a0e..402b011bd 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -20,6 +20,16 @@ input_type: data # campaign's merged catalogues. Required, and set by the run config. # run: smk-g6 +# How ngmix treats a stamp's neighbours (its BLEND_HANDLING): `uberseg` zeroes +# their weight from the coadd segmentation map, whose stamps tile_detect writes +# into the sexcat as SEG_VIGNET; `noisefill` is ngmix's own default, without +# the segmentation stamps. Chosen per campaign: flipping it on an existing run +# dir reruns tile_detect and everything after it for every finished tile, and +# reclaimed exposures make that destructive -- a comparison needs its own +# `run:`. A campaign that ran as noisefill resumes unchanged with +# `blend_handling: noisefill` in its run config. +blend_handling: uberseg + # Default settings per input_type, can be overridden by a user-defined run config file. # - tile_detection: both run SExtractor on the tile. `unions_catalogue` (data) # also fetches the UNIONS per-tile catalogue and joins the detections to it, diff --git a/workflow/config/cfis/config_tile_Ng_template.ini b/workflow/config/cfis/config_tile_Ng_template.ini index 5d4c510e5..3931e48d7 100644 --- a/workflow/config/cfis/config_tile_Ng_template.ini +++ b/workflow/config/cfis/config_tile_Ng_template.ini @@ -122,3 +122,9 @@ MAG_ZP = 30.0 # via ShapePipe's getexpanded (env-expanded, not plain getint). ID_OBJ_MIN = $NGMIX_ROW_MIN ID_OBJ_MAX = $NGMIX_ROW_MAX + +# Neighbour treatment: noisefill, or uberseg, which reads the tile catalogue's +# SEG_VIGNET column. The workflow sets SP_BLEND_HANDLING to uberseg under its +# default blend_handling: uberseg, and leaves it unset under noisefill +# (workflow/config.yaml). +BLEND_HANDLING = ${SP_BLEND_HANDLING:-noisefill} diff --git a/workflow/config/cfis/config_tile_Sx.ini b/workflow/config/cfis/config_tile_Sx.ini index 43a59c493..6ff1d048d 100644 --- a/workflow/config/cfis/config_tile_Sx.ini +++ b/workflow/config/cfis/config_tile_Sx.ini @@ -112,6 +112,14 @@ BKG_FROM_HEADER = False # OBJECTS, -OBJECTS, SEGMENTATION, APERTURES CHECKIMAGE = BACKGROUND, SEGMENTATION +# Write the SEGMENTATION check image's stamps, cut on each VIGNET's grid, as +# the int32 SEG_VIGNET column ngmix's UberSeg blend handling reads (needs the +# SEGMENTATION check image; SExtractor also computes X/Y_IMAGE_DBL, which +# centre the stamps and are left out of the sexcat). The workflow sets +# SP_SEG_VIGNET under its default blend_handling: uberseg, and leaves +# it unset under noisefill (workflow/rules/tile.smk). +SEG_VIGNET = ${SP_SEG_VIGNET:-False} + # File name suffix for the output sextractor files (optional) SUFFIX = sexcat diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index 107bce4f9..8dc767443 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -54,6 +54,16 @@ catalogues; unpaired rows leave the catalogue. The catalogue is SExtractor run by MegaPipe on the same DR6 image with the same configuration, so the join is exact; on other pixels (a DR5 image) it fails loudly. +``blend_handling: uberseg``, the default, gives ngmix UberSeg's hard mask +instead of its own noise-fill default; the mask reads one coadd segmentation +stamp per object from the sexcat's SEG_VIGNET column. tile_detect cuts it from +SExtractor's SEGMENTATION check image on the grid of each object's VIGNET, +before the join, which relabels every stamp to the catalogue's NUMBER (-1 on +the footprints of rows that leave, 0 on sky). The switch is two prologue +variables and nothing else (blend_env, below), so ``blend_handling: +noisefill`` runs exactly the chain the committed inis run with neither +variable set. + There is no `tile_mask` rule, and there will not be one (PR #847). ShapePipe generates no masks: tiles have no instrument flag image of their own, so tile_detect runs SExtractor with FLAG_IMAGE = False against @@ -472,6 +482,26 @@ rule tile_merge_headers: shell: sp_shell("tile_merge_headers", "config_tile_Mh_exp.ini") +# --- blend handling ------------------------------------------------------- + +# What blend_handling: uberseg (the default) exports into a stage's prologue; +# nothing under noisefill, whose prologues (rerun triggers) are those the +# committed inis need no variable for. The inis read each variable with the +# module's own noise-fill default: +# SEG_VIGNET = ${SP_SEG_VIGNET:-False} in config_tile_Sx.ini, +# BLEND_HANDLING = ${SP_BLEND_HANDLING:-noisefill} in +# config_tile_Ng_template.ini. Being in the prologue, a flip reruns +# tile_detect, and with it the tile's whole shape chain. +BLEND_ENV = { + "tile_detect": {"SP_SEG_VIGNET": "True"}, + "tile_ngmix": {"SP_BLEND_HANDLING": "uberseg"}, +} if BLEND_HANDLING == "uberseg" else {} + + +def blend_env(stage): + return BLEND_ENV.get(stage, {}) + + # The UNIONS per-tile catalogue (CFIS..r.cat), from a local mirror or # from vos, the way tile_get_images fetches the image; tile_detect joins its # SExtractor detections to it. Reads only tile_numbers.txt, which unit_pre @@ -509,15 +539,28 @@ if TILE_DETECTION == "unions_catalogue": def detect_env(tile): - """tile_detect's prologue exports: the catalogue to join, if any. + """tile_detect's prologue exports: the catalogue to join, if any, and + SP_SEG_VIGNET under blend_handling: uberseg. - Empty under tile_detection: sextractor (image simulations), so that no - value left exported in the submitting shell reaches the job. + SP_MATCH_CATALOGUE is empty under tile_detection: sextractor (image + simulations), so that no value left exported in the submitting shell + reaches the job. """ - if TILE_DETECTION != "unions_catalogue": - return {"SP_MATCH_CATALOGUE": ""} - gic = f"{tile_dir(tile)}/output/run_sp_tile_Gic/get_images_runner/output" - return {"SP_MATCH_CATALOGUE": f"{gic}/CFIS_cat{unit_num(tile)}.cat"} + match = "" + if TILE_DETECTION == "unions_catalogue": + gic = f"{tile_dir(tile)}/output/run_sp_tile_Gic/get_images_runner/output" + match = f"{gic}/CFIS_cat{unit_num(tile)}.cat" + return {"SP_MATCH_CATALOGUE": match, **blend_env("tile_detect")} + + +# tile_detect's extra memory under uberseg: add_seg_vignet reads the 400 MB +# SEGMENTATION check image and adds an int32 SEG_VIGNET column the size of +# VIGNET, which the join then holds twice. MEASURED on DR6 tile 202.301 +# (36,022 detections; sextractor_runner as tile_detect runs it, without the +# post-processing): peak RSS 1.53 GiB under noisefill, 2.93 GiB under uberseg. +# Scaled to the 1.84 GiB worst tile above, uberseg needs ~3.5 GiB; 7000 MB is +# 2x that. +DETECT_SEG_MEM_MB = 3000 if BLEND_HANDLING == "uberseg" else 0 # SExtractor object detection on the tile; under unions_catalogue joined to @@ -549,11 +592,11 @@ rule tile_detect: # 4000 MB is 2.1x the worst tile, and an OOM retries at 8000. The runtime # is ~7x the slowest tile, for /scratch I/O on nibi (a 400 MB image and # weight in, a ~430 MB catalogue out). nibi bills max(cores, mem_GB/4), so - # the job bills 1 core-equivalent. + # the job bills 1 core-equivalent; under uberseg, DETECT_SEG_MEM_MB below. threads: 1 retries: 1 resources: - mem_mb = lambda wc, attempt: 4000 * attempt, + mem_mb = lambda wc, attempt: (4000 + DETECT_SEG_MEM_MB) * attempt, runtime = lambda wc, attempt: 20 * attempt shell: sp_shell("tile_detect", "config_tile_Sx.ini") @@ -614,6 +657,13 @@ rule tile_vignets: sp_shell("tile_vignets", f"config_tile_PiViVi_{PSF_MODEL}.ini", check_args=' --run-dir "$SP_LOCAL" --unit {wildcards.tile}') +# Each ngmix chunk's extra memory under uberseg: Tile_cat holding the tile's +# SEG_VIGNET column, a view into the loaded table like VIGNET. Measured on DR6 +# 202.301 (36065 rows, 51x51): Tile_cat peak and resident RSS 955 MiB with the +# column, 595 without (+360 MiB, the column's size). 500 covers 186.307's +# ~37.9k rows (~375 MiB) with a third to spare. +NGMIX_SEG_MEM_MB = 500 if BLEND_HANDLING == "uberseg" else 0 + # ngmix shape measurement — N chunks per tile (D4). Each chunk LOOKS UP its own # CLOSED catalogue-row range in the file tile_vignets materialised at the top of this # group job (TILE_NGMIX_RANGES); the ranges are knowable only at EXECUTION time, @@ -644,7 +694,8 @@ rule tile_ngmix: f"{TILE_DIR}/logs/tile_ngmix_{{chunk}}.json" params: pre = lambda wc: unit_pre("tile_ngmix", wc.tile, - env={"SP_NGMIX_CHUNK": wc.chunk, "NGMIX_N_CHUNKS": NGMIX_CHUNKS}, + env={"SP_NGMIX_CHUNK": wc.chunk, "NGMIX_N_CHUNKS": NGMIX_CHUNKS, + **blend_env("tile_ngmix")}, # Two steps, not `eval "$(...)"`: a command substitution inside eval # discards the script's exit status, so a missing sexcat would fall # through to shapepipe_run with an unset range and fail as something @@ -803,7 +854,10 @@ rule tile_ngmix: # margin over the worst tile measured, and the first thing to give would # be the page cache holding the store. Take it only with a measurement # of cache behaviour under pressure, not on the arithmetic alone. - mem_mb = lambda wc, attempt: 5000 * attempt, + # + # Under uberseg each chunk also holds the tile's SEG_VIGNET column, + # int32 and the size of VIGNET: NGMIX_SEG_MEM_MB below. + mem_mb = lambda wc, attempt: (5000 + NGMIX_SEG_MEM_MB) * attempt, # 120 on the FIRST attempt, and ATTEMPT-SCALED after it. MEASURED # two ways, and the second is why the margin is thinner than it looks: # * alone (job 20795277, tile 198.305 chunk 1): ~76 min for 3547 diff --git a/workflow/scripts/run_config.py b/workflow/scripts/run_config.py index 669a3ec60..7cf9c18c3 100644 --- a/workflow/scripts/run_config.py +++ b/workflow/scripts/run_config.py @@ -28,6 +28,10 @@ # never mention `$run` must still set it. REQUIRED = ("run", "tile_list", "inputs.tiles", "inputs.exposures", "outputs.run_dir", "outputs.index_db") +# ngmix's BLEND_HANDLING values, which `blend_handling` takes. A copy, because +# the Snakefile validates it outside the container; tests/unit holds it equal +# to ngmix_package.ngmix.BLEND_HANDLINGS. +BLEND_HANDLINGS = ("noisefill", "uberseg") def merge(base, over):