1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 """
19 Auxiliary functions for running and postprocessing Omega pipeline code.
20 """
21
22 from __future__ import division
23
24 import os
25 import sys
26 import numpy
27 import re
28
29 from socket import getfqdn
30 from lal import LIGOTimeGPS, lalStrainUnit
31
32 from pylal import git_version
33 from glue.ligolw import lsctables
34 from glue import segments
35 from glue.lal import Cache,CacheEntry
36
37 __author__ = 'Duncan Macleod <duncan.macleod@ligo.org>'
38 __version__ = git_version.id
39 __date__ = git_version.date
40
41
42
43
44
45 _comment = re.compile('[#%]')
46 _delim = re.compile('[\t\,\s]+')
47
48 -def trigger(line, columns=lsctables.SnglBurst.__slots__, virgo=False, ifo=None, channel=None):
49 """
50 Convert a line from an Omega-format ASCII file into a SnglBurst object.
51 """
52
53 if isinstance(line, str):
54 dat = map(float, _delim.split(line.rstrip()))
55 else:
56 dat = map(float, line)
57
58
59 if virgo:
60 start, stop, peak, freq, band, cln, cle, snr = dat
61 start = LIGOTimeGPS(peak-duration/2)
62 stop = LIGOTimeGPS(peak+duration/2)
63 peak = LIGOTimeGPS(peak)
64 av_freq = freq
65 av_band = band
66 err_freq = 0
67 clusters = False
68 elif len(dat)==11:
69 peak, freq, duration, band, amplitude, cls, cle, cln, av_freq,\
70 av_band, err_freq = dat
71 start = LIGOTimeGPS(peak-duration/2)
72 stop = LIGOTimeGPS(peak+duration/2)
73 peak = LIGOTimeGPS(peak)
74 snr = (2*amplitude)**(1/2)
75 clusters = True
76 elif len(dat)==8:
77 peak, freq, duration, band, amplitude, cls, cle, cln = dat
78 start = LIGOTimeGPS(peak-duration/2)
79 stop = LIGOTimeGPS(peak+duration/2)
80 peak = LIGOTimeGPS(peak)
81 av_freq = freq
82 av_band = band
83 err_freq = 0
84 snr = (2*amplitude)**(1/2)
85 clusters = True
86 elif len(dat)==5:
87 peak, freq, duration, band, amplitude = dat
88 start = LIGOTimeGPS(peak-duration/2)
89 stop = LIGOTimeGPS(peak+duration/2)
90 av_freq = freq
91 av_band = band
92 err_freq = 0
93 snr = (2*amplitude)**(1/2)
94 clusters = False
95 elif len(dat)==4:
96 peak, av_freq, amplitude, chisq = dat
97 snr = (amplitude)**(1/2)
98 peak = LIGOTimeGPS(peak)
99 else:
100 raise ValueError("Wrong number of columns in ASCII line. Cannot read.")
101
102
103 t = lsctables.SnglBurst()
104
105
106 if 'ifo' in columns: t.ifo = ifo
107 if 'channel' in columns: t.channel = channel
108 if 'search' in columns: t.search = 'omega'
109
110
111 if 'start_time' in columns: t.start_time = start.gpsSeconds
112 if 'start_time_ns' in columns: t.start_time_ns = start.gpsNanoSeconds
113 if 'peak_time' in columns: t.peak_time = peak.gpsSeconds
114 if 'peak_time_ns' in columns: t.peak_time_ns = peak.gpsNanoSeconds
115 if 'stop_time' in columns: t.stop_time = stop.gpsSeconds
116 if 'stop_time_ns' in columns: t.stop_time_ns = stop.gpsNanoSeconds
117
118
119 if 'ms_start_time' in columns: t.ms_start_time = start.gpsSeconds
120 if 'ms_start_time_ns' in columns: t.ms_start_time_ns = start.gpsNanoSeconds
121 if 'ms_stop_time' in columns: t.ms_stop_time = stop.gpsSeconds
122 if 'ms_stop_time_ns' in columns: t.ms_stop_time_ns = stop.gpsNanoSeconds
123
124
125 if 'duration' in columns: t.duration = duration
126 if 'ms_duration' in columns: t.ms_duration = duration
127
128
129 if 'central_freq' in columns: t.central_freq = freq
130 if 'peak_frequency' in columns: t.peak_frequency = av_freq
131 if 'peak_frequency_eror' in columns: t.peak_frequency_error = err_freq
132 if 'bandwidth' in columns: t.bandwidth = av_band
133 if 'flow' in columns: t.flow = freq-band/2
134 if 'fhigh' in columns: t.fhigh = freq+band/2
135
136
137 if 'ms_bandwidth' in columns: t.ms_bandwidth = band
138 if 'ms_flow' in columns: t.ms_flow = freq-band/2
139 if 'ms_fhigh' in columns: t.ms_fhigh = freq+band/2
140
141
142 if 'snr' in columns: t.snr = snr
143 if 'ms_snr' in columns: t.ms_snr = snr
144 if 'amplitude' in columns: t.amplitude = amplitude
145
146
147 if 'cluster_size' in columns or 'param_one_value' in columns:
148 t.param_one_name = 'cluster_size'
149 if clusters:
150 t.param_one_value = cls
151 else:
152 t.param_one_value = numpy.NaN
153 if 'cluster_norm_energy' in columns or 'param_two_value' in columns:
154 t.param_two_name = 'cluster_norm_energy'
155 if clusters:
156 t.param_two_value = cle
157 else:
158 t.param_two_value = numpy.NaN
159 if 'cluster_number' in columns or 'param_three_value' in columns:
160 t.param_three_name = 'cluster_number'
161 if clusters:
162 t.param_three_value = cln
163 else:
164 t.param_three_value = numpy.NaN
165
166 return t
167
168
169
170
171
172 -def fromfile(fobj, start=None, end=None, ifo=None, channel=None,\
173 columns=None, virgo=False):
174
175 """
176 Load triggers from an Omega format text file into a SnglBurstTable object.
177 Use start and end to restrict the returned triggers, and give ifo and
178 channel to fill those columns in the table.
179
180 If columns is given as a list, only those columns in the table will be
181 filled. This is advisable to speed up future operations on this table.
182
183 Arguments :
184
185 fname : file or str
186 file object or filename path to read with numpy.loadtext
187
188 Keyword arguments :
189
190 start : float
191 minimum peak time for returned triggers
192 end : float
193 maximum peak time for returned triggers
194 ifo : str
195 name of IFO to fill in table
196 channel : str
197 name of channel to fill in table
198 columns : iterable
199 list of columnnames to populate in table
200 virgo : [ True | False ]
201 fobj written in Virgo OmegaOnline format
202 """
203
204
205 if columns==None:
206 columns = lsctables.SnglBurst.__slots__
207 usercolumns = False
208 else:
209 usercolumns = True
210 if start or end:
211 if not start:
212 start = -numpy.inf
213 if not end:
214 end = numpy.inf
215 span = segments.segment(start, end)
216 if 'peak_time' not in columns: columns.append('peak_time')
217 if 'peak_time_ns' not in columns: columns.append('peak_time_ns')
218 check_time = True
219 else:
220 check_time = False
221
222
223 if 'snr' in columns and not 'amplitude' in columns:
224 columns.append('amplitude')
225
226
227
228
229
230 out = lsctables.New(lsctables.SnglBurstTable, columns=columns)
231 append = out.append
232
233
234 if usercolumns:
235
236 if isinstance(columns[0], str):
237 columns = map(str.lower, columns)
238 if isinstance(columns[0], unicode):
239 columns = map(unicode.lower, columns)
240
241 for c in out.columnnames:
242 if c.lower() not in columns:
243 idx = out.columnnames.index(c)
244 out.columnnames.pop(idx)
245 out.columntypes.pop(idx)
246
247
248
249
250
251
252 if hasattr(fobj, 'readline'):
253 fh = fobj
254 else:
255 fh = open(fobj, 'r')
256
257
258 for i,line in enumerate(fh):
259 if _comment.match(line): continue
260 t = trigger(line, columns=columns, virgo=virgo, channel=channel, ifo=ifo)
261 if not check_time or (check_time and float(t.get_peak()) in span):
262 append(t)
263
264
265 if not hasattr(fobj, 'readline'):
266 fh.close()
267
268 return out
269
270
271
272
273
274 -def tofrequencyseries(bursttable, fcol='peak_frequency', pcol=None,\
275 name="", epoch=LIGOTimeGPS(), deltaF=0, f0=0,\
276 unit=lalStrainUnit):
277 """
278 Returns a numpy.array and REAL8FrequencySeries built from these
279 OmegaSpectrum triggers. The array holds the discrete frequencies at
280 which the sectrum was calculated and the series holds the data and
281 associated metadata.
282
283 If pcol is not given, the series data is the square of the SNR of each
284 'trigger'.
285 """
286
287 freq = bursttable.getColumnByName('peak_frequency')
288 if pcol:
289 data = bursttable.getColumnByName('pcol')
290 else:
291 data = bursttable.getColumnByName('snr')**2
292 freq,data = list(map(numpy.asarray, zip(*sorted(zip(freq, data)))))
293
294 if int(epoch) == 0 and len(bursttable) != 0:
295 epoch = LIGOTimeGPS(float(bursttable[0].get_time()))
296 if deltaF == 0 and len(bursttable) > 1:
297 deltaF = freq[1]-freq[0]
298 if f0 == 0 and len(bursttable) != 0:
299 f0 = freq[0]
300
301 series = seriesutils.fromarray(data, name=name, epoch=epoch, deltaT=deltaF,\
302 f0=f0, unit=unit, frequencyseries=True)
303 return freq,series
304
305
306
307
308
309 -def get_cache(start, end, ifo, channel, mask='DOWNSELECT', checkfilesexist=False,\
310 **kwargs):
311 """
312 Returns a glue.lal.Cache contatining CacheEntires for all omega online
313 trigger files between the given start and end time for the given ifo.
314 """
315 cache = Cache()
316
317
318 host = { 'G1':'atlas', 'H1':'ligo-wa', 'H2':'ligo-wa', 'L1':'ligo-la'}
319 if (not kwargs.has_key('directory') and not re.search(host[ifo],getfqdn())):
320 sys.stderr.write("warning: Omega online files are not available for "+\
321 "IFO=%s on this host." % ifo)
322 sys.stderr.flush()
323 return cache
324
325 span = segments.segment(start,end)
326 if ifo == 'G1':
327 if channel:
328 kwargs.setdefault('directory', '/home/omega/online/%s/segments' % channel.replace(':','_'))
329 else:
330 kwargs.setdefault('directory', '/home/omega/online/G1/segments')
331 kwargs.setdefault('epoch', 0)
332 else:
333 kwargs.setdefault('directory',\
334 '/home/omega/online/%s/archive/S6/segments' % ifo)
335 kwargs.setdefault('epoch', 931211808)
336 kwargs.setdefault('duration', 64)
337 kwargs.setdefault('overlap', 8)
338
339
340 append = cache.append
341 splitext = os.path.splitext
342 isfile = os.path.isfile
343 intersects = span.intersects
344 segment = segments.segment
345 from_T050017 = CacheEntry.from_T050017
346 basedir = kwargs['directory']
347 basetime = kwargs['epoch']
348 triglength = kwargs['duration']
349 overlap = kwargs['overlap']
350
351
352 start_time = int(start-numpy.mod(start-basetime,triglength-overlap))
353 t = start_time
354
355
356 while t<end:
357 if ifo == 'G1':
358 trigfile = '%s/%.5d/%.10d-%10.d/%s-OMEGA_TRIGGERS_%s-%.10d-%d.txt'\
359 % (basedir, t/100000, t, t+triglength, ifo, mask, t, triglength)
360 else:
361 trigfile = '%s/%.10d-%10.d/%s-OMEGA_TRIGGERS_%s-%.10d-%d.txt'\
362 % (basedir, t, t+triglength, ifo, mask, t, triglength)
363 if intersects(segment(t, t+triglength))\
364 and (not checkfilesexist or isfile(trigfile)):
365 append(from_T050017(trigfile))
366 t+=triglength-overlap
367
368 cache.sort(key=lambda e: e.path)
369
370 return cache
371
372
373
374
375
376 -def fromfiles(filelist, start=None, end=None, columns=None, verbose=False,\
377 virgo=False, channel=None):
378 """
379 Read omega triggers from a list of ASCII filepaths.
380 """
381
382 if verbose:
383 sys.stdout.write("Extracting Omega triggers from %d files... "\
384 % len(filelist))
385 sys.stdout.flush()
386 delete = '\b\b\b'
387 num = len(filelist)/100
388
389 out = lsctables.New(lsctables.SnglBurstTable, columns=columns)
390 extend = out.extend
391
392 for i,fp in enumerate(filelist):
393 with open(fp, "r") as f:
394 extend(fromfile(f, start=start, end=end, columns=columns,\
395 virgo=virgo, channel=channel))
396 if verbose:
397 progress = int((i+1)/num)
398 sys.stdout.write('%s%.2d%%' % (delete, progress))
399 sys.stdout.flush()
400
401 if verbose:
402 sys.stdout.write("\n")
403
404 return out
405
406 -def fromlalcache(cache, start=None, end=None, columns=None, verbose=False,\
407 virgo=False,channel=None):
408 """
409 Read omega triggers from a glue.lal.Cache list of ASCII CacheEntries.
410 """
411 return fromfiles(cache.pfnlist(), start=start, end=end, columns=columns,\
412 verbose=verbose, virgo=virgo, channel=channel)
413