From b01f7c98a2f1355aa37be1d8e06922661f2a8e12 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Thu, 26 Mar 2020 09:13:04 -0400 Subject: [PATCH 1/8] replace mag_err with mag_aca_err --- mica/stats/update_guide_stats.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mica/stats/update_guide_stats.py b/mica/stats/update_guide_stats.py index cd5c7c8c..102e139b 100755 --- a/mica/stats/update_guide_stats.py +++ b/mica/stats/update_guide_stats.py @@ -112,7 +112,7 @@ ('var', 'int'), ('pos_err', 'int'), ('mag_aca', 'float'), - ('mag_err', 'int'), + ('mag_aca_err', 'int'), ('mag_band', 'int'), ('pos_catid', 'int'), ('aspq1', 'int'), From 36a70f45fc3fade98c4227b395274d8482cb814f Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Thu, 26 Mar 2020 17:55:23 -0400 Subject: [PATCH 2/8] Set initialization val to 99.0 for mag cols --- mica/stats/update_guide_stats.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/mica/stats/update_guide_stats.py b/mica/stats/update_guide_stats.py index 102e139b..4bc4faac 100755 --- a/mica/stats/update_guide_stats.py +++ b/mica/stats/update_guide_stats.py @@ -451,7 +451,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 'mag' in col: + table[col] = 99.0 + for col in np.dtype(GUIDE_COLS['obs']).names: if col in obsid_info: table[col][:] = obsid_info[col] From 0c6e2ce946efc3863c0cbc79861855152a3c1a35 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Thu, 26 Mar 2020 17:55:51 -0400 Subject: [PATCH 3/8] Add mag percentile cols --- mica/stats/update_guide_stats.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/mica/stats/update_guide_stats.py b/mica/stats/update_guide_stats.py index 4bc4faac..12a55984 100755 --- a/mica/stats/update_guide_stats.py +++ b/mica/stats/update_guide_stats.py @@ -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 @@ -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'), @@ -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]) @@ -327,7 +330,8 @@ 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)]) @@ -456,7 +460,7 @@ def table_gui_stats(obsid_info, gui_stats, star_info, catalog, temp): 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 'mag' in col: + if 'aoacmag' in col: table[col] = 99.0 for col in np.dtype(GUIDE_COLS['obs']).names: From 065aaf0881e8b4136cdd317759c27433b7463ad4 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Fri, 27 Mar 2020 12:11:26 -0400 Subject: [PATCH 4/8] Update script to make it easier to run chunks --- mica/stats/update_guide_stats.py | 49 +++++++++++++++++--------------- 1 file changed, 26 insertions(+), 23 deletions(-) diff --git a/mica/stats/update_guide_stats.py b/mica/stats/update_guide_stats.py index 12a55984..c8aa6ccb 100755 --- a/mica/stats/update_guide_stats.py +++ b/mica/stats/update_guide_stats.py @@ -135,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") @@ -149,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 @@ -344,12 +344,12 @@ def calc_gui_stats(data, star_info): 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: @@ -358,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 @@ -404,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: @@ -413,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( @@ -497,13 +496,15 @@ 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() @@ -511,7 +512,7 @@ def _save_gui_stats(t): 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( @@ -540,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(): From 47867c909d9965a8a82f103746884dbd446fe40e Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Fri, 27 Mar 2020 12:11:47 -0400 Subject: [PATCH 5/8] Update version of guide stats --- mica/stats/update_guide_stats.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mica/stats/update_guide_stats.py b/mica/stats/update_guide_stats.py index c8aa6ccb..e30ab3f9 100755 --- a/mica/stats/update_guide_stats.py +++ b/mica/stats/update_guide_stats.py @@ -33,7 +33,7 @@ if not len(logger.handlers): logger.addHandler(logging.StreamHandler()) -STAT_VERSION = 0.5 +STAT_VERSION = 0.6 GUIDE_COLS = { 'obs': [ From ebd1a4d52aea94872da6dcd584061061d08a2042 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Mon, 6 Apr 2020 12:58:31 -0400 Subject: [PATCH 6/8] TABLE_FILE no longer defined in update_gude_stats --- mica/stats/tests/test_guide_stats.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mica/stats/tests/test_guide_stats.py b/mica/stats/tests/test_guide_stats.py index 1990cb98..6c54e476 100644 --- a/mica/stats/tests/test_guide_stats.py +++ b/mica/stats/tests/test_guide_stats.py @@ -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') From 55137549640b02b1b85676e0b8a38236ba5a651d Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Mon, 6 Apr 2020 14:07:10 -0400 Subject: [PATCH 7/8] Update for new _save_guide_stats api in test --- mica/stats/tests/test_guide_stats.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/mica/stats/tests/test_guide_stats.py b/mica/stats/tests/test_guide_stats.py index 6c54e476..9c95a77a 100644 --- a/mica/stats/tests/test_guide_stats.py +++ b/mica/stats/tests/test_guide_stats.py @@ -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) From 0f11c8d9beb95d1fe9724d8cf7b397e324c97806 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 8 Apr 2020 12:28:14 -0400 Subject: [PATCH 8/8] Update denominator to be the possible kalman samples --- mica/stats/update_guide_stats.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mica/stats/update_guide_stats.py b/mica/stats/update_guide_stats.py index e30ab3f9..e4a522e9 100755 --- a/mica/stats/update_guide_stats.py +++ b/mica/stats/update_guide_stats.py @@ -336,8 +336,8 @@ def calc_gui_stats(data, star_info): 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