pylal.lalapps_cbc

1 # Copyright (C) 2006 Kipp Cannon 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 # 21 # Preamble 22 # 23 # ============================================================================= 24 # 25 26 27 """ 28 Inspiral injection identification library. Contains code providing the 29 capacity to search a list of sngl_inspiral candidates for events 30 matching entries in a sim_inspiral list of software injections, recording the 31 matches as inspiral <--> injection coincidences using the standard coincidence 32 infrastructure. Also, any pre-recorded inspiral <--> inspiral coincidences are 33 checked for cases where all the inspiral events in a coincidence match an 34 injection, and these are recorded as coinc <--> injection coincidences, 35 again using the standard coincidence infrastructure. 36 """ 37 38 39 import bisect 40 import sys 41 42 43 from glue.ligolw import lsctables 44 from glue.ligolw.utils import coincs as ligolw_coincs 45 from glue.ligolw.utils import process as ligolw_process 46 from pylal import git_version 47 from lalinspiral import thinca 48 from pylal import ligolw_rinca 49 from lalburst import timeslides as ligolw_tisi 50 from pylal import SimInspiralUtils 51 from pylal import SnglInspiralUtils 52 53 54 __author__ = "Kipp Cannon <kipp.cannon@ligo.org>" 55 __version__ = "git id %s" % git_version.id 56 __date__ = git_version.date 57 58 59 # 60 # ============================================================================= 61 # 62 # Speed Hacks 63 # 64 # ============================================================================= 65 # 66 67

68 -def sngl_inspiral___cmp__(self, other):

69 # compare self's end time to the LIGOTimeGPS instance other 70 return cmp(self.end_time, other.seconds) or cmp(self.end_time_ns, other.nanoseconds)

71 72 73 lsctables.SnglInspiral.__cmp__ = sngl_inspiral___cmp__ 74

75 -def sngl_ringdown___cmp__(self, other):

76 # compare self's end time to the LIGOTimeGPS instance other 77 return cmp(self.start_time, other.seconds) or cmp(self.start_time_ns, other.nanoseconds)

78 79 80 lsctables.SnglRingdown.__cmp__ = sngl_ringdown___cmp__ 81 82 # 83 # ============================================================================= 84 # 85 # Document Interface 86 # 87 # ============================================================================= 88 # 89 90 91 InspiralSICoincDef = lsctables.CoincDef(search = u"inspiral", search_coinc_type = 1, description = u"sim_inspiral<-->sngl_inspiral coincidences") 92 InspiralSCNearCoincDef = lsctables.CoincDef(search = u"inspiral", search_coinc_type = 2, description = u"sim_inspiral<-->coinc_event coincidences (nearby)") 93 94 RingdownSICoincDef = lsctables.CoincDef(search = u"ringdown", search_coinc_type = 1, description = u"sim_ringdown<-->sngl_ringdown coincidences") 95 RingdownSCNearCoincDef= lsctables.CoincDef(search = u"ringdown", search_coinc_type = 2, description = u"sim_ringdown<-->coinc_event coincidences (nearby)") 96

97 -class DocContents(object):

98 """ 99 A wrapper interface to the XML document. 100 """

101 - def __init__(self, xmldoc, bbdef, sbdef, scndef, process, sngl_type, sim_type, get_sngl_time):

