Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Set a temperature dependent mag limit for imposters #118

Merged
merged 3 commits into from
Sep 29, 2018
Merged

Set a temperature dependent mag limit for imposters #118

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6681955
Make calc_p_brightest into an AcqTable method
taldcroft Sep 26, 2018
0f708ac
Only include statistical outliers in imposters list
taldcroft Sep 26, 2018
92a830f
Fix regression tests
taldcroft Sep 28, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
115 changes: 64 additions & 51 deletions proseco/acq.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,8 @@
from astropy.table import Table

from chandra_aca.star_probs import acq_success_prob, prob_n_acq
from chandra_aca.transform import (pixels_to_yagzag, mag_to_count_rate)
from chandra_aca.transform import (pixels_to_yagzag, mag_to_count_rate,
snr_mag_for_t_ccd)
from mica.archive.aca_dark.dark_cal import get_dark_cal_image

from . import characteristics as CHAR
Expand Down Expand Up @@ -56,6 +57,13 @@ def get_acq_catalog(obsid=0, **kwargs):
acqs.set_attrs_from_kwargs(obsid=obsid, **kwargs)
acqs.set_stars()

# Only allow imposters that are statistical outliers and are brighter than
# this (temperature-dependent) threshold. See characterisics.py for more
# explanation.
acqs.imposters_mag_limit = snr_mag_for_t_ccd(acqs.t_ccd,
ref_mag=CHAR.imposter_mag_lim_ref_mag,
ref_t_ccd=CHAR.imposter_mag_lim_ref_t_ccd)

acqs.log(f'getting dark cal image at date={acqs.date} t_ccd={acqs.t_ccd:.1f}')

# If dark map not provided via input kwarg then get from dark cal archive
Expand Down Expand Up @@ -132,6 +140,7 @@ class AcqTable(ACACatalogTable):
cand_acqs = MetaAttribute(is_kwarg=False)
p_safe = MetaAttribute(is_kwarg=False)
_fid_set = MetaAttribute(is_kwarg=False, default=())
imposters_mag_limit = MetaAttribute(is_kwarg=False, default=20.0)

@classmethod
def empty(cls):
Expand Down Expand Up @@ -476,6 +485,52 @@ def get_initial_catalog(self):
for name, col in acqs_init.columns.items():
self[name] = col

def calc_p_brightest(self, acq, box_size, man_err=0, bgd=0):
"""
Calculate the probability that the `acq` star is the brightest
candidate in the search box.

This caches the spoiler and imposter stars in the acqs table (the row
corresponding to ``acq``). It is required that the first time this is
called that the box_size and man_err be the maximum, and this is checked.

:param acq: acq stars (AcqTable Row)
:param box_size: box size (float, arcsec)
:param man_err: maneuver error (float, arcsec, default=0)
:param bgd: assume background for imposters (float, e-sec, default=0)

:returns: probability that acq is the brightest (float)
"""
stars = self.stars
dark = self.dark
dither = self.dither

# Spoilers
ext_box_size = box_size + man_err
kwargs = dict(stars=stars, acq=acq, box_size=ext_box_size)
spoilers = get_intruders(acq, ext_box_size, 'spoilers',
n_sigma=2.0, # TO DO: put to characteristics
get_func=get_spoiler_stars, kwargs=kwargs)

# Imposters
ext_box_size = box_size + dither
kwargs = dict(star_row=acq['row'], star_col=acq['col'],
maxmag=acq['mag'] + acq['mag_err'], # + 1.5, TO DO: put to characteristics
box_size=ext_box_size,
dark=dark,
bgd=bgd, # TO DO deal with this
mag_limit=self.imposters_mag_limit
)
imposters = get_intruders(acq, ext_box_size, 'imposters',
n_sigma=1.0, # TO DO: put to characteristics
get_func=get_imposter_stars, kwargs=kwargs)

mags = np.concatenate([spoilers['mag'], imposters['mag']])
mag_errs = np.concatenate([spoilers['mag_err'], imposters['mag_err']])
prob = calc_p_brightest_compare(acq, mags, mag_errs)

return prob

def calc_p_safe(self, verbose=False):
"""
Calculate the probability of a safing action resulting from failure
Expand Down Expand Up @@ -725,7 +780,7 @@ def get_spoiler_stars(stars, acq, box_size):


def get_imposter_stars(dark, star_row, star_col, thresh=None,
maxmag=11.5, box_size=120, bgd=40, test=False):
maxmag=11.5, box_size=120, bgd=40, mag_limit=20.0, test=False):
"""
Note: current alg purposely avoids using the actual flight background
calculation because this is unstable to small fluctuations in values
Expand All @@ -737,9 +792,10 @@ def get_imposter_stars(dark, star_row, star_col, thresh=None,
:param star_row: row of acq star (float)
:param star_col: col of acq star (float)
:param thresh: PEA search hit threshold for a 2x2 block (e-/sec)
:param maxmag: Max mag (alternate way to specify ``thresh``)
:param maxmag: Max mag (alternate way to specify search hit ``thresh``)
:param box_size: box size (arcsec)
:param bgd: assumed flat background (float, e-/sec)
:param mag_limit: Max mag for imposter (using 6x6 readout)
:param test: hook for convenience in algorithm testing

:returns: numpy structured array of imposter stars
Expand Down Expand Up @@ -820,6 +876,9 @@ def get_imposter_stars(dark, star_row, star_col, thresh=None,

img, img_sum, mag, row, col = get_image_props(dark, c_row, c_col, bgd)

if mag > mag_limit:
continue

if pea_reject_image(img):
continue

Expand Down Expand Up @@ -937,52 +996,6 @@ def get_intruders(acq, box_size, name, n_sigma, get_func, kwargs):
return intruders


def calc_p_brightest(acq, box_size, stars, dark, man_err=0,
dither=20, bgd=0):
"""
Calculate the probability that the `acq` star is the brightest
candidate in the search box.

This caches the spoiler and imposter stars in the acqs table (the row
corresponding to ``acq``). It is required that the first time this is
called that the box_size and man_err be the maximum, and this is checked.

:param acq: acq stars (AcqTable Row)
:param box_size: box size (float, arcsec)
:param stars: StarsTable of stars in or near the ACA FOV
:param dark: dark current image (ndarray, e-/sec)
:param man_err: maneuver error (float, arcsec, default=0)
:param dither: dither (float, arsec, default=20)
:param bgd: assume background for imposters (float, e-sec, default=0)

:returns: probability that acq is the brightest (float)
"""
# Spoilers
ext_box_size = box_size + man_err
kwargs = dict(stars=stars, acq=acq, box_size=ext_box_size)
spoilers = get_intruders(acq, ext_box_size, 'spoilers',
n_sigma=2.0, # TO DO: put to characteristics
get_func=get_spoiler_stars, kwargs=kwargs)

# Imposters
ext_box_size = box_size + dither
kwargs = dict(star_row=acq['row'], star_col=acq['col'],
maxmag=acq['mag'] + acq['mag_err'], # + 1.5, TO DO: put to characteristics
box_size=ext_box_size,
dark=dark,
bgd=bgd, # TO DO deal with this
)
imposters = get_intruders(acq, ext_box_size, 'imposters',
n_sigma=1.0, # TO DO: put to characteristics
get_func=get_imposter_stars, kwargs=kwargs)

mags = np.concatenate([spoilers['mag'], imposters['mag']])
mag_errs = np.concatenate([spoilers['mag_err'], imposters['mag_err']])
prob = calc_p_brightest_compare(acq, mags, mag_errs)

return prob


def calc_p_on_ccd(row, col, box_size):
"""
Calculate the probability that star and initial tracked readout box
Expand Down Expand Up @@ -1085,8 +1098,8 @@ def __init__(self, acqs, acq, dither, stars, dark, t_ccd, date):
# man_err, independently because imposter prob is just a
# function of box_size not man_err). Technically also a
# function of dither, but that does not vary here.
p_brightest = calc_p_brightest(acq, box_size=box_size, stars=stars,
dark=dark, man_err=man_err, dither=dither)
p_brightest = acqs.calc_p_brightest(acq, box_size=box_size,
man_err=man_err)
self._p_brightest[box_size, man_err] = p_brightest

# Acquisition probability model value (function of box_size only)
Expand Down
15 changes: 15 additions & 0 deletions proseco/characteristics.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,3 +135,18 @@ def _get_fid_acq_stages():
return out

fid_acq_stages = _get_fid_acq_stages()

# For a given box size, set the mag limit to the magnitude where one expects
# (statistically) N imposters brighter than that limit. The idea is
# that statistically on-board you always have imposters fainter than this limit
# impacting acquisition, and thus they are already accounted for in the
# acquisition probability model. See:
#
# Here the limit is set to 2.5 imposters per 120" box at -10C based on
# the plot in cell 89 of:
#
# http://nbviewer.jupyter.org/url/cxc.cfa.harvard.edu/mta/ASPECT/ipynb/
# star_selection/imposter-stars.ipynb

imposter_mag_lim_ref_mag = 10.0
imposter_mag_lim_ref_t_ccd = -10.0
50 changes: 31 additions & 19 deletions proseco/tests/test_acq.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
from ..report import make_report
from ..acq import (get_p_man_err, bin2x2, CHAR,
get_imposter_stars,
get_image_props, calc_p_brightest,
get_image_props,
AcqTable, calc_p_on_ccd,
get_acq_catalog,
)
Expand All @@ -27,6 +27,18 @@
TEST_COLS = ('idx', 'slot', 'id', 'yang', 'zang', 'halfw')


def calc_p_brightest(acq, box_size, stars, dark, man_err=0, dither=20, bgd=0):
"""
Stub for original functional version of acq.calc_p_brightest, which was
turned into an AcqTable method.
"""
acqs = AcqTable()
acqs.stars = stars
acqs.dark = dark
acqs.dither = dither
return acqs.calc_p_brightest(acq, box_size, man_err, bgd)


def add_imposter(dark, acq, dyang, dzang, dmag):
"""
For testing, add an imposter (single hot pixel) at the specified delta location
Expand Down Expand Up @@ -391,9 +403,9 @@ def test_get_acq_catalog_19387():
' 3 3 37879992 318.47 -1565.92 60',
' 4 4 37882416 481.80 2204.44 60',
' 5 5 37880176 121.33 -1068.25 60',
' 6 -99 37881728 2046.89 1910.79 60',
' 7 6 37880376 -1356.71 1071.32 60',
' 8 7 38276824 -1822.26 -1813.66 60',
' 6 6 37881728 2046.89 1910.79 60',
' 7 7 37880376 -1356.71 1071.32 60',
' 8 -99 38276824 -1822.26 -1813.66 60',
' 9 -99 37882776 1485.00 127.97 60',
' 10 -99 37880152 -1542.43 970.39 60',
' 11 -99 37880584 -2005.80 2449.74 60',
Expand All @@ -412,8 +424,8 @@ def test_get_acq_catalog_19387():
' 3 3 37879992 318.47 -1565.92 60',
' 4 4 37882416 481.80 2204.44 60',
' 5 5 37880176 121.33 -1068.25 60',
' 7 6 37880376 -1356.71 1071.32 60',
' 8 7 38276824 -1822.26 -1813.66 60']
' 6 6 37881728 2046.89 1910.79 60',
' 7 7 37880376 -1356.71 1071.32 60']

