Python lalburst.SnglBurstUtils类代码示例

	def update_from(self, contents, filename = "", verbose = False):
		#
		# Add the live times from this file to the totals.
		#

		if verbose:
			print >>sys.stderr, "measuring live time ..."
		zero_lag_time_slides, background_time_slides = SnglBurstUtils.get_time_slides(contents.connection)
		self.zero_lag_live_time += SnglBurstUtils.time_slides_livetime(contents.seglists, zero_lag_time_slides.values(), verbose = verbose)
		self.background_live_time += SnglBurstUtils.time_slides_livetime(contents.seglists, background_time_slides.values(), verbose = verbose)

		#
		# Iterate over burst<-->burst coincidences.  Assume there
		# are no injections in this file.
		#

		if verbose:
			print >>sys.stderr, "retrieving sngl_burst<-->sngl_burst coincs ..."
		for id, likelihood, confidence, is_background in bb_id_likelihood_confidence_background(contents):
			record = coinc_detection_statistic(likelihood, confidence)
			if is_background:
				# remember the total number, but only keep
				# the highest ranked
				self.n_background_amplitudes += 1
				if len(self.background_amplitudes) < 1e7:
					heapq.heappush(self.background_amplitudes, record)
				else:
					heapq.heappushpop(self.background_amplitudes, record)
			else:
				self.zero_lag_amplitudes.append((record, filename, id))
		if verbose:
			print >>sys.stderr, "done"

	def finish(self):
		for instrument, data in sorted(self.dataA.items()):
			fig, axes = SnglBurstUtils.make_burst_plot(r"SNR in %s" % instrument, r"$(f_{\mathrm{recovered}} - f_{\mathrm{injected}})/ f_{\mathrm{injected}}$")

			axes.set_title("Cut-off Frequency Residual vs.\ SNR in %s\n(%d Found Injections)" % (instrument, data.found))
			axes.loglog()
			axes.plot(data.x, data.y, "k+")

			yield fig

		for instrument, data in sorted(self.dataB.items()):
			fig, axes = SnglBurstUtils.make_burst_plot(r"$f_{\mathrm{injected}}$ (Hz)", r"$f_{\mathrm{recovered}}$ (Hz)")
			axes.loglog()
			fig.set_size_inches(8, 8)
			axes.scatter(data.x, data.y, c = data.c, s = 16)
			xmin, xmax = min(data.x), max(data.x)
			ymin, ymax = min(data.y), max(data.y)
			xmin = ymin = min(xmin, ymin)
			xmax = ymax = max(xmax, ymax)
			axes.plot([xmin, xmax], [ymin, ymax], "k-")
			axes.set_xlim([xmin, xmax])
			axes.set_ylim([ymin, ymax])
			axes.set_title(r"Recovered Cut-off Frequency vs.\ Injected Cut-off Frequency in %s\\(%d Injections; red = high recovered SNR, blue = low recovered SNR)" % (instrument, data.found))

			yield fig

def measure_efficiency(filenames, threshold, live_time_program = "lalapps_power", upper_limit_scale = "E", tmp_path = None, verbose = False):
	# FIXME:  instruments are hard-coded.  bad bad bad.  sigh...
	if upper_limit_scale == "E":
		efficiency = EfficiencyData(("H1", "H2", "L1"), (lambda sim, instrument: sim.egw_over_rsquared), r"Equivalent Isotropic Energy ($M_{\odot} / \mathrm{pc}^{2}$)", 0.1)
	elif upper_limit_scale == "hrss":
		efficiency = EfficiencyData(("H1", "H2", "L1"), (lambda sim, instrument: sim.hrss), r"$h_{\mathrm{rss}}$", 0.1)
	else:
		raise ValueError("bad upper_limit_scale %s" % repr(upper_limit_scale))

	#
	# Iterate over injection files.
	#

	for n, filename in enumerate(filenames):
		#
		# Open the database file.
		#

		if verbose:
			print >>sys.stderr, "%d/%d: %s" % (n + 1, len(filenames), filename)
		working_filename = dbtables.get_connection_filename(filename, tmp_path = tmp_path, verbose = verbose)
		connection = sqlite3.connect(working_filename)
		connection.create_function("coinc_detection_statistic", 2, coinc_detection_statistic)
		database = SnglBurstUtils.CoincDatabase(connection, live_time_program)
		if verbose:
			SnglBurstUtils.summarize_coinc_database(database)

		#
		# Process database contents.
		#

		efficiency.add_contents(database, threshold)

		#
		# Done with this file.
		#

		database.connection.close()
		dbtables.discard_connection_filename(filename, working_filename, verbose = verbose)

	#
	# Compute efficiency from the data that have been collected
	#

	if verbose:
		print >>sys.stderr, "binning and smoothnig efficiency data ..."
	efficiency.finish(threshold)

	#
	# Done
	#

	return efficiency