102 # 103 # store the process row 104 # 105 106 self.process = process 107 108 # 109 # locate the sngl_inspiral and sim_inspiral tables 110 # 111 112 self.sngltable = sngl_type.get_table(xmldoc) 113 try: 114 self.simtable = sim_type.get_table(xmldoc) 115 except ValueError: 116 self.simtable = lsctables.SimInspiralTable.get_table(xmldoc) 117 print >>sys.stderr,"No SimRingdownTable, use SimInspiralTable instead!" 118 119 # 120 # construct the zero-lag time slide needed to cover the 121 # instruments listed in all the triggers, then determine 122 # its ID (or create it if needed) 123 # 124 # FIXME: in the future, the sim_inspiral table should 125 # indicate time slide at which the injection was done 126 # 127 128 self.tisi_id = ligolw_tisi.get_time_slide_id(xmldoc, {}.fromkeys(self.sngltable.getColumnByName("ifo"), 0.0), create_new = process) 129 130 # 131 # get coinc_definer row for sim_type <--> sngl_type 132 # coincs; this creates a coinc_definer table if the 133 # document doesn't have one 134 # 135 136 self.sb_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, sbdef.search, sbdef.search_coinc_type, create_new = True, description = sbdef.description) 137 138 # 139 # get coinc_def_id's for sngl_type <--> sngl_type, and 140 # the sim_type <--> coinc_event coincs. set all 141 # to None if this document does not contain any sngl_type 142 # <--> sngl_type coincs. 143 # 144 145 try: 146 bb_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, bbdef.search, bbdef.search_coinc_type, create_new = False) 147 except KeyError: 148 bb_coinc_def_id = None 149 self.scn_coinc_def_id = None 150 else: 151 self.scn_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, scndef.search, scndef.search_coinc_type, create_new = True, description = scndef.description) 152 153 # 154 # get coinc table, create one if needed 155 # 156 157 try: 158 self.coinctable = lsctables.CoincTable.get_table(xmldoc) 159 except ValueError: 160 self.coinctable = lsctables.New(lsctables.CoincTable) 161 xmldoc.childNodes[0].appendChild(self.coinctable) 162 self.coinctable.sync_next_id() 163 164 # 165 # get coinc_map table, create one if needed 166 # 167 168 try: 169 self.coincmaptable = lsctables.CoincMapTable.get_table(xmldoc) 170 except ValueError: 171 self.coincmaptable = lsctables.New(lsctables.CoincMapTable) 172 xmldoc.childNodes[0].appendChild(self.coincmaptable) 173 174 # 175 # index the document 176 # 177 # FIXME: type<-->type coincs should be organized by time 178 # slide ID, but since injections are only done at zero lag 179 # for now this is ignored. 180 # 181 182 # index sngl_type table 183 index = {} 184 for row in self.sngltable: 185 index[row.event_id] = row 186 # find IDs of type<-->type coincs 187 self.coincs = {} 188 for coinc in self.coinctable: 189 if coinc.coinc_def_id == bb_coinc_def_id: 190 self.coincs[coinc.coinc_event_id] = [] 191 # construct event list for each type<-->type coinc 192 for row in self.coincmaptable: 193 if row.coinc_event_id in self.coincs: 194 self.coincs[row.coinc_event_id].append(index[row.event_id]) 195 del index 196 # sort each event list by end/start time and convert to tuples 197 # for speed 198 199 for coinc_event_id, events in self.coincs.iteritems(): 200 events.sort(key=get_sngl_time) 201 self.coincs[coinc_event_id]=tuple(events) 202 # convert dictionary to a list 203 204 self.coincs = self.coincs.items() 205 206 # 207 # FIXME Is this true for inspirals too? 208 # sort sngl_type table by end/start time, and sort the coincs 209 # list by the end/start time of the first (earliest) event in 210 # each coinc (recall that the event tuple for each coinc 211 # has been time-ordered) 212 # 213 214 self.sngltable.sort(key=get_sngl_time) 215 self.coincs.sort(key=lambda(id,a): get_sngl_time(a[0])) 216 217 # 218 # set the window for type_near_time(). this window 219 # is the amount of time such that if an injection's end 220 # time and a inspiral event's end time differ by more than 221 # this it is *impossible* for them to match one another. 222 # 223 224 # FIXME I'll just make the windows 1.0 s 225 226 self.search_time_window = 1.0 227 self.coinc_time_window = 1.0

228 229

230 - def type_near_time(self, t):

231 """ 232 Return a list of the inspiral events whose peak times are 233 within self.search_time_window of t. 234 """ 235 return self.sngltable[bisect.bisect_left(self.sngltable, t - self.search_time_window):bisect.bisect_right(self.sngltable, t + self.search_time_window)]

236

237 - def coincs_near_time(self, t, get_time):

238 """ 239 Return a list of the (coinc_event_id, event list) tuples in 240 which at least one type event's end time is within 241 self.coinc_time_window of t. 242 """ 243 # FIXME: this test does not consider the time slide 244 # offsets that should be applied to the coinc, but for now 245 # injections are done at zero lag so this isn't a problem 246 # yet 247 return [(coinc_event_id, search_type) for coinc_event_id, search_type in self.coincs if (t - self.coinc_time_window <= get_time(search_type[-1])) and (get_time(search_type[0]) <= t + self.coinc_time_window)]

