From 642173ae9a557fa09c2e0297c40e525b0f1d3e60 Mon Sep 17 00:00:00 2001 From: Rahul Dhurkunde Date: Tue, 16 Jun 2026 06:32:32 -0700 Subject: [PATCH 1/6] Implement removal of loud coincs for the two detector case --- pycbc/events/coinc.py | 90 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 84 insertions(+), 6 deletions(-) diff --git a/pycbc/events/coinc.py b/pycbc/events/coinc.py index a36807e1366..7ecf0828e83 100644 --- a/pycbc/events/coinc.py +++ b/pycbc/events/coinc.py @@ -839,6 +839,7 @@ def __init__(self, num_templates, analysis_block, background_statistic, coinc_window_pad=.002, statistic_refresh_rate=None, return_background=False, + ifar_remove_threshold=None, **kwargs): """ Parameters @@ -893,6 +894,10 @@ class (in seconds), default not do do this self.timeslide_interval = timeslide_interval self.return_background = return_background self.coinc_window_pad = coinc_window_pad + self.ifar_remove_threshold = ifar_remove_threshold + # Set of integer chunk indices (gps_time // analysis_block) whose + # triggers are excluded from coincidence formation + self.loud_chunks = set() self.ifos = ifos if len(self.ifos) != 2: @@ -994,6 +999,7 @@ def from_cli(cls, args, num_templates, analysis_chunk, ifos): ifos=ifos, coinc_window_pad=args.coinc_window_pad, statistic_refresh_rate=args.statistic_refresh_rate, + ifar_remove_threshold=args.ifar_remove_threshold, **kwargs) @staticmethod @@ -1010,7 +1016,10 @@ def insert_args(parser): group.add_argument('--timeslide-interval', type=float, help="The interval between timeslides in seconds", default=0.1) group.add_argument('--ifar-remove-threshold', type=float, - help="NOT YET IMPLEMENTED", default=100.0) + help="If a zerolag coincidence has an inverse false alarm rate " + "(in years) above this threshold, the analysis chunks " + "containing its triggers are marked as loud and excluded " + "from background estimation", default=None) @staticmethod def verify_args(args, parser): @@ -1024,8 +1033,12 @@ def verify_args(args, parser): def background_time(self): """Return the amount of background time that the buffers contain""" time = 1.0 / self.timeslide_interval + # Dirty chunks are excluded from coincidence formation in both + # detectors, so they do not contribute to the background time + loud_time = len(self.loud_chunks) * self.analysis_block for ifo in self.singles: - time *= self.singles[ifo].filled_time * self.analysis_block + livetime = self.singles[ifo].filled_time * self.analysis_block + time *= max(livetime - loud_time, 0) return time def save_state(self, filename): @@ -1296,19 +1309,84 @@ def _find_coincs(self, results, valid_ifos): # (both zerolag and shifted are handled together) num_zerolag = 0 num_background = 0 + if len(cstat) > 0: offsets = numpy.concatenate(offsets) ctime0 = numpy.concatenate(ctimes[self.ifos[0]]).astype(numpy.float64) ctime1 = numpy.concatenate(ctimes[self.ifos[1]]).astype(numpy.float64) + good = None + if self.ifar_remove_threshold is not None and self.loud_chunks: + # Prune loud chunks older than the lookback time: their + # triggers have expired from the singles buffers, so they + # must no longer reduce the background time + min_end = max(ctime0.max(), ctime1.max()) - self.lookback_time + self.loud_chunks = { + c for c in self.loud_chunks + if (c + 1) * self.analysis_block > min_end + } + # Remove coincs with a trigger inside a loud chunk in + # either detector before clustering, so they enter neither + # the background nor the foreground + if self.loud_chunks: + loud = numpy.fromiter(self.loud_chunks, + dtype=numpy.int64) + chunk0 = (ctime0 // self.analysis_block).astype(numpy.int64) + chunk1 = (ctime1 // self.analysis_block).astype(numpy.int64) + keep = ~(numpy.isin(chunk0, loud) + | numpy.isin(chunk1, loud)) + good = numpy.flatnonzero(keep) + if len(good) < len(cstat): + logger.info("Removing %d coincs in loud chunks", + len(cstat) - len(good)) + logger.info("Clustering %s coincs", ppdets(self.ifos, "-")) - cidx = cluster_coincs(cstat, ctime0, ctime1, offsets, - self.timeslide_interval, - self.analysis_block + 2*self.time_window, - method='cython') + cluster_window = self.analysis_block + 2 * self.time_window + if good is None: + cidx = cluster_coincs(cstat, ctime0, ctime1, offsets, + self.timeslide_interval, + cluster_window, method='cython') + elif len(good): + cidx = good[cluster_coincs(cstat[good], ctime0[good], + ctime1[good], offsets[good], + self.timeslide_interval, + cluster_window, method='cython')] + else: + cidx = numpy.array([], dtype=numpy.int64) + offsets = offsets[cidx] zerolag_idx = (offsets == 0) bkg_idx = (offsets != 0) + if self.ifar_remove_threshold is not None: + # Mark the chunks containing the triggers of any loud + # zerolag candidate as loud. The candidate itself is still + # reported, but coincs involving these chunks are excluded + # from the background from now on. + new_loud = [] + for idx in cidx[zerolag_idx]: + ifar_val, _ = self.ifar(cstat[idx]) + if ifar_val <= self.ifar_remove_threshold: + continue + chunks = {int(ctime0[idx] // self.analysis_block), + int(ctime1[idx] // self.analysis_block)} + for chunk in chunks - self.loud_chunks: + self.loud_chunks.add(chunk) + new_loud.append(chunk) + logger.info( + "Dirty chunk [%d, %d): zerolag coinc with " + "IFAR %.2f above %.2f", + chunk * self.analysis_block, + (chunk + 1) * self.analysis_block, + ifar_val, self.ifar_remove_threshold + ) + if new_loud: + # Drop this update's background coincs involving the + # newly loud chunks before they enter the buffer + nd = numpy.array(new_loud, dtype=numpy.int64) + tc0 = (ctime0[cidx] // self.analysis_block).astype(numpy.int64) + tc1 = (ctime1[cidx] // self.analysis_block).astype(numpy.int64) + bkg_idx &= ~(numpy.isin(tc0, nd) | numpy.isin(tc1, nd)) + for ifo in self.ifos: single_expire[ifo] = numpy.concatenate(single_expire[ifo]) single_expire[ifo] = single_expire[ifo][cidx][bkg_idx] From 4dbabe9ece67f14354f661b01ba4eb85ad0225ba Mon Sep 17 00:00:00 2001 From: Rahul Dhurkunde Date: Thu, 25 Jun 2026 04:31:21 -0700 Subject: [PATCH 2/6] Keep zero lag coincs --- pycbc/events/coinc.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/pycbc/events/coinc.py b/pycbc/events/coinc.py index 7ecf0828e83..aa842c7b8d5 100644 --- a/pycbc/events/coinc.py +++ b/pycbc/events/coinc.py @@ -1033,7 +1033,7 @@ def verify_args(args, parser): def background_time(self): """Return the amount of background time that the buffers contain""" time = 1.0 / self.timeslide_interval - # Dirty chunks are excluded from coincidence formation in both + # Loud chunks are excluded from coincidence formation in both # detectors, so they do not contribute to the background time loud_time = len(self.loud_chunks) * self.analysis_block for ifo in self.singles: @@ -1324,19 +1324,18 @@ def _find_coincs(self, results, valid_ifos): c for c in self.loud_chunks if (c + 1) * self.analysis_block > min_end } - # Remove coincs with a trigger inside a loud chunk in - # either detector before clustering, so they enter neither - # the background nor the foreground + # Exclude background (timeslide) coincs in loud chunks; + # zerolag coincs are kept so loud signals/injections are + # always reported as candidates. if self.loud_chunks: - loud = numpy.fromiter(self.loud_chunks, - dtype=numpy.int64) chunk0 = (ctime0 // self.analysis_block).astype(numpy.int64) chunk1 = (ctime1 // self.analysis_block).astype(numpy.int64) - keep = ~(numpy.isin(chunk0, loud) - | numpy.isin(chunk1, loud)) - good = numpy.flatnonzero(keep) + loud = numpy.fromiter(self.loud_chunks, dtype=numpy.int64) + in_loud_block = (numpy.isin(chunk0, loud) + | numpy.isin(chunk1, loud)) + good = numpy.flatnonzero(~(in_loud_block & (offsets != 0))) if len(good) < len(cstat): - logger.info("Removing %d coincs in loud chunks", + logger.info("Removing %d background coincs in loud chunks", len(cstat) - len(good)) logger.info("Clustering %s coincs", ppdets(self.ifos, "-")) @@ -1373,7 +1372,7 @@ def _find_coincs(self, results, valid_ifos): self.loud_chunks.add(chunk) new_loud.append(chunk) logger.info( - "Dirty chunk [%d, %d): zerolag coinc with " + "Loud chunk [%d, %d): zerolag coinc with " "IFAR %.2f above %.2f", chunk * self.analysis_block, (chunk + 1) * self.analysis_block, From fabe216e3fcd04fd4de3eae0b3b3670ec9af6aa7 Mon Sep 17 00:00:00 2001 From: Rahul Dhurkunde Date: Mon, 29 Jun 2026 02:59:37 -0700 Subject: [PATCH 3/6] add _filter_loud_coincs method --- pycbc/events/coinc.py | 88 ++++++++++++++++++++++++++----------------- 1 file changed, 54 insertions(+), 34 deletions(-) diff --git a/pycbc/events/coinc.py b/pycbc/events/coinc.py index aa842c7b8d5..a6138813ec2 100644 --- a/pycbc/events/coinc.py +++ b/pycbc/events/coinc.py @@ -1029,12 +1029,49 @@ def verify_args(args, parser): parser.error(f"The single ifo ranking stat {args.sngl_ranking} " "requires --psd-variation.") + def _filter_loud_coincs(self, cstat, ctime0, ctime1, offsets): + """Remove background coincs that fall in loud chunks. + + Prunes stale loud chunks, then returns an index array selecting + only the coincs that are *not* in a loud chunk (background coincs + in loud chunks are excluded; zerolag coincs are always kept). + Returns slice(None) when no filtering is needed so the caller can + treat all cases uniformly. + """ + # Prune loud chunks older than the lookback time: their triggers + # have expired from the singles buffers, so they must no longer + # reduce the background time. + min_end = max(ctime0.max(), ctime1.max()) - self.lookback_time + self.loud_chunks = { + c for c in self.loud_chunks + if (c + 1) * self.analysis_block > min_end + } + if not self.loud_chunks: + return slice(None) + # Exclude background (timeslide) coincs in loud chunks; zerolag + # coincs are kept so loud signals/injections are always reported. + chunk0 = (ctime0 // self.analysis_block).astype(numpy.int64) + chunk1 = (ctime1 // self.analysis_block).astype(numpy.int64) + loud = numpy.fromiter(self.loud_chunks, dtype=numpy.int64) + in_loud_block = numpy.isin(chunk0, loud) | numpy.isin(chunk1, loud) + good = numpy.flatnonzero(~(in_loud_block & (offsets != 0))) + if len(good) < len(cstat): + logger.info( + "Removing %d background coincs in loud chunks", + len(cstat) - len(good), + ) + return good + @property def background_time(self): - """Return the amount of background time that the buffers contain""" + """Return the amount of background time that the buffers contain. + + A loud chunk is an analysis_block-length time segment identified as + containing a loud candidate (IFAR above ifar_remove_threshold). Loud + chunks are excluded from background coincidence formation in both + detectors, so they do not contribute to the background time. + """ time = 1.0 / self.timeslide_interval - # Loud chunks are excluded from coincidence formation in both - # detectors, so they do not contribute to the background time loud_time = len(self.loud_chunks) * self.analysis_block for ifo in self.singles: livetime = self.singles[ifo].filled_time * self.analysis_block @@ -1314,41 +1351,24 @@ def _find_coincs(self, results, valid_ifos): offsets = numpy.concatenate(offsets) ctime0 = numpy.concatenate(ctimes[self.ifos[0]]).astype(numpy.float64) ctime1 = numpy.concatenate(ctimes[self.ifos[1]]).astype(numpy.float64) - good = None + good = slice(None) if self.ifar_remove_threshold is not None and self.loud_chunks: - # Prune loud chunks older than the lookback time: their - # triggers have expired from the singles buffers, so they - # must no longer reduce the background time - min_end = max(ctime0.max(), ctime1.max()) - self.lookback_time - self.loud_chunks = { - c for c in self.loud_chunks - if (c + 1) * self.analysis_block > min_end - } - # Exclude background (timeslide) coincs in loud chunks; - # zerolag coincs are kept so loud signals/injections are - # always reported as candidates. - if self.loud_chunks: - chunk0 = (ctime0 // self.analysis_block).astype(numpy.int64) - chunk1 = (ctime1 // self.analysis_block).astype(numpy.int64) - loud = numpy.fromiter(self.loud_chunks, dtype=numpy.int64) - in_loud_block = (numpy.isin(chunk0, loud) - | numpy.isin(chunk1, loud)) - good = numpy.flatnonzero(~(in_loud_block & (offsets != 0))) - if len(good) < len(cstat): - logger.info("Removing %d background coincs in loud chunks", - len(cstat) - len(good)) + good = self._filter_loud_coincs(cstat, ctime0, ctime1, offsets) logger.info("Clustering %s coincs", ppdets(self.ifos, "-")) cluster_window = self.analysis_block + 2 * self.time_window - if good is None: - cidx = cluster_coincs(cstat, ctime0, ctime1, offsets, - self.timeslide_interval, - cluster_window, method='cython') - elif len(good): - cidx = good[cluster_coincs(cstat[good], ctime0[good], - ctime1[good], offsets[good], - self.timeslide_interval, - cluster_window, method='cython')] + good_cstat = cstat[good] + if len(good_cstat): + sub = cluster_coincs( + good_cstat, + ctime0[good], + ctime1[good], + offsets[good], + self.timeslide_interval, + cluster_window, + method='cython', + ) + cidx = numpy.arange(len(cstat), dtype=numpy.int64)[good][sub] else: cidx = numpy.array([], dtype=numpy.int64) From 10edbae2dcfeca537a6ec43425bed40e05aaabfd Mon Sep 17 00:00:00 2001 From: Rahul Dhurkunde Date: Mon, 29 Jun 2026 04:16:36 -0700 Subject: [PATCH 4/6] simplify loud chunk loop using continue guard --- pycbc/events/coinc.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/pycbc/events/coinc.py b/pycbc/events/coinc.py index a6138813ec2..9ba864e1591 100644 --- a/pycbc/events/coinc.py +++ b/pycbc/events/coinc.py @@ -1386,9 +1386,15 @@ def _find_coincs(self, results, valid_ifos): ifar_val, _ = self.ifar(cstat[idx]) if ifar_val <= self.ifar_remove_threshold: continue + # Both times are within the light travel time of each + # other (zerolag), so this set almost always has one + # element; two elements only if a trigger straddles a + # block boundary. chunks = {int(ctime0[idx] // self.analysis_block), int(ctime1[idx] // self.analysis_block)} - for chunk in chunks - self.loud_chunks: + for chunk in chunks: + if chunk in self.loud_chunks: + continue self.loud_chunks.add(chunk) new_loud.append(chunk) logger.info( From 29b871fb1c51c887a644fc423ad7134dcea39f53 Mon Sep 17 00:00:00 2001 From: Rahul Dhurkunde Date: Mon, 29 Jun 2026 07:15:14 -0700 Subject: [PATCH 5/6] Implement veto-window for triggers straddling near chunk boundaries --- pycbc/events/coinc.py | 47 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 41 insertions(+), 6 deletions(-) diff --git a/pycbc/events/coinc.py b/pycbc/events/coinc.py index 9ba864e1591..21527ae40ef 100644 --- a/pycbc/events/coinc.py +++ b/pycbc/events/coinc.py @@ -829,6 +829,36 @@ def data(self): return self.buffer[:self.index] +def chunk_indices_with_boundary(times, analysis_block, boundary_window): + """Return the set of chunk indices covering *times*, expanding into + neighbouring chunks when a time falls within *boundary_window* seconds + of a chunk edge. + + Parameters + ---------- + times : iterable of float + GPS trigger or injection times. + analysis_block : int or float + Chunk duration in seconds. + boundary_window : float + If a time is within this many seconds of a chunk boundary, the + adjacent chunk index is also included. + + Returns + ------- + set of int + """ + indices = set() + for t in times: + c = int(t // analysis_block) + indices.add(c) + if t - c * analysis_block < boundary_window: + indices.add(c - 1) + if (c + 1) * analysis_block - t < boundary_window: + indices.add(c + 1) + return indices + + class LiveCoincTimeslideBackgroundEstimator(object): """Rolling buffer background estimation.""" @@ -840,6 +870,7 @@ def __init__(self, num_templates, analysis_block, background_statistic, statistic_refresh_rate=None, return_background=False, ifar_remove_threshold=None, + boundary_veto_window=0.1, **kwargs): """ Parameters @@ -872,6 +903,10 @@ class (in seconds), default not do do this return_background: boolean If true, background triggers will also be included in the file output. + boundary_veto_window: float + If a loud trigger falls within this many seconds of a chunk + boundary, the neighbouring chunk is also flagged as loud. + Default 0.1 s. Applies to both IFAR-based and injection vetoes. kwargs: dict Additional options for the statistic to use. See stat.py for more details on statistic options. @@ -895,6 +930,7 @@ class (in seconds), default not do do this self.return_background = return_background self.coinc_window_pad = coinc_window_pad self.ifar_remove_threshold = ifar_remove_threshold + self.boundary_veto_window = boundary_veto_window # Set of integer chunk indices (gps_time // analysis_block) whose # triggers are excluded from coincidence formation self.loud_chunks = set() @@ -1386,12 +1422,11 @@ def _find_coincs(self, results, valid_ifos): ifar_val, _ = self.ifar(cstat[idx]) if ifar_val <= self.ifar_remove_threshold: continue - # Both times are within the light travel time of each - # other (zerolag), so this set almost always has one - # element; two elements only if a trigger straddles a - # block boundary. - chunks = {int(ctime0[idx] // self.analysis_block), - int(ctime1[idx] // self.analysis_block)} + chunks = chunk_indices_with_boundary( + [ctime0[idx], ctime1[idx]], + self.analysis_block, + self.boundary_veto_window, + ) for chunk in chunks: if chunk in self.loud_chunks: continue From d091ee20176ea644c250921d80b5c4afa8b51121 Mon Sep 17 00:00:00 2001 From: Rahul Dhurkunde Date: Wed, 1 Jul 2026 04:30:55 -0700 Subject: [PATCH 6/6] add ifar-remove-threshold to args --- test/test_live_coinc_compare.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_live_coinc_compare.py b/test/test_live_coinc_compare.py index fed36e8d9e8..e648bc845fc 100644 --- a/test/test_live_coinc_compare.py +++ b/test/test_live_coinc_compare.py @@ -79,6 +79,7 @@ def setUp(self, *args): store_background=True, coinc_window_pad=0.002, statistic_refresh_rate=None, + ifar_remove_threshold=None, ) # number of templates in the bank