Merge branch 'release/0.10.0'

examples/frequencyseries/index.rst

@@ -14,3 +14,4 @@
    coherence
    transfer_function
    rayleigh
+   inject

examples/frequencyseries/inject.py

@@ -0,0 +1,72 @@
+#!/usr/bin/env python
+
+
+# Copyright (C) Duncan Macleod (2013)
+#
+# This file is part of GWpy.
+#
+# GWpy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# GWpy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with GWpy.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Inject a known signal into a `FrequencySeries`
+
+It can often be useful to add some known signal to inherently random
+or noisy data. For example, one might want to investigate what
+would happen if a binary black hole merger signal occured at or near
+the time of a glitch. In LIGO data analysis, this procedure is referred
+to as an _injection_.
+
+In the example below we will create a stream of random, white Gaussian
+noise, then inject a loud, steady sinuosoid. We will do this in the
+frequency domain because it is much easier to model a sinusoid there.
+"""
+
+__author__ = "Alex Urban <alexander.urban@ligo.org>"
+__currentmodule__ = 'gwpy.timeseries'
+
+# First, we prepare one second of Gaussian noise:
+
+from numpy import random
+from gwpy.timeseries import TimeSeries
+noise = TimeSeries(random.normal(scale=.1, size=1024), sample_rate=1024)
+
+# To inject a signal in the frequency domain, we need to take an FFT:
+
+noisefd = noise.fft()
+
+# We can now easily inject a loud sinusoid of unit amplitude at, say,
+# 30 Hz. To do this, we use :meth:`~gwpy.types.series.Series.inject`.
+
+import numpy
+from gwpy.frequencyseries import FrequencySeries
+signal = FrequencySeries(numpy.array([1.]), f0=30, df=noisefd.df)
+injfd = noisefd.inject(signal)
+
+# We can then visualize the data before and after injection in the frequency
+# domain:
+
+from gwpy.plotter import FrequencySeriesPlot
+plot = FrequencySeriesPlot(numpy.abs(noisefd), numpy.abs(injfd), sep=True,
+                           sharex=True, sharey=True)
+plot.show()
+
+# Finally, for completeness we can visualize the effect before and after
+# injection back in the time domain:
+
+from gwpy.plotter import TimeSeriesPlot
+inj = injfd.ifft()
+plot = TimeSeriesPlot(noise, inj, sep=True, sharex=True, sharey=True)
+plot.show()
+
+# We can see why sinusoids are easier to inject in the frequency domain:
+# they only require adding at a single frequency.

examples/timeseries/index.rst

@@ -15,3 +15,4 @@
    statevector
    qscan
    pycbc-snr
+   inject

examples/timeseries/inject.py

@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+
+
+# Copyright (C) Duncan Macleod (2013)
+#
+# This file is part of GWpy.
+#
+# GWpy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# GWpy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with GWpy.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Inject a known signal into a `TimeSeries`
+
+It can often be useful to add some known signal to an inherently random
+or noisy timeseries. For example, one might want to examine what
+would happen if a binary black hole merger signal occured at or near
+the time of a glitch. In LIGO data analysis, this procedure is referred
+to as an _injection_.
+
+In the example below, we will create a stream of random, white Gaussian
+noise, then inject a simulation of GW150914 into it at a known time.
+"""
+
+__author__ = "Alex Urban <alexander.urban@ligo.org>"
+__currentmodule__ = 'gwpy.timeseries'
+
+# First, we prepare one second of Gaussian noise:
+
+from numpy import random
+from gwpy.timeseries import TimeSeries
+noise = TimeSeries(random.normal(scale=.1, size=16384), sample_rate=16384)
+
+# Then we can download a simulation of the GW150914 signal from LOSC:
+
+from astropy.utils.data import get_readable_fileobj
+source = 'https://losc.ligo.org/s/events/GW150914/P150914/'
+url = '%s/fig2-unfiltered-waveform-H.txt' % source
+with get_readable_fileobj(url) as f:
+    signal = TimeSeries.read(f, format='txt')
+signal.t0 = .5  # make sure this intersects with noise time samples
+
+# Note, since this simulation cuts off before a certain time, it is
+# important to taper its ends to zero to avoid ringing artifacts.
+# We can accomplish this using the
+# :meth:`~gwpy.timeseries.TimeSeries.taper` method.
+
+signal = signal.taper()
+
+# Since the time samples overlap, we can inject this into our noise data
+# using :meth:`~gwpy.types.series.Series.inject`:
+
+data = noise.inject(signal)
+
+# Finally, we can visualize the full process in the time domain:
+
+from gwpy.plotter import TimeSeriesPlot
+plot = TimeSeriesPlot(noise, signal, data, sep=True, sharex=True, sharey=True)
+plot.set_epoch(0)
+plot.show()
+
+# We can clearly see that the loud GW150914-like signal has been layered
+# on top of Gaussian noise with the correct amplitude and phase evolution.

