Merge branch '103-add-options-to-bayeswave_pipe-to-include-glitch-cleaning' into 'master'

Add glitch cleaning functionality for multiple detectors

See merge request !193

BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py

@@ -996,8 +996,8 @@ def condor_job_config(job_type, condor_job, config_parser):
     #
     if config_parser.getboolean('condor', 'osg-deploy'):
         # --- Force downstream jobs to run locally
-        if job_type in ['bayeswave_post', 'bayeswave_fpeak', 'megaplot.py',
-                        'megasky.py']:
+        if job_type in ['bayeswave_post', 'bayeswave_fpeak', 'bayeswave_clean_frame',
+                         'megaplot.py', 'megasky.py']:
             requires.append("(GLIDEIN_SITE=?=undefined)")
         else:
             try:
@@ -1605,12 +1605,14 @@ class bayeswave_clean_frameJob(pipeline.CondorDAGJob, pipeline.AnalysisJob):
         self.add_condor_cmd('initialdir', workdir)
         self.set_sub_file(os.path.join(workdir, 'bayeswave_clean_frame.sub'))
 
+
         self.set_stdout_file(os.path.join(workdir, 'logs',
-            'BayesWaveCleanFrame_$(macrooutputDir)-$(cluster)-$(process).out'))
+            'BayesWaveCleanFrame_$(macrooutdir)-$(macroifo)-$(cluster)-$(process).out'))
         self.set_stderr_file(os.path.join(workdir, 'logs',
-            'BayesWaveCleanFrame_$(macrooutputDir)-$(cluster)-$(process).err'))
+            'BayesWaveCleanFrame_$(macrooutdir)-$(macroifo)-$(cluster)-$(process).err'))
         self.set_log_file(os.path.join(workdir, 'logs',
-            'BayesWaveCleanFrame_$(macrooutputDir)-$(cluster)-$(process).log'))
+            'BayesWaveCleanFrame_$(macrooutdir)-$(macroifo)-$(cluster)-$(process).log'))
+
 
         if cp.has_option('condor','arch'):
             self.add_condor_cmd('+arch',cp.get('condor','arch'))
@@ -1633,28 +1635,16 @@ class bayeswave_clean_frameJob(pipeline.CondorDAGJob, pipeline.AnalysisJob):
         if cp.getboolean('condor', 'transfer-files'):
 
             # File transfers
-            transferstring='$(macrooutputDir),datafind'
+            transferstring='$(macrooutdir),datafind'
 
             if cp.getboolean('condor','copy-frames'):
                 transferstring += ',$(macroframes)'
 
             self.add_condor_cmd('transfer_input_files', transferstring)
-            self.add_condor_cmd('transfer_output_files', '$(macrooutputDir)')
+            self.add_condor_cmd('transfer_output_files', '$(macrooutdir)')
 
-
-        ifo_clean = ast.literal_eval(cp.get('bayeswave_clean_frame_options', 'ifo'))
-        self.add_opt('ifo', ifo_clean)
-   
-        channel_list = ast.literal_eval(cp.get('datafind', 'channel-list'))
-        frame_list = ast.literal_eval(cp.get('datafind', 'frtype-list'))
-
-        self.add_opt('channel', channel_list[ifo_clean])
-        self.add_opt('frame-type', frame_list[ifo_clean])
-         
-        
-        for item in cp.items('bayeswave_clean_frame_options'):
-            # Add any option and value which exists
-            self.add_opt(item[0], item[1])
+        clean_suffix = cp.get('bayeswave_clean_frame_options', 'clean-suffix')
+        self.add_opt('clean-suffix', clean_suffix)
 
 class bayeswave_clean_frameNode(pipeline.CondorDAGNode, pipeline.AnalysisNode):
 
@@ -1691,7 +1681,19 @@ class bayeswave_clean_frameNode(pipeline.CondorDAGNode, pipeline.AnalysisNode):
     
     def set_outdir(self, outdir):
         self.add_var_opt('outdir', outdir)
-  
+    
+    def set_ifo(self, ifo):
+        self.add_var_opt('ifo', ifo)
+   
+    def set_channel_name(self, channel):
+        self.add_var_opt('channel', channel)
+
+    def set_frame_type(self, frame_type):
+        self.add_var_opt('frame-type', frame_type)
+
+    def set_cache_file(self, cache_file):
+         self.add_var_opt('cachefile', cache_file)
+
 
 
 

BayesWaveUtils/scripts/bayeswave_pipe

@@ -41,6 +41,7 @@ from glue import pipeline
 
 from lalapps import inspiralutils
 from ligo import segments
+import lalframe.frread as fr
 
 from bayeswave_pipe import bayeswave_pipe_utils as pipe_utils
 
@@ -271,10 +272,12 @@ if cp.has_option('bayeswave_fpeak_options', 'fpeak-analysis'):
     if cp.has_option('bayeswave_fpeak_options','fpeak-flow'):
         fpeak_flow=cp.getfloat('bayeswave_fpeak_options','fpeak-flow')
 
-if cp.has_option('bayeswave_clean_frame_options', 'ifo'):
+if cp.has_option('bayeswave_clean_frame_options', 'ifos'):
     opts.bayeswave_clean_frame=True
 
-    clean_frame_ifo=cp.get('bayeswave_clean_frame_options', 'ifo')
+    clean_frame_ifos=ast.literal_eval(cp.get('bayeswave_clean_frame_options', 'ifos'))
+    clean_frame_srates = {}
+    channel_names = ast.literal_eval(cp.get('datafind','channel-list')) 
    
 
 if opts.bayeswave_clean_frame:
@@ -715,18 +718,14 @@ for ifo in ifo_list:
             os.chdir(curdiir)
         
         # Get the sampling rate of the frame
-        if opts.bayeswave_clean_frame:
+        if opts.bayeswave_clean_frame and ifo in clean_frame_ifos:
             try:
-                import lalframe.frread as fr
-                clean_ifo = ast.literal_eval(cp.get('bayeswave_clean_frame_options', 'ifo'))
                 trigtime = int(trigger_list.triggers[0].trigger_time)
-                channel_names = ast.literal_eval(cp.get('datafind','channel-list'))
-                example_ts = fr.read_timeseries(cache_files[clean_ifo], channel_names[clean_ifo], trigtime-1,2.0 )
-                clean_frame_srate = 1./example_ts.deltaT
-                print("Getting sampling rate of cleaned frame from input frame: %.2f Hz"%(clean_frame_srate))
+                example_ts = fr.read_timeseries(cache_files[ifo], channel_names[ifo], trigtime-1,2.0 )
+                clean_frame_srates[ifo] = 1./example_ts.deltaT
             except Exception as e:
-                clean_frame_srate = 16384.0
-                print("Setting clean frame srate to default: %.2f"%(clean_frame_srate))
+                clean_frame_srates[ifo] = 16384.0
+                print("Setting %s clean frame srate to default: %.2f"%(ifo, 16384.0))
             
         
 
@@ -1187,24 +1186,42 @@ for t,trigger in enumerate(trigger_list.triggers):
         # Add options for bayeswave_clean_frame node
         #
         if opts.bayeswave_clean_frame:
-            clean_frame_outdir = os.path.join(outputDir, 'clean_frame')
-            if not os.path.exists(clean_frame_outdir): 
-                os.makedirs(clean_frame_outdir)
-            bayeswave_clean_frame_node =  pipe_utils.bayeswave_clean_frameNode(bayeswave_clean_frame_job)
-            bayeswave_clean_frame_node.set_trigtime(trigger.trigger_time)
-            clean_frame_duration = 4.0*trigger.seglen
-            bayeswave_clean_frame_node.set_frame_start(int(trigger.trigger_time - clean_frame_duration/2.0))
-            bayeswave_clean_frame_node.set_frame_srate(clean_frame_srate)
-            bayeswave_clean_frame_node.set_seglen(trigger.seglen)
-            
+            bayeswave_clean_frame_nodes = []
             if cp.has_option('bayeswave_options', 'fullOnly'):
                 model_type = 'full'
             else:
                 model_type = 'glitch'
-            clean_ifo = ast.literal_eval(cp.get('bayeswave_clean_frame_options', 'ifo'))
-            glitch_param_file = os.path.join(outputDir, 'chains', '%s_params_%s.dat'%(model_type, clean_ifo))
-            bayeswave_clean_frame_node.set_glitch_param_file(glitch_param_file)
-            bayeswave_clean_frame_node.set_outdir(clean_frame_outdir)
+            clean_frame_ifos = ast.literal_eval(cp.get('bayeswave_clean_frame_options', 'ifos'))
+            for  ifo in clean_frame_ifos:
+                clean_frame_node =  pipe_utils.bayeswave_clean_frameNode(bayeswave_clean_frame_job)
+                clean_frame_node.set_trigtime(trigger.trigger_time)
+                clean_frame_node.set_ifo(ifo)
+                clean_frame_node.set_segment_start(segment_start)
+                 
+                trigtime = int(trigger.trigger_time)
+                example_ts = fr.read_timeseries(cache_files[ifo], channel_names[ifo], trigtime-1,2.0 )
+                clean_frame_srate = 1./example_ts.deltaT
+
+                for seg in frameSegs[ifo]:
+                    if(trigger.trigger_time in seg):
+                        clean_frame_start_time = seg[0]
+                        clean_frame_length = seg[1] - seg[0]
+                
+                
+                clean_frame_node.set_frame_start(clean_frame_start_time)
+                clean_frame_node.set_frame_length(clean_frame_length)
+                clean_frame_node.set_cache_file(cache_files[ifo])
+                clean_frame_node.set_frame_srate(clean_frame_srate)
+                clean_frame_node.set_seglen(trigger.seglen)            
+            
+                glitch_param_file = os.path.join(outputDir, 'chains', '%s_params_%s.dat.0'%(model_type, ifo))
+                clean_frame_node.set_glitch_param_file(glitch_param_file)
+                clean_frame_node.set_outdir(outputDir)
+
+                clean_frame_node.set_channel_name(channel_names[ifo])
+                clean_frame_node.set_frame_type(frtype_list[ifo])
+                
+                bayeswave_clean_frame_nodes.append(clean_frame_node)
             
         
 
@@ -1250,7 +1267,8 @@ for t,trigger in enumerate(trigger_list.triggers):
             if opts.fpeak_analysis:
                 bayeswave_fpeak_node.add_parent(bayeswave_node)
             if opts.bayeswave_clean_frame:
-                bayeswave_clean_frame_node.add_parent(bayeswave_node)
+                for clean_frame_node in bayeswave_clean_frame_nodes:
+                    clean_frame_node.add_parent(bayeswave_node)
         if not opts.skip_megapy:
             megasky_node.add_parent(bayeswave_node)
             megaplot_node.add_parent(bayeswave_post_node) 
@@ -1265,7 +1283,8 @@ for t,trigger in enumerate(trigger_list.triggers):
 
         dag.add_node(bayeswave_node)
         if opts.bayeswave_clean_frame:
-            dag.add_node(bayeswave_clean_frame_node)
+            for clean_frame_node in bayeswave_clean_frame_nodes:
+                dag.add_node(clean_frame_node)
 
         if not opts.skip_post and not opts.separate_post_dag:
             dag.add_node(bayeswave_post_node)

src/BayesWaveCleanFrame.c

 /*
- *  Copyright (C) 2018 Neil J. Cornish, Tyson B. Littenberg, James A. Clark, Jonah B. Kanner
+ *  Copyright (C) 2018 Neil J. Cornish, Tyson B. Littenberg, James A. Clark, Jonah B. Kanner,
+ *  Sudarshan Ghonge
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -25,6 +26,7 @@
 #include <math.h>
 #include <string.h>
 #include <fftw3.h>
+#include <sys/stat.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
@@ -91,7 +93,7 @@ struct CleanFrameData
 void print_usage();
 void parse_command_line_args(int argc, char **argv, struct CleanFrameData *data);
 
-static void output_frame(REAL8TimeSeries *timeData, REAL8TimeSeries *timeRes, REAL8TimeSeries *timeGlitch, CHAR *frameType, CHAR *ifo, CHAR *outdir);
+static void output_frame(REAL8TimeSeries *timeData, REAL8TimeSeries *timeRes, REAL8TimeSeries *timeGlitch, CHAR *frameType, CHAR *ifo, CHAR *cleandir);
 
 void printProgress (double percentage)
 {
@@ -181,6 +183,10 @@ int main(int argc, char *argv[])
   version = data->clean_suffix;
   outdir = data->outdir;
 
+  char cleandir[MAXSTRINGSIZE];
+  sprintf(cleandir,"%s/clean_frame", data->outdir);
+  mkdir(cleandir,S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+
   char outframeType[MAXSTRINGSIZE];
   char outframeChannel[MAXSTRINGSIZE];
   char outframeGlitchChannel[MAXSTRINGSIZE];
@@ -349,7 +355,7 @@ int main(int argc, char *argv[])
   
 
   // Output cleaned data!
-  output_frame(timeData, timeRes, timeGlitch, outframeType, data->ifo, outdir);
+  output_frame(timeData, timeRes, timeGlitch, outframeType, data->ifo, cleandir);
   XLALDestroyREAL8TimeSeries(timeData);
 
   fprintf(stdout,"\n");
@@ -368,7 +374,7 @@ static void output_frame(REAL8TimeSeries *timeData,
                          REAL8TimeSeries *timeGlitch,
                          CHAR *frameType,
                          CHAR *ifo,
-                         CHAR *outdir)
+                         CHAR *cleandir)
 {
   CHAR fname[2048];
   INT4 duration;
@@ -400,7 +406,7 @@ static void output_frame(REAL8TimeSeries *timeData,
 
   /* get frame filename */
   duration = gpsEnd - gpsStart;
-  snprintf( fname, FILENAME_MAX, "%s/%c-%s-%d-%d.gwf", outdir, ifo[0], frameType, gpsStart, duration );
+  snprintf( fname, FILENAME_MAX, "%s/%c-%s-%d-%d.gwf", cleandir, ifo[0], frameType, gpsStart, duration );
 
   /* define frame */
   frame = XLALFrameNew( &timeData->epoch, duration, "LIGO", 0, 1, detectorFlags );