pylal.dq.dqSegmentUtils

1 #!/usr/bin/env python 2 3 # Copyright (C) 2011 Duncan Macleod 4 # 5 # This program is free software; you can redistribute it and/or modify it 6 # under the terms of the GNU General Public License as published by the 7 # Free Software Foundation; either version 3 of the License, or (at your 8 # option) any later version. 9 # 10 # This program is distributed in the hope that it will be useful, but 11 # WITHOUT ANY WARRANTY; without even the implied warranty of 12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 13 # Public License for more details. 14 # 15 # You should have received a copy of the GNU General Public License along 16 # with this program; if not, write to the Free Software Foundation, Inc., 17 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 18 19 # ============================================================================== 20 # Preamble 21 # ============================================================================== 22 23 import os, sys, re, operator, math 24 from StringIO import StringIO 25 from glue import segments 26 from glue.ligolw import ligolw, lsctables, table, utils 27 from glue.ligolw.utils import segments as ligolw_segments 28 from glue.segmentdb import query_engine, segmentdb_utils 29 from glue.ligolw.utils import process as ligolw_process 30 from pylal.dq.dqTriggerUtils import def_get_time 31 32 from scipy.stats import poisson 33 34 LIGOTimeGPS = lsctables.LIGOTimeGPS 35 36 # Some boilerplate to make segmentlists picklable 37 import copy_reg 38 copy_reg.pickle(type(segments.segment(0, 1)), \ 39 lambda x:(segments.segment, (x[0], x[1]))) 40 copy_reg.pickle(type(segments.segmentlist([])), 41 lambda x:(segments.segmentlist, ([y for y in x], ))) 42 43 from glue import git_version 44 45 __author__ = "Andrew P Lundgren <andrew.lundgren@ligo.org>, Duncan Macleod <duncan.macleod@ligo.org>" 46 __version__ = "git id %s" % git_version.id 47 __date__ = git_version.date 48 49 """ 50 This module provides useful segment and veto tools for data quality investigations. 51 """ 52 53 # ============================================================================== 54 # Function to load segments from an xml file 55 # ============================================================================== 56

57 -def fromsegmentxml(file, dict=False, id=None):

58 59 """ 60 Read a glue.segments.segmentlist from the file object file containing an 61 xml segment table. 62 63 Arguments: 64 65 file : file object 66 file object for segment xml file 67 68 Keyword Arguments: 69 70 dict : [ True | False ] 71 returns a glue.segments.segmentlistdict containing coalesced 72 glue.segments.segmentlists keyed by seg_def.name for each entry in the 73 contained segment_def_table. Default False 74 id : int 75 returns a glue.segments.segmentlist object containing only those 76 segments matching the given segment_def_id integer 77 78 """ 79 80 # load xmldocument and SegmentDefTable and SegmentTables 81 xmldoc, digest = utils.load_fileobj(file, gz=file.name.endswith(".gz"), contenthandler = lsctables.use_in(ligolw.LIGOLWContentHandler)) 82 seg_def_table = lsctables.SegmentDefTable.get_table(xmldoc) 83 seg_table = lsctables.SegmentTable.get_table(xmldoc) 84 85 if dict: 86 segs = segments.segmentlistdict() 87 else: 88 segs = segments.segmentlist() 89 90 seg_id = {} 91 for seg_def in seg_def_table: 92 seg_id[int(seg_def.segment_def_id)] = str(seg_def.name) 93 if dict: 94 segs[str(seg_def.name)] = segments.segmentlist() 95 96 for seg in seg_table: 97 if dict: 98 segs[seg_id[int(seg.segment_def_id)]]\ 99 .append(segments.segment(seg.start_time, seg.end_time)) 100 continue 101 if id and int(seg.segment_def_id)==id: 102 segs.append(segments.segment(seg.start_time, seg.end_time)) 103 continue 104 segs.append(segments.segment(seg.start_time, seg.end_time)) 105 106 if dict: 107 for seg_name in seg_id.values(): 108 segs[seg_name] = segs[seg_name].coalesce() 109 else: 110 segs = segs.coalesce() 111 112 xmldoc.unlink() 113 114 return segs

115 116 # ============================================================================== 117 # Write to segment xml file 118 # ============================================================================== 119

120 -def tosegmentxml(file, segs):

121 122 """ 123 Write the glue.segments.segmentlist object segs to file object file in xml 124 format with appropriate tables. 125 """ 126 127 # generate empty document 128 xmldoc = ligolw.Document() 129 xmldoc.appendChild(ligolw.LIGO_LW()) 130 xmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.ProcessTable)) 131 xmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.ProcessParamsTable)) 132 133 # append process to table 134 process = ligolw_process.append_process(xmldoc,\ 135 program='pylal.dq.dqSegmentUtils',\ 136 version=__version__,\ 137 cvs_repository='lscsoft',\ 138 cvs_entry_time=__date__) 139 140 gpssegs = segments.segmentlist() 141 for seg in segs: 142 gpssegs.append(segments.segment(LIGOTimeGPS(seg[0]), LIGOTimeGPS(seg[1]))) 143 144 # append segs and seg definer 145 segments_tables = ligolw_segments.LigolwSegments(xmldoc) 146 segments_tables.add(ligolw_segments.LigolwSegmentList(active=gpssegs)) 147 # finalise 148 segments_tables.coalesce() 149 segments_tables.optimize() 150 segments_tables.finalize(process) 151 ligolw_process.set_process_end_time(process) 152 153 # write file 154 with utils.SignalsTrap(): 155 utils.write_fileobj(xmldoc, file, gz=False)

156 157 # ============================================================================== 158 # Function to load segments from a csv file 159 # ============================================================================== 160

161 -def fromsegmentcsv(csvfile):

162 163 """ 164 Read a glue.segments.segmentlist object from the file object file containin 165 a comma separated list of segments. 166 """ 167 168 def CSVLineToSeg(line): 169 tstart, tend = map(LIGOTimeGPS, line.split(',')) 170 return segments.segment(tstart, tend)

171 172 segs = segments.segmentlist([CSVLineToSeg(line) for line in csvfile]) 173 174 return segs.coalesce() 175 176 # ============================================================================== 177 # Function to parse a segment list for CBC analysable segments 178 # ============================================================================== 179

180 -def CBCAnalyzableSegs(seglist):

181 182 """ 183 Remove any segments shorter than 2064 seconds from seglist because ihope 184 won't analyze them. 185 """ 186 187 return segments.segmentlist([seg for seg in seglist if abs(seg) >= 2064])

188 189 # ============================================================================== 190 # Function to pad a list of segments given start and end paddings 191 # ============================================================================== 192

193 -def pad_segmentlist(seglist, start_pad, end_pad=None):

194 195 """ 196 Given a veto segmentlist, start pad, and end pad, pads and coalesces the 197 segments. Convention is to expand segments with positive padding, contract 198 with negative. Any segments that are not big enough to be contracted 199 appropriately are removed. 200 """ 201 202 if not end_pad: end_pad = start_pad 203 204 padded = lambda seg: segments.segment(seg[0]-start_pad, seg[1]+end_pad) 205 206 seglist = segments.segmentlist([padded(seg) for seg in seglist if \ 207 abs(seg) > start_pad+end_pad]) 208 209 return seglist.coalesce()

210 211 # ============================================================================== 212 # Function to crop a list of segments 213 # ============================================================================== 214

215 -def crop_segmentlist(seglist, end_chop=30):

216 217 """ 218 Given a segmentlist and time to chop, removes time from the end of each 219 segment (defaults to 30 seconds). 220 """ 221 222 chopped = lambda seg: segments.segment(seg[0], max(seg[0], seg[1] - end_chop)) 223 seglist = segments.segmentlist([chopped(seg) for seg in seglist]) 224 return seglist.coalesce()

225 226 # ============================================================================= 227 # Function to return segments in given gps time range 228 # ============================================================================= 229

230 -def grab_segments(start, end, flag,\ 231 segment_url='https://segdb.ligo.caltech.edu',\ 232 segment_summary=False):

233 234 """ 235 Returns a segmentlist containing the segments during which the given flag 236 was active in the given period. 237 238 Arguments: 239 240 start : int 241 GPS start time 242 end : int 243 GPS end time 244 flag : string 245 'IFO:NAME:VERSION' format string 246 247 Keyword arguments: 248 249 segment_url : string 250 url of segment database to query, default https://segdb.ligo.caltech.edu 251 segment_summary : [ True | False ] 252 also return the glue.segments.segmentlist defining the valid span of the 253 returned segments 254 """ 255 256 # set times 257 start = int(math.floor(start)) 258 end = int(math.ceil(end)) 259 260 # set query engine 261 connection = segmentdb_utils.setup_database(segment_url) 262 engine = query_engine.LdbdQueryEngine(connection) 263 264 # format flag name 265 if isinstance(flag, basestring): 266 flags = flag.split(',') 267 else: 268 flags = flag 269 270 segdefs = [] 271 for f in flags: 272 spec = f.split(':') 273 if len(spec) < 2 or len(spec) > 3: 274 raise AttributeError, "Included segements must be of the form "+\ 275 "ifo:name:version or ifo:name:*" 276 277 ifo = spec[0] 278 name = spec[1] 279 280 if len(spec) is 3 and spec[2] is not '*': 281 version = int(spec[2]) 282 if version < 1: 283 raise AttributeError, "Segment version numbers must be greater than zero" 284 else: 285 version = '*' 286 287 # expand segment definer 288 segdefs += segmentdb_utils.expand_version_number(engine, (ifo, name, version, \ 289 start, end, 0, 0)) 290 291 # query database and return 292 segs = segmentdb_utils.query_segments(engine, 'segment', segdefs) 293 segs = [s.coalesce() for s in segs] 294 if segment_summary: 295 segsums = segmentdb_utils.query_segments(engine, 'segment_summary', segdefs) 296 #segsums = reduce(operator.or_, segsums).coalesce() 297 segsums = [s.coalesce() for s in segsums] 298 segsummap = [segments.segmentlist() for f in flags] 299 for segdef,segsum in zip(segdefs, segsums): 300 try: 301 fidx = flags.index(':'.join(map(str, segdef[:3]))) 302 except ValueError: 303 fidx = flags.index(':'.join(segdef[:2])) 304 segsummap[fidx].extend(segsum) 305 if flag == flags[0]: 306 return segs[0],segsummap[0] 307 else: 308 return segs,segsummap 309 if flag == flags[0]: 310 return segs[0] 311 else: 312 return segs

313 314 # ============================================================================== 315 # Dump flags from segment database 316 # ============================================================================== 317

318 -def dump_flags(ifos=None, segment_url=None, match=None, unmatch=None,\ 319 latest=False):

320 321 """ 322 Returns the list of all flags defined in the database. 323 324 Keyword rguments: 325 ifo : [ str | list ] 326 list of ifos to query, or str for single ifo 327 segment_url : str 328 url of segment database, defaults to contents of S6_SEGMENT_SERVER 329 environment variable 330 match : [ str | regular pattern ] 331 regular expression to search against returned flag names, e.g, 'UPV' 332 unmatch : [ str | regular pattern ] 333 regular expression to negatively search against returned flag names 334 """ 335 336 if isinstance(ifos, str): 337 ifos = [ifos] 338 339 # get url 340 if not segment_url: 341 segment_url = os.getenv('S6_SEGMENT_SERVER') 342 343 # open connection to LDBD(W)Server 344 myClient = segmentdb_utils.setup_database(segment_url) 345 346 reply = StringIO(myClient.query(squery)) 347 xmldoc, digest = utils.load_fileobj(reply) 348 seg_def_table = lsctables.SegmentDefTable.get_table(xmldoc) 349 350 # sort table by ifo, name and version 351 seg_def_table.sort(key=lambda flag: (flag.ifos[0], flag.name, \ 352 flag.version), reverse=True) 353 354 flags = lsctables.New(type(seg_def_table)) 355 356 for row in seg_def_table: 357 358 # test re match 359 if match and not re.search(match, row.name): continue 360 361 # test re unmatch 362 if unmatch and re.search(unmatch, row.name): continue 363 364 # only append latest versions of multiple flags 365 flatest=True 366 if latest: 367 # get all flags with same ifo and name 368 vflags = [f for f in flags if row.name==f.name and\ 369 row.get_ifos()==f.get_ifos()] 370 # if later version in list, move on 371 for f in vflags: 372 if f.version>=row.version: 373 flatest=False 374 break 375 if not flatest: 376 continue 377 378 # append those flags matching ifos requirement 379 for ifo in ifos: 380 if ifo in row.get_ifos(): 381 flags.append(row) 382 break 383 384 return flags

385 386 # ============================================================================== 387 # Calculate Poisson safety 388 # ============================================================================== 389

390 -def poisson_safety(segs, injTable, livetime):

391 392 """ 393 Return a tuple containing the number of vetoed injections, the number 394 expected, and the Poisson safety probability based on the number of 395 injections vetoed relative to random chance according to Poisson statistics. 396 397 Arguments: 398 399 segs : glue.segments.segmentlist 400 list of segments to be tested 401 injTable : glue.ligolw.table.Table 402 table of injections 403 livetime : [ float ] 404 livetime of search 405 """ 406 407 deadtime = segs.__abs__() 408 409 get_time = def_get_time(injTable.tableName) 410 injvetoed = len([ inj for inj in injTable if get_time(inj) in segs ]) 411 412 injexp = len(injTable)*float(deadtime) / float(livetime) 413 414 prob = 1 - poisson.cdf(injvetoed-1, injexp) 415 416 return injvetoed, injexp, prob

417 418 # ============================================================================= 419 # Convert a data quality bit mask into segments 420 # ============================================================================= 421

422 -def _bits(i, n=8):

423 """ 424 Convert integer bit mask into binary bits. Returns a list of 0s or 1s 425 from 2^0 up to 2^n. 426 427 Example: 428 429 >>> _bits(295, n=8) 430 [1, 1, 1, 0, 0, 1, 0, 0] 431 """ 432 return [(0, 1)[int(i)>>j & 1] for j in xrange(n)]

433

434 -def DQSegments(time, data, dq_key):

435 436 """ 437 Returns a glue.segments.segmentlistdict of active segments for each bit 438 in a dq_key. 439 """ 440 441 segdict = segments.segmentlistdict() 442 for key in dq_key: 443 segdict[key] = segments.segmentlist() 444 445 # convert DQ bits into segments 446 for i in xrange(len(data)): 447 binary = _bits(data[i], len(dq_key)) 448 seg = segments.segment(time[i], time[i]-1) 449 for j, key in enumerate(dq_key): 450 if binary[j] == 1: 451 segdict[key].append(seg) 452 453 segdict = segdict.coalesce() 454 455 return segdict

456

Source Code for Module pylal.dq.dqSegmentUtils