Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Mag stats #219

Merged
merged 8 commits into from
Apr 14, 2020
Merged

Mag stats #219

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
b01f7c9
replace mag_err with mag_aca_err
jeanconn Mar 26, 2020
36a70f4
Set initialization val to 99.0 for mag cols
jeanconn Mar 26, 2020
0c6e2ce
Add mag percentile cols
jeanconn Mar 26, 2020
065aaf0
Update script to make it easier to run chunks
jeanconn Mar 27, 2020
47867c9
Update version of guide stats
jeanconn Mar 27, 2020
ebd1a4d
TABLE_FILE no longer defined in update_gude_stats
jeanconn Apr 6, 2020
5513754
Update for new _save_guide_stats api in test
jeanconn Apr 6, 2020
0f11c8d
Update denominator to be the possible kalman samples
jeanconn Apr 8, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 2 additions & 3 deletions mica/stats/tests/test_guide_stats.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
from .. import guide_stats
from .. import update_guide_stats

HAS_GS_TABLE = os.path.exists(update_guide_stats.TABLE_FILE)
HAS_GS_TABLE = os.path.exists(guide_stats.TABLE_FILE)


@pytest.mark.skipif('not HAS_GS_TABLE', reason='Test requires guide stats table')
Expand Down Expand Up @@ -44,9 +44,8 @@ def test_make_gui_stats():
# make a new table if the supplied file doesn't exist
fh, fn = tempfile.mkstemp(suffix='.h5')
os.unlink(fn)
update_guide_stats.TABLE_FILE = fn
obsid = 20001
obsid_info, gui, star_info, catalog, temp = update_guide_stats.calc_stats(obsid)
t = update_guide_stats.table_gui_stats(obsid_info, gui, star_info, catalog, temp)
update_guide_stats._save_gui_stats(t)
update_guide_stats._save_gui_stats(t, table_file=fn)
os.unlink(fn)
78 changes: 46 additions & 32 deletions mica/stats/update_guide_stats.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,6 @@
import argparse
import numpy as np
import tables
import tables3_api
from scipy.stats import scoreatpercentile
from astropy.table import Table
import logging
import warnings
Expand Down Expand Up @@ -35,7 +33,7 @@
if not len(logger.handlers):
logger.addHandler(logging.StreamHandler())

STAT_VERSION = 0.5
STAT_VERSION = 0.6

GUIDE_COLS = {
'obs': [
Expand Down Expand Up @@ -78,6 +76,11 @@
('aoacmag_min', 'float'),
('aoacmag_mean', 'float'),
('aoacmag_max', 'float'),
('aoacmag_5th', 'float'),
('aoacmag_16th', 'float'),
('aoacmag_50th', 'float'),
('aoacmag_84th', 'float'),
('aoacmag_95th', 'float'),
('aoacmag_std', 'float'),
('aoacyan_mean', 'float'),
('aoaczan_mean', 'float'),
Expand Down Expand Up @@ -112,7 +115,7 @@
('var', 'int'),
('pos_err', 'int'),
('mag_aca', 'float'),
('mag_err', 'int'),
('mag_aca_err', 'int'),
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @taldcroft and @javierggt ! It sounds like you also want a change in this PR somewhere to scale the mag_aca_err column to units of mag. That's fine. Though should that be in mica or in agasc.get_star (or get_stars) with a calibrate flag or some such?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well given that all the rest of this is uncalibrated it appears we cannot make a change here without having a one-off inconsistency. Just leave it. And don't bother with making an agasc issue. If one of us gets annoyed enough one day we'll do it, but it is a low priority and we don't need more low priority issues.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be less inconsistent to just UPPER CASE these?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be more consistent to upper case them, but would break all existing code that uses these. Not worth the trouble I'd say.

('mag_band', 'int'),
('pos_catid', 'int'),
('aspq1', 'int'),
Expand All @@ -132,12 +135,6 @@
}


SKA = os.environ['SKA']
DATA = os.path.join(SKA, 'data', 'guide_stats')
TABLE_FILE = os.path.join(DATA, 'guide_stats.h5')
SKIPPED_FILE = os.path.join(DATA, 'skipped.dat')


def get_options():
parser = argparse.ArgumentParser(
description="Update guide stats table")
Expand All @@ -146,6 +143,12 @@ def get_options():
help="check for missing observations in table and reprocess")
parser.add_argument("--obsid",
help="specific obsid to process. Not required in regular update mode")
parser.add_argument("--start",
help="start time for processing")
parser.add_argument("--stop",
help="stop time for processing")
parser.add_argument("--datafile",
default="gs.h5")
opt = parser.parse_args()
return opt

