Seed fix based on position for Pujol noise cancellation for image simulations

[fabianhervaspeters]

When metacal adds a fixnoise image, it is useful for the differently sheared image branches that the same noise realization is used in image simulations. This helps cancel noise between the sheared branches. Here, “branches” refers to image branches, not metacal branches.
This is the approach I used in the previous iteration, and it would be useful to add it to a simulation-mode run so that sigma_m converges more quickly. At the moment, the noise is fixed only at the tile level.

shapepipe/src/shapepipe/modules/ngmix_package/ngmix.py

Line 332 in 4f871d4

seed = int(''.join(re.findall(r'\d+', self._file_number_string)))

The code below, which was written for the previous version, defines the seed from the sky position using a 3-arcsecond grid. This means that even if the sheared versions have different detection ordering, the added noise is attached to the position.

#This is for simulation and ensures the seed is fixed per object FHP
#The logic is that the seed is fixed in boxes in 3 arcsec so that  
#different detections receive the same seed, the +1 on the CCD is in case of exposure overlap.
ra,dec,ccd_n = coord_list[0]
box_x = np.floor((ra*3600) / 3)+(ccd_n+1)
box_y = np.floor((dec*3600) / 3)+(ccd_n+2)
unique_seed = box_x  + box_y
np.random.seed(int(unique_seed))

In the new version, the self._rng should now be fixed inside the object loop for that to work

shapepipe/src/shapepipe/modules/ngmix_package/ngmix.py

Line 665 in 4f871d4

for i_tile, obj_id in enumerate(tile_cat.obj_id):

.

[cailmdaley]

Implemented in #803 (opt-in SEED_FROM_POSITION, default off; box math exactly as above, first-epoch coords). Two things surfaced during implementation that are worth knowing even if the PR evolves:

1. The prior has to be rebuilt from the per-object RNG, not just self._rng. In ngmix 2.x the guesser's initial guess comes from prior.sample(), which draws from the RNG the prior was constructed with — and the prior is built once per tile. So fixing self._rng inside the object loop is not enough by itself: the guess would stay on the shared tile stream and the cross-branch cancellation would silently fail (no error, just no cancellation). ngmix: position-seeded per-object RNG for Pujol noise cancellation (#796) #803 rebuilds the prior per object when the knob is on; verified byte-identical fits for the same object under different tile contexts and processing orders.

