Source Code for Module pylal.hacrutils

  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,  
 17   
 18  """ 
 19  Auxiliary functions for running and postprocessing HACR pipeline code. 
 20  """ 
 21   
 22  from __future__ import division 
 23   
 24  import MySQLdb 
 25  import sys 
 26  import re 
 27  import datetime 
 28   
 29  from lal import LIGOTimeGPS, GPSToUTC 
 30   
 31  from pylal import git_version 
 32   
 33  from glue.ligolw import lsctables 
 34   
 35  __author__  = 'Duncan Macleod <duncan.macleod@ligo.org>' 
 36  __version__ = git_version.id 
 37  __date__    = git_version.date 
 38   
 39  # ============================================================================= 
 40  # convert ascii line to HACR trigger 
 41  # ============================================================================= 
 42   
 43  _comment = re.compile('[#%]') 
 44  _delim   = re.compile('[\t\,\s]+') 
 45   
 46 -def trigger(line, columns=lsctables.SnglBurst.__slots__): 
 47      """ 
 48      Convert a line from an Omega-format ASCII file into a SnglBurst object. 
 49      """ 
 50   
 51      if isinstance(line, str): 
 52          dat = map(float, _delim.split(line.rstrip())) 
 53      else: 
 54          dat = map(float, line) 
 55   
 56      if len(dat)==8: 
 57          peak, peak_offset, freq, band, duration, n_pix, snr, totPower = dat 
 58          peak = LIGOTimeGPS(peak+peak_offset) 
 59          start    = LIGOTimeGPS(peak-duration/2) 
 60          stop     = LIGOTimeGPS(peak+duration/2) 
 61          av_freq  = freq 
 62          av_band  = band 
 63          err_freq = 0 
 64      else: 
 65          raise ValueError("Wrong number of columns in ASCII line. Cannot read.") 
 66   
 67      # set object 
 68      t = lsctables.SnglBurst() 
 69   
 70      # set times 
 71      if 'start_time' in columns:     t.start_time    = start.gpsSeconds 
 72      if 'start_time_ns' in columns:  t.start_time_ns = start.gpsNanoSeconds 
 73      if 'peak_time' in columns:      t.peak_time     = peak.gpsSeconds 
 74      if 'peak_time_ns' in columns:   t.peak_time_ns  = peak.gpsNanoSeconds 
 75      if 'stop_time' in columns:      t.stop_time     = stop.gpsSeconds 
 76      if 'stop_time_ns' in columns:   t.stop_time_ns  = stop.gpsNanoSeconds 
 77   
 78      # set ms times 
 79      if 'ms_start_time' in columns:     t.ms_start_time    = start.gpsSeconds 
 80      if 'ms_start_time_ns' in columns:  t.ms_start_time_ns = start.gpsNanoSeconds 
 81      if 'ms_stop_time' in columns:      t.ms_stop_time     = stop.gpsSeconds 
 82      if 'ms_stop_time_ns' in columns:   t.ms_stop_time_ns  = stop.gpsNanoSeconds 
 83   
 84      # set other times 
 85      if 'duration' in columns:     t.duration    = duration 
 86      if 'ms_duration' in columns:  t.ms_duration = duration 
 87   
 88      # set frequencies 
 89      if 'central_freq' in columns:         t.central_freq         = freq 
 90      if 'peak_frequency' in columns:       t.peak_frequency       = av_freq 
 91      if 'peak_frequency_eror' in columns:  t.peak_frequency_error = err_freq 
 92      if 'bandwidth' in columns:            t.bandwidth            = av_band 
 93      if 'flow' in columns:                 t.flow                 = freq-band/2 
 94      if 'fhigh' in columns:                t.fhigh                = freq+band/2 
 95   
 96      # set ms frequencies 
 97      if 'ms_bandwidth' in columns: t.ms_bandwidth = band 
 98      if 'ms_flow' in columns:      t.ms_flow      = freq-band/2 
 99      if 'ms_fhigh' in columns:     t.ms_fhigh     = freq+band/2 
100   
101      # set amplitude and snr 
102      if 'snr' in columns:       t.snr       = snr 
103      if 'ms_snr' in columns:    t.ms_snr    = snr 
104      #if 'amplitude' in columns: t.amplitude = totPower**(1/2) 
105   
106      # set other params 
107      if 'numPixels' in columns or 'param_two_value' in columns: 
108          t.param_two_name   = 'numPixels' 
109          t.param_two_value  = n_pix 
110      if 'totPower' in columns or 'param_three_value' in columns: 
111          t.param_three_name = 'cluster_number' 
112          t.param_two_name   = totPower 
113   
114      return t 
115   
116  # ============================================================================= 
117  # Find databases 
118  # ============================================================================= 
119   
120 -def find_databases(host, user="reader", passwd="readonly", match=None): 
121      """ 
122      Query for all databases on the given host 
123      """ 
124      connection = MySQLdb.connect(host=host, user=user, passwd=passwd) 
125      cursor = connection.cursor() 
126      cursor.execute("SHOW DATABASES") 
127      out = cursor.fetchall() 
128      connection.close() 
129      databases = [db[0] for db in out] 
130      if isinstance(match, str): 
131          match = re.compile(match)  
132      if match: 
133          databases = [db for db in databases if match.search(db)] 
134      return databases 
135   
136 -def find_geo_databases(start_time, end_time, host, user="reader",\ 
137                         passwd="readonly"): 
138      """ 
139      Query for all databases named geoYYYYMM whose year and month match the given 
140      [start_time, end_time) interval. 
141      """ 
142      # find all relevant databases 
143      match = "\Ageo\d\d\d\d\d\d\Z" 
144      alldbs = find_databases(host, user=user, passwd=passwd, match=match) 
145   
146      # format dates 
147      start_date = datetime.date(*GPSToUTC(int(start_time))[:3]) 
148      end_date   = datetime.date(*GPSToUTC(int(end_time))[:3]) 
149   
150      # pick out those databases whose dates look correct 
151      databases = [] 
152      d  = start_date 
153      dt = datetime.timedelta(days=28) 
154      while d<=end_date: 
155          db = 'geo%s' % d.strftime("%Y%m") 
156          if db not in databases: 
157              databases.append(db) 
158          d += dt 
159   
160      return databases 
161   
162  # ============================================================================= 
163  # Find databases 
164  # ============================================================================= 
165   
166 -def find_channels(host, database=None, user="reader", passwd="readonly",\ 
167                    match=None): 
168      """ 
169      Query for all channels in the given HACR database on the given host 
170      """ 
171      if not database: 
172          database = find_databases(host, user=user, passwd=passwd)[-1] 
173      connection = MySQLdb.connect(host=host, user=user, passwd=passwd,\ 
174                                   db=database) 
175      cursor = connection.cursor() 
176      query = "select channel from job where monitorName = 'chacr'" 
177      cursor.execute(query) 
178      out = cursor.fetchall() 
179      connection.close() 
180      channels = [ch[0] for ch in out] 
181      if isinstance(match, str):  
182          match = re.compile(match) 
183      if match: 
184          channels = [ch for ch in channels if match.search(ch)] 
185      return channels 
186   
187  # ============================================================================= 
188  # Find databases 
189  # ============================================================================= 
190   
191 -def get_triggers(start_time, end_time, channel, host, columns=None, pid=None,\ 
192                   user="reader", passwd="readonly", verbose=False): 
193      """ 
194      Query the given HACR MySQL host for HACR triggers in a given 
195      [start_time, stop_time) semi-open interval. Returns a LIGOLw SnglBurstTable. 
196      """ 
197      ifo = channel[0:2] 
198      if ifo != "G1": 
199          raise NotImplementedError("Access to non-GEO databases hasn't been "+\ 
200                                    "implemented yet.") 
201   
202      # 
203      # generate table 
204      # 
205   
206      out = lsctables.New(lsctables.SnglBurstTable, columns=columns) 
207      append = out.append 
208   
209      # remove unused names and types 
210      if columns: 
211          # map types 
212          if isinstance(columns[0], str): 
213              columns = map(str.lower, columns) 
214          if isinstance(columns[0], unicode): 
215              columns = map(unicode.lower, columns) 
216          # remove unused 
217          for c in out.columnnames: 
218              if c.lower() not in columns: 
219                  idx = out.columnnames.index(c) 
220                  out.columnnames.pop(idx) 
221                  out.columntypes.pop(idx) 
222      else: 
223          columns = lsctables.SnglBurst.__slots__ 
224       
225      # get the relevant databases 
226      databases  = find_geo_databases(start_time, end_time, host, user=user,\ 
227                                      passwd=passwd) 
228   
229      # 
230      # connect 
231      #  
232   
233      connection = MySQLdb.connect(host=host, user=user, passwd=passwd) 
234      cursor     = connection.cursor() 
235   
236      # 
237      # find triggers 
238      # 
239   
240      for db in databases: 
241          cursor.execute("use %s" % db) 
242   
243          # verify channel in database 
244          query = "select process_id, channel from job where monitorName = "+\ 
245                  "'chacr' and channel = '%s'" % channel 
246          numchan = int(cursor.execute(query)) 
247          result  = cursor.fetchall() 
248          if numchan == 0: 
249              raise AttributeError("%s not found in database %s.\n"\ 
250                                   % (channel, db)) 
251          elif numchan > 1: 
252              sys.stderr.write("warning: Multiple process_ids found for "+\ 
253                               "channel %s in monitor chacr. Triggers from all "\ 
254                               "processes will be amalgamated.\n" % channel) 
255   
256          # get process ids of relevance 
257          process_ids = [int(process[0]) for process in result\ 
258                         if pid==None or int(process[0]) == pid] 
259   
260          # loop over process ids 
261          for pid in process_ids: 
262              query = "select gps_start, gps_offset, freq_central, bandwidth, "\ 
263                      "duration, num_pixels, snr, totPower from mhacr where "\ 
264                      "process_id = %s AND gps_start >= %s AND gps_start < %s "\ 
265                      "order by gps_start asc" % (pid, start_time, end_time) 
266              ntrigs = cursor.execute(query) 
267              result = cursor.fetchall() 
268              if ntrigs: 
269                  for row in result: 
270                      append(trigger(row, columns=columns)) 
271   
272      connection.close() 
273      return out 
274