From 96a2f71dbba5bb317f7b01bb8fe81f11f89b387f Mon Sep 17 00:00:00 2001
From: Kipp Cannon <kipp.cannon@ligo.org>
Date: Wed, 14 Aug 2013 19:11:08 -0400
Subject: [PATCH] far.py: update for addtion of separate pdf arrays to
CoincParamsDistributions
also switch to using stock .filter() via custom filter thread
---
.../bin/gstlal_inspiral_calc_likelihood | 6 +-
.../gstlal_inspiral_create_prior_diststats | 4 +-
.../gstlal_inspiral_followups_from_gracedb | 4 +-
.../bin/gstlal_inspiral_plot_background | 4 +-
.../bin/gstlal_inspiral_plot_likelihoods | 8 +--
...lal_inspiral_plot_likelihoods_from_gracedb | 8 +--
gstlal-inspiral/bin/gstlal_llcbcnode | 8 +--
gstlal-inspiral/python/far.py | 58 +++++++++----------
8 files changed, 48 insertions(+), 52 deletions(-)
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_calc_likelihood b/gstlal-inspiral/bin/gstlal_inspiral_calc_likelihood
index 21119e8994..70af35007e 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_calc_likelihood
+++ b/gstlal-inspiral/bin/gstlal_inspiral_calc_likelihood
@@ -145,11 +145,11 @@ coincparamsdistributions, likelihood_seglists = Far.distribution_stats, Far.live
# calculate injections before writing to disk
if options.synthesize_injections != 0:
- coincparamsdistributions.raw_distributions.add_foreground_prior(n = options.synthesize_injections, instruments = Far.trials_table.get_sngl_ifos(), verbose = options.verbose)
+ coincparamsdistributions.distributions.add_foreground_prior(n = options.synthesize_injections, instruments = Far.trials_table.get_sngl_ifos(), verbose = options.verbose)
# add a uniform prior to background, by default 0 is added so it has no effect
if options.background_prior != 0:
- coincparamsdistributions.raw_distributions.add_background_prior(n = options.background_prior, instruments = Far.trials_table.get_sngl_ifos(), verbose = options.verbose)
+ coincparamsdistributions.distributions.add_background_prior(n = options.background_prior, instruments = Far.trials_table.get_sngl_ifos(), verbose = options.verbose)
if options.verbose:
print >>sys.stderr, "computing event densities ..."
@@ -161,7 +161,7 @@ if options.write_likelihood is not None:
utils.write_filename(xmldoc, options.write_likelihood, gz = (options.write_likelihood or "stdout").endswith(".gz"), verbose = options.verbose)
#FIXME should this come from somewhere else now?
-likelihood_ratio_func = snglcoinc.LikelihoodRatio(coincparamsdistributions.smoothed_distributions)
+likelihood_ratio_func = snglcoinc.LikelihoodRatio(coincparamsdistributions.distributions)
#
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_create_prior_diststats b/gstlal-inspiral/bin/gstlal_inspiral_create_prior_diststats
index 27e604eb1a..6cf63e8edc 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_create_prior_diststats
+++ b/gstlal-inspiral/bin/gstlal_inspiral_create_prior_diststats
@@ -100,9 +100,9 @@ options, filenames = parse_command_line()
coincparamsdistributions = far.DistributionsStats()
if options.background_prior != 0:
- coincparamsdistributions.raw_distributions.add_background_prior(options.background_prior, instruments = options.instrument, verbose = options.verbose)
+ coincparamsdistributions.distributions.add_background_prior(options.background_prior, instruments = options.instrument, verbose = options.verbose)
-coincparamsdistributions.raw_distributions.add_foreground_prior(verbose = options.verbose, n = options.synthesize_injection_count, instruments = options.instrument)
+coincparamsdistributions.distributions.add_foreground_prior(verbose = options.verbose, n = options.synthesize_injection_count, instruments = options.instrument)
# no segments, this file is independent of time
FAR = far.LocalRankingData(segments.segment(None, None), coincparamsdistributions)
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_followups_from_gracedb b/gstlal-inspiral/bin/gstlal_inspiral_followups_from_gracedb
index 777bf97613..b667dabb56 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_followups_from_gracedb
+++ b/gstlal-inspiral/bin/gstlal_inspiral_followups_from_gracedb
@@ -97,8 +97,8 @@ gracedb.writeLog(gid, "SNR vs time", filename = fname, filecontents = open(fname
likelihood_data, process_id = far.LocalRankingData.from_xml(utils.load_url(path, contenthandler = ligolw.LIGOLWContentHandler))
-counts = likelihood_data.distribution_stats.raw_distributions.background_rates
-inj = likelihood_data.distribution_stats.raw_distributions.injection_rates
+counts = likelihood_data.distribution_stats.distributions.background_rates
+inj = likelihood_data.distribution_stats.distributions.injection_rates
bgcol = (224/255.,224/255.,224/255.)
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_plot_background b/gstlal-inspiral/bin/gstlal_inspiral_plot_background
index 32d540959b..310a437da9 100644
--- a/gstlal-inspiral/bin/gstlal_inspiral_plot_background
+++ b/gstlal-inspiral/bin/gstlal_inspiral_plot_background
@@ -49,8 +49,8 @@ total_count = {'H1':0., 'L1':0., 'V1':0.}
for f in files:
Far = far.LocalRankingData.from_filenames([f], verbose = options.verbose)
dists, segs = Far.distribution_stats, Far.livetime_seg
- counts = dists.raw_distributions.background_rates
- inj = dists.raw_distributions.injection_rates
+ counts = dists.distributions.background_rates
+ inj = dists.distributions.injection_rates
likely = copy.deepcopy(inj)
for i, ifo in enumerate(['H1','L1', 'V1']):
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_plot_likelihoods b/gstlal-inspiral/bin/gstlal_inspiral_plot_likelihoods
index 8eec8c2ef0..ab7f6c5484 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_plot_likelihoods
+++ b/gstlal-inspiral/bin/gstlal_inspiral_plot_likelihoods
@@ -88,11 +88,11 @@ print >> sys.stderr, "smooth" in form
if "smooth" in form:
likelihood_data.distribution_stats.finish()
- counts = likelihood_data.distribution_stats.smoothed_distributions.background_rates
- inj = likelihood_data.distribution_stats.smoothed_distributions.injection_rates
+ counts = likelihood_data.distribution_stats.distributions.background_pdf
+ inj = likelihood_data.distribution_stats.distributions.injection_pdf
else:
- counts = likelihood_data.distribution_stats.raw_distributions.background_rates
- inj = likelihood_data.distribution_stats.raw_distributions.injection_rates
+ counts = likelihood_data.distribution_stats.distributions.background_rates
+ inj = likelihood_data.distribution_stats.distributions.injection_rates
bgcol = (224/255.,224/255.,224/255.)
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_plot_likelihoods_from_gracedb b/gstlal-inspiral/bin/gstlal_inspiral_plot_likelihoods_from_gracedb
index fed26ba6fc..47dd8d62b6 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_plot_likelihoods_from_gracedb
+++ b/gstlal-inspiral/bin/gstlal_inspiral_plot_likelihoods_from_gracedb
@@ -112,11 +112,11 @@ print >> sys.stderr, "smooth" in form
if "smooth" in form:
likelihood_data.distribution_stats.finish()
- counts = likelihood_data.distribution_stats.smoothed_distributions.background_rates
- inj = likelihood_data.distribution_stats.smoothed_distributions.injection_rates
+ counts = likelihood_data.distribution_stats.distributions.background_pdf
+ inj = likelihood_data.distribution_stats.distributions.injection_pdf
else:
- counts = likelihood_data.distribution_stats.raw_distributions.background_rates
- inj = likelihood_data.distribution_stats.raw_distributions.injection_rates
+ counts = likelihood_data.distribution_stats.distributions.background_rates
+ inj = likelihood_data.distribution_stats.distributions.injection_rates
bgcol = (224/255.,224/255.,224/255.)
diff --git a/gstlal-inspiral/bin/gstlal_llcbcnode b/gstlal-inspiral/bin/gstlal_llcbcnode
index e9a7584f8b..7a2407bf54 100755
--- a/gstlal-inspiral/bin/gstlal_llcbcnode
+++ b/gstlal-inspiral/bin/gstlal_llcbcnode
@@ -277,11 +277,11 @@ print >> sys.stderr, "smooth" in form
if "smooth" in form:
likelihood_data.distribution_stats.finish()
- counts = likelihood_data.distribution_stats.smoothed_distributions.background_rates
- inj = likelihood_data.distribution_stats.smoothed_distributions.injection_rates
+ counts = likelihood_data.distribution_stats.distributions.background_pdf
+ inj = likelihood_data.distribution_stats.distributions.injection_pdf
else:
- counts = likelihood_data.distribution_stats.raw_distributions.background_rates
- inj = likelihood_data.distribution_stats.raw_distributions.injection_rates
+ counts = likelihood_data.distribution_stats.distributions.background_rates
+ inj = likelihood_data.distribution_stats.distributions.injection_rates
bgcol = (224/255.,224/255.,224/255.)
diff --git a/gstlal-inspiral/python/far.py b/gstlal-inspiral/python/far.py
index 44a92f1ab0..0dc593873f 100644
--- a/gstlal-inspiral/python/far.py
+++ b/gstlal-inspiral/python/far.py
@@ -482,7 +482,7 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
P = 1.0
for name, value in params.items():
if name.endswith("_snr_chi"):
- P *= float(self.background_rates_interp[name](*value))
+ P *= float(self.background_pdf_interp[name](*value))
return P
def P_signal(self, params):
@@ -498,7 +498,7 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
for name, value in params.items():
if name.endswith("_snr_chi"):
- P *= float(self.injection_rates_interp[name](*value))
+ P *= float(self.injection_pdf_interp[name](*value))
return P
def add_background_prior(self, n = 1., transition = 10., instruments = None, prefactors_range = (1.0, 10.0), df = 40, verbose = False):
@@ -530,8 +530,6 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
self.add_foreground_prior(n = n, prefactors_range = prefactors_range, df = df, instruments = instruments, verbose = verbose)
def add_foreground_prior(self, n = 1., prefactors_range = (0.0, 0.10), df = 40, instruments = None, verbose = False):
- # FIXME: for maintainability, this should be modified to
- # use the .add_injection() method of self, but that will slow this down
pfs = numpy.linspace(prefactors_range[0], prefactors_range[1], 10)
for param, binarr in self.injection_rates.items():
new_binarr = rate.BinnedArray(binarr.bins)
@@ -581,7 +579,7 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
# convert signal (aka injection) (rho, chi^2/rho^2) PDFs
# into P(chi^2/rho^2 | rho)
- for name, pdf in self.injection_rates.items():
+ for name, pdf in self.injection_pdf.items():
if not name.endswith("_snr_chi"):
continue
bin_sizes = pdf.bins[1].upper() - pdf.bins[1].lower()
@@ -589,7 +587,7 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
nonzero = pdf.array[i] != 0
pdf.array[i] /= numpy.dot(numpy.compress(nonzero, pdf.array[i]), numpy.compress(nonzero, bin_sizes))
# rebuild the interpolator
- self.injection_rates_interp[name] = rate.InterpBinnedArray(pdf)
+ self.injection_pdf_interp[name] = rate.InterpBinnedArray(pdf)
@classmethod
def from_xml(cls, xml, name):
@@ -622,52 +620,51 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
#
+class FilterThread(snglcoinc.CoincParamsFilterThread):
+ # populate the bins with probabilities (not probability densities
+ # as in the default implementation)
+ def run(self):
+ with self.cpu:
+ if self.verbose:
+ with self.stderr:
+ print >>sys.stderr, "%s," % self.getName(),
+ rate.filter_array(self.binnedarray.array, self.filter)
+ self.binnedarray.array /= self.binnedarray.array.sum()
+
+
class DistributionsStats(object):
"""
A class used to populate a ThincaCoincParamsDistribution instance using
event parameter data.
"""
def __init__(self):
- self.raw_distributions = ThincaCoincParamsDistributions()
- self.smoothed_distributions = ThincaCoincParamsDistributions()
+ self.distributions = ThincaCoincParamsDistributions()
self.likelihood_pdfs = {}
self.target_length = 1000
def __add__(self, other):
out = type(self)()
- out.raw_distributions += self.raw_distributions
- out.raw_distributions += other.raw_distributions
- #FIXME do we also add the smoothed distributions??
+ out.distributions += self.distributions
+ out.distributions += other.distributions
return out
def add_single(self, event):
- self.raw_distributions.add_background(self.raw_distributions.coinc_params((event,), None))
+ self.distributions.add_background(self.distributions.coinc_params((event,), None))
def finish(self, verbose = False):
- self.smoothed_distributions = self.raw_distributions.copy()
- #self.smoothed_distributions.finish(verbose = verbose)
- # FIXME: should be the line above, we'll temporarily do
- # the following. the difference is that the above produces
- # PDFs while what follows produces probabilities in each
- # bin
if verbose:
print >>sys.stderr, "smoothing parameter distributions ...",
- for name, binnedarray in itertools.chain(self.smoothed_distributions.background_rates.items(), self.smoothed_distributions.injection_rates.items()):
- if verbose:
- print >>sys.stderr, "%s," % name,
- rate.filter_array(binnedarray.array, self.raw_distributions.filters[name])
- numpy.clip(binnedarray.array, 0.0, PosInf, binnedarray.array)
- binnedarray.array /= binnedarray.array.sum()
+ self.distributions.finish(verbose = verbose, filterthread = FilterThread)
if verbose:
print >>sys.stderr, "done"
def compute_single_instrument_background(self, instruments = None, verbose = False):
# initialize a likelihood ratio evaluator
- likelihood_ratio_evaluator = snglcoinc.LikelihoodRatio(self.smoothed_distributions)
+ likelihood_ratio_evaluator = snglcoinc.LikelihoodRatio(self.distributions)
# reduce typing
- background = self.smoothed_distributions.background_rates
- injections = self.smoothed_distributions.injection_rates
+ background = self.distributions.background_pdf
+ injections = self.distributions.injection_pdf
self.likelihood_pdfs.clear()
for param in background:
@@ -700,13 +697,12 @@ class DistributionsStats(object):
@classmethod
def from_xml(cls, xml, name):
self = cls()
- # FIXME: produce error if raw_distributions' binnings don't match the arrays in the file?
- self.raw_distributions, process_id = ThincaCoincParamsDistributions.from_xml(xml, name)
- self.smoothed_distributions = ThincaCoincParamsDistributions()
+ # FIXME: produce error if distributions' binnings don't match the arrays in the file?
+ self.distributions, process_id = ThincaCoincParamsDistributions.from_xml(xml, name)
return self, process_id
def to_xml(self, process, name):
- return self.raw_distributions.to_xml(process, name)
+ return self.distributions.to_xml(process, name)
#
--
GitLab