pylal.ligolw_miinjfind

1 # Copyright (C) 2012 Duncan M. Macleod, 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. 29 30 Contains code providing the capacity to search a list of multi_inspiral 31 candidates for events matching entries in a sim_inspiral list of 32 software injections, recording the matches as inspiral <--> injection 33 coincidences using the standard coincidence infrastructure. 34 """ 35 36 import bisect 37 import sys 38 39 from glue.ligolw import table 40 from glue.ligolw import lsctables 41 from glue.ligolw.utils import process as ligolw_process 42 43 from pylal import git_version, llwapp, \ 44 SimInspiralUtils, MultiInspiralUtils 45 from lalburst import timeslides as ligolw_tisi 46 from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS 47 48 49 __author__ = "Duncan M. Macleod <duncan.macleod@ligo.org>" 50 __credits__ = "Kipp Cannon <kipp.cannon@ligo.org>" 51 __version__ = git_version.id 52 __date__ = git_version.date 53 54 55 # 56 # ============================================================================= 57 # 58 # Speed Hacks 59 # 60 # ============================================================================= 61 # 62 63 64 lsctables.LIGOTimeGPS = LIGOTimeGPS 65 66

67 -def multi_inspiral___cmp__(self, other):

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

71 72 lsctables.MultiInspiral.__cmp__ = multi_inspiral___cmp__ 73 74 75 # 76 # ============================================================================= 77 # 78 # Document Interface 79 # 80 # ============================================================================= 81 # 82 83 84 MultiInspiralSICoincDef = lsctables.CoincDef( 85 search=u"inspiral", search_coinc_type=1, 86 description=u"sim_inspiral<-->multi_inspiral coincidences") 87

88 -class DocContents(object):

89 """A wrapper interface to the XML document. 90 """

91 - def __init__(self, xmldoc, sbdef, process):

92 # 93 # store the process row 94 # 95 96 self.process = process 97 98 # 99 # locate the multi_inspiral and sim_inspiral tables 100 # 101 102 self.multiinspiraltable = lsctables.MultiInspiralTable.get_table( 103 xmldoc) 104 self.siminspiraltable = lsctables.SimInspiralTable.get_table( 105 xmldoc) 106 107 # 108 # get out segment lists for programs that generated 109 # triggers (currently only used for time_slide vector 110 # construction) 111 # 112 113 search_summary = lsctables.SearchSummaryTable.get_table(xmldoc) 114 pids = set(self.multiinspiraltable.getColumnByName("process_id")) 115 seglists = search_summary.get_out_segmentlistdict(pids)\ 116 .coalesce() 117 118 # 119 # construct the zero-lag time slide needed to cover the 120 # instruments listed in all the triggers, then determine 121 # its ID (or create it if needed) 122 # 123 # FIXME: in the future, the sim_inspiral table should 124 # indicate time slide at which the injection was done 125 # 126 127 self.tisi_id = ligolw_tisi.get_time_slide_id( 128 xmldoc, 129 dict.fromkeys(seglists, 0.0), 130 create_new=process) 131 132 # 133 # get coinc_definer row for sim_inspiral <--> multi_inspiral 134 # coincs; this creates a coinc_definer table if the 135 # document doesn't have one 136 # 137 138 self.sb_coinc_def_id = llwapp.get_coinc_def_id( 139 xmldoc, 140 sbdef.search, 141 sbdef.search_coinc_type, 142 create_new=True, 143 description=sbdef.description) 144 145 # 146 # get coinc table, create one if needed 147 # 148 149 try: 150 self.coinctable = lsctables.CoincTable.get_table(xmldoc) 151 except ValueError: 152 self.coinctable = lsctables.New(lsctables.CoincTable) 153 xmldoc.childNodes[0].appendChild(self.coinctable) 154 self.coinctable.sync_next_id() 155 156 # 157 # get coinc_map table, create one if needed 158 # 159 160 try: 161 self.coincmaptable = lsctables.CoincMapTable.get_table(xmldoc) 162 except ValueError: 163 self.coincmaptable = lsctables.New(lsctables.CoincMapTable) 164 xmldoc.childNodes[0].appendChild(self.coincmaptable) 165 166 # 167 # sort multi_inspiral table by end time 168 # 169 170 self.multiinspiraltable.sort(lambda a, b: 171 cmp(a.end_time, b.end_time) or 172 cmp(a.end_time_ns, b.end_time_ns))

173

174 - def inspirals_near_endtime(self, t, dt):

175 """Return a list of the inspiral events whose peak times are 176 within self.inspiral_end_time_window of t. 177 """ 178 return self.multiinspiraltable[ 179 bisect.bisect_left(self.multiinspiraltable, t - dt): 180 bisect.bisect_right(self.multiinspiraltable, t + dt)]

181

182 - def sort_triggers_by_id(self):

183 """Sort the multi_inspiral table's rows by ID (tidy-up document 184 for output). 185 """ 186 self.multiinspiraltable.sort(lambda a, b: cmp(a.event_id, b.event_id))

187

188 - def new_coinc(self, coinc_def_id):

189 """Construct a new coinc_event row attached to the given 190 process, and belonging to the set of coincidences defined 191 by the given coinc_def_id. 192 """ 193 coinc = lsctables.Coinc() 194 coinc.process_id = self.process.process_id 195 coinc.coinc_def_id = coinc_def_id 196 coinc.coinc_event_id = self.coinctable.get_next_id() 197 coinc.time_slide_id = self.tisi_id 198 coinc.set_instruments(None) 199 coinc.nevents = 0 200 coinc.likelihood = None 201 self.coinctable.append(coinc) 202 return coinc

203 204 205 # 206 # ============================================================================= 207 # 208 # Add Process Information 209 # 210 # ============================================================================= 211 # 212 213 process_program_name = "ligolw_inspinjfind" 214

215 -def append_process(xmldoc, match_algorithm, time_window, loudest_by, comment):

216 """Convenience wrapper for adding process metadata to the document. 217 """ 218 process = ligolw_process.append_process(xmldoc, 219 program=process_program_name, 220 version=__version__, 221 cvs_repository=u"lscsoft", 222 cvs_entry_time=__date__, comment=comment) 223 params = [(u"--match-algorithm", u"lstring", match_algorithm), 224 (u"--time-window", u"lstring", time_window)] 225 if loudest_by: 226 params.append((u"--loudest-by", u"lstring", loudest_by)) 227 ligolw_process.append_process_params(xmldoc, process, params) 228 return process

229 230 231 # 232 # ============================================================================= 233 # 234 # Build sim_inspiral <--> multi_inspiral Coincidences 235 # 236 # ============================================================================= 237 # 238 239

240 -def find_multi_inspiral_matches(contents, sim, time_window, loudest_by=None):

241 """Scan the inspiral table for triggers matching sim. 242 """ 243 events = [inspiral for inspiral in 244 contents.inspirals_near_endtime(sim.get_end(), time_window)] 245 if loudest_by and len(events): 246 if hasattr(lsctables.MultiInspiral, "get_%s" % loudest_by): 247 rank = getattr(lsctables.MultiInspiral, "get_%s" % loudest_by) 248 else: 249 rank = lambda x: getattr(x, loudest_by) 250 return sorted(events, key=lambda mi: rank(mi))[-1:] 251 else: 252 return events

253 254

255 -def add_sim_inspiral_coinc(contents, sim, inspirals):

256 """Create a coinc_event in the coinc table, and add arcs in the 257 coinc_event_map table linking the sim_inspiral row and the list of 258 multi_inspiral rows to the new coinc_event row. 259 """ 260 coinc = contents.new_coinc(contents.sb_coinc_def_id) 261 if inspirals: 262 coinc.set_instruments( 263 lsctables.instrument_set_from_ifos(inspirals[0].ifos)) 264 coinc.nevents = len(inspirals) 265 coincmap = lsctables.CoincMap() 266 coincmap.coinc_event_id = coinc.coinc_event_id 267 coincmap.table_name = sim.simulation_id.table_name 268 coincmap.event_id = sim.simulation_id 269 contents.coincmaptable.append(coincmap) 270 for event in inspirals: 271 coincmap = lsctables.CoincMap() 272 coincmap.coinc_event_id = coinc.coinc_event_id 273 coincmap.table_name = event.event_id.table_name 274 coincmap.event_id = event.event_id 275 contents.coincmaptable.append(coincmap)

276 277 278 # 279 # ============================================================================= 280 # 281 # Library API 282 # 283 # ============================================================================= 284 # 285 286

287 -def ligolw_miinjfind(xmldoc, process, search, time_window, loudest_by=None, 288 verbose=False):

289 """Parse an XML document and find coincidences between entries 290 in sim_inspiral and multi_inspiral tables. 291 """ 292 if verbose: 293 sys.stderr.write("indexing ...\n") 294 295 sbdef = {"inspiral": MultiInspiralSICoincDef}[search] 296 297 contents = DocContents(xmldoc=xmldoc, sbdef=sbdef, process=process) 298 N = len(contents.siminspiraltable) 299 300 # 301 # Find sim_inspiral <--> multi_inspiral coincidences. 302 # 303 304 if verbose: 305 sys.stderr.write("constructing %s:\n" % sbdef.description) 306 for n, sim in enumerate(contents.siminspiraltable): 307 if verbose: 308 sys.stderr.write("\t%.1f%%\r" % (100.0 * n / N)) 309 sys.stderr.flush() 310 inspirals = find_multi_inspiral_matches(contents, sim, time_window, 311 loudest_by=loudest_by) 312 if inspirals: 313 add_sim_inspiral_coinc(contents, sim, inspirals) 314 if verbose: 315 sys.stderr.write("\t100.0%\n") 316 317 # 318 # Restore the original event order. 319 # 320 321 if verbose: 322 sys.stderr.write("finishing ...\n") 323 contents.sort_triggers_by_id() 324 325 # 326 # Done. 327 # 328 329 return xmldoc

330

Source Code for Module pylal.ligolw_miinjfind