Acq stats: comparing model to data and inferring results with MS filtering disabled


In [1]:
from __future__ import division

import sys

from pprint import pprint
from astropy.time import Time
from itertools import izip

import matplotlib.pyplot as plt
import tables
import numpy as np
from astropy.table import Table, join
%matplotlib inline

In [2]:
sys.path.insert(0, '/home/aldcroft/git/starcheck')
from starcheck import star_probs

In [3]:
tstart = Time('2014:001').cxcsec
filename = '/proj/sot/ska/data/acq_stats/acq_stats.h5'
with tables.openFile(filename) as h5:
    acq_stats = h5.root.data.readWhere('guide_tstart > {}'.format(tstart))
    acq_stats = Table(acq_stats)
acq_stats[:3]


Out[3]:
<Table masked=False length=3>
obsidobiacq_startguide_startguide_tstartone_shot_lengthrevisionslotidxtypeyangzanghalfwmagacqidstar_trackedspoiler_trackedimg_funcn_trak_intervmax_trak_cdymin_trak_cdymean_trak_cdymax_trak_cdzmin_trak_cdzmean_trak_cdzmax_trak_magmin_trak_magmean_trak_magcdycdzdydzion_raddef_pixmult_starsat_pixmag_obsyang_obszang_obsagasc_idcolor1radecepochpm_rapm_decvarpos_errmag_acamag_errmag_bandpos_catidaspq1aspq2aspq3acqq1acqq2acqq4n100_warm_fracccd_tempknown_badbad_comment
int64int64string168string168float64float64string120int64int64string40float64float64int64float64boolboolboolstring56int64float64float64float64float64float64float64float64float64float64float64float64float64float64boolboolboolboolfloat64float64float64int64float64float64float64float64int64int64int64int64float64int64int64int64int64int64int64int64int64int64float64float64boolstring120
5312102014:001:01:18:28.6852014:001:01:19:09.685504926416.86914.40329477311.001BOT-350.0-92.01208.057TrueTrueFalsestar10.8622305520160.7781506939510.8278042839020.5881069399090.5148868006110.5534652859688.1258.1258.1250.7781506939510.578296045783-6.38303387603-15.0257216855FalseFalseFalseFalse8.125-364.4-100.554451215040.0934997051954109.3238663748.464041282000.0-13-38-999978.0572748184214500999-99995863730.132767511564-16.0350036621False
5312102014:001:01:18:28.6852014:001:01:19:09.685504926416.86914.40329477311.012BOT-1490.0842.01208.069TrueTrueFalsestar10.5598019761420.4846866851040.5232331586920.1378519856050.09692627837260.1167202964728.06258.06258.06250.4846866851040.137851985605-6.67684443381-15.4669500247FalseFalseFalseFalse8.0625-1505.025833.954451216641.17810058594109.6657160248.123243112000.069-9999138.06908321381145009995284602710.132767511564-16.0350036621False
5312102014:001:01:18:28.6852014:001:01:19:09.685504926416.86914.40329477311.023BOT-619.01666.01209.261TrueTrueFalsestar10.163779796396-0.01214904779230.08508371629320.2765050033780.1355897341850.220691756559.18759.18759.18750.03755862298690.152129139578-7.12477668131-15.4519421955FalseFalseFalseFalse9.1875-634.11658.74451267201.11520028114109.2197940747.974353482000.0-9999-14-9999259.2614345550524500999-99992362360.132767511564-16.0350036621False

In [4]:
acq_probs = star_probs.acq_success_prob(date=acq_stats['guide_tstart'],
                                    t_ccd=acq_stats['ccd_temp'], 
                                    mag=acq_stats['mag_aca'], 
                                    color=acq_stats['color1'])
acq_stats['acq_prob'] = acq_probs

In [5]:
def acqid_colname(without_ms):
    return 'acqid_no_MS' if without_ms else 'acqid'

In [6]:
probs = acq_stats.copy()  # for compatibility with previous convention in notebook

In [7]:
probs[acqid_colname(True)] = probs['acqid'] | ((probs['img_func'] == 'star') 
                                               & ~probs['ion_rad'] & ~probs['sat_pix'])

In [8]:
# Group by obsid for later
pg = probs.group_by('obsid')

In [9]:
def get_probs_groups_dict():
    return {grp['obsid'][0].tolist(): grp for grp in pg.groups}

In [10]:
def get_n_expected(without_ms=False):
    obsids = pg['obsid'].groups.aggregate(lambda x: x[0])
    n_expected = pg['acq_prob'].groups.aggregate(np.sum)
    n_acqs = pg[acqid_colname(without_ms)].groups.aggregate(np.sum)
    return n_acqs, n_expected, obsids

