bin specific Chads extinction before marginalization

gstlal-inspiral/bin/gstlal_inspiral_workflow

@@ -131,9 +131,10 @@ elif args.command == "create":
 				triggers += inj_triggers
 			triggers = dag.add_and_simplify(triggers)
 			triggers = dag.calc_likelihood(triggers, dist_stats)
-			triggers = dag.merge_and_reduce(triggers)
 
 			pdfs = dag.calc_pdf(dist_stats)
+			pdfs = dag.extinct_bin(pdfs, triggers)
+			triggers = dag.merge_and_reduce(triggers)
 			pdfs = dag.marginalize_pdf(pdfs)
 
 			if config.filter.injections:

gstlal-inspiral/python/dags/layers/inspiral.py

@@ -2500,6 +2500,51 @@ def add_and_simplify_layer(config, dag, trigger_cache):
 	return remaining_cache
 
 
+def extinct_bin_layer(config, dag, pdf_cache, trigger_cache):
+	extinct_layer = Layer(
+		"gstlal_inspiral_extinct_bin",
+		name="gstlal_inspiral_extinct_bin",
+		requirements={
+			"request_cpus": 1,
+			"request_memory": 2000,
+			"request_disk": "1GB",
+			**config.condor.submit
+		},
+		transfer_files=config.condor.transfer_files,
+	)
+
+	output_cache = []
+	trigger_cache = trigger_cache.groupby('subtype')[""] # noninj triggers
+	trigger_cache = trigger_cache.groupby('bin')
+	for svd_bin, pdfs in pdf_cache.groupby("bin").items():
+		assert len(trigger_cache[svd_bin].files) == 1, "Must provide exactly one trigger file per bin to the extinct_bin layer"
+		trigger_file = trigger_cache[svd_bin].files[0]
+		output_filename = pdfs.name.filename(
+			config.ifos,
+			span=config.span,
+			svd_bin=f"{svd_bin}",
+		)
+		output_path = os.path.join(
+			pdfs.name.directory(root="rank", start = config.span[0]),
+			output_filename,
+		)
+		output_cache.append(output_path)
+
+		extinct_layer += Node(
+			inputs = [Option("input-trigger-file", trigger_file), Argument("input-pdfs", pdfs.files)],
+			outputs = Option("output", output_path),
+		)
+
+	dag.attach(extinct_layer)
+
+	pdf_cache = DataCache.from_files(
+		DataType.DIST_STAT_PDFS,
+		output_cache
+		)
+
+	return pdf_cache
+
+
 def merge_and_reduce_layer(config, dag, trigger_cache):
 	# FIXME: find better way of discovering SQL file
 	share_path = os.path.split(dagutil.which("gstlal_inspiral"))[0].replace("bin", "share/gstlal")
@@ -2788,4 +2833,5 @@ def layers():
 		"add_and_simplify": add_and_simplify_layer,
 		"merge_and_reduce": merge_and_reduce_layer,
 		"find_injections_lite": find_injections_lite_layer,
+		"extinct_bin": extinct_bin_layer,
 	}

gstlal-ugly/bin/Makefile.am

@@ -22,6 +22,7 @@ dist_bin_SCRIPTS = \
 	gstlal_inspiral_treebank \
 	gstlal_inspiral_treebank_dag \
 	gstlal_inspiral_bankviz \
+	gstlal_inspiral_extinct_bin \
 	gstlal_segments_trim \
 	gstlal_glitch_population \
 	gstlal_ifo_stat \

gstlal-ugly/bin/gstlal_inspiral_extinct_bin

+#!/usr/bin/env python3
+  
+from optparse import OptionParser
+from gstlal import far
+from ligo.lw import utils as ligolw_utils
+from ligo.lw import ligolw
+from ligo.lw import dbtables, lsctables
+from ligo.lw.utils import ligolw_sqlite
+from lalinspiral import thinca
+import os
+import sqlite3
+
+def parse_command_line():
+	parser = OptionParser()
+	parser.add_option("--input-trigger-file", help = "input trigger file from which to add zerolag")
+	parser.add_option("--output", help = "output")
+
+	options, urls = parser.parse_args()
+	return options, urls
+
+options, urls = parse_command_line()
+
+# load and marginalize the rankingstatpdf
+rankingstatpdf = far.marginalize_pdf_urls(urls, which = "RankingStatPDF")
+
+# read the trigeger xml file as a sqlite file
+with sqlite3.connect(":memory:") as connection:
+	@dbtables.use_in
+	class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
+		connection = connection
+	ligolw_sqlite.insert_from_url(options.input_trigger_file, contenthandler = LIGOLWContentHandler, preserve_ids = True)
+	del LIGOLWContentHandler
+
+	# add zerolag from trigger file to rankingstatpdf
+	xmldoc = dbtables.get_xml(connection)
+	coinc_def_id = lsctables.CoincDefTable.get_table(xmldoc).get_coinc_def_id(thinca.InspiralCoincDef.search, thinca.InspiralCoincDef.search_coinc_type, create_new = False)
+	xmldoc.unlink()
+	rankingstatpdf.collect_zero_lag_rates(connection, coinc_def_id)
+
+lrs = rankingstatpdf.noise_lr_lnpdf.centres()[0]
+bg = rankingstatpdf.noise_lr_lnpdf.array
+fg = rankingstatpdf.zero_lag_lr_lnpdf.array
+
+# zero out the beginning bins of each because they are notoriously bad and should just be ignored IMO
+bg[:10] = 0.
+fg[:10] = 0.
+
+# solve for the first nonzero bin after that
+#ix = (fg > 200).argmax()
+ix = (lrs > -20).argmax()
+
+print ("setting counts to zero below ix", ix, "at lr", lrs[ix])
+bg[:ix] = 0.
+fg[:ix] = 0.
+
+fg_ccdf = 1.0 - fg.cumsum() / fg.sum()
+bg_ccdf = 1.0 - bg.cumsum() / bg.sum()
+ 
+# Calculate an extincted PDF counts and a new ccdf from that
+survival_probability = 1.0 - fg_ccdf
+bg_extinct = bg * survival_probability
+bg_extinct /= bg_extinct.sum()
+bg_extinct *= fg.sum()
+#bg_ccdf_extinct = 1.0 - bg_extinct.cumsum() / bg_extinct.sum()
+
+rankingstatpdf.noise_lr_lnpdf.array = bg_extinct
+rankingstatpdf.noise_lr_lnpdf.normalize()
+rankingstatpdf.zero_lag_lr_lnpdf.array = fg
+rankingstatpdf.zero_lag_lr_lnpdf.normalize()
+
+
+xmldoc = ligolw.Document()
+xmldoc.appendChild(ligolw.LIGO_LW())
+far.gen_likelihood_control_doc(xmldoc, None, rankingstatpdf)
+ligolw_utils.write_filename(xmldoc, options.output)