Fix issues related to large p_acq change for small T_ccd change

docs/index.rst

@@ -311,6 +311,9 @@ The following environment variables are used by proseco:
   If this is a relative path then it is relative to ``<default_agasc_dir>``.
 - ``AGASC_SUPPLEMENT_ENABLED``: set to ``"False"`` to disable using the AGASC
    supplement. This is for testing and should not be used in production.
+- ``PROSECO_DISABLE_BOX_EDGE_DMAG``: If set to ``"True"`` then disable the
+  box edge dmag offset in the acq star selection. This is for testing and should not
+  be used in production.
 - ``PROSECO_DISABLE_OVERLAP_PENALTY``: if set to ``"True"`` then disable the
   overlap penalty in the acq star selection. This is for testing and should not
   be used in production.

proseco/acq.py

@@ -45,6 +45,10 @@
 OVERLAP_MAG_DEADBAND = 0.2  # overlap penalty applies for mag difference > deadband
 OVERLAP_PAD = 20  # arcsec, extra padding for overlap check
 
+# Constants in get_box_edge_dmag() function
+BOX_EDGE_DMAG_OFFSETS = np.array([-20.0, -10.0, -5.0, 0.0, 5.0, 12.5, 17.0, 20.0])
+BOX_EDGE_DMAG_VALUES = np.array([0.0, 0.03, 0.18, 0.75, 2.03, 5.12, 6.28, 10.0])
+
 
 def load_maxmags() -> dict:
     """
@@ -450,7 +454,7 @@ def update_idxs_halfws(self, idxs, halfws):
         """
         Update the rows of self to match the specified ``agasc_ids``
         and half widths.  These two input lists must match the length
-        of self and correspond to stars in self.cand_acqs.
+        of acqs and correspond to stars in self.cand_acqs.
 
         :param agasc_ids: list of AGASC IDs
         :param halfws: list of search box half widths
@@ -1282,6 +1286,22 @@ def get_spoiler_stars(stars, acq, box_size):
     return spoilers
 
 
