LAL  7.4.1.1-f1ed79c
series.py
Go to the documentation of this file.
1 # Copyright (C) 2008 Kipp Cannon
2 # 2015 Leo Singer
3 #
4 # Adapted from original pylal.series module to return SWIG lal datatypes
5 # instead of pylal datatypes.
6 #
7 # This program is free software; you can redistribute it and/or modify it
8 # under the terms of the GNU General Public License as published by the
9 # Free Software Foundation; either version 2 of the License, or (at your
10 # option) any later version.
11 #
12 # This program is distributed in the hope that it will be useful, but
13 # WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15 # Public License for more details.
16 #
17 # You should have received a copy of the GNU General Public License along
18 # with this program; if not, write to the Free Software Foundation, Inc.,
19 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
20 #
21 """
22 Code to assist in reading and writing LAL time- and frequency series data
23 encoded in LIGO Light-Weight XML format. The format recognized by the code
24 in this module is the same as generated by the array-related functions in
25 LAL's XML I/O code. The format is also very similar to the format used by
26 the DMT to store time- and frequency-series data in XML files,
27 """
28 
29 
30 from ligo.lw import ligolw
31 from ligo.lw import array as ligolw_array
32 from ligo.lw import param as ligolw_param
33 import lal
34 import numpy as np
35 
36 
37 Attributes = ligolw.sax.xmlreader.AttributesImpl
38 
39 
40 #
41 # =============================================================================
42 #
43 # XML I/O
44 #
45 # =============================================================================
46 #
47 
48 
49 def _build_series(series, dim_names, comment, delta_name, delta_unit):
50  elem = ligolw.LIGO_LW(Attributes({u"Name": str(series.__class__.__name__)}))
51  if comment is not None:
52  elem.appendChild(ligolw.Comment()).pcdata = comment
53  elem.appendChild(ligolw.Time.from_gps(series.epoch, u"epoch"))
54  elem.appendChild(ligolw_param.Param.from_pyvalue(u"f0", series.f0, unit=u"s^-1"))
55  delta = getattr(series, delta_name)
56  if np.iscomplexobj(series.data.data):
57  data = np.row_stack((np.arange(len(series.data.data)) * delta, series.data.data.real, series.data.data.imag))
58  else:
59  data = np.row_stack((np.arange(len(series.data.data)) * delta, series.data.data))
60  a = ligolw_array.Array.build(series.name, data, dim_names=dim_names)
61  a.Unit = str(series.sampleUnits)
62  dim0 = a.getElementsByTagName(ligolw.Dim.tagName)[0]
63  dim0.Unit = delta_unit
64  dim0.Start = series.f0
65  dim0.Scale = delta
66  elem.appendChild(a)
67  return elem
68 
69 
70 def _parse_series(elem, creatorfunc, delta_target_unit_string):
71  t, = elem.getElementsByTagName(ligolw.Time.tagName)
72  a, = elem.getElementsByTagName(ligolw.Array.tagName)
73  dims = a.getElementsByTagName(ligolw.Dim.tagName)
74  f0 = ligolw_param.get_param(elem, u"f0")
75 
76  if t.Type != u"GPS":
77  raise ValueError("epoch Type must be GPS")
78  epoch = t.pcdata
79 
80  # Target units: inverse seconds
81  inverse_seconds_unit = lal.Unit("s^-1")
82 
83  delta_target_unit = lal.Unit(delta_target_unit_string)
84 
85  # Parse units of f0 field
86  f0_unit = lal.Unit(str(f0.Unit))
87 
88  # Parse units of deltaF field
89  delta_unit = lal.Unit(str(dims[0].Unit))
90 
91  # Parse units of data
92  sample_unit = lal.Unit(str(a.Unit))
93 
94  # Initialize data structure
95  series = creatorfunc(
96  str(a.Name),
97  epoch,
98  f0.pcdata * float(f0_unit / inverse_seconds_unit),
99  dims[0].Scale * float(delta_unit / delta_target_unit),
100  sample_unit,
101  len(a.array.T)
102  )
103 
104  # Assign data
105  if np.iscomplexobj(series.data.data):
106  series.data.data = a.array[1] + 1j * a.array[2]
107  else:
108  series.data.data = a.array[1]
109 
110  # Done!
111  return series
112 
113 
114 def build_REAL4FrequencySeries(series, comment=None):
115  assert isinstance(series, lal.REAL4FrequencySeries)
116  return _build_series(series, (u"Frequency,Real", u"Frequency"), comment, 'deltaF', 's^-1')
117 
118 
120  return _parse_series(elem, lal.CreateREAL4FrequencySeries, "s^-1")
121 
122 
123 def build_REAL8FrequencySeries(series, comment=None):
124  assert isinstance(series, lal.REAL8FrequencySeries)
125  return _build_series(series, (u"Frequency,Real", u"Frequency"), comment, 'deltaF', 's^-1')
126 
127 
129  return _parse_series(elem, lal.CreateREAL8FrequencySeries, "s^-1")
130 
131 
132 def build_COMPLEX8FrequencySeries(series, comment=None):
133  assert isinstance(series, lal.COMPLEX8FrequencySeries)
134  return _build_series(series, (u"Frequency,Real,Imaginary", u"Frequency"), comment, 'deltaF', 's^-1')
135 
136 
138  return _parse_series(elem, lal.CreateCOMPLEX8FrequencySeries, "s^-1")
139 
140 
141 def build_COMPLEX16FrequencySeries(series, comment=None):
142  assert isinstance(series, lal.COMPLEX16FrequencySeries)
143  return _build_series(series, (u"Frequency,Real,Imaginary", u"Frequency"), comment, 'deltaF', 's^-1')
144 
145 
147  return _parse_series(elem, lal.CreateCOMPLEX16FrequencySeries, "s^-1")
148 
149 
150 def build_REAL4TimeSeries(series, comment=None):
151  assert isinstance(series, lal.REAL4TimeSeries)
152  return _build_series(series, (u"Time,Real", u"Time"), comment, 'deltaT', 's')
153 
154 
156  return _parse_series(elem, lal.CreateREAL4TimeSeries, "s")
157 
158 
159 def build_REAL8TimeSeries(series, comment=None):
160  assert isinstance(series, lal.REAL8TimeSeries)
161  return _build_series(series, (u"Time,Real", u"Time"), comment, 'deltaT', 's')
162 
163 
165  return _parse_series(elem, lal.CreateREAL8TimeSeries, "s")
166 
167 
168 def build_COMPLEX8TimeSeries(series, comment=None):
169  assert isinstance(series, lal.COMPLEX8TimeSeries)
170  return _build_series(series, (u"Time,Real,Imaginary", u"Time"), comment, 'deltaT', 's')
171 
172 
174  return _parse_series(elem, lal.CreateCOMPLEX8TimeSeries, "s")
175 
176 
177 def build_COMPLEX16TimeSeries(series, comment=None):
178  assert isinstance(series, lal.COMPLEX16TimeSeries)
179  return _build_series(series, (u"Time,Real,Imaginary", u"Time"), comment, 'deltaT', 's')
180 
181 
183  return _parse_series(elem, lal.CreateCOMPLEX16TimeSeries, "s")
184 
185 
186 #
187 # =============================================================================
188 #
189 # XML PSD I/O
190 #
191 # =============================================================================
192 #
193 
194 
195 def make_psd_xmldoc(psddict, xmldoc = None, root_name = u"psd"):
196  """
197  Construct an XML document tree representation of a dictionary of
198  frequency series objects containing PSDs. See also read_psd_xmldoc()
199  for a function to parse the resulting XML documents.
200 
201  If xmldoc is None (the default), then a new XML document is created and
202  the PSD dictionary added to it inside a LIGO_LW element. If xmldoc is
203  not None then the PSD dictionary is appended to the children of that
204  element inside a new LIGO_LW element. In both cases, the LIGO_LW
205  element's Name attribute is set to root_name. This will be looked for
206  by read_psd_xmldoc() when parsing the PSD document.
207  """
208  if xmldoc is None:
209  xmldoc = ligolw.Document()
210  lw = xmldoc.appendChild(ligolw.LIGO_LW(Attributes({u"Name": root_name})))
211  for instrument, psd in psddict.items():
212  fs = lw.appendChild(build_REAL8FrequencySeries(psd))
213  if instrument is not None:
214  fs.appendChild(ligolw_param.Param.from_pyvalue(u"instrument", instrument))
215  return xmldoc
216 
217 
218 def read_psd_xmldoc(xmldoc, root_name = u"psd"):
219  """
220  Parse a dictionary of PSD frequency series objects from an XML
221  document. See also make_psd_xmldoc() for the construction of XML
222  documents from a dictionary of PSDs. Interprets an empty frequency
223  series for an instrument as None.
224 
225  The XML document tree is searched for a LIGO_LW element whose Name
226  attribute is root_name (default is "psd"). If root_name is None all
227  REAL8Frequency series objects below xmldoc are included in the return
228  value.
229  """
230  if root_name is not None:
231  xmldoc, = (elem for elem in xmldoc.getElementsByTagName(ligolw.LIGO_LW.tagName) if elem.hasAttribute(u"Name") and elem.Name == root_name)
232  result = dict((ligolw_param.get_pyvalue(elem, u"instrument"), parse_REAL8FrequencySeries(elem)) for elem in xmldoc.getElementsByTagName(ligolw.LIGO_LW.tagName) if elem.hasAttribute(u"Name") and elem.Name == u"REAL8FrequencySeries")
233  # interpret empty frequency series as None
234  for instrument in result:
235  if len(result[instrument].data.data) == 0:
236  result[instrument] = None
237  return result
238 
239 
240 @ligolw_array.use_in
241 @ligolw_param.use_in
242 class PSDContentHandler(ligolw.LIGOLWContentHandler):
243  """A content handler suitable for reading PSD documents. Use like this:
244 
245  >>> from ligo.lw.utils import load_filename
246  >>> xmldoc = load_filename('psd.xml', contenthandler=PSDContentHandler)
247  >>> psds = read_psd_xmldoc(xmldoc)
248  """
def read_psd_xmldoc(xmldoc, root_name=u"psd")
Parse a dictionary of PSD frequency series objects from an XML document.
Definition: series.py:229
def parse_REAL4TimeSeries(elem)
Definition: series.py:155
def build_COMPLEX8FrequencySeries(series, comment=None)
Definition: series.py:132
def build_COMPLEX8TimeSeries(series, comment=None)
Definition: series.py:168
def build_REAL4FrequencySeries(series, comment=None)
Definition: series.py:114
def parse_COMPLEX16FrequencySeries(elem)
Definition: series.py:146
def build_REAL8TimeSeries(series, comment=None)
Definition: series.py:159
Attributes
Definition: series.py:37
def parse_COMPLEX16TimeSeries(elem)
Definition: series.py:182
def build_REAL4TimeSeries(series, comment=None)
Definition: series.py:150
def make_psd_xmldoc(psddict, xmldoc=None, root_name=u"psd")
Construct an XML document tree representation of a dictionary of frequency series objects containing ...
Definition: series.py:207
def parse_REAL8FrequencySeries(elem)
Definition: series.py:128
def parse_COMPLEX8TimeSeries(elem)
Definition: series.py:173
def parse_REAL4FrequencySeries(elem)
Definition: series.py:119
def build_REAL8FrequencySeries(series, comment=None)
Definition: series.py:123
def build_COMPLEX16TimeSeries(series, comment=None)
Definition: series.py:177
def parse_COMPLEX8FrequencySeries(elem)
Definition: series.py:137
def build_COMPLEX16FrequencySeries(series, comment=None)
Definition: series.py:141
def parse_REAL8TimeSeries(elem)
Definition: series.py:164