From 44baf31482d3625b7cfa71b545df5d52699072c4 Mon Sep 17 00:00:00 2001
From: Masaya Yamamoto <masaya.yamamoto@duke.edu>
Date: Mon, 4 Dec 2023 11:03:24 -0800
Subject: [PATCH 01/20] added functions to correct for extinction and v6 cuts

---
 des_y6utils/mdet.py | 207 +++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 205 insertions(+), 2 deletions(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index c8cf956..cbe39ee 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -7,6 +7,52 @@
 import healsparse
 
 
+def add_extinction_correction_columns(fold, fnew):
+    """This function changes and adds columns for extinction corrected magnitudes.
+    pgauss_band_flux_* are the columns that are the corrected fluxes. 
+    pgauss_band_flux_*_nodered are the columns that used to be pgauss_band_flux_*. 
+
+    Note: This function assumes you're working with fits files, not hdf5. 
+
+    Parameters
+    ----------
+    fold : original file path
+    fnew : new file path
+    
+    Returns
+    -------
+    None (saves the new file) 
+    """
+
+    import fitsio as fio
+    d = fio.read(fold)
+    bands = ['g', 'r', 'i', 'z']
+    fits = fio.FITS(fnew, 'rw')
+    
+    # compute dereddened fluxes and 
+    # replace the entries of pgauss_flux_band_* with the dereddened fluxes.
+    dustmap = _read_hsp_dustmap_local(
+                '/global/cfs/cdirs/des/y6-shear-catalogs/SFD_dust_4096.hsp'
+                )
+    dered = dustmap.get_values_pos(d["ra"], d["dec"])
+    flux_og = []
+    for ii,b in enumerate(bands):
+        flux = np.copy(d["pgauss_band_flux_"+b])
+        flux_og.append(flux)
+        mag_ = _compute_asinh_dered_mag(d["pgauss_band_flux_"+b], ii, dered)
+        flux_ = _compute_asinh_flux(mag_, ii)
+        d["pgauss_band_flux_"+b] = flux_
+    fits.write(d)
+
+    # make _nodered array with pgauss_band_flux_* entries, and add them to fits. 
+    for ii,b in enumerate(bands):    
+        fits[-1].insert_column('pgauss_band_flux_'+b+'_nodered', flux_og[ii])
+
+    fits.close()
+
+    return None
+
+
 def make_mdet_cuts(data, version, verbose=False):
     """A function to make the standard metadetection cuts.
 
@@ -40,6 +86,8 @@ def make_mdet_cuts(data, version, verbose=False):
         return _make_mdet_cuts_v4(data, verbose=verbose)
     elif str(version) == "5":
         return _make_mdet_cuts_v5(data, verbose=verbose)
+    elif str(version) == "6":
+        return _make_mdet_cuts_v6(data, verbose=verbose)
     else:
         raise ValueError("the mdet cut version '%r' is not recognized!" % version)
 
@@ -52,7 +100,7 @@ def _make_mdet_cuts_gauss(
     min_t_ratio=0.5,
     max_mfrac=0.1,
     max_s2n=np.inf,
-    mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v3.fits",
+    mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v1.hsp",
     max_t=100.0,
 ):
     """
@@ -61,7 +109,7 @@ def _make_mdet_cuts_gauss(
     msk : np.ndarray of bool
         A boolean array with the cuts. To cut the data, use `data[msk]`.
     """
-    msk = _make_mdet_cuts_raw_v345(
+    msk = _make_mdet_cuts_raw_v6(
         data,
         verbose=verbose,
         min_s2n=min_s2n,
@@ -258,6 +306,35 @@ def _make_mdet_cuts_v5(d, verbose=False):
 
     return msk
 
+def _make_mdet_cuts_v6(d, verbose=False):
+
+    msk = _make_mdet_cuts_raw_v6(
+        d,
+        verbose=verbose,
+        min_s2n=10,
+        n_terr=0,
+        min_t_ratio=0.5,
+        max_mfrac=0.1,
+        max_s2n=np.inf,
+        max_t=20.0,
+    )
+
+    size_sizeerr = (d['gauss_T_ratio']*d['gauss_psf_T']) * d['gauss_T_err']
+    size_s2n = (d['gauss_T_ratio']*d['gauss_psf_T']) / d['gauss_T_err']
+    msk_superspreader = ((size_sizeerr > 1) & (size_s2n < 10))
+    msk &= ~msk_superspreader
+
+    # apply the mask
+    hmap = _read_hsp_mask(
+        "y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v1.hsp"
+    )
+    in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
+    msk &= in_footprint
+    if verbose:
+        print("did mask cuts", np.sum(msk))
+
+    return msk
+
 
 def _compute_asinh_mags(flux, i):
     """This function and coefficients are from from Eli. Ask him.
@@ -291,6 +368,41 @@ def _compute_asinh_mags(flux, i):
     return mag
 
 
+def _compute_asinh_dered_mag(flux, i, map):
+
+    dered_fac = [3.186, 2.140, 1.196, 1.048]
+    mag = _compute_asinh_mags(flux, i)
+    dered_mag = mag - dered_fac[i]*map
+    
+    return dered_mag
+
+
+def _compute_asinh_flux(mag, i):
+
+    """This function does the inverse process of compute_asinh_mag.
+
+    Parameters
+    ----------
+    mag : float or np.ndarray
+        The magnitude.
+    i : int
+        The index of the band in griz (i.e., 0 for g, 1 for r, 2 for i, 3 for z).
+
+    Returns
+    -------
+    flux : float or np.ndarray
+        The asinh flux for the mag.
+    """
+
+    zp = 30.0
+    b_array = np.array([3.27e-12, 4.83e-12, 6.0e-12, 9.0e-12])
+    bscale = np.array(b_array) * 10.**(zp / 2.5)
+    A = (2.5 * np.log10(1.0 / b_array[i]) - mag) * 0.4 * np.log(10)
+    flux = 2*bscale[i] * np.sinh(A)
+
+    return flux
+
+
 @lru_cache
 def _read_hsp_mask(fname):
     mpth = _get_mask_path(fname)
@@ -301,6 +413,10 @@ def _read_hsp_mask(fname):
 def _read_hsp_mask_local(fname):
     return healsparse.HealSparseMap.read(fname)
 
+@lru_cache
+def _read_hsp_dustmap_local(fname):
+    return healsparse.HealSparseMap.read(fname)
+
 
 def _get_mask_path(fname):
     # get or make meds dir
@@ -528,3 +644,90 @@ def _make_mdet_cuts_raw_v345(
         print("did mdet cuts", np.sum(msk))
 
     return msk
+
+
+def _make_mdet_cuts_raw_v6(
+    d,
+    *,
+    min_s2n,
+    n_terr,
+    min_t_ratio,
+    max_mfrac,
+    max_s2n,
+    max_t,
+    verbose=False,
+):
+    """The raw v6 cuts come from analysis in Fall of 2023.
+
+      - We implemented i-band magnitude cut at 24.5. TO-DO: still need to test.
+      - We implemented galaxy size cut at T=20. 
+      - We updated a cut in the pgauss T-Terr plane. 
+    """
+
+    mag_g = _compute_asinh_mags(d["pgauss_band_flux_g"], 0)
+    mag_r = _compute_asinh_mags(d["pgauss_band_flux_r"], 1)
+    mag_i = _compute_asinh_mags(d["pgauss_band_flux_i"], 2)
+    mag_z = _compute_asinh_mags(d["pgauss_band_flux_z"], 3)
+
+    gmr = mag_g - mag_r
+    rmi = mag_r - mag_i
+    imz = mag_i - mag_z
+
+    msk = np.ones(d.shape[0]).astype(bool)
+
+    potential_flag_columns = [
+        "psfrec_flags",
+        "gauss_flags",
+        "pgauss_T_flags",
+        "pgauss_band_flux_flags_g",
+        "pgauss_band_flux_flags_r",
+        "pgauss_band_flux_flags_i",
+        "pgauss_band_flux_flags_z",
+        "mask_flags",
+    ]
+    for col in potential_flag_columns:
+        if col in d.dtype.names:
+            msk &= (d[col] == 0)
+            if verbose:
+                print("did cut " + col, np.sum(msk))
+
+    if "shear_bands" in d.dtype.names:
+        msk &= (d["shear_bands"] == "123")
+        if verbose:
+            print("did cut shear_bands", np.sum(msk))
+
+    if "pgauss_s2n" in d.dtype.names:
+        msk &= (d["pgauss_s2n"] > 5)
+        if verbose:
+            print("did cut pgauss_s2n", np.sum(msk))
+
+    # now do the rest
+    msk &= (
+        (d["gauss_s2n"] > min_s2n)
+        & (d["gauss_s2n"] < max_s2n)
+        & (d["mfrac"] < max_mfrac)
+        & (np.abs(gmr) < 5)
+        & (np.abs(rmi) < 5)
+        & (np.abs(imz) < 5)
+        & np.isfinite(mag_g)
+        & np.isfinite(mag_r)
+        & np.isfinite(mag_i)
+        & np.isfinite(mag_z)
+        & (mag_g < 26.5)
+        & (mag_r < 26.5)
+        & (mag_i < 24.5)
+        & (mag_z < 25.6)
+        & (d["pgauss_T"] < (1.6 - 3.1*d["pgauss_T_err"]))
+        & (
+            d["gauss_T_ratio"] >= np.maximum(
+                min_t_ratio,
+                (n_terr*d["gauss_T_err"]/d["gauss_psf_T"])
+            )
+        )
+        & ((d["gauss_T_ratio"] * d["gauss_psf_T"]) < max_t)
+    )
+
+    if verbose:
+        print("did mdet cuts", np.sum(msk))
+
+    return msk
\ No newline at end of file

From fc0c31b8f1b809aa805011d8d83d67c9b45eef58 Mon Sep 17 00:00:00 2001
From: Masaya Yamamoto <masaya.yamamoto@duke.edu>
Date: Mon, 4 Dec 2023 11:10:21 -0800
Subject: [PATCH 02/20] flake8

---
 des_y6utils/mdet.py | 71 +++++++++++++++++++++++----------------------
 1 file changed, 36 insertions(+), 35 deletions(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index cbe39ee..38fa14c 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -9,44 +9,44 @@
 
 def add_extinction_correction_columns(fold, fnew):
     """This function changes and adds columns for extinction corrected magnitudes.
-    pgauss_band_flux_* are the columns that are the corrected fluxes. 
-    pgauss_band_flux_*_nodered are the columns that used to be pgauss_band_flux_*. 
+    pgauss_band_flux_* are the columns that are the corrected fluxes.
+    pgauss_band_flux_*_nodered are the columns that used to be pgauss_band_flux_*.
 
-    Note: This function assumes you're working with fits files, not hdf5. 
+    Note: This function assumes you're working with fits files, not hdf5.
 
     Parameters
     ----------
     fold : original file path
     fnew : new file path
-    
+
     Returns
     -------
-    None (saves the new file) 
+    None (saves the new file)
     """
 
     import fitsio as fio
     d = fio.read(fold)
     bands = ['g', 'r', 'i', 'z']
     fits = fio.FITS(fnew, 'rw')
-    
-    # compute dereddened fluxes and 
+
+    # compute dereddened fluxes and
     # replace the entries of pgauss_flux_band_* with the dereddened fluxes.
     dustmap = _read_hsp_dustmap_local(
-                '/global/cfs/cdirs/des/y6-shear-catalogs/SFD_dust_4096.hsp'
-                )
+        '/global/cfs/cdirs/des/y6-shear-catalogs/SFD_dust_4096.hsp'
+    )
     dered = dustmap.get_values_pos(d["ra"], d["dec"])
     flux_og = []
-    for ii,b in enumerate(bands):
-        flux = np.copy(d["pgauss_band_flux_"+b])
+    for ii, b in enumerate(bands):
+        flux = np.copy(d["pgauss_band_flux_" + b])
         flux_og.append(flux)
-        mag_ = _compute_asinh_dered_mag(d["pgauss_band_flux_"+b], ii, dered)
+        mag_ = _compute_asinh_dered_mag(d["pgauss_band_flux_" + b], ii, dered)
         flux_ = _compute_asinh_flux(mag_, ii)
-        d["pgauss_band_flux_"+b] = flux_
+        d["pgauss_band_flux_" + b] = flux_
     fits.write(d)
 
-    # make _nodered array with pgauss_band_flux_* entries, and add them to fits. 
-    for ii,b in enumerate(bands):    
-        fits[-1].insert_column('pgauss_band_flux_'+b+'_nodered', flux_og[ii])
+    # make _nodered array with pgauss_band_flux_* entries, and add them to fits.
+    for ii, b in enumerate(bands):
+        fits[-1].insert_column('pgauss_band_flux_' + b + '_nodered', flux_og[ii])
 
     fits.close()
 
@@ -260,8 +260,8 @@ def _make_mdet_cuts_v4(d, verbose=False):
         max_t=np.inf,
     )
 
-    size_sizeerr = (d['gauss_T_ratio']*d['gauss_psf_T']) * d['gauss_T_err']
-    size_s2n = (d['gauss_T_ratio']*d['gauss_psf_T']) / d['gauss_T_err']
+    size_sizeerr = (d['gauss_T_ratio'] * d['gauss_psf_T']) * d['gauss_T_err']
+    size_s2n = (d['gauss_T_ratio'] * d['gauss_psf_T']) / d['gauss_T_err']
     msk_superspreader = ((size_sizeerr > 1) & (size_s2n < 10))
     msk &= ~msk_superspreader
 
@@ -290,8 +290,8 @@ def _make_mdet_cuts_v5(d, verbose=False):
         max_t=np.inf,
     )
 
-    size_sizeerr = (d['gauss_T_ratio']*d['gauss_psf_T']) * d['gauss_T_err']
-    size_s2n = (d['gauss_T_ratio']*d['gauss_psf_T']) / d['gauss_T_err']
+    size_sizeerr = (d['gauss_T_ratio'] * d['gauss_psf_T']) * d['gauss_T_err']
+    size_s2n = (d['gauss_T_ratio'] * d['gauss_psf_T']) / d['gauss_T_err']
     msk_superspreader = ((size_sizeerr > 1) & (size_s2n < 10))
     msk &= ~msk_superspreader
 
@@ -306,6 +306,7 @@ def _make_mdet_cuts_v5(d, verbose=False):
 
     return msk
 
+
 def _make_mdet_cuts_v6(d, verbose=False):
 
     msk = _make_mdet_cuts_raw_v6(
@@ -319,8 +320,8 @@ def _make_mdet_cuts_v6(d, verbose=False):
         max_t=20.0,
     )
 
-    size_sizeerr = (d['gauss_T_ratio']*d['gauss_psf_T']) * d['gauss_T_err']
-    size_s2n = (d['gauss_T_ratio']*d['gauss_psf_T']) / d['gauss_T_err']
+    size_sizeerr = (d['gauss_T_ratio'] * d['gauss_psf_T']) * d['gauss_T_err']
+    size_s2n = (d['gauss_T_ratio'] * d['gauss_psf_T']) / d['gauss_T_err']
     msk_superspreader = ((size_sizeerr > 1) & (size_s2n < 10))
     msk &= ~msk_superspreader
 
@@ -372,13 +373,12 @@ def _compute_asinh_dered_mag(flux, i, map):
 
     dered_fac = [3.186, 2.140, 1.196, 1.048]
     mag = _compute_asinh_mags(flux, i)
-    dered_mag = mag - dered_fac[i]*map
-    
+    dered_mag = mag - dered_fac[i] * map
+
     return dered_mag
 
 
 def _compute_asinh_flux(mag, i):
-
     """This function does the inverse process of compute_asinh_mag.
 
     Parameters
@@ -398,7 +398,7 @@ def _compute_asinh_flux(mag, i):
     b_array = np.array([3.27e-12, 4.83e-12, 6.0e-12, 9.0e-12])
     bscale = np.array(b_array) * 10.**(zp / 2.5)
     A = (2.5 * np.log10(1.0 / b_array[i]) - mag) * 0.4 * np.log(10)
-    flux = 2*bscale[i] * np.sinh(A)
+    flux = 2 * bscale[i] * np.sinh(A)
 
     return flux
 
@@ -413,6 +413,7 @@ def _read_hsp_mask(fname):
 def _read_hsp_mask_local(fname):
     return healsparse.HealSparseMap.read(fname)
 
+
 @lru_cache
 def _read_hsp_dustmap_local(fname):
     return healsparse.HealSparseMap.read(fname)
@@ -547,11 +548,11 @@ def _make_mdet_cuts_raw_v12(
         & (mag_r < 26.5)
         & (mag_i < 26.2)
         & (mag_z < 25.6)
-        & (d["pgauss_T"] < (1.9 - 2.8*d["pgauss_T_err"]))
+        & (d["pgauss_T"] < (1.9 - 2.8 * d["pgauss_T_err"]))
         & (
             d["wmom_T_ratio"] >= np.maximum(
                 min_t_ratio,
-                (1.0 + n_terr*d["wmom_T_err"]/d["wmom_psf_T"])
+                (1.0 + n_terr * d["wmom_T_err"] / d["wmom_psf_T"])
             )
         )
     )
@@ -630,11 +631,11 @@ def _make_mdet_cuts_raw_v345(
         & (mag_r < 26.5)
         & (mag_i < 26.2)
         & (mag_z < 25.6)
-        & (d["pgauss_T"] < (1.2 - 3.1*d["pgauss_T_err"]))
+        & (d["pgauss_T"] < (1.2 - 3.1 * d["pgauss_T_err"]))
         & (
             d["gauss_T_ratio"] >= np.maximum(
                 min_t_ratio,
-                (n_terr*d["gauss_T_err"]/d["gauss_psf_T"])
+                (n_terr * d["gauss_T_err"] / d["gauss_psf_T"])
             )
         )
         & ((d["gauss_T_ratio"] * d["gauss_psf_T"]) < max_t)
@@ -660,8 +661,8 @@ def _make_mdet_cuts_raw_v6(
     """The raw v6 cuts come from analysis in Fall of 2023.
 
       - We implemented i-band magnitude cut at 24.5. TO-DO: still need to test.
-      - We implemented galaxy size cut at T=20. 
-      - We updated a cut in the pgauss T-Terr plane. 
+      - We implemented galaxy size cut at T=20.
+      - We updated a cut in the pgauss T-Terr plane.
     """
 
     mag_g = _compute_asinh_mags(d["pgauss_band_flux_g"], 0)
@@ -717,11 +718,11 @@ def _make_mdet_cuts_raw_v6(
         & (mag_r < 26.5)
         & (mag_i < 24.5)
         & (mag_z < 25.6)
-        & (d["pgauss_T"] < (1.6 - 3.1*d["pgauss_T_err"]))
+        & (d["pgauss_T"] < (1.6 - 3.1 * d["pgauss_T_err"]))
         & (
             d["gauss_T_ratio"] >= np.maximum(
                 min_t_ratio,
-                (n_terr*d["gauss_T_err"]/d["gauss_psf_T"])
+                (n_terr * d["gauss_T_err"] / d["gauss_psf_T"])
             )
         )
         & ((d["gauss_T_ratio"] * d["gauss_psf_T"]) < max_t)
@@ -730,4 +731,4 @@ def _make_mdet_cuts_raw_v6(
     if verbose:
         print("did mdet cuts", np.sum(msk))
 
-    return msk
\ No newline at end of file
+    return msk

From 6c07ac4406cc10aa9b95652cf85601f58571c38c Mon Sep 17 00:00:00 2001
From: Masaya Yamamoto <masaya.yamamoto@duke.edu>
Date: Wed, 6 Dec 2023 07:02:05 -0800
Subject: [PATCH 03/20] added unit tests for asinh mag and changed part to add
 columns

---
 des_y6utils/mdet.py                 | 59 +++++++++++++++++++----------
 des_y6utils/tests/test_asinh_mag.py | 17 +++++++++
 2 files changed, 55 insertions(+), 21 deletions(-)
 create mode 100644 des_y6utils/tests/test_asinh_mag.py

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index 38fa14c..0110f4e 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -7,7 +7,7 @@
 import healsparse
 
 
-def add_extinction_correction_columns(fold, fnew):
+def add_extinction_correction_columns(fold, fnew, fdust):
     """This function changes and adds columns for extinction corrected magnitudes.
     pgauss_band_flux_* are the columns that are the corrected fluxes.
     pgauss_band_flux_*_nodered are the columns that used to be pgauss_band_flux_*.
@@ -27,26 +27,37 @@ def add_extinction_correction_columns(fold, fnew):
     import fitsio as fio
     d = fio.read(fold)
     bands = ['g', 'r', 'i', 'z']
-    fits = fio.FITS(fnew, 'rw')
-
-    # compute dereddened fluxes and
-    # replace the entries of pgauss_flux_band_* with the dereddened fluxes.
-    dustmap = _read_hsp_dustmap_local(
-        '/global/cfs/cdirs/des/y6-shear-catalogs/SFD_dust_4096.hsp'
-    )
-    dered = dustmap.get_values_pos(d["ra"], d["dec"])
-    flux_og = []
-    for ii, b in enumerate(bands):
-        flux = np.copy(d["pgauss_band_flux_" + b])
-        flux_og.append(flux)
-        mag_ = _compute_asinh_dered_mag(d["pgauss_band_flux_" + b], ii, dered)
-        flux_ = _compute_asinh_flux(mag_, ii)
-        d["pgauss_band_flux_" + b] = flux_
-    fits.write(d)
-
-    # make _nodered array with pgauss_band_flux_* entries, and add them to fits.
-    for ii, b in enumerate(bands):
-        fits[-1].insert_column('pgauss_band_flux_' + b + '_nodered', flux_og[ii])
+    with fio.FITS(fnew, 'rw') as fits:
+
+        # compute dereddened fluxes and
+        # replace the entries of pgauss_flux_band_* with the dereddened fluxes.
+        if os.path.exists(fdust):
+            # fdust : '/global/cfs/cdirs/des/y6-shear-catalogs/SFD_dust_4096.hsp'
+            dustmap = _read_hsp_dustmap_local(fdust)
+        else:
+            dustmap = _read_hsp_mask(fdust)
+
+        dered = dustmap.get_values_pos(d["ra"], d["dec"])
+        flux_og = []
+        for ii, b in enumerate(bands):
+            flux = np.copy(d["pgauss_band_flux_" + b])
+            flux_og.append(flux)
+            mag_ = _compute_asinh_dered_mag(d["pgauss_band_flux_" + b], ii, dered)
+            flux_ = _compute_asinh_flux(mag_, ii)
+            d["pgauss_band_flux_" + b] = flux_
+
+        # make _nodered array with pgauss_band_flux_* entries, and add them to fits.
+        new_dt = np.dtype(d.dtype.descr + [('pgauss_band_flux_g_nodered', 'f8'),
+                                           ('pgauss_band_flux_r_nodered', 'f8'),
+                                           ('pgauss_band_flux_i_nodered', 'f8'),
+                                           ('pgauss_band_flux_z_nodered', 'f8'),])
+        d_ = np.zeros(d.shape, dtype=new_dt)
+        for col in d.dtype.names:
+            d_[col] = d[col]
+        for ii, b in enumerate(bands):
+            d_['pgauss_band_flux_' + b + '_nodered'] = flux_og[ii]
+
+        fits.write(d_)
 
     fits.close()
 
@@ -414,6 +425,12 @@ def _read_hsp_mask_local(fname):
     return healsparse.HealSparseMap.read(fname)
 
 
+@lru_cache
+def _read_hsp_dustmap(fname):
+    mpth = _get_mask_path(fname)
+    return healsparse.HealSparseMap.read(mpth)
+
+
 @lru_cache
 def _read_hsp_dustmap_local(fname):
     return healsparse.HealSparseMap.read(fname)
diff --git a/des_y6utils/tests/test_asinh_mag.py b/des_y6utils/tests/test_asinh_mag.py
new file mode 100644
index 0000000..9a60ea7
--- /dev/null
+++ b/des_y6utils/tests/test_asinh_mag.py
@@ -0,0 +1,17 @@
+from ..mdet import _compute_asinh_dered_mag, _compute_asinh_flux, _compute_asinh_mags
+
+
+def test_compute_asinh_flux():
+    # this test ensures the function computes magnitudes and fluxes correctly.
+    flux_input = 10000
+    mag_g = _compute_asinh_mags(flux_input, 0)
+    flux_g = _compute_asinh_flux(mag_g, 0)
+
+    assert flux_input == flux_g
+
+def test_compute_asinh_dered_mag():
+    flux_input = 10000
+    mag_g = _compute_asinh_mags(flux_input, 0)
+    dered_mag_g = _compute_asinh_dered_mag(flux_input, 0, 1/3.186)
+
+    assert mag_g == dered_mag_g
\ No newline at end of file

From 25ad3e4648888dab6141e23ecc2f26a706ce76b5 Mon Sep 17 00:00:00 2001
From: Masaya Yamamoto <masaya.yamamoto@duke.edu>
Date: Wed, 6 Dec 2023 07:07:05 -0800
Subject: [PATCH 04/20] flake8

---
 des_y6utils/tests/test_asinh_mag.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/des_y6utils/tests/test_asinh_mag.py b/des_y6utils/tests/test_asinh_mag.py
index 9a60ea7..f4536b0 100644
--- a/des_y6utils/tests/test_asinh_mag.py
+++ b/des_y6utils/tests/test_asinh_mag.py
@@ -9,9 +9,10 @@ def test_compute_asinh_flux():
 
     assert flux_input == flux_g
 
+
 def test_compute_asinh_dered_mag():
     flux_input = 10000
     mag_g = _compute_asinh_mags(flux_input, 0)
-    dered_mag_g = _compute_asinh_dered_mag(flux_input, 0, 1/3.186)
+    dered_mag_g = _compute_asinh_dered_mag(flux_input, 0, 1 / 3.186)
 
-    assert mag_g == dered_mag_g
\ No newline at end of file
+    assert mag_g == dered_mag_g

From 241cf63198924df447f9fba54be9155b91ed7971 Mon Sep 17 00:00:00 2001
From: Masaya Yamamoto <masaya.yamamoto@duke.edu>
Date: Wed, 6 Dec 2023 07:37:03 -0800
Subject: [PATCH 05/20] fixed unit tests

---
 des_y6utils/tests/test_asinh_mag.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/des_y6utils/tests/test_asinh_mag.py b/des_y6utils/tests/test_asinh_mag.py
index f4536b0..e82c984 100644
--- a/des_y6utils/tests/test_asinh_mag.py
+++ b/des_y6utils/tests/test_asinh_mag.py
@@ -1,4 +1,5 @@
 from ..mdet import _compute_asinh_dered_mag, _compute_asinh_flux, _compute_asinh_mags
+import math
 
 
 def test_compute_asinh_flux():
@@ -7,7 +8,7 @@ def test_compute_asinh_flux():
     mag_g = _compute_asinh_mags(flux_input, 0)
     flux_g = _compute_asinh_flux(mag_g, 0)
 
-    assert flux_input == flux_g
+    assert math.isclose(flux_input, flux_g, rel_tol=1e-9)
 
 
 def test_compute_asinh_dered_mag():
@@ -15,4 +16,4 @@ def test_compute_asinh_dered_mag():
     mag_g = _compute_asinh_mags(flux_input, 0)
     dered_mag_g = _compute_asinh_dered_mag(flux_input, 0, 1 / 3.186)
 
-    assert mag_g == dered_mag_g
+    assert math.isclose(mag_g - 1.0, dered_mag_g, rel_tol=1e-9)

From 9a5a733c9faea8eae7a47de7791814f2346ef28e Mon Sep 17 00:00:00 2001
From: Masaya Yamamoto <masaya.yamamoto@duke.edu>
Date: Wed, 6 Dec 2023 14:13:31 -0800
Subject: [PATCH 06/20] rewrote how to apply dereddening and made changes to
 tests

---
 des_y6utils/mdet.py                 | 28 ++++++++++++++++++++++------
 des_y6utils/tests/test_asinh_mag.py | 12 ++++++++----
 2 files changed, 30 insertions(+), 10 deletions(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index 0110f4e..7f3108b 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -42,8 +42,7 @@ def add_extinction_correction_columns(fold, fnew, fdust):
         for ii, b in enumerate(bands):
             flux = np.copy(d["pgauss_band_flux_" + b])
             flux_og.append(flux)
-            mag_ = _compute_asinh_dered_mag(d["pgauss_band_flux_" + b], ii, dered)
-            flux_ = _compute_asinh_flux(mag_, ii)
+            flux_ = _compute_dered_flux(d["pgauss_band_flux_" + b], ii, dered)
             d["pgauss_band_flux_" + b] = flux_
 
         # make _nodered array with pgauss_band_flux_* entries, and add them to fits.
@@ -380,13 +379,30 @@ def _compute_asinh_mags(flux, i):
     return mag
 
 
-def _compute_asinh_dered_mag(flux, i, map):
+def _compute_dered_flux(flux, i, map):
+    """This function applies dereddening correction to fluxes.
+    Eli says we cannot apply the correction to the asinh mag.
+
+    Parameters
+    ----------
+    flux : float or np.ndarray
+        The flux.
+    i : int
+        The index of the band in griz (i.e., 0 for g, 1 for r, 2 for i, 3 for z).
+    map : float or np.ndarray
+        The E(B-V).
+
+    Returns
+    -------
+    dered_flux : float or np.ndarray
+        The dereddened flux.
+    """
 
     dered_fac = [3.186, 2.140, 1.196, 1.048]
-    mag = _compute_asinh_mags(flux, i)
-    dered_mag = mag - dered_fac[i] * map
+    dmag = dered_fac[i] * map
+    dered_flux = flux * 10**(dmag / 2.5)
 
-    return dered_mag
+    return dered_flux
 
 
 def _compute_asinh_flux(mag, i):
diff --git a/des_y6utils/tests/test_asinh_mag.py b/des_y6utils/tests/test_asinh_mag.py
index e82c984..747c2f1 100644
--- a/des_y6utils/tests/test_asinh_mag.py
+++ b/des_y6utils/tests/test_asinh_mag.py
@@ -1,4 +1,4 @@
-from ..mdet import _compute_asinh_dered_mag, _compute_asinh_flux, _compute_asinh_mags
+from ..mdet import _compute_dered_flux, _compute_asinh_flux, _compute_asinh_mags
 import math
 
 
@@ -11,9 +11,13 @@ def test_compute_asinh_flux():
     assert math.isclose(flux_input, flux_g, rel_tol=1e-9)
 
 
-def test_compute_asinh_dered_mag():
+def test_compute_dered_flux():
     flux_input = 10000
+    dered_flux = _compute_dered_flux(flux_input, 0, 2.5 / 3.186)
+
+    assert math.isclose(flux_input*10, dered_flux, rel_tol=1e-9)
+
     mag_g = _compute_asinh_mags(flux_input, 0)
-    dered_mag_g = _compute_asinh_dered_mag(flux_input, 0, 1 / 3.186)
+    dered_mag_g = _compute_asinh_mags(dered_flux, 0)
 
-    assert math.isclose(mag_g - 1.0, dered_mag_g, rel_tol=1e-9)
+    assert math.isclose(mag_g - 2.5, dered_mag_g, rel_tol=1e-8)

From 43bc46c7dbb4bada025910a0163c566baaee7f7b Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Tue, 12 Dec 2023 14:42:57 -0600
Subject: [PATCH 07/20] Update des_y6utils/mdet.py

---
 des_y6utils/mdet.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index 7f3108b..6c2be9d 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -749,7 +749,7 @@ def _make_mdet_cuts_raw_v6(
         & np.isfinite(mag_z)
         & (mag_g < 26.5)
         & (mag_r < 26.5)
-        & (mag_i < 24.5)
+        & (mag_i < 24.7)  # used to eliminate photo-z outliers motivated by cosmos2020
         & (mag_z < 25.6)
         & (d["pgauss_T"] < (1.6 - 3.1 * d["pgauss_T_err"]))
         & (

From 1c07c9dae15f1f4ce51e33d760dfbddf52215569 Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Fri, 15 Dec 2023 15:15:00 -0600
Subject: [PATCH 08/20] Apply suggestions from code review

---
 des_y6utils/mdet.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index 6c2be9d..b2de7aa 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -693,7 +693,7 @@ def _make_mdet_cuts_raw_v6(
 ):
     """The raw v6 cuts come from analysis in Fall of 2023.
 
-      - We implemented i-band magnitude cut at 24.5. TO-DO: still need to test.
+      - We implemented i-band magnitude cut at 24.7.
       - We implemented galaxy size cut at T=20.
       - We updated a cut in the pgauss T-Terr plane.
     """
@@ -750,6 +750,8 @@ def _make_mdet_cuts_raw_v6(
         & (mag_g < 26.5)
         & (mag_r < 26.5)
         & (mag_i < 24.7)  # used to eliminate photo-z outliers motivated by cosmos2020
+                          # the cut is 24.5 in cosmos mags with + 0.2 mags to correct
+                          # for the pgauss aperture
         & (mag_z < 25.6)
         & (d["pgauss_T"] < (1.6 - 3.1 * d["pgauss_T_err"]))
         & (

From 979c91230893f375e3d819c5eb091d0a1c2ae2b6 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Mon, 18 Dec 2023 13:03:46 -0600
Subject: [PATCH 09/20] REF refactor internal routines and unify some
 interfaces

---
 des_y6utils/__init__.py                       |   1 -
 des_y6utils/mdet.py                           | 128 ++--
 des_y6utils/piff/__init__.py                  |   9 -
 des_y6utils/piff/flag_good_regions.py         | 657 ------------------
 des_y6utils/piff/tests/__init__.py            |   0
 .../piff/tests/test_flag_good_regions.py      |  86 ---
 6 files changed, 58 insertions(+), 823 deletions(-)
 delete mode 100644 des_y6utils/piff/__init__.py
 delete mode 100644 des_y6utils/piff/flag_good_regions.py
 delete mode 100644 des_y6utils/piff/tests/__init__.py
 delete mode 100644 des_y6utils/piff/tests/test_flag_good_regions.py

diff --git a/des_y6utils/__init__.py b/des_y6utils/__init__.py
index e435846..d83bd64 100644
--- a/des_y6utils/__init__.py
+++ b/des_y6utils/__init__.py
@@ -12,4 +12,3 @@
 
 from . import mdet
 from . import viz
-from . import piff
diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index b2de7aa..af5acd0 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -7,60 +7,59 @@
 import healsparse
 
 
-def add_extinction_correction_columns(fold, fnew, fdust):
-    """This function changes and adds columns for extinction corrected magnitudes.
-    pgauss_band_flux_* are the columns that are the corrected fluxes.
-    pgauss_band_flux_*_nodered are the columns that used to be pgauss_band_flux_*.
-
-    Note: This function assumes you're working with fits files, not hdf5.
+def add_extinction_correction_columns(
+    data, dust_map_filename="SFD_dust_4096.hsp"
+):
+    """This function adds dereddened columns to the data and copies the existing columns
+    to columns with '_nodered' appended as a suffix.
 
     Parameters
     ----------
-    fold : original file path
-    fnew : new file path
+    data : np.ndarray
+        A structured array of data to be dereddened.
+    dust_map_filename : str
+        The path/name of the dust map. Assumed to be in HealSparse format.
 
     Returns
     -------
-    None (saves the new file)
+    new_data : np.ndarray
+        The new, dereddened data.
     """
 
-    import fitsio as fio
-    d = fio.read(fold)
     bands = ['g', 'r', 'i', 'z']
-    with fio.FITS(fnew, 'rw') as fits:
-
-        # compute dereddened fluxes and
-        # replace the entries of pgauss_flux_band_* with the dereddened fluxes.
-        if os.path.exists(fdust):
-            # fdust : '/global/cfs/cdirs/des/y6-shear-catalogs/SFD_dust_4096.hsp'
-            dustmap = _read_hsp_dustmap_local(fdust)
-        else:
-            dustmap = _read_hsp_mask(fdust)
-
-        dered = dustmap.get_values_pos(d["ra"], d["dec"])
-        flux_og = []
-        for ii, b in enumerate(bands):
-            flux = np.copy(d["pgauss_band_flux_" + b])
-            flux_og.append(flux)
-            flux_ = _compute_dered_flux(d["pgauss_band_flux_" + b], ii, dered)
-            d["pgauss_band_flux_" + b] = flux_
-
-        # make _nodered array with pgauss_band_flux_* entries, and add them to fits.
-        new_dt = np.dtype(d.dtype.descr + [('pgauss_band_flux_g_nodered', 'f8'),
-                                           ('pgauss_band_flux_r_nodered', 'f8'),
-                                           ('pgauss_band_flux_i_nodered', 'f8'),
-                                           ('pgauss_band_flux_z_nodered', 'f8'),])
-        d_ = np.zeros(d.shape, dtype=new_dt)
-        for col in d.dtype.names:
-            d_[col] = d[col]
-        for ii, b in enumerate(bands):
-            d_['pgauss_band_flux_' + b + '_nodered'] = flux_og[ii]
-
-        fits.write(d_)
-
-    fits.close()
-
-    return None
+
+    # compute dereddened fluxes and
+    # replace the entries of pgauss_flux_band_* with the dereddened fluxes.
+    if os.path.exists(dust_map_filename):
+        dustmap = _read_hsp_file_local(dust_map_filename)
+    else:
+        dustmap = _read_hsp_file(dust_map_filename)
+
+    dered = dustmap.get_values_pos(data["ra"], data["dec"])
+    flux_og = []
+    for ii, b in enumerate(bands):
+        flux = np.copy(data["pgauss_band_flux_" + b])
+        flux_og.append(flux)
+        flux_ = _compute_dered_flux(data["pgauss_band_flux_" + b], ii, dered)
+        data["pgauss_band_flux_" + b] = flux_
+
+    # make _nodered array with pgauss_band_flux_* entries, and add them to fits.
+    new_dt = np.dtype(
+        data.dtype.descr
+        + [
+            ('pgauss_band_flux_g_nodered', 'f8'),
+            ('pgauss_band_flux_r_nodered', 'f8'),
+            ('pgauss_band_flux_i_nodered', 'f8'),
+            ('pgauss_band_flux_z_nodered', 'f8'),
+        ]
+    )
+    new_data = np.zeros(data.shape, dtype=new_dt)
+    for col in data.dtype.names:
+        new_data[col] = data[col]
+    for ii, b in enumerate(bands):
+        new_data['pgauss_band_flux_' + b + '_nodered'] = flux_og[ii]
+
+    return new_data
 
 
 def make_mdet_cuts(data, version, verbose=False):
@@ -132,9 +131,9 @@ def _make_mdet_cuts_gauss(
 
     # apply the mask
     if os.path.exists(mask_name):
-        hmap = _read_hsp_mask_local(mask_name)
+        hmap = _read_hsp_file_local(mask_name)
     else:
-        hmap = _read_hsp_mask(mask_name)
+        hmap = _read_hsp_file(mask_name)
     in_footprint = hmap.get_values_pos(
         data["ra"], data["dec"], valid_mask=True
     )
@@ -173,9 +172,9 @@ def _make_mdet_cuts_wmom(
 
     # apply the mask
     if os.path.exists(mask_name):
-        hmap = _read_hsp_mask_local(mask_name)
+        hmap = _read_hsp_file_local(mask_name)
     else:
-        hmap = _read_hsp_mask(mask_name)
+        hmap = _read_hsp_file(mask_name)
     in_footprint = hmap.get_values_pos(
         data["ra"], data["dec"], valid_mask=True
     )
@@ -199,7 +198,7 @@ def _make_mdet_cuts_v1(d, verbose=False):
     )
 
     # apply the mask
-    hmap = _read_hsp_mask("y6-combined-hleda-gaiafull-hsmap16384-nomdet.fits")
+    hmap = _read_hsp_file("y6-combined-hleda-gaiafull-hsmap16384-nomdet.fits")
     in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
     msk &= in_footprint
     if verbose:
@@ -221,7 +220,7 @@ def _make_mdet_cuts_v2(d, verbose=False):
     )
 
     # apply the mask
-    hmap = _read_hsp_mask(
+    hmap = _read_hsp_file(
         "y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v2.fits"
     )
     in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
@@ -246,7 +245,7 @@ def _make_mdet_cuts_v3(d, verbose=False):
     )
 
     # apply the mask
-    hmap = _read_hsp_mask(
+    hmap = _read_hsp_file(
         "y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v2.fits"
     )
     in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
@@ -276,7 +275,7 @@ def _make_mdet_cuts_v4(d, verbose=False):
     msk &= ~msk_superspreader
 
     # apply the mask
-    hmap = _read_hsp_mask(
+    hmap = _read_hsp_file(
         "y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v2.fits"
     )
     in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
@@ -306,7 +305,7 @@ def _make_mdet_cuts_v5(d, verbose=False):
     msk &= ~msk_superspreader
 
     # apply the mask
-    hmap = _read_hsp_mask(
+    hmap = _read_hsp_file(
         "y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v3.fits"
     )
     in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
@@ -336,7 +335,7 @@ def _make_mdet_cuts_v6(d, verbose=False):
     msk &= ~msk_superspreader
 
     # apply the mask
-    hmap = _read_hsp_mask(
+    hmap = _read_hsp_file(
         "y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v1.hsp"
     )
     in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
@@ -431,28 +430,17 @@ def _compute_asinh_flux(mag, i):
 
 
 @lru_cache
-def _read_hsp_mask(fname):
-    mpth = _get_mask_path(fname)
-    return healsparse.HealSparseMap.read(mpth)
-
-
-@lru_cache
-def _read_hsp_mask_local(fname):
-    return healsparse.HealSparseMap.read(fname)
-
-
-@lru_cache
-def _read_hsp_dustmap(fname):
-    mpth = _get_mask_path(fname)
+def _read_hsp_file(fname):
+    mpth = _get_hsp_file_path(fname)
     return healsparse.HealSparseMap.read(mpth)
 
 
 @lru_cache
-def _read_hsp_dustmap_local(fname):
+def _read_hsp_file_local(fname):
     return healsparse.HealSparseMap.read(fname)
 
 
-def _get_mask_path(fname):
+def _get_hsp_file_path(fname):
     # get or make meds dir
     meds_dir = os.environ.get("MEDS_DIR", None)
     if meds_dir is None:
diff --git a/des_y6utils/piff/__init__.py b/des_y6utils/piff/__init__.py
deleted file mode 100644
index 412be8b..0000000
--- a/des_y6utils/piff/__init__.py
+++ /dev/null
@@ -1,9 +0,0 @@
-# flake8: noqa
-from .flag_good_regions import (
-    make_good_regions_for_piff_model_gal_grid,
-    make_good_regions_for_piff_model_star_and_gal_grid,
-    nanmad,
-    measure_t_grid_for_piff_model,
-    measure_star_t_for_piff_model,
-    map_star_t_to_grid,
-)
diff --git a/des_y6utils/piff/flag_good_regions.py b/des_y6utils/piff/flag_good_regions.py
deleted file mode 100644
index 5ed8647..0000000
--- a/des_y6utils/piff/flag_good_regions.py
+++ /dev/null
@@ -1,657 +0,0 @@
-import ngmix
-from ngmix.admom import AdmomFitter
-import galsim
-import copy
-import numpy as np
-from sklearn.preprocessing import PolynomialFeatures
-from sklearn.pipeline import make_pipeline
-from sklearn.linear_model import LinearRegression
-
-ALL_BAD = 2**0
-BAD_BOX = 2**1
-
-
-def _make_a_good_box_matts_hack(bad_msk, verbose=False):
-    flag = 0
-
-    b = dict(
-        xmax=bad_msk.shape[1]-1,
-        xmin=0,
-        ymax=bad_msk.shape[0]-1,
-        ymin=0,
-    )
-    if not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1]):
-        return b, flag
-
-    if np.all(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1]):
-        flag |= ALL_BAD
-        return b, flag
-
-    def _is_ok_range(b):
-        return (
-            (b["xmin"] <= b["xmax"])
-            and (b["ymin"] <= b["ymax"])
-        )
-
-    def _condition(b):
-        return (
-            (b["xmin"] <= b["xmax"])
-            and (b["ymin"] <= b["ymax"])
-            and np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
-        )
-
-    itr = 0
-    while _condition(b) and itr < 10_000:
-        new_b = copy.deepcopy(b)
-        curr_frac = np.sum(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
-        new_fracs = {}
-        for key, shift in zip(["ymax", "ymin", "xmax", "xmin"], [-1, 1, -1, 1]):
-            _new_b = copy.deepcopy(b)
-            _new_b[key] = b[key] + shift
-            if _is_ok_range(_new_b):
-                new_fracs[key] = np.sum(bad_msk[
-                    _new_b["ymin"]:_new_b["ymax"]+1,
-                    _new_b["xmin"]:_new_b["xmax"]+1]
-                )
-            else:
-                new_fracs[key] = np.nan
-
-        if np.all(np.isnan(list(new_fracs.values()))):
-            flag |= BAD_BOX
-            break
-
-        mval = np.nanmin(list(new_fracs.values()))
-        for key, shift in zip(["ymax", "ymin", "xmax", "xmin"], [-1, 1, -1, 1]):
-            _new_b = copy.deepcopy(b)
-            _new_b[key] = b[key] + shift
-            if new_fracs[key] == mval and _is_ok_range(_new_b):
-                new_b[key] = b[key] + shift
-                break
-
-        if verbose:
-            print("itr:", itr)
-            print("    curr # bad:", curr_frac)
-            print(
-                "    new # bad:",
-                " ".join("%s: %0.3f" % (k, v) for k, v in new_fracs.items()),
-            )
-            print("    new min val:", mval)
-            print(
-                "    edge|old|new: ",
-                " ".join("%s|%s|%s" % (k, b[k], new_b[k]) for k in b),
-            )
-        if new_b == b:
-            flag |= BAD_BOX
-            break
-
-        b = new_b
-        itr += 1
-
-    if (
-        b["xmin"] > b["xmax"]
-        or b["ymin"] > b["ymax"]
-        or np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
-    ):
-        flag |= BAD_BOX
-
-    return b, flag
-
-
-def _make_a_good_box(bad_msk, verbose=False):
-    """Maximum Empty Rectangle algorithm 1 from Naamad, Lee & Hsu, 1984
-​
-    https://www.sciencedirect.com/science/article/pii/0166218X84901240
-
-    This function was written by Mike Jarvis.
-    """
-
-    # Check for possible quick return
-    flag = 0
-    b = dict(
-        xmax=bad_msk.shape[1]-1,
-        xmin=0,
-        ymax=bad_msk.shape[0]-1,
-        ymin=0,
-    )
-    if not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1]):
-        return b, flag
-
-    if np.all(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1]):
-        flag |= ALL_BAD
-        return b, flag
-
-    # Algorithm MERAlg 1 (S,Al,Ar,Ab,At,MAXR)
-    # Input: Four boundary values of a rectangle A: Al, Ar, Ab, and At,
-    #        (left, right, bottom and top), and a set
-    #        S = {P1, P2, ... , Pn), Pi = (Xi, Yi) of points in A.
-    #
-    # Output: MAXR, the area of the MER defined by S and A.
-    #      ** Note: we change this to output the bounding box, not the area. **
-    #
-    # Method:
-    #
-    # 1. Let MGAP be the maximum gap in {Al, Ar, X1, X2, ... ,Xn}.
-    # 2. MAXR = MGAP * (At-Ab).
-    # 3. Sort S according to the Y coordinates of the points in descending order.
-    # 4. For i = 1 to n do steps 5-8.
-    # 5. Tl=Al, Tr=Ar.
-    # 6. For j=i+1 to n do step 7.
-    # 7. If Tl < Xj < Tr
-    #    Then do steps 7.1-7.2.
-    #       7.1. MAXR = MAX(MAXR, (Tr-Tl)*(Yi-Yj)).
-    #       7.2. If Xj>Xi
-    #            then Tr=Xj
-    #            else Tl=Xj
-    # 8. MAXR=MAX(MAXR, (Tr-Tl)*(Yi-Ab)).
-    # 9. For i = 1 to n do steps 10-12.
-    # 10. Ri = MIN(Ar U {Xj | (Xj,Yj) in S, Yj > Yi and Xj > Xi}).
-    # 11. Li = MAX(Al U {Xj | (Xj,Yj) in S, Yj > Yi and Xj < Xi}).
-    # 12. MAXR = MAX(MAXR, (Ri - Li) * (At - Yi)).
-
-    # Here is that algorithm translated to our case.
-    # Note: there are some modifications require to account for the fact that we have
-    # squares, not points as our bad area.
-    al = ab = -1
-    ar = bad_msk.shape[1]
-    at = bad_msk.shape[0]
-    y, x = np.where(bad_msk)
-
-    allx = np.sort(np.concatenate([[al, ar], x]))
-    gaps = np.diff(allx)
-    imgap = np.argmax(np.diff(allx))
-    maxr = (gaps[imgap]-1) * (at-ab-1)
-
-    # Keep track of the best answer yet.
-    # 5 numbers are maxr, tl, tr, tb, tt
-    # Note: the bounds we keep track of are the masked x,y just outside the box we want.
-    best = [maxr, allx[imgap], allx[imgap+1], ab, at]
-    if verbose:
-        print('initial best = ', best)
-
-    index = np.argsort(-y)
-    x = x[index]
-    y = y[index]
-
-    def update(tl, tr, tb, tt):
-        maxr = (tr-tl-1) * (tt-tb-1)
-        if maxr > best[0]:
-            best[:] = maxr, tl, tr, tb, tt
-            if verbose:
-                print('best => ', best)
-
-    n = len(x)
-    for i in range(n):
-        tl = al
-        tr = ar
-        for j in range(i+1, n):
-            if tl < x[j] < tr and y[j] < y[i]:
-                update(tl, tr, y[j], y[i])
-                if x[j] > x[i]:
-                    tr = x[j]
-                else:
-                    tl = x[j]
-        update(tl, tr, ab, y[i])
-    for i in range(n):
-        ri = np.min(x[(y > y[i]) & (x > x[i])], initial=ar)
-        li = np.max(x[(y > y[i]) & (x <= x[i])], initial=al)
-        update(li, ri, y[i], at)
-        ri = np.min(x[(y > y[i]) & (x >= x[i])], initial=ar)
-        li = np.max(x[(y > y[i]) & (x < x[i])], initial=al)
-        update(li, ri, y[i], at)
-
-    b = dict(
-        xmin=best[1]+1,
-        xmax=best[2]-1,
-        ymin=best[3]+1,
-        ymax=best[4]-1,
-    )
-    if best[0] == 0:
-        flag = BAD_BOX
-
-    if verbose:
-        print('final best = ', best)
-        print('b = ', b)
-
-    return b, flag
-
-
-def nanmad(x, axis=None):
-    """
-    median absolute deviation - scaled like a standard deviation
-
-        mad = 1.4826*median(|x-median(x)|)
-
-    Parameters
-    ----------
-    x: array-like
-        array to take MAD of
-    axis : {int, sequence of int, None}, optional
-        `axis` keyword for
-
-    Returns
-    -------
-    mad: float
-        MAD of array x
-    """
-    return 1.4826*np.nanmedian(np.abs(x - np.nanmedian(x, axis=axis)), axis=axis)
-
-
-def _rescale_b(b, grid_size):
-    b["xmin"] = b["xmin"] * grid_size
-    b["xmax"] = (b["xmax"] + 1) * grid_size
-    b["ymin"] = b["ymin"] * grid_size
-    b["ymax"] = (b["ymax"] + 1) * grid_size
-    return b
-
-
-def measure_t_grid_for_piff_model(
-    piff_mod, piff_kwargs, grid_size=128, seed=None
-):
-    """Make a bounding box of good regions for a given Piff model.
-
-    Parameters
-    ----------
-    piff_mod : piff Model
-        The Piff model read via `piff.read`.
-    piff_kwargs : dict, optional
-        Any extra keyword arguments to pass to `piff_mod.draw`. Typically these
-        might include things like `{"GI_COLOR": 0.61}`.
-    grid_size : int, optional
-        The grid size to construct the map. Must divide 2048 and 4096 evenly.
-        Default is 128.
-    seed : int, optional
-        The seed to use for the RNG fed to ngmix. Default of None will seed the
-        RNG on-the-fly.
-
-    Returns
-    -------
-    t_arr : array
-        The map of T over the model with one cell per `grid_size`.
-    """
-    piff_kwargs = piff_kwargs or {}
-    rng = np.random.RandomState(seed=seed)
-
-    if (
-        ((4096//grid_size) * grid_size != 4096)
-        or ((2048//grid_size) * grid_size != 2048)
-    ):
-        raise RuntimeError("The grid size must evenly divide 4096 and 2048.")
-
-    delta_to_pix = grid_size // 2
-    y = np.arange(4096//grid_size)*grid_size + delta_to_pix
-    x = np.arange(2048//grid_size)*grid_size + delta_to_pix
-
-    t_arr = np.zeros((y.shape[0], x.shape[0])) + np.nan
-    for i, yv in enumerate(y):
-        for j, xv in enumerate(x):
-            img = piff_mod.draw(
-                x=xv+1, y=yv+1, chipnum=list(piff_mod.wcs.keys())[0],
-                **piff_kwargs,
-            )
-            obs = ngmix.Observation(
-                image=img.array,
-                jacobian=ngmix.Jacobian(
-                    y=img.center.y - img.bounds.ymin,
-                    x=img.center.x - img.bounds.xmin,
-                    wcs=img.wcs.local(
-                        image_pos=galsim.PositionD(x=xv+1, y=yv+1)
-                    ).jacobian(),
-                )
-            )
-            # try 3 times just in case it fails at random
-            for _ in range(3):
-                try:
-                    res = AdmomFitter(
-                        rng=rng
-                    ).go(obs, ngmix.moments.fwhm_to_T(1))
-                    if res["flags"] == 0:
-                        t_arr[i, j] = res["T"]
-                except Exception:
-                    continue
-                else:
-                    break
-
-    return t_arr
-
-
-def _get_star_stamp_pos(s):
-    xint = int(np.floor(s.x - 1 + 0.5))
-    yint = int(np.floor(s.y - 1 + 0.5))
-    bbox = 17
-    bbox_2 = (bbox - 1)//2
-
-    return dict(
-        xstart=xint-bbox_2,
-        ystart=yint-bbox_2,
-        dim=bbox,
-        x=s.x - 1,
-        y=s.y - 1,
-    )
-
-
-def _get_star_piff_obs(piff_mod, s, piff_prop):
-
-    if piff_prop:
-        kwargs = {
-            piff_prop: s.data.properties[piff_prop]
-        }
-    else:
-        kwargs = {}
-    sres = _get_star_stamp_pos(s)
-
-    xv = sres["x"]+1
-    yv = sres["y"]+1
-    wcs = list(piff_mod.wcs.values())[0].local(
-        image_pos=galsim.PositionD(x=xv, y=yv)
-    ).jacobian()
-    img = galsim.ImageD(sres["dim"], sres["dim"], wcs=wcs)
-    cen = (
-        sres["x"] - sres["xstart"] + 1,
-        sres["y"] - sres["ystart"] + 1,
-    )
-    img = piff_mod.draw(
-        x=xv, y=yv, chipnum=list(piff_mod.wcs.keys())[0],
-        image=img, center=cen, **kwargs,
-    )
-    model_obs = ngmix.Observation(
-        image=img.array,
-        jacobian=ngmix.Jacobian(
-            y=cen[1]-1,
-            x=cen[0]-1,
-            wcs=wcs,
-        ),
-        meta=sres,
-    )
-    return model_obs
-
-
-def measure_star_t_for_piff_model(piff_mod, piff_prop=None, seed=None):
-    """Measure the model star sizes for a piff model.
-
-    Parameters
-    ----------
-    piff_mod : piff Model
-        The Piff model read via `piff.read`.
-    piff_prop : str, optional
-        The name of the piff property to use. Should be one of "GI_COLOR" or "IZ_COLOR"
-        for Y6.
-    seed : int, optional
-        The seed to use for the RNG fed to ngmix. Default of None will seed the
-        RNG on-the-fly.
-
-    Returns
-    -------
-    data : array
-        An array with columns x, y, t for the stars.
-    """
-    piff_prop = piff_prop or {}
-    rng = np.random.RandomState(seed=seed)
-
-    nd = len(piff_mod.stars)
-    data = np.zeros(nd, dtype=[("x", "f8"), ("y", "f8"), ("t", "f8")])
-    for col in data.dtype.names:
-        data[col] += np.nan
-
-    for i, s in enumerate(piff_mod.stars):
-
-        mobs = _get_star_piff_obs(piff_mod, s, piff_prop)
-
-        try:
-            res = ngmix.admom.AdmomFitter(
-                rng=rng
-            ).go(mobs, ngmix.moments.fwhm_to_T(1))
-            if res["flags"] == 0:
-                data["x"][i] = mobs.meta["x"]
-                data["y"][i] = mobs.meta["y"]
-                data["t"][i] = res["T"]
-        except Exception:
-            pass
-
-    return data
-
-
-def map_star_t_to_grid(star_data, grid_size=128, degree=2):
-    """Fit star T data to a polynomial and evaluate on a grid.
-
-    Parameters
-    ----------
-    star_data : array
-        The array pf star data from `measure_star_t_for_piff_model`.
-    grid_size : int, optional
-        The grid size to construct the map. Must divide 2048 and 4096 evenly.
-        Default is 128.
-    degree : int, optional
-        The degree of the polynomial.
-
-    Returns
-    -------
-    tg : array
-        The array of T values from evaluating the git on a grid.
-    """
-    # see this blog post
-    # https://towardsdatascience.com/polynomial-regression-with-scikit-learn-what-you-should-know-bed9d3296f2
-    polyreg = make_pipeline(PolynomialFeatures(degree), LinearRegression())
-    polyreg.fit(
-        np.array([star_data["x"], star_data["y"]]).T,
-        np.array(star_data["t"]),
-    )
-
-    if (
-        ((4096//grid_size) * grid_size != 4096)
-        or ((2048//grid_size) * grid_size != 2048)
-    ):
-        raise RuntimeError("The grid size must evenly divide 4096 and 2048.")
-
-    delta_to_pix = np.floor((grid_size-1)/2 + 0.5)
-    y, x = np.mgrid[0:4096:grid_size, 0:2048:grid_size] + delta_to_pix
-    tg = polyreg.predict(np.array([x.ravel(), y.ravel()]).T)
-    tg = tg.reshape(x.shape)
-
-    return tg
-
-
-def make_good_regions_for_piff_model_gal_grid(
-    piff_mod, piff_kwargs=None, grid_size=128, seed=None, verbose=False,
-    any_bad_thresh=25, flag_bad_thresh=15,
-):
-    """Make a bounding box of good regions for a given Piff model.
-
-    Parameters
-    ----------
-    piff_mod : piff Model
-        The Piff model read via `piff.read`.
-    piff_kwargs : dict, optional
-        Any extra keyword arguments to pass to `piff_mod.draw`. Typically these
-        might include things like `{"GI_COLOR": 0.61}`.
-    grid_size : int, optional
-        The grid size to construct the map. Must divide 2048 and 4096 evenly.
-        Default is 128.
-    seed : int, optional
-        The seed to use for the RNG fed to ngmix. Default of None will seed the
-        RNG on-the-fly.
-    verbose : bool, optional
-        If True, print some stats as the code is running. Default of False makes
-        the code silent.
-    any_bad_thresh : float, optional
-        The threshold used to figure out if any region of the CCD is bad. Models
-        with any regions where |t - t_mn| > any_bad_thresh * t_std are considered
-        bad. The default of 25 appears to account for PSF size variation while also
-        flagging weird models.
-    flag_bad_thresh : float, optional
-        The threshold used to mark good versus bad regions. Any region of the model
-        where |t - t_mn| > any_bad_thresh * t_std is marked as bad if any
-        region also exceeds the `any_bad_thresh`. Default of 15 accounts for
-        model size variation while also avoiding the very bad regions.
-
-    Returns
-    -------
-    data : dict
-        A dictionary with the following keys:
-
-            flags : int
-                The flags value. If non-zero, the entire model is flagged.
-            t_mn : float
-                The median T over the whole model.
-            t_std : float
-                The median absolute deviation over the whole model.
-            t_arr : array
-                The map of T over the model with one cell per `grid_size`.
-            bbox : dict
-                A dictionary with keys xmin, xmax, ymin, and ymax. Any model
-                whose center falls within this box is considered OK.
-    """
-    t_arr = measure_t_grid_for_piff_model(
-        piff_mod, piff_kwargs, grid_size=grid_size, seed=seed
-    )
-    flags = 0
-
-    b = dict(
-        xmax=t_arr.shape[1]-1,
-        xmin=0,
-        ymax=t_arr.shape[0]-1,
-        ymin=0,
-    )
-
-    if not np.any(np.isfinite(t_arr)):
-        flags |= ALL_BAD
-        t_mn = np.nan
-        t_std = np.nan
-
-    if flags == 0:
-        t_mn = np.nanmedian(t_arr)
-        t_std = nanmad(t_arr)
-
-        any_very_bad = (
-            (~np.isfinite(t_arr))
-            | (np.abs(t_arr - t_mn) > any_bad_thresh * t_std)
-        )
-        if np.any(any_very_bad):
-            some_bad = (
-                (~np.isfinite(t_arr))
-                | (np.abs(t_arr - t_mn) > flag_bad_thresh * t_std)
-            )
-            b, flag = _make_a_good_box(some_bad, verbose=verbose)
-            flags |= flag
-
-    return dict(
-        flags=flags,
-        t_mn=t_mn,
-        t_std=t_std,
-        t_arr=t_arr,
-        bbox=_rescale_b(b, grid_size),
-    )
-
-
-def make_good_regions_for_piff_model_star_and_gal_grid(
-    piff_mod, piff_kwargs=None, grid_size=128, seed=None, verbose=False,
-    any_bad_thresh=5, flag_bad_thresh=5, degree=2,
-):
-    """Make a bounding box of good regions for a given Piff model.
-
-    Parameters
-    ----------
-    piff_mod : piff Model
-        The Piff model read via `piff.read`.
-    piff_kwargs : dict, optional
-        Any extra keyword arguments to pass to `piff_mod.draw`. Typically these
-        might include things like `{"GI_COLOR": 0.61}`.
-    grid_size : int, optional
-        The grid size to construct the map. Must divide 2048 and 4096 evenly.
-        Default is 128.
-    seed : int, optional
-        The seed to use for the RNG fed to ngmix. Default of None will seed the
-        RNG on-the-fly.
-    verbose : bool, optional
-        If True, print some stats as the code is running. Default of False makes
-        the code silent.
-    any_bad_thresh : float, optional
-        The threshold used to figure out if any region of the CCD is bad. Models
-        with any regions where |t - t_mn| > any_bad_thresh * t_std are considered
-        bad. The default of 25 appears to account for PSF size variation while also
-        flagging weird models.
-    flag_bad_thresh : float, optional
-        The threshold used to mark good versus bad regions. Any region of the model
-        where |t - t_mn| > any_bad_thresh * t_std is marked as bad if any
-        region also exceeds the `any_bad_thresh`. Default of 15 accounts for
-        model size variation while also avoiding the very bad regions.
-    degree : int, optional
-        The degree of the polynomial.
-
-    Returns
-    -------
-    data : dict
-        A dictionary with the following keys:
-
-            flags : int
-                The flags value. If non-zero, the entire model is flagged.
-            t_star : float
-                The star T map.
-            t_gal : float
-                The gal T map.
-            bbox : dict
-                A dictionary with keys xmin, xmax, ymin, and ymax. Any model
-                whose center falls within this box is considered OK.
-    """
-    flags = 0
-
-    tg_arr = measure_t_grid_for_piff_model(
-        piff_mod, piff_kwargs, grid_size=grid_size, seed=seed
-    )
-    data = measure_star_t_for_piff_model(
-        piff_mod, piff_prop=list(piff_kwargs.keys())[0], seed=seed,
-    )
-
-    if np.all(np.isnan(data["t"])):
-        flags |= ALL_BAD
-        bad_msk = np.ones_like(tg_arr).astype(bool)
-
-    if flags == 0:
-        msk = np.isfinite(data["t"])
-
-        ts_arr = map_star_t_to_grid(
-            data[msk], grid_size=grid_size, degree=degree,
-        )
-    else:
-        ts_arr = None
-        bad_msk = np.ones_like(tg_arr).astype(bool)
-
-    b = dict(
-        xmax=tg_arr.shape[1]-1,
-        xmin=0,
-        ymax=tg_arr.shape[0]-1,
-        ymin=0,
-    )
-
-    if flags == 0 and (not np.any(np.isfinite(tg_arr))) or np.any(np.isnan(ts_arr)):
-        flags |= ALL_BAD
-        bad_msk = np.ones_like(tg_arr).astype(bool)
-
-    if flags == 0:
-        tdiff = tg_arr - ts_arr
-        t_mn = np.nanmedian(tdiff)
-        t_std = nanmad(tdiff)
-
-        any_very_bad = (
-            (~np.isfinite(tdiff))
-            | (np.abs(tdiff - t_mn) > any_bad_thresh * t_std)
-        )
-        if np.any(any_very_bad):
-            some_bad = (
-                (~np.isfinite(tg_arr))
-                | (np.abs(tdiff - t_mn) > flag_bad_thresh * t_std)
-            )
-            bad_msk = some_bad
-            b, flag = _make_a_good_box(some_bad, verbose=verbose)
-            flags |= flag
-        else:
-            bad_msk = np.zeros_like(tg_arr).astype(bool)
-    else:
-        bad_msk = np.ones_like(tg_arr).astype(bool)
-
-    return dict(
-        flags=flags,
-        t_star=ts_arr,
-        t_gal=tg_arr,
-        bad_msk=bad_msk,
-        bbox=_rescale_b(b, grid_size),
-    )
diff --git a/des_y6utils/piff/tests/__init__.py b/des_y6utils/piff/tests/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/des_y6utils/piff/tests/test_flag_good_regions.py b/des_y6utils/piff/tests/test_flag_good_regions.py
deleted file mode 100644
index 70849d8..0000000
--- a/des_y6utils/piff/tests/test_flag_good_regions.py
+++ /dev/null
@@ -1,86 +0,0 @@
-import numpy as np
-
-from ..flag_good_regions import _make_a_good_box
-
-
-def test_flag_good_regions_make_a_good_box_smoke():
-    bad_msk = np.zeros((12, 15)).astype(bool)
-    b, flag = _make_a_good_box(bad_msk)
-
-    assert flag == 0
-    assert b == {
-        "xmin": 0,
-        "xmax": 14,
-        "ymin": 0,
-        "ymax": 11,
-    }
-    assert not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
-
-
-def test_flag_good_regions_make_a_good_box():
-    bad_msk = np.zeros((12, 15)).astype(bool)
-    bad_msk[:, 0] = True
-    b, flag = _make_a_good_box(bad_msk)
-
-    assert flag == 0
-    assert b == {
-        "xmin": 1,
-        "xmax": 14,
-        "ymin": 0,
-        "ymax": 11,
-    }
-    assert not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
-
-
-def test_flag_good_regions_make_a_good_box_cen():
-    bad_msk = np.zeros((13, 15)).astype(bool)
-    bad_msk[6, 7] = True
-    b, flag = _make_a_good_box(bad_msk)
-    assert flag == 0
-    assert b == {
-        "xmin": 0,
-        "xmax": 6,
-        "ymin": 0,
-        "ymax": 12,
-    }
-    assert not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
-
-
-def test_flag_good_regions_make_a_good_box_cen_inv():
-    bad_msk = np.ones((13, 15)).astype(bool)
-    bad_msk[6, 7] = False
-    b, flag = _make_a_good_box(bad_msk)
-    assert flag == 0
-    assert b == {
-        "xmin": 7,
-        "xmax": 7,
-        "ymin": 6,
-        "ymax": 6,
-    }
-    assert not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
-
-
-def test_flag_good_regions_make_a_good_box_random():
-    rng = np.random.RandomState(seed=10)
-    some_ok = False
-    for _ in range(100):
-        bad_msk = rng.uniform(size=(12, 15)) > 0.8
-        b, flag = _make_a_good_box(bad_msk)
-        if flag == 0:
-            assert not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
-            some_ok = True
-
-    assert some_ok
-
-
-def test_flag_good_regions_make_a_good_box_all_bad():
-    bad_msk = np.ones((12, 15)).astype(bool)
-    b, flag = _make_a_good_box(bad_msk)
-
-    assert flag != 0
-    assert b == {
-        "xmin": 0,
-        "xmax": 14,
-        "ymin": 0,
-        "ymax": 11,
-    }

From 4376c56dbe4c58db272f6040dbde6b6735860abf Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Mon, 18 Dec 2023 13:06:50 -0600
Subject: [PATCH 10/20] BUG put back piff items

---
 des_y6utils/piff/__init__.py                  |   9 +
 des_y6utils/piff/flag_good_regions.py         | 657 ++++++++++++++++++
 des_y6utils/piff/tests/__init__.py            |   0
 .../piff/tests/test_flag_good_regions.py      |  86 +++
 4 files changed, 752 insertions(+)
 create mode 100644 des_y6utils/piff/__init__.py
 create mode 100644 des_y6utils/piff/flag_good_regions.py
 create mode 100644 des_y6utils/piff/tests/__init__.py
 create mode 100644 des_y6utils/piff/tests/test_flag_good_regions.py

diff --git a/des_y6utils/piff/__init__.py b/des_y6utils/piff/__init__.py
new file mode 100644
index 0000000..412be8b
--- /dev/null
+++ b/des_y6utils/piff/__init__.py
@@ -0,0 +1,9 @@
+# flake8: noqa
+from .flag_good_regions import (
+    make_good_regions_for_piff_model_gal_grid,
+    make_good_regions_for_piff_model_star_and_gal_grid,
+    nanmad,
+    measure_t_grid_for_piff_model,
+    measure_star_t_for_piff_model,
+    map_star_t_to_grid,
+)
diff --git a/des_y6utils/piff/flag_good_regions.py b/des_y6utils/piff/flag_good_regions.py
new file mode 100644
index 0000000..5ed8647
--- /dev/null
+++ b/des_y6utils/piff/flag_good_regions.py
@@ -0,0 +1,657 @@
+import ngmix
+from ngmix.admom import AdmomFitter
+import galsim
+import copy
+import numpy as np
+from sklearn.preprocessing import PolynomialFeatures
+from sklearn.pipeline import make_pipeline
+from sklearn.linear_model import LinearRegression
+
+ALL_BAD = 2**0
+BAD_BOX = 2**1
+
+
+def _make_a_good_box_matts_hack(bad_msk, verbose=False):
+    flag = 0
+
+    b = dict(
+        xmax=bad_msk.shape[1]-1,
+        xmin=0,
+        ymax=bad_msk.shape[0]-1,
+        ymin=0,
+    )
+    if not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1]):
+        return b, flag
+
+    if np.all(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1]):
+        flag |= ALL_BAD
+        return b, flag
+
+    def _is_ok_range(b):
+        return (
+            (b["xmin"] <= b["xmax"])
+            and (b["ymin"] <= b["ymax"])
+        )
+
+    def _condition(b):
+        return (
+            (b["xmin"] <= b["xmax"])
+            and (b["ymin"] <= b["ymax"])
+            and np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
+        )
+
+    itr = 0
+    while _condition(b) and itr < 10_000:
+        new_b = copy.deepcopy(b)
+        curr_frac = np.sum(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
+        new_fracs = {}
+        for key, shift in zip(["ymax", "ymin", "xmax", "xmin"], [-1, 1, -1, 1]):
+            _new_b = copy.deepcopy(b)
+            _new_b[key] = b[key] + shift
+            if _is_ok_range(_new_b):
+                new_fracs[key] = np.sum(bad_msk[
+                    _new_b["ymin"]:_new_b["ymax"]+1,
+                    _new_b["xmin"]:_new_b["xmax"]+1]
+                )
+            else:
+                new_fracs[key] = np.nan
+
+        if np.all(np.isnan(list(new_fracs.values()))):
+            flag |= BAD_BOX
+            break
+
+        mval = np.nanmin(list(new_fracs.values()))
+        for key, shift in zip(["ymax", "ymin", "xmax", "xmin"], [-1, 1, -1, 1]):
+            _new_b = copy.deepcopy(b)
+            _new_b[key] = b[key] + shift
+            if new_fracs[key] == mval and _is_ok_range(_new_b):
+                new_b[key] = b[key] + shift
+                break
+
+        if verbose:
+            print("itr:", itr)
+            print("    curr # bad:", curr_frac)
+            print(
+                "    new # bad:",
+                " ".join("%s: %0.3f" % (k, v) for k, v in new_fracs.items()),
+            )
+            print("    new min val:", mval)
+            print(
+                "    edge|old|new: ",
+                " ".join("%s|%s|%s" % (k, b[k], new_b[k]) for k in b),
+            )
+        if new_b == b:
+            flag |= BAD_BOX
+            break
+
+        b = new_b
+        itr += 1
+
+    if (
+        b["xmin"] > b["xmax"]
+        or b["ymin"] > b["ymax"]
+        or np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
+    ):
+        flag |= BAD_BOX
+
+    return b, flag
+
+
+def _make_a_good_box(bad_msk, verbose=False):
+    """Maximum Empty Rectangle algorithm 1 from Naamad, Lee & Hsu, 1984
+​
+    https://www.sciencedirect.com/science/article/pii/0166218X84901240
+
+    This function was written by Mike Jarvis.
+    """
+
+    # Check for possible quick return
+    flag = 0
+    b = dict(
+        xmax=bad_msk.shape[1]-1,
+        xmin=0,
+        ymax=bad_msk.shape[0]-1,
+        ymin=0,
+    )
+    if not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1]):
+        return b, flag
+
+    if np.all(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1]):
+        flag |= ALL_BAD
+        return b, flag
+
+    # Algorithm MERAlg 1 (S,Al,Ar,Ab,At,MAXR)
+    # Input: Four boundary values of a rectangle A: Al, Ar, Ab, and At,
+    #        (left, right, bottom and top), and a set
+    #        S = {P1, P2, ... , Pn), Pi = (Xi, Yi) of points in A.
+    #
+    # Output: MAXR, the area of the MER defined by S and A.
+    #      ** Note: we change this to output the bounding box, not the area. **
+    #
+    # Method:
+    #
+    # 1. Let MGAP be the maximum gap in {Al, Ar, X1, X2, ... ,Xn}.
+    # 2. MAXR = MGAP * (At-Ab).
+    # 3. Sort S according to the Y coordinates of the points in descending order.
+    # 4. For i = 1 to n do steps 5-8.
+    # 5. Tl=Al, Tr=Ar.
+    # 6. For j=i+1 to n do step 7.
+    # 7. If Tl < Xj < Tr
+    #    Then do steps 7.1-7.2.
+    #       7.1. MAXR = MAX(MAXR, (Tr-Tl)*(Yi-Yj)).
+    #       7.2. If Xj>Xi
+    #            then Tr=Xj
+    #            else Tl=Xj
+    # 8. MAXR=MAX(MAXR, (Tr-Tl)*(Yi-Ab)).
+    # 9. For i = 1 to n do steps 10-12.
+    # 10. Ri = MIN(Ar U {Xj | (Xj,Yj) in S, Yj > Yi and Xj > Xi}).
+    # 11. Li = MAX(Al U {Xj | (Xj,Yj) in S, Yj > Yi and Xj < Xi}).
+    # 12. MAXR = MAX(MAXR, (Ri - Li) * (At - Yi)).
+
+    # Here is that algorithm translated to our case.
+    # Note: there are some modifications require to account for the fact that we have
+    # squares, not points as our bad area.
+    al = ab = -1
+    ar = bad_msk.shape[1]
+    at = bad_msk.shape[0]
+    y, x = np.where(bad_msk)
+
+    allx = np.sort(np.concatenate([[al, ar], x]))
+    gaps = np.diff(allx)
+    imgap = np.argmax(np.diff(allx))
+    maxr = (gaps[imgap]-1) * (at-ab-1)
+
+    # Keep track of the best answer yet.
+    # 5 numbers are maxr, tl, tr, tb, tt
+    # Note: the bounds we keep track of are the masked x,y just outside the box we want.
+    best = [maxr, allx[imgap], allx[imgap+1], ab, at]
+    if verbose:
+        print('initial best = ', best)
+
+    index = np.argsort(-y)
+    x = x[index]
+    y = y[index]
+
+    def update(tl, tr, tb, tt):
+        maxr = (tr-tl-1) * (tt-tb-1)
+        if maxr > best[0]:
+            best[:] = maxr, tl, tr, tb, tt
+            if verbose:
+                print('best => ', best)
+
+    n = len(x)
+    for i in range(n):
+        tl = al
+        tr = ar
+        for j in range(i+1, n):
+            if tl < x[j] < tr and y[j] < y[i]:
+                update(tl, tr, y[j], y[i])
+                if x[j] > x[i]:
+                    tr = x[j]
+                else:
+                    tl = x[j]
+        update(tl, tr, ab, y[i])
+    for i in range(n):
+        ri = np.min(x[(y > y[i]) & (x > x[i])], initial=ar)
+        li = np.max(x[(y > y[i]) & (x <= x[i])], initial=al)
+        update(li, ri, y[i], at)
+        ri = np.min(x[(y > y[i]) & (x >= x[i])], initial=ar)
+        li = np.max(x[(y > y[i]) & (x < x[i])], initial=al)
+        update(li, ri, y[i], at)
+
+    b = dict(
+        xmin=best[1]+1,
+        xmax=best[2]-1,
+        ymin=best[3]+1,
+        ymax=best[4]-1,
+    )
+    if best[0] == 0:
+        flag = BAD_BOX
+
+    if verbose:
+        print('final best = ', best)
+        print('b = ', b)
+
+    return b, flag
+
+
+def nanmad(x, axis=None):
+    """
+    median absolute deviation - scaled like a standard deviation
+
+        mad = 1.4826*median(|x-median(x)|)
+
+    Parameters
+    ----------
+    x: array-like
+        array to take MAD of
+    axis : {int, sequence of int, None}, optional
+        `axis` keyword for
+
+    Returns
+    -------
+    mad: float
+        MAD of array x
+    """
+    return 1.4826*np.nanmedian(np.abs(x - np.nanmedian(x, axis=axis)), axis=axis)
+
+
+def _rescale_b(b, grid_size):
+    b["xmin"] = b["xmin"] * grid_size
+    b["xmax"] = (b["xmax"] + 1) * grid_size
+    b["ymin"] = b["ymin"] * grid_size
+    b["ymax"] = (b["ymax"] + 1) * grid_size
+    return b
+
+
+def measure_t_grid_for_piff_model(
+    piff_mod, piff_kwargs, grid_size=128, seed=None
+):
+    """Make a bounding box of good regions for a given Piff model.
+
+    Parameters
+    ----------
+    piff_mod : piff Model
+        The Piff model read via `piff.read`.
+    piff_kwargs : dict, optional
+        Any extra keyword arguments to pass to `piff_mod.draw`. Typically these
+        might include things like `{"GI_COLOR": 0.61}`.
+    grid_size : int, optional
+        The grid size to construct the map. Must divide 2048 and 4096 evenly.
+        Default is 128.
+    seed : int, optional
+        The seed to use for the RNG fed to ngmix. Default of None will seed the
+        RNG on-the-fly.
+
+    Returns
+    -------
+    t_arr : array
+        The map of T over the model with one cell per `grid_size`.
+    """
+    piff_kwargs = piff_kwargs or {}
+    rng = np.random.RandomState(seed=seed)
+
+    if (
+        ((4096//grid_size) * grid_size != 4096)
+        or ((2048//grid_size) * grid_size != 2048)
+    ):
+        raise RuntimeError("The grid size must evenly divide 4096 and 2048.")
+
+    delta_to_pix = grid_size // 2
+    y = np.arange(4096//grid_size)*grid_size + delta_to_pix
+    x = np.arange(2048//grid_size)*grid_size + delta_to_pix
+
+    t_arr = np.zeros((y.shape[0], x.shape[0])) + np.nan
+    for i, yv in enumerate(y):
+        for j, xv in enumerate(x):
+            img = piff_mod.draw(
+                x=xv+1, y=yv+1, chipnum=list(piff_mod.wcs.keys())[0],
+                **piff_kwargs,
+            )
+            obs = ngmix.Observation(
+                image=img.array,
+                jacobian=ngmix.Jacobian(
+                    y=img.center.y - img.bounds.ymin,
+                    x=img.center.x - img.bounds.xmin,
+                    wcs=img.wcs.local(
+                        image_pos=galsim.PositionD(x=xv+1, y=yv+1)
+                    ).jacobian(),
+                )
+            )
+            # try 3 times just in case it fails at random
+            for _ in range(3):
+                try:
+                    res = AdmomFitter(
+                        rng=rng
+                    ).go(obs, ngmix.moments.fwhm_to_T(1))
+                    if res["flags"] == 0:
+                        t_arr[i, j] = res["T"]
+                except Exception:
+                    continue
+                else:
+                    break
+
+    return t_arr
+
+
+def _get_star_stamp_pos(s):
+    xint = int(np.floor(s.x - 1 + 0.5))
+    yint = int(np.floor(s.y - 1 + 0.5))
+    bbox = 17
+    bbox_2 = (bbox - 1)//2
+
+    return dict(
+        xstart=xint-bbox_2,
+        ystart=yint-bbox_2,
+        dim=bbox,
+        x=s.x - 1,
+        y=s.y - 1,
+    )
+
+
+def _get_star_piff_obs(piff_mod, s, piff_prop):
+
+    if piff_prop:
+        kwargs = {
+            piff_prop: s.data.properties[piff_prop]
+        }
+    else:
+        kwargs = {}
+    sres = _get_star_stamp_pos(s)
+
+    xv = sres["x"]+1
+    yv = sres["y"]+1
+    wcs = list(piff_mod.wcs.values())[0].local(
+        image_pos=galsim.PositionD(x=xv, y=yv)
+    ).jacobian()
+    img = galsim.ImageD(sres["dim"], sres["dim"], wcs=wcs)
+    cen = (
+        sres["x"] - sres["xstart"] + 1,
+        sres["y"] - sres["ystart"] + 1,
+    )
+    img = piff_mod.draw(
+        x=xv, y=yv, chipnum=list(piff_mod.wcs.keys())[0],
+        image=img, center=cen, **kwargs,
+    )
+    model_obs = ngmix.Observation(
+        image=img.array,
+        jacobian=ngmix.Jacobian(
+            y=cen[1]-1,
+            x=cen[0]-1,
+            wcs=wcs,
+        ),
+        meta=sres,
+    )
+    return model_obs
+
+
+def measure_star_t_for_piff_model(piff_mod, piff_prop=None, seed=None):
+    """Measure the model star sizes for a piff model.
+
+    Parameters
+    ----------
+    piff_mod : piff Model
+        The Piff model read via `piff.read`.
+    piff_prop : str, optional
+        The name of the piff property to use. Should be one of "GI_COLOR" or "IZ_COLOR"
+        for Y6.
+    seed : int, optional
+        The seed to use for the RNG fed to ngmix. Default of None will seed the
+        RNG on-the-fly.
+
+    Returns
+    -------
+    data : array
+        An array with columns x, y, t for the stars.
+    """
+    piff_prop = piff_prop or {}
+    rng = np.random.RandomState(seed=seed)
+
+    nd = len(piff_mod.stars)
+    data = np.zeros(nd, dtype=[("x", "f8"), ("y", "f8"), ("t", "f8")])
+    for col in data.dtype.names:
+        data[col] += np.nan
+
+    for i, s in enumerate(piff_mod.stars):
+
+        mobs = _get_star_piff_obs(piff_mod, s, piff_prop)
+
+        try:
+            res = ngmix.admom.AdmomFitter(
+                rng=rng
+            ).go(mobs, ngmix.moments.fwhm_to_T(1))
+            if res["flags"] == 0:
+                data["x"][i] = mobs.meta["x"]
+                data["y"][i] = mobs.meta["y"]
+                data["t"][i] = res["T"]
+        except Exception:
+            pass
+
+    return data
+
+
+def map_star_t_to_grid(star_data, grid_size=128, degree=2):
+    """Fit star T data to a polynomial and evaluate on a grid.
+
+    Parameters
+    ----------
+    star_data : array
+        The array pf star data from `measure_star_t_for_piff_model`.
+    grid_size : int, optional
+        The grid size to construct the map. Must divide 2048 and 4096 evenly.
+        Default is 128.
+    degree : int, optional
+        The degree of the polynomial.
+
+    Returns
+    -------
+    tg : array
+        The array of T values from evaluating the git on a grid.
+    """
+    # see this blog post
+    # https://towardsdatascience.com/polynomial-regression-with-scikit-learn-what-you-should-know-bed9d3296f2
+    polyreg = make_pipeline(PolynomialFeatures(degree), LinearRegression())
+    polyreg.fit(
+        np.array([star_data["x"], star_data["y"]]).T,
+        np.array(star_data["t"]),
+    )
+
+    if (
+        ((4096//grid_size) * grid_size != 4096)
+        or ((2048//grid_size) * grid_size != 2048)
+    ):
+        raise RuntimeError("The grid size must evenly divide 4096 and 2048.")
+
+    delta_to_pix = np.floor((grid_size-1)/2 + 0.5)
+    y, x = np.mgrid[0:4096:grid_size, 0:2048:grid_size] + delta_to_pix
+    tg = polyreg.predict(np.array([x.ravel(), y.ravel()]).T)
+    tg = tg.reshape(x.shape)
+
+    return tg
+
+
+def make_good_regions_for_piff_model_gal_grid(
+    piff_mod, piff_kwargs=None, grid_size=128, seed=None, verbose=False,
+    any_bad_thresh=25, flag_bad_thresh=15,
+):
+    """Make a bounding box of good regions for a given Piff model.
+
+    Parameters
+    ----------
+    piff_mod : piff Model
+        The Piff model read via `piff.read`.
+    piff_kwargs : dict, optional
+        Any extra keyword arguments to pass to `piff_mod.draw`. Typically these
+        might include things like `{"GI_COLOR": 0.61}`.
+    grid_size : int, optional
+        The grid size to construct the map. Must divide 2048 and 4096 evenly.
+        Default is 128.
+    seed : int, optional
+        The seed to use for the RNG fed to ngmix. Default of None will seed the
+        RNG on-the-fly.
+    verbose : bool, optional
+        If True, print some stats as the code is running. Default of False makes
+        the code silent.
+    any_bad_thresh : float, optional
+        The threshold used to figure out if any region of the CCD is bad. Models
+        with any regions where |t - t_mn| > any_bad_thresh * t_std are considered
+        bad. The default of 25 appears to account for PSF size variation while also
+        flagging weird models.
+    flag_bad_thresh : float, optional
+        The threshold used to mark good versus bad regions. Any region of the model
+        where |t - t_mn| > any_bad_thresh * t_std is marked as bad if any
+        region also exceeds the `any_bad_thresh`. Default of 15 accounts for
+        model size variation while also avoiding the very bad regions.
+
+    Returns
+    -------
+    data : dict
+        A dictionary with the following keys:
+
+            flags : int
+                The flags value. If non-zero, the entire model is flagged.
+            t_mn : float
+                The median T over the whole model.
+            t_std : float
+                The median absolute deviation over the whole model.
+            t_arr : array
+                The map of T over the model with one cell per `grid_size`.
+            bbox : dict
+                A dictionary with keys xmin, xmax, ymin, and ymax. Any model
+                whose center falls within this box is considered OK.
+    """
+    t_arr = measure_t_grid_for_piff_model(
+        piff_mod, piff_kwargs, grid_size=grid_size, seed=seed
+    )
+    flags = 0
+
+    b = dict(
+        xmax=t_arr.shape[1]-1,
+        xmin=0,
+        ymax=t_arr.shape[0]-1,
+        ymin=0,
+    )
+
+    if not np.any(np.isfinite(t_arr)):
+        flags |= ALL_BAD
+        t_mn = np.nan
+        t_std = np.nan
+
+    if flags == 0:
+        t_mn = np.nanmedian(t_arr)
+        t_std = nanmad(t_arr)
+
+        any_very_bad = (
+            (~np.isfinite(t_arr))
+            | (np.abs(t_arr - t_mn) > any_bad_thresh * t_std)
+        )
+        if np.any(any_very_bad):
+            some_bad = (
+                (~np.isfinite(t_arr))
+                | (np.abs(t_arr - t_mn) > flag_bad_thresh * t_std)
+            )
+            b, flag = _make_a_good_box(some_bad, verbose=verbose)
+            flags |= flag
+
+    return dict(
+        flags=flags,
+        t_mn=t_mn,
+        t_std=t_std,
+        t_arr=t_arr,
+        bbox=_rescale_b(b, grid_size),
+    )
+
+
+def make_good_regions_for_piff_model_star_and_gal_grid(
+    piff_mod, piff_kwargs=None, grid_size=128, seed=None, verbose=False,
+    any_bad_thresh=5, flag_bad_thresh=5, degree=2,
+):
+    """Make a bounding box of good regions for a given Piff model.
+
+    Parameters
+    ----------
+    piff_mod : piff Model
+        The Piff model read via `piff.read`.
+    piff_kwargs : dict, optional
+        Any extra keyword arguments to pass to `piff_mod.draw`. Typically these
+        might include things like `{"GI_COLOR": 0.61}`.
+    grid_size : int, optional
+        The grid size to construct the map. Must divide 2048 and 4096 evenly.
+        Default is 128.
+    seed : int, optional
+        The seed to use for the RNG fed to ngmix. Default of None will seed the
+        RNG on-the-fly.
+    verbose : bool, optional
+        If True, print some stats as the code is running. Default of False makes
+        the code silent.
+    any_bad_thresh : float, optional
+        The threshold used to figure out if any region of the CCD is bad. Models
+        with any regions where |t - t_mn| > any_bad_thresh * t_std are considered
+        bad. The default of 25 appears to account for PSF size variation while also
+        flagging weird models.
+    flag_bad_thresh : float, optional
+        The threshold used to mark good versus bad regions. Any region of the model
+        where |t - t_mn| > any_bad_thresh * t_std is marked as bad if any
+        region also exceeds the `any_bad_thresh`. Default of 15 accounts for
+        model size variation while also avoiding the very bad regions.
+    degree : int, optional
+        The degree of the polynomial.
+
+    Returns
+    -------
+    data : dict
+        A dictionary with the following keys:
+
+            flags : int
+                The flags value. If non-zero, the entire model is flagged.
+            t_star : float
+                The star T map.
+            t_gal : float
+                The gal T map.
+            bbox : dict
+                A dictionary with keys xmin, xmax, ymin, and ymax. Any model
+                whose center falls within this box is considered OK.
+    """
+    flags = 0
+
+    tg_arr = measure_t_grid_for_piff_model(
+        piff_mod, piff_kwargs, grid_size=grid_size, seed=seed
+    )
+    data = measure_star_t_for_piff_model(
+        piff_mod, piff_prop=list(piff_kwargs.keys())[0], seed=seed,
+    )
+
+    if np.all(np.isnan(data["t"])):
+        flags |= ALL_BAD
+        bad_msk = np.ones_like(tg_arr).astype(bool)
+
+    if flags == 0:
+        msk = np.isfinite(data["t"])
+
+        ts_arr = map_star_t_to_grid(
+            data[msk], grid_size=grid_size, degree=degree,
+        )
+    else:
+        ts_arr = None
+        bad_msk = np.ones_like(tg_arr).astype(bool)
+
+    b = dict(
+        xmax=tg_arr.shape[1]-1,
+        xmin=0,
+        ymax=tg_arr.shape[0]-1,
+        ymin=0,
+    )
+
+    if flags == 0 and (not np.any(np.isfinite(tg_arr))) or np.any(np.isnan(ts_arr)):
+        flags |= ALL_BAD
+        bad_msk = np.ones_like(tg_arr).astype(bool)
+
+    if flags == 0:
+        tdiff = tg_arr - ts_arr
+        t_mn = np.nanmedian(tdiff)
+        t_std = nanmad(tdiff)
+
+        any_very_bad = (
+            (~np.isfinite(tdiff))
+            | (np.abs(tdiff - t_mn) > any_bad_thresh * t_std)
+        )
+        if np.any(any_very_bad):
+            some_bad = (
+                (~np.isfinite(tg_arr))
+                | (np.abs(tdiff - t_mn) > flag_bad_thresh * t_std)
+            )
+            bad_msk = some_bad
+            b, flag = _make_a_good_box(some_bad, verbose=verbose)
+            flags |= flag
+        else:
+            bad_msk = np.zeros_like(tg_arr).astype(bool)
+    else:
+        bad_msk = np.ones_like(tg_arr).astype(bool)
+
+    return dict(
+        flags=flags,
+        t_star=ts_arr,
+        t_gal=tg_arr,
+        bad_msk=bad_msk,
+        bbox=_rescale_b(b, grid_size),
+    )
diff --git a/des_y6utils/piff/tests/__init__.py b/des_y6utils/piff/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/des_y6utils/piff/tests/test_flag_good_regions.py b/des_y6utils/piff/tests/test_flag_good_regions.py
new file mode 100644
index 0000000..70849d8
--- /dev/null
+++ b/des_y6utils/piff/tests/test_flag_good_regions.py
@@ -0,0 +1,86 @@
+import numpy as np
+
+from ..flag_good_regions import _make_a_good_box
+
+
+def test_flag_good_regions_make_a_good_box_smoke():
+    bad_msk = np.zeros((12, 15)).astype(bool)
+    b, flag = _make_a_good_box(bad_msk)
+
+    assert flag == 0
+    assert b == {
+        "xmin": 0,
+        "xmax": 14,
+        "ymin": 0,
+        "ymax": 11,
+    }
+    assert not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
+
+
+def test_flag_good_regions_make_a_good_box():
+    bad_msk = np.zeros((12, 15)).astype(bool)
+    bad_msk[:, 0] = True
+    b, flag = _make_a_good_box(bad_msk)
+
+    assert flag == 0
+    assert b == {
+        "xmin": 1,
+        "xmax": 14,
+        "ymin": 0,
+        "ymax": 11,
+    }
+    assert not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
+
+
+def test_flag_good_regions_make_a_good_box_cen():
+    bad_msk = np.zeros((13, 15)).astype(bool)
+    bad_msk[6, 7] = True
+    b, flag = _make_a_good_box(bad_msk)
+    assert flag == 0
+    assert b == {
+        "xmin": 0,
+        "xmax": 6,
+        "ymin": 0,
+        "ymax": 12,
+    }
+    assert not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
+
+
+def test_flag_good_regions_make_a_good_box_cen_inv():
+    bad_msk = np.ones((13, 15)).astype(bool)
+    bad_msk[6, 7] = False
+    b, flag = _make_a_good_box(bad_msk)
+    assert flag == 0
+    assert b == {
+        "xmin": 7,
+        "xmax": 7,
+        "ymin": 6,
+        "ymax": 6,
+    }
+    assert not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
+
+
+def test_flag_good_regions_make_a_good_box_random():
+    rng = np.random.RandomState(seed=10)
+    some_ok = False
+    for _ in range(100):
+        bad_msk = rng.uniform(size=(12, 15)) > 0.8
+        b, flag = _make_a_good_box(bad_msk)
+        if flag == 0:
+            assert not np.any(bad_msk[b["ymin"]:b["ymax"]+1, b["xmin"]:b["xmax"]+1])
+            some_ok = True
+
+    assert some_ok
+
+
+def test_flag_good_regions_make_a_good_box_all_bad():
+    bad_msk = np.ones((12, 15)).astype(bool)
+    b, flag = _make_a_good_box(bad_msk)
+
+    assert flag != 0
+    assert b == {
+        "xmin": 0,
+        "xmax": 14,
+        "ymin": 0,
+        "ymax": 11,
+    }

From 97b464d52ca1d490bb9eb5dd2461265609eddd8a Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Mon, 18 Dec 2023 13:27:52 -0600
Subject: [PATCH 11/20] BUG put back piff

---
 des_y6utils/__init__.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/des_y6utils/__init__.py b/des_y6utils/__init__.py
index d83bd64..e435846 100644
--- a/des_y6utils/__init__.py
+++ b/des_y6utils/__init__.py
@@ -12,3 +12,4 @@
 
 from . import mdet
 from . import viz
+from . import piff

From 0b891b1ceb0004814e3d3dea8b40d599c69a7c6d Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Thu, 21 Mar 2024 15:18:21 -0500
Subject: [PATCH 12/20] ENH move to latest mask

---
 des_y6utils/mdet.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index af5acd0..b51a21e 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -109,7 +109,7 @@ def _make_mdet_cuts_gauss(
     min_t_ratio=0.5,
     max_mfrac=0.1,
     max_s2n=np.inf,
-    mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v1.hsp",
+    mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v2.hsp",
     max_t=100.0,
 ):
     """
@@ -336,7 +336,7 @@ def _make_mdet_cuts_v6(d, verbose=False):
 
     # apply the mask
     hmap = _read_hsp_file(
-        "y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v1.hsp"
+        "y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v2.hsp"
     )
     in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
     msk &= in_footprint

From b35cab68f88cfce60247cb6f7ec8a45e6cfe5f55 Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Sat, 23 Mar 2024 08:51:41 -0500
Subject: [PATCH 13/20] Update des_y6utils/mdet.py

---
 des_y6utils/mdet.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index b51a21e..fae904a 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -47,10 +47,10 @@ def add_extinction_correction_columns(
     new_dt = np.dtype(
         data.dtype.descr
         + [
-            ('pgauss_band_flux_g_nodered', 'f8'),
-            ('pgauss_band_flux_r_nodered', 'f8'),
-            ('pgauss_band_flux_i_nodered', 'f8'),
-            ('pgauss_band_flux_z_nodered', 'f8'),
+            ('pgauss_band_flux_g_nodered', 'f4'),
+            ('pgauss_band_flux_r_nodered', 'f4'),
+            ('pgauss_band_flux_i_nodered', 'f4'),
+            ('pgauss_band_flux_z_nodered', 'f4'),
         ]
     )
     new_data = np.zeros(data.shape, dtype=new_dt)

From 0356c4a4cefe828bf7b899dce523b4ea11ff0fa1 Mon Sep 17 00:00:00 2001
From: Matthew R Becker <becker.mr@gmail.com>
Date: Mon, 1 Apr 2024 17:27:51 -0400
Subject: [PATCH 14/20] BUG change flux errors too

---
 des_y6utils/mdet.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index fae904a..3dc0a59 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -43,6 +43,9 @@ def add_extinction_correction_columns(
         flux_ = _compute_dered_flux(data["pgauss_band_flux_" + b], ii, dered)
         data["pgauss_band_flux_" + b] = flux_
 
+        flux_err_ = _compute_dered_flux(data["pgauss_band_flux_err_" + b], ii, dered)
+        data["pgauss_band_flux_err_" + b] = flux_err_
+
     # make _nodered array with pgauss_band_flux_* entries, and add them to fits.
     new_dt = np.dtype(
         data.dtype.descr

From 3cf4aa7b76cda2c7a61092916933b3ff0617e261 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Mon, 1 Apr 2024 16:37:22 -0500
Subject: [PATCH 15/20] REF try to use less memory

---
 des_y6utils/mdet.py                 | 31 ++++++++++++++---------------
 des_y6utils/tests/test_asinh_mag.py |  4 ++--
 2 files changed, 17 insertions(+), 18 deletions(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index 3dc0a59..f5ce561 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -38,13 +38,11 @@ def add_extinction_correction_columns(
     dered = dustmap.get_values_pos(data["ra"], data["dec"])
     flux_og = []
     for ii, b in enumerate(bands):
-        flux = np.copy(data["pgauss_band_flux_" + b])
-        flux_og.append(flux)
-        flux_ = _compute_dered_flux(data["pgauss_band_flux_" + b], ii, dered)
-        data["pgauss_band_flux_" + b] = flux_
+        dered_fac = _compute_dered_flux_fac(ii, dered)
+        flux_og.append(data["pgauss_band_flux_" + b].copy())
 
-        flux_err_ = _compute_dered_flux(data["pgauss_band_flux_err_" + b], ii, dered)
-        data["pgauss_band_flux_err_" + b] = flux_err_
+        data["pgauss_band_flux_" + b] *= dered_fac
+        data["pgauss_band_flux_err_" + b] *= dered_fac
 
     # make _nodered array with pgauss_band_flux_* entries, and add them to fits.
     new_dt = np.dtype(
@@ -381,18 +379,21 @@ def _compute_asinh_mags(flux, i):
     return mag
 
 
-def _compute_dered_flux(flux, i, map):
-    """This function applies dereddening correction to fluxes.
+def _compute_dered_flux_fac(i, ebv_map_val):
+    """This function computes the dereddening factor for the flux.
+
     Eli says we cannot apply the correction to the asinh mag.
 
+    You use like this:
+
+        dered_flux = flux * _compute_dered_flux_fac(i, map)
+
     Parameters
     ----------
-    flux : float or np.ndarray
-        The flux.
     i : int
         The index of the band in griz (i.e., 0 for g, 1 for r, 2 for i, 3 for z).
-    map : float or np.ndarray
-        The E(B-V).
+    ebv_map_val : float or np.ndarray
+        The E(B-V) value(s) from the dust map.
 
     Returns
     -------
@@ -401,10 +402,8 @@ def _compute_dered_flux(flux, i, map):
     """
 
     dered_fac = [3.186, 2.140, 1.196, 1.048]
-    dmag = dered_fac[i] * map
-    dered_flux = flux * 10**(dmag / 2.5)
-
-    return dered_flux
+    dmag = dered_fac[i] * ebv_map_val
+    return 10**(dmag / 2.5)
 
 
 def _compute_asinh_flux(mag, i):
diff --git a/des_y6utils/tests/test_asinh_mag.py b/des_y6utils/tests/test_asinh_mag.py
index 747c2f1..22cade2 100644
--- a/des_y6utils/tests/test_asinh_mag.py
+++ b/des_y6utils/tests/test_asinh_mag.py
@@ -1,4 +1,4 @@
-from ..mdet import _compute_dered_flux, _compute_asinh_flux, _compute_asinh_mags
+from ..mdet import _compute_dered_flux_fac, _compute_asinh_flux, _compute_asinh_mags
 import math
 
 
@@ -13,7 +13,7 @@ def test_compute_asinh_flux():
 
 def test_compute_dered_flux():
     flux_input = 10000
-    dered_flux = _compute_dered_flux(flux_input, 0, 2.5 / 3.186)
+    dered_flux = flux_input * _compute_dered_flux_fac(0, 2.5 / 3.186)
 
     assert math.isclose(flux_input*10, dered_flux, rel_tol=1e-9)
 

From deb0ff47620bd0709e41b3da25d9e85ff60dca8a Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 25 Apr 2024 10:26:13 -0500
Subject: [PATCH 16/20] ENH finalize new mask

---
 .DS_Store           | Bin 0 -> 6148 bytes
 des_y6utils/mdet.py | 286 +++++++++++++++++++++++---------------------
 2 files changed, 151 insertions(+), 135 deletions(-)
 create mode 100644 .DS_Store

diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..ec3865d02b5f95538e9a0ad624efb479aee3ca3e
GIT binary patch
literal 6148
zcmeHK!AiqG5S?wSO({YTiXIod7OaJ8i<c1V4;aydN^MBdV9b^#wMZ%Csz2nP_&v_-
zZpBi2@FG%XVD`<<&SaUFu#*J<;f=#XfEoZesDy<I4nGL3lP*cYdI*Jn#tR(57&3?;
zUW#VNe`J8(T@^wYKngwh{{GTnk`9Ad1sL%j4C5%x8jTN8D3&%h%TC#;I=B9Vn)q2i
z8>e1>bWO7hrQ%?2_k+u5FzMB{PgRojqhv7B2~jw}kejP03Dv|?<0MRVu5T8cic{&;
zcBj*JyDppU<5^u!JDp}-wpyLptm5qKADo=`9^<D(y&6UZetcTCEY9H#jTOZ``=caQ
z$vsBQv*%fi%m6dM46Ho^_TY1>YhNiB$P6$8zh!{V2Z>7PTFedVqXP%}eWdjgAqm=a
zmmriDU5mLvTtN{g715*$d&LkY9sSbAxfXMSCLM&{8J}Z!7WRfB^zP`FIvj*+kVj^K
z8CYbXXr?vl|7WZ3|BFdHV+NRkwPHY&x_;Ngl5B5XDvo-sM7={Lp}5@OR|*=s6=N*5
d;yqL?=$B+5x)yVT=t1Eh0ZjuB%)p;A@CkmtPD20y

literal 0
HcmV?d00001

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index f5ce561..1666b13 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -110,7 +110,7 @@ def _make_mdet_cuts_gauss(
     min_t_ratio=0.5,
     max_mfrac=0.1,
     max_s2n=np.inf,
-    mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v2.hsp",
+    mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-extra-masks-v2.hsp",
     max_t=100.0,
 ):
     """
@@ -145,6 +145,125 @@ def _make_mdet_cuts_gauss(
     return msk
 
 
+def _make_mdet_cuts_v6(d, verbose=False):
+
+    msk = _make_mdet_cuts_raw_v6(
+        d,
+        verbose=verbose,
+        min_s2n=10,
+        n_terr=0,
+        min_t_ratio=0.5,
+        max_mfrac=0.1,
+        max_s2n=np.inf,
+        max_t=20.0,
+    )
+
+    # apply the mask
+    hmap = _read_hsp_file(
+        "y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-extra-masks-v2.hsp"
+    )
+    in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
+    msk &= in_footprint
+    if verbose:
+        print("did mask cuts", np.sum(msk))
+
+    return msk
+
+
+def _make_mdet_cuts_raw_v6(
+    d,
+    *,
+    min_s2n,
+    n_terr,
+    min_t_ratio,
+    max_mfrac,
+    max_s2n,
+    max_t,
+    verbose=False,
+):
+    """The raw v6 cuts come from analysis in Fall of 2023.
+
+      - We implemented i-band magnitude cut at 24.7.
+      - We implemented galaxy size cut at T=20.
+      - We updated a cut in the pgauss T-Terr plane.
+    """
+
+    mag_g = _compute_asinh_mags(d["pgauss_band_flux_g"], 0)
+    mag_r = _compute_asinh_mags(d["pgauss_band_flux_r"], 1)
+    mag_i = _compute_asinh_mags(d["pgauss_band_flux_i"], 2)
+    mag_z = _compute_asinh_mags(d["pgauss_band_flux_z"], 3)
+
+    gmr = mag_g - mag_r
+    rmi = mag_r - mag_i
+    imz = mag_i - mag_z
+
+    msk = np.ones(d.shape[0]).astype(bool)
+
+    potential_flag_columns = [
+        "psfrec_flags",
+        "gauss_flags",
+        "pgauss_T_flags",
+        "pgauss_band_flux_flags_g",
+        "pgauss_band_flux_flags_r",
+        "pgauss_band_flux_flags_i",
+        "pgauss_band_flux_flags_z",
+        "mask_flags",
+    ]
+    for col in potential_flag_columns:
+        if col in d.dtype.names:
+            msk &= (d[col] == 0)
+            if verbose:
+                print("did cut " + col, np.sum(msk))
+
+    if "shear_bands" in d.dtype.names:
+        msk &= (d["shear_bands"] == "123")
+        if verbose:
+            print("did cut shear_bands", np.sum(msk))
+
+    if "pgauss_s2n" in d.dtype.names:
+        msk &= (d["pgauss_s2n"] > 5)
+        if verbose:
+            print("did cut pgauss_s2n", np.sum(msk))
+
+    # now do the rest
+    msk &= (
+        (d["gauss_s2n"] > min_s2n)
+        & (d["gauss_s2n"] < max_s2n)
+        & (d["mfrac"] < max_mfrac)
+        & (np.abs(gmr) < 5)
+        & (np.abs(rmi) < 5)
+        & (np.abs(imz) < 5)
+        & np.isfinite(mag_g)
+        & np.isfinite(mag_r)
+        & np.isfinite(mag_i)
+        & np.isfinite(mag_z)
+        & (mag_g < 26.5)
+        & (mag_r < 26.5)
+        & (mag_i < 24.7)  # used to eliminate photo-z outliers motivated by cosmos2020
+                          # the cut is 24.5 in cosmos mags with + 0.2 mags to correct
+                          # for the pgauss aperture
+        & (mag_z < 25.6)
+        & (d["pgauss_T"] < (1.6 - 3.1 * d["pgauss_T_err"]))
+        & (
+            d["gauss_T_ratio"] >= np.maximum(
+                min_t_ratio,
+                (n_terr * d["gauss_T_err"] / d["gauss_psf_T"])
+            )
+        )
+        & ((d["gauss_T_ratio"] * d["gauss_psf_T"]) < max_t)
+    )
+
+    size_sizeerr = (d['gauss_T_ratio'] * d['gauss_psf_T']) * d['gauss_T_err']
+    size_s2n = (d['gauss_T_ratio'] * d['gauss_psf_T']) / d['gauss_T_err']
+    msk_superspreader = ((size_sizeerr > 1) & (size_s2n < 10))
+    msk &= ~msk_superspreader
+
+    if verbose:
+        print("did mdet cuts", np.sum(msk))
+
+    return msk
+
+
 def _make_mdet_cuts_wmom(
     data,
     verbose=False,
@@ -317,36 +436,6 @@ def _make_mdet_cuts_v5(d, verbose=False):
     return msk
 
 
-def _make_mdet_cuts_v6(d, verbose=False):
-
-    msk = _make_mdet_cuts_raw_v6(
-        d,
-        verbose=verbose,
-        min_s2n=10,
-        n_terr=0,
-        min_t_ratio=0.5,
-        max_mfrac=0.1,
-        max_s2n=np.inf,
-        max_t=20.0,
-    )
-
-    size_sizeerr = (d['gauss_T_ratio'] * d['gauss_psf_T']) * d['gauss_T_err']
-    size_s2n = (d['gauss_T_ratio'] * d['gauss_psf_T']) / d['gauss_T_err']
-    msk_superspreader = ((size_sizeerr > 1) & (size_s2n < 10))
-    msk &= ~msk_superspreader
-
-    # apply the mask
-    hmap = _read_hsp_file(
-        "y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v2.hsp"
-    )
-    in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
-    msk &= in_footprint
-    if verbose:
-        print("did mask cuts", np.sum(msk))
-
-    return msk
-
-
 def _compute_asinh_mags(flux, i):
     """This function and coefficients are from from Eli. Ask him.
 
@@ -465,25 +554,39 @@ def _download_fname_from_bnl(fpth):
     wget_res = subprocess.run("which wget", shell=True, capture_output=True)
     curl_res = subprocess.run("which curl", shell=True, capture_output=True)
 
-    bnl = "https://www.cosmo.bnl.gov/www/esheldon/data/y6-healsparse"
+    base_urls_to_try = [
+        "https://www.cosmo.bnl.gov/www/beckermr/des/y6/masks",
+        "https://www.cosmo.bnl.gov/www/esheldon/data/y6-healsparse",
+    ]
+
     if wget_res.returncode == 0:
-        subprocess.run(
-            "cd %s && wget %s/%s" % (
-                fdir, bnl, fname,
-            ),
-            shell=True,
-            check=True,
-            capture_output=True,
-        )
+        for bnl in base_urls_to_try:
+            ret = subprocess.run(
+                "cd %s && wget %s/%s" % (
+                    fdir, bnl, fname,
+                ),
+                shell=True,
+                capture_output=True,
+            )
+            if ret.returncode == 0:
+                break
+
+        ret.check_returncode()
+
     elif curl_res.returncode == 0:
-        subprocess.run(
-            "cd %s && curl -L %s/%s --output %s" % (
-                fdir, bnl, fname, fname,
-            ),
-            shell=True,
-            check=True,
-            capture_output=True,
-        )
+        for bnl in base_urls_to_try:
+            ret = subprocess.run(
+                "cd %s && curl -L %s/%s --output %s" % (
+                    fdir, bnl, fname, fname,
+                ),
+                shell=True,
+                capture_output=True,
+            )
+            if ret.returncode == 0:
+                break
+
+        ret.check_returncode()
+
     else:
         raise RuntimeError(
             "Could not download mask '%s' from BNL due "
@@ -670,90 +773,3 @@ def _make_mdet_cuts_raw_v345(
     return msk
 
 
-def _make_mdet_cuts_raw_v6(
-    d,
-    *,
-    min_s2n,
-    n_terr,
-    min_t_ratio,
-    max_mfrac,
-    max_s2n,
-    max_t,
-    verbose=False,
-):
-    """The raw v6 cuts come from analysis in Fall of 2023.
-
-      - We implemented i-band magnitude cut at 24.7.
-      - We implemented galaxy size cut at T=20.
-      - We updated a cut in the pgauss T-Terr plane.
-    """
-
-    mag_g = _compute_asinh_mags(d["pgauss_band_flux_g"], 0)
-    mag_r = _compute_asinh_mags(d["pgauss_band_flux_r"], 1)
-    mag_i = _compute_asinh_mags(d["pgauss_band_flux_i"], 2)
-    mag_z = _compute_asinh_mags(d["pgauss_band_flux_z"], 3)
-
-    gmr = mag_g - mag_r
-    rmi = mag_r - mag_i
-    imz = mag_i - mag_z
-
-    msk = np.ones(d.shape[0]).astype(bool)
-
-    potential_flag_columns = [
-        "psfrec_flags",
-        "gauss_flags",
-        "pgauss_T_flags",
-        "pgauss_band_flux_flags_g",
-        "pgauss_band_flux_flags_r",
-        "pgauss_band_flux_flags_i",
-        "pgauss_band_flux_flags_z",
-        "mask_flags",
-    ]
-    for col in potential_flag_columns:
-        if col in d.dtype.names:
-            msk &= (d[col] == 0)
-            if verbose:
-                print("did cut " + col, np.sum(msk))
-
-    if "shear_bands" in d.dtype.names:
-        msk &= (d["shear_bands"] == "123")
-        if verbose:
-            print("did cut shear_bands", np.sum(msk))
-
-    if "pgauss_s2n" in d.dtype.names:
-        msk &= (d["pgauss_s2n"] > 5)
-        if verbose:
-            print("did cut pgauss_s2n", np.sum(msk))
-
-    # now do the rest
-    msk &= (
-        (d["gauss_s2n"] > min_s2n)
-        & (d["gauss_s2n"] < max_s2n)
-        & (d["mfrac"] < max_mfrac)
-        & (np.abs(gmr) < 5)
-        & (np.abs(rmi) < 5)
-        & (np.abs(imz) < 5)
-        & np.isfinite(mag_g)
-        & np.isfinite(mag_r)
-        & np.isfinite(mag_i)
-        & np.isfinite(mag_z)
-        & (mag_g < 26.5)
-        & (mag_r < 26.5)
-        & (mag_i < 24.7)  # used to eliminate photo-z outliers motivated by cosmos2020
-                          # the cut is 24.5 in cosmos mags with + 0.2 mags to correct
-                          # for the pgauss aperture
-        & (mag_z < 25.6)
-        & (d["pgauss_T"] < (1.6 - 3.1 * d["pgauss_T_err"]))
-        & (
-            d["gauss_T_ratio"] >= np.maximum(
-                min_t_ratio,
-                (n_terr * d["gauss_T_err"] / d["gauss_psf_T"])
-            )
-        )
-        & ((d["gauss_T_ratio"] * d["gauss_psf_T"]) < max_t)
-    )
-
-    if verbose:
-        print("did mdet cuts", np.sum(msk))
-
-    return msk

From d7150ecfa0c70e51e5f447c52eaf3e8558d6b680 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 25 Apr 2024 10:26:48 -0500
Subject: [PATCH 17/20] BUG add git ignore

---
 .DS_Store  | Bin 6148 -> 0 bytes
 .gitignore |   2 ++
 2 files changed, 2 insertions(+)
 delete mode 100644 .DS_Store

diff --git a/.DS_Store b/.DS_Store
deleted file mode 100644
index ec3865d02b5f95538e9a0ad624efb479aee3ca3e..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 6148
zcmeHK!AiqG5S?wSO({YTiXIod7OaJ8i<c1V4;aydN^MBdV9b^#wMZ%Csz2nP_&v_-
zZpBi2@FG%XVD`<<&SaUFu#*J<;f=#XfEoZesDy<I4nGL3lP*cYdI*Jn#tR(57&3?;
zUW#VNe`J8(T@^wYKngwh{{GTnk`9Ad1sL%j4C5%x8jTN8D3&%h%TC#;I=B9Vn)q2i
z8>e1>bWO7hrQ%?2_k+u5FzMB{PgRojqhv7B2~jw}kejP03Dv|?<0MRVu5T8cic{&;
zcBj*JyDppU<5^u!JDp}-wpyLptm5qKADo=`9^<D(y&6UZetcTCEY9H#jTOZ``=caQ
z$vsBQv*%fi%m6dM46Ho^_TY1>YhNiB$P6$8zh!{V2Z>7PTFedVqXP%}eWdjgAqm=a
zmmriDU5mLvTtN{g715*$d&LkY9sSbAxfXMSCLM&{8J}Z!7WRfB^zP`FIvj*+kVj^K
z8CYbXXr?vl|7WZ3|BFdHV+NRkwPHY&x_;Ngl5B5XDvo-sM7={Lp}5@OR|*=s6=N*5
d;yqL?=$B+5x)yVT=t1Eh0ZjuB%)p;A@CkmtPD20y

diff --git a/.gitignore b/.gitignore
index b6e4761..8ba67bc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -127,3 +127,5 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+**/.DS_Store

From 968a7c3719349168e189fc8cf587dcbb85458a65 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 25 Apr 2024 10:35:14 -0500
Subject: [PATCH 18/20] STY please the flake8

---
 des_y6utils/mdet.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index 1666b13..101226e 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -771,5 +771,3 @@ def _make_mdet_cuts_raw_v345(
         print("did mdet cuts", np.sum(msk))
 
     return msk
-
-

From 4f49f26d7686424760bf1e3611365ece093c88d6 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 25 Apr 2024 10:39:29 -0500
Subject: [PATCH 19/20] REF use one location for default mask

---
 des_y6utils/mdet.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index 101226e..b080d45 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -7,6 +7,9 @@
 import healsparse
 
 
+DEFAULT_HSP_MASK = "y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-extra-masks-v2.hsp"  # noqa: E501
+
+
 def add_extinction_correction_columns(
     data, dust_map_filename="SFD_dust_4096.hsp"
 ):
@@ -110,7 +113,7 @@ def _make_mdet_cuts_gauss(
     min_t_ratio=0.5,
     max_mfrac=0.1,
     max_s2n=np.inf,
-    mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-extra-masks-v2.hsp",
+    mask_name=DEFAULT_HSP_MASK,
     max_t=100.0,
 ):
     """
@@ -159,9 +162,7 @@ def _make_mdet_cuts_v6(d, verbose=False):
     )
 
     # apply the mask
-    hmap = _read_hsp_file(
-        "y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-extra-masks-v2.hsp"
-    )
+    hmap = _read_hsp_file(DEFAULT_HSP_MASK)
     in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True)
     msk &= in_footprint
     if verbose:

From 7e2065be39865cbf6fc82d048827bd99b92fefdd Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 25 Apr 2024 11:04:41 -0500
Subject: [PATCH 20/20] REF no need for this one now

---
 des_y6utils/mdet.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py
index b080d45..55e1a47 100644
--- a/des_y6utils/mdet.py
+++ b/des_y6utils/mdet.py
@@ -556,7 +556,6 @@ def _download_fname_from_bnl(fpth):
     curl_res = subprocess.run("which curl", shell=True, capture_output=True)
 
     base_urls_to_try = [
-        "https://www.cosmo.bnl.gov/www/beckermr/des/y6/masks",
         "https://www.cosmo.bnl.gov/www/esheldon/data/y6-healsparse",
     ]