diff --git a/bin/live/pycbc_live_collated_dq_trigger_rates b/bin/live/pycbc_live_collated_dq_trigger_rates index 2a9d0982771..0086274950e 100644 --- a/bin/live/pycbc_live_collated_dq_trigger_rates +++ b/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') \ No newline at end of file diff --git a/companion.txt b/companion.txt index 06645ccbf10..0c77ceba7db 100644 --- a/companion.txt +++ b/companion.txt @@ -1,6 +1,3 @@ -# other tools which may be useful -gwpy>=0.8.1 - # HEALPix is very useful for some analysis. healpy diff --git a/pycbc/types/timeseries.py b/pycbc/types/timeseries.py index f81e26a58d4..38574558629 100644 --- a/pycbc/types/timeseries.py +++ b/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`.