pylal.ligolw

1 # Copyright (C) 2008--2012 Kipp Cannon, Drew G. Keppel 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 import bisect 28 import math 29 import sys 30 31 32 import lal 33 from glue import iterutils 34 from glue.ligolw import ligolw 35 from glue.ligolw import lsctables 36 from glue.ligolw.utils.coincs import get_coinc_def_id 37 from glue.ligolw.utils import process as ligolw_process 38 from glue.ligolw.utils import segments as ligolw_segments 39 from glue.ligolw.utils import search_summary as ligolw_search_summary 40 from glue import offsetvector 41 from pylal import git_version 42 from pylal import snglcoinc 43 from pylal.xlal import tools as xlaltools 44 from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS 45 from pylal.xlal.datatypes import snglinspiraltable 46 47 48 __author__ = "Kipp Cannon <kipp.cannon@ligo.org>" 49 __version__ = "git id %s" % git_version.id 50 __date__ = git_version.date 51 52 53 # 54 # ============================================================================= 55 # 56 # Speed Hacks 57 # 58 # ============================================================================= 59 # 60 61 62 # 63 # Use C row classes for memory and speed 64 # 65 66 67 lsctables.CoincMapTable.RowType = lsctables.CoincMap = xlaltools.CoincMap

68 69 70 # 71 # Construct a subclass of the C sngl_inspiral row class with the methods 72 # that are needed 73 # 74 75 76 -class SnglInspiral( snglinspiraltable.SnglInspiralTable):

77 __slots__ = () 78

79 - def get_end(self):

80 return LIGOTimeGPS(self.end_time, self.end_time_ns)

81

82 - def set_end(self, gps):

83 self.end_time, self.end_time_ns = gps.seconds, gps.nanoseconds

84

85 - def get_weighted_snr(self, fac):

86 return self.snr

87

88 - def get_new_snr(self, fac):

89 rchisq = self.chisq/(2*self.chisq_dof - 2) 90 nhigh = 2. 91 if rchisq > 1.: 92 return self.snr/((1+rchisq**(fac/nhigh))/2)**(1./fac) 93 else: 94 return self.snr

95

96 - def get_effective_snr(self, fac):

97 rchisq = self.chisq/(2*self.chisq_dof - 2) 98 return self.snr/( (1 + self.snr**2/fac) * rchisq )**(0.25)

99

100 - def get_snr_over_chi(self, fac):

101 return self.snr/self.chisq**(0.5)

102

103 - def __eq__(self, other):

104 return not ( 105 cmp(self.ifo, other.ifo) or 106 cmp(self.end_time, other.end_time) or 107 cmp(self.end_time_ns, other.end_time_ns) or 108 cmp(self.mass1, other.mass1) or 109 cmp(self.mass2, other.mass2) or 110 cmp(self.search, other.search) 111 )

112

113 - def __cmp__(self, other):

114 # compare self's end time to the LIGOTimeGPS instance 115 # other. allows bisection searches by GPS time to find 116 # ranges of triggers quickly 117 return cmp(self.end_time, other.seconds) or cmp(self.end_time_ns, other.nanoseconds)

118 119 120 # 121 # Use C LIGOTimeGPS type 122 # 123 124 125 lsctables.LIGOTimeGPS = LIGOTimeGPS

126 127 128 # 129 # ============================================================================= 130 # 131 # Add Process Information 132 # 133 # ============================================================================= 134 # 135 136 137 -def append_process(xmldoc, **kwargs):

138 process = ligolw_process.append_process( 139 xmldoc, 140 program = u"ligolw_thinca", 141 version = __version__, 142 cvs_repository = u"lscsoft", 143 cvs_entry_time = __date__, 144 comment = kwargs["comment"] 145 ) 146 147 params = [(u"--e-thinca-parameter", u"real_8", kwargs["e_thinca_parameter"])] 148 149 if kwargs["comment"] is not None: 150 params += [(u"--comment", u"lstring", kwargs["comment"])] 151 if kwargs["weighted_snr"] is not None: 152 params += [(u"--weighted-snr", u"lstring", kwargs["weighted_snr"])] 153 if kwargs["magic_number"] is not None: 154 params += [(u"--magic-number", u"real_8", kwargs["magic_number"])] 155 if kwargs["vetoes_name"] is not None: 156 params += [(u"--vetoes-name", u"lstring", kwargs["vetoes_name"])] 157 if kwargs["search_group"] is not None: 158 params += [(u"--search-group", u"lstring", kwargs["search_group"])] 159 if kwargs["trigger_program"] is not None: 160 params += [(u"--trigger-program", u"lstring", kwargs["trigger_program"])] 161 if kwargs["exact_match"] is not None: 162 params += [(u"--exact-match", None, None)] 163 if kwargs["depop_sngl_inspiral"] is not None: 164 params += [(u"--depop-sngl-inspiral", None, None)] 165 if kwargs["drop_veto_info"] is not None: 166 params += [(u"--drop-veto-info,", None,None)] 167 if kwargs["make_expr_tables"] is not None: 168 params += [(u"--make-expr-tables", None, None)] 169 if kwargs["verbose"] is not None: 170 params += [(u"--verbose", None, None)] 171 if kwargs["coinc_end_time_segment"] is not None: 172 params += [(u"--coinc-end-time-segment", u"lstring", kwargs["coinc_end_time_segment"])] 173 174 ligolw_process.append_process_params(xmldoc, process, params) 175 176 return process

177 178 179 # 180 # ============================================================================= 181 # 182 # CoincTables Customizations 183 # 184 # ============================================================================= 185 # 186 187 188 # 189 # The sngl_inspiral <--> sngl_inspiral coinc type. 190 # 191 192 193 InspiralCoincDef = lsctables.CoincDef(search = u"inspiral", search_coinc_type = 0, description = u"sngl_inspiral<-->sngl_inspiral coincidences")

194 195 196 # 197 # Custom snglcoinc.CoincTables subclass. 198 # 199 200 201 -class InspiralCoincTables(snglcoinc.CoincTables):

202 - def __init__(self, xmldoc, vetoes = None, program = u"inspiral", likelihood_func = None, likelihood_params_func = None):

203 snglcoinc.CoincTables.__init__(self, xmldoc) 204 205 # 206 # configure the likelihood ratio evaluator 207 # 208 209 if likelihood_func is None and likelihood_params_func is not None or likelihood_func is not None and likelihood_params_func is None: 210 raise ValueError("must provide both a likelihood function and a parameter function or neither") 211 self.likelihood_func = likelihood_func 212 self.likelihood_params_func = likelihood_params_func 213 214 # 215 # create a string uniquifier 216 # 217 218 self.uniquifier = {} 219 220 # 221 # find the coinc_inspiral table or create one if not found 222 # 223 224 try: 225 self.coinc_inspiral_table = lsctables.table.get_table(xmldoc, lsctables.CoincInspiralTable.tableName) 226 except ValueError: 227 self.coinc_inspiral_table = lsctables.New(lsctables.CoincInspiralTable) 228 xmldoc.childNodes[0].appendChild(self.coinc_inspiral_table) 229 230 # 231 # extract the coalesced out segment lists from the trigger generator 232 # 233 234 self.seglists = ligolw_search_summary.segmentlistdict_fromsearchsummary(xmldoc, program = program).coalesce() 235 if vetoes is not None: 236 self.seglists -= vetoes

237

238 - def append_coinc(self, process_id, time_slide_id, coinc_def_id, events, magic_number):

239 # 240 # populate the coinc_event and coinc_event_map tables 241 # 242 243 coinc = snglcoinc.CoincTables.append_coinc(self, process_id, time_slide_id, coinc_def_id, events) 244 245 # 246 # populate the coinc_inspiral table: 247 # 248 # - end_time is the end time of the first trigger in 249 # alphabetical order by instrument (!?) time-shifted 250 # according to the coinc's offset vector 251 # - mass is average of total masses 252 # - mchirp is average of mchirps 253 # - snr is root-sum-square of SNRs 254 # - false-alarm rates are blank 255 # 256 257 events = sorted(events, lambda a, b: cmp(a.ifo, b.ifo)) 258 259 coinc_inspiral = self.coinc_inspiral_table.RowType() 260 coinc_inspiral.coinc_event_id = coinc.coinc_event_id 261 coinc_inspiral.mass = sum(event.mass1 + event.mass2 for event in events) / len(events) 262 coinc_inspiral.mchirp = sum(event.mchirp for event in events) / len(events) 263 coinc_inspiral.snr = math.sqrt(sum(event.get_weighted_snr(fac = magic_number)**2 for event in events)) 264 # this will fail if chisq=0 for any trigger and you try to calculate effsnr or snr/chi 265 # if so, choose a different command line option for weighted-snr ! 266 coinc_inspiral.false_alarm_rate = None 267 coinc_inspiral.combined_far = None 268 coinc_inspiral.minimum_duration = min(event.template_duration for event in events) 269 coinc_inspiral.set_end(coinc_inspiral_end_time(events, self.time_slide_index[time_slide_id])) 270 coinc_inspiral.set_ifos(event.ifo for event in events) 271 self.coinc_inspiral_table.append(coinc_inspiral) 272 273 # 274 # record the instruments that were on at the time of the 275 # coinc. note that the start time of the coinc must be 276 # unslid to compare with the instrument segment lists 277 # 278 279 tstart = coinc_inspiral.get_end() 280 instruments = set([event.ifo for event in events]) 281 instruments |= set([instrument for instrument, segs in self.seglists.items() if tstart - self.time_slide_index[time_slide_id][instrument] in segs]) 282 coinc.set_instruments(instruments) 283 284 # 285 # if a likelihood ratio calculator is available, assign a 286 # likelihood ratio to the coinc 287 # 288 289 if self.likelihood_func is not None: 290 coinc.likelihood = self.likelihood_func(self.likelihood_params_func(events, self.time_slide_index[time_slide_id])) 291 292 # 293 # save memory by re-using strings 294 # 295 296 coinc.instruments = self.uniquifier.setdefault(coinc.instruments, coinc.instruments) 297 coinc_inspiral.ifos = self.uniquifier.setdefault(coinc_inspiral.ifos, coinc_inspiral.ifos) 298 299 # 300 # done 301 # 302 303 return coinc

304

305 # 306 # Custom function to compute the coinc_inspiral.end_time 307 # 308 309 -def coinc_inspiral_end_time(events, offset_vector):

310 """ 311 A custom function to compute the end_time of a coinc_inspiral trigger 312 @events: a tuple of sngl_inspiral triggers making up a single 313 coinc_inspiral trigger 314 @offset_vector: a dictionary of offsets to apply to different 315 detectors keyed by detector name 316 """ 317 events = sorted(events, lambda a, b: cmp(a.ifo, b.ifo)) 318 return events[0].get_end() + offset_vector[events[0].ifo]

319

320 321 # 322 # ============================================================================= 323 # 324 # Event List Management 325 # 326 # ============================================================================= 327 # 328 329 330 -class InspiralEventList(snglcoinc.EventList):

331 """ 332 A customization of the EventList class for use with the inspiral 333 search. 334 """

335 - def make_index(self):

336 """ 337 Sort events by end time so that a bisection search can 338 retrieve them. Note that the bisection search relies on 339 the __cmp__() method of the SnglInspiral row class having 340 previously been set to compare the event's end time to a 341 LIGOTimeGPS. 342 """ 343 self.sort(lambda a, b: cmp(a.end_time, b.end_time) or cmp(a.end_time_ns, b.end_time_ns))

344

345 - def set_dt(self, dt):

346 """ 347 If an event's end time differs by more than this many 348 seconds from the end time of another event then it is 349 *impossible* for them to be coincident. 350 """ 351 # add 1% for safety, and pre-convert to LIGOTimeGPS to 352 # avoid doing type conversion in loops 353 self.dt = LIGOTimeGPS(dt * 1.01)

354

355 - def get_coincs(self, event_a, offset_a, light_travel_time, threshold, comparefunc):

356 """ 357 The parameter 'threshold' holds the ethinca parameter (for metric-based coincidence tests) 358 or the time window (for exact match with a fixed time uncertainty window) 359 """ 360 # 361 # event_a's end time, with time shift applied 362 # 363 364 end = event_a.get_end() + offset_a - self.offset 365 366 # 367 # extract the subset of events from this list that pass 368 # coincidence with event_a (use bisection searches for the 369 # minimum and maximum allowed end times to quickly identify 370 # a subset of the full list) 371 # 372 373 return [event_b for event_b in self[bisect.bisect_left(self, end - self.dt) : bisect.bisect_right(self, end + self.dt)] if not comparefunc(event_a, offset_a, event_b, self.offset, light_travel_time, threshold)]

374

375 376 # 377 # ============================================================================= 378 # 379 # Coincidence Tests 380 # 381 # ============================================================================= 382 # 383 384 385 -def inspiral_max_dt(events, e_thinca_parameter):

386 """ 387 Given an e-thinca parameter and a list of sngl_inspiral events, 388 return the greatest \Delta t that can separate two events and they 389 still be considered coincident. 390 """ 391 # for each instrument present in the event list, compute the 392 # largest \Delta t interval for the events from that instrument, 393 # and return the sum of the largest two such \Delta t's. 394 return sum(sorted(max(xlaltools.XLALSnglInspiralTimeError(event, e_thinca_parameter) for event in events if event.ifo == instrument) for instrument in set(event.ifo for event in events))[-2:]) + 2. * lal.REARTH_SI / lal.C_SI

395

396 397 -def inspiral_max_dt_exact(events, e_thinca_parameter):

398 """ 399 Given an e-thinca parameter and a list of sngl_inspiral events, 400 return the greatest \Delta t that can separate two events and they 401 still be considered coincident, *if* I am doing exact match coincidence 402 """ 403 # for each instrument present in the event list, compute the 404 # largest \Delta t interval for the events from that instrument, 405 # and return the sum of the largest two such \Delta t's. 406 return sum(sorted(max( (e_thinca_parameter / event.Gamma0)**0.5 for event in events if event.ifo == instrument) for instrument in set(event.ifo for event in events))[-2:]) + 2. * lal.REARTH_SI / lal.C_SI

407

408 409 410 -def inspiral_coinc_compare(a, offseta, b, offsetb, light_travel_time, e_thinca_parameter):

411 """ 412 Returns False (a & b are coincident) if they pass the ellipsoidal 413 thinca test. Otherwise return True. 414 light_travel_time if given in units of seconds. 415 """ 416 if offseta: a.set_end(a.get_end() + offseta) 417 if offsetb: b.set_end(b.get_end() + offsetb) 418 try: 419 # FIXME: should it be "<" or "<="? 420 coincident = xlaltools.XLALCalculateEThincaParameter(a, b) <= e_thinca_parameter 421 except ValueError: 422 # ethinca test failed to converge == events are not 423 # coincident 424 coincident = False 425 if offseta: a.set_end(a.get_end() - offseta) 426 if offsetb: b.set_end(b.get_end() - offsetb) 427 return not coincident

428

429 430 -def inspiral_coinc_compare_exact(a, offseta, b, offsetb, light_travel_time, e_thinca_parameter):

431 """ 432 Returns False (a & b are coincident) if their component masses and spins 433 are equal and they pass the ellipsoidal thinca test. Otherwise return 434 True. light_travel_time if given in units of seconds. 435 """ 436 if inspiral_compare_masses_spins(a, b): 437 # Calculate metric dependent time-window in both detectors 438 twin_a = (e_thinca_parameter / a.Gamma0)**0.5 439 twin_b = (e_thinca_parameter / b.Gamma0)**0.5 440 return float(abs(a.get_end() + offseta - b.get_end() - offsetb)) > light_travel_time + twin_a + twin_b 441 else: 442 return True

443

444 -def inspiral_coinc_compare_exact_dt(a, offseta, b, offsetb, light_travel_time, delta_t):

445 """ 446 Returns False (a & b are coincident) if their component masses and spins 447 are equal and they have dt < (delta_t+light_travel_time) 448 after offsets are considered. Otherwise return True. light_travel_time 449 and delta_t are given in units of seconds. 450 """ 451 if inspiral_compare_masses_spins(a, b): 452 return float(abs(a.get_end() + offseta - b.get_end() - offsetb)) > light_travel_time + delta_t 453 else: 454 return True

455

456 -def inspiral_compare_masses_spins(a, b):

