diff --git a/gstlal-inspiral/bin/gstlal_inspiral_plot_background b/gstlal-inspiral/bin/gstlal_inspiral_plot_background
index d0666d96c94a7a5531c62c525dbc9f8c06d13daf..c8e02429cc1c49eff9d0e7b9a25b5fe39a8315c6 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_plot_background
+++ b/gstlal-inspiral/bin/gstlal_inspiral_plot_background
@@ -48,7 +48,6 @@ from optparse import OptionParser
import os
import sqlite3
import sys
-from tqdm import tqdm
from ligo.lw import dbtables
@@ -195,31 +194,13 @@ def parse_command_line():
#
-def load_distributions(filenames, verbose = False):
- # RankingStat objects go into a dictionary indexed by the template
- # bank ID set they are for, with RankinStat objects for the same
- # templates being summed. RankingStatPDF objects are summed (and
- # therefore are required to be fore the same templates). these
- # organizational choices reflect the current use cases for this
- # program
- rankingstats = {}
- rankingstatpdf = None
- for n, filename in enumerate(filenames, 1):
- if verbose:
- print("%d/%d:" % (n, len(filenames)), file=sys.stderr)
- this_rankingstat, this_rankingstatpdf = far.parse_likelihood_control_doc(far.ligolw_utils.load_filename(filename, verbose = verbose, contenthandler = far.RankingStat.LIGOLWContentHandler))
- if this_rankingstat is not None:
- template_ids = this_rankingstat.template_ids
- if template_ids in rankingstats:
- rankingstats[template_ids] += this_rankingstat
- else:
- rankingstats[template_ids] = this_rankingstat
- if this_rankingstatpdf is not None:
- if rankingstatpdf is None:
- rankingstatpdf = this_rankingstatpdf
- else:
- rankingstatpdf += this_rankingstatpdf
- for rankingstat in tqdm(rankingstats.values(), desc = "Density estimation", disable = not verbose):
+def load_distributions(filename, verbose = False):
+ # Summing up RankingStat objects and RankingStatPDF objects for
+ # the same template bank ID set is no longer done. It is assumed
+ # that the files given as input correspond to different SVD bins,
+ # otherwise their plots will have the same filenames.
+ rankingstat, rankingstatpdf = far.parse_likelihood_control_doc(far.ligolw_utils.load_filename(filename, verbose = verbose, contenthandler = far.RankingStat.LIGOLWContentHandler))
+ if rankingstat is not None:
if options.add_zerolag_to_background:
rankingstat.denominator.lnzerolagdensity = rankingstat.zerolag
rankingstat.finish()
@@ -228,13 +209,12 @@ def load_distributions(filenames, verbose = False):
seg = (0, 0)
if rankingstatpdf is not None:
seg = rankingstatpdf.segments.extent()
- elif rankingstats:
- rankingstat = list(rankingstats.values())[0]
+ elif rankingstat:
try:
seg = rankingstat.segmentlists.extent_all()
except ValueError:
pass
- return rankingstats, rankingstatpdf, seg
+ return rankingstat, rankingstatpdf, seg
def load_search_results(filenames, tmp_path = None, verbose = False):
@@ -337,14 +317,6 @@ options, filenames = parse_command_line()
os.makedirs(options.output_dir, exist_ok=True)
-#
-# load input
-#
-
-
-rankingstats, rankingstatpdf, seg = load_distributions(filenames, verbose = options.verbose)
-
-
if options.database:
timeslide_ln_lr, zerolag_ln_lr, timeslide_sngls, zerolag_sngls = load_search_results(options.database, tmp_path = options.tmp_space, verbose = options.verbose)
if options.scatter_log_lambdas is not None:
@@ -360,71 +332,78 @@ else:
#
-for bin_index, rankingstat in enumerate(sorted(rankingstats.values(), key = lambda rankingstat: sorted(rankingstat.template_ids))):
- # SNR and \chi^2
- for instrument in rankingstat.instruments:
- for snr_chi_type in ("background_pdf", "injection_pdf", "zero_lag_pdf", "LR"):
- 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], event_bank_chisq = options.event_bank_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 = datafind.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)
- if options.verbose:
- print("writing %s" % plotname, file=sys.stderr)
- fig.savefig(plotname)
+for num, filename in enumerate(filenames):
+ if verbose:
+ print("%d/%d:" % (num, len(filenames)), file=sys.stderr)
+ rankingstat, rankingstatpdf, seg = load_distributions(filename, verbose = options.verbose)
+
+ if rankingstat is not None:
+ # FIXME bin index is inferred from filename
+ # instead, store bin_index in rankingstat object and load here
+ bin_index = int(filename.split('-')[1][:4])
- # candidate rates
- fig = plotfar.plot_rates(rankingstat)
- plotname = datafind.T050017_filename("H1L1V1", "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_%04d_RATES" % (options.user_tag, bin_index), seg, options.output_format, path = options.output_dir)
- if options.verbose:
- print("writing %s" % plotname, file=sys.stderr)
- fig.savefig(plotname)
-
-
-# SNR PDFs
-if options.plot_snr_snr_pdfs:
- # assume PDFs are identical for all template bank bins, and pick
- # one set to plot
- rankingstat = list(rankingstats.values())[0]
- for (instruments, horizon_distances) in sorted(rankingstat.numerator.SNRPDF.snr_joint_pdf_cache.keys(), key = lambda x: sorted(x[1])):
- # they're stored as a frozen set of quantized key/value
- # pairs, need to unquantize them and get a dictionary back
- horizon_distances = rankingstat.numerator.SNRPDF.quantized_horizon_distances(horizon_distances)
- fig = plotfar.plot_snr_joint_pdf(rankingstat.numerator.SNRPDF, instruments, horizon_distances, rankingstat.min_instruments, options.max_snr, sngls = sngls)
- if fig is not None:
- plotname = datafind.T050017_filename(instruments, "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_SNR_PDF_%s" % (options.user_tag, "_".join(["%s_%s" % (k, horizon_distances[k]) for k in sorted(horizon_distances)]) ), seg, options.output_format, path = options.output_dir)
+ # SNR and \chi^2
+ for instrument in rankingstat.instruments:
+ for snr_chi_type in ("background_pdf", "injection_pdf", "zero_lag_pdf", "LR"):
+ 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], event_bank_chisq = options.event_bank_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 = datafind.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)
+ if options.verbose:
+ print("writing %s" % plotname, file=sys.stderr)
+ fig.savefig(plotname)
+
+ # candidate rates
+ fig = plotfar.plot_rates(rankingstat)
+ plotname = datafind.T050017_filename("H1L1V1", "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_%04d_RATES" % (options.user_tag, bin_index), seg, options.output_format, path = options.output_dir)
+ if options.verbose:
+ print("writing %s" % plotname, file=sys.stderr)
+ fig.savefig(plotname)
+
+ # SNR PDFs
+ if options.plot_snr_snr_pdfs:
+ # assume PDFs are identical for all template bank bins, and plot
+ # only the first one
+ if num == 0:
+ for (instruments, horizon_distances) in sorted(rankingstat.numerator.SNRPDF.snr_joint_pdf_cache.keys(), key = lambda x: sorted(x[1])):
+ # they're stored as a frozen set of quantized key/value
+ # pairs, need to unquantize them and get a dictionary back
+ horizon_distances = rankingstat.numerator.SNRPDF.quantized_horizon_distances(horizon_distances)
+ fig = plotfar.plot_snr_joint_pdf(rankingstat.numerator.SNRPDF, instruments, horizon_distances, rankingstat.min_instruments, options.max_snr, sngls = sngls)
+ if fig is not None:
+ plotname = datafind.T050017_filename(instruments, "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_SNR_PDF_%s" % (options.user_tag, "_".join(["%s_%s" % (k, horizon_distances[k]) for k in sorted(horizon_distances)]) ), seg, options.output_format, path = options.output_dir)
+ if options.verbose:
+ print("writing %s" % plotname, file=sys.stderr)
+ fig.savefig(plotname)
+
+ # ranking statistic PDFs and CCDFs
+ if rankingstatpdf is not None:
+ for Title, which, NAME in (("Noise", "noise", "NOISE"), ("Signal", "signal", "SIGNAL")):
+ fig = plotfar.plot_likelihood_ratio_pdf(rankingstatpdf, (options.min_log_lambda, options.max_log_lambda), Title, which = which)
+ plotname = datafind.T050017_filename("COMBINED", "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_%s_LIKELIHOOD_RATIO_PDF" % (options.user_tag, NAME), seg, options.output_format, path = options.output_dir)
if options.verbose:
print("writing %s" % plotname, file=sys.stderr)
fig.savefig(plotname)
-
-# ranking statistic PDFs and CCDFs
-if rankingstatpdf is not None:
- for Title, which, NAME in (("Noise", "noise", "NOISE"), ("Signal", "signal", "SIGNAL")):
- fig = plotfar.plot_likelihood_ratio_pdf(rankingstatpdf, (options.min_log_lambda, options.max_log_lambda), Title, which = which)
- plotname = datafind.T050017_filename("COMBINED", "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_%s_LIKELIHOOD_RATIO_PDF" % (options.user_tag, NAME), seg, options.output_format, path = options.output_dir)
+ if rankingstatpdf is not None and rankingstatpdf.zero_lag_lr_lnpdf.array.any():
+ fapfar = far.FAPFAR(rankingstatpdf.new_with_extinction())
+ if zerolag_ln_lr is not None:
+ xhi = max(zerolag_ln_lr + timeslide_ln_lr)[0]
+ xhi = 5. * math.ceil(xhi / 5.)
+ xhi = max(xhi, options.max_log_lambda)
+ else:
+ xhi = options.max_log_lambda
+ fig = plotfar.plot_likelihood_ratio_ccdf(fapfar, (options.min_log_lambda, xhi), observed_ln_likelihood_ratios = zerolag_ln_lr, is_open_box = True)
+ plotname = datafind.T050017_filename("COMBINED", "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_NOISE_LIKELIHOOD_RATIO_CCDF_OPEN_BOX" % options.user_tag, seg, options.output_format, path = options.output_dir)
if options.verbose:
print("writing %s" % plotname, file=sys.stderr)
fig.savefig(plotname)
-if rankingstatpdf is not None and rankingstatpdf.zero_lag_lr_lnpdf.array.any():
- fapfar = far.FAPFAR(rankingstatpdf.new_with_extinction())
- if zerolag_ln_lr is not None:
- xhi = max(zerolag_ln_lr + timeslide_ln_lr)[0]
- xhi = 5. * math.ceil(xhi / 5.)
- xhi = max(xhi, options.max_log_lambda)
- else:
- xhi = options.max_log_lambda
- fig = plotfar.plot_likelihood_ratio_ccdf(fapfar, (options.min_log_lambda, xhi), observed_ln_likelihood_ratios = zerolag_ln_lr, is_open_box = True)
- plotname = datafind.T050017_filename("COMBINED", "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_NOISE_LIKELIHOOD_RATIO_CCDF_OPEN_BOX" % options.user_tag, seg, options.output_format, path = options.output_dir)
- if options.verbose:
- print("writing %s" % plotname, file=sys.stderr)
- fig.savefig(plotname)
-
- fig = plotfar.plot_likelihood_ratio_ccdf(fapfar, (options.min_log_lambda, xhi), observed_ln_likelihood_ratios = timeslide_ln_lr, is_open_box = False)
- plotname = datafind.T050017_filename("COMBINED", "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_NOISE_LIKELIHOOD_RATIO_CCDF_CLOSED_BOX" % options.user_tag, seg, options.output_format, path = options.output_dir)
- if options.verbose:
- print("writing %s" % plotname, file=sys.stderr)
- fig.savefig(plotname)
+ fig = plotfar.plot_likelihood_ratio_ccdf(fapfar, (options.min_log_lambda, xhi), observed_ln_likelihood_ratios = timeslide_ln_lr, is_open_box = False)
+ plotname = datafind.T050017_filename("COMBINED", "GSTLAL_INSPIRAL_PLOT_BACKGROUND_%s_NOISE_LIKELIHOOD_RATIO_CCDF_CLOSED_BOX" % options.user_tag, seg, options.output_format, path = options.output_dir)
+ if options.verbose:
+ print("writing %s" % plotname, file=sys.stderr)
+ fig.savefig(plotname)