Remove gwpy

bin/live/pycbc_live_collated_dq_trigger_rates

@@ -25,13 +25,12 @@ import glob
 
 import numpy as np
 
-import gwdatafind
+from igwn_segments import segmentlist, segment
 
 import pycbc
+import pycbc.dq
 from pycbc.io import HFile
-
-from gwpy.timeseries import TimeSeriesDict
-from gwpy.segments import Segment, DataQualityFlag
+from pycbc.frame.frame import read_frame, query_and_read_frame
 
 
 def find_frames(ifo, frame_dir, frame_type, start, end):
@@ -70,7 +69,7 @@ parser.add_argument("--frame-directory",
                     help='Directory containing frame files. Required if '
                          'not using gwdatafind.')
 parser.add_argument("--frame-type", required=True)
-parser.add_argument("--analysis-flag-name", required=True)
+parser.add_argument("--analysis-flag-name")
 parser.add_argument("--dq-channel", required=True)
 parser.add_argument("--dq-ok-channel", required=True)
 parser.add_argument("--dq-thresh", required=True, type=float)
@@ -85,30 +84,36 @@ pycbc.init_logging(args.verbose)
 if not args.use_gwdatafind and args.frame_directory is None:
     raise ValueError('Must provide frame-directory if not using gwdatafind.')
 
-# Get observing segs
-ar_flag_name = args.analysis_flag_name.format(ifo=args.ifo)
-day_seg = Segment(args.gps_start_time, args.gps_end_time)
-observing_flag = DataQualityFlag.query(ar_flag_name, day_seg)
+if args.analysis_flag_name is not None:
+    # Get observing segs
+    ar_flag_name = args.analysis_flag_name.format(ifo=args.ifo)
+    observing_flag = pycbc.dq.query_flag(
+        args.ifo, ar_flag_name, args.gps_start_time, args.gps_end_time
+    )
+else:
+    # If we are not querying the database for segments, just assume that the
+    # full time is in observing mode
+    observing_flag = segmentlist(
+        [segment(args.gps_start_time, args.gps_end_time)]
+    )
 
 # shift observing flag to current time
-observing_flag = observing_flag.pad(args.replay_offset, args.replay_offset)
+observing_flag.protract(args.replay_offset)
 
-observing_segs = observing_flag.active
-livetime = observing_flag.livetime
+livetime = abs(observing_flag)
 logging.info(f'Found {livetime} seconds of observing time at {args.ifo}.')
 
 # for each segment, check how much time was dq flagged
 flagged_time = 0
 dq_channel = args.dq_channel.format(ifo=args.ifo)
 dq_ok_channel = args.dq_ok_channel.format(ifo=args.ifo)
-for seg in observing_segs:
+for seg in observing_flag:
     if args.use_gwdatafind:
-        frame_type = args.analysis_frame_type.format(ifo=args.ifo)
-        frames = gwdatafind.find_urls(
-            args.ifo[0],
-            frame_type,
-            seg[0],
-            seg[1],
+        ts_dq_channel, ts_dq_ok_channel = query_and_read_frame(
+            args.frame_type.format(ifo=args.ifo),
+            [dq_channel, dq_ok_channel],
+            start_time=seg[0],
+            end_time=seg[1]
         )
     else:
         frames = find_frames(
@@ -119,24 +124,31 @@ for seg in observing_segs:
             seg[1],
         )
 
-    tsdict = TimeSeriesDict.read(
-        frames,
-        channels=[dq_channel, dq_ok_channel],
-        start=seg[0],
-        end=seg[1],
-        pad=0,
-    )
+        ts_dq_channel, ts_dq_ok_channel = read_frame(
+            frames, [dq_channel, dq_ok_channel],
+            start_time=seg[0], end_time=seg[1]
+        )
 
-    dq_ok_flag = (tsdict[dq_ok_channel] == 1).to_dqflag()
-    dq_flag = (tsdict[dq_channel] <= args.dq_thresh).to_dqflag()
+    # build segmentlists:
+    # Is the channel actually working?
+    dq_ok_segs = ts_dq_ok_channel.bool_to_segmentlist()
 
-    valid_bad_dq = dq_flag & dq_ok_flag
+    # Is the channel above threshold?
+    dq_mask = (ts_dq_channel <= args.dq_thresh)
+    # If not, set it to zero
+    ts_dq_channel.data[dq_mask] = 0
+    dq_segs = ts_dq_channel.bool_to_segmentlist()
 
-    # pad flag outwards to match the padding used in the live script
+    # these are the times where dq channel is working, and
+    # the dq channel is above threshold
+    valid_bad_dq = dq_segs & dq_ok_segs
+
+    # pad / expand segments outwards by dq_padding seconds
     valid_bad_dq.protract(args.dq_padding)
     valid_bad_dq.coalesce()
 
-    flagged_time += (valid_bad_dq & observing_flag).livetime
+    # intersection with observing_flag and add livetime
+    flagged_time += abs(valid_bad_dq & observing_flag)
 
 logging.info(f'Found {flagged_time} seconds of dq flagged time at {args.ifo}.')
 
@@ -215,3 +227,5 @@ with HFile(args.output, 'w') as f:
     setting_str = ' '.join([str(s) for s in settings_to_hash])
     hash_object = hashlib.sha256(setting_str.encode())
     f.attrs['settings_hash'] = hash_object.hexdigest()
+
+logging.info('Done')

companion.txt

@@ -1,6 +1,3 @@
-# other tools which may be useful
-gwpy>=0.8.1
-
 # HEALPix is very useful for some analysis.
 healpy
 

pycbc/types/timeseries.py

@@ -22,6 +22,7 @@
 
 import numpy as _numpy
 from scipy.io.wavfile import write as write_wav
+from igwn_segments import segmentlist, segment
 
 from pycbc.types.array import Array, _convert, complex_same_precision_as, zeros
 from pycbc.types.utils import determine_epoch
@@ -1234,6 +1235,48 @@ def plot(self, **kwds):
             plot2 = pyplot.plot(self.sample_times, self.imag(), **kwds)
             return plot1, plot2
 
+    def bool_to_segmentlist(self, epoch=None):
+        """
+        Convert a boolean pycbc TimeSeries (Truth-like when condition holds) to
+        an igwn_segments.segmentlist of (start, end) in GPS seconds.
+        """
+
+        # Is the data truthlike?
+        # bools or numbers are OK, but we require finite values
+        arr = self.numpy()
+        if arr.dtype.kind != 'b' and not _numpy.issubdtype(arr.dtype, _numpy.number):
+            raise RuntimeError(
+                'To use bool_to_segmentlist, we require that the timeseries '
+                'is truthlike'
+            )
+
+        segs = segmentlist([])
+
+        if arr.size == 0:
+            return segs.coalesce()
+
+        # Convert to bool Treat NaNs as zeros
+        b = _numpy.nan_to_num(arr, nan=0).astype(bool)
+
+        # Pad with a leading/trailing False to detect edges at boundaries.
+        b = _numpy.concatenate(([False], b, [False]))
+
+        # Work out the transitions between true and false.
+        # starts = False to True transitions
+        starts = _numpy.flatnonzero((~b[:-1]) & b[1:])
+        # ends = True to False Transitions
+        ends = _numpy.flatnonzero(b[:-1] & (~b[1:])) 
+
+        # Convert indices to GPS times
+        starts_time = self.start_time + starts * self.delta_t
+        ends_time = self.start_time + ends * self.delta_t
+
+        # Convert to segments
+        for s, e in zip(starts_time, ends_time):
+            segs.append(segment(s, e))
+
+        return segs.coalesce()
+
 def load_timeseries(path, group=None):
     """Load a TimeSeries from an HDF5, ASCII or Numpy file. The file type is
     inferred from the file extension, which must be `.hdf`, `.txt` or `.npy`.