From 01fbf11d1584556867ca427d050f477e00a57773 Mon Sep 17 00:00:00 2001
From: Kipp Cannon <kipp.cannon@ligo.org>
Date: Sat, 25 Aug 2018 17:10:18 +0900
Subject: [PATCH] far.py: max-over-subset ln L

- redefine the ranking statistic to be the largest ln L for any subset of the triggers in a candidate to prevent a glitch in one detector from making a good candidate look bad
---
 gstlal-inspiral/python/far.py | 31 ++++++++++++++++++++++++++++---
 1 file changed, 28 insertions(+), 3 deletions(-)

diff --git a/gstlal-inspiral/python/far.py b/gstlal-inspiral/python/far.py
index a4621a322f..9d460ce57d 100644
--- a/gstlal-inspiral/python/far.py
+++ b/gstlal-inspiral/python/far.py
@@ -90,6 +90,29 @@ from gstlal.stats import inspiral_lr
 #
 
 
+def kwarggeniter(d, min_instruments):
+	d = tuple(sorted(d.items()))
+	return map(dict, itertools.chain(*(itertools.combinations(d, i) for i in range(min_instruments, len(d) + 1))))
+
+
+def kwarggen(segments, snrs, chi2s_over_snr2s, phase, dt, template_id, min_instruments):
+	# segments and template_id held fixed
+	for snrs, chi2s_over_snr2s, phase, dt in zip(
+		kwarggeniter(snrs, min_instruments),
+		kwarggeniter(chi2s_over_snr2s, min_instruments),
+		kwarggeniter(phase, min_instruments),
+		kwarggeniter(dt, min_instruments)
+	):
+		return {
+			"segments": segments,
+			"snrs": snrs,
+			"chi2s_over_snr2s": chi2s_over_snr2s,
+			"phase": phase,
+			"dt": dt,
+			"template_id": template_id
+		}
+
+
 class RankingStat(snglcoinc.LnLikelihoodRatioMixin):
 	ligo_lw_name_suffix = u"gstlal_inspiral_rankingstat"
 
@@ -112,12 +135,14 @@ class RankingStat(snglcoinc.LnLikelihoodRatioMixin):
 		self.denominator = inspiral_lr.LnNoiseDensity(template_ids = template_ids, instruments = instruments, delta_t = delta_t, min_instruments = min_instruments)
 		self.zerolag = inspiral_lr.LnLRDensity(template_ids = template_ids, instruments = instruments, delta_t = delta_t, min_instruments = min_instruments)
 
-	def __call__(self, *args, **kwargs):
+	def __call__(self, **kwargs):
 		# fast-path:  network SNR cut
 		if sum(snr**2. for snr in kwargs["snrs"].values()) < self.network_snrsq_threshold:
 			return NegInf
-		# full ln L ranking stat
-		return super(RankingStat, self).__call__(*args, **kwargs)
+		# full ln L ranking stat.  we define the ranking statistic
+		# to be the largest ln L from all allowed subsets of
+		# triggers
+		return max(super(RankingStat, self).__call__(**kwargs) for kwargs in kwarggen(**kwargs))
 
 	@property
 	def template_ids(self):
-- 
GitLab