diff --git a/gstlal-inspiral/bin/gstlal_inspiral b/gstlal-inspiral/bin/gstlal_inspiral
index 29ef3ad39bfee2884f87620fbc5d37d4b990ee92..7c70c1b8f169ebb96ad6befedbb61255d92ad7b8 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral
+++ b/gstlal-inspiral/bin/gstlal_inspiral
@@ -684,7 +684,7 @@ for output_file_number, (svd_bank, output_filename, likelihood_namedtuple, zero_
if coinc_params_distributions is None:
raise ValueError("\"%s\" does not contain parameter distribution data" % likelihood_namedtuple[0] if likelihood_namedtuple[0] is not None else likelihood_namedtuple[1])
else:
- coinc_params_distributions, seglists = far.ThincaCoincParamsDistributions(), segments.segmentlistdict((instrument, segments.segmentlist()) for instrument in detectors.channel_dict)
+ coinc_params_distributions, seglists = far.ThincaCoincParamsDistributions(instruments = set(detectors.channel_dict)), segments.segmentlistdict((instrument, segments.segmentlist()) for instrument in detectors.channel_dict)
if zero_lag_ranking_stats_file is None:
zero_lag_ranking_stats = None
@@ -692,6 +692,8 @@ for output_file_number, (svd_bank, output_filename, likelihood_namedtuple, zero_
_, zero_lag_ranking_stats, _ = far.parse_likelihood_control_doc(ligolw_utils.load_filename(zero_lag_ranking_stats_file, verbose = options.verbose, contenthandler = far.ThincaCoincParamsDistributions.LIGOLWContentHandler))
if zero_lag_ranking_stats is None:
raise ValueError("\"%s\" does not contain ranking statistic PDF data" % zero_lag_ranking_stat_file)
+ if zero_lag_ranking_stats.instruments != set(detectors.channel_dict):
+ raise ValueError("\"%s\" is for the instruments %s, we need %s" % (", ".join(sorted(zero_lag_ranking_stats.instruments)), ", ".join(sorted(detectors.channel_dict))))
#
@@ -705,7 +707,6 @@ for output_file_number, (svd_bank, output_filename, likelihood_namedtuple, zero_
filename = output_filename or "%s-%s_LLOID-%d-%d.xml.gz" % (lsctables.ifos_from_instrument_set(detectors.channel_dict.keys()).replace(",", ""), options.job_tag, int(detectors.seg[0]), int(abs(detectors.seg))),
process_params = process_params,
pipeline = pipeline,
- instruments = set(detectors.channel_dict),
seg = detectors.seg or segments.segment(LIGOTimeGPS(0), LIGOTimeGPS(2000000000)), # online data doesn't have a segment so make it all possible time
coincidence_threshold = options.coincidence_threshold,
coinc_params_distributions = coinc_params_distributions,
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_calc_likelihood b/gstlal-inspiral/bin/gstlal_inspiral_calc_likelihood
index 264ad8c70954f0be2a874fcc20ce7021a01e0b4f..a2ce0db137abf16fef852b5b087b158feacb98bb 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_calc_likelihood
+++ b/gstlal-inspiral/bin/gstlal_inspiral_calc_likelihood
@@ -171,7 +171,7 @@ ln_likelihood_ratio_func = snglcoinc.LnLikelihoodRatio(coincparamsdistributions)
if options.write_rates is not None:
- ranking_data = far.RankingData(None, seglists.keys())
+ ranking_data = far.RankingData(None, instruments = seglists.keys())
else:
ranking_data = None
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_calc_rank_pdfs b/gstlal-inspiral/bin/gstlal_inspiral_calc_rank_pdfs
index a4224e1a71aa904f349b7b0efb470429ac133b3a..225afa2dc68ff45042d061c69bc8675118457e94 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_calc_rank_pdfs
+++ b/gstlal-inspiral/bin/gstlal_inspiral_calc_rank_pdfs
@@ -140,7 +140,7 @@ coincparamsdistributions.finish(verbose = options.verbose)
#
-ranking_data = far.RankingData(coincparamsdistributions, instruments = seglists.keys(), sampler_coinc_params_distributions = sampler_coincparamsdistributions, nsamples = options.ranking_stat_samples, verbose = options.verbose)
+ranking_data = far.RankingData(coincparamsdistributions, sampler_coinc_params_distributions = sampler_coincparamsdistributions, nsamples = options.ranking_stat_samples, verbose = options.verbose)
#
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_create_prior_diststats b/gstlal-inspiral/bin/gstlal_inspiral_create_prior_diststats
index c22f76fbc817a5033a144a6d58637824c2e539b9..5be7c635588d7bcbe95b4a92834b7644f5bc574a 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_create_prior_diststats
+++ b/gstlal-inspiral/bin/gstlal_inspiral_create_prior_diststats
@@ -117,7 +117,7 @@ process = ligolw_process.register_to_xmldoc(xmldoc, u"gstlal_inspiral_create_pri
#
-coincparamsdistributions = far.ThincaCoincParamsDistributions()
+coincparamsdistributions = far.ThincaCoincParamsDistributions(instruments = options.instrument)
if options.background_prior > 0:
coincparamsdistributions.add_background_prior(n = dict((ifo, options.background_prior) for ifo in options.instrument), verbose = options.verbose)
diff --git a/gstlal-inspiral/bin/gstlal_ll_inspiral_create_prior_diststats b/gstlal-inspiral/bin/gstlal_ll_inspiral_create_prior_diststats
index c7b292cdebab5f7040236f75e5ccd990615c84c2..d2e523c27e0d240f8e772b1d189a7c641a2958e4 100755
--- a/gstlal-inspiral/bin/gstlal_ll_inspiral_create_prior_diststats
+++ b/gstlal-inspiral/bin/gstlal_ll_inspiral_create_prior_diststats
@@ -133,7 +133,7 @@ n_zerolag = dict(((ifo, float(abs(seg)) * options.num_templates / 50.) for ifo,
# Initialize an empty ThincaCoincParamsDistributions class
#
-diststats = far.ThincaCoincParamsDistributions()
+diststats = far.ThincaCoincParamsDistributions(instruments = set(segs))
#
# Add background, zero_lag and injection prior distributions in the SNR and chi-squared plane
@@ -179,7 +179,7 @@ process = ligolw_process.register_to_xmldoc(xmldoc, u"gstlal_ll_inspiral_create_
# Instantiate a new RankingData class from our distribution stats
#
-ranking_data = far.RankingData(diststats, instruments = options.instruments, process_id = process.process_id, nsamples = 100000, verbose = options.verbose)
+ranking_data = far.RankingData(diststats, process_id = process.process_id, nsamples = 100000, verbose = options.verbose)
#
# Simulate a measured coincident trigger likelihood histogram by just using the background one with a lower normalized count
diff --git a/gstlal-inspiral/python/far.py b/gstlal-inspiral/python/far.py
index b49295843ae294d8a6fdae7009b1776f53be9a16..e30f2b908ab9b335e73365d15080bc26c4c62c88 100644
--- a/gstlal-inspiral/python/far.py
+++ b/gstlal-inspiral/python/far.py
@@ -119,47 +119,62 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
# accidental_weight * (measured denominator)
numerator_accidental_weight = 0.
- # binnings (filter funcs look-up initialized in .__init__()
- snr_chi_binning = rate.NDBins((rate.ATanLogarithmicBins(3.6, 70., 600), rate.ATanLogarithmicBins(.001, 0.5, 300)))
+ # binnings (initialized in .__init__()
binnings = {
- "instruments": rate.NDBins((snglcoinc.InstrumentBins(),)),
- "H1_snr_chi": snr_chi_binning,
- "H2_snr_chi": snr_chi_binning,
- "H1H2_snr_chi": snr_chi_binning,
- "L1_snr_chi": snr_chi_binning,
- "V1_snr_chi": snr_chi_binning,
- "E1_snr_chi": snr_chi_binning,
- "E2_snr_chi": snr_chi_binning,
- "E3_snr_chi": snr_chi_binning,
- "E0_snr_chi": snr_chi_binning
}
- del snr_chi_binning
-
- def __init__(self, *args, **kwargs):
- super(ThincaCoincParamsDistributions, self).__init__(*args, **kwargs)
- self.horizon_history = horizonhistory.HorizonHistories()
- self.pdf_from_rates_func = {
- "instruments": self.pdf_from_rates_instruments,
- "H1_snr_chi": self.pdf_from_rates_snrchi2,
- "H2_snr_chi": self.pdf_from_rates_snrchi2,
- "H1H2_snr_chi": self.pdf_from_rates_snrchi2,
- "L1_snr_chi": self.pdf_from_rates_snrchi2,
- "V1_snr_chi": self.pdf_from_rates_snrchi2,
- "E1_snr_chi": self.pdf_from_rates_snrchi2,
- "E2_snr_chi": self.pdf_from_rates_snrchi2,
- "E3_snr_chi": self.pdf_from_rates_snrchi2,
- "E0_snr_chi": self.pdf_from_rates_snrchi2,
- }
+
+ def __init__(self, instruments = frozenset(("H1", "L1", "V1")), min_instruments = 2, **kwargs):
+ #
+ # check input
+ #
+
+ if min_instruments < 1:
+ raise ValueError("min_instruments=%d must be >= 1" % min_instruments)
+ if min_instruments > len(instruments):
+ raise ValueError("not enough instruments (%s) to satisfy min_instruments=%d" % (", ".join(sorted(instruments)), min_instruments))
+
+ # in the parent class this is a class attribute, but we use
+ # it as an instance attribute here
+ self.binnings = dict.fromkeys(("%s_snr_chi" % instrument for instrument in instruments), rate.NDBins((rate.ATanLogarithmicBins(3.6, 70., 600), rate.ATanLogarithmicBins(.001, 0.5, 300))))
+ self.binnings.update({
+ "instruments": rate.NDBins((snglcoinc.InstrumentBins(instruments),))
+ })
+
+ self.pdf_from_rates_func = dict.fromkeys(("%s_snr_chi" % instrument for instrument in instruments), self.pdf_from_rates_snrchi2)
+ self.pdf_from_rates_func.update({
+ "instruments": self.pdf_from_rates_instruments
+ })
+
+ # this can't be done until the binnings attribute has been
+ # populated
+ super(ThincaCoincParamsDistributions, self).__init__(**kwargs)
+
+ # record of horizon distances for all instruments in the
+ # network
+ self.horizon_history = horizonhistory.HorizonHistories((instrument, horizonhistory.NearestLeafTree()) for instrument in instruments)
+
+ # the minimum number of instruments required to form a
+ # candidate
+ self.min_instruments = min_instruments
+ assert min_instruments == 2 # FIXME: generalize
# set to True to include zero-lag histograms in background model
self.zero_lag_in_background = False
+ @property
+ def instruments(self):
+ return frozenset(self.horizon_history)
+
def __iadd__(self, other):
# NOTE: because we use custom PDF constructions, the stock
# .__iadd__() method for this class will not result in
# valid PDFs. the rates arrays *are* handled correctly by
# the .__iadd__() method, by fiat, so just remember to
# invoke .finish() to get the PDFs in shape afterwards
+ if self.instruments != other.instruments:
+ raise ValueError("incompatible instrument sets")
+ if self.min_instruments != other.min_instruments:
+ raise ValueError("incompatible minimum number of instruments")
super(ThincaCoincParamsDistributions, self).__iadd__(other)
self.horizon_history += other.horizon_history
return self
@@ -173,14 +188,10 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
#
#
- # 2D (snr, \chi^2) values. don't allow both H1 and H2 to
- # both contribute parameters to the same coinc. if both
- # have participated favour H1
+ # 2D (snr, \chi^2) values.
#
params = CoincParams(("%s_snr_chi" % event.ifo, (event.snr, event.chisq / event.snr**2)) for event in events)
- if "H2_snr_chi" in params and "H1_snr_chi" in params:
- del params["H2_snr_chi"]
#
# instrument combination
@@ -246,6 +257,8 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
print >>sys.stderr, "synthesizing signal-like (SNR, \\chi^2) distributions ..."
if df <= 0.:
raise ValueError("require df >= 0: %s" % repr(df))
+ if set(n) != self.instruments:
+ raise ValueError("n must provide a count for exactly each of %s" % ", ".join(sorted(self.instruments)))
pfs = numpy.logspace(numpy.log10(prefactors_range[0]), numpy.log10(prefactors_range[1]), 100)
for instrument, number_of_events in n.items():
binarr = rates_dict["%s_snr_chi" % instrument]
@@ -296,6 +309,8 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
# extrapolated background.
#
+ if set(n) != self.instruments:
+ raise ValueError("n must provide a count for exactly each of %s" % ", ".join(sorted(self.instruments)))
if verbose:
print >>sys.stderr, "adding tilt to (SNR, \\chi^2) background PDFs ..."
for instrument, number_of_events in n.items():
@@ -353,6 +368,7 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
zero_lag_coinc_counts = dict((instruments, self.zero_lag_rates["instruments"][instruments,]) for instruments in self.zero_lag_rates["instruments"].bins[0].centres()),
segs = segs,
delta_t = 0.005,
+ min_instruments = self.min_instruments,
verbose = verbose
).items():
self.background_rates["instruments"][instruments,] = count
@@ -369,7 +385,7 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
# instrument combinations. because the horizon distance is
# 0'ed when an instrument is off, this marginalization over
# horizon distance histories also accounts for duty cycles
- P = inspiral_extrinsics.P_instruments_given_signal(self.horizon_history)
+ P = inspiral_extrinsics.P_instruments_given_signal(self.horizon_history, min_instruments = self.min_instruments)
# populate binning from probabilities
for instruments, p in P.items():
@@ -490,13 +506,18 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
@classmethod
def from_xml(cls, xml, name):
- self = super(ThincaCoincParamsDistributions, cls).from_xml(xml, name)
- xml = self.get_xml_root(xml, name)
+ xml = cls.get_xml_root(xml, name)
+ self = super(ThincaCoincParamsDistributions, cls).from_xml(xml, name,
+ instruments = lsctables.instrument_set_from_ifos(ligolw_param.get_pyvalue(xml, u"instruments")),
+ min_instruments = ligolw_param.get_pyvalue(xml, u"min_instruments"),
+ )
self.horizon_history = horizonhistory.HorizonHistories.from_xml(xml, name)
return self
def to_xml(self, name):
xml = super(ThincaCoincParamsDistributions, self).to_xml(name)
+ xml.appendChild(ligolw_param.from_pyvalue(u"instruments", lsctables.ifos_from_instrument_set(self.instruments)))
+ xml.appendChild(ligolw_param.from_pyvalue(u"min_instruments", self.min_instruments))
xml.appendChild(self.horizon_history.to_xml(name))
return xml
@@ -665,7 +686,7 @@ class ThincaCoincParamsDistributions(snglcoinc.CoincParamsDistributions):
continue
params[key] = snr, chi2_over_snr2()
instruments.append(instrument)
- if coinc_only and len(instruments) < 2:
+ if coinc_only and len(instruments) < self.min_instruments:
continue
params["instruments"] = (frozenset(instruments),)
params.horizons = horizon_distance
@@ -736,7 +757,7 @@ class RankingData(object):
"ln_likelihood_ratio": rate.gaussian_window(8.)
}
- def __init__(self, coinc_params_distributions, instruments, sampler_coinc_params_distributions = None, process_id = None, nsamples = 1000000, min_instruments = 2, verbose = False):
+ def __init__(self, coinc_params_distributions, instruments = None, sampler_coinc_params_distributions = None, process_id = None, nsamples = 1000000, verbose = False):
self.background_likelihood_rates = {}
self.background_likelihood_pdfs = {}
self.signal_likelihood_rates = {}
@@ -745,16 +766,6 @@ class RankingData(object):
self.zero_lag_likelihood_pdfs = {}
self.process_id = process_id
- #
- # initialize binnings
- #
-
- instruments = tuple(instruments)
- for key in [frozenset(ifos) for n in range(min_instruments, len(instruments) + 1) for ifos in iterutils.choices(instruments, n)]:
- self.background_likelihood_rates[key] = rate.BinnedArray(self.binnings["ln_likelihood_ratio"])
- self.signal_likelihood_rates[key] = rate.BinnedArray(self.binnings["ln_likelihood_ratio"])
- self.zero_lag_likelihood_rates[key] = rate.BinnedArray(self.binnings["ln_likelihood_ratio"])
-
#
# bailout out used by .from_xml() class method to get an
# uninitialized instance
@@ -764,13 +775,22 @@ class RankingData(object):
return
#
- # now check input
+ # initialize binnings
#
- if min_instruments < 1:
- raise ValueError("min_instruments=%d must be >= 1" % min_instruments)
- if min_instruments > len(instruments):
- raise ValueError("not enough instruments to satisfy min_instruments")
+ if coinc_params_distributions is not None:
+ if instruments is not None and set(instruments) != coinc_params_distributions.instruments:
+ raise ValueError("instrument set does not match coinc_params_distributions (hint: don't supply instruments when initializing from ThincaCoincParamsDistributions)")
+ instruments = tuple(coinc_params_distributions.instruments)
+ elif instruments is None:
+ raise ValueError("must supply an instrument set when not initializing from a ThincaCoincParamsDistributions object")
+ else:
+ instruments = tuple(instruments)
+
+ for key in [frozenset(ifos) for n in range(coinc_params_distributions.min_instruments, len(instruments) + 1) for ifos in iterutils.choices(instruments, n)]:
+ self.background_likelihood_rates[key] = rate.BinnedArray(self.binnings["ln_likelihood_ratio"])
+ self.signal_likelihood_rates[key] = rate.BinnedArray(self.binnings["ln_likelihood_ratio"])
+ self.zero_lag_likelihood_rates[key] = rate.BinnedArray(self.binnings["ln_likelihood_ratio"])
#
# run importance-weighted random sampling to populate
diff --git a/gstlal-inspiral/python/inspiral.py b/gstlal-inspiral/python/inspiral.py
index 367834cc6e72f245e1c4427c8676edd5d21c1753..e4fd6bcfbb1612ab0455ab4acb0df590524edc60 100644
--- a/gstlal-inspiral/python/inspiral.py
+++ b/gstlal-inspiral/python/inspiral.py
@@ -516,7 +516,7 @@ class CoincsDocument(object):
class Data(object):
- def __init__(self, filename, process_params, pipeline, instruments, seg, coincidence_threshold, coinc_params_distributions, zero_lag_ranking_stats = None, marginalized_likelihood_file = None, likelihood_files_namedtuple = None, injection_filename = None, time_slide_file = None, comment = None, tmp_path = None, likelihood_snapshot_interval = None, thinca_interval = 50.0, sngls_snr_threshold = None, gracedb_far_threshold = None, gracedb_group = "Test", gracedb_search = "LowMass", gracedb_pipeline = "gstlal", gracedb_service_url = "https://gracedb.ligo.org/api/", replace_file = True, upload_auxiliary_data_to_gracedb = True, verbose = False):
+ def __init__(self, filename, process_params, pipeline, seg, coincidence_threshold, coinc_params_distributions, zero_lag_ranking_stats = None, marginalized_likelihood_file = None, likelihood_files_namedtuple = None, injection_filename = None, time_slide_file = None, comment = None, tmp_path = None, likelihood_snapshot_interval = None, thinca_interval = 50.0, sngls_snr_threshold = None, gracedb_far_threshold = None, gracedb_group = "Test", gracedb_search = "LowMass", gracedb_pipeline = "gstlal", gracedb_service_url = "https://gracedb.ligo.org/api/", replace_file = True, upload_auxiliary_data_to_gracedb = True, verbose = False):
#
# initialize
#
@@ -558,7 +558,7 @@ class Data(object):
# initialize document to hold coincs and segments
#
- self.coincs_document = CoincsDocument(filename, process_params, comment, instruments, seg, injection_filename = injection_filename, time_slide_file = time_slide_file, tmp_path = tmp_path, replace_file = replace_file, verbose = verbose)
+ self.coincs_document = CoincsDocument(filename, process_params, comment, coinc_params_distributions.instruments, seg, injection_filename = injection_filename, time_slide_file = time_slide_file, tmp_path = tmp_path, replace_file = replace_file, verbose = verbose)
#
# attach a StreamThinca instance to ourselves
diff --git a/gstlal-inspiral/python/lloidparts.py b/gstlal-inspiral/python/lloidparts.py
index 702c122053c62b7860285b0eb8ea4f3cef861f63..7867c173085b7ab365c4f870146e652d00d1188e 100644
--- a/gstlal-inspiral/python/lloidparts.py
+++ b/gstlal-inspiral/python/lloidparts.py
@@ -94,7 +94,6 @@ from gstlal import pipeio
from gstlal import pipeparts
from gstlal import reference_psd
from gstlal import simplehandler
-from gstlal.stats import horizonhistory
import lal
from lal import LIGOTimeGPS
@@ -340,15 +339,12 @@ class Handler(simplehandler.Handler):
"""
timestamp can be a float or a slice with float boundaries.
"""
- # retrieve the horizon history for this instrument or
- # initialize it if it doesn't exist yet
- try:
- horizon_history = self.dataclass.coinc_params_distributions.horizon_history[instrument]
- except KeyError:
- horizon_history = self.dataclass.coinc_params_distributions.horizon_history[instrument] = horizonhistory.NearestLeafTree()
- # NOTE: timestamps are floats here (or float slices).
- # whitener should be reporting PSDs with integer
- # timestamps, and, anyway, we don't need nanosecond
+ # retrieve the horizon history for this instrument
+ horizon_history = self.dataclass.coinc_params_distributions.horizon_history[instrument]
+ # NOTE: timestamps are floats here (or float slices), not
+ # LIGOTimeGPS. whitener should be reporting PSDs with
+ # integer timestamps so the timestamps are not naturally
+ # high precision, and, anyway, we don't need nanosecond
# precision for the horizon distance history.
horizon_history[timestamp] = horizon_distance
@@ -468,7 +464,7 @@ class Handler(simplehandler.Handler):
# set the horizon distance history to 0 at
# on-to-off transitions of whitened h(t)
if segtype == "whitehtsegments":
- if instrument in self.dataclass.coinc_params_distributions.horizon_history and self.dataclass.coinc_params_distributions.horizon_history[instrument]:
+ if self.dataclass.coinc_params_distributions.horizon_history[instrument]:
self._record_horizon_distance(instrument, slice(float(timestamp), None), 0.)
else:
self._record_horizon_distance(instrument, float(timestamp), 0.)