From 8b24b66b8e5b828b49d8770fbb5b44f08ece1b03 Mon Sep 17 00:00:00 2001 From: Kipp Cannon Date: Fri, 15 Nov 2019 13:51:24 +0900 Subject: [PATCH 1/5] lalapps_power_pipe: file() --> open() - for python3 --- lalapps/src/power/lalapps_power_pipe.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/lalapps/src/power/lalapps_power_pipe.py b/lalapps/src/power/lalapps_power_pipe.py index 92c9d6e511..4f263ba60e 100644 --- a/lalapps/src/power/lalapps_power_pipe.py +++ b/lalapps/src/power/lalapps_power_pipe.py @@ -1,5 +1,7 @@ # -# 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: -- GitLab From 50bc561e94cfbc4b61c641ec48ac722d820893d9 Mon Sep 17 00:00:00 2001 From: Kipp Cannon Date: Thu, 29 Aug 2019 18:39:10 +0900 Subject: [PATCH 2/5] lalapps_binj_pic: rewrite - this was badly bit-rotten, this patch gets it working again --- lalapps/src/power/lalapps_binj_pic.py | 198 ++++++++++++++++++-------- 1 file changed, 138 insertions(+), 60 deletions(-) diff --git a/lalapps/src/power/lalapps_binj_pic.py b/lalapps/src/power/lalapps_binj_pic.py index 66aeab918f..492b2f4dd6 100644 --- a/lalapps/src/power/lalapps_binj_pic.py +++ b/lalapps/src/power/lalapps_binj_pic.py @@ -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 " @@ -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) -- GitLab From 7192f73934ee4785d8ed051120bf368c1a2ad714 Mon Sep 17 00:00:00 2001 From: Kipp Cannon Date: Thu, 29 Aug 2019 12:34:44 +0900 Subject: [PATCH 3/5] lalapps: delete obsolete power programs - lalapps_power_plot_burca - lalapps_power_plot_burst - lalapps_power_veto - lalapps_power_plot_burca2 --- lalapps/.gitignore | 4 - lalapps/conda/meta.yaml.in.in | 4 - lalapps/src/power/Makefile.am | 6 +- lalapps/src/power/lalapps_power_plot_burca.py | 620 ------------------ .../src/power/lalapps_power_plot_burca2.py | 282 -------- lalapps/src/power/lalapps_power_plot_burst.py | 520 --------------- lalapps/src/power/lalapps_power_veto.py | 336 ---------- 7 files changed, 1 insertion(+), 1771 deletions(-) delete mode 100644 lalapps/src/power/lalapps_power_plot_burca.py delete mode 100644 lalapps/src/power/lalapps_power_plot_burca2.py delete mode 100644 lalapps/src/power/lalapps_power_plot_burst.py delete mode 100644 lalapps/src/power/lalapps_power_veto.py diff --git a/lalapps/.gitignore b/lalapps/.gitignore index 9d3090109e..4fc02b8571 100644 --- a/lalapps/.gitignore +++ b/lalapps/.gitignore @@ -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 diff --git a/lalapps/conda/meta.yaml.in.in b/lalapps/conda/meta.yaml.in.in index fbe20047f8..a822c9c8ac 100644 --- a/lalapps/conda/meta.yaml.in.in +++ b/lalapps/conda/meta.yaml.in.in @@ -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 diff --git a/lalapps/src/power/Makefile.am b/lalapps/src/power/Makefile.am index 2eab4b92b1..62d58bd954 100644 --- a/lalapps/src/power/Makefile.am +++ b/lalapps/src/power/Makefile.am @@ -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 diff --git a/lalapps/src/power/lalapps_power_plot_burca.py b/lalapps/src/power/lalapps_power_plot_burca.py deleted file mode 100644 index 275932c867..0000000000 --- a/lalapps/src/power/lalapps_power_plot_burca.py +++ /dev/null @@ -1,620 +0,0 @@ -# -# 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 -import matplotlib.cm -import numpy -from optparse import OptionParser -import sqlite3 -import sys - - -from ligo.lw import utils as ligolw_utils -import lal -from lal import rate -from lalburst import git_version -from lalburst import SnglBurstUtils -from lalburst import SimBurstUtils -from ligo import segments - - -__author__ = "Kipp Cannon " -__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 = "plotburca_", help = "Set the prefix for output filenames (default = \"plotburca_\").") - 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("--plot", metavar = "number", action = "append", default = None, help = "Generate the given plot number.") - parser.add_option("-s", "--skip", metavar = "number", default = 0, help = "Skip this many files at the start of the list (default = process all).") - parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.") - options, filenames = parser.parse_args() - - if options.plot: - options.plot = map(int, options.plot) - else: - options.plot = range(11) - - options.skip = int(options.skip) - - return options, (filenames or [None]) - - -# -# ============================================================================= -# -# Rate Contours -# -# ============================================================================= -# - - -class RateContours(object): - 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 add_contents(self, contents): - if self.tisi_rows is None: - # get a list of time slide dictionaries - self.tisi_rows = contents.time_slide_table.as_dict().values() - - # find the largest and smallest offsets - min_offset = min(offset for vector in self.tisi_rows for offset in vector.values()) - max_offset = max(offset for vector in self.tisi_rows for offset in vector.values()) - - # a guess at the time slide spacing: works if the - # time slides are distributed as a square grid over - # the plot area. (max - min)^2 gives the area of - # the time slide square in square seconds; dividing - # by the length of the time slide list gives the - # average area per time slide; taking the square - # root of that gives the average distance between - # adjacent time slides in seconds - time_slide_spacing = ((max_offset - min_offset)**2 / len(self.tisi_rows))**0.5 - - # use an average of 3 bins per time slide in each - # direction, but round to an odd integer - nbins = int(math.ceil((max_offset - min_offset) / time_slide_spacing * 3)) - - # construct the binning - self.counts = rate.BinnedArray(rate.NDBins((rate.LinearBins(min_offset, max_offset, nbins), rate.LinearBins(min_offset, max_offset, nbins)))) - - self.seglists |= contents.seglists - - for offsets in contents.connection.cursor().execute(""" -SELECT tx.offset, ty.offset FROM - coinc_event - JOIN time_slide AS tx ON ( - tx.time_slide_id == coinc_event.time_slide_id - ) - JOIN time_slide AS ty ON ( - ty.time_slide_id == coinc_event.time_slide_id - ) -WHERE - coinc_event.coinc_def_id == ? - AND tx.instrument == ? - AND ty.instrument == ? - """, (contents.bb_definer_id, self.x_instrument, self.y_instrument)): - try: - self.counts[offsets] += 1 - except IndexError: - # beyond plot boundaries - pass - - def finish(self): - livetime = rate.BinnedArray(self.counts.bins) - for offsets in self.tisi_rows: - self.seglists.offsets.update(offsets) - livetime[offsets[self.x_instrument], offsets[self.y_instrument]] += float(abs(self.seglists.intersection(self.seglists.keys()))) - zvals = self.counts.at_centres() / livetime.at_centres() - rate.filter_array(zvals, rate.gaussian_window(3, 3)) - xcoords, ycoords = self.counts.centres() - self.axes.contour(xcoords, ycoords, numpy.transpose(numpy.log(zvals))) - for offsets in self.tisi_rows: - if any(offsets): - # time slide vector is non-zero-lag - self.axes.plot((offsets[self.x_instrument],), (offsets[self.y_instrument],), "k+") - else: - # time slide vector is zero-lag - self.axes.plot((offsets[self.x_instrument],), (offsets[self.y_instrument],), "r+") - - self.axes.set_xlim([self.counts.bins().min[0], self.counts.bins().max[0]]) - self.axes.set_ylim([self.counts.bins().min[1], self.counts.bins().max[1]]) - self.axes.set_title(r"Coincident Event Rate vs. Instrument Time Offset (Logarithmic Rate Contours)") - - -# -# ============================================================================= -# -# Confidence Scatter -# -# ============================================================================= -# - -def incomplete_injection_coincs(contents): - # determine burst <--> burst coincidences for which at least one - # burst, *but not all*, was identified as an injection; these are - # places in the data where an injection was done, a coincident - # event was seen, but where, later, the injection was not found to - # match all events in the coincidence; these perhaps indicate - # power leaking from the injection into nearby tiles, or accidental - # coincidence of a marginal injection with near-by noise, etc, and - # so although they aren't "bang-on" reconstructions of injections - # they are nevertheless injections that are found and survive a - # coincidence cut. - for values in contents.connection.cursor().execute(""" -SELECT - bb_coinc_event.* -FROM - coinc_event AS bb_coinc_event -WHERE - bb_coinc_event.coinc_def_id == ? - AND EXISTS ( - SELECT - * - FROM - coinc_event AS sb_coinc_event - JOIN coinc_event_map AS a ON ( - a.coinc_event_id == sb_coinc_event.coinc_event_id - AND a.table_name == 'sngl_burst' - ) - JOIN coinc_event_map AS b ON ( - b.table_name == 'sngl_burst' - AND b.event_id == a.event_id - ) - WHERE - sb_coinc_event.coinc_def_id == ? - AND b.coinc_event_id == bb_coinc_event.coinc_event_id - AND sb_coinc_event.time_slide_id == bb_coinc_event.time_slide_id - ) -EXCEPT - SELECT - c.event_id - FROM - coinc_event_map AS c - JOIN coinc_event AS sc_coinc_event ON ( - c.coinc_event_id == sc_coinc_event.coinc_event_id - AND c.table_name == 'coinc_event' - ) - WHERE - sc_coinc_event.coinc_def_id == ? - AND sc_coinc_event.time_slide_id == bb_coinc_event.time_slide_id - """, (contents.bb_definer_id, contents.sb_definer_id, contents.scn_definer_id)): - yield contents.coinc_table.row_from_cols(values) - - -def magnitude_a(burst): - return burst.ms_hrss - - -def magnitude_b(burst): - gmst = lal.GreenwichMeanSiderealTime(burst.peak) - fplus, fcross = lal.ComputeDetAMResponse(lal.cached_detector_by_prefix[burst.ifo].response, SimBurstUtils.MW_CENTER_J2000_RA_RAD, SimBurstUtils.MW_CENTER_J2000_DEC_RAD, 0.0, gmst) - return (burst.ms_hrss**2.0 / (fplus**2.0 + fcross**2.0))**0.5 - - -class ConfidenceContours(object): - def __init__(self, x_instrument, y_instrument, magnitude, desc, min_magnitude, max_magnitude): - self.fig, self.axes = SnglBurstUtils.make_burst_plot("%s %s" % (x_instrument, desc), "%s %s" % (y_instrument, desc)) - self.fig.set_size_inches(6,6) - self.axes.loglog() - self.x_instrument = x_instrument - self.y_instrument = y_instrument - self.magnitude = magnitude - self.foreground_x = [] - self.foreground_y = [] - self.n_foreground = 0 - self.n_background = 0 - self.n_injections = 0 - self.foreground_bins = rate.BinnedArray(rate.NDBins((rate.LogarithmicBins(min_magnitude, max_magnitude, 1024), rate.LogarithmicBins(min_magnitude, max_magnitude, 1024)))) - self.background_bins = rate.BinnedArray(rate.NDBins((rate.LogarithmicBins(min_magnitude, max_magnitude, 1024), rate.LogarithmicBins(min_magnitude, max_magnitude, 1024)))) - self.coinc_injection_bins = rate.BinnedArray(rate.NDBins((rate.LogarithmicBins(min_magnitude, max_magnitude, 1024), rate.LogarithmicBins(min_magnitude, max_magnitude, 1024)))) - self.incomplete_coinc_injection_bins = rate.BinnedArray(rate.NDBins((rate.LogarithmicBins(min_magnitude, max_magnitude, 1024), rate.LogarithmicBins(min_magnitude, max_magnitude, 1024)))) - - 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 add_contents(self, contents): - if self.x_instrument not in contents.instruments or self.y_instrument not in contents.instruments: - # fast path: these are not the instruments we're - # looking for - return - for values in contents.connection.cursor().execute(""" -SELECT - coinc_event.*, - EXISTS ( - SELECT - * - FROM - time_slide - WHERE - time_slide.time_slide_id == coinc_event.time_slide_id - AND time_slide.offset != 0 - ) -FROM - coinc_event -WHERE - coinc_event.coinc_def_id == ? - """, (contents.bb_definer_id,)): - coinc = contents.coinc_table.row_from_cols(values) - not_zero_lag = values[-1] - if not_zero_lag: - self.n_foreground += 1 - x, y = self.coords(contents, coinc.coinc_event_id) - try: - self.foreground_bins[x, y] += 1 - #self.foreground_x.append(x) - #self.foreground_y.append(y) - except IndexError: - # not on plot axes - pass - else: - self.n_background += 1 - try: - self.background_bins[self.coords(contents, coinc.coinc_event_id)] += 1 - except IndexError: - # not on plot axes - pass - if contents.sce_definer_id is not None: - # this query assumes each injection can match at - # most 1 coinc - for values in contents.connection.cursor().execute(""" -SELECT - coinc_event.* -FROM - coinc_event AS sim_coinc - JOIN coinc_event_map ON ( - coinc_event_map.coinc_event_id == sim_coinc.coinc_event_id - ) - JOIN coinc_event ON ( - coinc_event_map.table_name == 'coinc_event' - AND coinc_event_map.event_id == coinc_event.coinc_event_id - ) -WHERE - sim_coinc.coinc_def_id == ? - """, (contents.sce_definer_id,)): - coinc = contents.coinc_event_table.row_from_cols(values) - self.n_injections += 1 - try: - self.coinc_injection_bins[self.coords(contents, coinc.coinc_event_id)] += 1 - except IndexError: - # not on plot axes - pass - for coinc in incomplete_injection_coincs(contents): - try: - self.incomplete_coinc_injection_bins[self.coords(contents, coinc.coinc_event_id)] += 1 - except IndexError: - # not on plot axes - pass - - def finish(self): - self.axes.set_title(r"\begin{center}Distribution of Coincident Events (%d Foreground, %d Background Events, %d Injections Found in Coincidence, Logarithmic Density Contours)\end{center}" % (self.n_foreground, self.n_background, self.n_injections)) - xcoords, ycoords = self.background_bins.centres() - - # prepare the data - rate.filter_array(self.foreground_bins.array, rate.gaussian_window(8, 8)) - rate.filter_array(self.background_bins.array, rate.gaussian_window(8, 8)) - rate.filter_array(self.coinc_injection_bins.array, rate.gaussian_window(8, 8)) - rate.filter_array(self.incomplete_coinc_injection_bins.array, rate.gaussian_window(8, 8)) - self.foreground_bins.logregularize() - self.background_bins.logregularize() - self.coinc_injection_bins.logregularize() - self.incomplete_coinc_injection_bins.logregularize() - - # plot background contours - max_density = math.log(self.background_bins.array.max()) - self.axes.contour(xcoords, ycoords, numpy.transpose(numpy.log(self.background_bins.array)), sorted(max_density - n for n in xrange(0, 10, 1)), cmap = matplotlib.cm.Greys) - - # plot foreground (zero-lag) contours - max_density = math.log(self.foreground_bins.array.max()) - self.axes.contour(xcoords, ycoords, numpy.transpose(numpy.log(self.foreground_bins.array)), sorted(max_density - n for n in xrange(0, 10, 1)), cmap = matplotlib.cm.Reds) - #self.axes.plot(self.foreground_x, self.foreground_y, "r+") - - # plot coincident injection contours - max_density = math.log(self.coinc_injection_bins.array.max()) - self.axes.contour(xcoords, ycoords, numpy.transpose(numpy.log(self.coinc_injection_bins.array)), sorted(max_density - n for n in xrange(0, 10, 1)), cmap = matplotlib.cm.Blues) - - # plot incomplete coincident injection contours - max_density = math.log(self.incomplete_coinc_injection_bins.array.max()) - self.axes.contour(xcoords, ycoords, numpy.transpose(numpy.log(self.incomplete_coinc_injection_bins.array)), sorted(max_density - n for n in xrange(0, 10, 1)), cmap = matplotlib.cm.Greens) - - # add diagonal line - lower = max(binning.min for binning in self.background_bins.bins) - upper = min(binning.max for binning in self.background_bins.bins) - self.axes.plot([lower, upper], [lower, upper], "k:") - - # fix axes limits - self.axes.set_xlim([self.background_bins.bins[0].min, self.background_bins.bins[0].max]) - self.axes.set_ylim([self.background_bins.bins[1].min, self.background_bins.bins[1].max]) - - -class ConfidenceContourProjection(object): - def __init__(self, x, y, magnitude, max_magnitude): - self.fig, self.axes = SnglBurstUtils.make_burst_plot("X", "Y") - self.fig.set_size_inches(6,6) - self.x = x - self.y = y - self.magnitude = magnitude - self.n_foreground = 0 - self.n_background = 0 - self.n_injections = 0 - max_magnitude = math.log10(max_magnitude) - self.foreground_bins = rate.BinnedArray(rate.NDBins((rate.LinearBins(-max_magnitude, max_magnitude, 1024), rate.LinearBins(-max_magnitude, max_magnitude, 1024)))) - self.background_bins = rate.BinnedArray(rate.NDBins((rate.LinearBins(-max_magnitude, max_magnitude, 1024), rate.LinearBins(-max_magnitude, max_magnitude, 1024)))) - self.coinc_injection_bins = rate.BinnedArray(rate.NDBins((rate.LinearBins(-max_magnitude, max_magnitude, 1024), rate.LinearBins(-max_magnitude, max_magnitude, 1024)))) - self.incomplete_coinc_injection_bins = rate.BinnedArray(rate.NDBins((rate.LinearBins(-max_magnitude, max_magnitude, 1024), rate.LinearBins(-max_magnitude, max_magnitude, 1024)))) - - def coords(self, contents, coinc_event_id): - mag = numpy.zeros(3, "Float64") - for burst in SnglBurstUtils.coinc_sngl_bursts(contents, coinc_event_id): - if burst.ifo == "H1": - mag[0] = math.log10(self.magnitude(burst)) - elif burst.ifo == "H2": - mag[1] = math.log10(self.magnitude(burst)) - elif burst.ifo == "L1": - mag[2] = math.log10(self.magnitude(burst)) - return numpy.inner(mag, self.x), numpy.inner(mag, self.y) - - def add_contents(self, contents): - for values in contents.connection.cursor().execute(""" -SELECT - coinc_event.*, - EXISTS ( - SELECT - * - FROM - time_slide - WHERE - time_slide.time_slide_id == coinc_event.time_slide_id - AND time_slide.offset != 0 - ) -FROM - coinc_event -WHERE - coinc_event.coinc_def_id == ? - """, (contents.bb_definer_id,)): - coinc = contents.coinc_table.row_from_cols(values) - not_zero_lag = values[-1] - x, y = self.coords(contents, coinc.coinc_event_id) - if not_zero_lag: - self.n_background += 1 - try: - self.background_bins[x, y] += 1 - except IndexError: - # not on plot axes - pass - else: - self.n_foreground += 1 - try: - self.foreground_bins[x, y] += 1 - #self.foreground_x.append(x) - #self.foreground_y.append(y) - except IndexError: - # not on plot axes - pass - if contents.sce_definer_id is not None: - # this query assumes each injection can match at - # most 1 coinc - for values in contents.connection.cursor().execute(""" -SELECT - coinc_event.* -FROM - coinc_event AS sim_coinc - JOIN coinc_event_map ON ( - coinc_event_map.coinc_event_id == sim_coinc.coinc_event_id - ) - JOIN coinc_event ON ( - coinc_event_map.table_name == 'coinc_event' - AND coinc_event_map.event_id == coinc_event.coinc_event_id - ) -WHERE - sim_coinc.coinc_def_id == ? - """, (contents.sce_definer_id,)): - coinc = contents.coinc_event_table.row_from_cols(values) - self.n_injections += 1 - try: - self.coinc_injection_bins[self.coords(contents, coinc.coinc_event_id)] += 1 - except IndexError: - # not on plot axes - pass - for coinc in incomplete_injection_coincs(contents): - try: - self.incomplete_coinc_injection_bins[self.coords(contents, coinc.coinc_event_id)] += 1 - except IndexError: - # not on plot axes - pass - - def finish(self): - self.axes.set_title(r"\begin{center}Distribution of Coincident Events (%d Foreground, %d Background Events, %d Injections Found in Coincidence, Logarithmic Density Contours)\end{center}" % (self.n_foreground, self.n_background, self.n_injections)) - xcoords, ycoords = self.background_bins.centres() - - # prepare the data - rate.filter_array(self.foreground_bins.array, rate.gaussian_window(8, 8)) - rate.filter_array(self.background_bins.array, rate.gaussian_window(8, 8)) - rate.filter_array(self.coinc_injection_bins.array, rate.gaussian_window(8, 8)) - rate.filter_array(self.incomplete_coinc_injection_bins.array, rate.gaussian_window(8, 8)) - self.foreground_bins.logregularize() - self.background_bins.logregularize() - self.coinc_injection_bins.logregularize() - self.incomplete_coinc_injection_bins.logregularize() - - # plot background contours - max_density = math.log(self.background_bins.array.max()) - self.axes.contour(xcoords, ycoords, numpy.transpose(numpy.log(self.background_bins.array)), [max_density - n for n in xrange(0, 10, 1)], cmap = matplotlib.cm.Greys) - - # plot foreground (zero-lag) contours - max_density = math.log(self.foreground_bins.array.max()) - self.axes.contour(xcoords, ycoords, numpy.transpose(numpy.log(self.foreground_bins.array)), [max_density - n for n in xrange(0, 10, 1)], cmap = matplotlib.cm.Reds) - #self.axes.plot(self.foreground_x, self.foreground_y, "r+") - - # plot coincident injection contours - max_density = math.log(self.coinc_injection_bins.array.max()) - self.axes.contour(xcoords, ycoords, numpy.transpose(numpy.log(self.coinc_injection_bins.array)), [max_density - n for n in xrange(0, 10, 1)], cmap = matplotlib.cm.Blues) - - # plot incomplete coincident injection contours - max_density = math.log(self.incomplete_coinc_injection_bins.array.max()) - self.axes.contour(xcoords, ycoords, numpy.transpose(numpy.log(self.incomplete_coinc_injection_bins.array)), [max_density - n for n in xrange(0, 10, 1)], cmap = matplotlib.cm.Greens) - - # fix axes limits - self.axes.set_xlim([self.background_bins.bins.min[0], self.background_bins.bins.max[0]]) - self.axes.set_ylim([self.background_bins.bins.min[1], self.background_bins.bins.max[1]]) - - -# -# ============================================================================= -# -# Rate vs. Confidence -# -# ============================================================================= -# - - -class RateVsConfidence(object): - def __init__(self, instrument): - self.fig, self.axes = SnglBurstUtils.make_burst_plot("%s Confidence" % instrument, "Coincident Event Rate (Hz)") - self.instrument = instrument - self.foreground = [] - self.background = [] - self.foreground_segs = segments.segmentlist() - self.background_segs = segments.segmentlist() - self.axes.loglog() - - def add_contents(self, contents): - for time_slide_id, not_zero_lag in contents.connection.cursor().execute("SELECT DISTINCT time_slide_id, EXISTS ( SELECT * FROM time_slide AS a WHERE a.time_slide_id == time_slide.time_slide_id AND a.offset != 0 ) FROM time_slide"): - if not_zero_lag: - bins = self.background - self.background_segs.append(contents.seglists[self.instrument]) - else: - bins = self.foreground - self.foreground_segs.append(contents.seglists[self.instrument]) - bins.extend(contents.connection.cursor().execute("SELECT sngl_burst.confidence FROM coinc_event JOIN coinc_event_map ON (coinc_event_map.coinc_event_id == coinc_event.coinc_event_id) JOIN sngl_burst ON (coinc_event_map.table_name == 'sngl_burst' AND coinc_event_map.event_id == sngl_burst.event_id) WHERE coinc_event.time_slide_id == ? AND coinc_event.coinc_def_id == ? AND sngl_burst.ifo == ?""", (time_slide_id, contents.bb_definer_id, self.instrument))) - - def finish(self): - self.axes.set_title("Cummulative Coincident Event Rate vs. Confidence in %s" % self.instrument) - self.background.sort() - self.foreground.sort() - background_y = numpy.arange(len(self.background), 0.0, -1.0, "Float64") / float(abs(self.background_segs)) - foreground_y = numpy.arange(len(self.foreground), 0.0, -1.0, "Float64") / float(abs(self.foreground_segs)) - self.axes.plot(self.background, background_y, "ko-") - self.axes.plot(self.foreground, foreground_y, "ro-", markeredgecolor = "r") - - -# -# ============================================================================= -# -# Main -# -# ============================================================================= -# - - -def new_plots(plots): - deltat_seg = segments.segment(-0.3, +0.3) - deltat_width = 0.03125 - l = [ - RateContours("H2", "H1"), - ConfidenceContours("H2", "H1", magnitude_a, "Confidence", 1, 10**10), - ConfidenceContours("H2", "L1", magnitude_a, "Confidence", 1, 10**10), - ConfidenceContours("L1", "H1", magnitude_a, "Confidence", 1, 10**10), - ConfidenceContours("H2", "H1", magnitude_b, r"Power / D.o.F. / ($F_{+}^{2} + F_{\times}^{2}$)", 1, 10**10), - ConfidenceContours("H2", "L1", magnitude_b, r"Power / D.o.F. / ($F_{+}^{2} + F_{\times}^{2}$)", 1, 10**10), - ConfidenceContours("L1", "H1", magnitude_b, r"Power / D.o.F. / ($F_{+}^{2} + F_{\times}^{2}$)", 1, 10**10), - ConfidenceContourProjection(numpy.array((-1/math.sqrt(2), +1/math.sqrt(2), 0), "Float64"), numpy.array((-1/math.sqrt(4), -1/math.sqrt(4), +1/math.sqrt(2)), "Float64"), magnitude_b, 10**5), - RateVsConfidence("H1"), - RateVsConfidence("H2"), - RateVsConfidence("L1") - ] - return [l[i] for i in plots] - - -options, filenames = parse_command_line() - - -plots = new_plots(options.plot) - - -for n, filename in enumerate(ligolw_utils.sort_files_by_size(filenames, options.verbose, reverse = True)[options.skip:]): - if options.verbose: - print("%d/%d: %s" % (n + 1, len(filenames) - options.skip, filename), file=sys.stderr) - - database = SnglBurstUtils.CoincDatabase(sqlite3.connect(filename), options.live_time_program) - if options.verbose: - SnglBurstUtils.summarize_coinc_database(database) - - for n, plot in zip(options.plot, plots): - if options.verbose: - print("adding to burca plot %d ..." % n, file=sys.stderr) - plot.add_contents(database) - - database.connection.close() - - -# delete the plots as we go to save memory -n = 0 -format = "%%s%%0%dd.%%s" % (int(math.log10(max(options.plot) or 1)) + 1) -while len(plots): - filename = format % (options.base, options.plot[n], options.format) - if options.verbose: - print("finishing plot %d ..." % options.plot[n], file=sys.stderr) - plots[0].finish() - if options.verbose: - print("writing %s ..." % filename, file=sys.stderr) - plots[0].fig.savefig(filename) - del plots[0] - n += 1 - -if options.verbose: - print("done.", file=sys.stderr) diff --git a/lalapps/src/power/lalapps_power_plot_burca2.py b/lalapps/src/power/lalapps_power_plot_burca2.py deleted file mode 100644 index b4da277c25..0000000000 --- a/lalapps/src/power/lalapps_power_plot_burca2.py +++ /dev/null @@ -1,282 +0,0 @@ -# -# 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 " -__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) diff --git a/lalapps/src/power/lalapps_power_plot_burst.py b/lalapps/src/power/lalapps_power_plot_burst.py deleted file mode 100644 index f183e52f85..0000000000 --- a/lalapps/src/power/lalapps_power_plot_burst.py +++ /dev/null @@ -1,520 +0,0 @@ -# -# 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 " -__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 -# -# ============================================================================= -# - - -class Delays(object): - def __init__(self, ifo, width, max): - self.fig, self.axes = SnglBurstUtils.make_burst_plot("Delay (s)", "Count / Delay") - self.ifo = ifo - self.nevents = 0 - # 21 bins per filter width - interval = segments.segment(0, max + 2) - self.bins = rate.BinnedDensity(rate.NDBins((rate.LinearBins(interval[0], interval[1], int(float(abs(interval)) / width) * 21),))) - self.axes.semilogy() - - def add_contents(self, contents, seglists): - self.nevents += contents[self.ifo].nevents - peaks = list(contents[self.ifo].peak_time) - peaks.sort() - for i in xrange(1, len(peaks)): - dt = float(peaks[i] - peaks[i - 1]) - try: - self.bins.count[dt,] += 1 - except IndexError: - # out of bounds - pass - - def finish(self): - self.axes.set_title("Time Between Triggers\n(%d Triggers)" % self.nevents) - - rate.filter_array(self.bins.array, rate.gaussian_window(21)) - xvals = self.bins.centres()[0] - yvals = self.bins.at_centres() - self.axes.plot(xvals, yvals, "k") - - self.axes.set_xlim((0, xvals[-1])) - self.axes.set_ylim((1, 10.0**(int(math.log10(yvals.max())) + 1))) - - -# -# ============================================================================= -# -# Rate vs. SNR -# -# ============================================================================= -# - - -class RateVsSNR(object): - def __init__(self, ifo): - self.fig, self.axes = SnglBurstUtils.make_burst_plot("SNR", "Trigger Rate (Hz)") - self.ifo = ifo - self.nevents = 0 - self.x = [] - self.seglist = segments.segmentlist() - self.axes.loglog() - - def add_contents(self, contents, seglists): - self.nevents += contents[self.ifo].nevents - self.x.extend(contents[self.ifo].snr) - self.seglist |= seglists[self.ifo] - - def finish(self): - self.axes.set_title("Cummulative Trigger Rate vs. SNR\n(%d Triggers)" % self.nevents) - self.x.sort() - self.y = numpy.arange(len(self.x), 0.0, -1.0) / float(abs(self.seglist)) - self.axes.plot(self.x, self.y, "ko-") - - -# -# ============================================================================= -# -# Rate vs. Confidence -# -# ============================================================================= -# - - -class RateVsConfidence(object): - 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 add_contents(self, contents, seglists): - self.nevents += contents[self.ifo].nevents - self.x.extend(contents[self.ifo].confidence) - self.seglist |= seglists[self.ifo] - - def finish(self): - self.axes.set_title("Cummulative Trigger Rate vs. Confidence\n(%d Triggers)" % self.nevents) - self.x.sort() - self.y = numpy.arange(len(self.x), 0.0, -1.0, "Float64") / float(abs(self.seglist)) - self.axes.plot(self.x, self.y, "ko-") - - -# -# ============================================================================= -# -# Time-Frequency Plane -# -# ============================================================================= -# - - -# moved from viz.py in pylal, but looks like it was ripped off of somebody -# else at some point. maybe rectfill() in matplotlib? -def tfplot(*args, **kwargs): - """ - tfplot(x, y, s=20, c='b', marker='o', cmap=None, norm=None, - vmin=None, vmax=None, alpha=1.0) - - Supported function signatures: - - TFPLOT(x, y) : make a scatter plot of x vs y - - TFPLOT(x, y, s) : make a scatter plot of x vs y with size in area - given by s - - TFPLOT(x, y, s, c) : make a scatter plot of x vs y with size in area - given by s and colors given by c - - TFPLOT(x, y, s, c, **kwargs) : control colormapping and scaling - with keyword args; see below - - Make a scatter plot of x versus y. s is a size in points^2 a scalar - or an array of the same length as x or y. c is a color and can be a - """ - shading = kwargs.get('shading', 'faceted') - cmap = kwargs.get('cmap', cm.get_cmap()) - norm = kwargs.get('norm', normalize()) - alpha = kwargs.get('alpha', 1.0) - vmin = kwargs.get('vmin', None) - vmax = kwargs.get('vmax', None) - a = kwargs.get('axes', gca()) - - try: - X, dX, Y, dY, C = args - except ValueError: - raise TypeError('Illegal arguments to rectfill; see help(rectfill)') - - Nx, = X.shape - verts = [( - (X[i,], Y[i,] ), - (X[i,]+dX[i,], Y[i,] ), - (X[i,]+dX[i,], Y[i,]+dY[i,]), - (X[i,], Y[i,]+dY[i,])) for i in range(Nx-1)] - C = array([C[i,] for i in range(Nx-1)]) - - if shading == 'faceted': - edgecolors = (0, 0, 0, 1), - else: - edgecolors = 'None' - - collection = PolyCollection(verts, edgecolors = edgecolors, antialiaseds = (0,), linewidths = (0.25,)) - collection.set_alpha(alpha) - collection.set_array(C) - if norm is not None: - assert isinstance(norm, normalize) - if cmap is not None: - assert isinstance(cmap, Colormap) - collection.set_cmap(cmap) - collection.set_norm(norm) - if norm is not None: - collection.set_clim(vmin, vmax) - minx = amin(X) - maxx = amax(X) - miny = amin(Y) - maxy = amax(Y) - corners = (minx, miny), (maxx, maxy) - a.update_datalim( corners ) - a.autoscale_view() - # add the collection last - a.add_collection(collection) - return collection - - -class TimeFrequencyPlane(object): - def __init__(self, ifo): - self.fig, self.axes = SnglBurstUtils.make_burst_plot("GPS Time (s)", "Frequency (Hz)") - self.ifo = ifo - self.nevents = 0 - self.seglist = segments.segmentlist() - - def add_contents(self, contents, seglists): - self.nevents += contents[self.ifo].nevents - tfplot(numpy.array(map(float, contents[self.ifo].start_time)), numpy.array(contents[self.ifo].duration), numpy.array(contents[self.ifo].lo_freq), numpy.array(contents[self.ifo].bandwidth), numpy.log(numpy.array(contents[self.ifo].confidence)), axes = self.axes) - self.seglist |= seglists[self.ifo] - - def finish(self): - self.axes.set_title("Time-Frequency Plane\n(%d Triggers)" % self.nevents) - 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) - - -# -# ============================================================================= -# -# Load Data -# -# ============================================================================= -# - - -summary = {} -seglists = segments.segmentlistdict() - - -for n, filename in enumerate(ligolw_utils.sort_files_by_size(filenames, options.verbose, reverse = True)): - 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) - seglists |= ligolw_search_summary.segmentlistdict_fromsearchsummary(xmldoc, options.livetime_program).coalesce() - xmldoc.unlink() - - -# -# ============================================================================= -# -# Plot -# -# ============================================================================= -# - - -def new_plots(ifo, plots): - l = ( - RateVsPeakFreq(ifo, segments.segment(options.frequency_range), 4), - Durations(ifo), - Delays(ifo, 0.25, 20), - RateVsSNR(ifo), - RateVsConfidence(ifo), - ConfidenceVsTime(ifo), - ConfidenceVsFrequencyScatter(ifo), - TimeFrequencyPlane(ifo) - ) - return [l[i] for i in plots] - - -for ifo in summary.keys(): - format = "%%s%s_%%0%dd.%%s" % (ifo, int(math.log10(max(options.plot) or 1)) + 1) - for plotnum, plot in zip(options.plot, new_plots(ifo, options.plot)): - filename = format % (options.base, plotnum, options.format) - if options.verbose: - print("adding to %s plot %d ..." % (ifo, plotnum), file=sys.stderr) - plot.add_contents(summary, seglists) - if options.verbose: - print("finishing %s plot %d ..." % (ifo, plotnum), file=sys.stderr) - plot.finish() - if options.verbose: - print("writing %s ..." % filename, file=sys.stderr) - plot.fig.savefig(filename) diff --git a/lalapps/src/power/lalapps_power_veto.py b/lalapps/src/power/lalapps_power_veto.py deleted file mode 100644 index 314f068db9..0000000000 --- a/lalapps/src/power/lalapps_power_veto.py +++ /dev/null @@ -1,336 +0,0 @@ -# -# Copyright (C) 2007 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 - - -from optparse import OptionParser -import sqlite3 -import sys - - -from ligo.lw import dbtables -from ligo.lw.utils import segments as ligolw_segments -from lalburst import git_version -from lalburst import SnglBurstUtils -from ligo import segments - - -__author__ = "Kipp Cannon " -__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, - usage = "%prog [options] --veto-segments-db=filename files ...", - description = "%prog probably does something..." - ) - parser.add_option("--no-vacuum", action = "store_true", default = False, help = "Don't vacuum the database after vetos. Vacuuming reduces the database size and improves the efficiency of indices, but can take a while so you might want to skip it.") - parser.add_option("-t", "--tmp-space", metavar = "path", help = "Path to a directory suitable for use as a work area while manipulating the database file. The database file will be worked on in this directory, and then moved to the final location when complete. This option is intended to improve performance when running in a networked environment, where there might be a local disk with higher bandwidth than is available to the filesystem on which the final output will reside.") - parser.add_option("--veto-segments-db", metavar = "filename", help = "Load veto segments from this SQLite DB file (required).") - parser.add_option("--veto-segments-name", metavar = "name", default = "sngl_burst_veto", help = "Use the segment lists named this to define the veto times (default = \"sngl_burst_veto\").") - parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.") - options, filenames = parser.parse_args() - - if options.veto_segments_db is None: - raise ValueError("missing required --veto-segments-db option") - - return options, (filenames or [None]) - - -# -# ============================================================================= -# -# Load Segments -# -# ============================================================================= -# - - -def load_segments(filename, name, verbose = False): - if verbose: - print("loading \"%s\" segments ... " % name, end=' ', file=sys.stderr) - connection = sqlite3.connect(filename) - segs = ligolw_segments.segmenttable_get_by_name(dbtables.get_xml(connection), name).coalesce() - connection.close() - if verbose: - print("done.", file=sys.stderr) - for ifo in segs: - print("loaded %d veto segment(s) for %s totalling %g s" % (len(segs[ifo]), ifo, float(abs(segs[ifo]))), file=sys.stderr) - return segs - - -# -# ============================================================================= -# -# Excess Power Veto -# -# ============================================================================= -# - - -def apply_excess_power_veto(contents, veto_segs, verbose = False): - if verbose: - SnglBurstUtils.summarize_coinc_database(contents) - - def sngl_burst_is_vetoed(ifo, start, start_ns, duration, veto_segs = veto_segs): - start = dbtables.lsctables.LIGOTimeGPS(start, start_ns) - return ifo in veto_segs and veto_segs[ifo].intersects_segment(segments.segment(start, start + duration)) - contents.connection.create_function("sngl_burst_is_vetoed", 4, sngl_burst_is_vetoed) - cursor = contents.connection.cursor() - - # - # Delete burst <--> burst coincs containing vetoed bursts - # - - if verbose: - print("applying excess power event veto strategy:", file=sys.stderr) - print("\tremoving vetoed burst <--> burst coincs ...", file=sys.stderr) - cursor.execute(""" -DELETE FROM - coinc_event -WHERE - coinc_def_id == ? - AND coinc_event_id IN ( - SELECT - coinc_event_id - FROM - coinc_event_map - JOIN sngl_burst ON ( - coinc_event_map.table_name == 'sngl_burst' - AND coinc_event_map.event_id == sngl_burst.event_id - ) - WHERE - sngl_burst_is_vetoed(sngl_burst.ifo, sngl_burst.start_time, sngl_burst.start_time_ns, sngl_burst.duration) - ) - """, (contents.bb_definer_id,)) - - # - # Delete sim <--> coinc coincs pointing to deleted coincs - # - - if contents.sc_definer_id is not None: - if verbose: - print("\tremoving vetoed sim <--> coinc coincs ...", file=sys.stderr) - cursor.execute(""" -DELETE FROM - coinc_event -WHERE - coinc_def_id == ? - AND coinc_event_id IN ( - SELECT - coinc_event_id - FROM - coinc_event_map - WHERE - table_name == 'coinc_event' - AND event_id NOT IN ( - SELECT - coinc_event_id - FROM - coinc_event - ) - ) - """, (contents.sc_definer_id,)) - - # - # Now that we no longer need to form links through the - # coinc_def_map table, delete vetoed burst rows from it - # - - if verbose: - print("\tremoving vetoed bursts from coinc_def_map table ...", file=sys.stderr) - cursor.execute(""" -DELETE FROM - coinc_event_map -WHERE - coinc_event_map.table_name == 'sngl_burst' - AND coinc_event_map.event_id IN ( - SELECT - event_id - FROM - sngl_burst - WHERE - sngl_burst_is_vetoed(ifo, start_time, start_time_ns, duration) - ) - """) - - # - # The coinc_def_map table no longer contains any rows pointing to - # vetoed bursts, update the event counts in sim <--> burst coincs - # and delete any for which this is now 0 - # - - if contents.sb_definer_id is not None: - if verbose: - print("\tupdating sim <--> burst event counts ...", file=sys.stderr) - cursor.execute(""" -UPDATE - coinc_event -SET - nevents = (SELECT COUNT(*) FROM coinc_event_map WHERE coinc_event_map.coinc_event_id == coinc_event.coinc_event_id AND coinc_event_map.table_name == 'sngl_burst') -WHERE - coinc_event.coinc_def_id == ? - """, (contents.sb_definer_id,)) - if verbose: - print("\tremoving empty sim <--> burst coincs ...", file=sys.stderr) - cursor.execute(""" -DELETE FROM - coinc_event -WHERE - coinc_def_id == ? - AND nevents == 0 - """, (contents.sb_definer_id,)) - - # - # We have removed all the rows from the coinc_event table that will - # be removed, so now remove multi_burst and coinc_event_map rows - # that point to non-existent coinc_events - # - - if verbose: - print("\ttrimming coinc_event_map table ...", file=sys.stderr) - cursor.execute(""" -DELETE FROM - coinc_event_map -WHERE - coinc_event_map.coinc_event_id NOT IN ( - SELECT - coinc_event_id - FROM - coinc_event - ) - """) - - if verbose: - print("\ttrimming multi_burst table ...", file=sys.stderr) - cursor.execute(""" -DELETE FROM - multi_burst -WHERE - multi_burst.coinc_event_id NOT IN ( - SELECT - coinc_event_id - FROM - coinc_event - ) - """) - - # - # Tables are now consistent, and contain no coincs involving a - # vetoed burst - # - - if verbose: - print("\tdone (excess power event vetos).", file=sys.stderr) - - -# -# ============================================================================= -# -# Main -# -# ============================================================================= -# - - -# -# Command line. -# - - -options, filenames = parse_command_line() - - -# -# Load veto segments. -# - - -veto_segs = load_segments(options.veto_segments_db, options.veto_segments_name, verbose = options.verbose) - - -# -# Iterate over files -# - - -for n, filename in enumerate(filenames): - # - # Open the database file. - # - - if options.verbose: - print("%d/%d: %s" % (n + 1, len(filenames), filename), end=' ', file=sys.stderr) - - working_filename = dbtables.get_connection_filename(filename, tmp_path = options.tmp_space, verbose = options.verbose) - connection = sqlite3.connect(working_filename) - connection.execute("PRAGMA synchronous = OFF;") - connection.execute("PRAGMA temp_store_directory = '%s';" % dbtables.tempfile.gettempdir()) - - # - # Apply vetoes - # - - apply_excess_power_veto(SnglBurstUtils.CoincDatabase(connection, "lalapps_power"), veto_segs, verbose = options.verbose) - - # - # Clean up - # - - if options.verbose: - print("committing ...", file=sys.stderr) - connection.commit() - if not options.no_vacuum: - if options.verbose: - print("vacuuming ...", file=sys.stderr) - connection.cursor().execute("VACUUM;") - connection.close() - del connection - dbtables.put_connection_filename(filename, working_filename, verbose = options.verbose) - if options.verbose: - print("done (%s)." % filename, file=sys.stderr) - - -# -# Done -# - - -if options.verbose: - print("done.", file=sys.stderr) -- GitLab From a41df6c87a4869ccb670b2734798ea663da874e5 Mon Sep 17 00:00:00 2001 From: Kipp Cannon Date: Tue, 3 Mar 2020 16:44:14 +0900 Subject: [PATCH 4/5] lalburst: remove six.h --- lalburst/python/lalburst/Makefile.am | 2 -- lalburst/python/lalburst/cs_gamma.c | 40 +++++++++++++++++----------- lalburst/python/lalburst/six.h | 1 - 3 files changed, 25 insertions(+), 18 deletions(-) delete mode 120000 lalburst/python/lalburst/six.h diff --git a/lalburst/python/lalburst/Makefile.am b/lalburst/python/lalburst/Makefile.am index 18c9f9732b..712468e403 100644 --- a/lalburst/python/lalburst/Makefile.am +++ b/lalburst/python/lalburst/Makefile.am @@ -6,8 +6,6 @@ include $(top_srcdir)/gnuscripts/lalsuite_vcs_info.am vcs_info_sources = git_version.py -noinst_HEADERS = six.h - if HAVE_PYTHON pymoduledir = $(pkgpythondir) diff --git a/lalburst/python/lalburst/cs_gamma.c b/lalburst/python/lalburst/cs_gamma.c index 84695c7b7b..27ea0cf348 100644 --- a/lalburst/python/lalburst/cs_gamma.c +++ b/lalburst/python/lalburst/cs_gamma.c @@ -46,7 +46,6 @@ #include #include #include -#include "six.h" #define CUSPS_PER_LOOP 1.0 /* c */ #define LOOP_RAD_POWER 50.0 /* Gamma */ @@ -395,20 +394,31 @@ static PyMethodDef cs_gammaMethods[] = { {NULL, NULL, 0, NULL} }; -static PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, - "cs_gamma", NULL, -1, cs_gammaMethods, - NULL, NULL, NULL, NULL -}; - -//They Python module initialization function. -PyMODINIT_FUNC PyInit_cs_gamma(void); /* silence no-previous-prototype warning */ -PyMODINIT_FUNC -PyInit_cs_gamma(void) +//Then Python module initialization function. +#if PY_MAJOR_VERSION < 3 +PyMODINIT_FUNC initcs_gamma(void); /* silence -Wmissing-prototypes */ +PyMODINIT_FUNC initcs_gamma(void) +#else +PyMODINIT_FUNC PyInit_cs_gamma(void); /* silence -Wmissing-prototypes */ +PyMODINIT_FUNC PyInit_cs_gamma(void) +#endif { - import_array(); - return PyModule_Create(&moduledef); -} +#if PY_MAJOR_VERSION < 3 + PyObject *module = Py_InitModule3("cs_gamma", cs_gammaMethods, NULL); +#else + static PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "cs_gamma", NULL, -1, cs_gammaMethods, + NULL, NULL, NULL, NULL + }; + PyObject *module = PyModule_Create(&moduledef); +#endif + import_array(); -SIX_COMPAT_MODULE(cs_gamma) +#if PY_MAJOR_VERSION < 3 + return; +#else + return module; +#endif +} diff --git a/lalburst/python/lalburst/six.h b/lalburst/python/lalburst/six.h deleted file mode 120000 index c6559f6900..0000000000 --- a/lalburst/python/lalburst/six.h +++ /dev/null @@ -1 +0,0 @@ -../../../gnuscripts/six.h \ No newline at end of file -- GitLab From 0ba57da0eda799b4fe78dc892de53f436f81a509 Mon Sep 17 00:00:00 2001 From: Kipp Cannon Date: Wed, 13 Nov 2019 18:25:04 +0900 Subject: [PATCH 5/5] lalapps_bucut, lalapps_binjfind: API updates - remove references to obsolete InterningRowBuilder type. --- lalapps/src/power/lalapps_binjfind.py | 18 ------------------ lalapps/src/power/lalapps_bucut.py | 11 +---------- 2 files changed, 1 insertion(+), 28 deletions(-) diff --git a/lalapps/src/power/lalapps_binjfind.py b/lalapps/src/power/lalapps_binjfind.py index 29d0ad040a..2cc1373cc2 100644 --- a/lalapps/src/power/lalapps_binjfind.py +++ b/lalapps/src/power/lalapps_binjfind.py @@ -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? diff --git a/lalapps/src/power/lalapps_bucut.py b/lalapps/src/power/lalapps_bucut.py index 0579ccf2ce..94eca569f6 100644 --- a/lalapps/src/power/lalapps_bucut.py +++ b/lalapps/src/power/lalapps_bucut.py @@ -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() -- GitLab