248

249 - def sort_triggers_by_id(self):

250 """ 251 Sort the sngl_table's rows by ID (tidy-up document 252 for output). 253 """ 254 self.sngltable.sort(lambda a, b: cmp(a.event_id, b.event_id))

255

256 - def new_coinc(self, coinc_def_id):

257 """ 258 Construct a new coinc_event row attached to the given 259 process, and belonging to the set of coincidences defined 260 by the given coinc_def_id. 261 """ 262 coinc = lsctables.Coinc() 263 coinc.process_id = self.process.process_id 264 coinc.coinc_def_id = coinc_def_id 265 coinc.coinc_event_id = self.coinctable.get_next_id() 266 coinc.time_slide_id = self.tisi_id 267 coinc.set_instruments(None) 268 coinc.nevents = 0 269 coinc.likelihood = None 270 self.coinctable.append(coinc) 271 return coinc

272 273 274 # 275 # ============================================================================= 276 # 277 # Add Process Information 278 # 279 # ============================================================================= 280 # 281 282 283 process_program_name = "lalapps_cbc_injfind" 284 285

286 -def append_process(xmldoc, match_algorithm, comment):

287 """ 288 Convenience wrapper for adding process metadata to the document. 289 """ 290 process = ligolw_process.append_process(xmldoc, program = process_program_name, version = __version__, cvs_repository = u"lscsoft", cvs_entry_time = __date__, comment = comment) 291 292 params = [(u"--match-algorithm", u"lstring", match_algorithm)] 293 ligolw_process.append_process_params(xmldoc, process, params) 294 295 return process

296 297 298 # 299 # ============================================================================= 300 # 301 # Injection <--> Inspiral/Ringdown Event Comparison Tests 302 # 303 # ============================================================================= 304 # 305

306 -def InspiralSnglCompare(sim, inspiral):

307 """ 308 Return False if the peak time of the sim is within 9 seconds of the inspiral event. 309 """ 310 return SnglInspiralUtils.CompareSnglInspiral(sim, inspiral, twindow = lsctables.LIGOTimeGPS(9))

311 312

313 -def InspiralNearCoincCompare(sim, inspiral):

314 """ 315 Return False if the peak time of the sim is within 9 seconds of the inspiral event. 316 """ 317 return SnglInspiralUtils.CompareSnglInspiral(sim, inspiral, twindow = lsctables.LIGOTimeGPS(9))

318

319 -def cmp_sngl_sim(sim, sngl, get_sim_time, get_sngl_time, twindow = lsctables.LIGOTimeGPS(9)):

320 tdiff = abs(get_sngl_time(sngl) - get_sim_time(sim)) 321 if tdiff < twindow: 322 return 0 323 else: 324 return cmp(get_sngl_time(sngl), get_sim_time(sim))

325 326 327 # 328 # ============================================================================= 329 # 330 # Build sim_type <--> sngl_type Coincidences 331 # 332 # ============================================================================= 333 # 334 335

336 -def find_sngl_type_matches(contents, sim, comparefunc, get_sim_time, get_sngl_time) :

337 """ 338 Scan the type table for triggers matching sim. 339 """ 340 return [type_table for type_table in contents.type_near_time(get_sim_time(sim)) if not comparefunc(sim, type_table, get_sim_time, get_sngl_time)]

341 342

343 -def add_sim_type_coinc(contents, sim, types):

344 """ 345 Create a coinc_event in the coinc table, and add arcs in the 346 coinc_event_map table linking the sim_type row and the list of 347 sngl_type rows to the new coinc_event row. 348 """ 349 coinc = contents.new_coinc(contents.sb_coinc_def_id) 350 coinc.set_instruments(set(event.ifo for event in types)) 351 coinc.nevents = len(types) 352 353 coincmap = lsctables.CoincMap() 354 coincmap.coinc_event_id = coinc.coinc_event_id 355 coincmap.table_name = sim.simulation_id.table_name 356 coincmap.event_id = sim.simulation_id 357 contents.coincmaptable.append(coincmap) 358 359 for event in types: 360 coincmap = lsctables.CoincMap() 361 coincmap.coinc_event_id = coinc.coinc_event_id 362 coincmap.table_name = event.event_id.table_name 363 coincmap.event_id = event.event_id 364 contents.coincmaptable.append(coincmap)