2. box_x + box_y collides along anti-diagonals — e.g. boxes (10,20) and (11,19) both give 30, so two distinct objects would share one noise stream. ngmix: position-seeded per-object RNG for Pujol noise cancellation (#796) #803 keeps your box math but combines the two indices with a Cantor pairing (collision-free) folded into the valid seed range.

With the knob off the module is byte-identical to current develop. The PR expects to rebase after #766, and the image-sims configs there should then set SEED_FROM_POSITION = True.

— Claude Fable on behalf of Cail

[fabianhervaspeters]

The way I compute the m-bias is, not sure where to put it

def sidney_style_selection_shear(gal_dic, m_dic):
    n_bins_des = 10
    npatch = 100
    key_ra = "X_WORLD"
    key_dec = "Y_WORLD"
    units = "deg"
    weight_type = "des"

    gal_pairs = [
        ("1p2z", "1m2z", "1"),
        ("1z2p", "1z2m", "2"),
    ]

    def safe_weights(weights):
        return np.nan_to_num(weights, nan=0)

    def weighted_masked_average(values, mask, weights):
        return np.average(values[mask], weights=safe_weights(weights[mask]))

    def shear_response(m, weights):
        r11 = (
            weighted_masked_average(m.p1["g1"], m.mask_dict["p1"], weights)
            - weighted_masked_average(m.m1["g1"], m.mask_dict["m1"], weights)
        ) / 0.02

        r22 = (
            weighted_masked_average(m.p2["g2"], m.mask_dict["p2"], weights)
            - weighted_masked_average(m.m2["g2"], m.mask_dict["m2"], weights)
        ) / 0.02

        return r11, r22

    out_dic = {}

    for gal_ind1, gal_ind2, comp in gal_pairs:
        common_tiles = set(gal_dic[gal_ind1]["TILE_ID"]) & set(gal_dic[gal_ind2]["TILE_ID"])
        common_tiles = np.array(list(common_tiles), dtype=float)

        tile_list = np.array([
            tile for tile in common_tiles
            if np.sum(gal_dic[gal_ind1]["TILE_ID"] == tile) > 5
        ])

        data1 = treecorr.Catalog(
            ra=gal_dic[gal_ind1][key_ra],
            dec=gal_dic[gal_ind1][key_dec],
            ra_units=units,
            dec_units=units,
            npatch=npatch,
        )

        data2 = treecorr.Catalog(
            ra=gal_dic[gal_ind2][key_ra],
            dec=gal_dic[gal_ind2][key_dec],
            ra_units=units,
            dec_units=units,
            npatch=npatch,
            patch_centers=data1.patch_centers,
        )

        gal_dic1 = gal_dic[gal_ind1]
        gal_dic2 = gal_dic[gal_ind2]

        if weight_type == "des":
            weights_p_all = get_w_des_fix(
                gal_dic1,
                m_dic[gal_ind1].R11,
                m_dic[gal_ind1].R22,
                n_bins_des,
            )
            weights_m_all = get_w_des_fix(
                gal_dic2,
                m_dic[gal_ind2].R11,
                m_dic[gal_ind2].R22,
                n_bins_des,
            )
        else:
            weights_p_all = m_dic[gal_ind1].ns["w"][m_dic[gal_ind1].mask_dict["ns"]]
            weights_m_all = m_dic[gal_ind2].ns["w"][m_dic[gal_ind2].mask_dict["ns"]]

        r11_1, r22_1 = shear_response(m_dic[gal_ind1], weights_p_all)
        r11_2, r22_2 = shear_response(m_dic[gal_ind2], weights_m_all)

        cat_e1_list = []
        cat_e2_list = []

        for patch in np.unique(data1.patch):
            cond1 = data1.patch == patch
            cond2 = data2.patch == patch

            weights1 = safe_weights(weights_p_all[~cond1])
            weights2 = safe_weights(weights_m_all[~cond2])

            mean_e1_1 = np.average(
                gal_dic1["NGMIX_ELL_NOSHEAR"][:, 0][~cond1],
                weights=weights1,
            )
            mean_e1_2 = np.average(
                gal_dic2["NGMIX_ELL_NOSHEAR"][:, 0][~cond2],
                weights=weights2,
            )

            mean_e2_1 = np.average(
                gal_dic1["NGMIX_ELL_NOSHEAR"][:, 1][~cond1],
                weights=weights1,
            )
            mean_e2_2 = np.average(
                gal_dic2["NGMIX_ELL_NOSHEAR"][:, 1][~cond2],
                weights=weights2,
            )

            cat_e1_list.append((mean_e1_1 - mean_e1_2) / (r11_1 + r11_2) / 0.025)
            cat_e2_list.append((mean_e2_1 - mean_e2_2) / (r22_1 + r22_2) / 0.025)

        cat_e1_list = np.array(cat_e1_list)
        cat_e2_list = np.array(cat_e2_list)

        cat_e1_mean = np.mean(cat_e1_list)
        cat_e2_mean = np.mean(cat_e2_list)

        jack_error_e1 = np.sqrt(np.sum((cat_e1_list - cat_e1_mean) ** 2))
        jack_error_e2 = np.sqrt(np.sum((cat_e2_list - cat_e2_mean) ** 2))

        out_dic[comp] = (
            cat_e1_mean - 1,
            cat_e2_mean - 1,
            jack_error_e1,
            jack_error_e2,
        )

    return out_dic

[cailmdaley]

thank you Fabian!