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
%matplotlib inline
In [2]:
SOTA2015_FIT_ALL = [9.0122372262846735,
6.5882697680183346,
-0.77727342693609758,
-2.1477049943999713,
0.67132336969927486,
0.33209623740442562]
In [3]:
SOTA2015_FIT_NO_1P5 = [9.6887121629584954,
9.1613039805725531,
-0.41919346946809299,
-2.3829996969445664,
0.54998935209070199,
0.47839261003413558,
0.0060869131558326916]
In [4]:
SOTA2015_FIT_ONLY_1P5 = [8.5417092977479818,
0.4448269654965103,
-3.5137851315167579,
-1.350542439869622,
1.527806122455764,
0.30973568451354011,
0.005414847154007865]
In [5]:
with tables.openFile('/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',
'warm_pix': 'n100_warm_frac',
'mag': 'mag_aca',
'known_bad': 'known_bad',
'color': 'color1'}
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
acqs['year'] = Time(acqs['tstart'], format='cxcsec').decimalyear.astype('f4')
acqs['quarter'] = (np.trunc((acqs['year'] - year_q0) * 4)).astype('f4')
acqs['fail'] = 1.0 - acqs['obc_id'].astype(np.int)
acqs['color_1p5'] = np.where(acqs['color'] == 1.5, 1, 0)
In [6]:
data_all = acqs.group_by('quarter')
data_all['mag10'] = data_all['mag'] - 10.0
data_all.sort('year')
ok = (data_all['year'] > 2007) & (data_all['mag'] > 6.0) & (data_all['mag'] < 11.0)
data_all = data_all[ok]
data_all = data_all.group_by('quarter')
data_mean = data_all.groups.aggregate(np.mean)
ok = np.ones(len(data_all), dtype=bool)
print('Filtering known bad obsids, start len = {}'.format(len(data_all)))
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 & (data_all['obsid'] != badid)
data_all = data_all[ok]
print('Filtering known bad obsids, end len = {}'.format(len(data_all)))
In [7]:
# SET DATA FILTERING for this run (TO DO: refactor into a function)
data_all = data_all[data_all['color'] != 1.5]
data_all
Out[7]:
In [8]:
def p_acq_fail(pars, x, data=None, mag=None, wp=None):
if data is None:
data = data_all
m10 = data['mag10'] if mag is None else mag - 10.0
if wp is None:
wp = data['warm_pix']
scl0, scl1, scl2 = pars[0:3]
off0, off1, off2 = pars[3:6]
p_bright_fail = pars[6]
scale = scl0 + scl1 * m10 + scl2 * m10**2
offset = off0 + off1 * m10 + off2 * m10**2 # + off_color1p5 * acq_data['color_1p5']
p_fail = offset + scale * wp
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
In [9]:
def fit_sota_model():
from sherpa import ui
data_id = 1
ui.set_method('simplex')
ui.set_stat('cash')
ui.load_user_model(p_acq_fail, '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_all['year']), np.array(data_all['fail'], dtype=np.float))
# Initial fit values from SOTA 2013 presentation (modulo typo there)
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, start_vals.next())
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 [10]:
def plot_fit_grouped(pars, group_col, group_bin, mask=None, log=False, colors='br', label=None):
if mask is None:
mask = np.ones(len(data_all), dtype=bool)
data = data_all[mask]
data['model'] = p_acq_fail(pars, None, data)
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)
fail_sigmas = np.sqrt(data_mean['fail'] * len_groups) / len_groups
plt.errorbar(data_mean[group_col], data_mean['fail'], yerr=fail_sigmas, fmt='.' + colors[0], label=label)
plt.plot(data_mean[group_col], data_mean['model'], '-' + colors[1])
if log:
ax = plt.gca()
ax.set_yscale('log')
In [11]:
def mag_filter(mag0, mag1):
ok = (data_all['mag'] > mag0) & (data_all['mag'] < mag1)
return ok
In [12]:
def wp_filter(wp0, wp1):
ok = (data_all['warm_pix'] > wp0) & (data_all['warm_pix'] < wp1)
return ok
In [13]:
print('Hang tight, this will take a few minutes')
fit = fit_sota_model()
In [14]:
print(fit)
In [15]:
parvals = [par.val for par in model.pars]
parvals
Out[15]:
In [25]:
plot_fit_grouped(parvals, 'mag', 0.25, wp_filter(0.10, 0.20), log=False, colors='cm', label='0.10 < WP < 0.2')
plot_fit_grouped(parvals, 'mag', 0.25, wp_filter(0.0, 0.10), log=False, colors='br', label='0 < WP < 0.10')
plt.legend(loc='best');
plt.ylim(0.0, 1.0);
plt.xlim(9, 11)
plt.grid();
In [17]:
plot_fit_grouped(parvals, 'warm_pix', 0.02, mag_filter(10, 10.6), log=True, colors='cm', label='10 < mag < 10.6')
plot_fit_grouped(parvals, 'warm_pix', 0.02, mag_filter(9, 10), log=True, colors='br', label='9 < mag < 10')
plt.legend(loc='best');
In [18]:
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(10, 10.6), colors='cm', label='10 < mag < 10.6')
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(9.5, 10), colors='br', label='9.5 < mag < 10')
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(9.0, 9.5), colors='gk', label='9.0 < mag < 9.5')
plt.legend(loc='best');
In [19]:
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(10, 10.6), colors='cm', label='10 < mag < 10.6', log=True)
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(9.5, 10), colors='br', label='9.5 < mag < 10', log=True)
plot_fit_grouped(parvals, 'year', 0.25, mag_filter(9.0, 9.5), colors='gk', label='9.0 < mag < 9.5', log=True)
plt.legend(loc='best');
In [20]:
ok = data_all['color'] == 1.5
if np.any(ok):
plot_fit_grouped(parvals, 'year', 0.5, mag_filter(10, 10.6) & ok, colors='cm', label='10 < mag < 10.6')
plot_fit_grouped(parvals, 'year', 0.5, mag_filter(9.5, 10) & ok, colors='br', label='9.5 < mag < 10')
plot_fit_grouped(parvals, 'year', 0.5, mag_filter(9.0, 9.5) & ok, colors='gk', label='9.0 < mag < 9.5')
plt.legend(loc='best');
In [21]:
model.pars
Out[21]:
In [22]:
parvals
Out[22]:
In [23]:
p_acq_fail(SOTA2015_FIT_ONLY_1P5, None, data=None, mag=np.array([9.0]), wp=0.18)
Out[23]:
In [ ]: