pylal.SnglInspiralUtils

1 # Copyright (C) 2006 Duncan A. Brown 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 copy 26 27 from pylal import SearchSummaryUtils 28 from glue.ligolw import ligolw 29 from glue.ligolw import table 30 from glue.ligolw import lsctables 31 from glue.ligolw import utils 32 from glue.ligolw.utils import ligolw_add 33 from glue import iterutils 34 from glue import segments 35 36 37 # 38 # ============================================================================= 39 # 40 # Input 41 # 42 # ============================================================================= 43 # 44

45 -class ExtractSnglInspiralTableLIGOLWContentHandler(ligolw.PartialLIGOLWContentHandler):

46 """ 47 LIGOLWContentHandler that will extract only the SnglInspiralTable from a document. 48 See glue.ligolw.LIGOLWContentHandler help for more info. 49 """

50 - def __init__(self,document):

51 def filterfunc(name,attrs): 52 return name==ligolw.Table.tagName and attrs.has_key('Name') and table.Table.TableName(attrs.get('Name'))==lsctables.SnglInspiralTable.tableName

53 ligolw.PartialLIGOLWContentHandler.__init__(self,document,filterfunc)

54 55

56 -class SnglInspiralID_old(object):

57 """ 58 Custom row ID thing for sngl_inspiral tables with int_8s event IDs. 59 """ 60 column_name = "event_id" 61

62 - def __init__(self, n = 0):

63 self.n = n

64

65 - def new(self, row):

66 self.n += 1 67 a = self.n // 100000 68 b = self.n % 100000 69 return lsctables.SnglInspiralID(a * 1000000000 + row.get_id_parts()[1] * 100000 + b)

70 71

72 -def ReadSnglInspiralFromFiles(fileList, verbose=False, filterFunc=None):

73 """ 74 Read the SnglInspiralTables from a list of files. 75 If filterFunc is not None, only keep triggers for which filterFunc 76 evaluates to True. Ex.: filterFunc=lambda sng: sng.snr >= 6.0 77 78 @param fileList: list of input files 79 @param verbose: print progress 80 """ 81 # NOTE: this function no longer carries out event ID mangling (AKA 82 # reassignment). Please adjust calling codes accordingly! 83 # This means that identical event IDs produced by lalapps_thinca in 84 # non-slide files having the same GPS start time will stay identical, 85 # affecting zerolag and injection runs made over the same data. 86 # 87 # In consequence, if the calling code is going to reconstruct coincs 88 # from the sngl event IDs, and if these include multiple injection 89 # runs, coinc finding should be done one file at a time - see the 90 # readCoincInspiralFromFiles function in CoincInspiralUtils.py 91 92 sngls = lsctables.New(lsctables.SnglInspiralTable, \ 93 columns=lsctables.SnglInspiralTable.loadcolumns) 94 95 lsctables.use_in(ExtractSnglInspiralTableLIGOLWContentHandler) 96 for i,file in enumerate(fileList): 97 if verbose: print str(i+1)+"/"+str(len(fileList))+": " 98 xmldoc = utils.load_filename(file, verbose=verbose, contenthandler=ExtractSnglInspiralTableLIGOLWContentHandler) 99 try: 100 sngl_table = lsctables.SnglInspiralTable.get_table(xmldoc) 101 if filterFunc is not None: 102 iterutils.inplace_filter(filterFunc, sngl_table) 103 except ValueError: #some xml files have no sngl table, that's OK 104 sngl_table = None 105 if sngl_table: sngls.extend(sngl_table) 106 xmldoc.unlink() #free memory 107 108 return sngls

109 110

111 -def ReadSnglInspiralSlidesFromFiles(fileList, shiftVector, vetoFile=None, 112 verbose=False):

113 """ 114 Function for reading time-slided single inspiral triggers 115 with automatic resliding the times, given a list of input files. 116 117 @param fileList: List of files containing single inspiral time-slided 118 triggers 119 @param shiftVector: Dictionary of time shifts to apply to triggers 120 keyed by IFO 121 @param vetoFile: segwizard formatted file used to veto all triggers 122 @param verbose: print progress 123 """ 124 # NOTE: This function also will not mangle (reassign) lalapps_thinca 125 # style event IDs (see the previous function). 126 # Hence it will fail if fed event ID-related options. 127 128 # read raw triggers 129 inspTriggers = ReadSnglInspiralFromFiles(fileList, verbose=verbose) 130 if inspTriggers: 131 # get the rings 132 segDict = SearchSummaryUtils.GetSegListFromSearchSummaries(fileList) 133 rings = segments.segmentlist(iterutils.flatten(segDict.values())) 134 rings.sort() 135 136 # perform the veto 137 if vetoFile is not None: 138 segList = segmentsUtils.fromsegwizard(open(vetoFile)) 139 inspTriggers = inspTriggers.veto(segList) 140 141 # now slide all the triggers within their appropriate ring 142 slideTriggersOnRingWithVector(inspTriggers, shiftVector, rings) 143 144 # return the re-slided triggers 145 return inspTriggers

146 147

148 -def ReadSnglInspiralsForPipelineStage(xmldoc, slideDict, stage):

149 """ 150 Collect the sngl_inspiral rows from a desired stage in the pipeline. 151 -- For the INSPIRAL and zerolag THINCA stages, the entire sngl_inspiral table is returned. 152 -- For the slide THINCA stages, return only the rows in the sngl_inspiral that 153 compose a coincident event from the desired time-slide 154 @param xmldoc: ligolw_xml doc 155 @param slideDict: dictionary of the desired time-slide (eg. {u'H1': 0.0, u'L1': 100.0}) 156 @param stage: the name of the stage (INSPIRAL_FIRST, THINCA_0_CAT_2, etc) 157 """ 158 159 sngls_tbl = lsctables.SnglInspiralTable.get_table(xmldoc) 160 if 'THINCA' in stage and slideDict: 161 # get the time-slides as a dictionary 162 time_slide_tbl = lsctables.TimeSlideTable.get_table(xmldoc) 163 time_slide_dict = time_slide_tbl.as_dict() 164 165 coinc_event_tbl = lsctables.CoincTable.get_table(xmldoc) 166 coinc_map_tbl = lsctables.CoincMapTable.get_table(xmldoc) 167 168 # get the time_slide_ids that have 169 time_slide_ids = set() 170 for tsid, offsets in time_slide_dict.items(): 171 if offsets.values() == slideDict.values(): 172 time_slide_ids.add( tsid ) 173 # get the coinc_event_ids for coincs from the desired slides 174 coinc_event_ids = set() 175 for coinc in coinc_event_tbl: 176 if coinc.time_slide_id in time_slide_ids: 177 coinc_event_ids.add( coinc.coinc_event_id ) 178 # get the inspiral event_ids associated with the above coinc events 179 event_ids = set() 180 for row in coinc_map_tbl: 181 if row.coinc_event_id in coinc_event_ids: 182 event_ids.add( row.event_id ) 183 184 sngls_tbl_eid = sngls_tbl.getColumnByName("event_id") 185 coinc_sngls_tbl = xmldoc.childNodes[0].insertBefore( lsctables.New(lsctables.SnglInspiralTable), sngls_tbl) 186 for idx, event_id in enumerate(event_ids): 187 coinc_sngls_tbl.insert( idx, sngls_tbl[sngls_tbl_eid.index(event_id)] ) 188 189 return coinc_sngls_tbl 190 else: 191 return sngls_tbl

192 193 194 # 195 # ============================================================================= 196 # 197 # Clustering 198 # 199 # ============================================================================= 200 # 201

202 -def CompareSnglInspiralByEndTime(a, b):

203 """ 204 Orders a and b by peak time. 205 """ 206 return cmp(a.get_end(), b.get_end())

207 208

209 -def CompareSnglInspiralBySnr(a, b):

210 """ 211 Orders a and b by peak time. 212 """ 213 return cmp(a.snr, b.snr)

214 215

216 -def CompareSnglInspiral(a, b, twindow = lsctables.LIGOTimeGPS(0)):

217 """ 218 Returns 0 if a and b are less than twindow appart. 219 """ 220 tdiff = abs(a.get_end() - b.get_end()) 221 if tdiff < twindow: 222 return 0 223 else: 224 return cmp(a.get_end(), b.get_end())

225 226 # 227 # ============================================================================= 228 # 229 # Time Sliding on Ring 230 # 231 # ============================================================================= 232 # 233

234 -def slideTriggersOnLines(triggerList, shifts):

235 """ 236 In-place modify trigger_list so that triggers are slid by appropriate value 237 of shifts. 238 239 @param triggerList: a SnglInspiralTable 240 @param shifts: a dictionary of the time-shifts keyed by IFO 241 """ 242 for trigger in triggerList: 243 end_time = trigger.get_end() 244 trigger.set_end( end_time + shifts[trigger.ifo] )

245

246 -def slideTimeOnRing(time, shift, ring):

247 """ 248 Return time after adding shift, but constrained to lie along the ring 249 """ 250 assert time in ring 251 # use ring[0] as an epoch, do arithmetic using floats relative to epoch 252 return ring[0] + (float(time - ring[0]) + shift) % float(abs(ring))

253 254

255 -def slideTriggersOnRings(triggerList, rings, shifts):

256 """ 257 In-place modify trigger_list so that triggers are slid by appropriate value 258 of shifts along their enclosing ring segment by the algorithm given in XXX. 259 This function calls the function slideTimeOnRing 260 261 @param triggerList: a SnglInspiralTable 262 @param rings: sorted segment list of possible rings 263 @param shifts: a dictionary of the time-shifts keyed by IFO 264 """ 265 for trigger in triggerList: 266 end_time = trigger.get_end() 267 trigger.set_end(slideTimeOnRing(end_time, shifts[trigger.ifo], rings[rings.find(end_time)]))

268

269 -def unslideTriggersOnRings(triggerList, rings, shifts):

270 """ 271 In-place modify trigger_list so that triggers are unslid by appropriate 272 value of shifts along their enclosing ring segment by the algorithm given in 273 XXX. 274 This function calls the function slideTriggersOnRing 275 276 @param triggerList: a SnglInspiralTable 277 @param rings: sorted segment list of possible rings 278 @param shifts: a dictionary of the time-shifts keyed by IFO 279 """ 280 slideTriggersOnRings(triggerList, rings, dict((ifo, -shift) for ifo, shift in shifts.items()))

281

282 -def slideTriggersOnRingWithVector(triggerList, shiftVector, rings):

283 """ 284 In-place modify trigger_list so that triggers are slid by 285 along their enclosing ring segment by the algorithm given in XXX. 286 Slide numbers are extracted from the event_id of each trigger, 287 and multiplied by the corresponding (ifo-keyed) entry in shift_vector 288 to get the total slide amount. 289 This function is called by ReadSnglInspiralSlidesFromFiles and 290 calls the function slideTimeOnRing 291 292 @param triggerList: a SnglInspiralTable 293 @param shiftVector: a dictionary of the unit time-shift vector, 294 keyed by IFO 295 @param rings: sorted segment list of possible rings 296 """ 297 for trigger in triggerList: 298 end_time = trigger.get_end() 299 trigger.set_end(slideTimeOnRing(end_time, trigger.get_slide_number() * shiftVector[trigger.ifo], rings[rings.find(end_time)]))

300

301 -def slideSegListDictOnRing(ring, seglistdict, shifts):

302 """ 303 Return seglistdict with segments that are slid by appropriate values of 304 shifts along the ring segment by the algorithm given in XXX. 305 306 @param ring: segment on which to cyclicly slide segments in 307 seglistdict 308 @param seglistdict: segments to be slid on ring 309 @param shifts: a dictionary of the time-shifts keyed by IFO 310 """ 311 # don't do this in loops 312 ring_duration = float(abs(ring)) 313 314 # automate multi-list arithmetic 315 ring = segments.segmentlistdict.fromkeys(seglistdict.keys(), segments.segmentlist([ring])) 316 317 # make a copy, and extract the segments that are in the ring. 318 seglistdict = seglistdict & ring 319 320 # apply the shift vector. the shift vector is first normalized so that 321 # all the shifts are non-negative and less than the duration of the ring. 322 seglistdict.offsets.update(dict((instrument, shift % ring_duration) for instrument, shift in shifts.items())) 323 324 # split the result into the pieces that are still in the ring, and the 325 # pieces that have fallen off the edge. both retain the shift vector in 326 # the offsets attribute. 327 extra = seglistdict - ring 328 seglistdict &= ring 329 330 # wrap the fallen-off pieces around the ring. 331 for instrument in extra.keys(): 332 extra.offsets[instrument] -= ring_duration 333 334 # return the union of the results. the shifts vector is retained in the 335 # offsets attribute of the result 336 return seglistdict | extra

337 338

339 -def compute_thinca_livetime(on_instruments, off_instruments, rings, vetoseglistdict, offsetvectors):

