...
 
Commits (12)
......@@ -831,6 +831,7 @@ for output_file_number, (svd_bank_url_dict, output_url, ranking_stat_output_url,
sngls_snr_threshold = options.singles_threshold,
tag = options.job_tag,
kafka_server = options.output_kafka_server,
cluster = options.data_source in ("lvshm", "framexmit"),# we cluster when running online within thinca.
verbose = options.verbose
)
if options.verbose:
......
......@@ -110,7 +110,7 @@ while true ; do
SERVER=$(cat ${REG}) || exit 1
ZEROLAG_COUNTS_URLS="${ZEROLAG_COUNTS_URLS} ${SERVER}${ZEROLAG_COUNTS_PATH}"
done || break
# NOTE we mix in previous samples
# NOTE comment this to unmix in previous samples
if [ -f ${OUTPUT} ]; then
cp -v ${OUTPUT} prev_${OUTPUT}
PREV_OUTPUT=prev_${OUTPUT}
......@@ -118,6 +118,8 @@ while true ; do
PREV_OUTPUT=""
fi
echo "Previous output: " ${PREV_OUTPUT}
#PREV_OUTPUT=""
# ENDNOTE
gstlal_inspiral_marginalize_likelihood --verbose --marginalize ranking-stat-pdf --density-estimate-zero-lag --output ${OUTPUT} ${RANKING_PDF_FILES} ${PREV_OUTPUT} ${ZEROLAG_COUNTS_URLS} || break
date +"%H:%M:%S" >&2
rm -vf ${RANKING_PDF_FILES}
......
......@@ -45,7 +45,7 @@ def schechter(mass, maxM, alpha):
parser = argparse.ArgumentParser(description = "Create analytic mass models for prior weighting of templates")
parser.add_argument("--template-bank", metavar='name', type=str, help='The input template bank file name.', required = True)
parser.add_argument("--output", metavar='name', type=str, help='The output file name', default = "inspiral_mass_model.h5")
parser.add_argument("--model", metavar='name', type=str, help='Mass model. Options are: ligo. If you want another one, submit a patch.')
parser.add_argument("--model", metavar='name', type=str, help='Mass model. Options are: ligo, narrow_bns. If you want another one, submit a patch.')
parser.add_argument("--verbose", help='Be verbose', action="store_true")
options = parser.parse_args()
......@@ -72,7 +72,14 @@ othernorm = schechter_norm(1, 200, 200., -2.35)
for row in sngl_inspiral_table:
assert row.template_id not in ids
tmplt_ids.append(int(row.template_id))
if options.model == "ligo":
if options.model == "narrow_bns":
mchirp = chirpmass(row.mass1, row.mass2)
sigma = 0.04
mean = 1.20
ids[row.template_id] = 1. / (2 * numpy.pi * sigma**2)**.5 * numpy.exp(-(mchirp - mean)**2 / 2. / sigma**2)
mchirps.append(mchirp)
probs.append(ids[row.template_id])
elif options.model == "ligo":
#
# BNS portion
......
......@@ -136,6 +136,8 @@ def parse_command_line():
parser.add_option("--add-zerolag-to-background", action = "store_true", help = "Add zerolag events to background before populating coincident parameter PDF histograms.")
parser.add_option("--user-tag", metavar = "user-tag", default = "ALL", help = "Set the adjustable component of the description fields in the output filenames (default = \"ALL\").")
parser.add_option("--plot-snr-snr-pdfs", action = "store_true", help = "Plot the full cache of snr-snr-pdfs.")
parser.add_option("--event-snr", metavar = "ifo:snr", action = "append", help = "SNR to plot on snr chisq plots. Pass as ifo:snr, e.g. H1:3.2. Can be passed multiple times, though only once for each ifo. Must also pass --event-chisq for ifo.")
parser.add_option("--event-chisq", metavar = "ifo:chisq", action = "append", help = "chisq to plot on snr chisq plots. Pass as ifo:chisq, e.g. H1:1.1. Can be passed multiple times, though only once for each ifo. Must also pass --event-snr for ifo.")
parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")
options, filenames = parser.parse_args()
......@@ -153,6 +155,22 @@ def parse_command_line():
options.user_tag = options.user_tag.upper()
options.event_snr_dict = {}
options.event_chisq_dict = {}
if options.event_snr or options.event_chisq:
for ifo_snr_str in options.event_snr:
ifo, snr = ifo_snr_str.split(':')
options.event_snr_dict[ifo] = float(snr)
for ifo_chisq_str in options.event_chisq:
if ifo_chisq_str.split(':')[0] != ifo:
continue
options.event_chisq_dict[ifo] = float(ifo_chisq_str.split(':')[1])
for ifo in options.event_snr_dict.keys():
if ifo not in options.event_chisq_dict.keys():
print ifo
raise ValueError("Each ifo provided to --event-snr or --event-chisq must be provided to both options")
return options, filenames
......@@ -331,7 +349,10 @@ for bin_index, rankingstat in enumerate(sorted(rankingstats.values(), key = lamb
# SNR and \chi^2
for instrument in rankingstat.instruments:
for snr_chi_type in ("background_pdf", "injection_pdf", "zero_lag_pdf", "LR"):
fig = plotfar.plot_snr_chi_pdf(rankingstat, instrument, snr_chi_type, options.max_snr, sngls = sngls)
if instrument in options.event_snr_dict.keys():
fig = plotfar.plot_snr_chi_pdf(rankingstat, instrument, snr_chi_type, options.max_snr, sngls = sngls, event_snr = options.event_snr_dict[instrument], event_chisq = options.event_chisq_dict[instrument])
else:
fig = plotfar.plot_snr_chi_pdf(rankingstat, instrument, snr_chi_type, options.max_snr, sngls = sngls)
if fig is None:
continue
plotname = inspiral_pipe.T050017_filename(instrument, "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_%04d_%s_SNRCHI2" % (options.user_tag, bin_index, snr_chi_type.upper()), seg, options.output_format, path = options.output_dir)
......
......@@ -129,7 +129,7 @@ class RankingStat(snglcoinc.LnLikelihoodRatioMixin):
pass
# network SNR threshold
network_snrsq_threshold = 5.5**2.
network_snrsq_threshold = 6.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)
......@@ -682,7 +682,7 @@ WHERE
# that have too small a count to have been well measured,
# and/or can't be modelled correctly by this fit anyway.
mode, = zl.argmax()
mask = (x < mode) | (zl.at_centres() < zl[mode,] - 9.)
mask = (x < mode) | (zl.at_centres() < zl[mode,] - 15.)
zl = numpy.ma.masked_array(zl.array, mask)
bg = numpy.ma.masked_array(bg, mask)
......
......@@ -473,7 +473,7 @@ class GracedBWrapper(object):
assert self.gracedb_client is not None, ".gracedb_client is None; did you forget to set .far_threshold?"
for gracedb_id in gracedb_ids:
if self.verbose:
print >>sys.stderr, "posting '%s' to gracedb ID %d ..." % (filename, gracedb_id)
print >>sys.stderr, "posting '%s' to gracedb ID %s ..." % (filename, gracedb_id)
for attempt in range(1, self.retries + 1):
try:
resp = self.gracedb_client.writeLog(gracedb_id, message, filename = filename, filecontents = contents, tagname = tag)
......@@ -509,15 +509,19 @@ class GracedBWrapper(object):
coinc_inspiral_index = last_coincs.coinc_inspiral_index
# NOTE if any are None, this becomes None.
# Pick the "best" coinc
# FIXME revisit depending on how clustering goes
best_coinc = [min((coinc_inspiral_index[coinc_event.coinc_event_id].combined_far, coinc_event) for coinc_event in last_coincs.coinc_event_index.values())]
# This appears to be a silly for loop since
# coinc_event_index will only have one value, but we're
# future proofing this at the point where it could have
# multiple clustered events
#for coinc_event in last_coincs.coinc_event_index.values():
for _, coinc_event in best_coinc:
# NOTE if any are None, this becomes None.
# best_coinc = [min((coinc_inspiral_index[coinc_event.coinc_event_id].combined_far, coinc_event) for coinc_event in last_coincs.coinc_event_index.values())]
# NOTE streamthinca currently records the max LR and max SNR
# triggers. Both will be uploaded if they are separate. Many
# times they are the same. NOTE NOTE NOTE FIXME FIXME FIXME.
# this loop would be a disaster if stream thinca doesn't
# cluster!
for coinc_event in last_coincs.coinc_event_index.values():
# revisit this "best coinc" if clustering is removed from streamthinca
#for _, coinc_event in best_coinc:
#
# continue if the false alarm rate is not low
# enough, or is nan, or there aren't enough
......
......@@ -198,6 +198,7 @@ class EyeCandy(object):
self.kafka_data["snr_history"] = ""
self.kafka_data["far_history"] = ""
self.kafka_data["ram_history"] = ""
self.kafka_data["coinc"] = ""
for instrument in instruments:
self.kafka_data["%s_snr_history" % instrument] = ""
self.kafka_data["%s/dqvector_on" % instrument] = ""
......@@ -218,6 +219,27 @@ class EyeCandy(object):
if last_coincs:
coinc_inspiral_index = last_coincs.coinc_inspiral_index
coinc_event_index = last_coincs.coinc_event_index
sngl_inspiral_index = last_coincs.sngl_inspiral_index
coinc_dict_list = []
for coinc_event_id in coinc_event_index:
coinc_dict = {}
for attr in ("combined_far", "snr", "false_alarm_rate"):
try:
coinc_dict[attr] = float(getattr(coinc_inspiral_index[coinc_event_id], attr))
except TypeError as e:
print >>sys.stderr, e, attr, getattr(coinc_inspiral_index[coinc_event_id], attr)
coinc_dict["end"] = float(coinc_inspiral_index[coinc_event_id].end)
for attr in ("likelihood",):
try:
coinc_dict[attr] = float(getattr(coinc_event_index[coinc_event_id], attr))
except TypeError as e:
print >>sys.stderr, e, attr, getattr(coinc_event_index[coinc_event_id], attr)
for sngl_row in sngl_inspiral_index[coinc_event_id]:
for attr in ("snr", "chisq", "mass1", "mass2", "spin1z", "spin2z"):
coinc_dict["%s_%s" % (sngl_row.ifo, attr)] = float(getattr(sngl_row, attr))
coinc_dict["%s_end" % sngl_row.ifo] = float(sngl_row.end)
coinc_dict_list.append(coinc_dict)
self.kafka_data["coinc"] += json.dumps(coinc_dict_list)
for coinc_inspiral in coinc_inspiral_index.values():
# latency in .minimum_duration
# FIXME: update when a proper column is available
......@@ -597,7 +619,7 @@ class Handler(simplehandler.Handler):
dumps of segment information, trigger files and background
distribution statistics.
"""
def __init__(self, mainloop, pipeline, coincs_document, rankingstat, horizon_distance_func, gracedbwrapper, zerolag_rankingstatpdf_url = None, rankingstatpdf_url = None, ranking_stat_output_url = None, ranking_stat_input_url = None, likelihood_snapshot_interval = None, sngls_snr_threshold = None, tag = "", kafka_server = "10.14.0.112:9092", verbose = False):
def __init__(self, mainloop, pipeline, coincs_document, rankingstat, horizon_distance_func, gracedbwrapper, zerolag_rankingstatpdf_url = None, rankingstatpdf_url = None, ranking_stat_output_url = None, ranking_stat_input_url = None, likelihood_snapshot_interval = None, sngls_snr_threshold = None, tag = "", kafka_server = "10.14.0.112:9092", cluster = False, verbose = False):
"""!
@param mainloop The main application's event loop
@param pipeline The gstreamer pipeline that is being
......@@ -621,6 +643,7 @@ class Handler(simplehandler.Handler):
# None to disable periodic snapshots, otherwise seconds
self.likelihood_snapshot_interval = likelihood_snapshot_interval
self.likelihood_snapshot_timestamp = None
self.cluster = cluster
self.gracedbwrapper = gracedbwrapper
# FIXME: detangle this
......@@ -1050,7 +1073,7 @@ class Handler(simplehandler.Handler):
for absent_instrument in self.absent_instruments:
self.stream_thinca.push(absent_instrument, (), buf_timestamp)
if self.stream_thinca.push(instrument, events, buf_timestamp):
flushed_sngls = self.stream_thinca.pull(self.rankingstat, fapfar = self.fapfar, zerolag_rankingstatpdf = self.zerolag_rankingstatpdf, coinc_sieve = lambda events, offset_vector: sum(event.snr**2. for event in events) < self.rankingstat.network_snrsq_threshold)
flushed_sngls = self.stream_thinca.pull(self.rankingstat, fapfar = self.fapfar, zerolag_rankingstatpdf = self.zerolag_rankingstatpdf, coinc_sieve = lambda events, offset_vector: sum(event.snr**2. for event in events) < self.rankingstat.network_snrsq_threshold, cluster = self.cluster)
self.coincs_document.commit()
# do GraceDB alerts and update eye candy
......@@ -1165,6 +1188,10 @@ class Handler(simplehandler.Handler):
buffer in order to close off open segment intervals before
writing to disk
"""
try:
self.snapshot_output_url("%s_LLOID" % self.tag, "xml.gz", verbose = self.verbose)
except TypeError as te:
print >>sys.stderr, "Warning: couldn't build output file on checkpoint, probably there aren't any triggers: %s" % te
# FIXME: the timestamp is used to close off open segments
# and so should *not* be the timestamp of the current
# buffer, necessarily, but rather a GPS time guaranteed to
......@@ -1173,10 +1200,6 @@ class Handler(simplehandler.Handler):
# warning" configuration using the GPS time of a trigger
# buffer would be an especially bad choice.
self.segmentstracker.flush_segments_to_disk(self.tag, timestamp)
try:
self.snapshot_output_url("%s_LLOID" % self.tag, "xml.gz", verbose = self.verbose)
except TypeError as te:
print >>sys.stderr, "Warning: couldn't build output file on checkpoint, probably there aren't any triggers: %s" % te
def __get_psd_xmldoc(self):
......
......@@ -340,7 +340,7 @@ def binned_rates_from_samples(samples):
"""
lo, hi = math.floor(samples.min()), math.ceil(samples.max())
nbins = int(math.sqrt(len(samples)) / 40.)
binnedarray = rate.BinnedArray(rate.NDBins((rate.LogarithmicPlusOverflowBins(lo if lo !=0. else samples.min(), hi, nbins),)))
binnedarray = rate.BinnedArray(rate.NDBins((rate.LogarithmicBins(lo, hi, nbins),)))
for sample in samples:
binnedarray[sample,] += 1.
rate.filter_array(binnedarray.array, rate.gaussian_window(5), use_fft = False)
......
......@@ -357,6 +357,12 @@ class LnSignalDensity(LnLRDensity):
except ValueError:
return NegInf
# FIXME NOTE
# Here we put in a penalty for single detector triggers.
# Motivated sort of by the coincidence probability.
if len(snrs) == 1:
lnP -= 5
# evaluate population model
lnP += self.population_model.lnP_template_signal(template_id, max(snrs.values()))
......
......@@ -263,7 +263,7 @@ class StreamThinca(object):
return self.time_slide_graph.push(instrument, events, t_complete)
def pull(self, rankingstat, fapfar = None, zerolag_rankingstatpdf = None, coinc_sieve = None, flush = False):
def pull(self, rankingstat, fapfar = None, zerolag_rankingstatpdf = None, coinc_sieve = None, flush = False, cluster = False):
# NOTE: rankingstat is not used to compute the ranking
# statistic, it supplies the detector livetime segment
# lists to determine which triggers are eligible for
......@@ -292,6 +292,8 @@ class StreamThinca(object):
flushed = []
flushed_unused = []
self.last_coincs.clear()
max_last_coinc_lr = None
max_last_coinc_snr = None
for node, events in self.time_slide_graph.pull(newly_reported = newly_reported, flushed = flushed, flushed_unused = flushed_unused, coinc_sieve = coinc_sieve, event_collector = self.backgroundcollector, flush = flush):
# construct row objects for coinc tables.
......@@ -310,6 +312,11 @@ class StreamThinca(object):
# some tasks for zero-lag candidates
if node.is_zero_lag:
# track biggest zero lag snr and lr
if (max_last_coinc_lr is not None and coinc.likelihood > max_last_coinc_lr[1].likelihood) or max_last_coinc_lr is None:
max_last_coinc_lr = (events, coinc, coincmaps, coinc_inspiral)
if (max_last_coinc_snr is not None and coinc_inspiral.snr > max_last_coinc_snr[3].snr) or max_last_coinc_snr is None:
max_last_coinc_snr = (events, coinc, coincmaps, coinc_inspiral)
# populate ranking statistic's zero-lag
# PDFs with triggers from all zero-lag
# candidates
......@@ -317,22 +324,44 @@ class StreamThinca(object):
for event in events:
rankingstat.zerolag.increment(event)
# latency goes in minimum_duration column. NOTE:
# latency is nonsense unless running live. FIXME:
# add a proper column for latency
coinc_inspiral.minimum_duration = gps_time_now - float(coinc_inspiral.end)
# finally, append coinc to tables
if not cluster:
self.coinc_tables.append_coinc(coinc, coincmaps, coinc_inspiral)
# add events to the zero-lag ranking
# statistic histogram
if zerolag_rankingstatpdf is not None and coinc.likelihood is not None:
zerolag_rankingstatpdf.zero_lag_lr_lnpdf.count[coinc.likelihood,] += 1
# latency goes in minimum_duration column. NOTE:
# latency is nonsense unless running live. FIXME:
# add a proper column for latency
self.last_coincs.add(events, coinc, coincmaps, coinc_inspiral)
coinc_inspiral.minimum_duration = gps_time_now - float(coinc_inspiral.end)
# finally, append coinc to tables
if cluster:
# If we have a max LR trigger, append it to last coincs and add
# events to the zero-lag ranking statistic histogram
if max_last_coinc_lr is not None:
events, coinc, coincmaps, coinc_inspiral = max_last_coinc_lr
self.coinc_tables.append_coinc(coinc, coincmaps, coinc_inspiral)
self.last_coincs.add(events, coinc, coincmaps, coinc_inspiral)
# it is only the clustered LR trigger that we want to count
if zerolag_rankingstatpdf is not None and coinc.likelihood is not None:
zerolag_rankingstatpdf.zero_lag_lr_lnpdf.count[coinc.likelihood,] += 1
# If we have a max SNR trigger, append it to last coincs
if max_last_coinc_snr is not None and (max_last_coinc_lr is not None and max_last_coinc_lr[1].coinc_event_id != max_last_coinc_snr[1].coinc_event_id):
events, coinc, coincmaps, coinc_inspiral = max_last_coinc_lr
self.coinc_tables.append_coinc(coinc, coincmaps, coinc_inspiral)
self.last_coincs.add(events, coinc, coincmaps, coinc_inspiral)
self.coinc_tables.append_coinc(coinc, coincmaps, coinc_inspiral)
self.last_coincs.add(events, coinc, coincmaps, coinc_inspiral)
# add selected singles to the noise model
......@@ -361,7 +390,17 @@ class StreamThinca(object):
# we record a lot of singles that aren't really used for
# any (retained) coincs.
self.sngl_inspiral_table.extend(newly_reported)
if not cluster:
self.sngl_inspiral_table.extend(newly_reported)
else:
sngls_in_coincs_ids = set()
if max_last_coinc_lr is not None:
for e in max_last_coinc_lr[0]:
sngls_in_coincs_ids.add(e.event_id)
if max_last_coinc_snr is not None:
for e in max_last_coinc_snr[0]:
sngls_in_coincs_ids.add(e.event_id)
self.sngl_inspiral_table.extend([sngl_trigger for sngl_trigger in newly_reported if sngl_trigger.event_id in sngls_in_coincs_ids])
# save all sngls above the requested sngls SNR threshold.
# all sngls that participated in coincs are already in the
......
......@@ -313,7 +313,7 @@ class HorizonDistance(object):
f_max + delta_f,
100., # reference frequency (Hz)
None, # LAL dictionary containing accessory parameters
lalsimulation.GetApproximantFromString("SEOBNRv4_ROM")
lalsimulation.GetApproximantFromString("IMRPhenomD")
)
assert hp.data.length > 0, "huh!? h+ has zero length!"
......