gwpy/frequencyseries/core.py

@@ -230,10 +230,9 @@ def ifft(self):
         # Undo normalization from TimeSeries.fft
         # The DC component does not have the factor of two applied
         # so we account for it here
-        dift = npfft.irfft(self.value * nout)
-        dift[1:] /= 2
+        dift = npfft.irfft(self.value * nout) / 2
         new = TimeSeries(dift, epoch=self.epoch, channel=self.channel,
-                         unit=self.unit * units.Hertz, dx=1/self.dx/nout)
+                         unit=self.unit, dx=1/self.dx/nout)
         return new
 
     def zpk(self, zeros, poles, gain, analog=True):

gwpy/plotter/histogram.py

@@ -186,23 +186,20 @@ class HistogramPlot(Plot):
     def __init__(self, *data, **kwargs):
         """Generate a new `HistogramPlot` from some ``data``.
         """
-        # extract histogram arguments
-        histargs = dict()
-        for key in ['bins', 'range', 'normed', 'weights', 'cumulative',
-                    'bottom', 'histtype', 'align', 'orientation', 'rwidth',
-                    'log', 'color', 'label', 'stacked', 'logbins']:
-            try:
-                histargs[key] = kwargs.pop(key)
-            except KeyError:
-                pass
+        # separate keyword arguments
+        axargs, histargs = self._parse_kwargs(kwargs)
+
         # generate Figure
         super(HistogramPlot, self).__init__(**kwargs)
-        # plot data
+
+        # create axes
         if data:
-            ax = self.gca()
+            ax = self.gca(**axargs)
             data = list(data)
         else:
             ax = None
+
+        # plot data
         while data:
             dataset = data.pop(0)
             if isinstance(dataset, Series):
@@ -212,6 +209,7 @@ def __init__(self, *data, **kwargs):
                 ax.hist_table(dataset, column, **histargs)
             else:
                 ax.hist(dataset, **histargs)
+
         if ax and histargs.get('logbins', False):
             if histargs.get('orientation', 'vertical') == 'vertical':
                 ax.set_xscale('log')

gwpy/plotter/series.py

@@ -121,7 +121,7 @@ def plot_series(self, series, **kwargs):
         return line
 
     @auto_refresh
-    def plot_mmm(self, mean_, min_=None, max_=None, **kwargs):
+    def plot_mmm(self, mean_, min_=None, max_=None, alpha=.1, **kwargs):
         """Plot a `Series` onto these axes, with shaded regions
 
         The ``mean_`` `Series` is plotted normally, while the ``min_``
@@ -139,9 +139,13 @@ def plot_mmm(self, mean_, min_=None, max_=None, **kwargs):
         max_ : `~gwpy.types.Series`
             second data set to shade to ``mean_``
 
+        alpha : `float`, `None`, optional
+            value for alpha channel for shading, default: ``0.1``
+
         **kwargs
             any other keyword arguments acceptable for
-            :meth:`~matplotlib.Axes.plot`
+            :meth:`~matplotlib.Axes.plot`, only used to draw the
+            lines
 
         Returns
         -------
@@ -163,15 +167,17 @@ def plot_mmm(self, mean_, min_=None, max_=None, **kwargs):
         meanline = self.plot_series(mean_, **kwargs)[0]
 
         # modify keywords for shading
-        kwargs.pop('label', None)
+        kwargs.update({
+            'label': '',
+        })
         color = kwargs.pop('color', meanline.get_color())
         linewidth = kwargs.pop('linewidth', meanline.get_linewidth()) / 2
 
         def _plot_shade(series):
-            line = self.plot(series, color=color, linewidth=linewidth,
-                             **kwargs)
+            line, = self.plot(series, color=color, linewidth=linewidth,
+                              **kwargs)
             coll = self.fill_between(series.xindex.value, series.value,
-                                     mean_.value, alpha=.1, color=color,
+                                     mean_.value, alpha=alpha, color=color,
                                      rasterized=kwargs.get('rasterized'))
             return line, coll
 

gwpy/plotter/table.py

@@ -25,18 +25,18 @@
 from __future__ import division
 
 import re
+from collections import OrderedDict
+
 from six import string_types
 
 import numpy
 from matplotlib import (collections, pyplot)
 from matplotlib.projections import register_projection
 
 from ..table import Table
-from ..time import LIGOTimeGPS
 from .core import Plot
 from .timeseries import (TimeSeriesAxes, TimeSeriesPlot)
 from .frequencyseries import FrequencySeriesPlot
-from .tex import float_to_latex
 
 __all__ = ['EventTableAxes', 'EventTablePlot']
 
@@ -279,52 +279,55 @@ def add_loudest(self, table, rank, x, y, *columns, **kwargs):
         y : `str`
             name of column to display on the Y-axis
 
-        color : `str`, optional
-            name of column by which to colour the data
-
         **kwargs
             any other arguments applicable to
             :meth:`~matplotlib.axes.Axes.text`
 
         Returns
         -------
-        out : `tuple`
-            (`collection`, `text`) tuple of items added to the `Axes`
+        collection : `~matplotlib.collections.PolyCollection`
+            the collection returned by :meth:`~matplotlib.axes.Axesscatter`
+            when marking the loudest event
+        text : `~matplotlib.text.Text`
+            the text added to these `Axes` with loudest event parameters
         """
