diff --git a/astra.yaml b/astra.yaml index 4237b1140..385d2151c 100644 --- a/astra.yaml +++ b/astra.yaml @@ -1329,8 +1329,8 @@ analyses: 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". Under every BLEND_HANDLING, defect pixels are replaced with - independent noise at the per-pixel background RMS when supplied, or + 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; `BLEND_HANDLING = uberseg` additionally zeroes the weight of neighbour-side pixels and leaves their image values untouched. The committed fill @@ -1355,7 +1355,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 @@ -1466,11 +1470,12 @@ 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. The 7 px interpolated radius arrives with - defect_fill's interpolate option (feat/defect-interpolation); until - then every defect is noise-filled and the 10 px radius applies. + 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_CENTRAL_DEFECT_RADIUS = 10; + EPOCH_INTERPOLATED_DEFECT_RADIUS = 7. default: fixed_radii options: disabled: 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..713e326f7 --- /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, defect, 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, + ) & ~defect + 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, defect, target): + """Replace the ``target`` pixels of every plane by a Clough-Tocher + interpolant of the clean pixels around them. + + @sc [decision:shape_measurement.defect_fill] shared-rotation-averaged-interpolant + The support is the clean pixels within ``SUPPORT_RADIUS`` (4 px) of the + target; no defect pixel enters it, so defect 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)``. + defect : numpy.ndarray of bool + Every defect pixel, shape ``(ny, nx)``; none of them supports the + interpolant. + 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) + defect = np.asarray(defect, 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(defect, 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 716afb9a0..e6f59c5f0 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -23,11 +23,20 @@ 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 # Neighbour treatments selectable with the BLEND_HANDLING option. BLEND_HANDLINGS = ("none", "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`). @@ -39,6 +48,12 @@ # @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') @@ -451,8 +466,7 @@ class Ngmix(object): Neighbour treatment. ``"none"`` (default) leaves neighbour pixels untouched; ``"uberseg"`` zeroes the weight of neighbour-side pixels from the coadd segmentation map and requires ``seg_cat_path``. Defect - pixels are noise-filled under both (see - :func:`prepare_ngmix_weights`). + 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"``. @@ -467,6 +481,14 @@ class Ngmix(object): Drop an epoch when more than this fraction of its stamp is masked (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 ----- @@ -478,8 +500,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``. """ @@ -503,6 +525,8 @@ def __init__( 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 @@ -521,6 +545,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). @@ -574,6 +603,10 @@ def __init__( 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 @@ -1055,6 +1088,10 @@ def process(self): 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) @@ -1101,6 +1138,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( @@ -1202,6 +1240,8 @@ def prepare_postage_stamps( 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, ): """Gather one object's epoch stamps, dropping epochs its defects spoil. @@ -1209,8 +1249,9 @@ def prepare_postage_stamps( 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 - noise-fills. Counting flags alone would keep epochs whose filled area - exceeds the cut. The default is 1/3; 10%, the DES Y3 and Y6 value, is + fills. 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. Parameters @@ -1234,6 +1275,14 @@ def prepare_postage_stamps( 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 ------- @@ -1316,14 +1365,17 @@ def prepare_postage_stamps( else None ) # Drop the epoch when too much of it would be zero-weighted and - # noise-filled (epoch-cut-on-defect-mask), or when a filled pixel - # would sit near the object (epoch-central-defect-veto). + # 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 has_central_defect(masked, epoch_central_defect_radius): + if central_defect_vetoes( + masked, epoch_central_defect_radius, defect_fill, + epoch_interpolated_defect_radius, + ): stamp.epoch_cuts["central_veto"] += 1 continue @@ -1746,6 +1798,59 @@ def has_central_defect(defect, radius): 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. @@ -1778,24 +1883,42 @@ def fill_defects(image, defect, noise): def prepare_ngmix_weights( gal, weight, flag, rng, bkg_rms=None, blend_handling="none", seg=None, object_number=None, - dilate_neighbour=0, + dilate_neighbour=0, defect_fill="noise", ): """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 replaced by an independent noise - realisation at their background RMS (:func:`fill_defects`), under either - ``blend_handling``. ``blend_handling`` decides only the neighbour side. + 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 side. @sc [decision:shape_measurement.defect_fill,decision:shape_measurement.blend_handling] defects-filled-neighbours-raw - Every pixel of ``defect_mask`` is zero-weighted and noise-filled whatever - ``blend_handling`` is, and the filled set equals the zero-weight defect - set. ``blend_handling`` acts on the neighbour side only: uberseg zeroes + Every pixel of ``defect_mask`` is zero-weighted and filled whatever + ``blend_handling`` is, and the filled set equals the defect set. + ``blend_handling`` acts on the neighbour side only: uberseg zeroes the weight of pixels nearer a neighbour's footprint 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. Neighbour pixels are never noise-filled. The noise - image covers the whole stamp and is independent of both. + an unsheared edge. Neighbour pixels are never filled. The noise image + covers the whole stamp and is independent of both. + + @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 ---------- @@ -1828,24 +1951,32 @@ 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. Returns ------- numpy.ndarray - Galaxy image with defect pixels noise-filled. + Galaxy image with defect pixels filled. numpy.ndarray Inverse-variance weight map for ngmix. numpy.ndarray Noise image: an independent realisation over the whole stamp, for - metacal's ``fixnoise``. + metacal's ``fixnoise``, interpolated where the galaxy image is. Raises ------ ValueError - If ``blend_handling`` is unknown, or ``"uberseg"`` lacks ``seg`` or - ``object_number``. + 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" @@ -1861,19 +1992,29 @@ def prepare_ngmix_weights( defect = defect_mask(weight, flag, bkg_rms) clean = ~defect + 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 = ( + ~(defect | 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 = ( - clean.astype(float) / sig_noise ** 2 + weighted.astype(float) / sig_noise ** 2 if sig_noise > 0 else np.zeros_like(gal, dtype=float) ) else: weight_map = np.zeros_like(gal, dtype=float) - weight_map[clean] = 1.0 / bkg_rms[clean] ** 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 @@ -1897,6 +2038,13 @@ 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, defect, noise_img_gal) + if interpolated.any(): + # One operator for the image and the noise image; a pixel whose + # support is degenerate keeps its noise fill. + filled = interpolate_defects([gal, noise_img], defect, interpolated) + done = interpolated & np.all(np.isfinite(filled), axis=0) + gal_filled[done] = filled[0][done] + noise_img = np.where(done, filled[1], noise_img) if blend_handling == "uberseg": weight_map = uberseg_weight( @@ -1910,7 +2058,7 @@ def make_ngmix_observation( gal, weight, flag, psf, wcs, rng, bkg_rms=None, centroid_source="wcs", offset=None, blend_handling="none", seg=None, object_number=None, - dilate_neighbour=0, + dilate_neighbour=0, defect_fill="noise", ): """Build an ngmix Observation for a single galaxy epoch. @@ -1964,6 +2112,9 @@ 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"``. Returns ------- @@ -1992,7 +2143,7 @@ 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, ) if centroid_source == "hsm": @@ -2223,7 +2374,7 @@ def make_runners(prior, flux_guess, rng): def do_ngmix_metacal( stamp, prior, flux_guess, rng, centroid_source="wcs", blend_handling="none", object_number=None, dilate_neighbour=0, - metacal_psf="fitgauss", + metacal_psf="fitgauss", defect_fill="noise", ): """Do Ngmix Metacal. @@ -2270,6 +2421,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 ------- @@ -2303,6 +2457,7 @@ 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, ) gal_obs_list.append(gal_obs) diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 2c73cfc6e..3099cd067 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -13,6 +13,7 @@ from shapepipe.modules.module_decorator import module_runner from shapepipe.modules.ngmix_package.ngmix import ( EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_INTERPOLATED_DEFECT_RADIUS, EPOCH_MASKED_FRACTION_CUT, Ngmix, ) @@ -135,8 +136,8 @@ def ngmix_runner( # weighted and untouched; "uberseg" 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 or - # invalid-RMS) are zero-weighted and noise-filled under both; see - # prepare_ngmix_weights. + # invalid-RMS) 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: @@ -168,6 +169,27 @@ def ngmix_runner( 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 @@ -225,6 +247,8 @@ def ngmix_runner( 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/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_defect_fill.py b/tests/module/test_ngmix_defect_fill.py index 0480c7ef7..cbbbaa3b0 100644 --- a/tests/module/test_ngmix_defect_fill.py +++ b/tests/module/test_ngmix_defect_fill.py @@ -2,12 +2,14 @@ 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 with noise at its background RMS, whatever -``BLEND_HANDLING`` is. 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. +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. """ import re @@ -15,6 +17,7 @@ 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 @@ -22,8 +25,14 @@ 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, @@ -461,8 +470,9 @@ def getboolean(self, _sec, key, fallback=False): def test_runner_threads_the_epoch_cut_options(tmp_path, monkeypatch): - """EPOCH_CENTRAL_DEFECT_RADIUS and EPOCH_MASKED_FRACTION_CUT reach Ngmix - as configured, and default to the module constants when absent. + """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. @@ -483,11 +493,15 @@ def process(self): for path in inputs: SqliteDict(path).close() - for options, radius, fraction in ( + for options, radius, fraction, fill, interpolated_radius in ( ({"EPOCH_CENTRAL_DEFECT_RADIUS": "7.5", - "EPOCH_MASKED_FRACTION_CUT": "0.1"}, 7.5, 0.1), + "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), + ngmix_module.EPOCH_MASKED_FRACTION_CUT, "noise", + EPOCH_INTERPOLATED_DEFECT_RADIUS), ): runner_module.ngmix_runner( inputs, {"output": str(tmp_path)}, "-001-001", @@ -495,6 +509,9 @@ def process(self): ) 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 ------- @@ -530,3 +547,187 @@ def test_hsm_centroid_ignores_raw_defect_values(): 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", ["none", "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", + ) 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