assert repr(acqs[TEST_COLS]).splitlines() == exp

Expand All @@ -436,9 +448,9 @@ def test_get_acq_catalog_21007():
' 1 1 189410928 -62.52 1763.04 160',
' 2 2 189409160 -2223.75 1998.69 160',
' 3 3 189417920 1482.94 243.72 160',
' 4 4 189015480 2222.47 -580.99 60',
' 5 5 189417752 1994.07 699.55 100',
' 6 6 189406216 -2311.90 -240.18 80',
' 4 4 189015480 2222.47 -580.99 100',
' 5 5 189417752 1994.07 699.55 60',
' 6 6 189406216 -2311.90 -240.18 120',
' 7 7 189416328 1677.88 137.11 60',
' 8 -99 189416496 333.11 -63.30 120',
' 9 -99 189410280 -495.21 1712.02 120',
Expand All @@ -457,9 +469,9 @@ def test_get_acq_catalog_21007():
' 1 1 189410928 -62.52 1763.04 160',
' 2 2 189409160 -2223.75 1998.69 160',
' 3 3 189417920 1482.94 243.72 160',
' 4 4 189015480 2222.47 -580.99 60',
' 5 5 189417752 1994.07 699.55 100',
' 6 6 189406216 -2311.90 -240.18 60',
' 4 4 189015480 2222.47 -580.99 80',
' 5 5 189417752 1994.07 699.55 60',
' 6 6 189406216 -2311.90 -240.18 120',
' 7 7 189416328 1677.88 137.11 60']

assert repr(acqs[TEST_COLS]).splitlines() == exp
Expand All @@ -481,9 +493,9 @@ def test_box_strategy_20603():
' 3 3 40114416 394.22 1204.43 160',
' 4 4 40112304 -1644.35 2032.47 80',
' 5 5 116923528 -2418.65 1088.40 160',
' 6 6 116791744 985.38 -1210.19 140',
' 7 -99 40108048 2.21 1619.17 120',
' 8 7 116785920 -673.94 -1575.87 80',
' 6 6 116791744 985.38 -1210.19 100',
' 7 7 40108048 2.21 1619.17 100',
' 8 -99 116785920 -673.94 -1575.87 120',
' 9 -99 116923744 -853.18 937.73 120',
' 10 -99 116792320 941.59 -1784.10 120',
' 11 -99 116918232 -2074.91 -1769.96 120',
Expand All @@ -501,8 +513,8 @@ def test_box_strategy_20603():
' 3 3 40114416 394.22 1204.43 140',
' 4 4 40112304 -1644.35 2032.47 160',
' 5 5 116923528 -2418.65 1088.40 160',
' 6 6 116791744 985.38 -1210.19 140',
' 8 7 116785920 -673.94 -1575.87 80']
' 6 6 116791744 985.38 -1210.19 160',
' 7 7 40108048 2.21 1619.17 160']

assert repr(acqs[TEST_COLS]).splitlines() == exp

Expand Down Expand Up @@ -660,8 +672,8 @@ def test_n_acq():
' 2 2 101 0.00 2000.00 160',
' 3 3 109 0.00 1500.00 160',
' 4 4 102 -2000.00 0.00 160',
' 5 -99 110 -1500.00 0.00 120',
' 6 5 103 0.00 -2000.00 140',
' 5 5 110 -1500.00 0.00 160',
' 6 -99 103 0.00 -2000.00 120',
' 7 -99 111 0.00 -1500.00 120',
' 8 -99 104 1000.00 1000.00 120',
' 9 -99 105 1000.00 -1000.00 120',
Expand Down
4 changes: 2 additions & 2 deletions proseco/tests/test_catalog.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,8 @@ def test_get_aca_catalog_20603():
' 5 6 40112304 BOT 6x6 -1644.35 2032.47 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 60',
' 0 10 116785920 ACQ 6x6 -673.94 -1575.87 20 1 80']
' 7 9 116791744 ACQ 6x6 985.38 -1210.19 20 1 160',
' 0 10 40108048 ACQ 6x6 2.21 1619.17 20 1 60']

repr(aca) # Apply default formats
assert aca[TEST_COLS].pformat(max_width=-1) == exp
Expand Down