Merge pull request #190 from sot/use-p2

Use P2 for calc_p_safe and optimization metric

proseco/acq.py

@@ -601,7 +601,14 @@ def calc_p_safe(self, verbose=False):
         Calculate the probability of a safing action resulting from failure
         to acquire at least two (2) acquisition stars.
 
+        This uses the probability of 2 or fewer stars => "conservative" p_fail at this
+        man_err.  This uses 2 stars instead of 1 or fewer (which is the actual criteria
+        for a safing action).  This allows for one star to be dropped for reasons not
+        reflected in the acq model probability and makes the optimization dig a bit deeper
+        in to the catalog beyond the brightest stars.
+
         :returns: p_safe (float)
+
         """
 
         p_no_safe = 1.0
@@ -612,14 +619,17 @@ def calc_p_safe(self, verbose=False):
 
             p_acqs = [acq['probs'].p_acqs(acq['halfw'], man_err) for acq in self]
 
-            p_n_cum = prob_n_acq(p_acqs)[1]
+            p_n_cum = prob_n_acq(p_acqs)[1]  # This returns (p_n, p_n_cum)
+
+            # Probability of 2 or fewer stars => conservative fail criteria
+            p2 = p_n_cum[2]
+
             if verbose:
                 self.log(f'man_err = {man_err}, p_man_err = {p_man_err}')
                 self.log('p_acqs =' + ' '.join(['{:.3f}'.format(val) for val in p_acqs]))
-                self.log('log10(p 1_or_fewer) = {:.2f}'.format(np.log10(p_n_cum[1])))
-            p_01 = p_n_cum[1]  # 1 or fewer => p_fail at this man_err
+                self.log('log10(p 2_or_fewer) = {:.2f}'.format(np.log10(p2)))
 
-            p_no_safe *= (1 - p_man_err * p_01)
+            p_no_safe *= (1 - p_man_err * p2)
 
         p_safe = 1 - p_no_safe
         self.p_safe = p_safe

proseco/tests/test_acq.py

@@ -502,7 +502,7 @@ def test_box_strategy_20603():
            '    4     4  40112304 -1644.35  2032.47   160',
            '    5     5 116923528 -2418.65  1088.40   160',
            '    6     6 116791744   985.38 -1210.19   160',
-           '    7     7  40108048     2.21  1619.17   160']
+           '    7     7  40108048     2.21  1619.17    60']
 
     assert repr(acqs[TEST_COLS]).splitlines() == exp
 
@@ -753,7 +753,7 @@ def test_acq_fid_catalog_probs_low_level():
 
     # Initial fid set is empty () and we check baseline p_safe
     assert acqs.fid_set == ()
-    assert np.allclose(np.log10(acqs.calc_p_safe()), -4.0,
+    assert np.allclose(np.log10(acqs.calc_p_safe()), -2.4,
                        rtol=0, atol=0.1)
 
     # This is the acq star spoiled by fid_id=2
@@ -777,7 +777,7 @@ def test_acq_fid_catalog_probs_low_level():
     # spoils an acq star.  This makes the p_safe value much worse.
     acqs.fid_set = (4, 3, 2)
     assert acqs.fid_set == (2, 3, 4)  # gets sorted when set
-    assert np.allclose(np.log10(acqs.calc_p_safe()), -2.6,
+    assert np.allclose(np.log10(acqs.calc_p_safe()), -1.3,
                        rtol=0, atol=0.1)
 
     # With fid_set = (1, 2, 4), the probability multiplier for catalog
@@ -790,7 +790,7 @@ def test_acq_fid_catalog_probs_low_level():
     # set does not spoil an acq star.
     for fid_set in ((1, 3, 4), ()):
         acqs.fid_set = fid_set
-        assert np.allclose(np.log10(acqs.calc_p_safe()), -4.0,
+        assert np.allclose(np.log10(acqs.calc_p_safe()), -2.4,
                            rtol=0, atol=0.1)
 
     # Check that p_acqs() method responds to fid_set in expected way
@@ -925,7 +925,7 @@ def test_acq_fid_catalog_one_cand_fid():
     aca = get_aca_catalog(**kwargs)
 
     assert aca.acqs['id'].tolist() == [1, 100, 101, 102, 103]
-    assert aca.acqs['halfw'].tolist() == [80, 160, 160, 160, 140]
+    assert aca.acqs['halfw'].tolist() == [80, 160, 160, 160, 120]
 
     assert aca.n_fid == 3
     assert aca.fids['id'].tolist() == [6]
@@ -976,7 +976,7 @@ def test_acq_fid_catalog_two_cand_fid(n_fid):
     aca = get_aca_catalog(**kwargs)
 
     assert aca.acqs['id'].tolist() == [1, 100, 101, 102, 103]
-    assert aca.acqs['halfw'].tolist() == [80, 160, 160, 160, 140]
+    assert aca.acqs['halfw'].tolist() == [80, 160, 160, 160, 120]
 
     assert aca.n_fid == n_fid
     assert aca.fids['id'].tolist() == [1, 2]

proseco/tests/test_catalog.py

@@ -43,14 +43,14 @@ def test_get_aca_catalog_20603():
            '---- --- --------- ---- --- -------- -------- --- --- -----',
            '   0   1         4  FID 8x8  2140.23   166.63   1   1    25',
            '   1   2         5  FID 8x8 -1826.28   160.17   1   1    25',
-           '   2   3  40112304  BOT 6x6 -1644.35  2032.47  20   1   160',
+           '   2   3  40112304  BOT 6x6 -1644.35  2032.47  20   1    80',
            '   3   4  40114416  BOT 6x6   394.22  1204.43  20   1   140',
            '   4   5 116791824  BOT 6x6   622.00  -953.60  20   1   160',
            '   5   6 116923528  BOT 6x6 -2418.65  1088.40  20   1   160',
            '   6   7  40113544  GUI 6x6   102.74  1133.37   1   1    25',
            '   6   8 116923496  ACQ 6x6 -1337.79  1049.27  20   1   120',
-           '   7   9 116791744  ACQ 6x6   985.38 -1210.19  20   1   160',
-           '   0  10  40108048  ACQ 6x6     2.21  1619.17  20   1    60']
+           '   7   9 116791744  ACQ 6x6   985.38 -1210.19  20   1   140',
+           '   0  10  40108048  ACQ 6x6     2.21  1619.17  20   1   140']
 
     repr(aca)  # Apply default formats
     assert aca[TEST_COLS].pformat(max_width=-1) == exp