1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25 import numpy
26 import itertools
27
28 from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS
29 from glue import segments
30 from glue.ligolw import table
31 from glue.ligolw import lsctables
32 from glue.ligolw import utils
33 from glue.ligolw import ilwd
34 from glue.ligolw import ligolw
35
36
37
38
39
40
41
42
43
45 """
46 Read the multiInspiral tables from a list of files
47 @param fileList: list of input files
48 """
49 if not fileList:
50 return multiInspiralTable(), None
51
52 multis = None
53
54 for thisFile in fileList:
55 doc = utils.load_filename(thisFile,
56 gz=(thisFile or "stdin").endswith(".gz"), contenthandler = lsctables.use_in(ligolw.LIGOLWContentHandler))
57
58 try:
59 multiInspiralTable = lsctables.MultiInspiralTable.get_table(doc)
60 if multis: multis.extend(multiInspiralTable)
61 else: multis = multiInspiralTable
62 except: multiInspiralTable = None
63 return multis
64
66 """
67 Read time-slid multiInspiral tables from a list of files
68 @param fileList: list of input files
69 """
70 if not fileList:
71 return multiInspiralTable(), None
72
73 multis = None
74 timeSlides = []
75
76 segmentDict = {}
77 for thisFile in fileList:
78
79 doc = utils.load_filename(thisFile,
80 gz=(thisFile or "stdin").endswith(".gz"), contenthandler = lsctables.use_in(ligolw.LIGOLWContentHandler))
81
82 timeSlideTable = lsctables.TimeSlideTable.get_table(doc)
83 slideMapping = {}
84 currSlides = {}
85
86
87 for slide in timeSlideTable:
88 currID = int(slide.time_slide_id)
89 if currID not in currSlides.keys():
90 currSlides[currID] = {}
91 currSlides[currID][slide.instrument] = slide.offset
92 elif slide.instrument not in currSlides[currID].keys():
93 currSlides[currID][slide.instrument] = slide.offset
94
95 for slideID,offsetDict in currSlides.items():
96 try:
97
98 offsetIndex = timeSlides.index(offsetDict)
99 slideMapping[slideID] = offsetIndex
100 except ValueError:
101
102 timeSlides.append(offsetDict)
103 slideMapping[slideID] = len(timeSlides) - 1
104
105
106 segmentMap = {}
107 timeSlideMapTable = lsctables.TimeSlideSegmentMapTable.get_table(doc)
108 for entry in timeSlideMapTable:
109 segmentMap[int(entry.segment_def_id)] = int(entry.time_slide_id)
110
111
112 segmentTable = lsctables.SegmentTable.get_table(doc)
113 for entry in segmentTable:
114 currSlidId = segmentMap[int(entry.segment_def_id)]
115 currSeg = entry.get()
116 if not segmentDict.has_key(slideMapping[currSlidId]):
117 segmentDict[slideMapping[currSlidId]] = segments.segmentlist()
118 segmentDict[slideMapping[currSlidId]].append(currSeg)
119 segmentDict[slideMapping[currSlidId]].coalesce()
120
121
122 try:
123 multiInspiralTable = lsctables.MultiInspiralTable.get_table(doc)
124
125 for multi in multiInspiralTable:
126 newID = slideMapping[int(multi.time_slide_id)]
127 multi.time_slide_id = ilwd.ilwdchar(\
128 "time_slide:time_slide_id:%d" % (newID))
129 if multis: multis.extend(multiInspiralTable)
130 else: multis = multiInspiralTable
131
132 except: raise
133
134 if not generate_output_tables:
135 return multis,timeSlides,segmentDict
136 else:
137
138 timeSlideTab = lsctables.New(lsctables.TimeSlideTable)
139
140 for slideID,offsetDict in enumerate(timeSlides):
141 for instrument in offsetDict.keys():
142 currTimeSlide = lsctables.TimeSlide()
143 currTimeSlide.instrument = instrument
144 currTimeSlide.offset = offsetDict[instrument]
145 currTimeSlide.time_slide_id = ilwd.ilwdchar(\
146 "time_slide:time_slide_id:%d" % (slideID))
147 currTimeSlide.process_id = ilwd.ilwdchar(\
148 "process:process_id:%d" % (0))
149 timeSlideTab.append(currTimeSlide)
150
151
152 timeSlideSegMapTab = lsctables.New(lsctables.TimeSlideSegmentMapTable)
153
154 for i in range(len(timeSlides)):
155 currMapEntry = lsctables.TimeSlideSegmentMap()
156 currMapEntry.time_slide_id = ilwd.ilwdchar(\
157 "time_slide:time_slide_id:%d" % (i))
158 currMapEntry.segment_def_id = ilwd.ilwdchar(\
159 "segment_def:segment_def_id:%d" % (i))
160 timeSlideSegMapTab.append(currMapEntry)
161
162
163 newSegmentTable = lsctables.New(lsctables.SegmentTable)
164
165 segmentIDCount = 0
166 for i in range(len(timeSlides)):
167 currSegList = segmentDict[i]
168 for seg in currSegList:
169 currSegment = lsctables.Segment()
170 currSegment.segment_id = ilwd.ilwdchar(\
171 "segment:segment_id:%d" %(segmentIDCount))
172 segmentIDCount += 1
173 currSegment.segment_def_id = ilwd.ilwdchar(\
174 "segment_def:segment_def_id:%d" % (i))
175 currSegment.process_id = ilwd.ilwdchar(\
176 "process:process_id:%d" % (0))
177 currSegment.set(seg)
178 currSegment.creator_db = -1
179 currSegment.segment_def_cdb = -1
180 newSegmentTable.append(currSegment)
181 return multis,timeSlides,segmentDict,timeSlideTab,newSegmentTable,\
182 timeSlideSegMapTab
183
184
185
186
187
188
189
190
191
192
194 """
195 Orders a and b by peak time.
196 """
197 return cmp(a.get_end(), b.get_end())
198
199
201 """
202 Orders a and b by peak time.
203 """
204 return cmp(a.snr, b.snr)
205
206
208 """
209 Returns 0 if a and b are less than twindow appart.
210 """
211 tdiff = abs(a.get_end() - b.get_end())
212 if tdiff < twindow:
213 return 0
214 else:
215 return cmp(a.get_end(), b.get_end())
216
217
219 """Cluster a MultiInspiralTable with a given ranking statistic and
220 clustering window.
221
222 This method separates the table rows into time bins, returning those
223 row that are
224 * loudest in their own bin, and
225 * louder than those events in the preceeding and following bins
226 that are within the clustering time window
227
228 @return: a new MultiInspiralTable containing those clustered events
229
230 @param mi_table:
231 MultiInspiralTable to cluster
232 @param dt:
233 width (seconds) of clustering window
234 @keyword loudest_by:
235 column by which to rank events, default: 'snr'
236
237 @type mi_table: glue.ligolw.lsctables.MultiInspiralTable
238 @type dt: float
239 @type loudest_by: string
240 @rtype: glue.ligolw.lsctables.MultiInspiralTable
241 """
242 cluster_table = table.new_from_template(mi_table)
243 if not len(mi_table):
244 return cluster_table
245
246
247 end_time = numpy.asarray(mi_table.get_end()).astype(float)
248 if hasattr(mi_table, "get_%s" % loudest_by):
249 stat = numpy.asarray(getattr(mi_table, "get_%s" % loudest_by)())
250 else:
251 stat = numpy.asarray(mi_table.get_column(loudest_by))
252
253
254 start = round(end_time.min())
255 end = round(end_time.max()+1)
256
257
258 num_bins = int((end-start)//dt + 1)
259 time_bins = []
260 loudest_stat = numpy.zeros(num_bins)
261 loudest_time = numpy.zeros(num_bins)
262 for n in range(num_bins):
263 time_bins.append([])
264
265
266 for i,(t,s) in enumerate(itertools.izip(end_time, stat)):
267 bin_ = int(float(t-start)//dt)
268 time_bins[bin_].append(i)
269 if s > loudest_stat[bin_]:
270 loudest_stat[bin_] = s
271 loudest_time[bin_] = t
272
273
274 for i,bin_ in enumerate(time_bins):
275 if len(bin_)<1:
276 continue
277 first = i==0
278 last = i==(num_bins-1)
279
280 prev = i-1
281 next_ = i+1
282 check_prev = (not first and len(time_bins[prev]) > 0)
283 check_next = (not last and len(time_bins[next_]) > 0)
284
285
286 idx = bin_[stat[bin_].argmax()]
287 s = stat[idx]
288 t = end_time[idx]
289
290
291
292 if (check_prev and (t - loudest_time[prev]) < dt and
293 s < loudest_stat[prev]):
294 continue
295
296
297 if (check_next and (loudest_time[next_] - t) < dt and
298 s < loudest_stat[next_]):
299 continue
300
301 loudest=True
302
303
304 if check_prev and not (t - loudest_time[prev]) < dt:
305 for idx2 in time_bins[prev]:
306 t2 = end_time[idx2]
307 if (t - end_time[idx2]) < dt and s < stat[idx2]:
308 loudest = False
309 break
310 if not loudest:
311 continue
312
313
314 if check_next and not (loudest_time[next_] - t) < dt:
315 for idx2 in time_bins[next_]:
316 if (end_time[idx2] - t) < dt and s < stat[idx2]:
317 loudest = False
318 break
319 if not loudest:
320 continue
321
322
323
324 cluster_table.append(mi_table[idx])
325
326 return cluster_table
327