From cee5509f117654d03be9d0376595fa5b8a874d8f Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 02:11:54 +0200 Subject: [PATCH] ngmix: position-seeded per-object RNG for Pujol noise cancellation (#796) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add an opt-in SEED_FROM_POSITION knob to the ngmix module (default False). When on, the object loop replaces the tile-level RNG with a per-object RandomState seeded from the object's sky position, and rebuilds the prior from that same RNG, so every stochastic step of the fit is a deterministic function of (ra, dec, ccd). Why: image-sim m-bias uses the Pujol estimator — the same scene under different applied shears, with shear response read from the branch difference of the SAME objects. Metacal's fixnoise adds a counter-noise realisation from an RNG seeded per tile, so 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 per object from position 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. Design and box math are Fabian's issue #796. The prior must be rebuilt per object, not just the RNG: the ngmix guesser draws its initial guess via prior.sample(), which consumes the RNG the prior was CONSTRUCTED with (ngmix 2.4.0 guessers.py / joint_prior.py). A per-object RNG alone would leave the guess drawing from the shared tile stream and break cancellation. The fixnoise counter-noise itself rides on the observation's attached noise image (use_noise_image=True), which prepare_ngmix_weights draws from the per-object RNG — so no fresh draw from the bootstrapper RNG is needed. Deviation from #796: Fabian combines the 3-arcsec position boxes as box_x + box_y, which collides along anti-diagonals — boxes (10,20) and (11,19) both give seed 30, so two distinct objects share a noise stream. This replaces the sum with a Cantor pairing (a bijection on the non-negative integers) after a zig-zag fold that maps signed box indices (Dec can be negative) onto non-negative ones, reduced mod 2**32 into the valid RandomState seed range. The box math is kept exactly as his: floor(ra*3600/3)+(ccd+1), floor(dec*3600/3)+(ccd+2), first-epoch coordinates (coord_list[0] convention), so all epochs of one object share the stream. Off by default => production path is byte-identical to develop (verified: metacal g for all five types matches unmodified develop bit-for-bit; the full fast suite passes, 326 tests). New property tests cover the seed's range, sub-box jitter invariance, adjacent-box distinctness, and the resolved anti-diagonal collision. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01J9qEVRCWEekT5ACNB3rjkq --- example/cfis/config_tile_Ng_template.ini | 9 ++ src/shapepipe/modules/ngmix_package/ngmix.py | 104 ++++++++++++++++++- src/shapepipe/modules/ngmix_runner.py | 12 +++ tests/module/test_ngmix.py | 78 ++++++++++++++ 4 files changed, 201 insertions(+), 2 deletions(-) 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