SVM-Based Glitch Detection and Removal

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]

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:

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 = []