340 """ 341 @on_instruments is an iterable of the instruments that must be on. 342 343 @off_instruments is an iterable of the instruments that must be off. 344 345 on_instruments and off_instruments must be disjoint. 346 347 @rings is a list of segments defining the analysis ring boundaries. They 348 can overlap, and do not need to be ordered. 349 350 @vetoseglistdict is a coalesced glue.segments.segmentlistdict object 351 providing the veto segments for whatever instruments have vetoes defined 352 for them. This can include veto lists for instruments other than those 353 listed in on_ and off_instruments (extra veto lists will be ignored), and 354 it need not provide lists for all instruments (instruments for which 355 there are no veto segment lists are assumed to be on at all times). 356 357 @offsetvectors is an iterable of dictionaries of instrument-offset pairs. 358 Each dictionary must contain entries for all instruments in the union of 359 on_instruments and off_instruments (it is allowed to name others as well, 360 but they will be ignored). An example of one dictionary of 361 instrument-offset pairs: {"H1": 0.0, "H2": 5.0, "L1": 10.0}. 362 363 The return value is a float giving the livetime in seconds. 364 """ 365 # local copies so they can be modified and iterated over more than once 366 # (in case generator expressions have been passed in) 367 on_instruments = set(on_instruments) 368 off_instruments = set(off_instruments) 369 370 # check that the on and off instruments are disjoint 371 if on_instruments & off_instruments: 372 raise ValueError, "on_instruments and off_instruments not disjoint" 373 374 # instruments that are not vetoed are assumed to be on 375 on_instruments &= set(vetoseglistdict.keys()) 376 377 # performance aid: only need offsets for instruments whose state is 378 # important 379 all_instruments = on_instruments | off_instruments 380 offsetvectors = tuple(dict((key, value) for key, value in offsetvector.items() if key in all_instruments) for offsetvector in offsetvectors) 381 382 # performance aid: if there are no offset vectors to consider, the 383 # livetime is trivial 384 if not offsetvectors: 385 return [] 386 387 # check that each offset vector provides values for all instruments of 388 # interest 389 for offsetvector in offsetvectors: 390 if not set(offsetvector.keys()).issuperset(all_instruments): 391 raise ValueError, "incomplete offset vector %s; missing instrument(s) %s" % (repr(offsetvector), ", ".join(all_instruments - set(offsetvector.keys()))) 392 393 # initialize the livetime sums 394 live_time = [0.0] * len(offsetvectors) 395 396 # the livetime is trivial if an instrument that must be off is never 397 # vetoed 398 if not set(vetoseglistdict.keys()).issuperset(off_instruments): 399 return live_time 400 401 # performance aid: don't need veto segment lists for instruments whose 402 # state is unimportant, nor veto segments that don't intersect the rings 403 coalesced_rings = segments.segmentlist(rings).coalesce() 404 vetoseglistdict = segments.segmentlistdict((key, segments.segmentlist(seg for seg in seglist if coalesced_rings.intersects_segment(seg))) for key, seglist in vetoseglistdict.items() if key in all_instruments) 405 406 # tot up the time when exactly the instruments that must be on are on 407 for ring in rings: 408 # don't do this in loops 409 ring = segments.segmentlist([ring]) 410 411 # performance aid: this is done in the loop, inside 412 # slideSegListDictOnRing(), but we can make that go faster by doing it 413 # here first 414 clipped_vetoseglistdict = segments.segmentlistdict((key, seglist & ring) for key, seglist in vetoseglistdict.items()) 415 416 # performance aid: if an instrument that must be vetoed is never 417 # vetoed in this ring, the livetime is zero 418 if not all(clipped_vetoseglistdict[key] for key in off_instruments): 419 continue 420 421 # iterate over offset vectors 422 for n, offsetvector in enumerate(offsetvectors): 423 # apply the offset vector to the vetoes, wrapping around the ring 424 slidvetoes = slideSegListDictOnRing(ring[0], clipped_vetoseglistdict, offsetvector) 425 426 # slidvetoes = times when instruments are vetoed, 427 # slidvetoes.union(on_instruments) = times when an instrument that 428 # must be on is vetoed 429 # 430 # ~slidvetoes = times when instruments are not vetoed, 431 # (~slidvetoes).union(off_instruments) = times when an instrument 432 # that must be off is not vetoed 433 live_time[n] += float(abs(ring - slidvetoes.union(on_instruments) - (~slidvetoes).union(off_instruments))) 434 435 # done 436 return live_time

437

Source Code for Module pylal.SnglInspiralUtils