Draft: WIP Positive/negative standard deviation refinement

I have changed the function find_transition in refine.py look at threshold crossings above and below the channel's mean by a given standard deviation. This makes the refinement more accurate by around 25ms in some cases (see images below) and does not effect the refinement in others. The tuple returned by find_transition is threshold, refined_gps from the first threshold crossing (or only threshold crossing or None if no crossings exist).

The number of standard deviations is still taken from config.py but now the sign doesn't matter, only the magnitude. Signs have been removed and comments added too config.py to indicate that.

find_transition now creates a list dictionary containing a list of first ind where threshold crosses and threshold, sorts that list so the first ind (crossing) comes first and then outputs ind and threshold from the beginning of the list. This list only contains entries if the ind exists.

eg: thresh_crossing_list = [[ind_that_comes_first, associated_threshold],[ind_that_comes_second, associated_threshold]]

eg: thresh_crossing_dict = {'thresh': {'lower': lower_threshold, 'upper': upper_threshold}, 'ind': {'lower': lower_ind, 'upper': upper_ind} }

The key of the minimum ind is retuned as first_crossing and will be upper, lower or None if no crossings exist. The threshold and refined_gps is then retuned using the first_crossing key.

I have removed the logic that was: if no threshold crossing exists, choose the max/min channel value for use in plotting the threshold. Do we think that is okay to do?

Closes #158

Old version, refined gps = .0688

New version, refined gps = .0429 (i.e. 25.9ms sooner)

Additionally, I have changed to input params of plot_indicators from taking the dictionary params to just taking channel. Previously the plotting looked at whether the upper or lower threshold was used in order to place the annotation "threshold" above or below the plotted line. This seemed tricky to do with the variable thresholds so I removed the dependence and put the annotation always below. This doesn't always look as smart if anyone can think of a better way of doing this Plotting will still look at whether threshold is above or below mean and plot "threshold" annotation in the correct place.

locklost/plugins/refine.py

 import logging
 import numpy as np
 import matplotlib.pyplot as plt
+from collections import defaultdict
 
 from gwpy.segments import Segment
 
@@ -9,7 +10,7 @@ from .. import data
 
 ##################################################
 
-def plot_indicators(event, params, refined_gps=None, threshold=None):
+def plot_indicators(event, channel, refined_gps, threshold):
     """Create graphs showing indicators of refined time.
 
     """
@@ -18,7 +19,7 @@ def plot_indicators(event, params, refined_gps=None, threshold=None):
 
     channels = [
         config.GRD_STATE_N_CHANNEL,
-        params['CHANNEL'],
+        channel,
     ]
     window = config.REFINE_PLOT_WINDOWS['WIDE']
     segment = Segment(*window).shift(t0)
@@ -27,6 +28,11 @@ def plot_indicators(event, params, refined_gps=None, threshold=None):
     y_grd, t_grd = bufs[0].yt()
     y_ind, t_ind = bufs[1].yt()
 
+    # calculate mean using first 5 seconds of y_ind data
+    samples = int(y_ind.sample_rate * 5)
+    locked_sample = bufs[1][0][:samples]
+    locked_mean = np.mean(locked_sample)
+
     plt.rcParams['text.usetex'] = False
 
     isc_loss_time = event.transition_gps - t0
@@ -43,7 +49,7 @@ def plot_indicators(event, params, refined_gps=None, threshold=None):
 
         # plot indicator
         color = 'blue'
-        label = params['CHANNEL']
+        label = channel
         # filter data based on window
         condition = np.logical_and(t0+window[0] <= t_ind, t_ind <= t0+window[1])
         idx = np.argwhere(condition)
@@ -52,7 +58,7 @@ def plot_indicators(event, params, refined_gps=None, threshold=None):
         ax1.plot(t-t0, y, label=label,
                  color=color, alpha=0.8, linewidth=2)
         ax1.grid(True)
-        ax1.set_ylabel(params['CHANNEL'])
+        ax1.set_ylabel(channel)
         ax1.yaxis.label.set_color(color)
         ax1.set_xlabel('Time [s] since {}'.format(t0))
 
@@ -65,7 +71,7 @@ def plot_indicators(event, params, refined_gps=None, threshold=None):
                 linestyle='--',
                 linewidth=2,
             )
-            if params['THRESHOLD'] > 0:
+            if locked_mean < threshold:
                 xyoffset = (8, 10)
                 valign = 'bottom'
             else:
@@ -159,6 +165,7 @@ def find_transition(channel, segment, std_threshold, minimum=None):
     returns `threshold`, `refined_gps` tuple
 
     """
+    threshold = None
     refined_gps = None
     buf = data.fetch(channel, segment)[0]
 
@@ -170,30 +177,43 @@ def find_transition(channel, segment, std_threshold, minimum=None):
     logging.info("locked mean/stddev: {}/{}".format(lock_mean, lock_stdd))
 
     # lock loss threshold is when moves past % threshold of std
-    threshold = lock_stdd * std_threshold + lock_mean
-    if std_threshold > 0:
-        threshold = min(threshold, max(buf.data))
-    else:
-        threshold = max(threshold, min(buf.data))
-    logging.info("threshold: {}".format(threshold))
+    upper_threshold = lock_stdd * abs(std_threshold) + lock_mean
+    lower_threshold = lock_stdd * -abs(std_threshold) + lock_mean
 
     # if the mean is less than the nominal full lock value then abort,
     # as this isn't a clean lock
     if minimum and lock_mean < minimum:
         logging.info("channel mean below minimum, unable to resolve time")
 
+    # look for positive and negative threshold crossings, add them to a
+    # dictionary, then choose the first crossing if one exists
     else:
-        if std_threshold > 0:
-            inds = np.where(buf.data > threshold)[0]
+        inds_above = np.where(buf.data > upper_threshold)[0]
+        inds_below = np.where(buf.data < lower_threshold)[0]
+        thresh_crossing_dict = defaultdict(lambda: defaultdict(int))
+        if inds_above.any():
+            thresh_crossing_dict['ind']['upper'] = np.min(inds_above)
+            thresh_crossing_dict['thresh']['upper'] = upper_threshold
+        if inds_below.any():
+            thresh_crossing_dict['ind']['lower'] = np.min(inds_below)
+            thresh_crossing_dict['thresh']['lower'] = lower_threshold
+        first_crossing = min(
+                thresh_crossing_dict['ind'],
+                key = lambda k: thresh_crossing_dict['ind'][k],
+                default = None
+                )
+        if first_crossing is None:
+            logging.info("no threshold crossings, unable to resolve time")
         else:
-            inds = np.where(buf.data < threshold)[0]
-
-        if inds.any():
-            ind = np.min(inds)
+            threshold = thresh_crossing_dict['thresh'][first_crossing]
+            ind = thresh_crossing_dict['ind'][first_crossing]
             refined_gps = buf.tarray[ind]
+            logging.info("threshold: {}".format(threshold))
             logging.info("refined time: {}".format(refined_gps))
-        else:
-            logging.info("no threshold crossings, unable to resolve time")
+
+    # add a threshold to be diplayed on graphs if refinement failed
+    if threshold is None:
+        threshold = lower_threshold
 
     return threshold, refined_gps
 
@@ -234,7 +254,7 @@ def refine_time(event):
     logging.info("plotting indicators...")
     plot_indicators(
         event,
-        params,
+        params['CHANNEL'],
         refined_gps=refined_gps,
         threshold=threshold,
     )