plotfar.py 17 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
# Copyright (C) 2011--2014  Kipp Cannon, Chad Hanna, Drew Keppel
# Copyright (C) 2013  Jacob Peoples
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

Ryan Everett's avatar
Ryan Everett committed
18 19 20
## @file

## @package plotfar
21 22 23 24 25 26 27 28 29 30 31 32 33 34

#
# =============================================================================
#
#                                   Preamble
#
# =============================================================================
#

import warnings
import math
import matplotlib
from matplotlib import figure
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
35
matplotlib.rcParams.update({
36 37 38 39 40 41 42 43
	"font.size": 10.0,
	"axes.titlesize": 10.0,
	"axes.labelsize": 10.0,
	"xtick.labelsize": 8.0,
	"ytick.labelsize": 8.0,
	"legend.fontsize": 8.0,
	"figure.dpi": 300,
	"savefig.dpi": 300,
44 45
	"text.usetex": True
})
46 47 48 49
import numpy
from gstlal import plotutil
from gstlal import far

50 51 52 53 54 55 56
def init_plot(figsize):
	fig = figure.Figure()
	FigureCanvas(fig)
	fig.set_size_inches(figsize)
	axes = fig.gca()

	return fig, axes
57

58

59 60 61 62 63 64
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"

65 66 67 68 69 70 71
	# also checks that which is an allowed value
	tag = {
		"background_pdf": "Noise",
		"injection_pdf": "Signal",
		"zero_lag_pdf": "Candidates",
		"LR": "LR"
	}[which]
72

73 74 75 76 77 78 79 80 81
	# sngls is a sequence of {instrument: (snr, chisq)} dictionaries,
	# obtain the co-ordinates for a sngls scatter plot for this
	# instrument from that.  need to handle case in which there are no
	# singles for this instrument
	if sngls is not None:
		sngls = numpy.array([sngl[instrument] for sngl in sngls if instrument in sngl])
		if not len(sngls):
			sngls = None

82 83
	if which == "background_pdf":
		# a ln PDF object
84
		binnedarray = rankingstat.denominator.densities["%s_%s" % (instrument, base)]
85 86 87
	elif which == "injection_pdf":
		# a ln PDF object.  numerator has only one, same for all
		# instruments
88
		binnedarray = rankingstat.numerator.densities["%s" % base]
89 90
	elif which == "zero_lag_pdf":
		# a ln PDF object
91
		binnedarray = rankingstat.zerolag.densities["%s_%s" % (instrument, base)]
92
	elif which == "LR":
93 94
		num = rankingstat.numerator.densities["%s" % base]
		den = rankingstat.denominator.densities["%s_%s" % (instrument, base)]
95 96 97 98 99 100 101 102
		assert num.bins == den.bins
		binnedarray = num.count.copy()
		with numpy.errstate(invalid = "ignore"):
			binnedarray.array[:,:] = num.at_centres() - den.at_centres()
		binnedarray.array[num.count.array == 0] = float("-inf")
	else:
		raise ValueError("invalid which (%s)" % which)

103 104 105 106
	# the last bin can have a centre at infinity, and its value is
	# always 0 anyway so there's no point in trying to include it
	x = binnedarray.bins[0].centres()[:-1]
	y = binnedarray.bins[1].centres()[:-1]
107
	z = binnedarray.at_centres()[:-1,:-1]
108
	if numpy.isnan(z).any():
109 110 111
		if numpy.isnan(z).all():
			warnings.warn("%s %s is all NaN, skipping" % (instrument, which))
			return None
112
		warnings.warn("%s %s contains NaNs" % (instrument, which))
113 114 115
		z = numpy.ma.masked_where(numpy.isnan(z), z)

	# the range of the plots
116
	xlo, xhi = rankingstat.snr_min, snr_max
117 118 119 120 121 122
	ylo, yhi = .0001, 1.

	x = x[binnedarray.bins[xlo:xhi, ylo:yhi][0]]
	y = y[binnedarray.bins[xlo:xhi, ylo:yhi][1]]
	z = z[binnedarray.bins[xlo:xhi, ylo:yhi]]

123 124
	fig, axes = init_plot((8., 8. / plotutil.golden_ratio))
	if which == "LR":
125 126
		norm = matplotlib.colors.Normalize(vmin = -80., vmax = +200.)
		levels = numpy.linspace(-80., +200, 141)
127
	elif which == "background_pdf":
128 129
		norm = matplotlib.colors.Normalize(vmin = -30., vmax = z.max())
		levels = 50
130
	elif which == "injection_pdf":
131 132
		norm = matplotlib.colors.Normalize(vmin = -60., vmax = z.max())
		levels = 50
133
	elif which == "zero_lag_pdf":
134 135
		norm = matplotlib.colors.Normalize(vmin = -30., vmax = z.max())
		levels = 50
136 137
	else:
		raise ValueError("invalid which (%s)" % which)
138 139

	mesh = axes.pcolormesh(x, y, z.T, norm = norm, cmap = "afmhot", shading = "gouraud")
140
	if which == "LR":
141 142 143 144
		cs = axes.contour(x, y, z.T, levels, norm = norm, colors = "k", linewidths = .5, alpha = .3)
		axes.clabel(cs, [-20., -10., 0., +10., +20.], fmt = "%g", fontsize = 8)
	else:
		axes.contour(x, y, z.T, levels, norm = norm, colors = "k", linestyles = "-", linewidths = .5, alpha = .3)
145
	if event_snr is not None and event_chisq is not None:
Kipp Cannon's avatar
Kipp Cannon committed
146
		axes.plot(event_snr, event_chisq / event_snr / event_snr, "ko", mfc = "None", mec = "g", ms = 14, mew=4)
147 148
	if sngls is not None:
		axes.plot(sngls[:,0], sngls[:,1] / sngls[:,0]**2., "b.", alpha = .2)
149 150 151 152 153 154
	axes.loglog()
	axes.grid(which = "both")
	#axes.set_xlim((xlo, xhi))
	#axes.set_ylim((ylo, yhi))
	fig.colorbar(mesh, ax = axes)
	axes.set_xlabel(r"$\mathrm{SNR}$")
155 156 157 158
	if bankchisq:
		label = "_{\mathrm{bank}}"
	else:
		label = ""
159
	if tag.lower() in ("signal",):
160
		axes.set_title(r"$\ln P(\chi%s^{2} / \mathrm{SNR}^{2} | \mathrm{SNR}, \mathrm{%s})$ in %s" % (label, tag.lower(), instrument))
161
	elif tag.lower() in ("noise", "candidates"):
162
		axes.set_title(r"$\ln P(\mathrm{SNR}, \chi%s^{2} / \mathrm{SNR}^{2} | \mathrm{%s})$ in %s" % (label, tag.lower(), instrument))
163
	elif tag.lower() in ("lr",):
164
		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))
165 166
	else:
		raise ValueError(tag)
167 168 169
	try:
		fig.tight_layout(pad = .8)
	except RuntimeError:
170 171 172 173
		if bankchisq:
			label = "bank"
		else:
			label = ""
174
		if tag.lower() in ("signal",):
175
			axes.set_title("ln P(chi%s^2 / SNR^2 | SNR, %s) in %s" % (label, tag.lower(), instrument))
176
		elif tag.lower() in ("noise", "candidates"):
177
			axes.set_title("ln P(SNR, chi%s^2 / SNR^2 | %s) in %s" % (label, tag.lower(), instrument))
178
		elif tag.lower() in ("lr",):
179
			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))
180
		fig.tight_layout(pad = .8)
181 182 183
	return fig


184
def plot_rates(rankingstat):
185 186 187 188 189 190 191 192 193 194 195 196
	fig = figure.Figure()
	FigureCanvas(fig)
	fig.set_size_inches((6., 6.))
	axes0 = fig.add_subplot(2, 2, 1)
	axes1 = fig.add_subplot(2, 2, 2)
	axes2 = fig.add_subplot(2, 2, 3)
	axes3 = fig.add_subplot(2, 2, 4)

	# singles counts
	labels = []
	sizes = []
	colours = []
197
	for instrument, count in sorted(rankingstat.denominator.triggerrates.counts.items()):
198 199 200 201
		labels.append("%s\n(%d)" % (instrument, count))
		sizes.append(count)
		colours.append(plotutil.colour_from_instruments((instrument,)))
	axes0.pie(sizes, labels = labels, colors = colours, autopct = "%.3g%%", pctdistance = 0.4, labeldistance = 0.8)
202
	axes0.set_title("Trigger Counts")
203

204
	# singles rates
205 206 207
	labels = []
	sizes = []
	colours = []
208 209 210 211
	for instrument, rate in sorted(rankingstat.denominator.triggerrates.densities.items()):
		labels.append("%s\n(%g Hz)" % (instrument, rate))
		sizes.append(rate)
		colours.append(plotutil.colour_from_instruments((instrument,)))
212
	axes1.pie(sizes, labels = labels, colors = colours, autopct = "%.3g%%", pctdistance = 0.4, labeldistance = 0.8)
213
	axes1.set_title(r"Mean Trigger Rates (per live-time)")
214

215
	# live time
216 217 218
	labels = []
	sizes = []
	colours = []
219 220 221 222 223
	seglists = rankingstat.denominator.triggerrates.segmentlistdict()
	for instruments in sorted(rankingstat.denominator.coinc_rates.all_instrument_combos, key = lambda instruments: sorted(instruments)):
		livetime = float(abs(seglists.intersection(instruments) - seglists.union(frozenset(seglists) - instruments)))
		labels.append("%s\n(%g s)" % (", ".join(sorted(instruments)), livetime))
		sizes.append(livetime)
224 225
		colours.append(plotutil.colour_from_instruments(instruments))
	axes2.pie(sizes, labels = labels, colors = colours, autopct = "%.3g%%", pctdistance = 0.4, labeldistance = 0.8)
226
	axes2.set_title("Live-time")
227

228
	# projected background counts
229 230 231
	labels = []
	sizes = []
	colours = []
232
	for instruments, count in sorted(rankingstat.denominator.candidate_count_model().items(), key = lambda (instruments, count): sorted(instruments)):
233
		labels.append("%s\n(%d)" % (", ".join(sorted(instruments)), count))
234 235 236
		sizes.append(count)
		colours.append(plotutil.colour_from_instruments(instruments))
	axes3.pie(sizes, labels = labels, colors = colours, autopct = "%.3g%%", pctdistance = 0.4, labeldistance = 0.8)
237 238
	axes3.set_title("Expected Background Candidate Counts\n(before $\ln \mathcal{L}$ cut)")

239 240
	fig.tight_layout(pad = .8)
	return fig
241

242

Kipp Cannon's avatar
Kipp Cannon committed
243
def plot_snr_joint_pdf(snrpdf, instruments, horizon_distances, min_instruments, max_snr, sngls = None):
244 245 246 247
	if len(instruments) < 1:
		raise ValueError("len(instruments) must be >= 1")

	# retrieve the PDF in binned array form (not the interpolator)
Kipp Cannon's avatar
Kipp Cannon committed
248
	binnedarray = snrpdf.get_snr_joint_pdf_binnedarray(instruments, horizon_distances, min_instruments)
249 250

	# the range of the axes
251
	xlo, xhi = far.RankingStat.snr_min, max_snr
252 253 254
	mask = binnedarray.bins[(slice(xlo, xhi),) * len(instruments)]

	# axes are in alphabetical order
255
	instruments = sorted(instruments)
256 257 258 259 260 261 262 263 264 265 266 267 268 269

	# sngls is a sequence of {instrument: (snr, chisq)} dictionaries,
	# digest into co-ordinate tuples for a sngls scatter plot
	if sngls is not None:
		# NOTE:  the PDFs are computed subject to the constraint
		# that the candidate is observed in precisely that set of
		# instruments, so we need to restrict ourselves, here, to
		# coincs that involve the combination of instruments in
		# question otherwise we'll be overlaying a scatter plot
		# that we don't believe to have been drawn from the PDF
		# we're showing.
		sngls = numpy.array([tuple(sngl[instrument][0] for instrument in instruments) for sngl in sngls if sorted(sngl) == instruments])

	x = [binning.centres() for binning in binnedarray.bins]
270 271 272 273 274
	z = binnedarray.array
	if numpy.isnan(z).any():
		warnings.warn("%s SNR PDF for %s contains NaNs" % (", ".join(instruments), ", ".join("%s=%g" % instdist for instdist in sorted(horizon_distances.items()))))
		z = numpy.ma.masked_where(numpy.isnan(z), z)

275 276
	x = [coords[m] for coords, m in zip(x, mask)]
	z = z[mask]
277 278 279 280

	# one last check for craziness to make error messages more
	# meaningful
	assert not numpy.isnan(z).any()
281
	assert not (z < 0.).any()
282

283 284 285 286
	# plot the natural logarithm of the PDF
	with numpy.errstate(divide = "ignore"):
		z = numpy.log(z)

287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334
	if len(instruments) == 1:
		# 1D case
		fig, axes = init_plot((5., 5. / plotutil.golden_ratio))

		# FIXME:  hack to allow all-0 PDFs to be plotted.  remove
		# when we have a version of matplotlib that doesn't crash,
		# whatever version of matplotlib that is
		if numpy.isinf(z).all():
			z[:] = -60.
			z[0] = -55.

		axes.semilogx(x[0], z, color = "k")
		ylo, yhi = -40., max(0., z.max())
		if sngls is not None and len(sngls) == 1:
			axes.axvline(sngls[0, 0])
		axes.set_xlim((xlo, xhi))
		axes.set_ylim((ylo, yhi))
		axes.grid(which = "both", linestyle = "-", linewidth = 0.2)
		axes.set_xlabel(r"$\mathrm{SNR}_{\mathrm{%s}}$" % instruments[0])
		axes.set_ylabel(r"$\ln P(\mathrm{SNR}_{\mathrm{%s}})$" % instruments[0])

	elif len(instruments) == 2:
		# 2D case
		fig, axes = init_plot((5., 4.))

		# FIXME:  hack to allow all-0 PDFs to be plotted.  remove
		# when we have a version of matplotlib that doesn't crash,
		# whatever version of matplotlib that is
		if numpy.isinf(z).all():
			z[:,:] = -60.
			z[0,0] = -55.

		norm = matplotlib.colors.Normalize(vmin = -40., vmax = max(0., z.max()))

		mesh = axes.pcolormesh(x[0], x[1], z.T, norm = norm, cmap = "afmhot", shading = "gouraud")
		axes.contour(x[0], x[1], z.T, 50, norm = norm, colors = "k", linestyles = "-", linewidths = .5, alpha = .3)

		if sngls is not None and len(sngls) == 1:
			axes.plot(sngls[0, 0], sngls[0, 1], "ko", mfc = "None", mec = "g", ms = 14, mew=4)
		elif sngls is not None:
			axes.plot(sngls[:,0], sngls[:,1], "b.", alpha = .2)

		axes.loglog()
		axes.grid(which = "both", linestyle = "-", linewidth = 0.2)
		fig.colorbar(mesh, ax = axes)
		# co-ordinates are in alphabetical order
		axes.set_xlabel(r"$\mathrm{SNR}_{\mathrm{%s}}$" % instruments[0])
		axes.set_ylabel(r"$\mathrm{SNR}_{\mathrm{%s}}$" % instruments[1])