def measure_threshold(filenames, n_survivors, live_time_program = "lalapps_power", tmp_path = None, open_box = False, verbose = False):
	#
	# Initialize the book-keeping object.
	#

	rate_vs_threshold_data = RateVsThresholdData()

	#
	# Iterate over non-injection files.
	#

	for n, filename in enumerate(filenames):
		#
		# Open the database file.
		#

		if verbose:
			print >>sys.stderr, "%d/%d: %s" % (n + 1, len(filenames), filename)
		working_filename = dbtables.get_connection_filename(filename, tmp_path = tmp_path, verbose = verbose)
		database = SnglBurstUtils.CoincDatabase(sqlite3.connect(working_filename), live_time_program)
		if verbose:
			SnglBurstUtils.summarize_coinc_database(database)

		#
		# Process database contents.
		#

		rate_vs_threshold_data.update_from(database, filename = filename, verbose = verbose)

		#
		# Done with this file.
		#

		database.connection.close()
		dbtables.discard_connection_filename(filename, working_filename, verbose = verbose)

	#
	# Determine likelihood threshold.
	#

	if verbose:
		print >>sys.stderr, "finishing rate vs. threshold measurement ..."
	rate_vs_threshold_data.finish(n_survivors, open_box = open_box, verbose = verbose)

	#
	# Done.
	#

	return rate_vs_threshold_data

	def __init__(self, ifo):
		self.fig, self.axes = SnglBurstUtils.make_burst_plot("Confidence", "Trigger Rate (Hz)")
		self.ifo = ifo
		self.nevents = 0
		self.x = []
		self.seglist = segments.segmentlist()
		self.axes.loglog()

def plot_rate_upper_limit(rate_data):
	print >>sys.stderr, "plotting rate upper limit ..."

	#
	# contour plot in frequency-energy plane
	#

	fig, axes = SnglBurstUtils.make_burst_plot("Centre Frequency (Hz)", rate_data.amplitude_lbl)
	axes.loglog()

	xcoords, ycoords = rate_data.rate_array.centres()
	# zvals = log_{10}(R_{0.9} / 1 Hz)
	zvals = numpy.log10(numpy.clip(rate_data.rate_array.array, 0, 1) / 1.0)
	cset = axes.contour(xcoords, ycoords, numpy.transpose(zvals), 19)
	cset.clabel(inline = True, fontsize = 5, fmt = r"$%g$", colors = "k")

	#axes.set_ylim((3e-17, 3e-10))

	axes.set_title(r"%g\%% Confidence Rate Upper Limit ($\log_{10} R_{%g} / 1\,\mathrm{Hz}$)" % (100 * rate_data.confidence, rate_data.confidence))

	#print >>sys.stderr, "writing lalapps_excesspowerfinal_upper_limit_1.pdf ..."
	#fig.savefig("lalapps_excesspowerfinal_upper_limit_1.pdf")
	print >>sys.stderr, "writing lalapps_excesspowerfinal_upper_limit_1.png ..."
	fig.savefig("lalapps_excesspowerfinal_upper_limit_1.png")

	#
	# rate vs. energy curve at a sample frequency
	#

	fig, axes = SnglBurstUtils.make_burst_plot(rate_data.amplitude_lbl, r"Rate Upper Limit (Hz)")
	axes.loglog()

	zvals = rate_data.rate_array[110, :]
	max = min(zvals) * 1e4
	ycoords = numpy.compress(zvals <= max, ycoords)
	zvals = numpy.compress(zvals <= max, zvals)
	max = min(ycoords) * 1e6
	zvals = numpy.compress(ycoords <= max, zvals)
	ycoords = numpy.compress(ycoords <= max, ycoords)
	axes.plot(ycoords, zvals, "k-")

	axes.set_title(r"%g\%% Confidence Rate Upper Limit at $110\,\mathrm{Hz}$" % (100 * rate_data.confidence))

	#print >>sys.stderr, "writing lalapps_excesspowerfinal_upper_limit_2.pdf ..."
	#fig.savefig("lalapps_excesspowerfinal_upper_limit_2.pdf")
	print >>sys.stderr, "writing lalapps_excesspowerfinal_upper_limit_2.png ..."
	fig.savefig("lalapps_excesspowerfinal_upper_limit_2.png")

def process_file(filename, products, live_time_program, tmp_path = None, veto_segments_name = None, verbose = False):
	#
	# connect to database and summarize contents
	#

	working_filename = dbtables.get_connection_filename(filename, tmp_path = tmp_path, verbose = verbose)
	contents = SnglBurstUtils.CoincDatabase(sqlite3.connect(working_filename), live_time_program, search = "StringCusp", veto_segments_name = veto_segments_name)
	if verbose:
		SnglBurstUtils.summarize_coinc_database(contents, filename = working_filename)

	#
	# augment summary with extra stuff we need.  the filename
	# is recorded for dumping debuggin information related to
	# missed injections.  if burca was run with the
	# --coincidence-segments option then the value is copied
	# into a segmentlistdict to facilitate the computation of
	# livetime
	#

	contents.filename = filename

	contents.coincidence_segments = ligolwprocess.get_process_params(contents.xmldoc, "lalapps_burca", "--coincidence-segments")
	if contents.coincidence_segments:
		# as a side-effect, this enforces the rule that
		# burca has been run on the input file exactly once
		contents.coincidence_segments, = contents.coincidence_segments
		# FIXME:  remove when LAL accepts unicode
		contents.coincidence_segments = contents.coincidence_segments.encode('utf-8')
		contents.coincidence_segments = segments.segmentlistdict.fromkeys(contents.seglists, segmentsUtils.from_range_strings(contents.coincidence_segments.split(","), boundtype = dbtables.lsctables.LIGOTimeGPS).coalesce())
	else:
		contents.coincidence_segments = None

	#
	# process contents
	#

	for n, product in enumerate(products):
		if verbose:
			print("%s: adding to product %d ..." % (working_filename, n), file=sys.stderr)
		product.add_contents(contents, verbose = verbose)

	#
	# close
	#

	contents.connection.close()
	dbtables.discard_connection_filename(filename, working_filename, verbose = verbose)

	def __init__(self, x_instrument, y_instrument):
		self.fig, self.axes = SnglBurstUtils.make_burst_plot("%s Offset (s)" % x_instrument, "%s Offset (s)" % y_instrument)
		self.fig.set_size_inches(6,6)
		self.x_instrument = x_instrument
		self.y_instrument = y_instrument
		self.tisi_rows = None
		self.seglists = segments.segmentlistdict()
		self.counts = None

	def coords(self, contents, coinc_event_id):
		x = y = None
		for burst in SnglBurstUtils.coinc_sngl_bursts(contents, coinc_event_id):
			if burst.ifo == self.x_instrument:
				x = self.magnitude(burst)
			elif burst.ifo == self.y_instrument:
				y = self.magnitude(burst)
		return x, y

def diagnostic_plot(z, bins, title, ylabel, filename):
	print >>sys.stderr, "generating %s ..." % filename
	fig, axes = SnglBurstUtils.make_burst_plot("Centre Frequency (Hz)", ylabel)
	axes.loglog()
	xcoords, ycoords = bins.centres()
	cset = axes.contour(xcoords, ycoords, numpy.transpose(z), 9)
	cset.clabel(inline = True, fontsize = 5, fmt = r"$%g$", colors = "k")
	axes.set_title(title)
	fig.savefig(filename)

	def __init__(self, instruments):
		self.fig, self.axes = SnglBurstUtils.make_burst_plot(r"$f_{\mathrm{injected}}$ (Hz)", r"$f_{\mathrm{recovered}}$ (Hz)")
		self.axes.loglog()
		self.fig.set_size_inches(8, 8)
		self.instruments = set(instruments)
		self.matches = 0
		self.x = []
		self.y = []
		self.c = []

def plot_multi_Efficiency_hrss_vs_freq(efficiencies):
	fig, axes = SnglBurstUtils.make_burst_plot("Frequency (Hz)", r"$h_{\mathrm{rss}}$")
	axes.loglog()

	e = efficiencies[0]
	xcoords, ycoords = e.efficiency.centres()
	zvals = e.efficiency.ratio()
	error = e.error
	for n, e in enumerate(efficiencies[1:]):
		error += e.error
		other_xcoords, other_ycoords = e.efficiency.centres()
		if (xcoords != other_xcoords).any() or (ycoords != other_ycoords).any():
			# binnings don't match, can't compute product of
			# efficiencies
			raise ValueError("binning mismatch")
		zvals *= e.efficiency.ratio()
	error /= len(efficiencies)

	nfound = numpy.array([len(e.found_x) for e in efficiencies], dtype = "double")
	ninjected = numpy.array([len(e.injected_x) for e in efficiencies], dtype = "double")

	# the model for guessing the ninjected in the coincidence case is
	# to assume that the injections done in the instrument with the
	# most injections were done into all three with a probability given
	# by the ratio of the number actually injected into each instrument
	# to the number injected into the instrument with the most
	# injected, and then to assume that these are independent random
	# occurances and that to be done in coincidence an injection must
	# be done in all three instruments.

	ninjected_guess = (ninjected / ninjected.max()).prod() * ninjected.min()

	# the model for guessing the nfound in the coincidence case is to
	# assume that each injection is found or missed in each instrument
	# at random, and to be found in the coincidence case it must be
	# found in all three.

	nfound_guess = (nfound / ninjected).prod() * ninjected_guess

	ninjected_guess = int(round(ninjected_guess))
	nfound_guess = int(round(nfound_guess))

	instruments = r" \& ".join(sorted("+".join(sorted(e.instruments)) for e in efficiencies))

	cset = axes.contour(xcoords, ycoords, numpy.transpose(zvals), (.1, .2, .3, .4, .5, .6, .7, .8, .9))
	cset.clabel(inline = True, fontsize = 5, fmt = r"$%%g \pm %g$" % error, colors = "k")

	axes.set_title(r"%s Estimated Injection Detection Efficiency ($\sim$%d of $\sim$%d Found)" % (instruments, nfound_guess, ninjected_guess))

	return fig

	def add_contents(self, contents):
		# this outer loop assumes each injection can only be found
		# in at most one coinc, otherwise the "found" count is
		# wrong.
		for values in contents.connection.cursor().execute("""
SELECT
	sim_burst.*,
	coinc_event.coinc_event_id
FROM
	coinc_event
	JOIN coinc_event_map ON (
		coinc_event_map.coinc_event_id == coinc_event.coinc_event_id
	)
	JOIN sim_burst ON (
		coinc_event_map.table_name == 'sim_burst'
		AND coinc_event_map.event_id == sim_burst.simulation_id
	)
WHERE
	coinc_def_id == ?
		""", (contents.sb_definer_id,)):
			sim = contents.sim_burst_table.row_from_cols(values)
			coinc_event_id = values[-1]
			sim_peak = sim.time_at_instrument(self.instrument)
			self.found += 1
			bursts = tuple(SnglBurstUtils.coinc_sngl_bursts(contents, coinc_event_id))
			coinc_dt = 0
			for burst in bursts:
				dt = float(burst.peak - sim_peak)
				if burst.ifo == self.instrument:
					try:
						self.offsets.count[dt,] += 1.0
					except IndexError:
						# outside plot range
						pass
				coinc_dt += dt * burst.ms_snr
			coinc_dt /= sum(burst.ms_snr for burst in bursts)
			try:
				self.coinc_offsets.count[coinc_dt,] += 1.0
			except IndexError:
				# outside plot range
				pass