365 366 367 # 368 # ============================================================================= 369 # 370 # Build sim_type <--> coinc Coincidences 371 # 372 # ============================================================================= 373 # 374 375

376 -def find_coinc_matches(coincs, sim, comparefunc):

377 """ 378 Return a list of the coinc_event_ids of the inspiral<-->inspiral coincs 379 in which all inspiral events match sim. 380 """ 381 # FIXME: this test does not consider the time slide offsets that 382 # should be applied to the coinc, but for now injections are done 383 # at zero lag so this isn't a problem yet 384 return [coinc_event_id for coinc_event_id, inspirals in coincs if True not in [bool(comparefunc(sim, inspiral)) for inspiral in inspirals]]

385

386 -def find_near_coinc_matches(coincs, sim, comparefunc, get_sim_time, get_sngl_time):

387 """ 388 Return a list of the coinc_event_ids of the type<-->type coincs 389 in which at least one type event matches sim. 390 """ 391 # FIXME: this test does not consider the time slide offsets that 392 # should be applied to the coinc, but for now injections are done 393 # at zero lag so this isn't a problem yet 394 return [coinc_event_id for coinc_event_id, coinc_types in coincs if False in [bool(comparefunc(sim, coinc_type, get_sim_time, get_sngl_time)) for coinc_type in coinc_types]]

395 396

397 -def add_sim_coinc_coinc(contents, sim, coinc_event_ids, coinc_def_id):

398 """ 399 Create a coinc_event in the coinc table, and add arcs in the 400 coinc_event_map table linking the sim_type row and the list of 401 coinc_event rows to the new coinc_event row. 402 """ 403 coinc = contents.new_coinc(coinc_def_id) 404 coinc.nevents = len(coinc_event_ids) 405 406 coincmap = lsctables.CoincMap() 407 coincmap.coinc_event_id = coinc.coinc_event_id 408 coincmap.table_name = sim.simulation_id.table_name 409 coincmap.event_id = sim.simulation_id 410 contents.coincmaptable.append(coincmap) 411 412 for coinc_event_id in coinc_event_ids: 413 coincmap = lsctables.CoincMap() 414 coincmap.coinc_event_id = coinc.coinc_event_id 415 coincmap.table_name = coinc_event_id.table_name 416 coincmap.event_id = coinc_event_id 417 contents.coincmaptable.append(coincmap)

418 419 420 # 421 # ============================================================================= 422 # 423 # Library API 424 # 425 # ============================================================================= 426 # 427 428

429 -def lalapps_cbc_injfind(xmldoc, process, search, snglcomparefunc, nearcoinccomparefunc, verbose = False):

430 # 431 # Analyze the document's contents. 432 # 433 434 if verbose: 435 print >>sys.stderr, "indexing ..." 436 437 bbdef = {"inspiral": thinca.InspiralCoincDef, "ringdown": ligolw_rinca.RingdownCoincDef}[search] 438 sbdef = {"inspiral": InspiralSICoincDef, "ringdown": RingdownSICoincDef}[search] 439 scndef = {"inspiral": InspiralSCNearCoincDef, "ringdown": RingdownSCNearCoincDef}[search] 440 sngl_type = {"inspiral": lsctables.SnglInspiralTable, "ringdown": lsctables.SnglRingdownTable}[search] 441 sim_type = {"inspiral": lsctables.SimInspiralTable, "ringdown": lsctables.SimRingdownTable}[search] 442 get_sngl_time = {"inspiral": lsctables.SnglInspiral.get_end, "ringdown": lsctables.SnglRingdown.get_start}[search] 443 444 contents = DocContents(xmldoc = xmldoc, bbdef = bbdef, sbdef = sbdef, scndef = scndef, process = process, sngl_type = sngl_type, sim_type = sim_type, get_sngl_time = get_sngl_time) 445 N = len(contents.simtable) 446 447 if isinstance(contents.simtable, lsctables.SimInspiralTable): 448 get_sim_time = lsctables.SimInspiral.get_end 449 elif isinstance(contents.simtable, lsctables.SimRingdownTable): 450 get_sim_time = lsctables.SimRingdown.get_start 451 else: 452 raise ValueError, "Unknown sim table" 453 454 # 455 # Find sim_type <--> sngl_type coincidences. 456 # 457 458 if verbose: 459 print >>sys.stderr, "constructing %s:" % sbdef.description 460 for n, sim in enumerate(contents.simtable): 461 if verbose: 462 print >>sys.stderr, "\t%.1f%%\r" % (100.0 * n / N), 463 types = find_sngl_type_matches(contents, sim, snglcomparefunc, get_sim_time, get_sngl_time) 464 if types: 465 add_sim_type_coinc(contents, sim, types) 466 if verbose: 467 print >>sys.stderr, "\t100.0%" 468 469 # 470 # Find sim_type <--> coinc_event coincidences. 471 # 472 473 if contents.scn_coinc_def_id: 474 if verbose: 475 print >>sys.stderr, "constructing %s:" % (scndef.description) 476 for n, sim in enumerate(contents.simtable): 477 if verbose: 478 print >>sys.stderr, "\t%.1f%%\r" % (100.0 * n / N), 479 coincs = contents.coincs_near_time(get_sim_time(sim), get_sngl_time) 480 coinc_event_ids = find_near_coinc_matches(coincs, sim, nearcoinccomparefunc, get_sim_time, get_sngl_time) 481 if coinc_event_ids: 482 add_sim_coinc_coinc(contents, sim, coinc_event_ids, contents.scn_coinc_def_id) 483 if verbose: 484 print >>sys.stderr, "\t100.0%" 485 486 # 487 # Find sim_type <--> sim_type mappings 488 # if both a sim_inspiral and sim_ring table 489 # are present. 490 # 491 492 if verbose: 493 print >>sys.stderr, "Checking for both SimInspiralTable and SimRingdownTable" 494 if isinstance(contents.simtable, lsctables.SimRingdownTable): 495 try: 496 insp_simtable = lsctables.SimInspiralTable.get_table(xmldoc) 497 except ValueError: 498 print >>sys.stderr,"No SimInspiralTable, only SimRingdownTable present" 499 insp_simtable = None 500 if insp_simtable is not None: 501 if verbose: 502 print >> sys.stderr, "found SimInspiralTable, creating maps to SimRingdownTable" 503 # create an index of the sim_ringdown geocent_start_times 504 sim_ring_time_map = dict([ [(str(row.process_id), row.geocent_start_time), str(row.simulation_id)] for row in contents.simtable]) 505 # create an index of the sim_ringdown's coinc_event_ids 506 sim_ring_ceid_map = dict([ [str(row.event_id), row.coinc_event_id] for row in contents.coincmaptable if row.table_name == "sim_ringdown" ]) 507 # cycle over the sim_inspiral table, creating the maps 508 for this_sim_insp in insp_simtable: 509 if (str(this_sim_insp.process_id), this_sim_insp.geocent_end_time) in sim_ring_time_map: 510 this_sim_ring_id = sim_ring_time_map[( str(this_sim_insp.process_id), this_sim_insp.geocent_end_time )] 511 else: 512 continue 513 if str(this_sim_ring_id) in sim_ring_ceid_map: 514 this_ceid = sim_ring_ceid_map[ str(this_sim_ring_id) ] 515 else: 516 continue 517 # add to the coinc_event_map table 518 new_mapping = lsctables.CoincMap() 519 new_mapping.table_name = "sim_inspiral" 520 new_mapping.event_id = this_sim_insp.simulation_id 521 new_mapping.coinc_event_id = this_ceid 522 contents.coincmaptable.append(new_mapping) 523 524 # 525 # Restore the original event order. 526 # 527 528 if verbose: 529 print >>sys.stderr, "finishing ..." 530 contents.sort_triggers_by_id() 531 532 # 533 # Done. 534 # 535 536 return xmldoc

537

Source Code for Module pylal.lalapps_cbc_injfind