457 """ 458 Returns True if a and b have identical masses and spins. Returns False 459 if a and b have differing masses and spins. 460 """ 461 # define mchirp, eta tuple 462 a_masses = (a.mchirp, a.eta) 463 b_masses = (b.mchirp, b.eta) 464 if (a_masses != b_masses): 465 return False 466 467 try: 468 # check for spin columns (from events in sngl_inspiral table) 469 a_spins = (a.spin1x, a.spin1y, a.spin1z, a.spin2x, a.spin2y, a.spin2z) 470 b_spins = (b.spin1x, b.spin1y, b.spin1z, b.spin2x, b.spin2y, b.spin2z) 471 except: 472 # use spin correction terms for older templates 473 a_spins = (a.beta, a.chi) 474 b_spins = (b.beta, b.chi) 475 476 if (a_spins == b_spins): 477 return True 478 else: 479 return False

480

481 # 482 # ============================================================================= 483 # 484 # Compare Functions 485 # 486 # ============================================================================= 487 # 488 489 490 -def default_ntuple_comparefunc(events, offset_vector):

491 """ 492 Default ntuple test function. Accept all ntuples. 493 """ 494 return False

495

496 497 # 498 # ============================================================================= 499 # 500 # Library API 501 # 502 # ============================================================================= 503 # 504 505 506 -def replicate_threshold(threshold, instruments):

507 """ 508 From a single threshold and a list of instruments, return a 509 dictionary whose keys are every instrument pair (both orders), and 510 whose values are all the same single threshold. 511 512 Example: 513 514 >>> replicate_threshold(6, ["H1", "H2"]) 515 {("H1", "H2"): 6, ("H2", "H1"): 6} 516 """ 517 instruments = sorted(instruments) 518 thresholds = dict((pair, threshold) for pair in iterutils.choices(instruments, 2)) 519 instruments.reverse() 520 thresholds.update(dict((pair,threshold) for pair in iterutils.choices(instruments, 2))) 521 return thresholds

522

523 -def get_vetoes(xmldoc, vetoes_name, verbose = False):

524 if not ligolw_segments.has_segment_tables(xmldoc): 525 if verbose: 526 print >>sys.stderr, "warning: no segment definitions found, vetoes will not be applied" 527 vetoes = None 528 elif not ligolw_segments.has_segment_tables(xmldoc, name = vetoes_name): 529 if verbose: 530 print >>sys.stderr, "warning: document contains segment definitions but none named \"%s\", vetoes will not be applied" % options.vetoes_name 531 vetoes = None 532 else: 533 vetoes = ligolw_segments.segmenttable_get_by_name(xmldoc, vetoes_name).coalesce() 534 return vetoes

535

536 -def ligolw_thinca( 537 xmldoc, 538 process_id, 539 coinc_definer_row, 540 event_comparefunc, 541 thresholds, 542 ntuple_comparefunc = default_ntuple_comparefunc, 543 magic_number = None, 544 veto_segments = None, 545 trigger_program = u"inspiral", 546 likelihood_func = None, 547 likelihood_params_func = None, 548 verbose = False, 549 max_dt_func = None 550 ):

551 if not max_dt_func: 552 err_msg = "Must supply max_dt_func keyword argument to " 553 err_msg += "ligolw_thinca function." 554 raise ValueError(err_msg) 555 # 556 # prepare the coincidence table interface. 557 # 558 559 if verbose: 560 print >>sys.stderr, "indexing ..." 561 coinc_tables = InspiralCoincTables( 562 xmldoc, 563 vetoes = veto_segments, 564 program = trigger_program, 565 likelihood_func = likelihood_func, 566 likelihood_params_func = likelihood_params_func 567 ) 568 coinc_def_id = get_coinc_def_id( 569 xmldoc, 570 coinc_definer_row.search, 571 coinc_definer_row.search_coinc_type, 572 create_new = True, 573 description = coinc_definer_row.description 574 ) 575 sngl_index = dict((row.event_id, row) for row in lsctables.table.get_table(xmldoc, lsctables.SnglInspiralTable.tableName)) 576 577 # 578 # build the event list accessors, populated with events from those 579 # processes that can participate in a coincidence. apply vetoes by 580 # removing events from the lists that fall in vetoed segments 581 # 582 583 eventlists = snglcoinc.EventListDict(InspiralEventList, lsctables.SnglInspiralTable.get_table(xmldoc)) 584 if veto_segments is not None: 585 for eventlist in eventlists.values(): 586 iterutils.inplace_filter((lambda event: event.ifo not in veto_segments or event.get_end() not in veto_segments[event.ifo]), eventlist) 587 588 # 589 # set the \Delta t parameter on all the event lists 590 # 591 592 max_dt = max_dt_func(lsctables.table.get_table(xmldoc, lsctables.SnglInspiralTable.tableName), thresholds) 593 if verbose: 594 print >>sys.stderr, "event bisection search window will be %.16g s" % max_dt 595 for eventlist in eventlists.values(): 596 eventlist.set_dt(max_dt) 597 598 # 599 # replicate the ethinca parameter for every possible instrument 600 # pair 601 # 602 603 thresholds = replicate_threshold(thresholds, set(eventlists)) 604 605 # 606 # construct offset vector assembly graph 607 # 608 609 time_slide_graph = snglcoinc.TimeSlideGraph(coinc_tables.time_slide_index, verbose = verbose) 610 611 # 612 # retrieve all coincidences, apply the final n-tuple compare func 613 # and record the survivors 614 # 615 616 for node, coinc in time_slide_graph.get_coincs(eventlists, event_comparefunc, thresholds, verbose = verbose): 617 ntuple = tuple(sngl_index[id] for id in coinc) 618 if not ntuple_comparefunc(ntuple, node.offset_vector): 619 coinc_tables.append_coinc( 620 process_id, 621 node.time_slide_id, 622 coinc_def_id, 623 ntuple, 624 magic_number 625 ) 626 627 # 628 # remove time offsets from events 629 # 630 631 del eventlists.offsetvector 632 633 # 634 # done 635 # 636 637 return xmldoc

638

639 640 # 641 # ============================================================================= 642 # 643 # GraceDB Utilities 644 # 645 # ============================================================================= 646 # 647 648 649 # 650 # Device to extract sngl_inspiral coincs from a source XML document tree. 651 # 652 653 654 -class sngl_inspiral_coincs(object):

655 """ 656 Dictionary-like device to extract XML document trees containing 657 individual sngl_inspiral coincs from a source XML document tree 658 containing several. 659 660 An instance of the class is initialized with an XML document tree. 661 The coinc event ID of a sngl_inspiral<-->sngl_inspiral coinc in 662 the document can then be used like a dictionary key to retrieve a 663 newly-constructed XML document containing that coinc by itself. 664 The output document trees are complete, self-describing, documents 665 with all metadata about the event from the source document 666 preserved. 667 668 Example: 669 670 >>> coincs = sngl_inspiral_coincs(xmldoc) 671 >>> print(coincs.coinc_def_id) 672 coinc_definer:coinc_def_id:0 673 >>> coincs.keys() 674 [<glue.ligolw.ilwd.cached_ilwdchar_class object at 0x41a4328>] 675 >>> coinc_id = coincs.keys()[0] 676 >>> print(coinc_id) 677 coinc_event:coinc_event_id:83763 678 >>> coincs[coinc_id].write() 679 <?xml version='1.0' encoding='utf-8'?> 680 <!DOCTYPE LIGO_LW SYSTEM "http://ldas-sw.ligo.caltech.edu/doc/ligolwAPI/html/ligolw_dtd.txt"> 681 <LIGO_LW> 682 <Table Name="process:table"> 683 <Column Type="lstring" Name="process:comment"/> 684 <Column Type="lstring" Name="process:node"/> 685 ... 686 687 The XML documents returned from this class share references to the 688 row objects in the original document. Modifications to the row 689 objects in the tables returned by this class will affect both the 690 original document and all other documents returned by this class. 691 However, each retrieval constructs a new document from scratch, 692 they are not cached nor re-used, therefore this operation can be 693 time consuming if it needs to be performed repeatedly but the table 694 objects and document trees can be edited without affecting each 695 other. 696 697 If the source document is modified after this class has been 698 instantiated, the behaviour is undefined. 699 700 To assist with memory clean-up, it is helpful to invoke the 701 .unlink() method on the XML trees returned by this class when they 702 are no longer needed. 703 """

704 - def __init__(self, xmldoc):

705 """ 706 Initialize an instance of the class. xmldoc is the source 707 XML document tree from which the 708 sngl_inspiral<-->sngl_inspiral coincs will be extracted. 709 """ 710 # 711 # find all tables 712 # 713 714 self.process_table = lsctables.table.get_table(xmldoc, lsctables.ProcessTable.tableName) 715 self.process_params_table = lsctables.table.get_table(xmldoc, lsctables.ProcessParamsTable.tableName) 716 self.search_summary_table = lsctables.table.get_table(xmldoc, lsctables.SearchSummaryTable.tableName) 717 self.sngl_inspiral_table = lsctables.table.get_table(xmldoc, lsctables.SnglInspiralTable.tableName) 718 self.coinc_def_table = lsctables.table.get_table(xmldoc, lsctables.CoincDefTable.tableName) 719 self.coinc_event_table = lsctables.table.get_table(xmldoc, lsctables.CoincTable.tableName) 720 self.coinc_inspiral_table = lsctables.table.get_table(xmldoc, lsctables.CoincInspiralTable.tableName) 721 self.coinc_event_map_table = lsctables.table.get_table(xmldoc, lsctables.CoincMapTable.tableName) 722 self.time_slide_table = lsctables.table.get_table(xmldoc, lsctables.TimeSlideTable.tableName) 723 724 # 725 # index the process, process params, search_summary, 726 # sngl_inspiral and time_slide tables 727 # 728 729 self.process_index = dict((row.process_id, row) for row in self.process_table) 730 self.process_params_index = {} 731 for row in self.process_params_table: 732 self.process_params_index.setdefault(row.process_id, []).append(row) 733 self.search_summary_index = dict((row.process_id, row) for row in self.search_summary_table) 734 self.sngl_inspiral_index = dict((row.event_id, row) for row in self.sngl_inspiral_table) 735 self.time_slide_index = {} 736 for row in self.time_slide_table: 737 self.time_slide_index.setdefault(row.time_slide_id, []).append(row) 738 739 # 740 # find the sngl_inspiral<-->sngl_inspiral coincs 741 # 742 743 self.coinc_def, = (row for row in self.coinc_def_table if row.search == InspiralCoincDef.search and row.search_coinc_type == InspiralCoincDef.search_coinc_type) 744 self.coinc_event_index = dict((row.coinc_event_id, row) for row in self.coinc_event_table if row.coinc_def_id == self.coinc_def.coinc_def_id) 745 self.coinc_inspiral_index = dict((row.coinc_event_id, row) for row in self.coinc_inspiral_table) 746 assert set(self.coinc_event_index) == set(self.coinc_inspiral_index) 747 self.coinc_event_map_index = dict((coinc_event_id, []) for coinc_event_id in self.coinc_event_index) 748 for row in self.coinc_event_map_table: 749 try: 750 self.coinc_event_map_index[row.coinc_event_id].append(row) 751 except KeyError: 752 continue

753 754 @property

755 - def coinc_def_id(self):

756 """ 757 The coinc_def_id of the sngl_inspiral<-->sngl_inspiral 758 coincs in the source XML document. 759 """ 760 return self.coinc_def.coinc_def_id

761

762 - def sngl_inspirals(self, coinc_event_id):

763 """ 764 Return a list of the sngl_inspiral rows that participated 765 in the coincidence given by coinc_event_id. 766 """ 767 return [self.sngl_inspiral_index[event_id] for event_id in self.coinc_event_map_index[coinc_event_id]]

768

769 - def offset_vector(self, time_slide_id):

770 """ 771 Return the offsetvector given by time_slide_id. 772 """ 773 return offsetvector.offsetvector((row.instrument, row.offset) for row in self.time_slide_index[time_slide_id])

774

775 - def __getitem__(self, coinc_event_id):

776 """ 777 Construct and return an XML document containing the 778 sngl_inspiral<-->sngl_inspiral coinc carrying the given 779 coinc_event_id. 780 """ 781 newxmldoc = ligolw.Document() 782 newxmldoc.appendChild(ligolw.LIGO_LW()) 783 784 # when making these, we can't use table.new_from_template() 785 # because we need to ensure we have a Table subclass, not a 786 # DBTable subclass 787 new_process_table = newxmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.ProcessTable, self.process_table.columnnames)) 788 new_process_params_table = newxmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.ProcessParamsTable, self.process_params_table.columnnames)) 789 new_search_summary_table = newxmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.SearchSummaryTable, self.search_summary_table.columnnames)) 790 new_sngl_inspiral_table = newxmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.SnglInspiralTable, self.sngl_inspiral_table.columnnames)) 791 new_coinc_def_table = newxmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.CoincDefTable, self.coinc_def_table.columnnames)) 792 new_coinc_event_table = newxmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.CoincTable, self.coinc_event_table.columnnames)) 793 new_coinc_inspiral_table = newxmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.CoincInspiralTable, self.coinc_inspiral_table.columnnames)) 794 new_coinc_event_map_table = newxmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.CoincMapTable, self.coinc_event_map_table.columnnames)) 795 new_time_slide_table = newxmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.TimeSlideTable, self.time_slide_table.columnnames)) 796 797 new_coinc_def_table.append(self.coinc_def) 798 coinc_event = self.coinc_event_index[coinc_event_id] 799 new_coinc_event_table.append(coinc_event) 800 new_coinc_inspiral_table.append(self.coinc_inspiral_index[coinc_event_id]) 801 map(new_coinc_event_map_table.append, self.coinc_event_map_index[coinc_event_id]) 802 map(new_time_slide_table.append, self.time_slide_index[coinc_event.time_slide_id]) 803 for row in new_coinc_event_map_table: 804 new_sngl_inspiral_table.append(self.sngl_inspiral_index[row.event_id]) 805 806 for process_id in set(new_sngl_inspiral_table.getColumnByName("process_id")) | set(new_coinc_event_table.getColumnByName("process_id")) | set(new_time_slide_table.getColumnByName("process_id")): 807 # process row is required 808 new_process_table.append(self.process_index[process_id]) 809 try: 810 map(new_process_params_table.append, self.process_params_index[process_id]) 811 except KeyError: 812 # process_params rows are optional 813 pass 814 try: 815 new_search_summary_table.append(self.search_summary_index[process_id]) 816 except KeyError: 817 # search_summary rows are optional 818 pass 819 820 return newxmldoc

821

822 - def __iter__(self):

823 """ 824 Iterate over the coinc_event_id's in the source document. 825 """ 826 return iter(self.coinc_event_index)

827

828 - def __nonzero__(self):

829 return bool(self.coinc_event_index)

830

831 - def keys(self):

832 """ 833 A list of the coinc_event_id's of the 834 sngl_inspiral<-->sngl_inspiral coincs available in the 835 source XML document. 836 """ 837 return self.coinc_event_index.keys()

838

839 - def items(self):

840 """ 841 Yield a sequence of (coinc_event_id, XML tree) tuples, one 842 for each sngl_inspiral<-->sngl_inspiral coinc in the source 843 document. 844 845 NOTE: to allow this to work more easily with very large 846 documents, instead of returning the complete sequence as a 847 pre-constructed list this method is implemented as a 848 generator. 849 """ 850 for coinc_event_id in self: 851 yield (coinc_event_id, self[coinc_event_id])

852 853 iteritems = items 854

855 - def column_index(self, table_name, column_name):

856 """ 857 Return a dictionary mapping coinc_event_id to the values in 858 the given column in the given table. 859 860 Example: 861 862 >>> print(coincs.column_index("coinc_event", "likelihood")) 863 864 Only columns in the coinc_event and coinc_inspiral tables 865 can be retrieved this way. 866 """ 867 if not lsctables.table.CompareTableNames(table_name, lsctables.CoincTable.tableName): 868 return dict(zip(self.coinc_event_table.getColumnByName("coinc_event_id"), self.coinc_event_table.getColumnByName(column_name))) 869 elif not lsctables.table.CompareTableNames(table_name, lsctables.CoincInspiralTable.tableName): 870 return dict(zip(self.coinc_inspiral_table.getColumnByName("coinc_event_id"), self.coinc_inspiral_table.getColumnByName(column_name))) 871 raise ValueError(table_name)

872

Source Code for Module pylal.ligolw_sstinca