In [35]:
from __future__ import division

from pprint import pprint

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

In [36]:
filename = 'https://drive.google.com/uc?export=download&id=0B52P_nYJgKHRMnNkaUwtbER6X2c'
acq_stats = Table.read(filename, format='fits')
acq_stats.remove_columns(['tstart', 'tstop', 'slot', 'cat_pos', 'type', 'obi', 'revision', 'mag', 'obc_id'])
acq_stats[:3]


Out[36]:
obsididxagasc_idyangzanghalfwmag_obsyang_obszang_obscolor
114781046540360302.0873.0120.08.625260.65907.0251.29115
114781146537584-477.01784.0120.010.0-518.951818.4750.3842
114781246536768872.01398.0120.09.3125831.21432.1751.11265

In [37]:
filename = 'https://drive.google.com/uc?export=download&id=0B52P_nYJgKHRQkZzWkdXNXRDeDQ'
# filename = 'starcheck_acq_probs.fits.gz'  # local version
acq_probs = Table.read(filename)
acq_probs['acquired'] = np.where(acq_probs['obc_id'] == 'ID', 1, 0)
acq_probs[:2]


Downloading https://drive.google.com/uc?export=download&id=0B52P_nYJgKHRQkZzWkdXNXRDeDQ
|===========================================| 518k/518k (100.00%)         0s
Out[37]:
obsididxdateccd_tempn100magn_warncolor_0p7color_1p5spoilerbad_histacq_probobc_idold_acq_probacquired
1214842011:360:06:30:45.352-17.61797428130.08837398591576.798000000.985ID0.98461
1214852011:360:06:30:45.352-17.61797428130.08837398591578.552000000.985ID0.98461

In [38]:
probs = join(acq_probs, acq_stats, keys=('obsid', 'idx'))
probs[:3]


Out[38]:
obsididxdateccd_tempn100magn_warncolor_0p7color_1p5spoilerbad_histacq_probobc_idold_acq_probacquiredagasc_idyangzanghalfwmag_obsyang_obszang_obscolor
1214842011:360:06:30:45.352-17.61797428130.08837398591576.798000000.985ID0.98461281575096160.0-1989.0120.06.8125154.1-1979.5751.07015
1214852011:360:06:30:45.352-17.61797428130.08837398591578.552000000.985ID0.98461282111880-518.0-1338.0120.08.625-525.025-1329.20.3315
1214862011:360:06:30:45.352-17.61797428130.08837398591576.693000000.985ID0.984612822189522112.01394.0120.06.68752106.1751404.6750.884

In [39]:
len(acq_probs), len(probs)


Out[39]:
(28508, 28476)

In [40]:
fudge_color_1p5 = False
if fudge_color_1p5:
    ok = probs['color_1p5'] == 1
    probs['acq_prob'][ok] -= 0.15

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

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

In [43]:
def calc_acq_probs(probs):
    n_slots = len(probs)
    acq_probs = np.zeros(9, dtype=np.float64)
    for ii in xrange(2 ** n_slots):
        total_prob = 1.0
        n_acq = 0

        for slot in xrange(n_slots):
            if ii & (2 ** slot):
                prob = probs[slot]
                n_acq += 1
            else:
                prob = 1 - probs[slot]
            total_prob *= prob

        acq_probs[n_acq] += total_prob
    return acq_probs

In [44]:
def cum_prob(acq_probs):
    """
    """
    n_expected = 0.0
    cum_probs = np.zeros(9, dtype=np.float64)
    for n_acq in xrange(9):
        n_expected += n_acq * acq_probs[n_acq]
        for n_acq2 in xrange(n_acq, 9):
            cum_probs[n_acq] += acq_probs[n_acq2]

    cum_probs = np.clip(cum_probs, 0, 1)
    cum_probs_fewer = np.clip(1.0 - cum_probs, 1e-15, 1.0)

    return {'n_expected': n_expected,
            'probs_n_or_more': cum_probs,
            'probs_fewer_than_n': cum_probs_fewer,
            'log10_probs_fewer_than_n': np.log10(cum_probs_fewer)}

In [45]:
def cum_prob_by_obsid(obsid):
    pok = pgd[obsid]
    print('Actual stars: {}'.format(np.sum(pok['acquired'])))
    print('CURRENT')
    acq_probs = calc_acq_probs(pok['old_acq_prob'])
    pprint(cum_prob(acq_probs))
    print()
    print('NEW')
    acq_probs = calc_acq_probs(pok['acq_prob'])
    pprint(cum_prob(acq_probs))

