pylal.MultiInspiralUtils

1 # Copyright (C) 2006 Sukanta Bose 2 # 3 # This program is free software; you can redistribute it and/or modify it 4 # under the terms of the GNU General Public License as published by the 5 # Free Software Foundation; either version 2 of the License, or (at your 6 # option) any later version. 7 # 8 # This program is distributed in the hope that it will be useful, but 9 # WITHOUT ANY WARRANTY; without even the implied warranty of 10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 # Public License for more details. 12 # 13 # You should have received a copy of the GNU General Public License along 14 # with this program; if not, write to the Free Software Foundation, Inc., 15 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 16 17 # 18 # ============================================================================= 19 # 20 # Preamble 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 # Input 40 # 41 # ============================================================================= 42 # 43

44 -def ReadMultiInspiralFromFiles(fileList):

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 # extract the multi inspiral table 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

65 -def ReadMultiInspiralTimeSlidesFromFiles(fileList,generate_output_tables=False):

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 # Extract the time slide table 82 timeSlideTable = lsctables.TimeSlideTable.get_table(doc) 83 slideMapping = {} 84 currSlides = {} 85 # NOTE: I think some of this is duplicated in the glue definition of the 86 # time slide table. Probably should move over to that 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 # Is the slide already in the list and where? 98 offsetIndex = timeSlides.index(offsetDict) 99 slideMapping[slideID] = offsetIndex 100 except ValueError: 101 # If not then add it 102 timeSlides.append(offsetDict) 103 slideMapping[slideID] = len(timeSlides) - 1 104 105 # Get the mapping table 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 # Extract the segment table 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 # extract the multi inspiral table 122 try: 123 multiInspiralTable = lsctables.MultiInspiralTable.get_table(doc) 124 # Remap the time slide IDs 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 # except: multiInspiralTable = None 132 except: raise 133 134 if not generate_output_tables: 135 return multis,timeSlides,segmentDict 136 else: 137 # Make a new time slide table 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 # Make a new mapping table 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 # Make a new segment table 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 # Clustering 189 # 190 # ============================================================================= 191 # 192

193 -def CompareMultiInspiralByEndTime(a, b):

194 """ 195 Orders a and b by peak time. 196 """ 197 return cmp(a.get_end(), b.get_end())

198 199

200 -def CompareMultiInspiralBySnr(a, b):

201 """ 202 Orders a and b by peak time. 203 """ 204 return cmp(a.snr, b.snr)

205 206

207 -def CompareMultiInspiral(a, b, twindow = LIGOTimeGPS(0)):

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

218 -def cluster_multi_inspirals(mi_table, dt, loudest_by="snr"):

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 # get data 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 # get times 254 start = round(end_time.min()) 255 end = round(end_time.max()+1) 256 257 # generate bins 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 # bin triggers 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 # loop over all bins 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 # pick loudest event in bin 286 idx = bin_[stat[bin_].argmax()] 287 s = stat[idx] 288 t = end_time[idx] 289 290 # trigger was loudest in it's bin, search loudest event 291 # in previous bin 292 if (check_prev and (t - loudest_time[prev]) < dt and 293 s < loudest_stat[prev]): 294 continue 295 296 # Same for the next bin 297 if (check_next and (loudest_time[next_] - t) < dt and 298 s < loudest_stat[next_]): 299 continue 300 301 loudest=True 302 303 # trigger was loudest in it's bin, search previous bin 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 # if still loudest, check the next bin 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 # this was the loudest trigger in its vicinity, 323 # keep it and move to the next bin 324 cluster_table.append(mi_table[idx]) 325 326 return cluster_table

327

Source Code for Module pylal.MultiInspiralUtils