Commit 7192f739 authored by Kipp Cannon's avatar Kipp Cannon
Browse files

lalapps: delete obsolete power programs

- lalapps_power_plot_burca
- lalapps_power_plot_burst
- lalapps_power_veto
- lalapps_power_plot_burca2
parent 50bc561e
......@@ -175,12 +175,8 @@ src/power/lalapps_power_online_pipe
src/power/lalapps_power_pipe
src/power/lalapps_power_plot_binj
src/power/lalapps_power_plot_binjtf
src/power/lalapps_power_plot_burca
src/power/lalapps_power_plot_burca2
src/power/lalapps_power_plot_burst
src/power/lalapps_power_plot_burstrate
src/power/lalapps_power_plot_detresponse
src/power/lalapps_power_veto
src/power/lalapps_simburst_to_frame
src/power/lalapps_xml_plotlalseries
src/pulsar/CreateEphemeris/lalapps_create_solar_system_ephemeris
......
......@@ -151,10 +151,6 @@ test:
- lalapps_power_pipe --help
- lalapps_power_plot_binj --help
- lalapps_power_plot_binjtf --help
- lalapps_power_plot_burca --help
- lalapps_power_plot_burca2 --help
- lalapps_power_plot_burst --help
- lalapps_power_veto --help
- lalapps_pulsar_frequency_evolution --help
- lalapps_pulsar_parameter_estimation_nested --help
- lalapps_randombank --help
......
......@@ -30,11 +30,7 @@ pybin_scripts = \
lalapps_power_likelihood_pipe \
lalapps_power_pipe \
lalapps_power_plot_binjtf \
lalapps_power_plot_burca2 \
lalapps_power_plot_burst \
lalapps_power_plot_binj \
lalapps_power_plot_burca \
lalapps_power_veto
lalapps_power_plot_binj
pkgpython_PYTHON = power.py
endif
......
This diff is collapsed.
#
# Copyright (C) 2006 Kipp Cannon
#
# 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.
#
# =============================================================================
#
# Preamble
#
# =============================================================================
#
from __future__ import print_function
try:
from fpconst import PosInf, NegInf
except ImportError:
# fpconst is not part of the standard library and might not be
# available
PosInf = float("+inf")
NegInf = float("-inf")
import math
import numpy
from optparse import OptionParser
import sys
from lalburst import git_version
from lalburst import burca_tailor
from lalburst import SnglBurstUtils
__author__ = "Kipp Cannon <kipp.cannon@ligo.org>"
__version__ = "git id %s" % git_version.id
__date__ = git_version.date
#
# =============================================================================
#
# Command Line
#
# =============================================================================
#
def parse_command_line():
parser = OptionParser(
version = "Name: %%prog\n%s" % git_version.verbose_msg
)
parser.add_option("-b", "--base", metavar = "base", default = "plotburca2", help = "Set the prefix for output filenames (default = \"plotburca2\").")
parser.add_option("-f", "--format", metavar = "format", default = "png", help = "Set the output image format (default = \"png\").")
parser.add_option("-l", "--live-time-program", metavar = "program", default = "lalapps_power", help = "Set the name, as it appears in the process table, of the program whose search summary entries define the search live time (default = \"lalapps_power\").")
parser.add_option("--with-zero-lag", action = "store_true", help = "Also plot the parameter distributions for zero lag events.")
parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")
options, filenames = parser.parse_args()
return options, (filenames or [None])
#
# =============================================================================
#
# Coincidence Parameter Distribution Plot
#
# =============================================================================
#
def add_to_plot(axes, red, black, blue, gmst, plottype = "LR"):
if plottype == "P":
if blue is not None:
y = blue[:,gmst]
l, = axes.plot(numpy.compress(y > 0, blue.centres()[0]), numpy.compress(y > 0, y), "b")
y = black[:,gmst]
l, = axes.plot(numpy.compress(y > 0, black.centres()[0]), numpy.compress(y > 0, y), "k")
#l.set_dashes([8, 4])
y = red[:,gmst]
l, = axes.plot(numpy.compress(y > 0, red.centres()[0]), numpy.compress(y > 0, y), "r")
elif plottype == "LR":
y1 = red[:,gmst]
y2 = black[:,gmst]
l, = axes.plot(numpy.compress((y1 > 0) & (y2 > 0), black.centres()[0]), numpy.compress((y1 > 0) & (y2 > 0), y1 / y2), "k")
else:
raise ValueError(plottype)
ymin, ymax = axes.get_ylim()
if ymin < 1e-5 * ymax:
ymin = 1e-5 * ymax
axes.set_ylim((ymin, ymax))
axes.set_xlim((-2.0, +2.0))
def add_to_dtplot(axes, red, black, blue, gmst, plottype = "LR"):
x = numpy.arange(len(black[:,gmst]), dtype = "double")
if plottype == "P":
if blue is not None:
y = blue[:,gmst]
l, = axes.plot(numpy.compress(y > 0, x), numpy.compress(y > 0, y), "b")
y = black[:,gmst]
l, = axes.plot(numpy.compress(y > 0, x), numpy.compress(y > 0, y), "k")
#l.set_dashes([8, 4])
y = red[:,gmst]
l, = axes.plot(numpy.compress(y > 0, x), numpy.compress(y > 0, y), "r")
elif plottype == "LR":
y1 = red[:,gmst]
y2 = black[:,gmst]
l, = axes.plot(numpy.compress((y1 > 0) & (y2 > 0), x), numpy.compress((y1 > 0) & (y2 > 0), y1 / y2), "k")
else:
raise ValueError(plottype)
ymin, ymax = axes.get_ylim()
if ymin < 1e-5 * ymax:
ymin = 1e-5 * ymax
axes.set_ylim((ymin, ymax))
axes.set_xlim((0, len(x) - 1))
if black.bins[0].scale > 150:
xticks = (NegInf, -0.02, -0.01, -0.005, -0.002, 0, 0.002, 0.005, 0.01, 0.02, PosInf)
else:
xticks = (NegInf, -0.05, -0.02, -0.01, -0.005, 0, 0.005, 0.01, 0.02, 0.05, PosInf)
xticklabels = [r"$-\infty$"] + [r"$%.3g$" % (1000 * x) for x in xticks[1:-1]] + [r"$\infty$"]
axes.set_xticks([black.bins[x,:][0] for x in xticks])
axes.set_xticklabels(xticklabels)
#
# distributions is a burca_tailor.BurcaCoincParamsDistributions instance
#
def plot_coinc_params(distributions, gmst, plottype = "LR", with_zero_lag = False):
#
# Create a figure.
#
fig = SnglBurstUtils.figure.Figure()
SnglBurstUtils.FigureCanvas(fig)
#
# How many instrument pairs are there?
#
pairs = set(tuple(name.split("_")[:2]) for name in distributions.denominator.densities if "_" in name)
#
# How many plots in each row?
#
n_horiz = len(pairs)
#
# How many rows?
#
n_vert = 5
#
# Each of the first len(pairs) sub plot's aspect ratio is the golden
# ratio.
#
size = 4.0
fig.set_size_inches(n_horiz * size, n_vert / ((1 + math.sqrt(5)) / 2) * size)
#
# Plot layout.
#
vlabel_allowance = (n_vert * .05) / fig.figheight.get()
hlabel_allowance = (n_horiz * .1) / fig.figwidth.get()
border = .07 / fig.figwidth.get()
width = 1.0 / n_horiz - hlabel_allowance - 2 * border
height = 1.0 / n_vert
#
# Iterate over instrument pairs.
#
for i, pair in enumerate(pairs):
#
# Construct the axes for this instrument pair.
#
left = float(i) / n_horiz + hlabel_allowance + border
dt_axes = fig.add_axes((left, 0 * height + vlabel_allowance + border, width, height - vlabel_allowance - 2 * border))
df_axes = fig.add_axes((left, 1 * height + vlabel_allowance + border, width, height - vlabel_allowance - 2 * border))
dh_axes = fig.add_axes((left, 2 * height + vlabel_allowance + border, width, height - vlabel_allowance - 2 * border))
dband_axes = fig.add_axes((left, 3 * height + vlabel_allowance + border, width, height - vlabel_allowance - 2 * border))
ddur_axes = fig.add_axes((left, 4 * height + vlabel_allowance + border, width, height - vlabel_allowance - 2 * border))
dt_axes.semilogy()
df_axes.semilogy()
dh_axes.semilogy()
dband_axes.semilogy()
ddur_axes.semilogy()
dt_axes.set_xlabel(r"$t_{\mathrm{%s}} - t_{\mathrm{%s}}$ (ms)" % pair)
df_axes.set_xlabel(r"$(f_{\mathrm{%s}} - f_{\mathrm{%s}}) / \left< f \right>$" % pair)
dh_axes.set_xlabel(r"$({h_{\mathrm{rss}}}_{\mathrm{%s}} - {h_{\mathrm{rss}}}_{\mathrm{%s}}) / \left< h_{\mathrm{rss}} \right>$" % pair)
dband_axes.set_xlabel(r"$(\Delta f_{\mathrm{%s}} - \Delta f_{\mathrm{%s}}) / \left< \Delta f \right>$" % pair)
ddur_axes.set_xlabel(r"$(\Delta t_{\mathrm{%s}} - \Delta t_{\mathrm{%s}}) / \left< \Delta t \right>$" % pair)
#
# Plot the data on them.
#
prefix = "%s_%s_" % pair
for axes, suffix in zip((df_axes, dh_axes, dband_axes, ddur_axes), ("df", "dh", "dband", "ddur")):
red = distributions.numerator.densities[prefix + suffix]
black = distributions.denominator.densities[prefix + suffix]
if with_zero_lag:
blue = distributions.candidates.densities[prefix + suffix]
else:
blue = None
add_to_plot(axes, red, black, blue, gmst, plottype = plottype)
suffix = "dt"
red = distributions.numerator.densities[prefix + suffix]
black = distributions.denominator.densities[prefix + suffix]
if with_zero_lag:
blue = distributions.candidates.densities[prefix + suffix]
else:
blue = None
add_to_dtplot(dt_axes, red, black, blue, gmst, plottype = plottype)
#
# Done.
#
return fig
#
# =============================================================================
#
# Main
#
# =============================================================================
#
options, filenames = parse_command_line()
filenames.sort()
distributions, seglists = burca_tailor.EPGalacticCoreCoincParamsDistributions.from_filenames(filenames, u"lalapps_burca_tailor", verbose = options.verbose)
if options.verbose:
print("applying filters ...", file=sys.stderr)
distributions.finish()
for gmst in numpy.arange(0.0, 2 * math.pi, 2 * math.pi / 10):
filename = "%s-%%s-%d-%d-%s.%s" % (options.base, int(seglists.extent_all()[0]), int(abs(seglists.extent_all())), ("%.2f" % gmst).replace(".", "_"), options.format)
if options.verbose:
print("writing %s ..." % (filename % "P"), file=sys.stderr)
plot_coinc_params(distributions, gmst, plottype = "P", with_zero_lag = options.with_zero_lag).savefig(filename % "P")
if options.verbose:
print("writing %s ..." % (filename % "LR"), file=sys.stderr)
plot_coinc_params(distributions, gmst, plottype = "LR").savefig(filename % "LR")
if options.verbose:
print("done.", file=sys.stderr)
#
# Copyright (C) 2006 Kipp Cannon
#
# 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.
#
# =============================================================================
#
# Preamble
#
# =============================================================================
#
from __future__ import print_function
import math
from optparse import OptionParser
import numpy
import sys
from ligo.lw import ligolw
from ligo.lw import lsctables
from ligo.lw import utils as ligolw_utils
from ligo.lw.utils import search_summary as ligolw_search_summary
from lal import rate
from lal.utils import CacheEntry
from lalburst import git_version
from lalburst import SnglBurstUtils
from ligo import segments
lsctables.use_in(ligolw.LIGOLWContentHandler)
__author__ = "Kipp Cannon <kipp.cannon@ligo.org>"
__version__ = "git id %s" % git_version.id
__date__ = git_version.date
#
# =============================================================================
#
# Command Line
#
# =============================================================================
#
def parse_command_line():
parser = OptionParser(
version = "Name: %%prog\n%s" % git_version.verbose_msg
)
parser.add_option("-b", "--base", metavar = "base", default = "plotburst_", help = "set the prefix for output filenames (default = plotburst_)")
parser.add_option("-f", "--format", metavar = "format", default = "png", help = "set the output image format (default = png)")
parser.add_option("-i", "--input-cache", metavar = "filename", default = None, help = "get file names from this LAL cache file")
parser.add_option("--plot", metavar = "number", action = "append", default = None, help = "only generate the given plot number")
parser.add_option("--frequency-range", metavar = "low,high", default = "0,2200", help = "set the peak frequency range for the plots (default = 0,2200)")
parser.add_option("--livetime-program", metavar = "name", default = "lalapps_power", help = "set the name of the program whose search_summary rows will define the livetime (default = \"lalapps_power\").")
parser.add_option("-v", "--verbose", action = "store_true", help = "be verbose")
options, filenames = parser.parse_args()
if options.input_cache:
filenames.extend([c.path for c in map(CacheEntry, file(options.input_cache))])
if options.plot:
options.plot = map(int, options.plot)
else:
options.plot = range(8)
options.frequency_range = map(float, options.frequency_range.split(","))
return options, (filenames or [None])
options, filenames = parse_command_line()
#
# =============================================================================
#
# Input
#
# =============================================================================
#
class Summary(object):
def __init__(self):
self.nevents = 0
self.start_time = []
self.duration = []
self.peak_time = []
self.peak_freq = []
self.bandwidth = []
self.lo_freq = []
self.snr = []
self.confidence = []
def snglburst_append(self, row):
global summary
if row.ifo not in summary:
summary[row.ifo] = Summary()
summary[row.ifo].nevents += 1
summary[row.ifo].start_time.append(lsctables.LIGOTimeGPS(row.start_time, row.start_time_ns))
summary[row.ifo].duration.append(row.duration)
summary[row.ifo].peak_time.append(lsctables.LIGOTimeGPS(row.peak_time, row.peak_time_ns))
summary[row.ifo].peak_freq.append(row.peak_frequency)
summary[row.ifo].bandwidth.append(row.bandwidth)
summary[row.ifo].lo_freq.append(row.central_freq - row.bandwidth / 2.0)
summary[row.ifo].snr.append(row.snr)
summary[row.ifo].confidence.append(row.confidence)
lsctables.SnglBurstTable.append = snglburst_append
#
# =============================================================================
#
# Confidence vs. Time
#
# =============================================================================
#
class ConfidenceVsTime(object):
def __init__(self, ifo):
self.fig, self.axes = SnglBurstUtils.make_burst_plot("GPS Time (s)", "Confidence")
self.ifo = ifo
self.nevents = 0
self.x = []
self.y = []
self.seglist = segments.segmentlist()
self.axes.semilogy()
def add_contents(self, contents, seglists):
self.nevents += contents[self.ifo].nevents
self.x.extend(map(float, contents[self.ifo].peak_time))
self.y.extend(contents[self.ifo].confidence)
self.seglist |= seglists[self.ifo]
def finish(self):
self.axes.set_title("Trigger Confidence vs. Time\n(%d Triggers)" % self.nevents)
self.axes.plot(self.x, self.y, "k+")
for seg in ~self.seglist & segments.segmentlist([segments.segment(self.axes.get_xlim())]):
self.axes.axvspan(float(seg[0]), float(seg[1]), facecolor = "k", alpha = 0.2)
#
# =============================================================================
#
# Confidence vs. Peak Frequency
#
# =============================================================================
#
class ConfidenceVsFrequencyScatter(object):
def __init__(self, ifo):
self.fig, self.axes = SnglBurstUtils.make_burst_plot("Peak Frequency (Hz)", "Confidence")
self.ifo = ifo
self.nevents = 0
self.x = []
self.y = []
self.axes.semilogy()
def add_contents(self, contents, seglists):
self.nevents += contents[self.ifo].nevents
self.x.extend(contents[self.ifo].peak_freq)
self.y.extend(contents[self.ifo].confidence)
def finish(self):
self.axes.set_title("Trigger Confidence vs. Peak Frequency\n(%d Triggers)" % self.nevents)
self.axes.plot(self.x, self.y, "k+")
self.axes.set_xlim((min(self.x), max(self.x)))
#
# =============================================================================
#
# Rate vs. Peak Frequency
#
# =============================================================================
#
class RateVsPeakFreq(object):
def __init__(self, ifo, interval, width):
self.fig, self.axes = SnglBurstUtils.make_burst_plot("Peak Frequency (Hz)", "Trigger Rate Spectral Density (triggers / s / Hz)")
self.ifo = ifo
self.nevents = 0
# 21 bins per filter width
bins = int(float(abs(interval)) / width) * 21
binning = rate.NDBins((rate.LinearBins(interval[0], interval[1], bins),))
self.rate = rate.BinnedDensity(binning)
def add_contents(self, contents, seglists):
self.nevents += contents[self.ifo].nevents
for f in contents[self.ifo].peak_freq:
try:
self.rate.count[f,] += 1.0
except IndexError:
raise ValueError("trigger peak frequency %g Hz outside plot range [%g Hz, %g Hz]" % (f, self.rate.bins[0].min, self.rate.bins[0].max))
def finish(self):
self.axes.set_title("Trigger Rate vs. Peak Frequency\n(%d Triggers)" % self.nevents)
# 21 bins per filter width
rate.filter_array(self.rate.array, rate.gaussian_window(21))
xvals = self.rate.centres()[0]
self.axes.plot(xvals, self.rate.at_centres(), "k")
self.axes.semilogy()
self.axes.set_xlim((min(xvals), max(xvals)))
#
# =============================================================================
#
# Trigger Duration Histogram
#
# =============================================================================
#
class Durations(object):
def __init__(self, ifo):
self.fig, self.axes = SnglBurstUtils.make_burst_plot("Duration (s)", "Trigger Count")
self.ifo = ifo
self.nevents = 0
self.bins = {}
def add_contents(self, contents, seglists):
self.nevents += contents[self.ifo].nevents
for dt in contents[self.ifo].duration:
if dt not in self.bins:
self.bins[dt] = 0
self.bins[dt] += 1
def finish(self):
self.axes.set_title("Trigger Durations\n(%d Triggers)" % self.nevents)
data = self.bins.items()
data.sort()
self.axes.plot([d[0] for d in data], [d[1] for d in data], "ko-")
#
# =============================================================================
#
# Time Between Triggers Histogram