Merge branch 'gwastro:master' into offline_gracedb_prep_scripts

bin/all_sky_search/pycbc_add_statmap

@@ -361,8 +361,14 @@ fg_fars_out = np.sum(isincombo_mask * fg_fars, axis=0)
 fg_fars_exc_out = np.sum(isincombo_mask * fg_fars_exc, axis=0)
 
 # Apply any limits as appropriate
-fg_fars_out = significance.apply_far_limit(fg_fars_out, significance_dict, combo=fg_coinc_type)
-fg_fars_exc_out = significance.apply_far_limit(fg_fars_exc_out, significance_dict, combo=fg_coinc_type)
+fg_fars_out = significance.apply_far_limit(
+    fg_fars_out,
+    significance_dict,
+    combo=fg_coinc_type)
+fg_fars_exc_out = significance.apply_far_limit(
+    fg_fars_exc_out,
+    significance_dict,
+    combo=fg_coinc_type)
 
 fg_ifar = conv.sec_to_year(1. / fg_fars_out)
 fg_ifar_exc = conv.sec_to_year(1. / fg_fars_exc_out)
@@ -562,6 +568,7 @@ while True:
     final_combined_fg = final_combined_fg + \
                             combined_fg_data.select(where_combined)
     combined_fg_data = combined_fg_data.remove(where_combined)
+    fg_coinc_type = np.delete(fg_coinc_type, where_combined)
     n_triggers -= 1
 
     logging.info('Removing background triggers at time {} within window '
@@ -605,6 +612,17 @@ while True:
             sep_bg_data[key].data['decimation_factor'],
             bg_t_y,
             **significance_dict[ifo_combo_key])
+        fg_far = significance.apply_far_limit(
+            fg_far,
+            significance_dict,
+            combo=key,
+        )
+        bg_far = significance.apply_far_limit(
+            bg_far,
+            significance_dict,
+            combo=key,
+        )
+
         sep_bg_data[key].data['ifar'] = 1. / bg_far
         sep_fg_data[key].data['ifar'] = 1. / fg_far
         sep_fg_data[key].data['fap'] = 1 - \
@@ -631,9 +649,15 @@ while True:
     isincombo_mask = np.array([list(is_in_combo_time[ct])
                                for ct in all_ifo_combos])
     fg_fars = np.array([list(far[ct]) for ct in all_ifo_combos])
+    fg_fars_out = np.sum(isincombo_mask * fg_fars, axis=0)
+    fg_fars_out = significance.apply_far_limit(
+        fg_fars_out,
+        significance_dict,
+        combo=fg_coinc_type,
+    )
     # Combine the FARs with the mask to obtain the new ifars
     combined_fg_data.data['ifar'] = conv.sec_to_year(
-        1. / np.sum(isincombo_mask * fg_fars, axis=0))
+        1. / fg_fars_out)
     fg_time -= args.cluster_window
     combined_fg_data.data['fap'] = 1 - \
         np.exp(-conv.sec_to_year(fg_time) / combined_fg_data.data['ifar'])

bin/all_sky_search/pycbc_calculate_dqflag

@@ -26,6 +26,8 @@ parser.add_argument('--flag', type=str, required=True,
                     help="name of the data quality flag")
 parser.add_argument('--ifo', type=str, required=True,
                     help="interferometer to analyze")
+parser.add_argument('--sample-frequency', required=True, type=int,
+                    help="timeseries sampling frequency in Hz")
 parser.add_argument("--output-file", required=True,
                     help="name of output file")
 parser.add_argument("--output-channel", required=False,
@@ -35,13 +37,14 @@ args = parser.parse_args()
 pycbc.init_logging(args.verbose)
 
 
-def make_dq_ts(dq_segs, start, end, ifo, flag):
+def make_dq_ts(dq_segs, start, end, ifo, flag, freq):
     """ Create a data quality timeseries
     """
     logging.info('Creating data quality timeseries for flag %s', 
                  flag)
     # Keep just times which belong to science_segments
-    dq_times = numpy.arange(start, end)
+    dur = end - start
+    dq_times = numpy.linspace(start, end, int(freq*dur), endpoint=False)
     dq = numpy.zeros(len(dq_times))
     # Identify times within segments for the chosen flag
     indexes, _ = veto.indices_within_segments(dq_times, [dq_segs], ifo=ifo,
@@ -51,7 +54,7 @@ def make_dq_ts(dq_segs, start, end, ifo, flag):
     else:
         logging.warning('Veto definer segment list is empty for flag %s-%s',
                         ifo, flag)
-    return TimeSeries(dq, epoch=start, delta_t=1.)
+    return TimeSeries(dq, epoch=start, delta_t=1.0/freq)
 
 
 ifo = args.ifo
@@ -63,7 +66,7 @@ if len(flag_list)>1:
 
 # Create data quality time series
 dq = make_dq_ts(args.dq_segments, args.gps_start_time, args.gps_end_time, 
-                ifo, flag)
+                ifo, flag, args.sample_frequency)
 
 dq.save(args.output_file, group=args.output_channel)
 

bin/all_sky_search/pycbc_coinc_findtrigs

@@ -181,7 +181,7 @@ if args.randomize_template_order:
     shuffle(template_ids)
     template_ids = template_ids[tmin:tmax]
 else:
-    template_ids = np.array([range(tmin, tmax)])
+    template_ids = numpy.array([range(tmin, tmax)])
 
 original_bank_len = len(template_ids)
 

bin/all_sky_search/pycbc_coinc_statmap

@@ -257,10 +257,22 @@ bg_far_exc, fg_far_exc = significance.get_far(
     background_time_exc,
     **significance_dict[ifo_combo])
 
-fg_far = significance.apply_far_limit(fg_far, significance_dict, combo=ifo_combo)
-bg_far = significance.apply_far_limit(bg_far, significance_dict, combo=ifo_combo)
-fg_far_exc = significance.apply_far_limit(fg_far_exc, significance_dict, combo=ifo_combo)
-bg_far_exc = significance.apply_far_limit(bg_far_exc, significance_dict, combo=ifo_combo)
+fg_far = significance.apply_far_limit(
+    fg_far,
+    significance_dict,
+    combo=ifo_combo)
+bg_far = significance.apply_far_limit(
+    bg_far,
+    significance_dict,
+    combo=ifo_combo)
+fg_far_exc = significance.apply_far_limit(
+    fg_far_exc,
+    significance_dict,
+    combo=ifo_combo)
+bg_far_exc = significance.apply_far_limit(
+    bg_far_exc,
+    significance_dict,
+    combo=ifo_combo)
 
 f['background/ifar'] = conv.sec_to_year(1. / bg_far)
 f['background_exc/ifar'] = conv.sec_to_year(1. / bg_far_exc)
@@ -421,6 +433,17 @@ while numpy.any(ifar_foreground >= background_time):
         return_counts=True,
         **significance_dict[ifo_combo])
 
+    fg_far = significance.apply_far_limit(
+        fg_far,
+        significance_dict,
+        combo=ifo_combo
+    )
+    bg_far = significance.apply_far_limit(
+        bg_far,
+        significance_dict,
+        combo=ifo_combo,
+    )
+
     # Update the ifar_foreground criteria depending on whether foreground
     # triggers are being removed via inclusive or exclusive background.   
     if args.hierarchical_removal_against == 'inclusive':
@@ -435,6 +458,11 @@ while numpy.any(ifar_foreground >= background_time):
             exc_zero_trigs.decimation_factor,
             background_time_exc,
             **significance_dict[ifo_combo])
+        fg_far_exc = significance.apply_far_limit(
+            fg_far_exc,
+            significance_dict,
+            combo=ifo_combo
+        )
         ifar_foreground = 1. / fg_far_exc
     # ifar_foreground has been updated and the code can continue.
 

bin/all_sky_search/pycbc_coinc_statmap_inj

@@ -37,12 +37,21 @@ if args.verbose:
     log_level = logging.INFO
     logging.basicConfig(format='%(asctime)s : %(message)s', level=log_level)
 
-ifo_key = ''.join(args.ifos)
-significance_dict = significance.digest_significance_options([ifo_key], args)
 
 window = args.cluster_window
 logging.info("Loading coinc zerolag triggers")
 zdata = pycbc.io.MultiifoStatmapData(files=args.zero_lag_coincs, ifos=args.ifos)
+
+if 'ifos' in zdata.attrs:
+    ifos = zdata.attrs['ifos'].split(' ')
+    logging.info('using ifos from file {}'.format(args.zero_lag_coincs[0]))
+else:
+    ifos = args.ifos
+    logging.info('using ifos from command line input')
+
+ifo_key = ''.join(ifos)
+significance_dict = significance.digest_significance_options([ifo_key], args)
+
 zdata = zdata.cluster(window)
 
 f = h5py.File(args.output_file, "w")
@@ -51,7 +60,7 @@ f.attrs['num_of_ifos'] = zdata.attrs['num_of_ifos']
 f.attrs['pivot'] = zdata.attrs['pivot']
 f.attrs['fixed'] = zdata.attrs['fixed']
 f.attrs['timeslide_interval'] = zdata.attrs['timeslide_interval']
-f.attrs['ifos'] = ' '.join(sorted(args.ifos))
+f.attrs['ifos'] = ' '.join(sorted(ifos))
 
 # Copy over the segment for coincs and singles
 for key in zdata.seg.keys():
@@ -90,7 +99,11 @@ if len(zdata) > 0:
         background_time,
         **significance_dict[ifo_key])
 
-    fg_far_exc = significance.apply_far_limit(fg_far_exc, significance_dict, combo=ifo_key)
+    fg_far_exc = significance.apply_far_limit(
+        fg_far_exc,
+        significance_dict,
+        combo=ifo_key,
+    )
 
     ifar_exc = 1. / fg_far_exc
     fap_exc = 1 - numpy.exp(- coinc_time / ifar_exc)