In [11]:
def get_n_or_fewer_expected(n_or_fewer, without_ms=False):
    lte_probs = []
    for i0, i1 in izip(pg.groups.indices[:-1], pg.groups.indices[1:]):
        n_acq_probs, n_or_fewer_probs = star_probs.prob_n_acq(pg['acq_prob'][i0:i1])
        lte_probs.append(n_or_fewer_probs[n_or_fewer])
    return np.array(lte_probs)

In [12]:
def plot_n_expected(n_acqs, n_expected):
    plt.clf()
    plt.plot(n_expected, n_acqs + np.random.uniform(-0.3, 0.3, size=len(n_acqs)), '.', alpha=0.3)
    plt.plot([0,8], [0,8], '--m')
    plt.ylim(-0.1, 8.5)
    plt.xlim(-0.1, 8.5)
    plt.grid()

In [13]:
def plot_acq_prob_vs_actual(data=probs, verbose=False, without_ms=False):
    bins = np.array([0, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0])
    p_actuals = []
    p_acq_probs = []
    for p0, p1 in zip(bins[:-1], bins[1:]):
        ok = (data['acq_prob'] >= p0) & (data['acq_prob'] < p1)
        n = np.sum(ok)
        if n > 10:
            pok = data[acqid_colname(without_ms)][ok]
            p_actuals.append(np.sum(pok) / n)
            pok = data['acq_prob'][ok]
            p_acq_probs.append(np.mean(pok))
            if verbose:
                print(p0, p1, len(pok))
    plt.clf()
    plt.plot(p_acq_probs, p_actuals, '.-')
    plt.grid()
    plt.plot([0, 1], [0, 1], '--r')
    plt.title('Actual vs predicted ACQ success{}'
              .format(' (without MS)' if without_ms else ''))
    plt.xlabel('Mean predicted acq probability')
    plt.ylabel('Actual acq success rate')

Mean actual success rate vs. predicted within bins


In [14]:
plot_acq_prob_vs_actual();



In [15]:
plot_acq_prob_vs_actual(without_ms=True);


Actual number of stars per obsid vs. expected

The mean predicted number of acq stars correlates with the actual when MS filtering is removed.

This demonstrates that there has to be some clustering of properties by star field that impact star acquisition success. That property is now known to be Multiple Stars Suspected (MS) filtering.

Label as "bad fields" the obsids where the actual was at least 2 stars less than expected.


In [16]:
def plot_n_expected_per_obsid(without_ms=False):
    n_acqs, n_expected, obsids = get_n_expected(without_ms)
    plt.plot(n_expected, n_acqs + np.random.uniform(-0.1, 0.1, size=len(n_acqs)), '.', alpha=0.5)
    plt.xlim(0, 8.5)
    plt.ylim(0, 8.5)
    plt.grid()
    plt.plot([0,8.5], [0, 8.5], '--r')
    plt.plot([2,8.5], [0, 6.5], '--m')
    plt.xlabel('Predicted acq stars')
    plt.ylabel('Actual acq stars')
    plt.title('Acq stars per obsid');

In [17]:
plot_n_expected_per_obsid()



In [18]:
plot_n_expected_per_obsid(without_ms=True)


Why mean expected number of stars is not ideal

Here we consider two extreme cases that both have 4.0 expected stars.

In [15]:
def get_acq_probs_table(probs):
    n_acq_probs, n_or_fewer_probs = star_probs.prob_n_acq(probs)
    dat = Table([np.arange(9), n_acq_probs * 100, n_or_fewer_probs * 100], 
                names=['n_stars', 'P(n)', 'P(n_or_fewer)'])
    dat['P(n)'].format = '{0:.1f}'
    dat['P(n_or_fewer)'].format = '{0:.1f}'
    return dat

Eight faint stars with P_acq = 0.5


In [16]:
get_acq_probs_table([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5])


Out[16]:
<Table masked=False length=9>
n_starsP(n)P(n_or_fewer)
int64float64float64
00.40.4
13.13.5
210.914.5
321.936.3
427.363.7
521.985.5
610.996.5
73.199.6
80.4100.0

Four bright stars and four very faint stars


In [17]:
print_acq_probs_table([0.985, 0.985, 0.985, 0.985, 0.015, 0.015, 0.015, 0.015])


Out[17]:
<Table masked=False length=9>
n_starsP(n)P(n_or_fewer)
int64float64float64
00.00.0
10.00.0
20.10.1
35.45.5
488.994.5
55.499.9
60.1100.0
70.0100.0
80.0100.0

Actual number of stars per obsid vs. P(2 or fewer)


In [19]:
def plot_lte_per_obsid(n_lte, without_ms=False):
    n_acqs, n_expected, obsids = get_n_expected(without_ms)
    lte_probs = get_n_or_fewer_expected(n_lte, without_ms)
    plt.semilogx(lte_probs, 
             n_acqs + np.random.uniform(-0.1, 0.1, size=len(n_acqs)), 
             '.', alpha=0.5)
    # plt.xlim(0, 1)
    plt.ylim(0, 8.5)
    plt.grid()
    plt.xlabel('Probability for {} or fewer stars'.format(n_lte))
    plt.ylabel('Actual acq stars')
    plt.title('P({} or fewer stars) vs. N actual acq stars'
             .format(n_lte));
    print('Expected # of {} or fewer: {:.2f}'.format(n_lte, np.sum(lte_probs)))
    print('Actual # of {} or fewer: {:d}'.format(n_lte, np.count_nonzero(n_acqs <= n_lte)))

In [20]:
plot_lte_per_obsid(2, without_ms=False)


Expected # of 2 or fewer: 0.34
Actual # of 2 or fewer: 1

In [21]:
plot_lte_per_obsid(3, without_ms=False)


Expected # of 3 or fewer: 3.67
Actual # of 3 or fewer: 5

In [22]:
plot_lte_per_obsid(2, without_ms=True)


Expected # of 2 or fewer: 0.34
Actual # of 2 or fewer: 0

In [23]:
plot_lte_per_obsid(1, without_ms=False)


Expected # of 1 or fewer: 0.02
Actual # of 1 or fewer: 0

In [24]:
lte_2_probs = get_n_or_fewer_expected(2, without_ms=False)
n_acqs, n_expected, obsids = get_n_expected(without_ms=False)

In [25]:
ok = lte_2_probs > 0.01
Table([obsids[ok], n_acqs[ok], n_expected[ok], lte_2_probs[ok]],
      names=['obsid', 'n_acq_actual', 'n_expected', 'P(lte2)'])


Out[25]:
<Table masked=False length=5>
obsidn_acq_actualn_expectedP(lte2)
int64int64float64float64
1566254.548407547910.0240487472506
1611654.936089098190.0212924052593
1627954.175505311870.0386672047219
1699854.690982141230.0190188543562
1763154.651372035450.0423979520903

In [26]:
import webbrowser

In [27]:
def open_mica(obsid):
    webbrowser.open_new_tab('http://kadi.cfa.harvard.edu/mica/?obsid_or_date={}'.format(obsid))

In [28]:
open_mica(15662)

In [29]:
open_mica(17631)

In [30]:
twos = n_acqs == 2
print(obsids[twos])


obsid
-----
16311

In [31]:
pg.groups[twos]


