Commit f068099f authored by Kipp Cannon's avatar Kipp Cannon
Browse files

Merge branch 'burstupdate' into 'master'

Burst Update

See merge request lscsoft/lalsuite!1405
parents 23a9d225 0ba57da0
......@@ -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
......
......@@ -34,6 +34,7 @@ import sys
from PIL import Image
from glue.text_progress_bar import ProgressBar
from ligo.lw import ligolw
from ligo.lw import lsctables
from ligo.lw import utils as ligolw_utils
......@@ -49,7 +50,15 @@ class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
pass
lsctables.SimBurst = lsctables.SimBurstTable.RowType = lalmetaio.SimBurst
class SimBurst(lalmetaio.SimBurst):
def __init__(self, **kwargs):
super(SimBurst, self).__init__()
for key, value in kwargs.items():
setattr(self, key, value)
@property
def time_geocent_gps_ns(self):
return self.time_geocent_gps.gpsNanoSeconds
lsctables.SimBurst = lsctables.SimBurstTable.RowType = SimBurst
__author__ = "Kipp Cannon <kipp.cannon@ligo.org>"
......@@ -72,31 +81,38 @@ def parse_command_line():
usage = "%prog [options] filename",
description = "Convert an image into a LIGO Light-Weight XML file containing a list of sine-Gaussian burst injections. When injected into data, the injections will cause a waterfall plot to display the image."
)
parser.add_option("-l", "--f-low", metavar = "Hz", type = "float", default = 64.0, help = "Set the low-frequency limit of the tiling (default = 64).")
parser.add_option("-d", "--delta-f", metavar = "Hz", type = "float", default = 16.0, help = "Set the frequency spacing of the tiling (default = 16).")
parser.add_option("-t", "--delta-t", metavar = "s", type = "float", default = 1.0 / 16, help = "Set the time spacing of the tiling (default = 1/16).")
parser.add_option("-H", "--height", metavar = "pixels", type = "int", default = 64, help = "Set the number of tiles in the frequency domain (default = 64).")
parser.add_option("-l", "--f-low", metavar = "Hz", type = "float", default = 64.0, help = "Set the low-frequency limit of the tiling in Hertz (default = 64).")
parser.add_option("-d", "--delta-f", metavar = "Hz", type = "float", default = 16.0, help = "Set the bandwidth of the pixels in Hertz (default = 16). Must be > 0. The product of --delta-f and --delta-t must be at least 2/pi.")
parser.add_option("-t", "--delta-t", metavar = "s", type = "float", default = 1.0 / 16, help = "Set the duration of the pixels in seconds (default = 1/16). Must be > 0. The product of --delta-f and --delta-t must be at least 2/pi.")
parser.add_option("-f", "--overlap-fraction", metavar = "fraction", type = "float", default = 0.25, help = "Set the fraction by which adjacent tiles overlap (default = 0.25). The pixels centres are spaced (1 - --overlap-fraction) * --delta-f apart in the frequency domain. The value must be in the range [0, 1).")
parser.add_option("-H", "--height", metavar = "pixels", type = "int", default = 64, help = "Set the number of pixles in the frequency domain (default = 64). The image will be scaled to this vertical size, and the number of pixels in the time domain (horizontal size) will be fixed by the image aspect ratio.")
parser.add_option("-o", "--output", metavar = "filename", help = "Set the name of the output file (default = stdout).")
parser.add_option("-s", "--gps-start-time", metavar = "seconds", help = "Set the start time of the tiling in GPS seconds (required).")
parser.add_option("-f", "--overlap-fraction", metavar = "fraction", type = "float", default = 0.25, help = "Set the fraction by which adjacent tiles overlap (default = 0.25). The value must be in the range [0, 1).")
parser.add_option("--sample-rate", metavar = "Hz", type = "int", default = 16384, help = "Set the sample rate of the data into which the injections will be placed (default = 16384). This information is required in order to normalize each pixel accurately. If the wrong value is used, the result will be the addition of noise to the image.")
parser.add_option("-n", "--hrss-scale", metavar = "hrss", type = "float", default = 1e-20, help = "Set the hrss corresponding to \"white\" (default = 1e-20).")
parser.add_option("--sample-rate", metavar = "Hz", type = "int", default = 16384, help = "Set the sample rate in Hertz of the data into which the injections will be placed (default = 16384). This information is required in order to normalize each pixel accurately. If the wrong value is used, the result will be the addition of noise to the image. The highest frequency pixel must have a centre frequency < 1/2 this frequency.")
parser.add_option("-n", "--hrss-scale", metavar = "hrss", type = "float", default = 1e-20, help = "Set the single-pixel hrss (root-sum-square strain) corresponding to \"white\" (default = 1e-20).")
parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")
parser.add_option("-T", "--time-slide-xml", metavar = "filename", help = "Associate injections with the first time slide ID in this XML file (required).")
options, filenames = parser.parse_args()
# save for the process_params table
options.options_dict = dict(options.__dict__)
if options.gps_start_time is None:
raise ValueError("missing required option --gps-start-time")
if not (0 <= options.overlap_fraction < 1.0):
raise ValueError("--overlap-fraction must be in [0, 1)")
if options.delta_f * options.delta_t < 2 / math.pi:
if not 0. < options.delta_f:
raise ValueError("--delta-f must be > 0")
if not 0. < options.delta_t:
raise ValueError("--delta-t must be > 0")
if options.delta_f * options.delta_t < 2. / math.pi:
raise ValueError("the product of --delta-t and --delta-f must be >= 2/pi")
if not (0. <= options.overlap_fraction < 1.):
raise ValueError("--overlap-fraction must be in [0, 1)")
if options.f_low + options.height * (1. - options.overlap_fraction) * options.delta_f >= options.sample_rate / 2.:
raise ValueError("total waveform bandwidth exceeds Nyquist frequency")
if options.time_slide_xml is None:
raise ValueError("missing required option time-slide-xml")
# for the process_params table
options.options_dict = dict(options.__dict__)
# type-cast
options.gps_start_time = lsctables.LIGOTimeGPS(options.gps_start_time)
return options, filenames
......@@ -114,71 +130,133 @@ def parse_command_line():
options, filenames = parse_command_line()
if options.verbose:
print("time-frequency tiles have %g degrees of freedom" % (2 * options.delta_t * options.delta_f), file=sys.stderr)
#
# use the time-slide file to start the output document
#
xmldoc = ligolw_utils.load_filename(options.time_slide_xml, verbose = options.verbose, contenthandler = LIGOLWContentHandler)
#
# add our metadata
#
xmldoc = ligolw.Document()
xmldoc.appendChild(ligolw.LIGO_LW())
process = ligolw_process.register_to_xmldoc(xmldoc, u"lalapps_binj_pic", options.options_dict, version = git_version.version)
time_slide_table = xmldoc.childNodes[-1].appendChild(lsctables.TimeSlideTable.get_table(ligolw_utils.load_filename(options.time_slide_xml, verbose = options.verbose, contenthandler = LIGOLWContentHandler)))
#
# use whatever time slide vector comes first in the table (lazy)
#
time_slide_table = lsctables.TimeSlideTable.get_table(xmldoc)
time_slide_id = time_slide_table[0].time_slide_id
if options.verbose:
print("associating injections with time slide ID \"%s\": %s" % (time_slide_id, time_slide_table.as_dict()[time_slide_id]), file=sys.stderr)
for row in time_slide_table:
row.process_id = process.process_id
sim_burst_tbl = xmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.SimBurstTable, ["process_id", "simulation_id", "time_slide_id", "waveform", "waveform_number", "ra", "dec", "psi", "time_geocent_gps", "time_geocent_gps_ns", "duration", "frequency", "bandwidth", "egw_over_rsquared", "pol_ellipse_angle", "pol_ellipse_e"]))
print("associating injections with time slide (%d) %s" % (time_slide_id, time_slide_table.as_dict()[time_slide_id]), file = sys.stderr)
#
# find or add a sim_burst table
#
try:
lsctables.SimBurstTable.get_table(xmldoc)
except ValueError:
# no sim_burst table in document
pass
else:
raise ValueError("%s contains a sim_burst table. this program isn't smart enough to deal with that." % options.time_slide_xml)
sim_burst_tbl = xmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.SimBurstTable, ["process:process_id", "simulation_id", "time_slide:time_slide_id", "waveform", "waveform_number", "ra", "dec", "psi", "time_geocent_gps", "time_geocent_gps_ns", "duration", "frequency", "bandwidth", "egw_over_rsquared", "pol_ellipse_angle", "pol_ellipse_e"]))
#
# populate the sim_burst table with pixels
#
for filename in filenames:
if options.verbose:
print("time-frequency tiles have %g degrees of freedom" % (2. * options.delta_t * options.delta_f), file = sys.stderr)
for n, filename in enumerate(filenames, 1):
if options.verbose:
print("loading %s ..." % filename, file=sys.stderr)
print("%d/%d: loading %s ..." % (n, len(filenames), filename), file = sys.stderr)
img = Image.open(filename)
width, height = img.size
width, height = int(round(width / float(height) * options.height)), options.height
width, height = int(round(width * options.height / float(height))), options.height
if options.verbose:
print("converting to %dx%d grayscale ... " % (width, height), file=sys.stderr)
print("converting to %dx%d grayscale ... " % (width, height), file = sys.stderr)
img = img.resize((width, height)).convert("L")
progress = ProgressBar("computing pixels", max = width * height) if options.verbose else None
for i in range(width):
for j in range(height):
# new row
row = lsctables.SimBurst()
row.process_id = process.process_id
row.simulation_id = sim_burst_tbl.get_next_id()
row.time_slide_id = time_slide_id
# source orientation
row.ra = row.dec = row.psi = 0
if progress is not None:
progress.increment()
# amplitude. hrss column is ignored by waveform
# generation code. it is included for convenience,
# to record the desired pixel brightness. because
# band- and time-limited white-noise burst
row.waveform = "BTLWNB"
row.waveform_number = 0
row.pol_ellipse_e = row.pol_ellipse_angle = 0
# waveforms are random, the waveform's final
# ampltiude (in the egw_over_rsquared column) is
# determined by generating the burst at a canonical
# amplitude, measuring its hrss, then rescaling to
# achieve the desired value. this process requires
# the final sample rate to be known.
hrss = options.hrss_scale * img.getpixel((i, height - 1 - j)) / 255.0
if hrss == 0.:
# don't generate injections for black
# pixels
continue
# create and initialize table row
row = lsctables.SimBurst(
# metadata
process_id = process.process_id,
simulation_id = sim_burst_tbl.get_next_id(),
time_slide_id = time_slide_id,
# source orientation
ra = 0.,
dec = 0.,
psi = 0.,
# band- and time-limited white-noise burst
waveform = "BTLWNB",
waveform_number = 0,
pol_ellipse_e = 0.,
pol_ellipse_angle = 0.,
# time-frequency co-ordinates
time_geocent_gps = options.gps_start_time + i * options.delta_t * (1.0 - options.overlap_fraction),
frequency = options.f_low + j * options.delta_f * (1.0 - options.overlap_fraction),
bandwidth = options.delta_f,
duration = options.delta_t,
# brightness. hrss is ignored, set
# egw_over_rsquared to 1, to be re-scaled
# after measuring waveform's real
# brightness
hrss = hrss,
egw_over_rsquared = 1.
)
# generate waveform, then scale egw_over_rsquared
# to get desired brightness
row.egw_over_rsquared *= hrss / lalsimulation.MeasureHrss(*lalburst.GenerateSimBurst(row, 1.0 / options.sample_rate))
# put row into table
sim_burst_tbl.append(row)
del progress
# time-frequency co-ordinates
row.frequency = options.f_low + j * options.delta_f * (1.0 - options.overlap_fraction)
gps = options.gps_start_time + i * options.delta_t * (1.0 - options.overlap_fraction)
row.time_geocent_gps, row.time_geocent_gps_ns = gps.seconds, gps.nanoseconds
row.bandwidth = options.delta_f
row.duration = options.delta_t
# amplitude. hrss column is ignored by waveform
# generation code. it is included for convenience.
# the waveform is normalized by generating the
# burst, measuring its hrss, then rescaling to
# achieve the desired value. this process requires
# the sample rate to be known.
row.hrss = options.hrss_scale * img.getpixel((i, height - 1 - j)) / 255.0
if row.hrss != 0:
row.egw_over_rsquared = 1
row.egw_over_rsquared *= row.hrss / lalsimulation.MeasureHrss(*lalburst.GenerateSimBurst(row, 1.0 / 16384))
sim_burst_tbl.append(row)
if options.verbose:
print("generating sim_burst table ... %d injections\r" % len(sim_burst_tbl), end=' ', file=sys.stderr)
if options.verbose:
print(file=sys.stderr)
#
# write output
#
ligolw_utils.write_filename(xmldoc, options.output, gz = (options.output or "stdout").endswith(".gz"), verbose = options.verbose)
......@@ -82,23 +82,6 @@ def parse_command_line():
return options, (filenames or [None])
#
# =============================================================================
#
# Input
#
# =============================================================================
#
#
# Use interning row builder to save memory.
#
binjfind.lsctables.table.TableStream.RowBuilder = binjfind.lsctables.table.InterningRowBuilder
#
# =============================================================================
#
......@@ -151,7 +134,6 @@ for n, filename in enumerate(filenames):
if options.verbose:
print("%d/%d:" % (n + 1, len(filenames)), end=' ', file=sys.stderr)
xmldoc = ligolw_utils.load_filename(filename, verbose = options.verbose, contenthandler = ligolw.LIGOLWContentHandler)
binjfind.lsctables.table.InterningRowBuilder.strings.clear()
#
# have we already procesed it?
......
......@@ -129,17 +129,9 @@ def parse_command_line():
#
@lsctables.use_in
class ContentHandler(ligolw.LIGOLWContentHandler):
pass
lsctables.use_in(ContentHandler)
#
# Use interning row builder to save memory.
#
lsctables.table.TableStream.RowBuilder = lsctables.table.InterningRowBuilder
#
......@@ -505,7 +497,6 @@ else:
for filename in filenames:
xmldoc = ligolw_utils.load_filename(filename, verbose = options.verbose, contenthandler = ContentHandler)
lsctables.table.InterningRowBuilder.strings.clear()
xmldoc = ligolw_bucut(xmldoc, options, keep_this_sngl_burst, veto_segments = veto_segments, del_non_coincs = options.coinc_only, del_skipped_injections = options.inj_made_only, program = options.program, verbose = options.verbose)
ligolw_utils.write_filename(xmldoc, filename, verbose = options.verbose, gz = (filename or "stdout").endswith(".gz"))
xmldoc.unlink()
#
# Copyright (C) 2006 Kipp C. Cannon
# Copyright (C) 2005--2010,2013,2015,2016,2019 Kipp Cannon
# Copyright (C) 2004--2006 Saikat Ray-Majumder
# Copyright (C) 2003--2005 Duncan Brown
#
# 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
......@@ -121,7 +123,7 @@ def parse_config_file(options):
seglistdict = segments.segmentlistdict()
tiling_phase = {}
for ifo in config.get("pipeline", "ifos").split():
seglistdict[ifo] = segmentsUtils.fromsegwizard(file(config.get("pipeline", "seglist_%s" % ifo)), coltype = LIGOTimeGPS).coalesce()
seglistdict[ifo] = segmentsUtils.fromsegwizard(open(config.get("pipeline", "seglist_%s" % ifo)), coltype = LIGOTimeGPS).coalesce()
try:
offset = config.getfloat("pipeline", "tiling_phase_%s" % ifo)
except NoOptionError:
......
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()