inspiral_extrinsics: add NumeratorSNRCHIPDF

- subclass of rate.BinnedLnPDF to implement 1-parameter family of 1-D PDFs P(\chi^2 | snr, signal)

gstlal-inspiral/python/stats/inspiral_extrinsics.py

 # Copyright (C) 2016,2017  Kipp Cannon
+# Copyright (C) 2011--2014  Kipp Cannon, Chad Hanna, Drew Keppel
+# Copyright (C) 2013  Jacob Peoples
 #
 # This program is free software; you can redistribute it and/or modify it
 # under the terms of the GNU General Public License as published by the
@@ -47,6 +49,7 @@ from glue.ligolw import array as ligolw_array
 from glue.ligolw import param as ligolw_param
 from glue.ligolw import utils as ligolw_utils
 from glue.text_progress_bar import ProgressBar
+from gstlal import stats as gstlalstats
 import lal
 from lal import rate
 from lalburst import snglcoinc
@@ -63,7 +66,8 @@ from gstlal import paths as gstlal_config_paths
 __all__ = [
 	"instruments_rate_given_noise",
 	"P_instruments_given_signal",
-	"SNRPDF"
+	"SNRPDF",
+	"NumeratorSNRCHIPDF"
 ]
 
 
@@ -915,3 +919,100 @@ def lnP_dt_dphi(params, coincidence_window_extension, model = "noise"):
 
 	else:
 		raise ValueError("invalid data model '%s'" % mode)
+
+
+#
+# =============================================================================
+#
+#                               SNR, \chi^2 PDF
+#
+# =============================================================================
+#
+
+
+class NumeratorSNRCHIPDF(rate.BinnedLnPDF):
+	"""
+	Reports ln P(chi^2/rho^2 | rho, signal)
+	"""
+	# NOTE:  by happy coincidence numpy's broadcasting rules allow the
+	# inherited .at_centres() method to be used as-is
+	def __init__(self, *args, **kwargs):
+		super(rate.BinnedLnPDF, self).__init__(*args, **kwargs)
+		self.volume = self.bins[1].upper() - self.bins[1].lower()
+		# FIXME: instead of .shape and repeat() use .broadcast_to()
+		# when we can rely on numpy >= 1.8
+		#self.volume = numpy.broadcast_to(self.volume, self.bins.shape)
+		self.volume.shape = (1, len(self.volume))
+		self.volume = numpy.repeat(self.volume, len(self.bins[0]), 0)
+		self.norm = numpy.zeros((len(self.bins[0]), 1))
+
+	def __getitem__(self, coords):
+		return numpy.log(super(rate.BinnedLnPDF, self).__getitem__(coords)) - self.norm[self.bins(*coords)[0]]
+
+	def marginalize(self, *args, **kwargs):
+		raise NotImplementedError
+
+	def __iadd__(self, other):
+		super(rate.BinnedLnPDF, self).__iadd__(other)
+		self.norm = numpy.log(numpy.exp(self.norm) + numpy.exp(other.norm))
+		return self
+
+	def normalize(self):
+		# replace the vector with a new one so that we don't
+		# interfere with any copies that might have been made
+		with numpy.errstate(divide = "ignore"):
+			self.norm = numpy.log(self.array.sum(axis = 1))
+		self.norm.shape = (len(self.norm), 1)
+
+	@staticmethod
+	def add_signal_model(lnpdf, n, prefactors_range, df, inv_snr_pow = 4., snr_min = 4., progressbar = None):
+		if df <= 0.:
+			raise ValueError("require df >= 0: %s" % repr(df))
+		pfs = numpy.linspace(prefactors_range[0], prefactors_range[1], 100)
+		if progressbar is not None:
+			progressbar.max = len(pfs)
+
+		# FIXME:  except for the low-SNR cut, the slicing is done
+		# to work around various overflow and loss-of-precision
+		# issues in the extreme parts of the domain of definition.
+		# it would be nice to identify the causes of these and
+		# either fix them or ignore them one-by-one with a comment
+		# explaining why it's OK to ignore the ones being ignored.
+		# for example, computing snrchi2 by exponentiating the sum
+		# of the logs of the terms might permit its evaluation
+		# everywhere on the domain.  can ncx2pdf() be made to work
+		# everywhere?
+		snrindices, rcossindices = lnpdf.bins[snr_min:1e10, 1e-10:1e10]
+		snr, dsnr = lnpdf.bins[0].centres()[snrindices], lnpdf.bins[0].upper()[snrindices] - lnpdf.bins[0].lower()[snrindices]
+		rcoss, drcoss = lnpdf.bins[1].centres()[rcossindices], lnpdf.bins[1].upper()[rcossindices] - lnpdf.bins[1].lower()[rcossindices]
+
+		snr2 = snr**2.
+		snrchi2 = numpy.outer(snr2, rcoss) * df
+
+		arr = numpy.zeros_like(lnpdf.array)
+		for pf in pfs:
+			if progressbar is not None:
+				progressbar.increment()
+			arr[snrindices, rcossindices] += gstlalstats.ncx2pdf(snrchi2, df, numpy.array([pf * snr2]).T)
+
+		# convert to counts by multiplying by bin volume, and also
+		# multiply by an SNR powr law
+		arr[snrindices, rcossindices] *= numpy.outer(dsnr / snr**inv_snr_pow, drcoss)
+
+		# normalize to a total count of n
+		arr *= n / arr.sum()
+
+		# add to lnpdf
+		lnpdf.array += arr
+
+	def to_xml(self, *args, **kwargs):
+		elem = super(rate.BinnedLnPDF, self).to_xml(*args, **kwargs)
+		elem.appendChild(ligolw_array.Array.build("norm", self.norm))
+		return elem
+
+	@classmethod
+	def from_xml(cls, xml, name):
+		xml = cls.get_xml_root(xml, name)
+		self = super(rate.BinnedLnPDF, cls).from_xml(xml, name)
+		self.norm = ligolw_array.get_array(xml, "norm").array
+		return self