Out[31]:
<Table masked=False length=8>
obsidobiacq_startguide_startguide_tstartone_shot_lengthrevisionslotidxtypeyangzanghalfwmagacqidstar_trackedspoiler_trackedimg_funcn_trak_intervmax_trak_cdymin_trak_cdymean_trak_cdymax_trak_cdzmin_trak_cdzmean_trak_cdzmax_trak_magmin_trak_magmean_trak_magcdycdzdydzion_raddef_pixmult_starsat_pixmag_obsyang_obszang_obsagasc_idcolor1radecepochpm_rapm_decvarpos_errmag_acamag_errmag_bandpos_catidaspq1aspq2aspq3acqq1acqq2acqq4n100_warm_fracccd_tempknown_badbad_commentacq_probacqid_no_MS
int64int64string168string168float64float64string120int64int64string40float64float64int64float64boolboolboolstring56int64float64float64float64float64float64float64float64float64float64float64float64float64float64boolboolboolboolfloat64float64float64int64float64float64float64float64int64int64int64int64float64int64int64int64int64int64int64int64int64int64float64float64boolstring120float64bool
1631112014:189:23:01:18.7602014:189:23:01:56.685521247783.8698.659720140661.0011ACQ-1663.0-1185.012010.154FalseTrueFalsestar10.0976189206281-0.01534555167150.02734271394140.1633690561170.04092718670340.10532741312510.12510.12510.1250.09731364663620.040927186703410.3021134882-7.59605033039FalseFalseTrueFalse10.125-1661.0-1185.2753312228961.47729992867184.4081978130.378027182000.0-2217-99994710.153566360574500999-9999-9999-99990.160091439733-14.7281799316False0.695861079607True
1631112014:189:23:01:18.7602014:189:23:01:56.685521247783.8698.659720140661.0112ACQ29.01703.012010.24FalseTrueFalsestar1-1.86381106038-2.26874320965-2.083570757090.5675757021730.3140413104860.41747496416410.510.187510.3515625-2.180391701040.4200603542178.02342523425-7.21660013528FalseFalseTrueFalse10.312529.451703.8252609650481.49175012112185.1206673929.682199282000.0-9999-18-99994810.240335464564500999-9999-9999-99990.160091439733-14.7281799316False0.626839632142True
1631112014:189:23:01:18.7602014:189:23:01:56.685521247783.8698.659720140661.0213ACQ-1849.0675.012010.262FalseTrueFalsestar11.589461526591.588914365741.589187946171.534533114031.517354572761.525943843399.93759.93759.93751.588914365741.5173545727611.793438569-6.12048311992FalseFalseTrueFalse9.9375-1845.25677.353313545441.05995059013184.9960506130.267009732000.0-9999-9999-99994210.2624740601545004781271271270.160091439733-14.7281799316False0.608061764417True
1631112014:189:23:01:18.7602014:189:23:01:56.685521247783.8698.659720140661.034BOT1992.01284.012010.258FalseTrueFalsestar1-0.579422132961-0.877140580925-0.722012485002-0.840992309763-1.20343737046-1.0123486207310.187510.010.1136363636-0.682325881264-1.03906033059.52123109695-8.67421350726FalseFalseTrueFalse10.1251993.91282.752609654160.419050037861184.7984221429.200373312000.0-4413-99992610.257882118234500999-9999-9999-99990.160091439733-14.7281799316False0.611991700597True
1631112014:189:23:01:18.7602014:189:23:01:56.685521247783.8698.659720140661.045BOT489.02016.012010.253FalseTrueFalsestar1-1.27671028545-1.56147330637-1.376193759660.256143801436-0.1274227633040.05348058891110.37510.2510.34375-1.561473306370.05016093508228.6420315682-7.58606336112FalseFalseTrueFalse10.375489.7752016.0752609660480.612849771976185.1696238729.5335812000.0-3018-99993110.252649307344500999-99992532530.160091439733-14.7281799316False0.61644807734True
1631112014:189:23:01:18.7602014:189:23:01:56.685521247783.8698.659720140661.056BOT-1434.0-379.01206.46TrueTrueFalsestar10.05582818614790.009488984769180.02741712215420.2914801929440.2228432578940.2474049211476.43756.43756.43750.02920763861180.23507565733510.233893536-7.40213947118FalseFalseFalseFalse6.4375-1432.075-379.03313682400.271149814129184.6317571730.249046052000.077-119-999966.4603681564314500999-9999-9999-99990.160091439733-14.7281799316False0.985True
1631112014:189:23:01:18.7602014:189:23:01:56.685521247783.8698.659720140661.069ACQ-230.0872.012010.055FalseTrueFalsestar11.985252301511.819121178541.91900822892-3.08439794902-3.19389938464-3.12548887110.062510.010.02678571431.94302141046-3.1026982353512.1472649673-10.7395378339FalseFalseTrueFalse10.0-226.025869.5752609678481.5184.8941031329.822740972000.0-9999-12-99994510.055201530564500999-9999-99994060.160091439733-14.7281799316False0.458421956971True
1631112014:189:23:01:18.7602014:189:23:01:56.685521247783.8698.659720140661.0710ACQ-1397.0443.01208.709TrueTrueFalsestar10.1605083209830.04623429522990.116650008012-0.0508174201224-0.116586432105-0.08189825384428.68758.6258.6406250.159281581753-0.067936944161610.3637988547-7.70541481293FalseFalseFalseFalse8.625-1395.0442.9753313505440.866150021553184.879941230.167953232000.0-54-9999-9999158.708631515524500999-9999-9999-99990.160091439733-14.7281799316False0.985True

In [32]:
threes = n_acqs == 3
print(obsids[threes])


obsid
-----
15756
16620
17322
17632

In [33]:
low_n_exp = n_expected < 4.5
Table([obsids[low_n_exp], n_expected[low_n_exp], n_acqs[low_n_exp]],
      names=['obsid', 'n_expected', 'n_acq_actual'])


Out[33]:
<Table masked=False length=4>
obsidn_expectedn_acq_actual
int64float64int64
161054.473231835785
162794.175505311875
163744.486593916454
173784.159997960975

In [ ]: