Commit ebe7c676 authored by Chad Hanna's avatar Chad Hanna Committed by Patrick Godwin

gstlal_inspiral_create_prior_diststats, gstlal_inspiral_rerank_pipe, far.py...

gstlal_inspiral_create_prior_diststats, gstlal_inspiral_rerank_pipe, far.py inspiral_pipe.py, inspiral_lr.py: add IDQ likelihood
parent eef57455
Pipeline #79396 passed with stages
in 35 minutes and 33 seconds
......@@ -82,6 +82,7 @@ def parse_command_line():
parser.add_option("--svd-file", metavar = "filename", help = "The SVD file to read the template ids from")
parser.add_option("--mass-model-file", metavar = "filename", help = "The mass model file to read from (hdf5 format)")
parser.add_option("--dtdphi-file", metavar = "filename", help = "dtdphi snr ratio pdfs to read from (hdf5 format). Default passed by gstlal_inspiral_pipe, but not when run as a standalone program.")
parser.add_option("--idq-file", metavar = "filename", help = "idq glitch file (hdf5 format)")
parser.add_option("--psd-xml", type = "string", help = "Specify a PSD to use for computing template bandwidth. Required if df is bandwidth")
options, filenames = parser.parse_args()
......@@ -157,7 +158,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, dtdphi_file = options.dtdphi_file, horizon_factors = horizon_factors)
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, dtdphi_file = options.dtdphi_file, horizon_factors = horizon_factors, idq_file = options.idq_file)
if options.background_prior > 0:
rankingstat.denominator.add_noise_model(number_of_events = options.background_prior)
......
......@@ -78,6 +78,9 @@ def parse_command_line():
# dtdphi option
parser.add_option("--dtdphi-file", metavar = "filename", action = "append", help = "dtdphi snr ratio pdfs to read from (hdf5 format)")
# idq option
parser.add_option("--idq-file", metavar = "filename", action = "append", help = "idq glitch file (hdf5 format)")
# trigger generation options
parser.add_option("--web-dir", metavar = "directory", help = "Set the web directory like /home/USER/public_html")
parser.add_option("--coincidence-threshold", metavar = "value", type = "float", default = 0.005, help = "Set the coincidence window in seconds (default = 0.005). The light-travel time between instruments will be added automatically in the coincidence test.")
......@@ -247,7 +250,7 @@ if __name__ == '__main__':
model_node, model_file = inspiral_pipe.mass_model_layer(dag, jobs, [], instruments, options, boundary_seg, reference_psd)
# marginalize jobs
marg_nodes = inspiral_pipe.marginalize_layer(dag, jobs, [], lloid_output, lloid_diststats, options, boundary_seg, instrument_set, model_node, model_file, reference_psd, svd_dtdphi_map)
marg_nodes = inspiral_pipe.marginalize_layer(dag, jobs, [], lloid_output, lloid_diststats, options, boundary_seg, instrument_set, model_node, model_file, reference_psd, svd_dtdphi_map, options.idq_file)
# calc rank PDF jobs
rankpdf_nodes, rankpdf_zerolag_nodes = inspiral_pipe.calc_rank_pdf_layer(dag, jobs, marg_nodes, options, boundary_seg, instrument_set)
......
......@@ -132,8 +132,8 @@ class RankingStat(snglcoinc.LnLikelihoodRatioMixin):
# network SNR threshold
network_snrsq_threshold = 49.0
def __init__(self, template_ids = None, instruments = frozenset(("H1", "L1", "V1")), population_model_file = None, dtdphi_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, dtdphi_file = dtdphi_file, min_instruments = min_instruments, horizon_factors = horizon_factors)
def __init__(self, template_ids = None, instruments = frozenset(("H1", "L1", "V1")), population_model_file = None, dtdphi_file = None, min_instruments = 1, delta_t = 0.005, horizon_factors = None, idq_file = None):
self.numerator = inspiral_lr.LnSignalDensity(template_ids = template_ids, instruments = instruments, delta_t = delta_t, population_model_file = population_model_file, dtdphi_file = dtdphi_file, min_instruments = min_instruments, horizon_factors = horizon_factors, idq_file = idq_file)
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)
......@@ -211,9 +211,17 @@ class RankingStat(snglcoinc.LnLikelihoodRatioMixin):
return self.numerator.population_model_file
@property
def dtdphi_file(self, **kwargs):
def dtdphi_file(self):
return self.numerator.dtdphi_file
@property
def idq_file(self):
return self.numerator.idq_file
@property
def horizon_factors(self):
return self.numerator.horizon_factors
@property
def segmentlists(self):
return self.denominator.segmentlists
......@@ -403,6 +411,7 @@ class DatalessRankingStat(RankingStat):
self.numerator = inspiral_lr.DatalessLnSignalDensity(*args, **kwargs)
kwargs.pop("population_model_file", None)
kwargs.pop("dtdphi_file", None)
kwargs.pop("idq_file", None)
self.denominator = inspiral_lr.DatalessLnNoiseDensity(*args, **kwargs)
def finish(self):
......
......@@ -593,7 +593,7 @@ def merge_cluster_layer(dag, jobs, parent_nodes, db, db_cache, sqlfile, input_fi
)
def marginalize_layer(dag, jobs, svd_nodes, lloid_output, lloid_diststats, options, boundary_seg, instrument_set, model_node, model_file, ref_psd, svd_dtdphi_map):
def marginalize_layer(dag, jobs, svd_nodes, lloid_output, lloid_diststats, options, boundary_seg, instrument_set, model_node, model_file, ref_psd, svd_dtdphi_map, idq_file = None):
instruments = "".join(sorted(instrument_set))
margnodes = {}
......@@ -626,6 +626,14 @@ def marginalize_layer(dag, jobs, svd_nodes, lloid_output, lloid_diststats, optio
# FIXME we keep this here in case we someday want to have a
# mass bin dependent prior, but it really doesn't matter for
# the time being.
prior_input_files = {
"svd-file": svd_file,
"mass-model-file": model_file,
"dtdphi-file": svd_dtdphi_map[bin_key],
"psd-xml": ref_psd
}
if idq_file is not None:
prior_input_files["idq-file"] = idq_file
priornode = dagparts.DAGNode(jobs['createPriorDistStats'], dag,
parent_nodes = parent_nodes,
opts = {
......@@ -635,12 +643,7 @@ def marginalize_layer(dag, jobs, svd_nodes, lloid_output, lloid_diststats, optio
"coincidence-threshold":options.coincidence_threshold,
"df": "bandwidth"
},
input_files = {
"svd-file": svd_file,
"mass-model-file": model_file,
"dtdphi-file": svd_dtdphi_map[bin_key],
"psd-xml": ref_psd
},
input_files = prior_input_files,
output_files = {"write-likelihood": rankfile('CREATE_PRIOR_DIST_STATS', job=jobs['createPriorDistStats'])}
)
# Create a file that has the priors *and* all of the diststats
......
......@@ -80,6 +80,24 @@ __all__ = [
TYPICAL_HORIZON_DISTANCE = 150.
# utility function to make an idq interpolator from an h5 file
def idq_interp(idq_file):
# set a default function to simply return 0
out = {}
for ifo in ("H1", "L1", "V1", "K1"):
out[ifo] = lambda x: numpy.zeros(len(x))
if idq_file is not None:
import h5py
from scipy.interpolate import interp1d
f = h5py.File(idq_file, "r")
for ifo in f:
# Pad the boundaries - FIXME remove
t = [0.] + list(numpy.array(f[ifo]["time"])) + [2000000000.]
d = [0.] + list(numpy.array(f[ ifo]["data"])) + [0.]
out[ifo] = interp1d(t, d)
return out
#
# =============================================================================
#
......@@ -292,9 +310,10 @@ class LnLRDensity(snglcoinc.LnLRDensity):
class LnSignalDensity(LnLRDensity):
def __init__(self, *args, **kwargs):
population_model_file = kwargs.pop("population_model_file", None)
dtdphi_file = kwargs.pop("dtdphi_file", None)
self.population_model_file = kwargs.pop("population_model_file", None)
self.dtdphi_file = kwargs.pop("dtdphi_file", None)
self.horizon_factors = kwargs.pop("horizon_factors", None)
self.idq_file = kwargs.pop("idq_file", None)
super(LnSignalDensity, self).__init__(*args, **kwargs)
# install SNR, chi^2 PDF (one for all instruments)
self.densities = {
......@@ -304,13 +323,14 @@ class LnSignalDensity(LnLRDensity):
# record of horizon distances for all instruments in the
# network
self.horizon_history = horizonhistory.HorizonHistories((instrument, horizonhistory.NearestLeafTree()) for instrument in self.instruments)
self.population_model_file = population_model_file
self.dtdphi_file = dtdphi_file
# source population model
# FIXME: introduce a mechanism for selecting the file
self.population_model = inspiral_intrinsics.SourcePopulationModel(self.template_ids, filename = self.population_model_file)
self.InspiralExtrinsics = inspiral_extrinsics.InspiralExtrinsics(self.min_instruments, filename = self.dtdphi_file)
# initialize an IDQ likelihood of glitch interpolator
self.idq_glitch_lnl = idq_interp(self.idq_file)
def set_horizon_factors(self, horizon_factors):
self.horizon_factors = horizon_factors
......@@ -373,6 +393,14 @@ class LnSignalDensity(LnLRDensity):
# evalute the (snr, \chi^2 | snr) PDFs (same for all
# instruments)
interp = self.interps["snr_chi"]
# Evaluate the IDQ glitch probability # FIXME put in denominator
for ifo, seg in segments.items():
# FIXME don't just use the last segment, somehow include whole template duration?
t = float(seg[1])
# NOTE choose max over +-1 seconds because the sampling is only at 1 Hz.
lnP -= max(self.idq_glitch_lnl[ifo]([t-1., t, t+1.]))
return lnP + sum(interp(snrs[instrument], chi2_over_snr2) for instrument, chi2_over_snr2 in chi2s_over_snr2s.items())
def __iadd__(self, other):
......@@ -386,6 +414,10 @@ class LnSignalDensity(LnLRDensity):
raise ValueError("incompatible dtdphi files")
if self.dtdphi_file is None and other.dtdphi_file is not None:
self.dtdphi_file = other.dtdphi_file
if self.idq_file is not None and other.idq_file is not None and other.idq_file != self.idq_file:
raise ValueError("incompatible idq files")
if self.idq_file is None and other.idq_file is not None:
self.idq_file = other.idq_file
if self.horizon_factors is not None and other.horizon_factors is not None and other.horizon_factors != self.horizon_factors:
# require that the horizon factors be the same within 1%
for k in self.horizon_factors:
......@@ -404,10 +436,12 @@ class LnSignalDensity(LnLRDensity):
new.horizon_history = self.horizon_history.copy()
new.population_model_file = self.population_model_file
new.dtdphi_file = self.dtdphi_file
new.idq_file = self.idq_file
# okay to use references because read-only data
new.population_model = self.population_model
new.InspiralExtrinsics = self.InspiralExtrinsics
new.horizon_factors = self.horizon_factors
new.idq_glitch_lnl = self.idq_glitch_lnl
return new
def local_mean_horizon_distance(self, gps, window = segments.segment(-32., +2.)):
......@@ -553,6 +587,7 @@ class LnSignalDensity(LnLRDensity):
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))
xml.appendChild(ligolw_param.Param.from_pyvalue(u"dtdphi_file", self.dtdphi_file))
xml.appendChild(ligolw_param.Param.from_pyvalue(u"idq_file", self.idq_file))
return xml
@classmethod
......@@ -562,6 +597,12 @@ class LnSignalDensity(LnLRDensity):
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.dtdphi_file = ligolw_param.get_pyvalue(xml, u"dtdphi_file")
# FIXME this should be removed once there are not legacy files
try:
self.idq_file = ligolw_param.get_pyvalue(xml, u"idq_file")
except ValueError as e:
print >> sys.stderr, "idq file not found", e
self.idq_file = None
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?
......@@ -571,6 +612,7 @@ class LnSignalDensity(LnLRDensity):
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, filename = self.dtdphi_file)
self.idq_glitch_lnl = idq_interp(self.idq_file)
return self
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment