Commit 8702adfb authored by Ryan Magee's avatar Ryan Magee

Bankchisq: Histograms added, other python QOL changes. Still not incorporated in LR

parent 1ae71845
Pipeline #141436 passed with stages
in 37 minutes and 46 seconds
......@@ -138,6 +138,7 @@ def parse_command_line():
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("--bank-chisq", action = "store_true", help = "If provided, plots bank chisq background instead of chisq.")
parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")
options, filenames = parser.parse_args()
......@@ -350,9 +351,9 @@ for bin_index, rankingstat in enumerate(sorted(rankingstats.values(), key = lamb
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])
fig = plotfar.plot_snr_chi_pdf(rankingstat, instrument, snr_chi_type, options.max_snr, bankchisq = options.bank_chisq, 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)
fig = plotfar.plot_snr_chi_pdf(rankingstat, instrument, snr_chi_type, options.max_snr, bankchisq = options.bank_chisq, sngls = sngls)
if fig is None:
continue
plotname = dagparts.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)
......
......@@ -60,6 +60,15 @@ for bank in svd_bank.read_banks(form.getlist("url")[0], svd_bank.DefaultContentH
tmps = bank.bank_fragments[0].mix_matrix.shape[1] / 2
row = bank.sngl_inspiral_table
#bc = bank.bank_correlation_matrix
#
#sngl_inspiral = bank.sngl_inspiral_table
#
#print "<div>"
#pyplot.pcolormesh(abs(bc))
#pyplot.savefig(sys.stdout, format="svg")
#print "</div>"
for n in range(tmps):
mc, s1z, s2z = row[n].mchirp, row[n].spin1z, row[n].spin2z
......@@ -79,6 +88,7 @@ for bank in svd_bank.read_banks(form.getlist("url")[0], svd_bank.DefaultContentH
maxrate = max(maxrate, frag.rate)
mint = min(t[0], mint)
pyplot.xlabel('time (s) mc:%.2f s1z:%.2f s2z:%.2f' % (mc, s1z, s2z))
#pyplot.xlabel('time (s) mchirp = %f, number %d' % (sngl_inspiral[n].mchirp, n))
#pyplot.xlim([-.25, 0])
pyplot.subplot(212)
......@@ -86,8 +96,10 @@ for bank in svd_bank.read_banks(form.getlist("url")[0], svd_bank.DefaultContentH
l = len(bank.autocorrelation_bank[n,:])
ts = l / float(maxrate)
t = numpy.linspace(-ts, 0, l)
pyplot.plot(numpy.real(bank.autocorrelation_bank[n,:]))
pyplot.plot(numpy.imag(bank.autocorrelation_bank[n,:]))
#pyplot.plot(numpy.real(bank.autocorrelation_bank[n,:]))
#pyplot.plot(numpy.imag(bank.autocorrelation_bank[n,:]))
pyplot.plot(numpy.real(bank.bank_correlation_matrix[n,:]))
pyplot.plot(numpy.imag(bank.bank_correlation_matrix[n,:]))
pyplot.savefig(sys.stdout, format="svg")
print "</div>"
......
......@@ -509,7 +509,7 @@ SELECT
coinc_inspiral.ifos,
coinc_event.instruments,
(SELECT
group_concat(sngl_inspiral.ifo || ":" || sngl_inspiral.snr || ":" || sngl_inspiral.chisq || ":" || sngl_inspiral.mass1 || ":" || sngl_inspiral.mass2 || ":" || sngl_inspiral.spin1z || ":" || sngl_inspiral.spin2z, " ")
group_concat(sngl_inspiral.ifo || ":" || sngl_inspiral.snr || ":" || sngl_inspiral.chisq || ":" || sngl_inspiral.bank_chisq || ":" || sngl_inspiral.mass1 || ":" || sngl_inspiral.mass2 || ":" || sngl_inspiral.spin1z || ":" || sngl_inspiral.spin2z, " ")
FROM
sngl_inspiral
JOIN coinc_event_map ON (
......@@ -550,7 +550,7 @@ SELECT
coinc_inspiral.ifos,
coinc_event.instruments,
(SELECT
group_concat(sngl_inspiral.ifo || ":" || sngl_inspiral.snr || ":" || sngl_inspiral.chisq || ":" || sngl_inspiral.mass1 || ":" || sngl_inspiral.mass2 || ":" || sngl_inspiral.spin1z || ":" || sngl_inspiral.spin2z, " ")
group_concat(sngl_inspiral.ifo || ":" || sngl_inspiral.snr || ":" || sngl_inspiral.chisq || ":" || sngl_inspiral.bank_chisq || ":" || sngl_inspiral.mass1 || ":" || sngl_inspiral.mass2 || ":" || sngl_inspiral.spin1z || ":" || sngl_inspiral.spin2z, " ")
FROM
sngl_inspiral
JOIN coinc_event_map ON (
......@@ -616,10 +616,13 @@ SELECT distinct_ifos.ifos, count(*) FROM coinc_inspiral JOIN distinct_ifos ON (d
"S_2z",
"H1 SNR",
"H1 chisq",
"H1 bankchisq",
"L1 SNR",
"L1 chisq",
"L1 bankchisq",
"V1 SNR",
"V1 chisq"
"V1 chisq",
"V1 bankchisq",
]
for rank, values in enumerate(candidates, 1):
if values[0] is None:
......@@ -648,12 +651,12 @@ SELECT distinct_ifos.ifos, count(*) FROM coinc_inspiral JOIN distinct_ifos ON (d
row.extend([truncate(float(v), 3) for v in values[5:7]] + list(values[7:9]))
# populate row with trigger data
ifodict = {"H1": ['-', '-'], "L1": ['-', '-'], "V1": ['-', '-']}
ifodict = {"H1": ['-', '-', '-'], "L1": ['-', '-', '-'], "V1": ['-', '-', '-']}
template_params = ['-'] * 4
for ifo_row in values[9].split():
ifodata = ifo_row.split(":")
ifodict[ifodata[0]] = [truncate(float(v), 3) for v in ifodata[1:3]]
template_params = [truncate(float(v), 3) for v in ifodata[3:]]
ifodict[ifodata[0]] = [truncate(float(v), 3) for v in ifodata[1:4]]
template_params = [truncate(float(v), 3) for v in ifodata[4:]]
# add template params + SNR/chi^2 for all ifos
row.extend(template_params)
......@@ -1187,6 +1190,26 @@ WHERE
axes.set_xlim((self.snr_min, None))
yield fig, "chi2_vs_rho_%s" % instrument, True
fig, axes = create_plot(r"$\rho$", r"$\chi_\mathrm{bank}^{2}$")
axes.set_title(r"$\chi_\mathrm{bank}^{2}$ vs.\ $\rho$ in %s (Closed Box)" % instrument)
axes.loglog(self.injections[instrument].snr, self.injections[instrument].bankveto, 'r.', label = "Inj")
axes.loglog(self.background[instrument].snr, self.background[instrument].bankveto, "kx", label = "Background")
axes.legend(loc = "upper left")
axes.set_xlim((self.snr_min, None))
yield fig, "bankchi2_vs_rho_%s" % instrument, False
fig, axes = create_plot(r"$\rho$", r"$\chi_\mathrm{bank}^{2}$")
axes.set_title(r"$\chi_\mathrm{bank}^{2}$ vs.\ $\rho$ in %s" % instrument)
axes.loglog(self.injections[instrument].snr, self.injections[instrument].bankveto, 'r.', label = "Inj")
axes.loglog(self.background[instrument].snr, self.background[instrument].bankveto, "kx", label = "Background")
axes.loglog(self.zerolag[instrument].snr, self.zerolag[instrument].bankveto, "bx", label = "Zero-lag")
axes.legend(loc = "upper left")
axes.set_xlim((self.snr_min, None))
yield fig, "bankchi2_vs_rho_%s" % instrument, True
# the following plots require background and injections
if self.background[instrument] and self.injections[instrument]:
# for the following plots, sort injections by SNR -- descending
......@@ -1195,6 +1218,7 @@ WHERE
inj_etas = numpy.array(self.injections[instrument].eta)[sortable]
inj_snrs = numpy.array(self.injections[instrument].snr)[sortable]
inj_chisq = numpy.array(self.injections[instrument].chi2)[sortable]
inj_bankchisq = numpy.array(self.injections[instrument].bankveto)[sortable]
# sort background triggers by SNR -- ascending
sortable = numpy.array(self.background[instrument].snr).argsort()
......@@ -1202,6 +1226,7 @@ WHERE
bg_etas = numpy.array(self.background[instrument].eta)[sortable]
bg_snrs = numpy.array(self.background[instrument].snr)[sortable]
bg_chisq = numpy.array(self.background[instrument].chi2)[sortable]
bg_bankchisq = numpy.array(self.background[instrument].bankveto)[sortable]
bg_chieffs = numpy.array(self.background[instrument].chieff)[sortable]
# mass dependence of chisq
......@@ -1213,6 +1238,15 @@ WHERE
axes.semilogy()
yield fig, "chi2_vs_mc_vs_rho_%s" % instrument, False
# mass dependence of bankchisq
fig, axes = create_plot(r"$M_\mathrm{chirp}$", r"$\chi_\mathrm{bank}^{2}/\rho^2$")
coll = axes.scatter(inj_mcs, inj_bankchisq / inj_snrs**2, c = inj_snrs, label = "Injections", vmax = 20, linewidth = 0.2, alpha = 0.8)
fig.colorbar(coll, ax = axes).set_label("SNR in %s" % instrument)
axes.scatter(bg_mcs, bg_bankchisq / bg_snrs**2, edgecolor = 'k', c = bg_snrs, marker = 'x', label = "Background")
axes.legend(loc = "upper left")
axes.semilogy()
yield fig, "bankchi2_vs_mc_vs_rho_%s" % instrument, False
# background parameter distributions
fig, axes = create_plot(r"$M_\mathrm{chirp}$", r"$\eta$")
axes.set_title(r"Background Events in %s (Closed Box)" % instrument)
......
......@@ -56,7 +56,12 @@ def init_plot(figsize):
return fig, axes
def plot_snr_chi_pdf(rankingstat, instrument, which, snr_max, event_snr = None, event_chisq = None, sngls = None):
def plot_snr_chi_pdf(rankingstat, instrument, which, snr_max, bankchisq = False, event_snr = None, event_chisq = None, sngls = None):
if bankchisq:
base = "snr_bankchi"
else:
base = "snr_chi"
# also checks that which is an allowed value
tag = {
"background_pdf": "Noise",
......@@ -76,17 +81,17 @@ def plot_snr_chi_pdf(rankingstat, instrument, which, snr_max, event_snr = None,
if which == "background_pdf":
# a ln PDF object
binnedarray = rankingstat.denominator.densities["%s_snr_chi" % instrument]
binnedarray = rankingstat.denominator.densities["%s_%s" % (instrument, base)]
elif which == "injection_pdf":
# a ln PDF object. numerator has only one, same for all
# instruments
binnedarray = rankingstat.numerator.densities["snr_chi"]
binnedarray = rankingstat.numerator.densities["%s" % base]
elif which == "zero_lag_pdf":
# a ln PDF object
binnedarray = rankingstat.zerolag.densities["%s_snr_chi" % instrument]
binnedarray = rankingstat.zerolag.densities["%s_%s" % (instrument, base)]
elif which == "LR":
num = rankingstat.numerator.densities["snr_chi"]
den = rankingstat.denominator.densities["%s_snr_chi" % instrument]
num = rankingstat.numerator.densities["%s" % base]
den = rankingstat.denominator.densities["%s_%s" % (instrument, base)]
assert num.bins == den.bins
binnedarray = num.count.copy()
with numpy.errstate(invalid = "ignore"):
......@@ -147,24 +152,31 @@ def plot_snr_chi_pdf(rankingstat, instrument, which, snr_max, event_snr = None,
#axes.set_ylim((ylo, yhi))
fig.colorbar(mesh, ax = axes)
axes.set_xlabel(r"$\mathrm{SNR}$")
axes.set_ylabel(r"$\chi^{2} / \mathrm{SNR}^{2}$")
if bankchisq:
label = "_{\mathrm{bank}}"
else:
label = ""
if tag.lower() in ("signal",):
axes.set_title(r"$\ln P(\chi^{2} / \mathrm{SNR}^{2} | \mathrm{SNR}, \mathrm{%s})$ in %s" % (tag.lower(), instrument))
axes.set_title(r"$\ln P(\chi%s^{2} / \mathrm{SNR}^{2} | \mathrm{SNR}, \mathrm{%s})$ in %s" % (label, tag.lower(), instrument))
elif tag.lower() in ("noise", "candidates"):
axes.set_title(r"$\ln P(\mathrm{SNR}, \chi^{2} / \mathrm{SNR}^{2} | \mathrm{%s})$ in %s" % (tag.lower(), instrument))
axes.set_title(r"$\ln P(\mathrm{SNR}, \chi%s^{2} / \mathrm{SNR}^{2} | \mathrm{%s})$ in %s" % (label, tag.lower(), instrument))
elif tag.lower() in ("lr",):
axes.set_title(r"$\ln P(\chi^{2} / \mathrm{SNR}^{2} | \mathrm{SNR}, \mathrm{signal} ) / P(\mathrm{SNR}, \chi^{2} / \mathrm{SNR}^{2} | \mathrm{noise})$ in %s" % instrument)
axes.set_title(r"$\ln P(\chi%s^{2} / \mathrm{SNR}^{2} | \mathrm{SNR}, \mathrm{signal} ) / P(\mathrm{SNR}, \chi%s^{2} / \mathrm{SNR}^{2} | \mathrm{noise})$ in %s" % (label, label, instrument))
else:
raise ValueError(tag)
try:
fig.tight_layout(pad = .8)
except RuntimeError:
if bankchisq:
label = "bank"
else:
label = ""
if tag.lower() in ("signal",):
axes.set_title("ln P(chi^2 / SNR^2 | SNR, %s) in %s" % (tag.lower(), instrument))
axes.set_title("ln P(chi%s^2 / SNR^2 | SNR, %s) in %s" % (label, tag.lower(), instrument))
elif tag.lower() in ("noise", "candidates"):
axes.set_title("ln P(SNR, chi^2 / SNR^2 | %s) in %s" % (tag.lower(), instrument))
axes.set_title("ln P(SNR, chi%s^2 / SNR^2 | %s) in %s" % (label, tag.lower(), instrument))
elif tag.lower() in ("lr",):
axes.set_title("ln P(chi^2 / SNR^2 | SNR, signal) / P(SNR, chi^2 / SNR^2 | noise) in %s" % instrument)
axes.set_title("ln P(chi%s^2 / SNR^2 | SNR, signal) / P(SNR, chi%s^2 / SNR^2 | noise) in %s" % (label, label, instrument))
fig.tight_layout(pad = .8)
return fig
......
......@@ -114,9 +114,14 @@ class LnLRDensity(snglcoinc.LnLRDensity):
snr_min = 4.0
chi2_over_snr2_min = 1e-5
chi2_over_snr2_max = 1e20
bankchi2_over_snr2_min = 1e-5
bankchi2_over_snr2_max = 1e20
# SNR, \chi^2 binning definition
snr_chi_binning = rate.NDBins((rate.ATanLogarithmicBins(2.6, 26., 300), rate.ATanLogarithmicBins(.001, 0.2, 280)))
snr_bankchi_binning = rate.NDBins((rate.ATanLogarithmicBins(2.6, 26., 300), rate.ATanLogarithmicBins(.001, 0.2, 280)))
chi_bankchi_binning = rate.NDBins((rate.ATanLogarithmicBins(.001, 0.2, 280), rate.ATanLogarithmicBins(.001, 0.2, 280)))
chi_ratio_binning = rate.NDBins((rate.ATanLogarithmicBins(.001, 0.2, 280), rate.ATanLogarithmicBins(0.1, 10, 280)))
def __init__(self, template_ids, instruments, delta_t, min_instruments = 2):
#
......@@ -144,6 +149,11 @@ class LnLRDensity(snglcoinc.LnLRDensity):
self.densities = {}
for instrument in instruments:
self.densities["%s_snr_chi" % instrument] = rate.BinnedLnPDF(self.snr_chi_binning)
self.densities["%s_snr_bankchi" % instrument] = rate.BinnedLnPDF(self.snr_bankchi_binning)
self.densities["%s_chi_bankchi" % instrument] = rate.BinnedLnPDF(self.chi_bankchi_binning)
self.densities["%s_chi_ratio" % instrument] = rate.BinnedLnPDF(self.chi_ratio_binning)
#FIXME finish method is broken and won't make correct kernel for chi/bankchi or
# chi/ratio histograms. This should be fine until we incorporate in LR.
def __iadd__(self, other):
if type(other) != type(self):
......@@ -201,6 +211,9 @@ class LnLRDensity(snglcoinc.LnLRDensity):
if event.snr < self.snr_min:
return
self.densities["%s_snr_chi" % event.ifo].count[event.snr, event.chisq / event.snr**2.] += 1.0
self.densities["%s_snr_bankchi" % event.ifo].count[event.snr, event.bank_chisq / event.snr**2.] += 1.0
self.densities["%s_chi_bankchi" % event.ifo].count[event.chisq / event.snr**2., event.bank_chisq / event.snr**2.] += 1.0
self.densities["%s_chi_ratio" % event.ifo].count[event.chisq / event.snr**2., event.bank_chisq / event.chisq] += 1.0
def copy(self):
new = type(self)(self.template_ids, self.instruments, self.delta_t, self.min_instruments)
......@@ -921,7 +934,9 @@ class LnNoiseDensity(LnLRDensity):
# requested events to this portion of the model
arr *= number_of_events / arr.sum()
for lnpdf in self.densities.values():
# FIXME Handle this properly.
for key, lnpdf in self.densities.items():
if "chi_bankchi" not in key and "chi_ratio" not in key:
# add in the noise model
lnpdf.array += arr
# re-normalize
......
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