made slope thresholding a command-line option

pemcoupling/coupfunc.py

@@ -152,7 +152,8 @@ class CoupFunc(PEMFrequencySeries):
             injection_type='broadband',
             peak_search_window=0, smooth_params=(0, 0), smooth_log=False,
             notch_windows=[], broadband_freqs=None, comb_freq=None,
-            line_freq=None, injection_name=None, verbose=False):
+            line_freq=None, injection_name=None, slope_threshold=None,
+            verbose=False):
         """
         Calculates coupling factors from sensor spectra and DARM spectra.
 
@@ -277,19 +278,11 @@ class CoupFunc(PEMFrequencySeries):
                 notch_mask[(lo < freqs) & (freqs < hi)] = False
             freq_mask = (freq_mask & notch_mask)
         # Slope threshold
-        slope_window = 10
-        slope_thresh = 0.01
-        moving_avg = np.convolve(sens_inj_asd, np.ones((slope_window,)) / (slope_window), mode='same')
-        inj_slope = np.zeros_like(sens_inj_asd)
-        inj_slope[slope_window / 2: -slope_window / 2] = moving_avg[slope_window:] - moving_avg[:-slope_window]
-        inj_slope /= (slope_window * 2)
-#         inj_slope = np.zeros_like(sens_inj_asd)
-#         inj_slope[slope_window: -slope_window] = sens_inj_asd[slope_window * 2:] - sens_inj_asd[:-slope_window * 2]
-#         inj_slope /= (slope_window * 2)
-        inj_slope_below_thresh = np.abs(inj_slope) / sens_inj_asd < slope_thresh
-        freq_mask = (freq_mask & inj_slope_below_thresh)
-#         print(inj_slope[(freqs > 100) & (freqs < 102)])
-#         print(inj_slope_below_thresh[(freqs > 100) & (freqs < 102)])
+        if slope_threshold is not None:
+            slope_window = 0.03
+            inj_slope = pem_utils.asd_slope(sens_inj_asd, window=slope_window)
+            inj_slope_below_thresh = (inj_slope < slope_threshold) & (inj_slope > 1.0 / slope_threshold)
+            freq_mask = (freq_mask & inj_slope_below_thresh)
         #######################################################################
         # COMPUTE COUPLING FACTORS WHERE APPLICABLE
         # sens_ratio = sens_inj / np.maximum(sens_bkg, sens_bkg_ref)

pemcoupling/injection.py

@@ -31,6 +31,7 @@ class Injection(object):
             calibration_units={},
             darm_threshold=2,
             sensor_threshold=2,
+            slope_threshold=None,
             notch_windows=[],
             peak_search_window=0,
             broadband_freqs=None,
@@ -60,6 +61,7 @@ class Injection(object):
         self.darm_threshold = darm_threshold
         self.sensor_threshold = sensor_threshold
         self.sensor_threshold_autoscale = sensor_threshold_autoscale
+        self.slope_threshold = slope_threshold
         self.notch_windows = notch_windows
         self.peak_search_window = peak_search_window
         self.broadband_freqs = broadband_freqs
@@ -357,6 +359,7 @@ class Injection(object):
             'darm_threshold',
             'sensor_threshold',
             'sensor_threshold_autoscale',
+            'slope_threshold',
             'injection_type',
             'smooth_dict',
             'notch_windows',

pemcoupling/main.py

@@ -116,7 +116,7 @@ def main(args):
         try:
             sect_cfg = config_dict[section]
         except:
-            logging.wargnin("invalid config section " + args.section)
+            logging.warning("invalid config section " + args.section)
             continue
         logging.info("Handling conflicts between config options "
                      "and command line options.")
@@ -154,6 +154,7 @@ def main(args):
             'sensor_threshold_autoscale',
             sect_cfg.get('sensor_threshold_autoscale'))
         peak_search_window = sect_cfg.get('peak_search_window')
+        slope_threshold = arg_dict.get('slope_threshold')
         # Injection types and frequency ranges
         broadband_freqs_custom = None
         comb_freq_custom = None
@@ -317,6 +318,7 @@ def main(args):
                 darm_threshold=darm_threshold,
                 sensor_threshold=sensor_threshold,
                 sensor_threshold_autoscale=sensor_threshold_autoscale,
+                slope_threshold=slope_threshold,
                 notch_windows=darm_notch_data,
                 peak_search_window=peak_search_window,
                 broadband_freqs=broadband_freqs,

pemcoupling/parse.py

@@ -183,6 +183,9 @@ def parse_args(script, args=None):
             "--line", nargs="?", type=float, const=True, default=False,
             help="measure coupling only near this frequency (frequency " +
             "is inferred from injection name if a value is not given here)")
+        analysis_args.add_argument(
+            "--slope-threshold", type=float,
+            help="slope threshold for coupling function")
         plot_args.add_argument(
             "-r", "--ratio-plot", action="store_true", default=False,
             help="create ratio plot")

pemcoupling/utils.py

@@ -454,8 +454,19 @@ def near_line(line_freq, freqs, delta_f):
     return np.abs(freqs - line_freq) < (0.5 * delta_f)
 
 
-def in_freq_band(band, freqs, pad=0.1):
-    fmin = (1 - pad) * band[0]
-    fmax = (1 + pad) * band[1]
+def in_freq_band(band, freqs, pad=0):
+    fmin = (1.0 - pad) * band[0]
+    fmax = (1.0 + pad) * band[1]
     out = ((freqs > fmin) & (freqs < fmax))
     return out
+
+
+def asd_slope(asd, window):
+    slope = np.ones_like(asd)
+    for i in range(1, slope.size - 1):
+        lo = max(0, i - int(window * i))
+        hi = min(slope.size, i + int(window * i))
+        avg_below = asd[lo:i].mean()
+        avg_above = asd[i:hi].mean()
+        slope[i] = avg_above / avg_below
+    return slope