diff --git a/astra.yaml b/astra.yaml index 24dd350d3..bf019aa51 100644 --- a/astra.yaml +++ b/astra.yaml @@ -637,9 +637,10 @@ analyses: catalogue_neighbour_marking: label: Neighbour pixels in the catalogue path's VIGNET rationale: >- - ngmix masks a neighbour only where the tile VIGNET is -1e30 (flag - 2**10: zero weight, noise-filled under noisefill), and SExtractor - writes -1e30 on neighbours' footprints and off the image. The + ngmix's noisefill masks a neighbour only where the tile VIGNET is + -1e30 (zero weight, noise-filled; the off-tile rows and columns are + defects instead), and SExtractor writes -1e30 on neighbours' + footprints and off the image. The unions_catalogue converter reproduces that from the catalogue's own r-band segmentation map (CFIS..r.seg.fits.fz, on the tile's pixel grid), so both tile_detection options mask the same way. @@ -1376,22 +1377,26 @@ 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 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, or by + the unions_catalogue converter from the catalogue's segmentation + map; catalogue_neighbour_marking). 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 @@ -1413,7 +1418,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 @@ -1431,11 +1440,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 @@ -1445,10 +1454,18 @@ 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 + 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. + noisefill (the default; the committed config sets no key) gives + weight 0 to the pixels marked -1e30 in the tile VIGNET on other + detections' footprints (about 93% of neighbour-footprint pixels on + a SExtractor-mode sims tile; every one in catalogue mode, from the + segmentation map), excluding the off-tile + rows and columns, which are defects, and replaces them with noise; + unmarked + neighbour pixels keep their weight and light. uberseg ignores the + markers and 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: esheldon/meds get_uberseg returns a weight map (a @@ -1469,15 +1486,16 @@ analyses: 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 + default: noisefill 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: @@ -1508,11 +1526,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, @@ -1523,16 +1543,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 @@ -1542,14 +1564,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 @@ -1565,12 +1589,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/ngmix_package/__init__.py b/src/shapepipe/modules/ngmix_package/__init__.py index e30d36bda..38f22b09a 100644 --- a/src/shapepipe/modules/ngmix_package/__init__.py +++ b/src/shapepipe/modules/ngmix_package/__init__.py @@ -36,8 +36,8 @@ PIXEL_SCALE : float, optional Pixel scale in arcsec. Optional override; when omitted (or non-positive) it is read from the image WCS so it cannot drift from the pixels. Only - sets the centroid-prior width and noise window -- the fit Jacobian is - built per object from the full WCS. + sets the centroid-prior width -- the fit Jacobian is built per object + from the full WCS. SAVE_BATCH : int, optional Save the output catalogue in batches of this size; default is ``-1`` (no batch saving) 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 6c05ef62e..471e62b2c 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 @@ -22,11 +23,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') @@ -515,7 +545,7 @@ class Tile_cat(): 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). + leaves ``self.seg`` unset, which only ``"uberseg"`` needs. """ def __init__( @@ -550,7 +580,7 @@ def get_data(self, cat_path): # 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. + # uberseg unavailable; ``"noisefill"`` does not read it. self.seg = None if self.seg_cat_path: seg_cat = file_io.FITSCatalogue( @@ -603,9 +633,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. @@ -622,6 +658,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 @@ -677,11 +716,11 @@ def pixel_scale_from_wcs(f_wcs_file): """Representative pixel scale (arcsec) from the tile's image WCS. The ngmix fit builds each object's Jacobian from the full per-epoch WCS, - so this scalar only sets the centroid-prior width and the noise-window - scale (see :func:`get_prior`, :func:`get_noise`). A single value read from - the astrometry is therefore sufficient -- and, unlike a hard-coded config - constant, it cannot silently drift from the pixels it describes (this - mirrors SExtractor's ``PIXEL_SCALE 0`` convention). + so this scalar only sets the centroid-prior width (see :func:`get_prior`). + A single value read from the astrometry is therefore sufficient -- and, + unlike a hard-coded config constant, it cannot silently drift from the + pixels it describes (this mirrors SExtractor's ``PIXEL_SCALE 0`` + convention). Parameters ---------- @@ -751,15 +790,34 @@ 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``. + 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, + zeroes the weight of neighbour-side pixels from the coadd + segmentation map and requires ``seg_cat_path``. Defect pixels are + filled under both (see :func:`prepare_ngmix_weights`). seg_cat_path : str, optional Path to the coadd-frame segmentation VIGNET catalogue (see :class:`Tile_cat`). Required when ``blend_handling="uberseg"``. 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 ----- @@ -771,8 +829,8 @@ 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, or ``"uberseg"`` + is selected without ``seg_cat_path``. """ @@ -794,6 +852,10 @@ def __init__( 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 @@ -812,6 +874,11 @@ def __init__( f"Unknown BLEND_HANDLING '{blend_handling}'; expected one of" + f" {BLEND_HANDLINGS}" ) + if defect_fill not in DEFECT_FILLS: + raise ValueError( + f"Unknown DEFECT_FILL '{defect_fill}'; expected one of" + + f" {DEFECT_FILLS}" + ) # Fail fast at construction (not deep in the per-epoch loop) when # uberseg is requested without its segmentation input (shapepipe#776). @@ -863,6 +930,12 @@ def __init__( 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 @@ -874,8 +947,8 @@ def __init__( # Pixel scale: an explicit positive PIXEL_SCALE overrides; otherwise # derive it from the image WCS so it can never drift from the pixels # (mirrors SExtractor's ``PIXEL_SCALE 0`` convention). Only the - # centroid-prior width and noise window use it -- the fit Jacobian is - # built per object from the full WCS. + # centroid-prior width uses it -- the fit Jacobian is built per + # object from the full WCS. if pixel_scale is None or pixel_scale <= 0: self._pixel_scale = pixel_scale_from_wcs( self._vignet_cat.f_wcs_file @@ -1281,6 +1354,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( @@ -1314,10 +1389,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 @@ -1358,6 +1441,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( @@ -1426,6 +1510,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() @@ -1447,11 +1538,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. - @sc [decision:shape_measurement.central_defect_veto,decision:shape_measurement.epoch_masked_fraction_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``. + + 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 @@ -1518,19 +1673,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 @@ -1564,6 +1730,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) @@ -1573,9 +1740,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 @@ -1586,6 +1751,52 @@ 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; SExtractor writes it, and in + catalogue mode the converter paints it from the catalogue's segmentation + map (``read_ext_sexcat._extract_vignets``). 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. @@ -1805,9 +2016,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 ---------- @@ -1869,33 +2080,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. @@ -1905,49 +2343,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." + + " Set SEG_VIGNET_PATH on the ngmix run (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) ) @@ -1960,34 +2439,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. @@ -2031,7 +2510,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). @@ -2041,6 +2521,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 ------- @@ -2069,16 +2554,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) @@ -2299,7 +2786,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. @@ -2324,9 +2811,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 @@ -2346,6 +2833,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 ------- @@ -2379,6 +2869,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 75d26f60a..08113aa84 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( @@ -63,7 +69,7 @@ def ngmix_runner( # Pixel scale -- optional override. When absent (or non-positive) it is # derived from the image WCS inside Ngmix, so it cannot drift from the - # pixels. Only the centroid-prior width and noise window use it. + # pixels. Only the centroid-prior width uses it. if config.has_option(module_config_sec, "PIXEL_SCALE"): pixel_scale = config.getfloat(module_config_sec, "PIXEL_SCALE") else: @@ -138,10 +144,13 @@ 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") else: @@ -155,6 +164,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 @@ -216,6 +266,10 @@ def ngmix_runner( 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/read_ext_sexcat_package/read_ext_sexcat.py b/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py index ff4f58cc8..4b1d66012 100644 --- a/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py +++ b/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py @@ -48,8 +48,9 @@ def _build_ldac_imhead(img_header): # SExtractor's VIGNET value for pixels that are not the object's: off the -# image and on a neighbour's segmentation footprint. ngmix flags these pixels -# (tile VIGNET == -1e30) and noise-fills them at zero weight. +# image and on a neighbour's segmentation footprint. ngmix splits these pixels +# (tile VIGNET == -1e30, ``ngmix.split_tile_markers``): the off-image rows and +# columns are defects, the footprint pixels its neighbour mask. BIG = -1e30 # The label of a footprint no catalogue object claims in the relabelled 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 ed87b4be2..a69aa0e69 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 hypothesis import given from hypothesis import strategies as st @@ -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), ) psf = dict( n_epoch=1, g_psf=[.01, -.01], g_psf_err=[.001, .001], @@ -641,7 +644,9 @@ 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), ) @@ -658,6 +663,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) diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py new file mode 100644 index 000000000..f068c1290 --- /dev/null +++ b/tests/module/test_ngmix_defect_fill.py @@ -0,0 +1,1268 @@ +"""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.read_ext_sexcat_package import read_ext_sexcat + + +# --- 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] + + +# --- 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, 0.186, 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, 0.186, 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, 0.186, 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) + + +# --- Catalogue mode: converter-built tile VIGNETs --------------------------- +# +# In catalogue mode the tile VIGNET comes from read_ext_sexcat, which paints +# -1e30 off the image and on other footprints of the DR6 segmentation map. +# 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 _converter_stamps(seg, number, x, y): + """Converter tile VIGNETs on a unit image, and each stamp's off-image + mask.""" + relabelled, _ = read_ext_sexcat.relabel_seg(seg, number, x, y) + vignets = read_ext_sexcat._extract_vignets( + np.ones(seg.shape, np.float32), x, y, N_STAMP, seg=relabelled, + number=number, + ) + half = N_STAMP // 2 + ny, nx = seg.shape + off_image = [] + for xi, yi in zip(x, y): + rows = int(np.rint(yi)) - 1 - half + np.arange(N_STAMP) + cols = int(np.rint(xi)) - 1 - half + np.arange(N_STAMP) + off_image.append( + ((rows < 0) | (rows >= ny))[:, None] + | ((cols < 0) | (cols >= nx))[None, :] + ) + return vignets, off_image + + +def test_dr6_converter_stamps_split_into_off_image_and_neighbours(): + """On the real 202.301 patch, every converter stamp splits into its + off-image pixels (off-tile) and its other -1e30 pixels (neighbours). + + Failure modes: the converter's 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 + number = np.array(objects["NUMBER"]) + x, y = np.array(objects["X_IMAGE"]), np.array(objects["Y_IMAGE"]) + vignets, off_image = _converter_stamps(seg, number, 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 + number, x, y = np.array([1, 5]), np.array([15.0, 30.0]), np.array([30.0, 22.0]) + vignets, off_image = _converter_stamps(seg, number, 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 0c80486b2..5c20789f6 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). """ @@ -226,7 +226,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 +239,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 +256,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 +275,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(): @@ -422,3 +406,4 @@ def test_runner_missing_seg_vignet_file_raises(tmp_path): "NGMIX_RUNNER", _RecordingLogger(), ) + 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/universes/committed.yaml b/universes/committed.yaml index bc34d7601..f90c5874b 100644 --- a/universes/committed.yaml +++ b/universes/committed.yaml @@ -54,9 +54,9 @@ analyses: psf_likelihood_noise: psf_noise_1em5 megacam_ccd_flip: megapipe_flip defect_fill: noise - blend_handling: none + blend_handling: noisefill epoch_masked_fraction_cut: one_third - central_defect_veto: disabled + central_defect_veto: fixed_radii catalogue_assembly: decisions: star_galaxy_classification: deferred_downstream