diff --git a/example/cfis/config_tile_Ng_template.ini b/example/cfis/config_tile_Ng_template.ini index 5bfbc2c7..6214765f 100644 --- a/example/cfis/config_tile_Ng_template.ini +++ b/example/cfis/config_tile_Ng_template.ini @@ -83,5 +83,14 @@ PIXEL_SCALE = 0.186 # noisy for stars and flagged as incorrect — see #767). CENTROID_SOURCE = wcs +# SEED_FROM_POSITION (optional): FALSE (default) seeds one RNG per tile. +# TRUE seeds a per-object RNG from the object's sky position instead, so that +# metacal's fixnoise counter-noise and the fit guesses are identical for the +# same object across Pujol image-simulation shear branches and cancel in the +# branch difference — shrinking the m-bias error. This is FOR IMAGE +# SIMULATIONS ONLY (Pujol noise cancellation, ngmix#796); leave it off (or +# omit) for real data, where it has no benefit. +# SEED_FROM_POSITION = FALSE + ID_OBJ_MIN = X ID_OBJ_MAX = X diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 371264e2..12d8f245 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -119,6 +119,64 @@ def get_prior(pixel_scale, rng, T_range=None, F_range=None): ) +def position_seed(ra, dec, ccd): + """Deterministic RNG seed from an object's sky position (ngmix#796). + + For image-simulation m-bias with the Pujol estimator, the same scene is + simulated under different applied shears ("image branches") and the shear + response is read from the branch difference of the SAME objects. Metacal's + ``fixnoise`` adds a counter-noise realisation drawn from an RNG; if that RNG + is seeded per tile, an object gets *different* added noise in different + branches (detection order differs), and the noise fails to cancel in the + difference, inflating ``sigma_m``. Seeding the per-object RNG from sky + position instead makes the same object draw the same added noise (and the + same fit guesses) in every branch, so both cancel and the m-bias error + shrinks. Off in production; a knob for the sim path only. + + Box math (kept exactly as Fabian's issue #796):: + + box_x = floor(ra * 3600 / 3) + (ccd + 1) + box_y = floor(dec * 3600 / 3) + (ccd + 2) + + The 3-arcsec boxes make the seed robust to detection-order changes and to + small centroid jitter within a box: the same physical object lands in the + same box across branches and so gets the same seed. The ``+ccd`` offsets + disambiguate exposure overlap. + + **Deviation from #796.** Fabian combines the boxes as ``box_x + box_y``, + which collides along anti-diagonals — e.g. boxes ``(10, 20)`` and + ``(11, 19)`` both give seed 30, so two distinct objects share a noise + stream. We instead combine them with a Cantor pairing (a bijection on + non-negative integers), after a zig-zag fold that maps signed box indices + (dec can be negative) onto non-negative ones. The result is reduced mod + ``2**32`` to land in the valid ``numpy.random.RandomState`` seed range + ``[0, 2**32)``. The box math is untouched; only the collision-prone sum is + replaced. + + Parameters + ---------- + ra, dec : float + Object sky position in degrees (first-epoch coordinates; all epochs of + one object share them). + ccd : int + CCD number of the first epoch. + + Returns + ------- + int + Seed in ``[0, 2**32)`` for ``numpy.random.RandomState``. + """ + box_x = int(np.floor((ra * 3600) / 3) + (ccd + 1)) + box_y = int(np.floor((dec * 3600) / 3) + (ccd + 2)) + # Zig-zag fold signed -> non-negative (0,-1,1,-2,... -> 0,1,2,3,...) so the + # Cantor pairing, which is a bijection on the non-negative integers, stays + # injective for southern (dec<0) boxes. + zx = 2 * box_x if box_x >= 0 else -2 * box_x - 1 + zy = 2 * box_y if box_y >= 0 else -2 * box_y - 1 + cantor = (zx + zy) * (zx + zy + 1) // 2 + zy + return cantor % (2 ** 32) + + class Tile_cat(): """Tile_cat. @@ -186,6 +244,10 @@ def __init__( self.wcs = [] self.ra = [] self.dec = [] + # CCD number of the first epoch, used only to build the per-object + # position seed (see :func:`position_seed`, ngmix#796). ``None`` when + # position seeding is off, so the field costs nothing on the hot path. + self.ccd = None self.bkg_sub = bkg_sub self.megacam_flip = megacam_flip @@ -273,6 +335,13 @@ class Ngmix(object): (robust for galaxies); ``"wcs"`` uses the catalog sky position projected through the WCS (better for stars, whose HSM moments are noisy). See :func:`make_ngmix_observation`. + seed_from_position : bool, optional + If ``True``, replace the tile-level RNG with a per-object RNG seeded + from the object's sky position (:func:`position_seed`) inside the + object loop, so metacal's ``fixnoise`` counter-noise and the fit + guesses cancel across Pujol image-simulation branches (ngmix#796). The + default ``False`` leaves the production path byte-identical. See + :func:`position_seed` for the physics and the seed construction. Raises ------ @@ -294,6 +363,7 @@ def __init__( id_obj_min=-1, id_obj_max=-1, centroid_source="hsm", + seed_from_position=False, ): if len(input_file_list) not in {6, 7}: @@ -325,6 +395,7 @@ def __init__( self._id_obj_min = id_obj_min self._id_obj_max = id_obj_max self._centroid_source = centroid_source + self._seed_from_position = seed_from_position self._w_log = w_log @@ -332,6 +403,11 @@ def __init__( seed = int(''.join(re.findall(r'\d+', self._file_number_string))) self._rng = np.random.RandomState(seed) self._w_log.info(f'Random generator initialisation seed = {seed}') + if self._seed_from_position: + self._w_log.info( + 'SEED_FROM_POSITION on: per-object RNG seeded from sky position' + ' for Pujol noise cancellation (image sims, ngmix#796)' + ) @classmethod def MegaCamFlip(self, vign, ccd_nb): @@ -684,6 +760,25 @@ def process(self): n_no_epoch += 1 continue + # Position-seeded per-object RNG for Pujol noise cancellation in + # image sims (ngmix#796): the same object gets the same fixnoise + # counter-noise and fit guesses in every shear branch, so both + # cancel in the branch difference. The prior is rebuilt from the + # same per-object RNG because the guesser draws its initial guess + # via prior.sample() (ngmix guessers.py), which consumes the RNG the + # prior was CONSTRUCTED with — so a per-object rng alone would leave + # the guess drawing from the shared tile stream and break + # cancellation. Off in production, where the single tile-level + # self._rng and the tile-level prior carry the whole loop. + if self._seed_from_position: + obj_rng = np.random.RandomState( + position_seed(stamp.ra[0], stamp.dec[0], stamp.ccd) + ) + obj_prior = get_prior(self._pixel_scale, obj_rng) + else: + obj_rng = self._rng + obj_prior = prior + try: flux_guess = ( tile_cat.flux[i_tile] @@ -692,9 +787,9 @@ def process(self): ) res, psf_res, psf_orig_res = do_ngmix_metacal( stamp, - prior, + obj_prior, flux_guess, - self._rng, + obj_rng, centroid_source=self._centroid_source, ) except Exception as ee: @@ -856,6 +951,11 @@ def prepare_postage_stamps(vignet, obj_id, i_tile, tile_cat): stamp.wcs.append(epoch_wcs) 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 + # for the position seed (ngmix#796). All epochs of one object share the + # ra/dec above, so first-epoch CCD pins one deterministic seed stream. + if stamp.ccd is None: + stamp.ccd = int(ccd_n) return stamp diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 31028c7f..99bb9573 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -88,6 +88,17 @@ def ngmix_runner( else: centroid_source = "wcs" + # Seed the per-object RNG from sky position instead of per tile, so + # metacal's fixnoise counter-noise (and the fit guesses) cancel across + # Pujol image-simulation shear branches (ngmix#796). Default False leaves + # the production path byte-identical. + if config.has_option(module_config_sec, "SEED_FROM_POSITION"): + seed_from_position = config.getboolean( + module_config_sec, "SEED_FROM_POSITION" + ) + else: + seed_from_position = False + # Initialise class instance ngmix_inst = Ngmix( input_file_list, @@ -101,6 +112,7 @@ def ngmix_runner( id_obj_min=id_obj_min, id_obj_max=id_obj_max, centroid_source=centroid_source, + seed_from_position=seed_from_position, ) # Process ngmix shape measurement and metacalibration diff --git a/tests/module/test_ngmix.py b/tests/module/test_ngmix.py index 37db4e3f..3fcb7a51 100644 --- a/tests/module/test_ngmix.py +++ b/tests/module/test_ngmix.py @@ -135,6 +135,84 @@ def test_metacal_is_reproducible_with_fixed_seed(): npt.assert_array_equal(_metacal_noshear_g(42), _metacal_noshear_g(42)) +# --- position_seed: the Pujol noise-cancellation seed (ngmix#796) --- +# The seed is a pure function of (ra, dec, ccd) that must be (a) invariant to +# sub-box centroid jitter, (b) distinct for neighbouring boxes, (c) in the +# valid RandomState range, and (d) free of the anti-diagonal collision of +# Fabian's original ``box_x + box_y`` sum. These properties are what make the +# per-object noise/guess draws cancel across image-simulation shear branches. + +# In-survey positions: RA in [0, 360), Dec in [-40, 40], MegaCam CCDs 0..39. +_ra = st.floats(min_value=0.0, max_value=359.999, allow_nan=False) +_dec = st.floats(min_value=-40.0, max_value=40.0, allow_nan=False) +_ccd = st.integers(min_value=0, max_value=39) + + +@given(_ra, _dec, _ccd) +def test_position_seed_in_valid_randomstate_range(ra, dec, ccd): + """Seed always lands in [0, 2**32), the RandomState-acceptable range.""" + from shapepipe.modules.ngmix_package.ngmix import position_seed + + seed = position_seed(ra, dec, ccd) + assert 0 <= seed < 2 ** 32 + + +@given(_ra, _dec, _ccd, st.floats(min_value=-1.4, max_value=1.4)) +def test_position_seed_invariant_to_sub_box_jitter(ra, dec, ccd, jitter_arcsec): + """Centroid jitter < half the 3-arcsec box leaves the seed unchanged. + + An object re-detected with a slightly different centroid across branches + stays in the same box and so must get the same seed — the robustness that + lets the added noise cancel in the branch difference. We nudge the position + to a box centre first so a bounded jitter cannot cross a boundary. + """ + from shapepipe.modules.ngmix_package.ngmix import position_seed + + # Snap ra/dec to their box centres (1.5 arcsec into the 3-arcsec box), so a + # |jitter| < 1.5 arcsec cannot leave the box. + ra_c = (np.floor(ra * 3600 / 3) * 3 + 1.5) / 3600 + dec_c = (np.floor(dec * 3600 / 3) * 3 + 1.5) / 3600 + base = position_seed(ra_c, dec_c, ccd) + jittered = position_seed( + ra_c + jitter_arcsec / 3600, dec_c + jitter_arcsec / 3600, ccd + ) + assert base == jittered + + +@given(_ra, _dec, _ccd) +def test_position_seed_distinguishes_adjacent_boxes(ra, dec, ccd): + """Neighbouring 3-arcsec boxes get different seeds (in RA and in Dec). + + Both positions are snapped to box centres and the neighbour is one full box + (3 arcsec) away, so the shift genuinely lands in the next box rather than + grazing a boundary where ``3 / 3600``'s float error could leave ``floor`` + in the same box. + """ + from shapepipe.modules.ngmix_package.ngmix import position_seed + + ra_c = (np.floor(ra * 3600 / 3) * 3 + 1.5) / 3600 + dec_c = (np.floor(dec * 3600 / 3) * 3 + 1.5) / 3600 + base = position_seed(ra_c, dec_c, ccd) + assert position_seed(ra_c + 3 / 3600, dec_c, ccd) != base + assert position_seed(ra_c, dec_c + 3 / 3600, ccd) != base + + +def test_position_seed_resolves_fabian_anti_diagonal_collision(): + """The Cantor pairing avoids the ``box_x + box_y`` anti-diagonal collision. + + Fabian's #796 formula sums the box indices, so boxes ``(10, 20)`` and + ``(11, 19)`` collide (both -> 30) and two distinct objects would share a + noise stream. The box centres below realise exactly that pair; our pairing + must separate them. + """ + from shapepipe.modules.ngmix_package.ngmix import position_seed + + # Box centres for boxes (box_x, box_y) at ccd=0: box_x = floor(ra*1200)+1. + ra_a, dec_a = (10 - 1 + 0.5) / 1200, (20 - 2 + 0.5) / 1200 + ra_b, dec_b = (11 - 1 + 0.5) / 1200, (19 - 2 + 0.5) / 1200 + assert position_seed(ra_a, dec_a, 0) != position_seed(ra_b, dec_b, 0) + + # The two PSF families carry distinct ellipticity AND size (shapepipe#749): # the original image PSF (-> *_psf_orig) and the metacal reconvolution kernel # (-> *_psf_reconv) are independent fits, so a regression to the old aliasing