page_foreground output to hdf summary file (#3781)

* add hdf output to pycbc_page_foreground

* more outputs in hdf summary file

bin/hdfcoinc/pycbc_page_foreground

@@ -182,3 +182,8 @@ if args.output_file.endswith('.html'):
 elif args.output_file.endswith('.xml') or args.output_file.endswith('.xml.gz'):
     fortrigs.to_coinc_xml_object(args.output_file)
 
+elif args.output_file.endswith('.hdf'):
+    fortrigs.to_coinc_hdf_object(args.output_file)
+
+else:
+    raise NotImplementedError("Output file format not recognised")

pycbc/io/hdf.py

@@ -25,6 +25,8 @@
 from pycbc import events, conversions, pnutils
 from pycbc.events import ranking, veto
 from pycbc.events.stat import sngl_statistic_dict
+from pycbc.events import mean_if_greater_than_zero
+from pycbc.pnutils import mass1_mass2_to_mchirp_eta
 
 class HFile(h5py.File):
     """ Low level extensions to the capabilities of reading an hdf5 File
@@ -986,6 +988,89 @@ def to_coinc_xml_object(self, file_name):
 
         ligolw_utils.write_filename(outdoc, file_name)
 
+    def to_coinc_hdf_object(self, file_name):
+        ofd = h5py.File(file_name,'w')
+
+        # Some fields are special cases:
+        logging.info("Outputting end times")
+        time = self.get_end_time()
+        ofd.create_dataset('end_time', data=time, dtype=np.float32)
+
+        logging.info("Getting IFARs and outputting as FARs")
+        far = 1. / self.get_coincfile_array('ifar')
+        ofd.create_dataset('false_alarm_rate', data=far, dtype=np.float32)
+
+        far_exc = 1. / self.get_coincfile_array('ifar_exc')
+        ofd.create_dataset('false_alarm_rate_exclusive', data=far_exc,
+                           dtype=np.float32)
+
+        logging.info("Outputting p values")
+        fap = self.get_coincfile_array('fap')
+        ofd.create_dataset('p_value', data=fap,
+                           dtype=np.float32)
+
+        fap_exc = self.get_coincfile_array('fap_exc')
+        ofd.create_dataset('p_value_exclusive', data=fap_exc,
+                           dtype=np.float32)
+
+        # Coinc fields
+        logging.info("Outputting coinc info")
+        for field in ['stat']:
+            logging.info(field)
+            vals = self.get_coincfile_array(field)
+            ofd.create_dataset(field, data=vals, dtype=np.float32)
+
+        # Bank fields
+        logging.info("Outputting associated bank values")
+        for field in ['mass1','mass2','spin1z','spin2z']:
+            logging.info(field)
+            vals = self.get_bankfile_array(field)
+            ofd.create_dataset(field, data=vals, dtype=np.float32)
+
+        mass1 = self.get_bankfile_array('mass1')
+        mass2 = self.get_bankfile_array('mass2')
+        mchirp, _ = mass1_mass2_to_mchirp_eta(mass1, mass2)
+
+        ofd.create_dataset('chirp_mass', data=mchirp, dtype=np.float32)
+
+        sngl_fields = [f for f in self.sngl_files[self.ifos[0]].columns
+                       if not ('template' in f)]
+        # Single-detector fields
+        logging.info("Outputting single-detector info")
+        for field in sngl_fields:
+            if field in ['gating','search','snr']: continue
+            vals_valid = self.get_snglfile_array_dict(field)
+            for ifo in self.ifos:
+                vals = vals_valid[ifo][0]
+                valid = vals_valid[ifo][1]
+                vals[np.logical_not(valid)] = -1.
+                ofd.create_dataset(ifo + '_'+field, data=vals,
+                                   dtype=np.float32)
+
+
+        logging.info("Outputting single-detector and network SNRs")
+        snr_vals_valid = self.get_snglfile_array_dict('snr')
+        network_snr_sq = np.zeros_like(snr_vals_valid[self.ifos[0]][0])
+        for ifo in self.ifos:
+            vals = snr_vals_valid[ifo][0]
+            valid = snr_vals_valid[ifo][1]
+            vals[np.logical_not(valid)] = -1.
+            ofd.create_dataset(ifo + '_snr', data=vals,
+                               dtype=np.float32)
+            network_snr_sq[valid] += vals[valid] ** 2.0
+        network_snr = np.sqrt(network_snr_sq)
+        ofd.create_dataset('network_snr', data=network_snr, dtype=np.float32)
+
+        logging.info("Outputting template duration")
+        temp_dur_vals_valid = self.get_snglfile_array_dict('template_duration')
+        temp_dur = np.zeros_like(snr_vals_valid[self.ifos[0]][0])
+        for ifo in self.ifos:
+            vals = temp_dur_vals_valid[ifo][0]
+            valid = temp_dur_vals_valid[ifo][1]
+            temp_dur[valid] = vals[valid]
+        ofd.create_dataset('template_duration', data=temp_dur, dtype=np.float32)
+        ofd.close()
+
 class ReadByTemplate(object):
     # default assignment to {} is OK for a variable used only in __init__
     def __init__(self, filename, bank=None, segment_name=None, veto_files=None,