pylal.KW_veto

1 #!/usr/bin/env python 2 # 3 # Copyright (C) 2009 Tomoki Isogai 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 %prog --channel_name=channel_name [--segment_file=File | --gps_start_time=GPSTime --gps_end_time=GPSTime] [options] 21 22 Tomoki Isogai (isogait@carleton.edu) 23 24 This program gets all the KW triggers that fall into a specified segment list above a certain threshold for the channels specified. 25 Supported ifos are H1, H2, L1, and V1. 26 If you are running at the right cluster specified below, the code will use the following path as KW daily dump directory as default. 27 28 For S5 LIGO triggers: 29 /archive/home/lindy/public_html/triggers/s5/ at CIT 30 31 For post S5, trigger locations are: 32 H1: 33 /archive/home/lindy/public_html/triggers/s6/ at LHO 34 H2: 35 /archive/home/lindy/public_html/triggers/s6/ at LLO 36 V1: 37 /archive/home/mabizoua/public_html/KW/ at CIT 38 39 If you want to use triggers at other location, you can specify --KW_location. 40 41 For post S5 triggers, channel name must follow the notation: 42 (ifo)_(channel name)_(min freq)_(max freq) in capital letters. 43 For example, 44 H1_LSC-DARM_ERR_64_1024 45 V1_Pr_B1_ACp_40_1250 46 47 For S5 LIGO triggers, channel name follows the notation: 48 s5_(ifo)_(channel name) in small letters. 49 For example, 50 s5_l1_pobi 51 52 You can omit the frequency part if there is no ambiguity, but the code will give you an error if there are several possibilities. 53 "ls" the above directory to see channel names and max/min frequency available. 54 55 If --out_format is given, a file will be created for each channel with the name 56 (channel)-(start_time)-(duration)_KWtrigs.(specified extention) 57 If --name_tag is given, (name_tag)_ will be added to the name as prefix. 58 If --out_format is not given, the code prints out the result in stdout. 59 You can specify --order_by to sort the output. Supported options are 'GPSTime asc' for ascending time, 'GPSTime desc' for descending time, 'KWSignificance asc' for ascending KW Significance, and 'KWSignificance desc' for descending KW Significance. Default is 'GPSTime asc'. 60 """ 61 62 # ============================================================================= 63 # 64 # PREAMBLE 65 # 66 # ============================================================================= 67 68 69 from __future__ import division 70 71 import os 72 import sys 73 import shutil 74 import re 75 import optparse 76 77 import sqlite3 78 79 from glue.segments import segment, segmentlist 80 from glue import segmentsUtils 81 82 from pylal import git_version 83 from pylal import KW_veto_utils as utils 84 85 __author__ = "Tomoki Isogai <isogait@carleton.edu>" 86 __date__ = "7/10/2009" 87 __version__ = "2.0" 88

89 -def parse_commandline():

90 """ 91 Parse the options given on the command-line. 92 """ 93 parser = optparse.OptionParser(usage = __doc__,version=git_version.verbose_msg) 94 95 parser.add_option("-K", "--KW_location", default=None, 96 help="Location of KW trigger folder if you are not using the folder specified in --help.") 97 parser.add_option("-c", "--channel_name", action="append", default=[], 98 help="Channel names you want KW triggers from. See --help for the naming format. Can be provided multiple times to specify more than one channel. At least one channel is required.") 99 parser.add_option("-m", "--min_thresh", type="int", default=0, 100 help="Minimum KW significance threshold for KW triggers. (Default: 0)") 101 parser.add_option("-S", "--segment_file", default=None, 102 help="Segment file on which KW triggers are retrieved. This option or --gps_start_time and --gps_end_time are required.") 103 parser.add_option("-s", "--gps_start_time", type="int", 104 help="GPS start time on which the KW triggers are retrieved. Required unless --segment_file is specified.") 105 parser.add_option("-e", "--gps_end_time", type="int", 106 help="GPS end time on which KW triggers are retrieved. Required unless --segment_file is specified.") 107 parser.add_option("-n", "--name_tag", default=None, 108 help="If given, this will be added as prefix in the output file name.") 109 parser.add_option("-f", "--out_format", default='stdout', 110 help="Save in this format if specified, otherwise output in stdout. See --help for the supported format and output file name.") 111 parser.add_option("-o","--order_by",default="GPSTime asc", 112 help="Order of the output. See -help for supported order. (Default: GPSTime asc)") 113 parser.add_option("-l", "--scratch_dir", default=".", 114 help="Scratch directory to be used for database engine. Specify local scratch directory for better performance and less fileserver load.") 115 parser.add_option("-v", "--verbose", action="store_true", default=False, 116 help="Run verbosely") 117 118 opts, args = parser.parse_args() 119 120 ############################ Sanity Checks ################################ 121 122 # at least 1 channel is necessary 123 if len(opts.channel_name) < 1: 124 parser.error("Error: at least 1 channel for --channel_name is required.") 125 126 # check if necessary input exists 127 if opts.segment_file is not None: 128 if not os.path.isfile(opts.segment_file): 129 parser.error("Error: segment file %s not found"%opts.segment_file) 130 else: 131 for t in ('gps_start_time','gps_end_time'): 132 if getattr(opts,t) is None: 133 parser.error("Error: either --segment_file or --gps_start/end_time must be specified.") 134 135 if opts.channel_name == []: # 136 parser.error("Error: --channel_name is a required parameter") 137 138 # check if scratch directory exists 139 if not os.path.exists(opts.scratch_dir): 140 parser.error("Error: %s does not exist"%opts.scratch_dir) 141 142 # standardize the format (starting from period .) and see it's supported 143 if opts.out_format[0] != ".": 144 opts.out_format = "." + opts.out_format 145 if opts.out_format not in (\ 146 ".stdout",".pickle",".pickle.gz",".mat",".txt",".txt.gz",".db"): 147 parser.error("Error: %s is not a supported format"%opts.out_format) 148 149 # check order_by is valid and avoid SQL injection attack 150 if opts.order_by not in ("GPSTime asc","GPSTime desc","KWSignificance asc","KWSignificance desc"): 151 parser.error("Error: %s is not valid. See -help for the supported order."%opts.order_by) 152 153 # show parameters 154 if opts.verbose: 155 print >> sys.stderr, "running KW_veto_getKW.py..." 156 print >> sys.stderr, git_version.verbose_msg 157 print >> sys.stderr, "" 158 print >> sys.stderr, "******************* PARAMETERS ********************" 159 print >> sys.stderr, "KW trigger directory:" 160 if opts.KW_location is None: 161 print >> sys.stderr, "/archive/home/lindy/putlic_html/triggers/s5 for S5 LIGO triggers" 162 print >> sys.stderr, "/archive/home/lindy/public_html/triggers/s6 for post S5 LIGO triggers" 163 print >> sys.stderr, "/archive/home/mabizoua/public_html/KW for Virgo triggers" 164 else: 165 print >> sys.stderr, opts.KW_location; 166 print >> sys.stderr,'channels:' 167 print >> sys.stderr, opts.channel_name; 168 print >> sys.stderr,'minimum threshold:' 169 print >> sys.stderr, opts.min_thresh; 170 if opts.segment_file is not None: 171 print >> sys.stderr, 'segment file:' 172 print >> sys.stderr, opts.segment_file; 173 else: 174 print >> sys.stderr, 'start/end GPS time:' 175 print >> sys.stderr, "%d - %d"%(opts.gps_start_time,opts.gps_end_time) 176 print >> sys.stderr, 'name tag:' 177 print >> sys.stderr, opts.name_tag; 178 print >> sys.stderr,'output format:' 179 print >> sys.stderr, opts.out_format[1:]; 180 print >> sys.stderr, 'order by:' 181 print >> sys.stderr, opts.order_by; 182 print >> sys.stderr,'scratch directory:' 183 print >> sys.stderr, opts.scratch_dir; 184 print >> sys.stderr, "" 185 186 return opts

187

188 -def find_file_from_channel(channel,daily_dir):

189 """ 190 From the files in daily_dir, find the file that contains KW triggers 191 for the channel, and return the full channels name and the file path. 192 """ 193 # this function standardize the channel names 194 def norm(channel_name): 195 return channel_name.upper().replace("-","_")

196 197 all_files = \ 198 [f for f in os.listdir(daily_dir) if os.path.splitext(f)[1] == '.trg'] 199 candidates = [f for f in all_files if re.match(norm(channel), norm(f)) != None] 200 if len(candidates) == 1: 201 full_name = norm(os.path.splitext(candidates[0])[0]) 202 return os.path.join(daily_dir,candidates[0]), full_name 203 elif len(candidates) < 1: 204 print >> sys.stderr, "Warning: no match found for %s in %s. See --help for channel name format. Attempting to ignore..."%(channel, daily_dir) 205 return None, None 206 # When there are more than two possibilities see if one of their name 207 # mathes perfectly. If so, warn user and use it, otherwise give an error 208 # and show the possibilities. 209 else: 210 refined_candidates = [f for f in candidates if re.match(norm(os.path.splitext(f)[0]),norm(channel)) != None] 211 if len(refined_candidates) == 1: 212 print >> sys.stderr, """ 213 Warning: Multiple possible files with the channel name %s: 214 %s 215 Using %s and ignoring... 216 """%(channel,", ".join(candidates),refined_candidates[0]) 217 full_name = norm(os.path.splitext(refined_candidates[0])[0]) 218 return os.path.join(daily_dir,refined_candidates[0]), full_name 219 else: 220 print >> sys.stderr, """ 221 Error: Multiple possible files with the channel name %s: 222 %s 223 Please change the channels name so that it is unique. 224 See --help for channel name format. 225 """%(channel, ", ".join(candidates)) 226 sys.exit(1) 227

228 -def get_trigs(channel, segs, min_thresh, trigs_loc=None,name_tag=None,\ 229 scratch_dir=".",verbose=True):

230 """ 231 Get time and KW significance of KW triggers for a particular channel that 232 occured in the specified segments and above specified KW significance 233 threshold. 234 ifo has to be one of H1, H2, L1, V1. 235 """ 236 if verbose: print >> sys.stderr, "getting data for %s..."%channel 237 238 ## initialize SQLite database 239 start_time = segs[0][0] 240 end_time = segs[-1][1] 241 duration = end_time - start_time 242 prefix = "%s-%d-%d-KWtrigs"%(channel, start_time, duration) 243 if name_tag != None: 244 prefix = name_tag + "_" + prefix 245 dbname = prefix+".db" 246 # if the name already exists, rename the old one to avoid collision 247 utils.rename(dbname) 248 249 global KW_working_filename # so that it can be erased when error occurs 250 KW_working_filename = utils.get_connection_filename(\ 251 dbname,tmp_path=scratch_dir,verbose=verbose) 252 253 KWconnection = sqlite3.connect(KW_working_filename) 254 KWcursor = KWconnection.cursor() 255 256 ## create a table for retrieved triggers 257 KWcursor.execute('create table KWtrigs (GPSTime double, KWSignificance double, frequency int)') 258 259 ## determine the KW trigger file we need 260 ifo = channel.split("_")[0].upper() 261 262 # Virgo case 263 if trigs_loc is None and ifo == "V1": 264 trigs_loc = "/archive/home/mabizoua/public_html/KW/" 265 # LIGO case 266 elif trigs_loc is None and ifo in ("H0","H1","H2","L0","L1"): 267 trigs_loc = "/archive/home/lindy/public_html/triggers/s6/" 268 # for S5, first two letters were not ifo but 's5' 269 elif trigs_loc is None and ifo == "S5": 270 trigs_loc = "/archive/home/lindy/public_html/triggers/s5/" 271 elif trigs_loc is None: 272 print >> sys.stderr, "Error: channel name %s is unsupported. See --help for name format."%channel 273 sys.exit(1) 274 275 # sanity check 276 if not os.path.exists(trigs_loc): 277 print >> sys.stderr, "Error: KW daily dump %s not found."%trigs_loc 278 279 ## select daily dump directories in the folder 280 daily_dump_dirs = \ 281 [d for d in os.listdir(trigs_loc) if \ 282 re.match(r"(?P<start>\d{9,10}?)_(?P<end>\d{9,10}?)",d) != None \ 283 and len(d) < 22] 284 # sort by time 285 daily_dump_dirs.sort() 286 287 # get the necessary daily dump folder 288 analyze_dumps = [d for d in daily_dump_dirs \ 289 if int(d.split("_")[0]) < end_time and \ 290 int(d.split("_")[1]) > start_time] 291 292 # for S5, some channels are not in certain dumps: do a little trick to keep 293 # the name 294 full_channel_name = "" 295 for dump in analyze_dumps: 296 trigs_file, tmp_name = \ 297 find_file_from_channel(channel,os.path.join(trigs_loc,dump)) 298 if full_channel_name == "" and tmp_name != None: 299 full_channel_name = tmp_name 300 if trigs_file != None: 301 if verbose: print "retreiving data from %s..."%trigs_file 302 303 for line in open(trigs_file): 304 # get central time and KW significance 305 trig = line.split() 306 t = float(trig[2]) 307 s = float(trig[7]) 308 f = float(trig[3]) 309 310 # check if KW trig is in the given segment and if its significance 311 # is above the minimum specified 312 # FIXME: check inf/nan just in case 313 if t in segs and s > min_thresh: 314 # insert into the database 315 # micro second for GPS time is accurate enough 316 KWcursor.execute("insert into KWtrigs values (?, ?, ?)", ("%.3f"%t, "%.2f"%s,"%d"%f)) 317 318 if full_channel_name == "": # means there is no KW trigger 319 full_channel_name = channel # better than nothing... 320 321 # commit the insertions to the database 322 KWconnection.commit() 323 324 return dbname, KW_working_filename, KWconnection, KWcursor, full_channel_name

325 326 327 # ============================================================================= 328 # 329 # Main 330 # 331 # ============================================================================= 332

333 -def main():

334 ## parse the command line 335 opts = parse_commandline() 336 337 ## get the segments on which triggers are retrieved 338 # case 1: segment file is given 339 if opts.segment_file is not None: 340 # read segment file 341 if opts.segment_file.endswith(".txt"): 342 seg_list = utils.read_segfile(opts.segment_file) 343 elif opts.segment_file.endswith(".xml") or opts.segment_file.endswith(".xml.gz"): 344 seg_list = utils.read_segfile_xml(opts.segment_file,opts.verbose) 345 # case 2: start and end GPS time are given 346 else: 347 seg_list = segmentlist([segment(opts.gps_start_time,opts.gps_end_time)]) 348 349 ## loop over each channels and get KW triggers 350 for chan in opts.channel_name: 351 # wrap in try/finally clause so that the code can erase temporary database 352 # when it encounters an error and had to stop in the middle 353 try: 354 ## get triggers 355 dbname, KW_working_filename, KWconnection, KWcursor, full_chan_name = \ 356 get_trigs(chan,seg_list,opts.min_thresh,trigs_loc=opts.KW_location,name_tag=opts.name_tag,scratch_dir=opts.scratch_dir,verbose=opts.verbose) 357 358 # save/display the result 359 outname = dbname.replace(chan, full_chan_name).replace(".db", opts.out_format) 360 utils.save_db(KWcursor, "KWtrigs", outname, KW_working_filename, 361 order_by=opts.order_by, verbose=opts.verbose) 362 363 # close the connection to database 364 KWconnection.close() 365 366 finally: 367 # erase temporal database 368 if globals().has_key('KW_working_filename'): 369 db = globals()['KW_working_filename'] 370 if opts.verbose: 371 print >> sys.stderr, "removing temporary workspace '%s'..." % db 372 os.remove(db)

373 374 if __name__ == "__main__": 375 main() 376

Source Code for Module pylal.KW_veto_getKW