+def force_even_and_on_ccd(val: int, offset: int) -> int:
+    """Force val to be even and on the CCD (0-1024).
+
+    Parameters
+    ----------
+    val : int
+        Input value.
+    offset : int
+        If ``val`` is odd then add ``sign(offset)`` to make even.
+    """
+    if val % 2 == 1:
+        val += np.sign(offset)
+    val = np.clip(val, 0, 1024)
+    return val
+
+
 def get_imposter_stars(
     dark,
     star_row,
@@ -1320,26 +1340,10 @@ def get_imposter_stars(
     box_col = int(box_size.col)
 
     # Make sure box is within CCD
-    box_r0 = np.clip(acq_row - box_row, 0, 1024)
-    box_r1 = np.clip(acq_row + box_row, 0, 1024)
-    box_c0 = np.clip(acq_col - box_col, 0, 1024)
-    box_c1 = np.clip(acq_col + box_col, 0, 1024)
-
-    # Make sure box has even number of pixels on each edge.  Increase
-    # box by one if needed.
-    #
-    # TO DO: Test the clipping and shrinking code
-    #
-    if (box_r1 - box_r0) % 2 == 1:
-        if box_r1 == 1024:
-            box_r0 -= 1
-        else:
-            box_r1 += 1
-    if (box_c1 - box_c0) % 2 == 1:
-        if box_c1 == 1024:
-            box_c0 -= 1
-        else:
-            box_c1 += 1
+    box_r0 = force_even_and_on_ccd(acq_row - box_row, -1)
+    box_r1 = force_even_and_on_ccd(acq_row + box_row, 1)
+    box_c0 = force_even_and_on_ccd(acq_col - box_col, -1)
+    box_c1 = force_even_and_on_ccd(acq_col + box_col, 1)
 
     # Get bgd-subtracted dark current image corresponding to the search box
     # and bin in 2x2 blocks.
@@ -1517,14 +1521,59 @@ def get_intruders(acq, box_size, name, n_sigma, get_func, kwargs):
 
     colnames = ["yang", "zang", "mag", "mag_err"]
     if len(intruders) == 0:
-        intruders = {name: np.array([], dtype=np.float64) for name in colnames}
+        out = {name: np.array([], dtype=np.float64) for name in colnames}
     else:
-        ok = (np.abs(intruders["yang"] - acq["yang"]) < box_size.y) & (
-            np.abs(intruders["zang"] - acq["zang"]) < box_size.z
+        dys = intruders["yang"] - acq["yang"]
+        dzs = intruders["zang"] - acq["zang"]
+        box_margin = (
+            0.0 if os.environ.get("PROSECO_DISABLE_BOX_EDGE_DMAG") == "True" else 20.0
         )
-        intruders = {name: intruders[name][ok] for name in ["mag", "mag_err"]}
+        ok = (np.abs(dys) < box_size.y + box_margin) & (
+            np.abs(dzs) < box_size.z + box_margin
+        )
+        intruders = intruders[ok]
+        mags = []
+        for dy, dz, mag in zip(dys, dzs, intruders["mag"]):
+            mag += get_box_edge_dmag(dy, box_size.y)  # noqa: PLW2901
+            mag += get_box_edge_dmag(dz, box_size.z)  # noqa: PLW2901
+            mags.append(mag)
+        out = {
+            "yang": intruders["yang"],
+            "zang": intruders["zang"],
+            "mag": np.array(mags),
+            "mag_err": intruders["mag_err"],
+        }
+
+    return out
+
+
+def get_box_edge_dmag(dyz, box_size):
+    """
+    Calculate the magnitude penalty for an intruder near the edge of a search box.
+
+    This uses linear interpolation of an empirical relationship determined in the
+    ``pr409-intruders-at-search-box-edge.ipynb`` notebook.
 
-    return intruders
+    Parameters
+    ----------
+    dyz : float
+        Offset from the box center to the intruder (arcsec).
+    box_size : float
+        Size of the box (arcsec).
+
+    Returns
+    -------
+    dmag : float
+        Magnitude penalty for being near the box edge.
+    """
+    if os.environ.get("PROSECO_DISABLE_BOX_EDGE_DMAG") == "True":
+        return 0.0
+
+    d_box_edge = abs(dyz) - box_size
+    # Note: np.interp is at least 5x faster than scipy.interpolate.interp1d (even with
+    # that returning a lambda func).
+    dmag = np.interp(x=d_box_edge, xp=BOX_EDGE_DMAG_OFFSETS, fp=BOX_EDGE_DMAG_VALUES)
+    return dmag
 
 
 def calc_p_on_ccd(row, col, box_size):

proseco/tests/conftest.py

@@ -17,6 +17,11 @@ def disable_fid_offsets(monkeypatch):
     monkeypatch.setenv("PROSECO_ENABLE_FID_OFFSET", "False")
 
 
+@pytest.fixture()
+def disable_box_edge_dmag(monkeypatch):
+    monkeypatch.setenv("PROSECO_DISABLE_BOX_EDGE_DMAG", "True")
+
+
 @pytest.fixture(autouse=True)
 def use_fixed_chandra_models(monkeypatch):
     monkeypatch.setenv("CHANDRA_MODELS_DEFAULT_VERSION", "3.48")

proseco/tests/test_acq.py

@@ -207,7 +207,7 @@ def get_test_cand_acqs():
     return CACHE["cand_acqs"].copy()
 
 
-def test_calc_p_brightest_same_bright():
+def test_calc_p_brightest_same_bright(disable_box_edge_dmag):
     """
     Test for an easy situation of three spoiler/imposters with exactly
     the same brightness as acq star so that p_brighter is always 0.5.
@@ -239,7 +239,7 @@ def test_calc_p_brightest_same_bright():
     assert np.allclose(probs, [0.25, 0.25, 0.25, 0.3334, 0.5, 1.0], rtol=0, atol=0.01)
 
 
-def test_calc_p_brightest_same_bright_asymm_dither():
+def test_calc_p_brightest_same_bright_asymm_dither(disable_box_edge_dmag):
     """
     Test for an easy situation of three spoiler/imposters with exactly
     the same brightness as acq star so that p_brighter is always 0.5.
@@ -308,7 +308,7 @@ def test_calc_p_brightest_same_bright_asymm_dither():
     assert np.allclose(probs, [0.333, 0.333, 0.333, 0.5, 1.0, 1.0], rtol=0, atol=0.01)
 
 
-def test_calc_p_brightest_1mag_brighter():
+def test_calc_p_brightest_1mag_brighter(disable_box_edge_dmag):
     """
     Test for the situation of spoiler/imposter that is 1 mag brighter.
     """
@@ -541,7 +541,9 @@ def test_get_acq_catalog_21007(proseco_agasc_1p7, disable_overlap_penalty):
     assert info == exp
 
 
-def test_box_strategy_20603(proseco_agasc_1p7, disable_overlap_penalty):
+def test_box_strategy_20603(
+    proseco_agasc_1p7, disable_overlap_penalty, disable_box_edge_dmag
+):
     """Test for PR #32 that doesn't allow p_acq to be reduced below 0.1.
 
     The idx=8 (mag=10.50) star was previously selected with 160 arsec box.

proseco/tests/test_catalog.py

@@ -109,12 +109,12 @@ def test_get_aca_catalog_20603(proseco_agasc_1p7):
         "   1   2         5  FID 8x8 -1829.63   156.96   1   1    25",
         "   2   3 116791824  BOT 6x6   622.00  -953.60  28   1   160",
         "   3   4  40114416  BOT 6x6   394.22  1204.43  24   1   140",
-        "   4   5  40112304  BOT 6x6 -1644.35  2032.47  12   1    80",
+        "   4   5  40112304  BOT 6x6 -1644.35  2032.47  16   1   100",
         "   5   6  40113544  GUI 6x6   102.74  1133.37   1   1    25",
         "   0   7 116923496  ACQ 6x6 -1337.79  1049.27  20   1   120",
         "   1   8 116923528  ACQ 6x6 -2418.65  1088.40  28   1   160",
         "   5   9 116791744  ACQ 6x6   985.38 -1210.19  28   1   160",
-        "   6  10  40108048  ACQ 6x6     2.21  1619.17  24   1   140",
+        "   6  10  40108048  ACQ 6x6     2.21  1619.17  16   1   100",
     ]
 
     assert aca[TEST_COLS].pformat(max_width=-1) == exp
@@ -653,16 +653,22 @@ def test_bad_obsid():
     assert "ValueError: text does not have OBSID" in aca.exception
 
 
-def test_bad_pixel_dark_current():
+@pytest.mark.parametrize("no_box_edge_dmag", [True, False])
+def test_bad_pixel_dark_current(monkeypatch, no_box_edge_dmag):
     """
     Test avoidance of bad_pixels = [[-245, 0, 454, 454]]
 
     - Put a bright star near this bad column and confirm it is not picked at
       all.
     - Put a bright star at col = 454 - 20 / 5 - 105 / 5 (dither and search box
-      size) and confirm it is picked with search box = 100 arcsec.
+      size) and confirm it is picked with appropriate search box size:
+       - 100 arcsec if the box edge dmag offset is disabled
+       - 80 arcsec if the box edge dmag offset is enabled (default). This means that
+         the bad pixel imposter is leaking into the search box as intended.
 
     """
+    if no_box_edge_dmag:
+        monkeypatch.setenv("PROSECO_DISABLE_BOX_EDGE_DMAG", "True")
     dark = DARK40.copy()
     stars = StarsTable.empty()
     stars.add_fake_constellation(mag=np.linspace(8.0, 8.1, 4), n_stars=4)
@@ -687,7 +693,8 @@ def test_bad_pixel_dark_current():
     exp_ids = [2, 100, 101, 102, 103]
     assert sorted(aca.guides["id"]) == sorted(exp_ids + [4])
     assert aca.acqs["id"].tolist() == exp_ids
-    assert aca.acqs["halfw"].tolist() == [100, 160, 160, 160, 160]
+    exp_halfw = 100 if no_box_edge_dmag else 80
+    assert aca.acqs["halfw"].tolist() == [exp_halfw, 160, 160, 160, 160]
 
 
 def test_fid_trap_effect():

proseco/tests/test_jupiter.py

@@ -272,7 +272,7 @@ def test_jupiter_midline():
 
 
 @pytest.mark.parametrize("col_dist_arcsec", np.arange(-300, 305, 20))
-def test_jupiter_acquisition(col_dist_arcsec):
+def test_jupiter_acquisition(col_dist_arcsec, disable_box_edge_dmag):
     """
     Test how jupiter is handled during acquisition.
 

ruff.toml

@@ -2,7 +2,11 @@ extend = "ruff-base.toml"
 
 # These are files to exclude for this project.
 extend-exclude = [
-  "**/*.ipynb",  # commonly not ruff-compliant
+  "analysis/pr388-backstop-vs-proseco-star-positions.ipynb",
+  "analysis/pr388-validate-drift-model-fid-positions.ipynb",
+  "analysis/maxmag/interpolate-maxmags-contour-plot.ipynb",
+  "analysis/maxmag/maxmags-demo.ipynb",
+  "analysis/maxmag/validate-maxmags-recent-loads.ipynb",
 ]
 
 # These are rules that commonly cause many ruff warnings. Code will be improved by