gstlal_inspiral, gstlal_inspiral_create_prior_diststats, far.py,...

gstlal_inspiral, gstlal_inspiral_create_prior_diststats, far.py, inspiral_lr.py, svd_bank.py: add support for per template corrections to the horizon distance used in the likelihood ratio using the sigma values computed during the SVD decomposition

gstlal-inspiral/bin/gstlal_inspiral

@@ -691,8 +691,10 @@ for output_file_number, (svd_bank_url_dict, output_url, ranking_stat_output_url,
 	# assume all instruments have the same templates, just extract them
 	# from one of the instruments at random
 	sngl_inspiral_table = banks.values()[0][0].sngl_inspiral_table.copy()
+	horizon_factors = {}
 	for bank in banks.values()[0]:
 		sngl_inspiral_table.extend(bank.sngl_inspiral_table)
+		horizon_factors.update(bank.horizon_factors)
 	@bottle.route("/template_bank.xml.gz")
 	def get_template_bank_xml(sngl_inspiral_table = sngl_inspiral_table):
 		xmldoc = ligolw.Document()
@@ -704,6 +706,7 @@ for output_file_number, (svd_bank_url_dict, output_url, ranking_stat_output_url,
 		output.close()
 		return outstr
 	template_ids = frozenset(row.template_id for row in sngl_inspiral_table)
+	assert set(horizon_factors) == template_ids, "horizon factors are not assigned for each template id"
 	assert len(template_ids) == len(sngl_inspiral_table), "template IDs are not unique within the template bank"
 	@bottle.route("/template_ids.txt")
 	def get_template_ids(template_ids = sorted(template_ids)):
@@ -783,8 +786,9 @@ for output_file_number, (svd_bank_url_dict, output_url, ranking_stat_output_url,
 				rankingstat.template_ids = template_ids
 			elif rankingstat.template_ids != template_ids:
 				raise ValueError("\"%s\" is for the wrong templates" % options.ranking_stat_input)
+			rankingstat.numerator.set_horizon_factors(horizon_factors)
 	if rankingstat is None:
-		rankingstat = far.RankingStat(template_ids = template_ids, instruments = all_instruments, delta_t = options.coincidence_threshold, min_instruments = options.min_instruments)
+		rankingstat = far.RankingStat(template_ids = template_ids, instruments = all_instruments, delta_t = options.coincidence_threshold, min_instruments = options.min_instruments, horizon_factors = horizon_factors)
 		rankingstat.numerator.add_signal_model(df = 40)
 
 
@@ -811,7 +815,7 @@ for output_file_number, (svd_bank_url_dict, output_url, ranking_stat_output_url,
 			verbose = options.verbose
 		),
 		rankingstat = rankingstat,
-		horizon_distance_func = svd_bank.make_horizon_distance_func(banks),
+		horizon_distance_func = banks.values()[0][0].horizon_distance_func,# they should all be the same
 		gracedbwrapper = inspiral.GracedBWrapper(
 			instruments = rankingstat.instruments,
 			far_threshold = options.gracedb_far_threshold,

gstlal-inspiral/bin/gstlal_inspiral_create_prior_diststats

@@ -94,10 +94,12 @@ def parse_command_line():
 		raise ValueError("unrecognized arguments after options: %s" % " ".join(filenames))
 
 	template_ids = []
+	horizon_factors = {}
 	for n, bank in enumerate(svd_bank.read_banks(options.svd_file, contenthandler = LIGOLWContentHandler, verbose = options.verbose)):
 		template_ids += [row.template_id for row in bank.sngl_inspiral_table]
+		horizon_factors.update(bank.horizon_factors)
 
-	return options, process_params, filenames, template_ids
+	return options, process_params, filenames, template_ids, horizon_factors
 
 
 #
@@ -114,7 +116,7 @@ def parse_command_line():
 #
 
 
-options, process_params, filenames, template_ids = parse_command_line()
+options, process_params, filenames, template_ids, horizon_factors = parse_command_line()
 
 
 #
@@ -132,7 +134,7 @@ process = ligolw_process.register_to_xmldoc(xmldoc, u"gstlal_inspiral_create_pri
 #
 
 
-rankingstat = far.RankingStat(template_ids = template_ids, instruments = options.instrument, delta_t = options.coincidence_threshold, min_instruments = options.min_instruments, population_model_file = options.mass_model_file)
+rankingstat = far.RankingStat(template_ids = template_ids, instruments = options.instrument, delta_t = options.coincidence_threshold, min_instruments = options.min_instruments, population_model_file = options.mass_model_file, horizon_factors = horizon_factors)
 
 if options.background_prior > 0:
 	rankingstat.denominator.add_noise_model(number_of_events = options.background_prior)

gstlal-inspiral/python/far.py

@@ -129,10 +129,10 @@ class RankingStat(snglcoinc.LnLikelihoodRatioMixin):
 		pass
 
 	# network SNR threshold
-	network_snrsq_threshold = 6.0**2.
+	network_snrsq_threshold = 5.0**2.
 
-	def __init__(self, template_ids = None, instruments = frozenset(("H1", "L1", "V1")), population_model_file = None, min_instruments = 1, delta_t = 0.005):
-		self.numerator = inspiral_lr.LnSignalDensity(template_ids = template_ids, instruments = instruments, delta_t = delta_t, population_model_file = population_model_file, min_instruments = min_instruments)
+	def __init__(self, template_ids = None, instruments = frozenset(("H1", "L1", "V1")), population_model_file = None, min_instruments = 1, delta_t = 0.005, horizon_factors = None):
+		self.numerator = inspiral_lr.LnSignalDensity(template_ids = template_ids, instruments = instruments, delta_t = delta_t, population_model_file = population_model_file, min_instruments = min_instruments, horizon_factors = horizon_factors)
 		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)
 

gstlal-inspiral/python/stats/inspiral_lr.py

@@ -39,6 +39,7 @@ from scipy import interpolate
 from scipy import stats
 import sys
 import warnings
+import json
 
 
 from glue.ligolw import ligolw
@@ -292,6 +293,7 @@ class LnLRDensity(snglcoinc.LnLRDensity):
 class LnSignalDensity(LnLRDensity):
 	def __init__(self, *args, **kwargs):
 		population_model_file = kwargs.pop("population_model_file", None)
+		self.horizon_factors = kwargs.pop("horizon_factors", None)
 		super(LnSignalDensity, self).__init__(*args, **kwargs)
 
 		# install SNR, chi^2 PDF (one for all instruments)
@@ -309,6 +311,9 @@ class LnSignalDensity(LnLRDensity):
 		self.population_model = inspiral_intrinsics.SourcePopulationModel(self.template_ids, filename = self.population_model_file)
 		self.InspiralExtrinsics = inspiral_extrinsics.InspiralExtrinsics(self.min_instruments)
 
+	def set_horizon_factors(self, horizon_factors):
+		self.horizon_factors = horizon_factors
+
 	def __call__(self, segments, snrs, chi2s_over_snr2s, phase, dt, template_id):
 		assert frozenset(segments) == self.instruments
 		if len(snrs) < self.min_instruments:
@@ -341,7 +346,12 @@ class LnSignalDensity(LnLRDensity):
 		horizon = sorted(horizons.values())[-self.min_instruments] / TYPICAL_HORIZON_DISTANCE
 		if not horizon:
 			return NegInf
-		lnP = 3. * math.log(horizon) + math.log(len(self.template_ids))
+		# horizon factors adjusts the computed horizon factor according
+		# to the sigma values computed at the time of the SVD. This
+		# gives a good approximation to the horizon distance for each
+		# template without computing them each explicitly. Only one
+		# template has its horizon calculated explicitly.
+		lnP = 3. * math.log(horizon * self.horizon_factors[template_id]) + math.log(len(self.template_ids))
 
 		# Add P(instruments | horizon distances)
 		try:
@@ -378,6 +388,10 @@ class LnSignalDensity(LnLRDensity):
 			raise ValueError("incompatible mass model file names")
 		if self.population_model_file is None and other.population_model_file is not None:
 			self.population_model_file = other.population_model_file
+		if self.horizon_factors is not None and other.horizon_factors is not None and other.horizon_factors != self.horizon_factors:
+			raise ValueError("incompatible horizon_factors")
+		if self.horizon_factors is None and other.horizon_factors is not None:
+			self.horizon_factors = other.horizon_factors
 
 		return self
 
@@ -391,6 +405,7 @@ class LnSignalDensity(LnLRDensity):
 		# okay to use references because read-only data
 		new.population_model = self.population_model
 		new.InspiralExtrinsics = self.InspiralExtrinsics
+		new.horizon_factors = self.horizon_factors
 		return new
 
 	def local_mean_horizon_distance(self, gps, window = segments.segment(-32., +2.)):
@@ -534,6 +549,7 @@ class LnSignalDensity(LnLRDensity):
 		xml = super(LnSignalDensity, self).to_xml(name)
 		xml.appendChild(self.horizon_history.to_xml(u"horizon_history"))
 		xml.appendChild(ligolw_param.Param.from_pyvalue(u"population_model_file", self.population_model_file))
+		xml.appendChild(ligolw_param.Param.from_pyvalue(u"horizon_factors", json.dumps(self.horizon_factors) if self.horizon_factors is not None else None))
 		return xml
 
 	@classmethod
@@ -542,6 +558,13 @@ class LnSignalDensity(LnLRDensity):
 		self = super(LnSignalDensity, cls).from_xml(xml, name)
 		self.horizon_history = horizonhistory.HorizonHistories.from_xml(xml, u"horizon_history")
 		self.population_model_file = ligolw_param.get_pyvalue(xml, u"population_model_file")
+		self.horizon_factors = ligolw_param.get_pyvalue(xml, u"horizon_factors")
+		if self.horizon_factors is not None:
+			# FIXME, how do we properly decode the json, I assume something in ligolw can do this?
+			self.horizon_factors = self.horizon_factors.replace("\\","").replace('\\"','"')
+			self.horizon_factors = json.loads(self.horizon_factors)
+			self.horizon_factors = dict((int(k), v) for k, v in self.horizon_factors.items())
+			assert set(self.template_ids) == set(self.horizon_factors)
 		self.population_model = inspiral_intrinsics.SourcePopulationModel(self.template_ids, filename = self.population_model_file)
 		self.InspiralExtrinsics = inspiral_extrinsics.InspiralExtrinsics(self.min_instruments)
 		return self

gstlal-inspiral/python/svd_bank.py

@@ -388,6 +388,9 @@ def read_banks(filename, contenthandler, verbose = False):
 		bank.autocorrelation_mask = ligolw_array.get_array(root, 'autocorrelation_mask').array
 		bank.sigmasq = ligolw_array.get_array(root, 'sigmasq').array
 
+		# prepare the horizon distance factors
+		bank.horizon_factors = dict((row.template_id, sigmasq**.5) for row, sigmasq in zip(bank.sngl_inspiral_table, bank.sigmasq))
+
 		# attach a reference to the psd
 		bank.processed_psd = processed_psd
 
@@ -418,6 +421,15 @@ def read_banks(filename, contenthandler, verbose = False):
 			bank.bank_fragments.append(frag)
 
 		banks.append(bank)
+	template_id, func = horizon_distance_func(banks)
+	horizon_norm = None
+	for bank in banks:
+		if template_id in bank.horizon_factors:
+			assert horizon_norm is None
+			horizon_norm = bank.horizon_factors[template_id]
+	for bank in banks:
+		bank.horizon_distance_func = func
+		bank.horizon_factors = dict((tid, f / horizon_norm) for (tid, f) in bank.horizon_factors.items())
 	xmldoc.unlink()
 	return banks
 
@@ -449,27 +461,21 @@ def svdbank_templates_mapping(filenames, contenthandler, verbose = False):
 	return mapping
 
 
-def make_horizon_distance_func(banks):
+def horizon_distance_func(banks):
 	"""
 	Takes a dictionary of objects returned by read_banks keyed by instrument
 	"""
 	# span is [15 Hz, 0.85 * Nyquist frequency]
 	# find the Nyquist frequency for the PSD to be used for each
 	# instrument.  require them to all match
-	sngl_inspiral_table = banks.values()[0][0].sngl_inspiral_table.copy()
-	for bank in banks.values()[0]:
-		sngl_inspiral_table.extend(bank.sngl_inspiral_table)
-	nyquists = set(max(rate for bank in banklist for rate in bank.get_rates()) // 2 for instrument, banklist in banks.items())
+	nyquists = set((max(bank.get_rates()) for bank in banks))
 	assert len(nyquists) == 1, "all banks must have the same Nyquist frequency to define a consistent horizon distance function (got %s)" % ", ".join("%g" % rate for rate in sorted(nyquists))
-	# assume default 32 s PSD.  this is not required to be correct, but
+	# assume default 4 s PSD.  this is not required to be correct, but
 	# for best accuracy it should not be larger than the true value and
 	# for best performance it should not be smaller than the true
 	# value.
-	deltaF = 1. / 32.
-	# FIXME (from Chad) What is the 5/3 for???
-	# pick (m1, m2) from the median template ranked by Mchirp^(5/3)
-	# to provide the canonical waveform model.  See Maggiore equation
-	# (4.3).
-	assert len(sngl_inspiral_table) > 0, "no templates:  must have templates to define horizon distance function"
-	median_row = sorted(sngl_inspiral_table, key = lambda row: row.mchirp**(5./3.))[len(sngl_inspiral_table) // 2]
-	return reference_psd.HorizonDistance(15.0, 0.85 * max(nyquists), deltaF, median_row.mass1, median_row.mass2)
+	deltaF = 1. / 4.
+	# use the minimum template id as the cannonical horizon function
+	template_id, m1, m2 = min((row.template_id, row.mass1, row.mass2) for bank in banks for row in bank.sngl_inspiral_table)
+
+	return template_id, reference_psd.HorizonDistance(15.0, 0.85 * max(nyquists), deltaF, m1, m2)