Expand Down Expand Up @@ -315,8 +318,8 @@ def calc_gui_stats(data, star_info):
stats['star_tracked'] = np.any(dr < 5.0)
stats['spoiler_tracked'] = np.any(dr > 5.0)
deltas = {'dy': dy, 'dz': dz, 'dr': dr}
stats['dr_5th'] = scoreatpercentile(deltas['dr'], 5)
stats['dr_95th'] = scoreatpercentile(deltas['dr'], 95)
stats['dr_5th'] = np.percentile(deltas['dr'], 5)
stats['dr_95th'] = np.percentile(deltas['dr'], 95)
for ax in deltas:
stats['{}_mean'.format(ax)] = np.mean(deltas[ax])
stats['{}_std'.format(ax)] = np.std(deltas[ax])
Expand All @@ -327,25 +330,26 @@ def calc_gui_stats(data, star_info):
stats['aoacmag_mean'] = np.mean(mag)
stats['aoacmag_max'] = np.max(mag)
stats['aoacmag_std'] = np.std(mag)

for perc in [5, 16, 50, 84, 95]:
stats[f'aoacmag_{perc}th'] = np.percentile(mag, perc)
stats['aoacyan_mean'] = np.mean(kal['AOACYAN{}'.format(slot)])
stats['aoaczan_mean'] = np.mean(kal['AOACZAN{}'.format(slot)])

for dist in ['0.3', '1', '3', '5']:
stats['f_within_{}'.format(dist)] = np.count_nonzero(dr < float(dist)) / stats['n_kalman']
stats['f_outside_5'] = np.count_nonzero(dr > 5) / stats['n_kalman']
stats['f_within_{}'.format(dist)] = np.count_nonzero(dr < float(dist)) / len(kal)
stats['f_outside_5'] = np.count_nonzero(dr > 5) / len(kal)

gui_stats[slot] = stats

return gui_stats


def _get_obsids_to_update(check_missing=False):
def _get_obsids_to_update(check_missing=False, table_file=None, start=None, stop=None):
if check_missing:
last_tstart = '2007:271'
last_tstart = start if start is not None else '2007:271'
kadi_obsids = events.obsids.filter(start=last_tstart)
try:
h5 = tables.open_file(TABLE_FILE, 'r')
h5 = tables.open_file(table_file, 'r')
tbl = h5.root.data[:]
h5.close()
except:
Expand All @@ -354,13 +358,13 @@ def _get_obsids_to_update(check_missing=False):
obsids = [o.obsid for o in kadi_obsids if o.obsid not in tbl['obsid']]
else:
try:
h5 = tables.open_file(TABLE_FILE, 'r')
h5 = tables.open_file(table_file, 'r')
tbl = h5.get_node('/', 'data')
last_tstart = tbl.cols.kalman_tstart[tbl.colindexes['kalman_tstart'][-1]]
h5.close()
except:
last_tstart = '2002:012'
kadi_obsids = events.obsids.filter(start=last_tstart)
last_tstart = start if start is not None else '2002:012'
kadi_obsids = events.obsids.filter(start=last_tstart, stop=stop)
# Skip the first obsid (as we already have it in the table)
obsids = [o.obsid for o in kadi_obsids][1:]
return obsids
Expand Down Expand Up @@ -400,7 +404,7 @@ def calc_stats(obsid):
if not manvr or not dwell:
raise ValueError("No manvr or dwell for {}".format(obsid))
if not manvr.get_next():
raise ValueError("No *next* manvr so can't calculate dwell".format(obsid))
raise ValueError("No *next* manvr so can't calculate dwell")
if not manvr.guide_start:
raise ValueError("No guide transition for {}".format(obsid))
if not manvr.kalman_start:
Expand All @@ -409,7 +413,6 @@ def calc_stats(obsid):