335

336 337 338
	else:
		# FIXME:  figure out how to plot 3+D PDFs
		return None
339

340
	axes.set_title(r"$\ln P(%s | \{%s\}, \mathrm{signal})$" % (", ".join("\mathrm{SNR}_{\mathrm{%s}}" % instrument for instrument in instruments), ", ".join("{D_{\mathrm{H}}}_{\mathrm{%s}}=%.3g" % item for item in sorted(horizon_distances.items()))))
341
	fig.tight_layout(pad = .8)
342
	return fig
343

344

345 346 347
def plot_likelihood_ratio_pdf(rankingstatpdf, (xlo, xhi), title, which = "noise"):
	fig, axes = init_plot((8., 8. / plotutil.golden_ratio))

Kipp Cannon's avatar
Kipp Cannon committed
348 349 350 351
	if rankingstatpdf.zero_lag_lr_lnpdf.array.any():
		extincted = rankingstatpdf.new_with_extinction()
	else:
		extincted = None
352 353 354

	if which == "noise":
		lnpdf = rankingstatpdf.noise_lr_lnpdf
Kipp Cannon's avatar
Kipp Cannon committed
355
		extinctedlnpdf = extincted.noise_lr_lnpdf if extincted is not None else None
356 357 358
		zerolag_lnpdf = rankingstatpdf.zero_lag_lr_lnpdf
	elif which == "signal":
		lnpdf = rankingstatpdf.signal_lr_lnpdf
Kipp Cannon's avatar
Kipp Cannon committed
359
		extinctedlnpdf = extincted.signal_lr_lnpdf if extincted is not None else None
360
		zerolag_lnpdf = None
361
	else:
362 363 364
		raise ValueError("invalid which (%s)" % which)

	axes.semilogy(lnpdf.bins[0].centres(), numpy.exp(lnpdf.at_centres()), color = "r", label = "%s model without extinction" % title)
Kipp Cannon's avatar
Kipp Cannon committed
365 366
	if extincted is not None:
		axes.semilogy(extinctedlnpdf.bins[0].centres(), numpy.exp(extinctedlnpdf.at_centres()), color = "k", label = "%s model with extinction" % title)
367 368
	if zerolag_lnpdf is not None:
		axes.semilogy(zerolag_lnpdf.bins[0].centres(), numpy.exp(zerolag_lnpdf.at_centres()), color = "k", linestyle = "--", label = "Observed candidate density")
369 370

	axes.grid(which = "major", linestyle = "-", linewidth = 0.2)
371
	axes.set_title(r"Ranking Statistic Distribution Density Model for %s" % title)
372
	axes.set_xlabel(r"$\ln \mathcal{L}$")
373 374 375 376 377
	axes.set_ylabel(r"$P(\ln \mathcal{L} | \mathrm{%s})$" % which)
	yhi = math.exp(lnpdf[xlo:xhi,].max())
	ylo = math.exp(lnpdf[xlo:xhi,].min())
	if zerolag_lnpdf is not None:
		yhi = max(yhi, math.exp(zerolag_lnpdf[xlo:xhi,].max()))
378 379 380
	ylo = max(yhi * 1e-40, ylo)
	axes.set_ylim((10**math.floor(math.log10(ylo) - .5), 10**math.ceil(math.log10(yhi) + .5)))
	axes.set_xlim((xlo, xhi))
381
	axes.legend(loc = "lower left", handlelength = 3)
382 383
	fig.tight_layout(pad = .8)
	return fig
384

Kipp Cannon's avatar
Kipp Cannon committed
385

386
def plot_likelihood_ratio_ccdf(fapfar, (xlo, xhi), observed_ln_likelihood_ratios = None, is_open_box = False, ln_likelihood_ratio_markers = None):
387 388
	assert xlo < xhi

389
	fig, axes = init_plot((8., 8. / plotutil.golden_ratio))
390

391
	x = numpy.linspace(xlo, xhi, int(math.ceil(xhi - xlo)) * 8)
Kipp Cannon's avatar
Kipp Cannon committed
392
	axes.semilogy(x, fapfar.fap_from_rank(x), color = "k")
393

394
	ylo = fapfar.fap_from_rank(xhi)
395 396 397
	ylo = 10**math.floor(math.log10(ylo))
	yhi = 10.

398 399 400 401
	if observed_ln_likelihood_ratios is not None:
		observed_ln_likelihood_ratios = numpy.array(observed_ln_likelihood_ratios)
		x = observed_ln_likelihood_ratios[:,0]
		y = observed_ln_likelihood_ratios[:,1]
402 403
		axes.semilogy(x, y, color = "k", linestyle = "", marker = "+", label = r"Candidates" if is_open_box else r"Candidates (time shifted)")
		axes.legend(loc = "upper right")
404

405 406
	if ln_likelihood_ratio_markers is not None:
		for ln_likelihood_ratio in ln_likelihood_ratio_markers:
407
			axes.axvline(ln_likelihood_ratio)
408

409
	axes.set_xlim((xlo, xhi))
410
	axes.set_ylim((ylo, yhi))
411
	axes.grid(which = "major", linestyle = "-", linewidth = 0.2)
412
	axes.set_title(r"False Alarm Probability vs.\ Log Likelihood Ratio")
413
	axes.set_xlabel(r"$\ln \mathcal{L}$")
414
	axes.set_ylabel(r"$P(\mathrm{one\ or\ more\ candidates} \geq \ln \mathcal{L} | \mathrm{noise})$")
415 416
	fig.tight_layout(pad = .8)
	return fig
417

Kipp Cannon's avatar
Kipp Cannon committed
418

419
def plot_horizon_distance_vs_time(rankingstat, (tlo, thi), masses = (1.4, 1.4), tref = None):
420
	fig, axes = init_plot((8., 8. / plotutil.golden_ratio))
421 422 423 424 425
	axes.xaxis.set_major_locator(matplotlib.ticker.MultipleLocator(1800.))
	axes.yaxis.set_minor_locator(matplotlib.ticker.MultipleLocator(5.))
	axes.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(50.))

	yhi = 1.
426
	for instrument, history in rankingstat.numerator.horizon_history.items():
427 428 429 430
		x = numpy.array([t for t in history.keys() if tlo <= t < thi])
		y = list(map(history.__getitem__, x))
		if tref is not None:
			x -= float(tref)
Kipp Cannon's avatar
Kipp Cannon committed
431
		axes.plot(x, y, color = plotutil.colour_from_instruments([instrument]), label = "%s" % instrument)
432 433 434 435
		try:
			yhi = max(max(y), yhi)
		except ValueError:
			pass
436
	if tref is not None:
Kipp Cannon's avatar
Kipp Cannon committed
437
		axes.set_xlabel("Time From GPS %.2f (s)" % float(tref))
438 439
	else:
		axes.set_xlim((math.floor(tlo), math.ceil(thi)))
Kipp Cannon's avatar
Kipp Cannon committed
440
		axes.set_xlabel("GPS Time (s)")
441
	axes.set_ylim((0., math.ceil(yhi / 10.) * 10.))
Kipp Cannon's avatar
Kipp Cannon committed
442 443
	axes.set_ylabel("Horizon Distance (Mpc)")
	axes.set_title(r"Horizon Distance for $%.3g\,\mathrm{M}_{\odot}$--$%.3g\,\mathrm{M}_{\odot}$ vs.\ Time" % masses)
444
	axes.grid(which = "major", linestyle = "-", linewidth = 0.2)
445
	axes.legend(loc = "lower left")
446 447
	fig.tight_layout(pad = .8)
	return fig