ngmix: position-seeded per-object RNG for Pujol noise cancellation (#796)

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

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,13 +395,19 @@ 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
 
         # Initiatlise random generator
         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
 

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

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