Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALApps 10.1.0.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
lalapps_string_meas_likelihood.py
Go to the documentation of this file.
1##python
2# Copyright (C) 2010--2018 Kipp Cannon
3#
4# This program is free software; you can redistribute it and/or modify it
5# under the terms of the GNU General Public License as published by the
6# Free Software Foundation; either version 2 of the License, or (at your
7# option) any later version.
8#
9# This program is distributed in the hope that it will be useful, but
10# WITHOUT ANY WARRANTY; without even the implied warranty of
11# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
12# Public License for more details.
13#
14# You should have received a copy of the GNU General Public License along
15# with this program; if not, write to the Free Software Foundation, Inc.,
16# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17
18
19#
20# =============================================================================
21#
22# Preamble
23#
24# =============================================================================
25#
26
27
28from __future__ import print_function
29
30
31import math
32from optparse import OptionParser
33import sqlite3
34import string
35import sys
36
37
38from lal.utils import CacheEntry
39
40
41from igwn_ligolw import dbtables
42from igwn_ligolw import ligolw
43from igwn_ligolw import utils as ligolw_utils
44from igwn_ligolw.utils import process as ligolw_process
45from igwn_ligolw.utils import search_summary as ligolw_search_summary
46from lalburst import SnglBurstUtils
47from lalburst import git_version
48from lalburst import stringutils
49import igwn_segments as segments
50
51
52__author__ = "Kipp Cannon <kipp.cannon@ligo.org>"
53__version__ = "git id %s" % git_version.id
54__date__ = git_version.date
55
56
57# characters allowed to appear in the description string
58T010150_letters = set(string.ascii_lowercase + string.ascii_uppercase + string.digits + "_+#")
59
60
61#
62# =============================================================================
63#
64# Command Line
65#
66# =============================================================================
67#
68
69
71 parser = OptionParser(
72 version = "Name: %%prog\n%s" % git_version.verbose_msg,
73 usage = "%prog [options] [filename ...]",
74 description = "%prog analyzes a collection of sqlite3 database files containing lalburst_coinc outputs of string-cusp coincidence events, and measures probability distributions for a variety of parameters computed from those coincidences. The distributions are written to a likelihood data file in XML format, which can later be used by to assign likelihoods to the coincidences. The files to be processed can be named on the command line and/or provided by a LAL cache file."
75 )
76 parser.add_option("-o", "--output", metavar = "filename", default = None, help = "Set the name of the likelihood data file to write (default = stdout).")
77 parser.add_option("-c", "--input-cache", metavar = "filename", help = "Also process the files named in this LAL cache. See lalapps_path2cache for information on how to produce a LAL cache file.")
78 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.")
79 parser.add_option("--T010150", metavar = "description", default = None, help = "Write the output to a file whose name is compatible with the file name format described in LIGO-T010150-00-E, \"Naming Convention for Frame Files which are to be Processed by LDAS\". The description string will be used to form the second field in the file name.")
80 parser.add_option("--injection-reweight", metavar = "off|astrophysical", default = "off", help = "Set the weight function to be applied to the injections (default = \"off\"). When \"off\", the injections are all given equal weight and so the injection population is whatever was injected. When set to \"astrophysical\", the injections are reweighted to simulate an amplitude^{-4} distribution.")
81 parser.add_option("--injection-reweight-cutoff", metavar = "amplitude", default = 1e-20, type = "float", help = "When using the astrophysical injection reweighting, do not allow the weight assigned to arbitrarily low-amplitude injections to grow without bound, instead clip the weight assigned to injections to the weight given to injections with this amplitude (default = 1e-20, 0 = disabled). This option is ignored when astrophysical reweighting is not being performed.")
82 parser.add_option("--vetoes-name", metavar = "name", help = "Set the name of the segment lists to use as vetoes (default = do not apply vetoes).")
83 parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")
84 options, filenames = parser.parse_args()
85
86 paramdict = options.__dict__.copy()
87
88 if options.T010150 is not None:
89 if options.output is not None:
90 raise ValueError("cannot set both --T010150 and --output")
91 if options.T010150 == "":
92 options.T010150 = "STRING_LIKELIHOOD"
93 elif set(options.T010150) - T010150_letters:
94 raise ValueError("invalid characters in description \"%s\"" % options.T010150)
95
96 if options.input_cache:
97 filenames += [CacheEntry(line).path for line in file(options.input_cache)]
98
99 if not filenames:
100 raise ValueError("no input files!")
101
102 if options.injection_reweight not in ("off", "astrophysical"):
103 raise ValueError("--injection-reweight \"%s\" not recognized" % options.injections_reweight)
104
105 return options, filenames, paramdict
106
107
108#
109# =============================================================================
110#
111# Injection ReWeight Functions
112#
113# =============================================================================
114#
115
116
117def get_injection_weight_func(contents, reweight_type, amplitude_cutoff):
118 if reweight_type == "off":
119 def weight_func(sim, amplitude_cutoff = amplitude_cutoff):
120 # amplitude cut-off is not used
121 return 1.0
122 elif reweight_type == "astrophysical":
123 population, = ligolw_process.get_process_params(contents.xmldoc, "lalapps_binj", "--population")
124 if population != "string_cusp":
125 raise ValueError("lalapps_binj was not run with --population=\"string_cusp\"")
126 def weight_func(sim, amplitude_cutoff = amplitude_cutoff):
127 # the "string_cusp" injection population is uniform
128 # in log A, meaning the number of injections
129 # between log A and log A + d log A is independent
130 # of A. that corresponds to a distribution density
131 # in A of P(A) dA \propto A^{-1} dA. we want P(A)
132 # dA \propto A^{-4} dA, so each physical injection
133 # needs to be treated as if it was A^{-3} virtual
134 # injections, then the density of virtual
135 # injections will be P(A) dA \propto A^{-4} dA.
136 # the factor of 10^{21} simply renders the
137 # amplitudes closer to 1; since double-precision
138 # arithmetic is used this should have no affect on
139 # the results, but it might help make the numbers a
140 # little easier for humans to look at should anyone
141 # have occasion to do so.
142 return (max(sim.amplitude, amplitude_cutoff) * 1e21)**-3
143 else:
144 raise ValueError(reweight_type)
145 return weight_func
146
147
148#
149# =============================================================================
150#
151# Main
152#
153# =============================================================================
154#
155
156
157#
158# Command line.
159#
160
161
162options, filenames, paramdict = parse_command_line()
163
164
165#
166# Clear the statistics book-keeping object.
167#
168
169
170# FIXME: don't hard-code instruments
171distributions = stringutils.StringCoincParamsDistributions(["H1", "L1", "V1"])
172segs = segments.segmentlistdict()
173
174
175#
176# Start output document
177#
178
179
180xmldoc = ligolw.Document()
181xmldoc.appendChild(ligolw.LIGO_LW())
182process = ligolw_process.register_to_xmldoc(xmldoc, program = u"lalapps_string_meas_likelihood", paramdict = paramdict, version = __version__, cvs_repository = "lscsoft", cvs_entry_time = __date__, comment = u"")
183
184
185#
186# Iterate over files
187#
188
189
190for n, filename in enumerate(filenames):
191 #
192 # Open the database file.
193 #
194
195 if options.verbose:
196 print("%d/%d: %s" % (n + 1, len(filenames), filename), file=sys.stderr)
197
198 working_filename = dbtables.get_connection_filename(filename, tmp_path = options.tmp_space, verbose = options.verbose)
199 connection = sqlite3.connect(str(working_filename))
200 if options.tmp_space is not None:
201 dbtables.set_temp_store_directory(connection, options.tmp_space, verbose = options.verbose)
202
203 #
204 # Summarize the database.
205 #
206
207 contents = SnglBurstUtils.CoincDatabase(connection, live_time_program = "StringSearch", search = "StringCusp", veto_segments_name = options.vetoes_name)
208 if options.verbose:
209 SnglBurstUtils.summarize_coinc_database(contents)
210 if not contents.seglists and options.verbose:
211 print("\twarning: no segments found", file=sys.stderr)
212 segs |= contents.seglists
213
214 #
215 # Record statistics. Assume all files with sim_burst tables are
216 # the outputs of injection runs, and others aren't.
217 #
218
219 if contents.sim_burst_table is None:
220 # iterate over burst<-->burst coincs
221 for is_background, events, offsetvector in contents.get_noninjections():
222 params = distributions.coinc_params([event for event in events if event.ifo not in contents.vetoseglists or event.peak not in contents.vetoseglists[event.ifo]], offsetvector)
223 if is_background:
224 distributions.denominator.increment(**params)
225 else:
226 distributions.candidates.increment(**params)
227 else:
228 weight_func = get_injection_weight_func(contents, options.injection_reweight, options.injection_reweight_cutoff)
229 # iterate over burst<-->burst coincs matching injections
230 # "exactly"
231 for sim, events, offsetvector in contents.get_injections():
232 distributions.numerator.increment(
233 weight = weight_func(sim),
234 **distributions.coinc_params([event for event in events if event.ifo not in contents.vetoseglists or event.peak not in contents.vetoseglists[event.ifo]], offsetvector)
235 )
236
237 #
238 # Clean up.
239 #
240
241 contents.xmldoc.unlink()
242 connection.close()
243 dbtables.discard_connection_filename(filename, working_filename, verbose = options.verbose)
244
245
246#
247# Output.
248#
249
250
251ligolw_search_summary.append_search_summary(xmldoc, process, ifos = segs.keys(), inseg = segs.extent_all(), outseg = segs.extent_all())
252xmldoc.childNodes[-1].appendChild(distributions.to_xml(u"string_cusp_likelihood"))
253process.set_end_time_now()
254
255
256def T010150_basename(instruments, description, seg):
257 start = int(math.floor(seg[0]))
258 duration = int(math.ceil(seg[1] - start))
259 return "%s-%s-%d-%d" % ("+".join(sorted(instruments)), description, start, duration)
260if options.T010150:
261 filename = "%s.xml.gz" % T010150_basename(segs.keys(), options.T010150, segs.extent_all())
262else:
263 filename = options.output
264ligolw_utils.write_filename(xmldoc, filename, verbose = options.verbose)
def get_injection_weight_func(contents, reweight_type, amplitude_cutoff)
def T010150_basename(instruments, description, seg)