logger.info("Found obsid manvr at {}".format(manvr.start))
logger.info("Found dwell at {}".format(dwell.start))
guide_start = manvr.guide_start
starcheck = get_starcheck_catalog_at_date(manvr.guide_start)
if starcheck is None or 'cat' not in starcheck or not len(starcheck['cat']):
raise ValueError('No starcheck catalog found for {}'.format(
Expand Down Expand Up @@ -451,7 +454,14 @@ def table_gui_stats(obsid_info, gui_stats, star_info, catalog, temp):
logger.info("arranging stats into tabular data")
cols = (GUIDE_COLS['obs'] + GUIDE_COLS['cat'] + GUIDE_COLS['stat']
+ GUIDE_COLS['agasc'] + GUIDE_COLS['temp'] + GUIDE_COLS['bad'])

# Initialize all values to zero
table = Table(np.zeros((1, 8), dtype=cols).flatten())
# Set all columns with mag info to 99.0 initial value instead of zero
for col in table.dtype.names:
if 'aoacmag' in col:
table[col] = 99.0
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is somewhat reasonable to "initialize" the numpy table with values of 99 for all the aoacmag_* cols if we want 99.0 for that kind of missing data. Not sure if other columns should get that treatment too.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is probably the only controversial thing in here. As with everything, hard to know how much to rummage around and improve things and how much to just back away slowly.


for col in np.dtype(GUIDE_COLS['obs']).names:
if col in obsid_info:
table[col][:] = obsid_info[col]
Expand Down Expand Up @@ -486,21 +496,23 @@ def table_gui_stats(obsid_info, gui_stats, star_info, catalog, temp):
return table


def _save_gui_stats(t):
if not os.path.exists(TABLE_FILE):
def _save_gui_stats(t, table_file=None):
if table_file is None:
return
if not os.path.exists(table_file):
cols = (GUIDE_COLS['obs'] + GUIDE_COLS['cat'] + GUIDE_COLS['stat']
+ GUIDE_COLS['agasc'] + GUIDE_COLS['temp'] + GUIDE_COLS['bad'])
desc, byteorder = tables.descr_from_dtype(np.dtype(cols))
filters = tables.Filters(complevel=5, complib='zlib')
h5 = tables.open_file(TABLE_FILE, 'a')
h5 = tables.open_file(table_file, 'a')
tbl = h5.create_table('/', 'data', desc, filters=filters,
expectedrows=1e6)
tbl.cols.obsid.create_index()
tbl.cols.kalman_tstart.create_csindex()
tbl.cols.agasc_id.create_index()
h5.close()
del h5
h5 = tables.open_file(TABLE_FILE, 'a')
h5 = tables.open_file(table_file, 'a')
tbl = h5.get_node('/', 'data')
have_obsid_coord = tbl.get_where_list(
'(obsid == {}) & (obi == {})'.format(
Expand Down Expand Up @@ -529,23 +541,25 @@ def update(opt):
if opt.obsid:
obsids = [int(opt.obsid)]
else:
obsids = _get_obsids_to_update(check_missing=opt.check_missing)
obsids = _get_obsids_to_update(table_file=opt.datafile, check_missing=opt.check_missing,
start=opt.start, stop=opt.stop)
for obsid in obsids:
t = time.localtime()

logger.info("Processing obsid {}".format(obsid))
try:
obsid_info, gui_stats, star_info, guide_catalog, temp = calc_stats(obsid)
except Exception as e:
open(SKIPPED_FILE, 'a').write("{}: {}\n".format(obsid, e))
open(os.path.splitext(opt.datafile)[0] + '_skipped.dat', 'a').write(
"{}: {}\n".format(obsid, e))
logger.info("Skipping obsid {}: {}".format(obsid, e))
continue
if not len(gui_stats):
open(SKIPPED_FILE, 'a').write("{}: No stats\n".format(obsid))
open(os.path.splitext(opt.datafile)[0] + '_skipped.dat', 'a').write(
"{}: No stats\n".format(obsid))
logger.info("Skipping obsid {}, no stats determined".format(obsid))
continue
t = table_gui_stats(obsid_info, gui_stats, star_info, guide_catalog, temp)
_save_gui_stats(t)
_save_gui_stats(t, opt.datafile)


def main():
Expand Down