Trigger timeseries ranking (gwastro#4491)

* allow all rankings in pycbc_plot_trigger_timeseries

* remove uinrelated chnage

* Update bin/minifollowups/pycbc_plot_trigger_timeseries

Co-authored-by: Tito Dal Canton <[email protected]>

* help message fix

* Some more safety about unwanted events going away from the plot

* more logging

---------

Co-authored-by: Tito Dal Canton <[email protected]>

bin/minifollowups/pycbc_plot_trigger_timeseries

@@ -19,6 +19,7 @@ import argparse, logging, pycbc.version, pycbc.results, sys
 from pycbc.types import MultiDetOptionAction
 from pycbc.events import ranking
 from pycbc.io import HFile
+import h5py
 import matplotlib; matplotlib.use('Agg'); import pylab
 import numpy
 
@@ -39,8 +40,9 @@ parser.add_argument('--times', nargs='+', type=float,
 parser.add_argument('--special-trigger-ids', nargs='+', type=int,
     action=MultiDetOptionAction, metavar="IFO:GPS_TIME",
     help="The set of special trigger ids to plot a star at")
-parser.add_argument('--plot-type', choices=['snr', 'newsnr'], default='snr',
-    help="Which plot to make; an 'snr' or a newsnr' plot.")
+parser.add_argument('--plot-type',
+    choices=ranking.sngls_ranking_function_dict, default='snr',
+    help="Which single-detector ranking statistic to plot.")
 parser.add_argument('--output-file')
 parser.add_argument('--log-y-axis', action='store_true')
 
@@ -49,73 +51,80 @@ pycbc.init_logging(args.verbose)
 
 any_data = False
 
-# organize the single detector triggers files by ifo in a dict
 fig = pylab.figure()
+
+min_rank = numpy.inf
 for ifo in args.single_trigger_files.keys():
+    logging.info("Getting %s triggers", ifo)
     t = args.times[ifo]
 
-    if args.special_trigger_ids:
-        id_loud = args.special_trigger_ids[ifo]
-
     data = HFile(args.single_trigger_files[ifo], 'r')
 
     # Identify trigger indices within window
-    select_func = lambda *inputs: ((inputs[0] < (t + args.window)) & 
-                                   (inputs[0] > (t - args.window)))
-    times, snr, chisq, chisq_dof = data.select(select_func, ifo + '/end_time',
-                                               ifo + '/snr', ifo + '/chisq',
-                                               ifo + '/chisq_dof')
-                                 
-    # center times on the trigger/chosen time
-    times = times - t
-
-    if args.plot_type == 'snr':
-        val = snr
-
-        if args.special_trigger_ids:
-            top = data[ifo]['snr'][id_loud]
-
-    if args.plot_type == 'newsnr':
-        rchisq = chisq / (2 * chisq_dof - 2)
-        if len(rchisq) > 0:
-            val = ranking.newsnr(snr, rchisq)
-        else:
-            val = numpy.array([])
-
-        # Get the newnsr of the loudest trigger so we can plot a star there
-        if args.special_trigger_ids:
-            rchisq = data[ifo]['chisq'][id_loud] / \
-                (data[ifo]['chisq_dof'][id_loud] * 2 - 2)
-            top = ranking.newsnr(data[ifo]['snr'][id_loud], rchisq)
-
-    if type(val) in [numpy.float32, numpy.float64] or len(val) > 0:
-        any_data = True
-
-    pylab.scatter(times, val, color=pycbc.results.ifo_color(ifo), marker='x',
+    keep = abs(data[ifo + '/end_time'][:] - t) < args.window
+    if not any(keep):
+        # No triggers in this window, add to the legend and continue
+        # Make sure it isnt on the plot
+        pylab.scatter(-2 * args.window, 0,
+                      color=pycbc.results.ifo_color(ifo),
+                      marker='x',
+                      label=ifo)
+        continue
+
+    any_data = True
+    logging.info("Keeping %d triggers in the window", sum(keep))
+
+
+    data_dict = {k: data[ifo][k][keep] for k in data[ifo]
+                 if isinstance(data[ifo][k], h5py.Dataset)
+                 and data[ifo][k].size == keep.size}
+
+    logging.info("Getting %s", args.plot_type)
+    rank = ranking.get_sngls_ranking_from_trigs(data_dict, args.plot_type)
+
+    if all(rank == 1):
+        # This is the default value to say "downranked beyond consideration"
+        # so skip these events
+        pylab.scatter(-2 * args.window, 0,
+                      color=pycbc.results.ifo_color(ifo),
+                      marker='x',
+                      label=ifo)
+        continue
+
+    pylab.scatter(data_dict['end_time'] - t, rank,
+                  color=pycbc.results.ifo_color(ifo), marker='x',
                   label=ifo)
 
+    # Minimum rank which hasn't been set to the default of 1
+    min_rank = min(min_rank, rank[rank > 1].min())
+
     if args.special_trigger_ids:
-        any_data = True
-        pylab.scatter([0], [top], marker='*', s=50, color='yellow')
+        special_idx = numpy.flatnonzero(keep) == args.special_trigger_ids[ifo]
+        if not any(special_idx):
+            logging.info("IDX %d not in kept list", args.special_trigger_ids[ifo])
+            continue
+
+        pylab.scatter((data_dict['end_time'] - t)[special_idx],
+                      rank[special_idx], marker='*', s=50, color='yellow')
 
 if args.log_y_axis and any_data:
     pylab.yscale('log')
 
+if not numpy.isinf(min_rank):
+    pylab.ylim(ymin=min_rank)
+
 pylab.xlabel('time (s)')
 pylab.ylabel(args.plot_type)
 
-try:
-    if len(val) > 0:
-        pylab.ylim(ymin=val.min())
-except TypeError:
-    pylab.ylim(ymin=val)
-
 pylab.xlim(xmin=-args.window, xmax=args.window)
 pylab.legend()
 pylab.grid()
+
+logging.info("Saving figure")
 pycbc.results.save_fig_with_metadata(fig, args.output_file,
             cmd = ' '.join(sys.argv),
             title = 'Single Detector Trigger Timeseries (%s)' % args.plot_type,
             caption = 'Time series showing the single detector triggers'
                       ' centered around the time of the trigger of interest.',
          )
+logging.info("Done!")