31from igwn_ligolw
import lsctables
32from igwn_ligolw.utils
import process
as ligolw_process
33from igwn_ligolw.utils
import search_summary
as ligolw_search_summary
34import igwn_segments
as segments
35from .
import snglcluster
38__author__ =
"Kipp Cannon <kipp.cannon@ligo.org>"
39from .git_version
import date
as __date__
40from .git_version
import version
as __version__
52process_program_name =
"lalburst_cluster"
56 return ligolw_process.register_to_xmldoc(
58 program = process_program_name,
60 "cluster_algorithm": cluster_algorithm
62 version = __version__,
63 cvs_repository =
"lscsoft",
64 cvs_entry_time = __date__,
91 for colname
in (
"peak_frequency",
"ms_start_time",
"ms_start_time_ns",
"ms_duration",
"ms_flow",
"ms_bandwidth",
"ms_hrss",
"ms_snr",
"ms_confidence"):
93 sngl_burst_table.getColumnByName(colname)
95 sngl_burst_table.appendColumn(colname)
102 for row
in sngl_burst_table:
103 row.peak_frequency = row.central_freq
104 row.ms_period = row.period
105 row.ms_band = row.band
106 row.ms_hrss = row.amplitude
108 row.ms_confidence = row.confidence
127 For speed, convert peak times to floats relative to epoch.
129 if not len(sngl_burst_table):
131 offset = sngl_burst_table[0].peak
132 for row
in sngl_burst_table:
133 row.peak_time = float(row.peak - offset)
140 Restore peak times to absolute LIGOTimeGPS values.
142 for row
in sngl_burst_table:
143 row.peak = offset + row.peak_time
148 Sort key for grouping excess power triggers near triggers with
149 which they might cluster.
151 return (a.ifo, a.channel, a.search, a.start)
156 Returns True if a's and b's (ifo, channel, seach) are different or
157 if the periods they span are disjoint. Knowing excess power
159 then if for a pair of events this function returns
False, we know
160 the result will also be
False for all other events farther apart
in
161 the list. This
is used to terminate the scan
for events to
164 return (a.ifo, a.channel, a.search) != (b.ifo, b.channel, b.search) or a.period.disjoint(b.period)
167def ExcessPowerTestFunc(a, b):
169 Return
False if a
and b cluster. To cluster, two events must be
170 from the same channel of the same instrument,
and their
171 time-frequency tiles must be non-disjoint.
173 return (a.ifo, a.channel, a.search) != (b.ifo, b.channel, b.search) or a.period.disjoint(b.period) or a.band.disjoint(b.band)
176def ExcessPowerClusterFunc(a, b):
178 Modify a
in place to be a cluster constructed
from a
and b. The
179 cluster
's time-frequency tile is the smallest tile that contains
180 the original two tiles, and the
"most signficiant" contributor
for
181 the cluster
is the tile whose boundaries are the SNR^{2} weighted
182 average boundaries of the two contributing tiles. The
"most
183 signficiant" contributor's h_{rss}, SNR, and confidence, are copied
184 verbatim from whichever of the two contributing tiles has the
185 highest confidence. The modified event a
is returned.
188 # In the special case of the two events being the exact same
189 # time-frequency tile, simply preserve the one with the highest
190 # confidence and discard the other.
193 if a.period == b.period and a.band == b.band:
194 if b.ms_confidence > a.ms_confidence:
202 if b.ms_confidence > a.ms_confidence:
203 a.ms_hrss = b.ms_hrss
205 a.ms_confidence = b.ms_confidence
206 a.ms_period = snglcluster.weighted_average_seg(a.ms_period, a.snr**2.0, b.ms_period, b.snr**2.0)
207 a.ms_band = snglcluster.weighted_average_seg(a.ms_band, a.snr**2.0, b.ms_band, b.snr**2.0)
215 a.peak_time = (a.snr**2.0 * a.peak_time + b.snr**2.0 * b.peak_time) / (a.snr**2.0 + b.snr**2.0)
216 a.peak_frequency = (a.snr**2.0 * a.peak_frequency + b.snr**2.0 * b.peak_frequency) / (a.snr**2.0 + b.snr**2.0)
227 a.amplitude += b.amplitude
228 a.snr = math.sqrt(a.snr**2.0 + b.snr**2.0)
234 a.confidence = a.ms_confidence
241 a.band = snglcluster.smallest_enclosing_seg(a.band, b.band)
248 a.period = snglcluster.smallest_enclosing_seg(a.period, b.period)
259 Modify a in place to be a cluster constructed from a and b. The
260 cluster's time-frequency tile is the smallest tile that contains
261 the original two tiles, and the
"most signficiant" contributor
for
262 the cluster
is the tile whose boundaries are the SNR^{2} weighted
263 average boundaries of the two contributing tiles. The
"most
264 signficiant" contributor's h_{rss}, SNR, and confidence, are copied
265 verbatim from whichever of the two contributing tiles has the
266 highest confidence. The modified event a
is returned.
269 # In the special case of the two events being the exact same
270 # time-frequency tile, simply preserve the one with the highest
271 # confidence and discard the other.
274 if a.period == b.period and a.band == b.band:
283 if b.ms_snr > a.ms_snr:
285 a.ms_period = snglcluster.weighted_average_seg(a.ms_period, a.snr**2.0, b.ms_period, b.snr**2.0)
286 a.ms_band = snglcluster.weighted_average_seg(a.ms_band, a.snr**2.0, b.ms_band, b.snr**2.0)
294 a.peak_time = (a.snr**2.0 * a.peak_time + b.snr**2.0 * b.peak_time) / (a.snr**2.0 + b.snr**2.0)
295 a.peak_frequency = (a.snr**2.0 * a.peak_frequency + b.snr**2.0 * b.peak_frequency) / (a.snr**2.0 + b.snr**2.0)
306 a.amplitude += b.amplitude
307 a.snr = math.sqrt(a.snr**2.0 + b.snr**2.0)
314 a.band = snglcluster.smallest_enclosing_seg(a.band, b.band)
321 a.period = snglcluster.smallest_enclosing_seg(a.period, b.period)
351 Run the clustering algorithm on the list of burst candidates. The
352 return value
is the tuple (xmldoc, changed), where xmldoc
is the
353 input document,
and changed
is a boolean that
is True if the
354 contents of the sngl_burst table were altered,
and False if the
355 triggers were
not modified by the clustering process.
357 If the document does
not contain a sngl_burst table, then the
358 document
is not modified (including no modifications to the process
363 # Extract live time segment and sngl_burst table
367 sngl_burst_table = lsctables.SnglBurstTable.get_table(xmldoc)
371 print(
"document does not contain a sngl_burst table, skipping ...", file=sys.stderr)
373 seglists = ligolw_search_summary.segmentlistdict_fromsearchsummary_out(xmldoc, program = program).coalesce()
380 print(
"pre-processing ...", file=sys.stderr)
381 preprocess_output =
prefunc(sngl_burst_table)
387 table_changed = snglcluster.cluster_events(sngl_burst_table, testfunc, clusterfunc, sortkeyfunc = sortkeyfunc, bailoutfunc = bailoutfunc, verbose = verbose)
394 print(
"post-processing ...", file=sys.stderr)
395 postfunc(sngl_burst_table, preprocess_output)
402 process.instruments = seglists.keys()
403 ligolw_search_summary.append_search_summary(xmldoc, process, inseg = seglists.extent_all(), outseg = seglists.extent_all(), nevents = len(sngl_burst_table))
409 return xmldoc, table_changed
def ExcessPowerSortKeyFunc(a)
Sort key for grouping excess power triggers near triggers with which they might cluster.
def ExcessPowerBailoutFunc(a, b)
Returns True if a's and b's (ifo, channel, seach) are different or if the periods they span are disjo...
def ExcessPowerPostFunc(sngl_burst_table, offset)
Restore peak times to absolute LIGOTimeGPS values.
def OmegaClusterFunc(a, b)
Modify a in place to be a cluster constructed from a and b.
def append_process(xmldoc, cluster_algorithm, comment)
def add_ms_columns(xmldoc)
def add_ms_columns_to_table(sngl_burst_table)
def bucluster(xmldoc, program, process, prefunc, postfunc, testfunc, clusterfunc, sortkeyfunc=None, bailoutfunc=None, verbose=False)
Run the clustering algorithm on the list of burst candidates.
def ExcessPowerPreFunc(sngl_burst_table)
For speed, convert peak times to floats relative to epoch.