In [46]:
def get_n_expected():
    n_acqs = []
    n_expected = []
    obsids = []
    for grp in pg.groups:
        obsids.append(grp['obsid'][0])
        n_acqs.append(np.sum(grp['acquired']))
        acq_probs = calc_acq_probs(grp['acq_prob'])
        out = cum_prob(acq_probs)
        n_expected.append(out['n_expected'])
    return np.array(n_acqs), np.array(n_expected), np.array(obsids)

In [47]:
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 [48]:
def plot_acq_prob_vs_actual(data=probs, verbose=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['acquired'][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.xlabel('Mean predicted acq probability')
    plt.ylabel('Actual acq success rate')

Mean actual success rate vs. predicted within bins


In [49]:
plot_acq_prob_vs_actual();



In [50]:
n_acqs, n_expected, obsids = get_n_expected()

Actual number of stars per obsid vs. expected

The mean predicted number of acq stars does not correlate with the actual.

This demonstrates that there has to be some clustering of properties by star field that impact star acquisition success.

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


In [51]:
plt.plot(n_expected, n_acqs + np.random.uniform(-0.1, 0.1, size=len(n_acqs)), '.', alpha=0.3)
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');


Colors of stars in bad fields

Collect the stars in bad fields (actual acquired < expected - 2) and plot histogram of colors. Do the same for good fields.

Result: No significant difference


In [52]:
bad_mask = n_expected - n_acqs > 2  # These arrays are per-obsid (per-field)

In [53]:
colors = []
for grp in pg.groups[bad_mask].groups:
    colors.append(grp['color'])
bad_colors = np.concatenate(colors)

In [54]:
colors = []
for grp in pg.groups[~bad_mask].groups:
    colors.append(grp['color'])
colors = np.concatenate(colors)
good_colors = colors[~np.isnan(colors)]

In [55]:
plt.hist(good_colors, bins=np.arange(-0.2, 1.6, 0.1), normed=True, alpha=0.5, histtype='stepfilled');
plt.hist(bad_colors, bins=np.arange(-0.2, 1.6, 0.1), normed=True, alpha=0.5, facecolor='r', histtype='stepfilled');


Colors of bad stars

Do the simpler thing and just look at colors of stars that were acquired vs. those not acquired.

Result: No strong difference, except perhaps 40% more color=1.5 stars


In [56]:
id_colors = probs['color'][probs['acquired'].astype(np.bool)]
no_id_colors = probs['color'][~probs['acquired'].astype(np.bool)]

In [57]:
plt.hist(id_colors, bins=np.arange(-0.2, 1.6, 0.1), normed=True, alpha=0.5, histtype='stepfilled');
plt.hist(no_id_colors, bins=np.arange(-0.2, 1.6, 0.1), normed=True, alpha=0.5, facecolor='r', histtype='stepfilled');



In [58]:
print(obsids[bad_mask])


[12911 13363 13516 13518 13994 14033 14955 15378 15756 16311 16413 16425
 16426 52512 52600 52650 52782 53405 53520 53656 53976 54116 54190 54245
 54455 54649 54745 54767 54775 54936]

In [59]:
import webbrowser

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

In [61]:
open_mica(12911)

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


[16311]

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


[14033 15756 52650 54775]

In [77]:
delta = n_expected - n_acqs
isort = np.argsort(delta)[::-1]
isort5 = isort[:5]
delta[isort5]


Out[77]:
array([ 4.30906546,  3.64365345,  3.39872472,  3.30658321,  3.20969516])

In [78]:
n_acqs[isort5]


Out[78]:
array([3, 2, 4, 4, 3])

In [79]:
n_expected[isort5]


Out[79]:
array([ 7.30906546,  5.64365345,  7.39872472,  7.30658321,  6.20969516])

In [81]:
obsids[isort5]


Out[81]:
array([54775, 16311, 54936, 52600, 14033])

In [72]:
open_mica(54775)  # Worst case no apparent star problems

In [73]:
open_mica(16311)  # Warnings, bad histories.  52650 is a repeat attitude with only 3 good acqs

In [74]:
open_mica(54936)  # Nothing unusual, looks clean apart from one spoiler

In [75]:
open_mica(52600)  # Looks clean, warm CCD

In [76]:
open_mica(14033)  # Known issues (bad histories, bad field stars etc)

In [82]:
def get_t_warm(obsid):
    ok = probs['obsid'] == obsid
    row = probs[ok][0]
    return row['ccd_temp'], row['n100']

In [84]:
for obsid in obsids[isort5]:
    print obsid, get_t_warm(obsid)


54775 (-17.141580463699999, 0.096301258532)
16311 (-14.638885349500001, 0.16147280099399999)
54936 (-17.845648749599999, 0.085070262840999999)
52600 (-14.4402241442, 0.165395076625)
14033 (-17.877394094, 0.088915601275700001)

In [ ]: