Start adding support for idq based reweighting.

Includes argument for whether to calculate dq penalty,
and stat that utilizes dq penalty

bin/pycbc_live

@@ -29,7 +29,7 @@ from pycbc.filter import LiveBatchMatchedFilter, compute_followup_snr_series
 from pycbc.filter import followup_event_significance
 from pycbc.strain import StrainBuffer
 from pycbc.events.ranking import newsnr
-from pycbc.events.coinc import LiveCoincTimeslideBackgroundEstimator as Coincer, MultiRingBuffer
+from pycbc.events.coinc import LiveCoincTimeslideBackgroundEstimator as Coincer
 from pycbc.events.single import LiveSingle
 from pycbc.io.live import CandidateForGraceDB
 from pycbc.io.hdf import recursively_save_dict_contents_to_group
@@ -817,8 +817,11 @@ parser.add_argument('--idq-state-channel', action=MultiDetMultiColonOptionAction
                     help="Channel containing information about the state of idq."
                          "Used to determine when idq is usable for statistics. ")
 parser.add_argument('--idq-threshold', type=float,
-                    help='Threshold used to veto triggers at times of '
-                    'low iDQ False Alarm Probability')
+                    help="Threshold used to flag triggers at times of "
+                         "low iDQ False Alarm Probability")
+parser.add_argument('--idq-reweighting', type=str,
+                    help="Path to file with idq reweighting. If not given idq is"
+                         "used to veto. Ignored if no iDQ channel given.")
 parser.add_argument('--data-quality-channel',
                     action=MultiDetMultiColonOptionAction,
                     help="Channel containing data quality information. Used "
@@ -1230,7 +1233,7 @@ with ctx:
                 results.pop(ifo)
 
             # Veto single detector triggers that fail the DQ vector
-            # and record idq value for each trigger
+            # and process dq data for each trigger
             for ifo in results:
                 if data_reader[ifo].dq:
                     logging.info("Checking %s's DQ vector", ifo)
@@ -1246,14 +1249,22 @@ with ctx:
                             results[ifo][key] = results[ifo][key][idx]
                 if data_reader[ifo].idq:
                     logging.info("Reading %s's iDQ information", ifo)
+                    start = data_reader[ifo].start_time
                     times = results[ifo]['end_time']
-                    idq_vals = data_reader[ifo].idq.lookup_idq(times)
-                    idq_state = data_reader[ifo].idq_state.lookup_idq(times)
-                    # idq state will be 1 if idq is working, 0 if not
-                    # if idq state is 0, set val to nan
-                    # statistic must implement case for nan values
-                    idq_vals[idq_state==0] = numpy.nan
-                    results[ifo]['idq'] = idq_vals
+                    if args.idq_reweighting:
+                        # FIXME implement the following function
+                        results[ifo]['dq_penalty'] = data_reader[ifo].idq.dq_penalty(
+                                start, valid_pad, times)
+                    else:
+                        # use idq as a veto
+                        idx = data_reader[ifo].idq.indices_of_flag(
+                                start, valid_pad, times,
+                                padding=data_reader[ifo].dq_padding)
+                        logging.info('Keeping %d/%d %s triggers after iDQ',
+                                    len(idx), len(times), ifo)
+                        for key in results[ifo]:
+                            if len(results[ifo][key]):
+                                results[ifo][key] = results[ifo][key][idx]
 
             # Calculate and add the psd variation for the results
             if args.psd_variation:

pycbc/events/ranking.py

@@ -127,15 +127,15 @@ def newsnr_sgveto_psdvar_scaled_threshold(snr, bchisq, sgchisq, psd_var_val,
     else:
         return nsnr[0]
 
-def newsnr_sgveto_idq(snr, brchisq, sgchisq, idq):
+def newsnr_sgveto_idq(snr, brchisq, sgchisq, dq_penalty):
     """ Effective SNR with downweighting based on iDQ False Alarm Probability
     """
     nsnr_sg = newsnr_sgveto(snr, brchisq, sgchisq)
     # nan values in idq indicate idq state was not ok
     # replace nan with 1 so no downweighting occurs
     idq = numpy.nan_to_num(idq, nan=1.0)
-    nsnr_sqr = nsnr_sg**2 + 2*numpy.log(idq)
-    nsnr_sqr[nsnr_sqr<0] = 0
+    nsnr_sqr = nsnr_sg**2 - 2*numpy.log(dq_penalty)
+    nsnr_sqr[nsnr_sqr < 0] = 0
     nsnr = numpy.sqrt(nsnr_sqr)
     if hasattr(snr, '__len__'):
         return nsnr
@@ -301,7 +301,7 @@ def get_newsnr_sgveto_psdvar_scaled_threshold(trigs):
                      trigs['psd_var_val'][:])
     return numpy.array(nsnr_sg_psdt, ndmin=1, dtype=numpy.float32)
 
-def get_newsnr_sgveto_idq(trigs):
+def get_newsnr_sgveto_dq(trigs):
     """
     Calculate newsnr re-weigthed by the sine-gaussian veto
     using the log of the iDQ false alarm probability
@@ -316,11 +316,11 @@ def get_newsnr_sgveto_idq(trigs):
         Array of newsnr values
     """
     dof = 2. * trigs['chisq_dof'][:] - 2.
-    if'idq' in trigs.keys():
+    if'dq_penalty' in trigs.keys():
         nsnr_sg_idq = newsnr_sgveto_idq(trigs['snr'][:],
                             trigs['chisq'][:] / dof,
                             trigs['sg_chisq'][:],
-                            trigs['idq'][:])
+                            trigs['dq_penalty'][:])
     else:
         nsnr_sg_idq = newsnr_sgveto(trigs['snr'][:],
                             trigs['chisq'][:] / dof,
@@ -336,7 +336,7 @@ def get_newsnr_sgveto_idq(trigs):
     'newsnr_sgveto_psdvar_threshold': get_newsnr_sgveto_psdvar_threshold,
     'newsnr_sgveto_psdvar_scaled': get_newsnr_sgveto_psdvar_scaled,
     'newsnr_sgveto_psdvar_scaled_threshold': get_newsnr_sgveto_psdvar_scaled_threshold,
-    'newsnr_sgveto_idq': get_newsnr_sgveto_idq,
+    'newsnr_sgveto_idq': get_newsnr_sgveto_dq,
 }