diff --git a/astra.yaml b/astra.yaml index 87c844909..4237b1140 100644 --- a/astra.yaml +++ b/astra.yaml @@ -1329,15 +1329,11 @@ 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". 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, 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 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 @@ -1377,11 +1373,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 @@ -1391,9 +1387,10 @@ 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. + Covers only pixels shared with a neighbour; defect_fill does not + depend on it: defects are noise-filled under every BLEND_HANDLING. + none (the default; the committed config + sets no key) leaves neighbours fully weighted and untouched. uberseg zeroes the weight of pixels nearer a neighbour's coadd segmentation footprint than the target's (DILATE_NEIGHBOUR, default 1) and leaves the image untouched, as official uberseg does: @@ -1418,7 +1415,7 @@ analyses: default: none options: none: - label: No neighbour treatment (BLEND_HANDLING = noisefill) + label: No neighbour treatment (BLEND_HANDLING = none) description: >- Neighbour pixels keep their full weight and image values, so the fit sees all neighbour light, which biases shapes toward @@ -1454,10 +1451,10 @@ 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 + 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 + 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 @@ -1469,16 +1466,17 @@ 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 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. + Values: + EPOCH_CENTRAL_DEFECT_RADIUS = 10. + 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 @@ -1488,14 +1486,13 @@ 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 flag + bit, including the tile-coverage bit 2**10 set where the tile vignet + is off-image), zero-weight pixels and invalid-RMS pixels, the set + defect_fill noise-fills. 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 @@ -1511,12 +1508,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 66323822c..012b3264b 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/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 62502c4cd..716afb9a0 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 @@ -25,7 +26,18 @@ from shapepipe.pipeline import file_io # Neighbour treatments selectable with the BLEND_HANDLING option. -BLEND_HANDLINGS = ("noisefill", "uberseg") +BLEND_HANDLINGS = ("none", "uberseg") + +# 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 # @sc [decision:shape_measurement.metacal_scheme] METACAL_TYPES = ('noshear', '1p', '1m', '2p', '2m') @@ -197,7 +209,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__( @@ -232,7 +244,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; ``"none"`` does not read it. self.seg = None if self.seg_cat_path: seg_cat = file_io.FITSCatalogue( @@ -287,7 +299,7 @@ def __init__( self.flags = [] 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 "none". 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. @@ -304,6 +316,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 @@ -359,11 +374,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 ---------- @@ -432,16 +447,26 @@ class Ngmix(object): propagated on the vignette. ``"hsm"`` re-centers on the 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``. + blend_handling : {"none", "uberseg"}, optional + 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`). 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 masked 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 masked + (see :func:`prepare_postage_stamps`); the default is + ``EPOCH_MASKED_FRACTION_CUT``. Notes ----- @@ -472,10 +497,12 @@ def __init__( id_obj_max=-1, bkg_sub=True, centroid_source="wcs", - blend_handling="noisefill", + blend_handling="none", 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, ): # Base count = catalogue + vignets, excluding the f_wcs headers (passed @@ -545,6 +572,8 @@ 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._w_log = w_log @@ -556,8 +585,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 @@ -993,6 +1022,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 for i_tile, obj_id in enumerate(tile_cat.obj_id): @@ -1022,10 +1053,14 @@ 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, ) + 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 @@ -1139,6 +1174,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}" + ) vignet_cat.close() # Put all results together @@ -1158,11 +1200,47 @@ 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, ): - """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 + 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 + the alternative to test. - @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``. + + 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 @@ -1231,17 +1309,23 @@ def prepare_postage_stamps( 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 - 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 + # noise-filled (epoch-cut-on-defect-mask), or when a filled pixel + # would sit near the object (epoch-central-defect-veto). + 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): + 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 @@ -1284,9 +1368,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 @@ -1516,9 +1598,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 ---------- @@ -1580,33 +1662,163 @@ 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 flag, or (when a + background RMS map is given) a non-finite or non-positive RMS. This one + set is zero-weighted and noise-filled by :func:`prepare_ngmix_weights` + and counted by the epoch cuts in :func:`prepare_postage_stamps`. 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 + 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 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 (:func:`defect_mask`). + noise : numpy.ndarray + Noise realisation on the stamp grid. + + Returns + ------- + numpy.ndarray + ``image`` with ``defect`` pixels taken from ``noise``. + """ + return np.where(defect, noise, image) + + def prepare_ngmix_weights( gal, weight, flag, rng, bkg_rms=None, - blend_handling="noisefill", seg=None, object_number=None, + blend_handling="none", seg=None, object_number=None, dilate_neighbour=0, ): - """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 replaced by an independent noise + realisation at their background RMS (:func:`fill_defects`), under either + ``blend_handling``. ``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 + 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. Parameters ---------- gal : numpy.ndarray + Background-subtracted galaxy stamp. weight : numpy.ndarray + Exposure weight stamp; zero marks a defect. flag : numpy.ndarray + 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. - 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`). + 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 : {"none", "uberseg"}, optional + Neighbour treatment. ``"none"`` (default) leaves neighbour + pixels weighted and untouched. ``"uberseg"`` 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. @@ -1620,12 +1832,18 @@ def prepare_ngmix_weights( 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 noise-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``. + + Raises + ------ + ValueError + If ``blend_handling`` 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 blend_handling not in BLEND_HANDLINGS: @@ -1633,32 +1851,39 @@ def prepare_ngmix_weights( 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) + clean = ~defect 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 + clean.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[clean] = 1.0 / bkg_rms[clean] ** 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) ) @@ -1671,33 +1896,20 @@ 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) - 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, + blend_handling="none", seg=None, object_number=None, dilate_neighbour=0, ): """Build an ngmix Observation for a single galaxy epoch. @@ -1740,9 +1952,9 @@ def make_ngmix_observation( Sub-pixel ``[row, col]`` coadd-centroid offset propagated from the stamp extractor. Required for ``centroid_source="wcs"`` (ignored for ``"hsm"``). - blend_handling : {"noisefill", "uberseg"}, optional + blend_handling : {"none", "uberseg"}, optional Neighbour treatment passed through to :func:`prepare_ngmix_weights`; - the default ``"noisefill"`` is the historical behaviour. + the default ``"none"`` leaves neighbour pixels untouched. seg : numpy.ndarray, optional Segmentation map on the stamp grid. Required for ``blend_handling="uberseg"`` (ignored otherwise). @@ -1785,11 +1997,12 @@ def make_ngmix_observation( 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) @@ -2009,7 +2222,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, + blend_handling="none", object_number=None, dilate_neighbour=0, metacal_psf="fitgauss", ): """Do Ngmix Metacal. @@ -2033,10 +2246,10 @@ def do_ngmix_metacal( coadd-centroid offset in ``stamp.offsets``, propagated from the stamp extractor); ``"hsm"`` uses the adaptive-moment centroid from the stamp pixels — see that function. - blend_handling : {"noisefill", "uberseg"}, optional + blend_handling : {"none", "uberseg"}, optional Neighbour treatment passed through to - :func:`make_ngmix_observation`; the default ``"noisefill"`` is the - historical behaviour. ``"uberseg"`` consumes ``stamp.segs`` and + :func:`make_ngmix_observation`; the default ``"none"`` leaves + neighbour pixels untouched. ``"uberseg"`` consumes ``stamp.segs`` and ``object_number``. object_number : int, optional Central object's segmentation label — its SExtractor ``NUMBER`` diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 63f03872b..2c73cfc6e 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -11,7 +11,11 @@ from sqlitedict import SqliteDict from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.ngmix_package.ngmix import Ngmix +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_CENTRAL_DEFECT_RADIUS, + EPOCH_MASKED_FRACTION_CUT, + Ngmix, +) @module_runner( @@ -52,7 +56,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: @@ -127,14 +131,16 @@ 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: "none" (default) leaves neighbour pixels + # 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. if config.has_option(module_config_sec, "BLEND_HANDLING"): blend_handling = config.get(module_config_sec, "BLEND_HANDLING") else: - blend_handling = "noisefill" + blend_handling = "none" # DILATE_NEIGHBOUR (optional): binary-dilation iterations enlarging the # uberseg neighbour mask, to absorb the few-pixel coadd-vs-epoch seg-overlay @@ -144,6 +150,24 @@ def ngmix_runner( else: dilate_neighbour = 1 + # EPOCH_CENTRAL_DEFECT_RADIUS (optional, pixels): drop an epoch when a + # masked pixel 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 masked (flagged, zero-weight or invalid-RMS). + 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 + # 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 @@ -199,6 +223,8 @@ 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, ) # Process ngmix shape measurement and metacalibration 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_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py new file mode 100644 index 000000000..0480c7ef7 --- /dev/null +++ b/tests/module/test_ngmix_defect_fill.py @@ -0,0 +1,532 @@ +"""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 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. +""" + +import re +from types import SimpleNamespace + +import numpy as np +import numpy.testing as npt +from astropy.io import fits +from astropy.wcs import WCS +from hypothesis import given +from hypothesis import strategies as st +from sqlitedict import SqliteDict +from shapepipe.modules.ngmix_package import ngmix as ngmix_module + +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_CENTRAL_DEFECT_RADIUS, + Ngmix, + make_ngmix_observation, + prepare_ngmix_weights, + prepare_postage_stamps, + uberseg_weight, +) + + +# --- 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(["none", "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))) + 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) + + +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 and EPOCH_MASKED_FRACTION_CUT reach Ngmix + as configured, and default to the module constants 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 in ( + ({"EPOCH_CENTRAL_DEFECT_RADIUS": "7.5", + "EPOCH_MASKED_FRACTION_CUT": "0.1"}, 7.5, 0.1), + ({}, EPOCH_CENTRAL_DEFECT_RADIUS, + ngmix_module.EPOCH_MASKED_FRACTION_CUT), + ): + 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 + + +# --- 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 + ) diff --git a/tests/module/test_ngmix_uberseg.py b/tests/module/test_ngmix_uberseg.py index 0c80486b2..4d6ebb32a 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.""" +def test_none_ignores_seg_and_dilate_kwargs(): + """Under BLEND_HANDLING = none, 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,15 @@ 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-neighbours-raw). + """ npix = 41 gal, flag, weight = _gal_flag_weight(npix=npix) seg, centre, neigh = two_object_seg(npix=npix, sep=12) @@ -294,13 +274,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 +405,11 @@ def test_runner_missing_seg_vignet_file_raises(tmp_path): "NGMIX_RUNNER", _RecordingLogger(), ) + + +def test_retired_noisefill_name_is_rejected(): + """Do not silently interpret the retired neighbour option.""" + with pytest.raises(ValueError, match="noisefill"): + prepare_ngmix_weights(np.ones((5, 5)), np.ones((5, 5)), + np.zeros((5, 5)), np.random.RandomState(3), + blend_handling="noisefill") 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 63b6e795b..4e3e02bec 100644 --- a/universes/committed.yaml +++ b/universes/committed.yaml @@ -54,7 +54,7 @@ analyses: defect_fill: noise blend_handling: none epoch_masked_fraction_cut: one_third - central_defect_veto: disabled + central_defect_veto: fixed_radii catalogue_assembly: decisions: star_galaxy_classification: deferred_downstream