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]:
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]
Out[37]:
In [38]:
probs = join(acq_probs, acq_stats, keys=('obsid', 'idx'))
probs[:3]
Out[38]:
In [39]:
len(acq_probs), len(probs)
Out[39]:
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')
In [49]:
plot_acq_prob_vs_actual();
In [50]:
n_acqs, n_expected, obsids = get_n_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');
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');
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])
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])
In [64]:
threes = n_acqs == 3
print(obsids[threes])
In [77]:
delta = n_expected - n_acqs
isort = np.argsort(delta)[::-1]
isort5 = isort[:5]
delta[isort5]
Out[77]:
In [78]:
n_acqs[isort5]
Out[78]:
In [79]:
n_expected[isort5]
Out[79]:
In [81]:
obsids[isort5]
Out[81]:
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)
In [ ]: