Pygrb offline workflow

bin/pycbc_multi_inspiral

@@ -70,7 +70,8 @@ parser.add_argument("--cluster-method", choices=["template", "window"],
                          "chirp length.")
 parser.add_argument("--cluster-window", type=float, default=0,
                     help="Length of clustering window in seconds.")
-parser.add_argument("--bank-veto-bank-file", type=str)
+parser.add_argument("--bank-veto-bank-file", type=str, help="Path to the "
+                    "bank file used to compute the the bank chi-square veto.")
 parser.add_argument("--chisq-bins", default=0)
 # Commenting out options which are not yet implemented
 # parser.add_argument("--chisq-threshold", type=float, default=0) 

bin/pygrb/pycbc_pygrb_efficiency

@@ -98,24 +98,20 @@ parser.add_argument("-g", "--glitch-check-factor", action="store",
                     "exclusion efficiencies this value is multiplied " +
                     "to the offsource around the injection trigger to " +
                     "determine if it is just a loud glitch.")
+parser.add_argument("--injection-set-name", action="store", type=str,
+                    default="", help="Name of the injection set to be " +
+                    "used in the plot title.")
 parser.add_argument("-C", "--cluster-window", action="store", type=float,
                     default=0.1, help="The cluster window used " +
                     "to cluster triggers in time.")
-ppu.pygrb_add_missed_injs_input_opt(parser)
 ppu.pygrb_add_injmc_opts(parser)
 ppu.pygrb_add_bestnr_opts(parser)
 opts = parser.parse_args()
 
 init_logging(opts.verbose, format="%(asctime)s: %(levelname)s: %(message)s")
 
 # Check options
-if (opts.found_file is None) and (opts.missed_file is None):
-    do_injections = False
-elif (opts.found_file) and opts.missed_file:
-    do_injections = True
-else:
-    err_msg = "Must provide both found and missed file if running injections."
-    parser.error(err_msg)
+do_injections = opts.found_missed_file is not None
 
 if not opts.newsnr_threshold:
     opts.newsnr_threshold = opts.snr_threshold
@@ -124,9 +120,8 @@ if not opts.newsnr_threshold:
 outdir = os.path.split(os.path.abspath(opts.background_output_file))[0]
 trig_file = opts.offsource_file
 onsource_file = opts.onsource_file
-found_file = opts.found_file
-missed_file = opts.missed_file
-inj_set_name = os.path.split(os.path.abspath(missed_file))[1].split('INSPINJ')[1].split('_')[1]
+found_missed_file = opts.found_missed_file
+inj_set_name = opts.injection_set_name
 chisq_index = opts.chisq_index
 chisq_nhigh = opts.chisq_nhigh
 wf_err = opts.waveform_error
@@ -376,7 +371,7 @@ if do_injections:
     logging.info("%d found injections analysed.", len(found_injs))
 
     # Missed injections (ones not recovered at all)
-    missed_injs = ppu.load_injections(missed_file, vetoes, sim_table=True)
+    missed_injs = ppu.load_injections(found_missed_file, vetoes, sim_table=True)
 
     # Process missed injections 'missed_inj'
     missed_inj = {}

bin/pygrb/pycbc_pygrb_minifollowups

@@ -35,8 +35,9 @@ import pycbc.workflow.minifollowups as mini
 import pycbc.version
 import pycbc.events
 from pycbc.results import layout
-from pycbc.results import pygrb_postprocessing_utils as ppu
+from pycbc.results.pygrb_postprocessing_utils import extract_ifos
 from pycbc.workflow.plotting import PlotExecutable
+from pycbc.workflow.grb_utils import build_veto_filelist, build_segment_filelist
 
 __author__ = "Francesco Pannarale <[email protected]>"
 __version__ = pycbc.version.git_verbose_msg
@@ -72,10 +73,10 @@ def make_timeseries_plot(workflow, trig_file, snr_type, central_time,
                           tags=tags+extra_tags).create_node()
     node.add_input_opt('--trig-file', trig_file)
     # Pass the veto files
-    veto_files = ppu.build_veto_filelist(workflow)
+    veto_files = build_veto_filelist(workflow)
     node.add_input_list_opt('--veto-files', veto_files)
     # Pass the segment files
-    seg_files = ppu.build_segment_filelist(workflow)
+    seg_files = build_segment_filelist(workflow)
     node.add_input_list_opt('--seg-files', seg_files)
     # Other shared tuning values
     for opt in ['chisq-index', 'chisq-nhigh', 'null-snr-threshold',
@@ -155,7 +156,7 @@ num_events = min(num_events, len(fp['BestNR'][:]))
 
 # Determine ifos used in the analysis
 trig_file = resolve_url_to_file(os.path.abspath(args.trig_file))
-ifos = ppu.extract_ifos(os.path.abspath(args.trig_file))
+ifos = extract_ifos(os.path.abspath(args.trig_file))
 num_ifos = len(ifos)
 
 # (Loudest) off/on-source events are on time-slid data so the

bin/pygrb/pycbc_pygrb_page_tables

@@ -48,6 +48,38 @@ __program__ = "pycbc_pygrb_page_tables"
 # =============================================================================
 # Functions
 # =============================================================================
+# Function to load trigger data
+def load_data(input_file, vetoes, ifos, opts, injections=False):
+    """Load data from a trigger/injection file"""
+
+    # Initialize the dictionary
+    data = {}
+    data['trig_time'] = None
+    data['coherent'] = None
+    data['reweighted_snr'] = None
+
+    if input_file:
+        if injections:
+            logging.info("Loading injections...")
+            # This will eventually become ppu.load_injections()
+            trigs_or_injs = ppu.load_triggers(input_file, vetoes)
+        else:
+            logging.info("Loading triggers...")
+            trigs_or_injs = ppu.load_triggers(input_file, vetoes)
+
+        data['trig_time'] = trigs_or_injs['network/end_time_gc'][:]
+        num_trigs_or_injs = len(data['trig_time'])
+        data['coherent'] = trigs_or_injs['network/coherent_snr'][:]
+        data['reweighted_snr'] = trigs_or_injs['network/reweighted_snr'][:]
+
+        if injections:
+            logging.info("%d injections found.", num_trigs_or_injs)
+        else:
+            logging.info("%d triggers found.", num_trigs_or_injs)
+
+    return data
+
+
 def get_column(table, column_name):
     """Wrapper for glue.ligolw.lsctables.MultiInspiralTable.get_column
     method. Easier for h5py switch. Note: Must still replace
@@ -190,7 +222,6 @@ parser.add_argument("-F", "--offsource-file", action="store", required=True,
 # As opposed to offsource-file and trig-file, this only contains onsource
 parser.add_argument("--onsource-file", action="store", default=None,
                     help="Location of on-source trigger file.")
-ppu.pygrb_add_missed_injs_input_opt(parser)
 ppu.pygrb_add_injmc_opts(parser)
 ppu.pygrb_add_bestnr_opts(parser)
 parser.add_argument("--num-loudest-off-trigs", action="store",
@@ -200,6 +231,7 @@ parser.add_argument("--quiet-found-injs-output-file", default=None, #required=Tr
                     help="Quiet-found injections html output file.")
 parser.add_argument("--missed-found-injs-output-file", default=None, #required=True,
                     help="Missed-found injections html output file.")
+# TODO: group hdf5 files into a single one
 parser.add_argument("--quiet-found-injs-h5-output-file", default=None, #required=True,
                     help="Quiet-found injections h5 output file.")
 parser.add_argument("--loudest-offsource-trigs-output-file", default=None, #required=True,
@@ -241,9 +273,9 @@ if output_files.count(None) == len(output_files):
 if opts.quiet_found_injs_output_file or opts.missed_found_injs_output_file or\
     opts.quiet_found_injs_h5_output_file:
     do_injections = True
-    if (opts.found_file is None) and (opts.missed_file is None):
-        err_msg = "Must provide both found and missed injections file "
-        err_msg += "locations if processing injections."
+    if opts.found_missed_file is None:
+        err_msg = "Must provide found-missed injections file "
+        err_msg += "location if processing injections."
         parser.error(err_msg)
 else:
     do_injections = False
@@ -260,8 +292,7 @@ if not opts.newsnr_threshold:
 # Store options used multiple times in local variables
 trig_file = opts.offsource_file
 onsource_file = opts.onsource_file
-found_file = opts.found_file
-missed_file = opts.missed_file
+found_missed_file = opts.found_missed_file
 chisq_index = opts.chisq_index
 chisq_nhigh = opts.chisq_nhigh
 wf_err = opts.waveform_error
@@ -316,8 +347,7 @@ ifos, vetoes = ppu.extract_ifos_and_vetoes(trig_file, opts.veto_files,
 
 # Load triggers, time-slides, and segment dictionary
 logging.info("Loading triggers.")
-trigs = ppu.load_xml_table(trig_file, lsctables.MultiInspiralTable.tableName)
-logging.info("%d triggers loaded.", len(trigs))
+trig_data = load_data(trig_file, vetoes, ifos, opts)
 logging.info("Loading timeslides.")
 slide_dict = ppu.load_time_slides(trig_file)
 logging.info("Loading segments.")
@@ -587,8 +617,8 @@ else:
 # =======================
 if do_injections:
 
-    found_trig_table = ppu.load_injections(found_file, vetoes)
-    found_inj_table = ppu.load_injections(found_file, vetoes, sim_table=True)
+    found_trig_table = ppu.load_injections(found_missed_file, vetoes)
+    found_inj_table = ppu.load_injections(found_missed_file, vetoes, sim_table=True)
 
     logging.info("Missed/found injections/triggers loaded.")
 
@@ -618,7 +648,7 @@ if do_injections:
     logging.info("%d found injections analysed.", len(found_injs['mchirp']))
 
     # Missed injections (ones not recovered at all)
-    missed_inj_table = ppu.load_injections(missed_file, vetoes, sim_table=True)
+    missed_inj_table = ppu.load_injections(found_missed_file, vetoes, sim_table=True)
     missed_injs = lsctable_to_dict(missed_inj_table, ifos, opts)
 
     # Avoids a problem with formatting in the non-static html output file

bin/pygrb/pycbc_pygrb_plot_chisq_veto

@@ -76,8 +76,12 @@ def load_data(input_file, vetoes, ifos, opts, injections=False):
         if snr_type in ['coherent', 'null', 'reweighted']:
             data[snr_type] = trigs_or_injs['network/%s_snr' % snr_type][:]
         elif snr_type == 'single':
-            att = opts.ifo[0].lower()
-            data[snr_type] = trigs_or_injs['%s/snr_%s' % (opts.ifo, att)][:]
+            key = opts.ifo + '/snr_' + opts.ifo.lower()
+            # TODO: Undoes current L710 in EventManagerCoherent.write_to_hdf
+            # and is ugly hardcoded!
+            if opts.ifo.lower() != 'h1':
+                key = key[:-1]
+            data[snr_type] = trigs_or_injs[key][:]
 
         # Calculate coincident SNR
         elif snr_type == 'coincident':
@@ -88,20 +92,38 @@ def load_data(input_file, vetoes, ifos, opts, injections=False):
                      'bank': 'bank_chisq',
                      'auto': 'cont_chisq'}
 
-        data[veto_type] = \
-            trigs_or_injs['%s/%s' % (opts.ifo, veto_tags[veto_type])][:]
+        # FIXME: network data contains my_network_chisq
+        chisq_key = 'network/my_network_chisq' if injections \
+            else opts.ifo + '/' + veto_tags[veto_type]
+        # Will not work until single IFO data is in found_missed_file
+        data[veto_type] = trigs_or_injs[chisq_key][:]
 
         # Floor single IFO chi-square at 0.005
         numpy.putmask(data[veto_type], data[veto_type] == 0, 0.005)
 
-        data['dof'] = numpy.unique(
-            trigs_or_injs['%s/%s_dof' % (opts.ifo, veto_tags[veto_type])][:])
+        if not injections:
+            dof_key = '%s/%s_dof' % (opts.ifo, veto_tags[veto_type])
+            data['dof'] = numpy.unique(
+                trigs_or_injs[dof_key][:])
 
-        logging.info("%d triggers found.", num_trigs_or_injs)
+        label = "injections" if injections else "triggers"
+
+        logging.info("{0} {1} found.".format(num_trigs_or_injs, label))
 
     return data
 
 
+# Function to calculate chi-square weight for the reweighted SNR
+def new_snr_chisq(snr, new_snr, chisq_dof, chisq_index=4.0, chisq_nhigh=3.0):
+    """Returns the chi-square value needed to weight SNR into new SNR"""
+
+    chisqnorm = (snr/new_snr)**chisq_index
+    if chisqnorm <= 1:
+        return 1E-20
+
+    return chisq_dof * (2*chisqnorm - 1)**(chisq_nhigh/chisq_index)
+
+
 # Function that produces the contrours to be plotted
 def calculate_contours(trig_data, opts, veto_type, new_snrs=None):
     """Generate the contours for the veto plots"""
@@ -130,9 +152,9 @@ def calculate_contours(trig_data, opts, veto_type, new_snrs=None):
     # Loop over SNR values and calculate chisq variable needed
     for j, snr in enumerate(snr_vals):
         for i, new_snr in enumerate(new_snrs):
-            contours[i][j] = plu.new_snr_chisq(snr, new_snr, dof,
-                                               opts.chisq_index,
-                                               opts.chisq_nhigh)
+            contours[i][j] = new_snr_chisq(snr, new_snr, dof,
+                                           opts.chisq_index,
+                                           opts.chisq_nhigh)
 
     return contours, snr_vals, cont_value
 
@@ -161,14 +183,17 @@ parser = ppu.pygrb_initialize_plot_parser(description=__doc__,
 parser.add_argument("-t", "--trig-file", action="store",
                     default=None, required=True,
                     help="The location of the trigger file")
+parser.add_argument("--found-missed-file",
+                    help="The hdf injection results file", required=False)
 parser.add_argument("-z", "--zoom-in", default=False, action="store_true",
                     help="Output file a zoomed in version of the plot.")
 parser.add_argument("-y", "--y-variable", required=True,
                     choices=['standard', 'bank', 'auto'],
                     help="Quantity to plot on the vertical axis.")
-parser.add_argument("--snr-type", default=None,
-                    choices=['coherent', 'single', 'coincident', 'null',
-                    'reweighted'], help="SNR value to plot on x-axis. Can be " +
+# TODO: get 'single' to work
+parser.add_argument("--snr-type", default='coherent',
+                    choices=['coherent', 'coincident', 'null', 'reweighted',
+                    'single'], help="SNR value to plot on x-axis. Can be " +
                     "coherent, coincident, null, reweighted, or single IFO SNR")
 ppu.pygrb_add_bestnr_opts(parser)
 opts = parser.parse_args()
@@ -177,7 +202,7 @@ init_logging(opts.verbose, format="%(asctime)s: %(levelname)s: %(message)s")
 
 # Check options
 trig_file = os.path.abspath(opts.trig_file)
-found_file = os.path.abspath(opts.found_file) if opts.found_file else None
+found_missed_file = os.path.abspath(opts.found_missed_file) if opts.found_missed_file else None
 zoom_in = opts.zoom_in
 veto_type = opts.y_variable
 ifo = opts.ifo
@@ -186,13 +211,14 @@ snr_type = opts.snr_type
 # otherwise the single IFO SNR is used
 if snr_type == 'single':
     if ifo is None:
-        err_msg = "--ifo must be given when using --use-sngl-ifo-snr"
+        err_msg = "--ifo must be given to plot single IFO SNR veto"
         parser.error(err_msg)
 
 # Veto is intended as a single IFO quantity. Network chisq will be obsolete.
-if ifo is None:
-    err_msg = "--ifo must be given to plot single IFO SNR veto"
-    parser.error(err_msg)
+# TODO: fix vetoes
+#if ifo is None:
+#    err_msg = "--ifo must be given to plot single IFO SNR veto"
+#    parser.error(err_msg)
 
 # Prepare plot title and caption
 veto_labels = {'standard': "Chi Square Veto",
@@ -211,7 +237,7 @@ if opts.plot_caption is None:
                          "Black line: veto line.  " +
                          "Gray shaded region: Vetoed area.  " +
                          "Yellow lines: contours of new SNR.")
-    if found_file:
+    if found_missed_file:
         opts.plot_caption = ("Red crosses: injections triggers.  ") +\
                             opts.plot_caption
 
@@ -235,7 +261,7 @@ if ifo and ifo not in ifos:
 trig_data = load_data(trig_file, vetoes, ifos, opts)
 
 # Extract (or initialize) injection data
-inj_data = load_data(found_file, vetoes, ifos, opts, injections=True)
+inj_data = load_data(found_missed_file, vetoes, ifos, opts, injections=True)
 
 # Sanity checks
 if trig_data[snr_type] is None and inj_data[snr_type] is None:
@@ -259,14 +285,14 @@ if snr_type == 'coincident':
 
 # Determine the minumum and maximum SNR value we are dealing with
 x_min = opts.sngl_snr_threshold
-x_max = 1.1*plu.axis_max_value(trig_data[snr_type], inj_data[snr_type], found_file)
+x_max = 1.1*plu.axis_max_value(trig_data[snr_type], inj_data[snr_type], found_missed_file)
 
 # Determine y-axis minimum value and label
 y_label = "%s Single %s" % (ifo, veto_labels[veto_type])
 y_min = 1
 
 # Determine the maximum bank veto value we are dealing with
-y_max = plu.axis_max_value(trig_data[veto_type], inj_data[veto_type], found_file)
+y_max = plu.axis_max_value(trig_data[veto_type], inj_data[veto_type], found_missed_file)
 
 # Determine contours for plots
 conts = None