-        ylim = self.get_ylim()
+        # get loudest row
         idx = table[rank].argmax()
         row = table[idx]
-        disp = "Loudest event:"
+
+        # mark loudest row with star
+        coll = self.scatter([float(row[x])], [float(row[y])],
+                            marker='*', zorder=1000, facecolor='gold',
+                            edgecolor='black', s=200)
+
+        # get text
         columns = [x, y, rank] + list(columns)
-        scat = []
-        for i, column in enumerate(columns):
-            if not column or column in columns[:i]:
-                continue
-            if i:
-                disp += ','
-            val = row[column]
-            if i < 2:
-                scat.append([float(val)])
-            column = get_column_string(column)
-            if pyplot.rcParams['text.usetex'] and column.endswith('Time'):
-                disp += (r" %s$= %s$" % (column, LIGOTimeGPS(float(val))))
-            elif pyplot.rcParams['text.usetex']:
-                disp += (r" %s$=$ %s" % (column, float_to_latex(val, '%.3g')))
-            else:
-                disp += " %s = %.2g" % (column, val)
-        disp = disp.rstrip(',')
-        pos = kwargs.pop('position', [0.5, 1.00])
-        kwargs.setdefault('transform', self.axes.transAxes)
-        kwargs.setdefault('verticalalignment', 'bottom')
-        kwargs.setdefault('horizontalalignment', 'center')
-        args = pos + [disp]
-        self.scatter(*scat, marker='*', zorder=1000, facecolor='gold',
-                     edgecolor='black', s=200)
-        self.text(*args, **kwargs)
-        if self.get_title():
-            pos = self.title.get_position()
-            self.title.set_position((pos[0], pos[1] + 0.05))
-        self.set_ylim(*ylim)
+        loudtext = _loudest_text(row, columns)
+
+        # get position for new text
+        try:
+            pos = kwargs.pop('position')
+        except KeyError:  # user didn't specify, set default and shunt title
+            pos = [0.5, 1.00]
+            tpos = self.title.get_position()
+            self.title.set_position((tpos[0], tpos[1] + 0.05))
+
+        # parse text kwargs
+        text_kw = {  # defaults
+            'transform': self.axes.transAxes,
+            'verticalalignment': 'bottom',
+            'horizontalalignment': 'center',
+        }
+        text_kw.update(kwargs)
+        if 'ha' in text_kw:  # handle short versions or alignment params
+            text_kw['horizontalalignment'] = text_kw.pop('ha')
+        if 'va' in text_kw:
+            text_kw['verticalalignment'] = text_kw.pop('va')
+
+        # add text
+        text = self.text(pos[0], pos[1], loudtext, **text_kw)
+
+        return coll, text
 
 
 register_projection(EventTableAxes)
@@ -665,3 +668,22 @@ def get_column_string(column):
             # escape underscore
             words[i] = re.sub(r'(?<!\\)_', r'\_', words[i])
     return ' '.join(words)
+
+
+def _loudest_text(row, columns):
+    """Format the text for `EventTableAxes.add_loudest`
+    """
+    coltxt = []
+    for col in OrderedDict.fromkeys(columns):  # unique list of cols
+        # format column name
+        colstr = get_column_string(col)
+
+        # format row value
+        try:
+            valstr = '{0:.2f}'.format(row[col]).rstrip('.0')
+        except ValueError:  # not float()able
+            valstr = str(row[col])
+
+        coltxt.append('{col} = {val}'.format(col=colstr, val=valstr))
+
+    return 'Loudest event: {0}'.format(', '.join(coltxt))