diff --git a/nitrates/config.py b/nitrates/config.py index 4483365..75158e6 100644 --- a/nitrates/config.py +++ b/nitrates/config.py @@ -52,6 +52,9 @@ # Table of bright known sources from the Trans Monitor bright_source_table_fname = os.path.join(dir, "data", "bright_src_cat.fits") +# Support Vector Model used in glitch_btis +SVM_model = os.path.join(dir, "data", "svm_original_model.pkl") + EBINS0 = [15.0, 24.0, 35.0, 48.0, 64.0, 84.0, 120.0, 171.5, 245.0] EBINS1 = [24.0, 35.0, 48.0, 64.0, 84.0, 120.0, 171.5, 245.0, 350.0] diff --git a/nitrates/data/svm_original_model.pkl b/nitrates/data/svm_original_model.pkl new file mode 100644 index 0000000..6ba5c9f Binary files /dev/null and b/nitrates/data/svm_original_model.pkl differ diff --git a/nitrates/data_prep/do_data_setup.py b/nitrates/data_prep/do_data_setup.py index 8143dd9..f74cee1 100644 --- a/nitrates/data_prep/do_data_setup.py +++ b/nitrates/data_prep/do_data_setup.py @@ -11,6 +11,7 @@ from datetime import datetime import argparse import logging, traceback +import pickle from ..lib.time_funcs import met2astropy, utc2met, met2utc_str, apy_time2met from ..lib.sqlite_funcs import ( @@ -45,11 +46,13 @@ get_btis_for_glitches, check_if_in_GTI, find_and_remove_cr_glitches, + SVM_dpi_eval, ) from ..data_scraping.db_ql_funcs import get_gainoff_fname from ..data_scraping.api_funcs import get_sao_file from ..HeasoftTools.bat_tool_funcs import bateconvert from ..lib.search_config import Config +from ..config import SVM_model def query_data_metslice(conn, met0, met1, table_name="SwiftQLevent"): sql = """SELECT * FROM %s @@ -148,10 +151,17 @@ def evfnames2write( else: gti_tot = all_gtis[0] + with open(SVM_model, "rb") as f: + clf = pickle.load(f) + glitch_btis = get_btis_for_glitches( ev_data0, gti_tot["START"][0], gti_tot["STOP"][-1] ) - for bti in glitch_btis: + + # new function added to evaluate glitches with a ML classifier + realglitch_btis = SVM_dpi_eval(ev_data0, glitch_btis, clf) + + for bti in realglitch_btis: logging.info("Found glitch bti: ") logging.info(bti) if len(gti_pnt) >= 1: diff --git a/nitrates/lib/gti_funcs.py b/nitrates/lib/gti_funcs.py index babb845..2f40164 100644 --- a/nitrates/lib/gti_funcs.py +++ b/nitrates/lib/gti_funcs.py @@ -1,7 +1,6 @@ from astropy.table import Table, vstack import numpy as np -import logging - +import logging, pickle, os def union_gtis(gti_tabs): """Union of overlapping time intervals. @@ -126,8 +125,14 @@ def mk_gti_bl(times, GTI, time_pad=0.0): bl = bl | bl_ return bl - -def get_btis_for_glitches(evdata, tstart, tstop, tbin_size=16e-3): +# SNR threshold editing enabled as arguments for glitch energy parameters +def get_btis_for_glitches(evdata, tstart, tstop, tbin_size=16e-3, lowE_snr_thresh=6.0, snr_ratio_thresh=2.0): + """ Finds where heat pipe glitches occur depending on SNR thresholds. + + Deafult set to lowE = 6 and the ratio = 2...algorithm was + previously hard-coded 10 and 3. + Returns list of possible glitch BTIS expanded to 64ms. + """ bins = np.arange(tstart, tstop + tbin_size / 2.0, tbin_size) ebl = evdata["ENERGY"] <= 25.0 ebl2 = evdata["ENERGY"] > 50.0 @@ -138,9 +143,9 @@ def get_btis_for_glitches(evdata, tstart, tstop, tbin_size=16e-3): stds2 = (h2 - np.mean(h2)) / np.std(h2) # bl_bad = (stds>10.0)&(stds2<2.5) - bl_lowE_highSNR = stds > 10.0 + bl_lowE_highSNR = stds > lowE_snr_thresh # bl_highE_lowSNR = (stds2<2.5)|((stds/stds2)>5) - bl_highE_lowSNR = (stds / np.abs(stds2)) > 3 + bl_highE_lowSNR = (stds / np.abs(stds2)) > snr_ratio_thresh logging.debug("N_lowE_highSNR: " + str(np.sum(bl_lowE_highSNR))) if np.sum(bl_lowE_highSNR) > 0: logging.debug("LowE highSNRs: ") @@ -156,13 +161,49 @@ def get_btis_for_glitches(evdata, tstart, tstop, tbin_size=16e-3): t1 = t0 + tbin_size bad_twind = ( - bins[:-1][bl_bad][i] - tbin_size / 2.0, - bins[:-1][bl_bad][i] + 1.5 * tbin_size, + bins[:-1][bl_bad][i] - (tbin_size / 2.0) - tbin_size, + bins[:-1][bl_bad][i] + (1.5 * tbin_size) + tbin_size, ) bad_twinds.append(bad_twind) return bad_twinds +def SVM_dpi_eval(evtable, possibad_twinds, clf, threshold=0.45): + """ Evaluates glitch probability and returns those with >.45 + glitch probability to the rest of the pipeline. + + Makes a detector plane image (DPI) out of counts from 64ms of event + data and determines glitch pattern likelihood. + + glitch_btis greater than .45 glitch probability get returned as + a list of true btis as 'realglitch_btis' + """ + xbins = np.arange(286 + 1) - 0.5 + ybins = np.arange(173 + 1) - 0.5 + + real_bads = [] + + for tmin, tmax in possibad_twinds: + bl = (evtable['TIME'] >= tmin) & (evtable['TIME'] <= tmax) + binned_events = evtable[bl] + + if len(binned_events) == 0: + continue + + dpi = np.histogram2d(binned_events['DETX'], binned_events['DETY'], bins=[xbins, ybins])[0] + flatdpi = dpi.ravel().reshape(1, -1) + + # class 0 are glitches, 1 are non-glitch + probabilities = clf.predict_proba(flatdpi)[0] + glitch_prob = probabilities[0] + logging.debug("Glitch probability: %.4f" % glitch_prob) + + if glitch_prob > threshold: + real_bads.append((tmin, tmax)) + + return(real_bads) + + def find_cr_glitch_times(ev_data, tmin, tmax, tbin_size=5e-5, emin=50, max_cnts=20): bad_times = []