Fit values here were computed 2018-Mar-08
This version introduces a dependence on the search box size. Search box sizes of 160 or 180 arcsec (required for at least 3 star slots) were used in normal operations starting in the MAR2017 products. This followed two PMSTA anomalies.
In [1]:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from astropy.table import Table
from astropy.time import Time
import tables
from scipy import stats
import tables3_api
%matplotlib inline
In [2]:
SOTA2017_FIT_NO_1P5 = [3.8981465963928441, 5.5208216663935739, 2.0187091292966395, # offsets
-2.2103221221745111, 0.37783433000968347, 0.10035462978065751, # scales
0.0038541777023636111] # brighter than 8.5 mag
In [3]:
SOTA2017_FIT_ONLY_1P5 = [4.4598955940740002, 7.3654868182850661, 4.380944461070051,
-1.4766615762918867, 0.53879889036008366, -0.36463411364645115,
0.0020022525242344045]
In [4]:
SOTA2015_FIT_ALL = [3.9438714542029976, 5.4601129927961134, 1.6582423213669775,
-2.0646518576907495, 0.36414269305801689, -0.0075143036207362852,
0.003740065500207244]
In [5]:
SOTA2015_FIT_NO_1P5 = [4.092016310373646, 6.5415918325159641, 1.8191919043258409,
-2.2301709573082413, 0.30337711472920426, 0.10116735012955963,
0.0043395964215468185]
In [6]:
SOTA2015_FIT_ONLY_1P5 = [4.786710417762472, 4.839392687262392, 1.8646719319052267,
-1.4926740399312248, 0.76412972998935347, -0.20229644263097146,
0.0016270748026844457]
In [7]:
with tables.open_file('/proj/sot/ska/data/acq_stats/acq_stats.h5', 'r') as h5:
cols = h5.root.data.cols
names = {'tstart': 'guide_tstart',
'obsid': 'obsid',
'obc_id': 'acqid',
'halfwidth': 'halfw',
'warm_pix': 'n100_warm_frac',
'mag': 'mag_aca',
'known_bad': 'known_bad',
'color': 'color1',
'img_func': 'img_func',
'ion_rad': 'ion_rad',
'sat_pix': 'sat_pix',
'agasc_id': 'agasc_id',
't_ccd': 'ccd_temp',
'slot': 'slot'}
acqs = Table([getattr(cols, h5_name)[:] for h5_name in names.values()],
names=list(names.keys()))
year_q0 = 1999.0 + 31. / 365.25 # Jan 31 approximately
In [8]:
acqs['year'] = Time(acqs['tstart'], format='cxcsec').decimalyear.astype('f4')
acqs['quarter'] = (np.trunc((acqs['year'] - year_q0) * 4)).astype('f4')
acqs['color_1p5'] = np.where(acqs['color'] == 1.5, 1, 0)
In [9]:
# Filter for year and mag (previously used data through 2007:001)
ok = (acqs['year'] > 2014) & (acqs['mag'] > 6.0) & (acqs['mag'] < 10.6)
In [10]:
# Filter known bad obsids
print('Filtering known bad obsids, start len = {}'.format(np.count_nonzero(ok)))
bad_obsids = [
# Venus
2411,2414,6395,7306,7307,7308,7309,7311,7312,7313,7314,7315,7317,7318,7406,583,
7310,9741,9742,9743,9744,9745,9746,9747,9749,9752,9753,9748,7316,15292,16499,
16500,16501,16503,16504,16505,16506,16502,
]
for badid in bad_obsids:
ok = ok & (acqs['obsid'] != badid)
print('Filtering known bad obsids, end len = {}'.format(np.count_nonzero(ok)))
In [11]:
data_all = acqs[ok]
data_all.sort('year')
data_all['mag10'] = data_all['mag'] - 10.0
In [12]:
# Adjust probability (in probit space) for box size. See:
# https://github.com/sot/skanb/blob/master/pea-test-set/fit_box_size_acq_prob.ipynb
b1 = 0.96
b2 = -0.30
box0 = (data_all['halfwidth'] - 120) / 120 # normalized version of box, equal to 0.0 at nominal default
data_all['box_delta'] = b1 * box0 + b2 * box0**2
In [13]:
# Create 'fail' column, rewriting history as if the OBC always
# ignore the MS flag in ID'ing acq stars. Set ms_disabled = False
# to not do this
obc_id = data_all['obc_id']
obc_id_no_ms = (data_all['img_func'] == 'star') & ~data_all['sat_pix'] & ~data_all['ion_rad']
data_all['fail'] = np.where(obc_id | obc_id_no_ms, 0.0, 1.0)
In [14]:
data_all = data_all.group_by('quarter')
data_mean = data_all.groups.aggregate(np.mean)
In [15]:
def p_fail(pars, m10, wp, box_delta=0.0):
"""
Acquisition probability model
:param pars: 7 parameters (3 x offset, 3 x scale, p_fail for bright stars)
:param m10: mag - 10
:param wp: warm pixel fraction
:param box: search box half width (arcsec)
"""
scl0, scl1, scl2 = pars[0:3]
off0, off1, off2 = pars[3:6]
p_bright_fail = pars[6]
# Make sure inputs have same dimensions
m10, wp, box_delta = np.broadcast_arrays(m10, wp, box_delta)
scale = scl0 + scl1 * m10 + scl2 * m10**2
offset = off0 + off1 * m10 + off2 * m10**2
p_fail = offset + scale * wp + box_delta
p_fail = stats.norm.cdf(p_fail) # probit transform
p_fail[m10 < -1.5] = p_bright_fail # For stars brighter than 8.5 mag use a constant
return p_fail
def p_acq_fail(data=None):
"""
Sherpa fit function wrapper to ensure proper use of data in fitting.
"""
if data is None:
data = data_all
m10 = data['mag10']
wp = data['warm_pix']
box_delta = data['box_delta']
def sherpa_func(pars, x):
return p_fail(pars, m10, wp, box_delta)
return sherpa_func
In [16]:
def fit_sota_model(data_mask=None):
from sherpa import ui
data = data_all if data_mask is None else data_all[data_mask]
data_id = 1
ui.set_method('simplex')
ui.set_stat('cash')
ui.load_user_model(p_acq_fail(data), 'model')
ui.add_user_pars('model', ['scl0', 'scl1', 'scl2', 'off0', 'off1', 'off2', 'p_bright_fail'])
ui.set_model(data_id, 'model')
ui.load_arrays(data_id, np.array(data['year']), np.array(data['fail'], dtype=np.float))
# Initial fit values from fit of all data
start_vals = iter(SOTA2015_FIT_ALL) # Offset
fmod = ui.get_model_component('model')
for name in ('scl', 'off'):
for num in (0, 1, 2):
comp_name = name + str(num)
setattr(fmod, comp_name, next(start_vals))
comp = getattr(fmod, comp_name)
comp.min = -100000
comp.max = 100000
# ui.freeze(comp)
fmod.p_bright_fail = 0.025
fmod.p_bright_fail.min = 0.0
fmod.p_bright_fail.max = 1.0
# ui.freeze(fmod.p_bright_fail)
ui.fit(data_id)
# conf = ui.get_confidence_results()
return ui.get_fit_results()
In [17]:
def plot_fit_grouped(pars, group_col, group_bin, mask=None, log=False, colors='br', label=None, probit=False):
data = data_all if mask is None else data_all[mask]
data['model'] = p_acq_fail(data)(pars, None)
group = np.trunc(data[group_col] / group_bin)
data = data.group_by(group)
data_mean = data.groups.aggregate(np.mean)
len_groups = np.diff(data.groups.indices)
data_fail = data_mean['fail']
model_fail = data_mean['model']
# Possibly plot the data and model probabilities in probit space
if probit:
data_fail = stats.norm.ppf(data_fail)
model_fail = stats.norm.ppf(model_fail)
fail_sigmas = np.sqrt(data_fail * len_groups) / len_groups
plt.errorbar(data_mean[group_col], data_fail, yerr=fail_sigmas, fmt='.' + colors[0], label=label)
plt.plot(data_mean[group_col], model_fail, '-' + colors[1])
if log:
ax = plt.gca()
ax.set_yscale('log')
In [18]:
def mag_filter(mag0, mag1):
ok = (data_all['mag'] > mag0) & (data_all['mag'] < mag1)
return ok
In [19]:
def wp_filter(wp0, wp1):
ok = (data_all['warm_pix'] > wp0) & (data_all['warm_pix'] < wp1)
return ok
In [20]:
def plot_fit_all(parvals, mask=None, probit=False):
if mask is None:
mask = np.ones(len(data_all), dtype=bool)
plt.figure()
plot_fit_grouped(parvals, 'mag', 0.25, wp_filter(0.20, 0.25) & mask, log=False,
colors='gk', label='0.20 < WP < 0.25')
plot_fit_grouped(parvals, 'mag', 0.25, wp_filter(0.10, 0.20) & mask, log=False,
colors='cm', label='0.10 < WP < 0.20')
plt.legend(loc='upper left');
plt.ylim(0.001, 1.0);
plt.xlim(9, 11)
plt.grid()
plt.figure()
plot_fit_grouped(parvals, 'mag', 0.25, wp_filter(0.20, 0.25) & mask, probit=True, colors='gk', label='0.20 < WP < 0.25')
plot_fit_grouped(parvals, 'mag', 0.25, wp_filter(0.10, 0.20) & mask, probit=True, colors='cm', label='0.10 < WP < 0.20')
plt.legend(loc='upper left');
# plt.ylim(0.001, 1.0);
plt.xlim(9, 11)
plt.grid()
plt.figure()
plot_fit_grouped(parvals, 'warm_pix', 0.02, mag_filter(10.3, 10.6) & mask, log=False, colors='gk', label='10.3 < mag < 10.6')
plot_fit_grouped(parvals, 'warm_pix', 0.02, mag_filter(10, 10.3) & mask, log=False, colors='cm', label='10 < mag < 10.3')
plot_fit_grouped(parvals, 'warm_pix', 0.02, mag_filter(9, 10) & mask, log=False, colors='br', label='9 < mag < 10')
plt.legend(loc='best')
plt.grid()
plt.figure()
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(10.3, 10.6) & mask, colors='gk', label='10.3 < mag < 10.6')
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(10, 10.3) & mask, colors='cm', label='10 < mag < 10.3')
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(9.5, 10) & mask, colors='br', label='9.5 < mag < 10')
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(9.0, 9.5) & mask, colors='gk', label='9.0 < mag < 9.5')
plt.legend(loc='best')
plt.grid()
plt.figure()
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(10.3, 10.6) & mask, colors='gk', label='10.3 < mag < 10.6', probit=True)
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(10, 10.3) & mask, colors='cm', label='10 < mag < 10.3', probit=True)
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(9.5, 10) & mask, colors='br', label='9.5 < mag < 10', probit=True)
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(9.0, 9.5) & mask, colors='gk', label='9.0 < mag < 9.5', probit=True)
plt.legend(loc='best')
plt.grid();
In [21]:
# fit = fit_sota_model(data_all['color'] == 1.5, ms_disabled=True)
mask_no_1p5 = data_all['color'] != 1.5
In [22]:
print('Hang tight, this could take a few minutes')
fit_no_1p5 = fit_sota_model(mask_no_1p5)
In [23]:
plot_fit_all(fit_no_1p5.parvals, mask=mask_no_1p5)
In [24]:
plot_fit_all(SOTA2017_FIT_NO_1P5, mask=mask_no_1p5)
In [25]:
print('Hang tight, this could take a few minutes')
mask_1p5 = data_all['color'] == 1.5
fit_1p5 = fit_sota_model(mask_1p5)
In [26]:
plot_fit_all(fit_1p5.parvals, mask=mask_1p5)
In [27]:
mag = np.linspace(9, 11, 30)
for wp in (0.1, 0.2, 0.3):
plt.plot(mag, p_fail(fit_no_1p5.parvals, mag-10, wp), 'r',
label='2015 model' if wp == 0.1 else None)
plt.plot(mag, p_fail(SOTA2017_FIT_NO_1P5, mag-10, wp), 'b',
label='2017 model' if wp == 0.1 else None)
plt.grid()
plt.xlabel('Mag')
plt.ylim(0, 1)
plt.title('Failure prob vs. mag for Wp=(0.1, 0.2, 0.3)')
plt.legend(loc='upper left')
plt.ylabel('Prob');
In [28]:
mag = np.linspace(9, 11, 30)
for wp in (0.1, 0.2, 0.3):
plt.plot(mag, stats.norm.ppf(p_fail(fit_no_1p5.parvals, mag-10, wp)), 'r',
label='2018 model' if wp == 0.1 else None)
plt.plot(mag, stats.norm.ppf(p_fail(SOTA2017_FIT_NO_1P5, mag-10, wp)), 'b',
label='2017 model' if wp == 0.1 else None)
plt.grid()
plt.xlabel('Mag')
# plt.ylim(0, 1)
plt.title('Failure prob vs. mag for Wp=(0.1, 0.2, 0.3)')
plt.legend(loc='upper left')
plt.ylabel('Prob');
In [29]:
for mag in (10.0, 10.25, 10.5):
wp = np.linspace(0, 0.4, 30)
plt.plot(wp, p_fail(fit_no_1p5.parvals, mag-10, wp), 'r',
label='2015 model' if mag == 10.0 else None)
plt.plot(wp, p_fail(SOTA2017_FIT_NO_1P5, mag-10, wp), 'b',
label='2017 model' if mag == 10.0 else None)
plt.grid()
plt.xlabel('Warm pix frac')
plt.ylim(0, 1)
plt.title('Failure prob vs. Wp for mag=(10.0, 10.25, 10.5)')
plt.ylabel('Fail prob');
In [30]:
plt.hist(data_all['warm_pix'], bins=100)
plt.grid()
plt.xlabel('Warm pixel fraction');
In [31]:
plt.hist(data_all['mag'], bins=np.arange(6, 11.1, 0.1))
plt.grid()
plt.xlabel('Mag_aca')
Out[31]:
In [32]:
plt.plot(data_all['year'], data_all['warm_pix'])
plt.ylim(0, None)
plt.grid();
In [33]:
ok = (data_all['mag'] > 10.3) & (data_all['mag'] < 10.6)
da = data_all[ok]
In [34]:
from Chandra.Time import DateTime
yr0 = DateTime('2017:079').frac_year
yr1 = DateTime('2017-10-01T00:00:00').frac_year
yr2 = DateTime('2018:060').frac_year
In [ ]:
In [35]:
ok1 = (yr0 < da['year']) & (da['year'] < yr1)
ok2 = (yr1 < da['year']) & (da['year'] < yr2)
plt.figure()
plt.hist(da['mag'][ok1], facecolor='C0')
plt.hist(da['mag'][ok2], facecolor='C1', alpha=0.5)
Out[35]:
In [36]:
print(np.mean(da['fail'][ok1]))
print(np.mean(da['fail'][ok2]))
In [37]:
print(np.count_nonzero(ok1))
print(np.count_nonzero(ok2))
In [38]:
print(np.mean(da['mag'][ok1]))
print(np.mean(da['mag'][ok2]))
In [39]:
print(len(set(da['agasc_id'][ok1])))
print(len(set(da['agasc_id'][ok2])))
In [40]:
print(np.mean(da['warm_pix'][ok1]))
print(np.mean(da['warm_pix'][ok2]))
In [41]:
print(np.mean(da['t_ccd'][ok1]))
print(np.mean(da['t_ccd'][ok2]))
In [42]:
from scipy.stats import binom
In [43]:
n = np.count_nonzero(ok1)
k = np.count_nonzero(da['fail'][ok1])
print(binom.ppf(0.01, n, k/n) / n)
print(binom.ppf(0.99, n, k/n) / n)
In [44]:
n = np.count_nonzero(ok2)
k = np.count_nonzero(da['fail'][ok2])
print(binom.ppf(0.01, n, k/n) / n)
print(binom.ppf(0.99, n, k/n) / n)
In [45]:
print(np.mean(da['halfwidth'][ok1]))
print(np.mean(da['halfwidth'][ok2]))
In [46]:
print(np.mean(da['ion_rad'][ok1]))
print(np.mean(da['ion_rad'][ok2]))
In [47]:
print(np.mean(da['sat_pix'][ok1]))
print(np.mean(da['sat_pix'][ok2]))
In [48]:
plot_fit_grouped(fit_no_1p5.parvals, 'year', 0.10, mag_filter(10.3, 10.6) & mask_no_1p5,
colors='gk', label='10.3 < mag < 10.6')
plt.xlim(2016.0, None)
y0, y1 = plt.ylim()
x = DateTime('2017-10-01T00:00:00').frac_year
plt.plot([x, x], [y0, y1], '--r', alpha=0.5)
plt.grid();
In [49]:
plot_fit_grouped(fit_no_1p5.parvals, 'year', 0.10, mag_filter(10.0, 10.3) & mask_no_1p5,
colors='gk', label='10.0 < mag < 10.3')
plt.xlim(2016.0, None)
y0, y1 = plt.ylim()
x = DateTime('2017-10-01T00:00:00').frac_year
plt.plot([x, x], [y0, y1], '--r', alpha=0.5)
plt.grid();
In [ ]: