Merge 6418c95 into 106e7be

bin/gwdetchar-omega

@@ -32,6 +32,7 @@ use('agg')  # nopep8
 
 from gwpy.utils import gprint
 from gwpy.time import tconvert
+from gwpy.segments import Segment
 from gwpy.table import EventTable
 from gwpy.timeseries import TimeSeriesDict
 from gwpy.detector import (Channel, ChannelList)
@@ -187,7 +188,7 @@ def get_widths(x0, xdata):
         yield width
 
 
-def eventgram(time, data, search=0.5, frange=(0, numpy.inf),
+def eventgram(time, data, search=0.25, frange=(0, numpy.inf),
               qrange=(4, 96), snrthresh=5.5, mismatch=0.2):
     """Create an eventgram with the Q-plane that has the most significant
     tile.
@@ -304,11 +305,11 @@ gprint('Launching Omega scans...')
 for block in blocks[:]:
     gprint('Processing block %s' % block.name)
     chans = [c.name for c in block.channels]
-    # read in fftlength seconds of data
-    # centered on gps
+    # read in `duration` seconds of data, centered on gps
     duration = block.duration
     fftlength = block.fftlength
-    data = TimeSeriesDict.get(chans, gps-256-fftlength/4, gps+256+fftlength/4,
+    pad = max(fftlength/4, 1)
+    data = TimeSeriesDict.get(chans, gps-duration/2-pad, gps+duration/2+pad,
                               frametype=block.frametype, nproc=args.nproc,
                               verbose=args.verbose)
     # compute qscans
@@ -321,61 +322,54 @@ for block in blocks[:]:
         if block.resample:
             series = series.resample(block.resample)
 
-        # filter the timeseries
+        # highpass and whiten the timeseries
         corner = c.frange[0] / 1.5
         wn = 2 * corner * series.dt.decompose().value
         hpfilt = butter(12, wn, btype='highpass', analog=False, output='sos')
-        hpseries = series.filter(hpfilt, filtfilt=True)
-        asd = series.asd(fftlength, fftlength/2, method='lal_median_mean')
-        wseries = hpseries.whiten(fftlength, fftlength/2, window='hann',
-                                  asd=asd, detrend='linear')
+        series = series.filter(hpfilt, filtfilt=True)
+        series = series.whiten(fftlength=fftlength, overlap=fftlength/2,
+                               window='hann', method='lal_median_mean')
 
-        # crop the timeseries
-        wseries = wseries.crop(gps-duration/2, gps+duration/2)
-        hpseries = hpseries.crop(gps-duration/2, gps+duration/2)
+        # crop filter-corrupted boundaries
+        series = series.crop(gps-duration/2, gps+duration/2)
 
-        # compute eventgrams
+        # compute eventgram
         try:
-            table = eventgram(gps, wseries, frange=c.frange, qrange=c.qrange,
+            table = eventgram(gps, series, frange=c.frange, qrange=c.qrange,
                               snrthresh=c.snrthresh, mismatch=c.mismatch)
         except UnboundLocalError:
             if args.verbose:
                 gprint('Channel is misbehaved, removing it from the analysis')
-            del series, hpseries, wseries, asd
+            del series
             block.channels.remove(c)
             continue
+
+        # if insignificant, skip this analysis
         if table.Z < table.engthresh and not c.always_plot:
             if args.verbose:
                 gprint('Channel not significant at white noise false alarm '
                        'rate %s Hz' % far)
-            del series, hpseries, wseries, asd, table
+            del series
             block.channels.remove(c)
             continue
         Q = table.q
-        rtable = eventgram(gps, hpseries, frange=table.frange, qrange=(Q, Q),
-                           snrthresh=c.snrthresh, mismatch=c.mismatch)
 
-        # compute Q-transforms
+        # compute Q-transform spectrograms
         tres = min(c.pranges) / 500
-        fres = c.frange[0] / 5
-        qscan = wseries.q_transform(qrange=(Q, Q), frange=c.frange,
-                                    tres=tres, fres=fres, gps=gps,
-                                    search=0.25, whiten=False)
-        rqscan = hpseries.q_transform(qrange=(Q, Q), frange=c.frange,
-                                      tres=tres, fres=fres, gps=gps,
-                                      search=0.25, whiten=False)
+        fres = min(1. / fftlength, c.frange[0] / 32)
+        outseg = Segment(gps - max(c.pranges)/2, gps + max(c.pranges)/2)
+        qscan = series.q_transform(qrange=(Q, Q), frange=c.frange, tres=tres,
+                                   fres=fres, gps=gps, search=0.25,
+                                   whiten=False, outseg=outseg)
 
         # prepare plots
         if args.verbose:
-            gprint('Plotting omega scans for channel %s...' % c.name)
-        # work out figure size
+            gprint('Plotting omega scans for whitened channel %s...' % c.name)
         figsize = [8, 5]
-        for span, png1, png2, png3, png4, png5, png6, png7, png8, png9 in zip(
-            c.pranges, c.plots['qscan_whitened'],
-            c.plots['qscan_autoscaled'], c.plots['qscan_raw'],
-            c.plots['timeseries_raw'], c.plots['timeseries_highpassed'],
-            c.plots['timeseries_whitened'], c.plots['eventgram_raw'],
-            c.plots['eventgram_whitened'], c.plots['eventgram_autoscaled']
+        for span, png1, png2, png3, png4, png5 in zip(
+            c.pranges, c.plots['qscan_whitened'], c.plots['qscan_autoscaled'],
+            c.plots['timeseries_whitened'], c.plots['eventgram_whitened'],
+            c.plots['eventgram_autoscaled']
         ):
             # plot whitened qscan
             plot.omega_plot(qscan, gps, span, c.name, str(png1), qscan=True,
@@ -384,29 +378,15 @@ for block in blocks[:]:
             # plot autoscaled, whitened qscan
             plot.omega_plot(qscan, gps, span, c.name, str(png2), qscan=True,
                             colormap=args.colormap, figsize=figsize)
-            # plot raw qscan
-            plot.omega_plot(rqscan, gps, span, c.name, str(png3), qscan=True,
-                            clim=(0, 25), colormap=args.colormap,
-                            figsize=figsize)
-            # plot raw timeseries
-            plot.omega_plot(series, gps, span, c.name, str(png4),
-                            ylabel='Amplitude', figsize=figsize)
-            # plot highpassed timeseries
-            plot.omega_plot(hpseries, gps, span, c.name, str(png5),
-                            ylabel='Highpassed Amplitude', figsize=figsize)
             # plot whitened timeseries
-            plot.omega_plot(wseries, gps, span, c.name, str(png6),
+            plot.omega_plot(series, gps, span, c.name, str(png3),
                             ylabel='Whitened Amplitude', figsize=figsize)
-            # plot raw eventgram
-            plot.omega_plot(rtable, gps, span, c.name, str(png7),
-                            eventgram=True, clim=(0, 25),
-                            colormap=args.colormap, figsize=figsize)
             # plot whitened eventgram
-            plot.omega_plot(table, gps, span, c.name, str(png8),
+            plot.omega_plot(table, gps, span, c.name, str(png4),
                             eventgram=True, clim=(0, 25),
                             colormap=args.colormap, figsize=figsize)
             # plot autoscaled whitened eventgram
-            plot.omega_plot(table, gps, span, c.name, str(png9),
+            plot.omega_plot(table, gps, span, c.name, str(png5),
                             eventgram=True, colormap=args.colormap,
                             figsize=figsize)
 
@@ -417,11 +397,44 @@ for block in blocks[:]:
         c.t = table.tc
         c.f = table.fc
 
+        # repeat the analysis with highpassed data
+        if args.verbose:
+            gprint('Repeating for raw data from channel %s...' % c.name)
+        series = data[c.name].astype('float64')
+        if block.resample:
+            series = series.resample(block.resample)
+        for span, png in zip(c.pranges, c.plots['timeseries_raw']):
+            plot.omega_plot(series, gps, span, c.name, str(png),
+                            ylabel='Amplitude', figsize=figsize)
+        series = series.filter(hpfilt, filtfilt=True).crop(gps-duration/2,
+                                                           gps+duration/2)
+        table = eventgram(gps, series, frange=table.frange, qrange=(Q, Q),
+                          snrthresh=c.snrthresh, mismatch=c.mismatch)
+        qscan = series.q_transform(qrange=(Q, Q), frange=c.frange, tres=tres,
+                                   fres=fres, gps=gps, search=0.25,
+                                   whiten=False, outseg=outseg)
+        for span, png1, png2, png3 in zip(
+            c.pranges, c.plots['qscan_raw'], c.plots['timeseries_highpassed'],
+            c.plots['eventgram_raw']
+        ):
+            # plot raw qscan
+            plot.omega_plot(qscan, gps, span, c.name, str(png1), qscan=True,
+                            clim=(0, 25), colormap=args.colormap,
+                            figsize=figsize)
+            # plot highpassed timeseries
+            plot.omega_plot(series, gps, span, c.name, str(png2),
+                            ylabel='Highpassed Amplitude', figsize=figsize)
+            # plot raw eventgram
+            plot.omega_plot(table, gps, span, c.name, str(png3),
+                            eventgram=True, clim=(0, 25),
+                            colormap=args.colormap, figsize=figsize)
+
         # update HTML output
         html.write_qscan_page(ifo, gps, blocks, **htmlv)
 
         # delete intermediate data products
-        del qscan, rqscan, table, rtable, series, hpseries, wseries, asd
+        del data[c.name]
+        del (qscan, table, series)
 
     # delete data
     del data

gwdetchar/omega/html.py

@@ -35,7 +35,7 @@
 
 from glue import markup
 from gwpy.time import tconvert
-from gwpy.plotter.colors import GW_OBSERVATORY_COLORS
+from gwpy.plot.colors import GW_OBSERVATORY_COLORS
 from ..io.html import (JQUERY_JS, BOOTSTRAP_CSS, BOOTSTRAP_JS)
 from .. import __version__
 

gwdetchar/omega/plot.py

@@ -24,7 +24,7 @@
 from matplotlib import cm
 from matplotlib import rcParams
 
-from gwpy.plotter.colors import GW_OBSERVATORY_COLORS
+from gwpy.plot.colors import GW_OBSERVATORY_COLORS
 
 __author__ = 'Alex Urban <alexander.urban@ligo.org>'
 __credits__ = 'Duncan Macleod <duncan.macleod@ligo.org>'
@@ -51,12 +51,11 @@ def omega_plot(series, gps, span, channel, output, colormap='viridis',
     # construct plot
     if qscan:
         # plot Q-transform
-        plot = series.crop(gps-span/2, gps+span/2).plot(figsize=figsize)
+        plot = series.crop(gps-span/2, gps+span/2).imshow(figsize=figsize)
     elif eventgram:
         # plot eventgram
-        plot = series.plot('central_time', 'central_freq', 'duration',
-                           'bandwidth', color='energy',
-                           figsize=figsize)
+        plot = series.tile('central_time', 'central_freq', 'duration',
+                           'bandwidth', color='energy', figsize=figsize)
     else:
         # set color by IFO
         ifo = channel[:2]
@@ -65,19 +64,14 @@ def omega_plot(series, gps, span, channel, output, colormap='viridis',
                            figsize=figsize)
     ax = plot.gca()
     # set time axis units
-    ax.set_epoch(gps)
-    ax.set_xscale('auto-gps')
-    if span <= 1.:
-        ax.set_xlabel('Time [milliseconds]')
-    else:
-        ax.set_xlabel('Time [seconds]')
+    ax.set_xscale('auto-gps', epoch=gps)
+    ax.set_xlabel('Time [seconds]')
     # set y-axis properties
     chan = channel.replace('_', r'\_')
     if (qscan or eventgram):
         ax.set_yscale('log')
         ax.set_title('%s at %.3f with $Q=%.1f$' % (chan, gps, series.q))
-        plot.add_colorbar(cmap=colormap, clim=clim,
-                          label='Normalized energy')
+        ax.colorbar(cmap=colormap, clim=clim, label='Normalized energy')
     else:
         ax.set_yscale('linear')
         ax.set_title('%s at %.3f' % (chan, gps))

requirements.txt

@@ -8,7 +8,7 @@ numpy >= 1.10
 scipy >= 0.16
 matplotlib >= 2.0.0
 astropy >= 1.2
-gwpy >= 0.7
+gwpy >= 0.12.0
 lalsuite
 lscsoft-glue
 dqsegdb

setup.py

@@ -58,7 +58,7 @@
     'scipy>=0.16',
     'matplotlib>=2.0.0',
     'astropy>=1.2',
-    'gwpy>=0.5',
+    'gwpy>=0.12.0',
     'lscsoft-glue',
     'gwtrigfind',
 ]