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]:
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')
In [14]:
plot_acq_prob_vs_actual();
In [15]:
plot_acq_prob_vs_actual(without_ms=True);
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)
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
In [16]:
get_acq_probs_table([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5])
Out[16]:
In [17]:
print_acq_probs_table([0.985, 0.985, 0.985, 0.985, 0.015, 0.015, 0.015, 0.015])
Out[17]:
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)
In [21]:
plot_lte_per_obsid(3, without_ms=False)
In [22]:
plot_lte_per_obsid(2, without_ms=True)
In [23]:
plot_lte_per_obsid(1, without_ms=False)
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]:
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])
In [31]:
pg.groups[twos]
Out[31]:
In [32]:
threes = n_acqs == 3
print(obsids[threes])
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]:
In [ ]: