implemented weighted shear response

Author: martinkilbinger

config/calibration/mask_v1.X.6.yaml

@@ -0,0 +1,110 @@
+# Config file for masking and calibration.
+# Standard cuts without halo masks, type v1.X.6.
+
+# General parameters (can also given on command line)
+params:
+  input_path: unions_shapepipe_comprehensive_struc_2024_v1.X.c.hdf5
+  cmatrices: True
+  sky_regions: False
+  verbose: True
+
+# Masks
+## Using columns in 'dat' group (ShapePipe flags)
+dat:
+  # SExtractor flags
+  - col_name: FLAGS
+    label: SE FLAGS
+    kind: equal
+    value: 0
+
+  # Duplicate objects
+  - col_name: overlap
+    label: tile overlap
+    kind: equal
+    value: True
+
+  # ShapePipe flags
+  - col_name: IMAFLAGS_ISO
+    label: SP mask
+    kind: equal
+    value: 0
+
+  # Number of epochs
+  - col_name: N_EPOCH
+    label: r"$n_{\rm epoch}$"
+    kind: greater_equal
+    value: 2
+
+  # Magnitude range
+  - col_name: mag
+    label: mag range
+    kind: range
+    value: [15, 30]
+
+  # ngmix flags
+  - col_name: NGMIX_MOM_FAIL
+    label: "ngmix moments failure"
+    kind: equal
+    value: 0
+
+  # invalid PSF ellipticities
+  - col_name: NGMIX_ELL_PSFo_NOSHEAR_0
+    label: "bad PSF ellipticity comp 1"
+    kind: not_equal
+    value: -10
+  - col_name: NGMIX_ELL_PSFo_NOSHEAR_1
+    label: "bad PSF ellipticity comp 2"
+    kind: not_equal
+    value: -10
+
+## Using columns in 'dat_ext' group (post-processing flags)
+dat_ext:
+
+  # Stars
+  - col_name: 4_Stars
+    label: "Stars"
+    kind: equal
+    value: False
+
+  # Manual mask
+  - col_name: 8_Manual
+    label: "manual mask"
+    kind: equal
+    value: False
+
+  # r-band footprint
+  - col_name: 64_r
+    label: "r-band imaging"
+    kind: equal
+    value: False
+
+  # Maximask
+  - col_name: 1024_Maximask
+    label: "maximask"
+    kind: equal
+    value: False
+
+  # Rough pointing coverage
+  - col_name: npoint3
+    label: r"$n_{\rm pointing}$"
+    kind: greater_equal
+    value: 3
+
+# Metacal parameters
+metacal:
+  # Ellipticity dispersion
+  sigma_eps_prior: 0.34
+
+  # Signal-to-noise range
+  gal_snr_min: 10
+  gal_snr_max: 1000
+
+  # Relative-size (hlr / hlr_psf) range
+  gal_rel_size_min: 0.707
+  gal_rel_size_max: 3
+
+  # Correct relative size for ellipticity?
+  gal_size_corr_ell: False
+
+  # Weight for global response matrix, None for unweighted mean
+  global_R_weight: w

notebooks/calibrate_comprehensive_cat.py

@@ -38,8 +38,8 @@
 dat, dat_ext = obj.read_cat(load_into_memory=False)
 
 # %%
-#n_test = -1
-n_test = 100000
+n_test = -1
+#n_test = 100000
 if n_test > 0:
     print(f"MKDEBUG testing only first {n_test} objects")
     dat = dat[:n_test]
@@ -234,7 +234,7 @@
     g1_uncal=g_uncorr[0],
     g2_uncal=g_uncorr[1],
     R=gal_metacal.R,
-    R_shear=R_shear,
+    R_shear=gal_metacal.R_shear_global,
     R_select=gal_metacal.R_selection,
     c=c,
     c_err=c_err,

sp_validation/basic.py

@@ -44,8 +44,6 @@ class metacal:
         masking type, one in 'gal', 'gal_mom', 'star'
     step : float, optional, default=0.01
         step h in finite differences
-    stat_operator : optional, default=np.mean
-        summary statistic for response matrices
     prefix : string, optional, default='NGMIX'
         to specify columns in input catalogue
     snr_min : float, optional, default=10
@@ -57,6 +55,9 @@ class metacal:
     rel_size_max : float, optional, default=3.0
         relative size maximum
     size_corr_ell : bool, optional, default=True
+    global_R_weight : str, optional,
+        weight column name for global response matrix; default is ``None``
+        (unweighted mean)
     sigma_eps : float, optional
         ellipticity dispersion (one component) for computation
         of weights; default is 0.34
@@ -74,21 +75,20 @@ def __init__(
         mask,
         masking_type='gal',
         step=0.01,
-        stat_operator=np.mean,
         prefix='NGMIX',
         snr_min=10,
         snr_max=500,
         rel_size_min=0.5,
         rel_size_max=3.0,
         size_corr_ell=True,
+        global_R_weight=None,
         sigma_eps=0.34,
         col_2d=True,
         verbose=False,
     ):
 
         self._masking_type = masking_type
         self._step = step
-        self._stat_operator = stat_operator
 
         # Cuts
         self._snr_min = snr_min
@@ -103,6 +103,8 @@ def __init__(
                 + f'rel_size_max={rel_size_max}, '
                 + f'size_corr_ell={size_corr_ell}'
             )
+            
+        self._global_R_weight = global_R_weight
 
         self._sigma_eps = sigma_eps
         self._col_2d = col_2d
@@ -204,7 +206,10 @@ def _read_data_ngmix(self, masked_data, m1, p1, m2, p2, ns):
         self._n_after_gal_mask = len(dict_tmp['flag'])
         if self._verbose:
             print(f"Number of objects on metacal input = {self._n_input}")
-            print(f"Number of objects after galaxy selection masking = {self._n_after_gal_mask}")
+            print(
+                "Number of objects after galaxy selection masking ="
+                + f" {self._n_after_gal_mask}"
+            )
 
         return m1, p1, m2, p2, ns
 
@@ -503,20 +508,20 @@ def _selection_response(self):
         h2 = 2 * self._step
 
         self.R11_s = (
-            self._stat_operator(self.ns['g1'][ma_p1])
-            - self._stat_operator(self.ns['g1'][ma_m1])
+            np.mean(self.ns['g1'][ma_p1])
+            - np.mean(self.ns['g1'][ma_m1])
         ) / h2
         self.R22_s = sign * (
-            self._stat_operator(self.ns['g2'][ma_p2])
-            - self._stat_operator(self.ns['g2'][ma_m2])
+            np.mean(self.ns['g2'][ma_p2])
+            - np.mean(self.ns['g2'][ma_m2])
         ) / h2
         self.R12_s = (
-            self._stat_operator(self.ns['g1'][ma_p2])
-            - self._stat_operator(self.ns['g1'][ma_m2])
+            np.mean(self.ns['g1'][ma_p2])
+            - np.mean(self.ns['g1'][ma_m2])
         ) / h2
         self.R21_s = (
-            self._stat_operator(self.ns['g2'][ma_p1])
-            - self._stat_operator(self.ns['g2'][ma_m1])
+            np.mean(self.ns['g2'][ma_p1])
+            - np.mean(self.ns['g2'][ma_m1])
         ) / h2
 
         self.R_selection = np.array([
@@ -530,8 +535,16 @@ def _total_response(self):
         ...
 
         """
-        R_ws = self._stat_operator(self.R_shear, 2)
-        self.R = R_ws + self.R_selection
+        if self._global_R_weight is None or self._global_R_weight == "None":
+            print("Computing unweighted response")
+            self.R_shear_global = np.mean(self.R_shear, axis=2)
+        else:
+            print("Computing response weighted by", self._global_R_weight)
+            # Get weights of masked no-shear objects
+            weights = self.ns[self._global_R_weight][self.mask_dict['ns']]
+            self.R_shear_global = np.average(self.R_shear, axis=2, weights=weights)
+
+        self.R = self.R_shear_global + self.R_selection
 
     def _return():
         """Add docstring.

sp_validation/calibration.py

@@ -28,7 +28,8 @@
 def get_calibrated_quantities(gal_metacal, shape_method='ngmix'):
     """Get Calibrated Quantities.
 
-    Return catalogue quantities for objects calibrated for multiplicative bias.
+    Return catalogue quantities for objects calibrated for multiplicative
+    bias.
 
     Parameters
     ----------
@@ -97,7 +98,9 @@ def get_calibrated_m_c(gal_metacal, shape_method='ngmix'):
 
     """
     # Get m-calibrated quantities
-    g_corr, g_uncorr, w, mask_metacal = get_calibrated_quantities(gal_metacal)
+    g_corr, g_uncorr, w, mask_metacal = get_calibrated_quantities(
+        gal_metacal
+    )
 
     # Additive bias
     c = np.zeros(2)
@@ -106,7 +109,8 @@ def get_calibrated_m_c(gal_metacal, shape_method='ngmix'):
     for comp in (0, 1):
         c[comp] = np.mean(g_uncorr[comp])
 
-        # MKDEBUG TODO: Use std of mean instead, which is consistent with jackknife
+        # MKDEBUG TODO: Use std of mean instead,
+        # which is consistent with jackknife
         c_err[comp] = np.std(g_uncorr[comp])
 
     # Shear estimate corrected for additive bias
@@ -269,19 +273,31 @@ def get_w_des(cat_gal, num_bins):
 
     for i in range(num_bins):
         for j in range(num_bins):
-            bin_mask = (df_gal['snr_log_bins'] == i) & (df_gal['size_ratio_log_bins'] == j)
+            bin_mask = (
+                (df_gal['snr_log_bins'] == i)
+                & (df_gal['size_ratio_log_bins'] == j)
+            )
             ngal = np.sum(bin_mask)
-            shape_noise = 0.5*(np.sum(df_gal[bin_mask]['e1_uncal']**2)/ngal+np.sum(df_gal[bin_mask]['e2_uncal']**2)/ngal)
-            response = 0.5*(np.average(df_gal[bin_mask]['R_g11'])+np.average(df_gal[bin_mask]['R_g22']))
-            df_gal.loc[bin_mask, 'w_des'] = response**2/shape_noise     
+            if ngal == 0:
+                print(f"Zero galaxies in snr/size_ratio bin ({i},{j})")
+            shape_noise = 0.5 * (
+                np.sum(df_gal[bin_mask]['e1_uncal'] ** 2) / ngal
+                + np.sum(df_gal[bin_mask]['e2_uncal'] ** 2) / ngal
+            )
+            response = 0.5 * (
+                np.average(df_gal[bin_mask]['R_g11'])
+                + np.average(df_gal[bin_mask]['R_g22'])
+            )
+            df_gal.loc[bin_mask, 'w_des'] = response ** 2 / shape_noise     
 
     return np.array(df_gal['w_des'])
 
 
 def get_alpha_leakage_per_object(cat_gal, num_bins, weight_type='des'):
     """
     Compute the leakage per object (Li et al. 2024)
-    Return an array of leakage coefficients obtained by binning in SNR and size.
+    Return an array of leakage coefficients obtained by binning in
+    SNR and size.
 
     Parameters
     ----------
@@ -293,13 +309,18 @@ def get_alpha_leakage_per_object(cat_gal, num_bins, weight_type='des'):
     Returns
     -------
     alpha_1 : np.array
-        Array containing the correction coefficient for the PSF leakage per object for the first component.
+        Array containing the correction coefficient for the PSF leakage
+        per object for the first component.
     alpha_2 : np.array
-        Array containing the correction coefficient for the PSF leakage per object for the second component.
+        Array containing the correction coefficient for the PSF leakage
+        per object for the second component.
     """
     assert weight_type in ['des', 'iv'], "weight_type must be either 'des' or 'iv'"
-    #Compute the size ratio
-    size_ratio = cat_gal['NGMIX_Tpsf_NOSHEAR']/(cat_gal['NGMIX_T_NOSHEAR']+cat_gal['NGMIX_Tpsf_NOSHEAR'])
+    # Compute the size ratio
+    size_ratio = (
+        cat_gal['NGMIX_Tpsf_NOSHEAR']
+        / (cat_gal['NGMIX_T_NOSHEAR'] + cat_gal['NGMIX_Tpsf_NOSHEAR'])
+    )
 
     df_gal = pd.DataFrame(np.array([
         np.array(cat_gal['e1'], dtype=np.float64),
@@ -327,7 +348,12 @@ def get_alpha_leakage_per_object(cat_gal, num_bins, weight_type='des'):
         mask_binr = df_gal['bin_R'].values == ibin_r
 
         #bin in snr
-        df_gal.loc[mask_binr, 'bin_snr'] = pd.qcut(df_gal.loc[mask_binr, 'snr'], n_bins_snr, labels=False, retbins=False)
+        df_gal.loc[mask_binr, 'bin_snr'] = pd.qcut(
+            df_gal.loc[mask_binr, 'snr'],
+            n_bins_snr,
+            labels=False,
+            retbins=False
+        )
 
     #group by bin
     df_gal_grouped = df_gal.groupby(['bin_R', 'bin_snr'])