LAL  7.2.4.1-d5ae9e8
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 six
35 import numpy as np
36 
37 
38 Attributes = ligolw.sax.xmlreader.AttributesImpl
39 
40 
41 #
42 # =============================================================================
43 #
44 # XML I/O
45 #
46 # =============================================================================
47 #
48 
49 
50 def _build_series(series, dim_names, comment, delta_name, delta_unit):
51  elem = ligolw.LIGO_LW(Attributes({u"Name": six.text_type(series.__class__.__name__)}))
52  if comment is not None:
53  elem.appendChild(ligolw.Comment()).pcdata = comment
54  elem.appendChild(ligolw.Time.from_gps(series.epoch, u"epoch"))
55  elem.appendChild(ligolw_param.Param.from_pyvalue(u"f0", series.f0, unit=u"s^-1"))
56  delta = getattr(series, delta_name)
57  if np.iscomplexobj(series.data.data):
58  data = np.row_stack((np.arange(len(series.data.data)) * delta, series.data.data.real, series.data.data.imag))
59  else:
60  data = np.row_stack((np.arange(len(series.data.data)) * delta, series.data.data))
61  a = ligolw_array.Array.build(series.name, data, dim_names=dim_names)
62  a.Unit = str(series.sampleUnits)
63  dim0 = a.getElementsByTagName(ligolw.Dim.tagName)[0]
64  dim0.Unit = delta_unit
65  dim0.Start = series.f0
66  dim0.Scale = delta
67  elem.appendChild(a)
68  return elem
69 
70 
71 def _parse_series(elem, creatorfunc, delta_target_unit_string):
72  t, = elem.getElementsByTagName(ligolw.Time.tagName)
73  a, = elem.getElementsByTagName(ligolw.Array.tagName)
74  dims = a.getElementsByTagName(ligolw.Dim.tagName)
75  f0 = ligolw_param.get_param(elem, u"f0")
76 
77  if t.Type != u"GPS":
78  raise ValueError("epoch Type must be GPS")
79  epoch = t.pcdata
80 
81  # Target units: inverse seconds
82  inverse_seconds_unit = lal.Unit("s^-1")
83 
84  delta_target_unit = lal.Unit(delta_target_unit_string)
85 
86  # Parse units of f0 field
87  f0_unit = lal.Unit(str(f0.Unit))
88 
89  # Parse units of deltaF field
90  delta_unit = lal.Unit(str(dims[0].Unit))
91 
92  # Parse units of data
93  sample_unit = lal.Unit(str(a.Unit))
94 
95  # Initialize data structure
96  series = creatorfunc(
97  str(a.Name),
98  epoch,
99  f0.pcdata * float(f0_unit / inverse_seconds_unit),
100  dims[0].Scale * float(delta_unit / delta_target_unit),
101  sample_unit,
102  len(a.array.T)
103  )
104 
105  # Assign data
106  if np.iscomplexobj(series.data.data):
107  series.data.data = a.array[1] + 1j * a.array[2]
108  else:
109  series.data.data = a.array[1]
110 
111  # Done!
112  return series
113 
114 
115 def build_REAL4FrequencySeries(series, comment=None):
116  assert isinstance(series, lal.REAL4FrequencySeries)
117  return _build_series(series, (u"Frequency,Real", u"Frequency"), comment, 'deltaF', 's^-1')
118 
119 
121  return _parse_series(elem, lal.CreateREAL4FrequencySeries, "s^-1")
122 
123 
124 def build_REAL8FrequencySeries(series, comment=None):
125  assert isinstance(series, lal.REAL8FrequencySeries)
126  return _build_series(series, (u"Frequency,Real", u"Frequency"), comment, 'deltaF', 's^-1')
127 
128 
130  return _parse_series(elem, lal.CreateREAL8FrequencySeries, "s^-1")
131 
132 
133 def build_COMPLEX8FrequencySeries(series, comment=None):
134  assert isinstance(series, lal.COMPLEX8FrequencySeries)
135  return _build_series(series, (u"Frequency,Real,Imaginary", u"Frequency"), comment, 'deltaF', 's^-1')
136 
137 
139  return _parse_series(elem, lal.CreateCOMPLEX8FrequencySeries, "s^-1")
140 
141 
142 def build_COMPLEX16FrequencySeries(series, comment=None):
143  assert isinstance(series, lal.COMPLEX16FrequencySeries)
144  return _build_series(series, (u"Frequency,Real,Imaginary", u"Frequency"), comment, 'deltaF', 's^-1')
145 
146 
148  return _parse_series(elem, lal.CreateCOMPLEX16FrequencySeries, "s^-1")
149 
150 
151 def build_REAL4TimeSeries(series, comment=None):
152  assert isinstance(series, lal.REAL4TimeSeries)
153  return _build_series(series, (u"Time,Real", u"Time"), comment, 'deltaT', 's')
154 
155 
157  return _parse_series(elem, lal.CreateREAL4TimeSeries, "s")
158 
159 
160 def build_REAL8TimeSeries(series, comment=None):
161  assert isinstance(series, lal.REAL8TimeSeries)
162  return _build_series(series, (u"Time,Real", u"Time"), comment, 'deltaT', 's')
163 
164 
166  return _parse_series(elem, lal.CreateREAL8TimeSeries, "s")
167 
168 
169 def build_COMPLEX8TimeSeries(series, comment=None):
170  assert isinstance(series, lal.COMPLEX8TimeSeries)
171  return _build_series(series, (u"Time,Real,Imaginary", u"Time"), comment, 'deltaT', 's')
172 
173 
175  return _parse_series(elem, lal.CreateCOMPLEX8TimeSeries, "s")
176 
177 
178 def build_COMPLEX16TimeSeries(series, comment=None):
179  assert isinstance(series, lal.COMPLEX16TimeSeries)
180  return _build_series(series, (u"Time,Real,Imaginary", u"Time"), comment, 'deltaT', 's')
181 
182 
184  return _parse_series(elem, lal.CreateCOMPLEX16TimeSeries, "s")
185 
186 
187 #
188 # =============================================================================
189 #
190 # XML PSD I/O
191 #
192 # =============================================================================
193 #
194 
195 
196 def make_psd_xmldoc(psddict, xmldoc = None, root_name = u"psd"):
197  """
198  Construct an XML document tree representation of a dictionary of
199  frequency series objects containing PSDs. See also read_psd_xmldoc()
200  for a function to parse the resulting XML documents.
201 
202  If xmldoc is None (the default), then a new XML document is created and
203  the PSD dictionary added to it inside a LIGO_LW element. If xmldoc is
204  not None then the PSD dictionary is appended to the children of that
205  element inside a new LIGO_LW element. In both cases, the LIGO_LW
206  element's Name attribute is set to root_name. This will be looked for
207  by read_psd_xmldoc() when parsing the PSD document.
208  """
209  if xmldoc is None:
210  xmldoc = ligolw.Document()
211  lw = xmldoc.appendChild(ligolw.LIGO_LW(Attributes({u"Name": root_name})))
212  for instrument, psd in psddict.items():
213  fs = lw.appendChild(build_REAL8FrequencySeries(psd))
214  if instrument is not None:
215  fs.appendChild(ligolw_param.Param.from_pyvalue(u"instrument", instrument))
216  return xmldoc
217 
218 
219 def read_psd_xmldoc(xmldoc, root_name = u"psd"):
220  """
221  Parse a dictionary of PSD frequency series objects from an XML
222  document. See also make_psd_xmldoc() for the construction of XML
223  documents from a dictionary of PSDs. Interprets an empty frequency
224  series for an instrument as None.
225 
226  The XML document tree is searched for a LIGO_LW element whose Name
227  attribute is root_name (default is "psd"). If root_name is None all
228  REAL8Frequency series objects below xmldoc are included in the return
229  value.
230  """
231  if root_name is not None:
232  xmldoc, = (elem for elem in xmldoc.getElementsByTagName(ligolw.LIGO_LW.tagName) if elem.hasAttribute(u"Name") and elem.Name == root_name)
233  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")
234  # interpret empty frequency series as None
235  for instrument in result:
236  if len(result[instrument].data.data) == 0:
237  result[instrument] = None
238  return result
239 
240 
241 @ligolw_array.use_in
242 @ligolw_param.use_in
243 class PSDContentHandler(ligolw.LIGOLWContentHandler):
244  """A content handler suitable for reading PSD documents. Use like this:
245 
246  >>> from ligo.lw.utils import load_filename
247  >>> xmldoc = load_filename('psd.xml', contenthandler=PSDContentHandler)
248  >>> psds = read_psd_xmldoc(xmldoc)
249  """
def read_psd_xmldoc(xmldoc, root_name=u"psd")
Parse a dictionary of PSD frequency series objects from an XML document.
Definition: series.py:230
def parse_REAL4TimeSeries(elem)
Definition: series.py:156
def build_COMPLEX8FrequencySeries(series, comment=None)
Definition: series.py:133
def build_COMPLEX8TimeSeries(series, comment=None)
Definition: series.py:169
def build_REAL4FrequencySeries(series, comment=None)
Definition: series.py:115
def parse_COMPLEX16FrequencySeries(elem)
Definition: series.py:147
def build_REAL8TimeSeries(series, comment=None)
Definition: series.py:160
Attributes
Definition: series.py:38
def parse_COMPLEX16TimeSeries(elem)
Definition: series.py:183
def build_REAL4TimeSeries(series, comment=None)
Definition: series.py:151
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:208
def parse_REAL8FrequencySeries(elem)
Definition: series.py:129
def parse_COMPLEX8TimeSeries(elem)
Definition: series.py:174
def parse_REAL4FrequencySeries(elem)
Definition: series.py:120
def build_REAL8FrequencySeries(series, comment=None)
Definition: series.py:124
def build_COMPLEX16TimeSeries(series, comment=None)
Definition: series.py:178
def parse_COMPLEX8FrequencySeries(elem)
Definition: series.py:138
def build_COMPLEX16FrequencySeries(series, comment=None)
Definition: series.py:142
def parse_REAL8TimeSeries(elem)
Definition: series.py:165