Cherry-pick for v2.8.3 (#5128)

* Make pycbc_multi_inspiral way less verbose (#5107)

* Make pycbc_multi_inspiral way less verbose

* Silence more bank veto logging

* Ensure that the PyGRB offsource segmentlist contains ints (#5108)

* Ensure that the offsource segmentlist has ints

* Cleaner soluziont

* Comment

* Change recursion parameter in resample (#5113)

* Expand multi_inspiral examples and tests (#5114)

* Cleanup and simplify examples

* Add 170817 cases with 1 and 2 detectors

* Handle different detectors in check

* Update test suite

* PyGRB: handling the 2 IFO case with the coherent search methods (#5115)

* PyGRB single IFO case fixes (#5118)

* Copy gating information as is when combining trigger files in PyGRB

* Missing verb in comment

* Apply coinc SNR cut even in single IFO PyGRB runs

* Comment update

* Comment update

* Missing help string

* Workaround for git.ligo.org authentication change (#5119)

* Workaround for git.ligo.org auth change

* Workaround for git.ligo.org auth change

* Forgot args

* Fix error

* Two possible URL schemes

* Fix regex

* Set version to 2.8.3

* updated for depcreated fillcolor value (#5123)

---------

Co-authored-by: Francesco Pannarale <[email protected]>
Co-authored-by: Ian Harry <[email protected]>
Co-authored-by: Alex Nitz <[email protected]>

Author: 4 people

bin/pycbc_multi_inspiral

@@ -624,7 +624,7 @@ with ctx:
                 if not inj_filter_rejector[ifo].template_segment_checker(
                     bank, t_num, stilde[ifo]
                 ):
-                    logging.info(
+                    logging.debug(
                         "Skipping segment %d/%d with template %d/%d as no "
                         "detectable injection is present",
                         s_num + 1,
@@ -648,7 +648,11 @@ with ctx:
             if not analyse_segment:
                 continue
             logging.info(
-                "Analyzing segment %d/%d", s_num + 1, len(segments[ifo])
+                "Analyzing template %d/%d, segment %d/%d",
+                t_num + 1,
+                n_bank,
+                s_num + 1,
+                len(segments[ifo])
             )
             # The following dicts with IFOs as keys are created to store
             # copies of the matched filtering results computed below.
@@ -664,7 +668,7 @@ with ctx:
             # - The list of normalized SNR values at the trigger locations
             snr = dict.fromkeys(args.instruments)
             for ifo in args.instruments:
-                logging.info(
+                logging.debug(
                     "Filtering template %d/%d, ifo %s", t_num + 1, n_bank, ifo
                 )
                 # The following lines unpack and store copies of the matched
@@ -696,10 +700,10 @@ with ctx:
 
             # Loop over (short) time-slides, starting with the zero-lag
             for slide in range(num_slides):
-                logging.info("Analyzing slide %d/%d", slide, num_slides)
+                logging.debug("Analyzing slide %d/%d", slide, num_slides)
                 # Loop over sky positions
                 for position_index in range(len(sky_positions)):
-                    logging.info(
+                    logging.debug(
                         "Analyzing sky position %d/%d",
                         position_index + 1,
                         len(sky_positions),
@@ -741,7 +745,7 @@ with ctx:
                                 args.instruments[0]
                             ]
                         )
-                    logging.info(
+                    logging.debug(
                         "Found %d coincident triggers", len(coinc_idx)
                     )
                     # Time delay is applied to indices to have them at the IFOs
@@ -753,9 +757,14 @@ with ctx:
                         % len(snr_dict[ifo])
                         for ifo in args.instruments
                     }
-                    # Calculate the coincident and coherent SNR.
-                    # First check there is enough data to compute the SNRs.
-                    if len(coinc_idx) != 0 and nifo > 1:
+                    # Calculate the quadrature sum of individual detector SNRs,
+                    # regardless of how many IFOs are available.
+                    # This accounts for a coherent SNR component and an
+                    # incoherent one.  It is referred to as coincident SNR
+                    # throughout the code, even in the single IFO case.
+                    # The calculation is performed only if triggers were
+                    # gathered so far.
+                    if len(coinc_idx) != 0:
                         # Find coinc SNR at trigger times and apply coinc SNR
                         # threshold (which depopulates coinc_idx accordingly)
                         (
@@ -768,32 +777,24 @@ with ctx:
                             args.coinc_threshold,
                             time_delay_idx[slide][position_index],
                         )
-                        logging.info(
+                        logging.debug(
                             "%d triggers above coincident SNR threshold",
                             len(coinc_idx),
                         )
                         if len(coinc_idx) != 0:
-                            logging.info(
+                            logging.debug(
                                 "With max coincident SNR = %.2f",
                                 max(rho_coinc),
                             )
-                    # If there is only one IFO, just take its triggers
-                    # and their SNRs
-                    elif len(coinc_idx) != 0 and nifo == 1:
-                        coinc_triggers = {
-                            args.instruments[0]: snr_dict[args.instruments[0]][
-                                coinc_idx_det_frame[args.instruments[0]]
-                            ]
-                        }
                     else:
                         coinc_triggers = {}
-                        logging.info("No coincident triggers were found")
+                        logging.debug("No coincident triggers were found")
                     # If there are triggers above coinc threshold and more
-                    # than 2 IFOs, then calculate the coherent statistics for
+                    # than 1 IFO, then calculate the coherent statistics for
                     # them and apply the cut on coherent SNR (with threshold
                     # equal to the coinc SNR one)
-                    if len(coinc_idx) != 0 and nifo > 2:
-                        logging.info("Calculating their coherent statistics")
+                    if len(coinc_idx) != 0 and nifo > 1:
+                        logging.debug("Calculating their coherent statistics")
                         # Plus and cross antenna pattern dictionaries
                         fp = {
                             ifo: antenna_pattern[ifo][position_index][0]
@@ -893,12 +894,12 @@ with ctx:
                                 project,
                                 rho_coinc,
                             )
-                        logging.info(
+                        logging.debug(
                             "%d triggers above coherent SNR threshold",
                             len(rho_coh),
                         )
                         if len(coinc_idx) != 0:
-                            logging.info(
+                            logging.debug(
                                 "With max coherent SNR = %.2f", max(rho_coh)
                             )
                             # Calculate the null SNR and apply the null SNR cut
@@ -918,11 +919,11 @@ with ctx:
                                 snrv=coinc_triggers,
                                 index=coinc_idx,
                             )
-                            logging.info(
+                            logging.debug(
                                 "%d triggers above null threshold", len(null)
                             )
                             if len(coinc_idx) != 0:
-                                logging.info(
+                                logging.debug(
                                     "With max null SNR = %.2f", max(null)
                                 )
                     # Now calculate the individual detector chi2 values
@@ -967,7 +968,7 @@ with ctx:
                             chisq, chisq_dof, coherent_ifo_trigs
                         )
                         # Calculate chisq reweighted SNR
-                        if nifo > 2:
+                        if nifo > 1:
                             reweighted_snr = ranking.newsnr(
                                 rho_coh,
                                 network_chisq_dict,
@@ -983,13 +984,6 @@ with ctx:
                                 null_grad=args.null_grad,
                                 null_step=args.null_step,
                             )
-                        elif nifo == 2:
-                            reweighted_snr = ranking.newsnr(
-                                rho_coinc,
-                                network_chisq_dict,
-                                q=args.chisq_index,
-                                n=args.chisq_nhigh,
-                            )
                         else:
                             rho_sngl = abs(
                                 snr_dict[args.instruments[0]][
@@ -1050,11 +1044,9 @@ with ctx:
                                 ifo_names,
                                 [ifo_out_vals[n] for n in ifo_names],
                             )
-                        if nifo > 2:
+                        if nifo > 1:
                             network_out_vals['coherent_snr'] = rho_coh
                             network_out_vals['null_snr'] = null
-                        elif nifo == 2:
-                            network_out_vals['coherent_snr'] = rho_coinc
                         else:
                             network_out_vals['coherent_snr'] = abs(
                                 snr_dict[args.instruments[0]][
@@ -1089,7 +1081,7 @@ with ctx:
                 cluster_window = int(template.chirp_length * sample_rate)
             # Cluster template events by slide
             for slide in range(num_slides):
-                logging.info("Clustering slide %d", slide)
+                logging.debug("Clustering slide %d", slide)
                 event_mgr.cluster_template_network_events(
                     'time_index', 'reweighted_snr', cluster_window, slide=slide
                 )

bin/pygrb/pycbc_grb_trig_combiner

@@ -137,7 +137,7 @@ def merge_hdf5_files(inputfiles, outputfile, verbose=False, **compression_kw):
     all_event_id_names = sngl_event_id_names + network_event_id_names
 
     # handle search datasets as a special case
-    # (they will the same in all files)
+    # (they will be the same in all files)
     search_datasets = set(filter(lambda x: "/search/" in x, dataset_names))
     gating_datasets = set(filter(lambda x: "/gating/" in x, dataset_names))
     once_only_datasets = search_datasets.union(gating_datasets)

examples/multi_inspiral/check_gw170817_trigs.py

@@ -1,37 +1,46 @@
 #!/usr/bin/env python
+
 # Read a pycbc_multi_inspiral HDF5 trigger file and check that it contains
-# triggers compatible with mock GW170817-like injections
+# triggers compatible with GW170817
 # 2022 Andrew Williamson, Tito Dal Canton
 
 import sys
 import logging
 import h5py
-from pycbc import init_logging
 
-init_logging(True)
-gw170817_time = 1187008882.43
-status = 0
+
 with h5py.File('GW170817_test_output.hdf', 'r') as f:
-    snrs = [
-        f['network/end_time_gc'][:],
-        f['network/coherent_snr'][:],
-        f['network/reweighted_snr'][:],
-        f['network/slide_id'][:]]
-# search for compatible trigs
-mask = (
-    (abs(gw170817_time - snrs[0]) < 0.1)
-    & (snrs[1] > 25)
-    & (snrs[2] > 25)
-    & (snrs[3] == 0)
+    end_time = f['network/end_time_gc'][:]
+    coh_snr = f['network/coherent_snr'][:]
+    rw_snr = f['network/reweighted_snr'][:]
+    slide_id = f['network/slide_id'][:]
+    present_detectors = ''.join(
+        d for d in 'HLV' if f'{d}1' in f.keys()
     )
+
+# basic sanity checks
+assert (end_time > 0).all()
+assert (coh_snr > 0).all()
+assert (rw_snr > 0).all()
+
+# search for trigs compatible with GW170817
+expected_time = 1187008882.43
+required_net_snr = {'H': 17, 'HL': 28, 'HLV': 28}
+mask = (
+    (abs(expected_time - end_time) < 0.1)
+    & (coh_snr > required_net_snr[present_detectors])
+    & (rw_snr > required_net_snr[present_detectors])
+    & (slide_id == 0)
+)
 n = mask.sum()
 if n > 0:
-    result = 'PASS'
+    print(
+        f'PASS: GW170817 found with max coherent SNR {max(coh_snr[mask]):.2f}, '
+        f'max reweighted SNR {max(rw_snr[mask]):.2f}'
+    )
     status = 0
 else:
-    result = 'FAIL'
+    print('FAIL: GW170817 not found')
     status = 1
-logging.info(
-    '%s: GW170817 found with coherent SNR = %.2f; reweighted SNR %.2f', result,
-    snrs[1][mask], snrs[2][mask])
+
 sys.exit(status)

examples/multi_inspiral/clean.sh

@@ -1 +1 @@
-rm -f bank_veto_bank.xml *.hdf *.gwf
+rm -f bank_veto_bank.xml *.hdf

examples/multi_inspiral/faceon_faceaway.sh

@@ -32,8 +32,8 @@ PAD=8
 BANK_FILE=gw170817_single_template.hdf
 BANK_VETO_FILE=bank_veto_bank.xml
 CHANNEL=SIMULATED_GW170817
-RA=5.016076270234897
-DEC=-0.408407044967
+RA="5.016076270234897 rad"
+DEC="-0.408407044967 rad"
 
 echo -e "\\n\\n>> [`date`] Getting template bank"
 wget -nv -nc ${CONFIG_URL}/${BANK_FILE}

examples/multi_inspiral/gw170817_h.sh

@@ -0,0 +1,62 @@
+#!/bin/bash
+
+set -e
+
+CONFIG_URL=https://github.com/gwastro/pycbc-config/raw/master/test/multi_inspiral
+BANK_FILE=gw170817_single_template.hdf
+echo -e "\\n\\n>> [`date`] Getting template bank"
+wget -nv -nc ${CONFIG_URL}/${BANK_FILE}
+
+EVENT=1187008882
+PAD=8
+START_PAD=111
+END_PAD=17
+GPS_START=$((EVENT - 192 - PAD))
+GPS_END=$((EVENT + 192 + PAD))
+TRIG_START=$((GPS_START + START_PAD))
+TRIG_END=$((GPS_END - END_PAD))
+OUTPUT=GW170817_test_output.hdf
+
+echo -e "\\n\\n>> [`date`] Running pycbc_multi_inspiral on GW170817 data"
+pycbc_multi_inspiral \
+    --verbose \
+    --projection left+right \
+    --instruments H1 \
+    --channel-name H1:GWOSC-16KHZ_R1_STRAIN \
+    --frame-type H1:GWOSC \
+    --trigger-time ${EVENT} \
+    --gps-start-time ${GPS_START} \
+    --gps-end-time ${GPS_END} \
+    --trig-start-time ${TRIG_START} \
+    --trig-end-time ${TRIG_END} \
+    --ra '3.44527994344 rad' \
+    --dec '-0.408407044967 rad' \
+    --bank-file ${BANK_FILE} \
+    --approximant IMRPhenomD \
+    --order -1 \
+    --low-frequency-cutoff 30 \
+    --sngl-snr-threshold 3.0 \
+    --chisq-bins "0.9*get_freq('fSEOBNRv4Peak',params.mass1,params.mass2,params.spin1z,params.spin2z)**(2./3.)" \
+    --pad-data 8 \
+    --strain-high-pass 25 \
+    --sample-rate 4096 \
+    --autogating-threshold 100 \
+    --autogating-cluster 0.5 \
+    --autogating-width 0.25 \
+    --autogating-taper 0.25 \
+    --autogating-pad 0 \
+    --cluster-method window \
+    --cluster-window 0.1 \
+    --segment-length 256 \
+    --segment-start-pad ${START_PAD} \
+    --segment-end-pad ${END_PAD} \
+    --psd-estimation median \
+    --psd-segment-length 32 \
+    --psd-segment-stride 8 \
+    --psd-num-segments 29 \
+    --do-shortslides \
+    --slide-shift 1 \
+    --output ${OUTPUT}
+
+echo -e "\\n\\n>> [`date`] Checking output files"
+python check_gw170817_trigs.py