# Copyright (C) 2012 Alex Nitz
# 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 3 of the License, or (at your
# self.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
#
# =============================================================================
#
"""This modules defines functions for clustering and thresholding timeseries to
produces event triggers
"""
import os.path
import copy
import itertools
import logging
import pickle
import numpy
import h5py
from pycbc.types import Array
from pycbc.scheme import schemed
from pycbc.detector import Detector
from . import coinc, ranking
from .eventmgr_cython import findchirp_cluster_over_window_cython
logger = logging.getLogger('pycbc.events.eventmgr')
[docs]
@schemed("pycbc.events.threshold_")
def threshold(series, value):
"""Return list of values and indices values over threshold in series.
"""
err_msg = "This function is a stub that should be overridden using the "
err_msg += "scheme. You shouldn't be seeing this error!"
raise ValueError(err_msg)
[docs]
@schemed("pycbc.events.threshold_")
def threshold_only(series, value):
"""Return list of values and indices whose values in series are
larger (in absolute value) than value
"""
err_msg = "This function is a stub that should be overridden using the "
err_msg += "scheme. You shouldn't be seeing this error!"
raise ValueError(err_msg)
# FIXME: This should be under schemed, but I don't understand that yet!
[docs]
def threshold_real_numpy(series, value):
arr = series.data
locs = numpy.where(arr > value)[0]
vals = arr[locs]
return locs, vals
[docs]
@schemed("pycbc.events.threshold_")
def threshold_and_cluster(series, threshold, window):
"""Return list of values and indices values over threshold in series.
"""
err_msg = "This function is a stub that should be overridden using the "
err_msg += "scheme. You shouldn't be seeing this error!"
raise ValueError(err_msg)
@schemed("pycbc.events.threshold_")
def _threshold_cluster_factory(series):
err_msg = "This class is a stub that should be overridden using the "
err_msg += "scheme. You shouldn't be seeing this error!"
raise ValueError(err_msg)
[docs]
class ThresholdCluster(object):
"""Create a threshold and cluster engine
Parameters
-----------
series : complex64
Input pycbc.types.Array (or subclass); it will be searched for
points above threshold that are then clustered
"""
def __new__(cls, *args, **kwargs):
real_cls = _threshold_cluster_factory(*args, **kwargs)
return real_cls(*args, **kwargs) # pylint:disable=not-callable
# The class below should serve as the parent for all schemed classes.
# The intention is that this class serves simply as the location for
# all documentation of the class and its methods, though that is not
# yet implemented. Perhaps something along the lines of:
#
# http://stackoverflow.com/questions/2025562/inherit-docstrings-in-python-class-inheritance
#
# will work? Is there a better way?
class _BaseThresholdCluster(object):
def threshold_and_cluster(self, threshold, window):
"""
Threshold and cluster the memory specified at instantiation with the
threshold and window size specified at creation.
Parameters
-----------
threshold : float32
The minimum absolute value of the series given at object initialization
to return when thresholding and clustering.
window : uint32
The size (in number of samples) of the window over which to cluster
Returns:
--------
event_vals : complex64
Numpy array, complex values of the clustered events
event_locs : uint32
Numpy array, indices into series of location of events
"""
pass
[docs]
def findchirp_cluster_over_window(times, values, window_length):
""" Reduce the events by clustering over a window using
the FindChirp clustering algorithm
Parameters
-----------
indices: Array
The list of indices of the SNR values
snr: Array
The list of SNR value
window_size: int
The size of the window in integer samples. Must be positive.
Returns
-------
indices: Array
The reduced list of indices of the SNR values
"""
assert window_length > 0, 'Clustering window length is not positive'
indices = numpy.zeros(len(times), dtype=numpy.int32)
tlen = len(times)
absvalues = numpy.asarray(abs(values))
times = numpy.asarray(times, dtype=numpy.int32)
k = findchirp_cluster_over_window_cython(times, absvalues, window_length,
indices, tlen)
return indices[0:k+1]
[docs]
def cluster_reduce(idx, snr, window_size):
""" Reduce the events by clustering over a window
Parameters
-----------
indices: Array
The list of indices of the SNR values
snr: Array
The list of SNR value
window_size: int
The size of the window in integer samples.
Returns
-------
indices: Array
The list of indices of the SNR values
snr: Array
The list of SNR values
"""
ind = findchirp_cluster_over_window(idx, snr, window_size)
return idx.take(ind), snr.take(ind)
class H5FileSyntSugar(object):
"""Convenience class that adds some syntactic sugar to h5py.File.
"""
def __init__(self, name, prefix=''):
self.f = h5py.File(name, 'w')
self.prefix = prefix
def __setitem__(self, name, data):
self.f.create_dataset(
self.prefix + '/' + name,
data=data,
compression='gzip',
compression_opts=9,
shuffle=True
)
[docs]
class EventManager(object):
def __init__(self, opt, column, column_types, **kwds):
self.opt = opt
self.global_params = kwds
self.event_dtype = [('template_id', int)]
for col, coltype in zip(column, column_types):
self.event_dtype.append((col, coltype))
self.events = numpy.array([], dtype=self.event_dtype)
self.accumulate = [self.events]
self.template_params = []
self.template_index = -1
self.template_events = numpy.array([], dtype=self.event_dtype)
self.write_performance = False
[docs]
def save_state(self, tnum_finished, filename):
"""Save the current state of the background buffers"""
from pycbc.io.hdf import dump_state
self.tnum_finished = tnum_finished
logger.info('Writing checkpoint file at template %s', tnum_finished)
fp = h5py.File(filename, 'w')
dump_state(self, fp, protocol=pickle.HIGHEST_PROTOCOL)
fp.close()
[docs]
@staticmethod
def restore_state(filename):
"""Restore state of the background buffers from a file"""
from pycbc.io.hdf import load_state
fp = h5py.File(filename, 'r')
try:
mgr = load_state(fp)
except Exception as e:
fp.close()
raise e
fp.close()
next_template = mgr.tnum_finished + 1
logger.info('Restoring with checkpoint at template %s', next_template)
return mgr.tnum_finished + 1, mgr
[docs]
@classmethod
def from_multi_ifo_interface(cls, opt, ifo, column, column_types, **kwds):
"""
To use this for a single ifo from the multi ifo interface requires
some small fixing of the opt structure. This does that. As we edit the
opt structure the process_params table will not be correct.
"""
opt = copy.deepcopy(opt)
opt_dict = vars(opt)
for arg, value in opt_dict.items():
if isinstance(value, dict):
setattr(opt, arg, getattr(opt, arg)[ifo])
return cls(opt, column, column_types, **kwds)
[docs]
def chisq_threshold(self, value, num_bins, delta=0):
remove = []
for i, event in enumerate(self.events):
xi = event['chisq'] / (event['chisq_dof'] +
delta * event['snr'].conj() * event['snr'])
if xi > value:
remove.append(i)
self.events = numpy.delete(self.events, remove)
[docs]
def newsnr_threshold(self, threshold):
""" Remove events with newsnr smaller than given threshold
"""
if not self.opt.chisq_bins:
raise RuntimeError('Chi-square test must be enabled in order to '
'use newsnr threshold')
nsnrs = ranking.newsnr(abs(self.events['snr']),
self.events['chisq'] / self.events['chisq_dof'])
remove_idxs = numpy.where(nsnrs < threshold)[0]
self.events = numpy.delete(self.events, remove_idxs)
[docs]
def keep_near_injection(self, window, injections):
from pycbc.events.veto import indices_within_times
if len(self.events) == 0:
return
inj_time = numpy.array(injections.end_times())
gpstime = self.events['time_index'].astype(numpy.float64)
gpstime = gpstime / self.opt.sample_rate + self.opt.gps_start_time
i = indices_within_times(gpstime, inj_time - window, inj_time + window)
self.events = self.events[i]
[docs]
def keep_loudest_in_interval(self, window, num_keep, statname="newsnr",
log_chirp_width=None):
if len(self.events) == 0:
return
e_copy = self.events.copy()
# Here self.events['snr'] is the complex SNR
e_copy['snr'] = abs(e_copy['snr'])
# Messy step because pycbc inspiral's internal 'chisq_dof' is 2p-2
# but stat.py / ranking.py functions use 'chisq_dof' = p
e_copy['chisq_dof'] = e_copy['chisq_dof'] / 2 + 1
statv = ranking.get_sngls_ranking_from_trigs(e_copy, statname)
# Convert trigger time to integer bin number
# NB time_index and window are in units of samples
wtime = (e_copy['time_index'] / window).astype(numpy.int32)
bins = numpy.unique(wtime)
if log_chirp_width:
from pycbc.conversions import mchirp_from_mass1_mass2
m1 = numpy.array([p['tmplt'].mass1 for p in self.template_params])
m2 = numpy.array([p['tmplt'].mass2 for p in self.template_params])
mc = mchirp_from_mass1_mass2(m1, m2)[e_copy['template_id']]
# convert chirp mass to integer bin number
imc = (numpy.log(mc) / log_chirp_width).astype(numpy.int32)
cbins = numpy.unique(imc)
keep = []
for b in bins:
if log_chirp_width:
for b2 in cbins:
bloc = numpy.where((wtime == b) & (imc == b2))[0]
bloudest = statv[bloc].argsort()[-num_keep:]
keep.append(bloc[bloudest])
else:
bloc = numpy.where((wtime == b))[0]
bloudest = statv[bloc].argsort()[-num_keep:]
keep.append(bloc[bloudest])
keep = numpy.concatenate(keep)
self.events = self.events[keep]
[docs]
def add_template_events(self, columns, vectors):
""" Add a vector indexed """
# initialize with zeros - since vectors can be None, look for the
# longest one that isn't
new_events = None
for v in vectors:
if v is not None:
new_events = numpy.zeros(len(v), dtype=self.event_dtype)
break
# they shouldn't all be None
assert new_events is not None
new_events['template_id'] = self.template_index
for c, v in zip(columns, vectors):
if v is not None:
if isinstance(v, Array):
new_events[c] = v.numpy()
else:
new_events[c] = v
self.template_events = numpy.append(self.template_events, new_events)
[docs]
def cluster_template_events(self, tcolumn, column, window_size):
""" Cluster the internal events over the named column
"""
cvec = self.template_events[column]
tvec = self.template_events[tcolumn]
if window_size == 0:
indices = numpy.arange(len(tvec))
else:
indices = findchirp_cluster_over_window(tvec, cvec, window_size)
self.template_events = numpy.take(self.template_events, indices)
[docs]
def new_template(self, **kwds):
self.template_params.append(kwds)
self.template_index += 1
[docs]
def add_template_params(self, **kwds):
self.template_params[-1].update(kwds)
[docs]
def finalize_template_events(self):
self.accumulate.append(self.template_events)
self.template_events = numpy.array([], dtype=self.event_dtype)
[docs]
def consolidate_events(self, opt, gwstrain=None):
self.events = numpy.concatenate(self.accumulate)
logger.info("We currently have %d triggers", len(self.events))
if opt.chisq_threshold and opt.chisq_bins:
logger.info("Removing triggers with poor chisq")
self.chisq_threshold(opt.chisq_threshold, opt.chisq_bins,
opt.chisq_delta)
logger.info("%d remaining triggers", len(self.events))
if opt.newsnr_threshold and opt.chisq_bins:
logger.info("Removing triggers with NewSNR below threshold")
self.newsnr_threshold(opt.newsnr_threshold)
logger.info("%d remaining triggers", len(self.events))
if opt.keep_loudest_interval:
logger.info("Removing triggers not within the top %s "
"loudest of a %s second interval by %s",
opt.keep_loudest_num, opt.keep_loudest_interval,
opt.keep_loudest_stat)
self.keep_loudest_in_interval\
(opt.keep_loudest_interval * opt.sample_rate,
opt.keep_loudest_num, statname=opt.keep_loudest_stat,
log_chirp_width=opt.keep_loudest_log_chirp_window)
logger.info("%d remaining triggers", len(self.events))
if opt.injection_window and hasattr(gwstrain, 'injections'):
logger.info("Keeping triggers within %s seconds of injection",
opt.injection_window)
self.keep_near_injection(opt.injection_window,
gwstrain.injections)
logger.info("%d remaining triggers", len(self.events))
self.accumulate = [self.events]
[docs]
def finalize_events(self):
self.events = numpy.concatenate(self.accumulate)
[docs]
def make_output_dir(self, outname):
path = os.path.dirname(outname)
if path != '':
if not os.path.exists(path) and path is not None:
os.makedirs(path)
[docs]
def write_events(self, outname):
"""Write the found events to a file. The only currently supported
format is HDF5, indicated by an .hdf or .h5 extension.
"""
self.make_output_dir(outname)
if outname.endswith(('.hdf', '.h5')):
self.write_to_hdf(outname)
return
raise ValueError('Unsupported event output file format')
[docs]
def write_to_hdf(self, outname):
self.events.sort(order='template_id')
th = numpy.array([p['tmplt'].template_hash for p in
self.template_params])
tid = self.events['template_id']
f = H5FileSyntSugar(outname, self.opt.channel_name[0:2])
if len(self.events):
f['snr'] = abs(self.events['snr'])
try:
# Precessing
f['u_vals'] = self.events['u_vals']
f['coa_phase'] = self.events['coa_phase']
f['hplus_cross_corr'] = self.events['hplus_cross_corr']
except Exception:
# Not precessing
f['coa_phase'] = numpy.angle(self.events['snr'])
f['chisq'] = self.events['chisq']
f['bank_chisq'] = self.events['bank_chisq']
f['bank_chisq_dof'] = self.events['bank_chisq_dof']
f['cont_chisq'] = self.events['cont_chisq']
f['end_time'] = self.events['time_index'] / \
float(self.opt.sample_rate) \
+ self.opt.gps_start_time
try:
# Precessing
template_sigmasq_plus = numpy.array(
[t['sigmasq_plus'] for t in self.template_params],
dtype=numpy.float32)
f['sigmasq_plus'] = template_sigmasq_plus[tid]
template_sigmasq_cross = numpy.array(
[t['sigmasq_cross'] for t in self.template_params],
dtype=numpy.float32)
f['sigmasq_cross'] = template_sigmasq_cross[tid]
# FIXME: I want to put something here, but I haven't yet
# figured out what it should be. I think we would also
# need information from the plus and cross correlation
# (both real and imaginary(?)) to get this.
f['sigmasq'] = template_sigmasq_plus[tid]
except Exception:
# Not precessing
f['sigmasq'] = self.events['sigmasq']
# Template durations should ideally be stored in the bank file.
# At present, however, a few plotting/visualization codes
# downstream in the offline search workflow rely on durations being
# stored in the trigger files instead.
template_durations = [p['tmplt'].template_duration for p in
self.template_params]
f['template_duration'] = numpy.array(template_durations,
dtype=numpy.float32)[tid]
# FIXME: Can we get this value from the autochisq instance?
cont_dof = self.opt.autochi_number_points
if self.opt.autochi_onesided is None:
cont_dof = cont_dof * 2
if self.opt.autochi_two_phase:
cont_dof = cont_dof * 2
if self.opt.autochi_max_valued_dof:
cont_dof = self.opt.autochi_max_valued_dof
f['cont_chisq_dof'] = numpy.repeat(cont_dof, len(self.events))
if 'chisq_dof' in self.events.dtype.names:
f['chisq_dof'] = self.events['chisq_dof'] / 2 + 1
else:
f['chisq_dof'] = numpy.zeros(len(self.events))
f['template_hash'] = th[tid]
if 'sg_chisq' in self.events.dtype.names:
f['sg_chisq'] = self.events['sg_chisq']
if self.opt.psdvar_segment is not None:
f['psd_var_val'] = self.events['psd_var_val']
if self.opt.trig_start_time:
f['search/start_time'] = numpy.array([self.opt.trig_start_time])
search_start_time = float(self.opt.trig_start_time)
else:
f['search/start_time'] = numpy.array([self.opt.gps_start_time +
self.opt.segment_start_pad])
search_start_time = float(self.opt.gps_start_time +
self.opt.segment_start_pad)
if self.opt.trig_end_time:
f['search/end_time'] = numpy.array([self.opt.trig_end_time])
search_end_time = float(self.opt.trig_end_time)
else:
f['search/end_time'] = numpy.array([self.opt.gps_end_time -
self.opt.segment_end_pad])
search_end_time = float(self.opt.gps_end_time -
self.opt.segment_end_pad)
if self.write_performance:
self.analysis_time = search_end_time - search_start_time
time_ratio = numpy.array(
[float(self.analysis_time) / float(self.run_time)])
temps_per_core = float(self.ntemplates) / float(self.ncores)
filters_per_core = float(self.nfilters) / float(self.ncores)
f['search/templates_per_core'] = \
numpy.array([float(temps_per_core) * float(time_ratio)])
f['search/filter_rate_per_core'] = \
numpy.array([filters_per_core / float(self.run_time)])
f['search/setup_time_fraction'] = \
numpy.array([float(self.setup_time) / float(self.run_time)])
f['search/run_time'] = numpy.array([float(self.run_time)])
if 'q_trans' in self.global_params:
qtrans = self.global_params['q_trans']
for key in qtrans:
if key == 'qtiles':
for seg in qtrans[key]:
for q in qtrans[key][seg]:
f['qtransform/%s/%s/%s' % (key, seg, q)] = \
qtrans[key][seg][q]
elif key == 'qplanes':
for seg in qtrans[key]:
f['qtransform/%s/%s' % (key, seg)] = qtrans[key][seg]
if 'gating_info' in self.global_params:
gating_info = self.global_params['gating_info']
for gate_type in ['file', 'auto']:
if gate_type in gating_info:
f['gating/' + gate_type + '/time'] = \
numpy.array([float(g[0]) for g in gating_info[gate_type]])
f['gating/' + gate_type + '/width'] = \
numpy.array([g[1] for g in gating_info[gate_type]])
f['gating/' + gate_type + '/pad'] = \
numpy.array([g[2] for g in gating_info[gate_type]])
f.f.close()
class EventManagerMultiDetBase(EventManager):
def __init__(self, opt, ifos, column, column_types, psd=None, **kwargs):
self.opt = opt
self.ifos = ifos
self.global_params = kwargs
if psd is not None:
self.global_params['psd'] = psd[ifos[0]]
# The events array does not like holding the ifo as string,
# so create a mapping dict and hold as an int
self.ifo_dict = {}
self.ifo_reverse = {}
for i, ifo in enumerate(ifos):
self.ifo_dict[ifo] = i
self.ifo_reverse[i] = ifo
self.event_dtype = [('template_id', int), ('event_id', int)]
for col, coltype in zip(column, column_types):
self.event_dtype.append((col, coltype))
self.events = numpy.array([], dtype=self.event_dtype)
self.event_id_map = {}
self.template_params = []
self.template_index = -1
self.template_event_dict = {}
self.coinc_list = []
self.write_performance = False
for ifo in ifos:
self.template_event_dict[ifo] = \
numpy.array([], dtype=self.event_dtype)
def add_template_events_to_ifo(self, ifo, columns, vectors):
""" Add a vector indexed """
# Just call through to the standard function
self.template_events = self.template_event_dict[ifo]
self.add_template_events(columns, vectors)
self.template_event_dict[ifo] = self.template_events
self.template_events = None
def write_gating_info_to_hdf(self, hf):
"""Write per-detector gating information to an h5py file object.
The information is laid out according to the following groups and
datasets: `/<detector>/gating/{file, auto}/{time, width, pad}` where
"file" and "auto" indicate respectively externally-provided gates and
internally-generated gates (autogating), and "time", "width" and "pad"
indicate the gate center times, total durations and padding durations
in seconds respectively.
"""
if 'gating_info' not in self.global_params:
return
gates = self.global_params['gating_info']
for ifo, gate_type in itertools.product(self.ifos, ['file', 'auto']):
if gate_type not in gates[ifo]:
continue
hf[f'{ifo}/gating/{gate_type}/time'] = numpy.array(
[float(g[0]) for g in gates[ifo][gate_type]]
)
hf[f'{ifo}/gating/{gate_type}/width'] = numpy.array(
[g[1] for g in gates[ifo][gate_type]]
)
hf[f'{ifo}/gating/{gate_type}/pad'] = numpy.array(
[g[2] for g in gates[ifo][gate_type]]
)
[docs]
class EventManagerCoherent(EventManagerMultiDetBase):
def __init__(self, opt, ifos, column, column_types, network_column,
network_column_types, segments, time_slides,
psd=None, **kwargs):
super(EventManagerCoherent, self).__init__(
opt, ifos, column, column_types, psd=None, **kwargs)
self.network_event_dtype = \
[(ifo + '_event_id', int) for ifo in self.ifos]
self.network_event_dtype.append(('template_id', int))
self.network_event_dtype.append(('event_id', int))
for col, coltype in zip(network_column, network_column_types):
self.network_event_dtype.append((col, coltype))
self.network_events = numpy.array([], dtype=self.network_event_dtype)
self.event_index = {}
for ifo in self.ifos:
self.event_index[ifo] = 0
self.event_index['network'] = 0
self.template_event_dict['network'] = numpy.array(
[], dtype=self.network_event_dtype)
self.segments = segments
self.time_slides = time_slides
[docs]
def cluster_template_network_events(self, tcolumn, column, window_size,
slide=0):
""" Cluster the internal events over the named column
Parameters
----------------
tcolumn
Indicates which column contains the time.
column
The named column to cluster.
window_size
The size of the window.
slide
Default is 0.
"""
slide_indices = (
self.template_event_dict['network']['slide_id'] == slide)
cvec = self.template_event_dict['network'][column][slide_indices]
tvec = self.template_event_dict['network'][tcolumn][slide_indices]
if not window_size == 0:
# cluster events over the window
indices = findchirp_cluster_over_window(tvec, cvec, window_size)
# if a slide_indices = 0
if any(~slide_indices):
indices = numpy.concatenate((
numpy.flatnonzero(~slide_indices),
numpy.flatnonzero(slide_indices)[indices]))
indices.sort()
# get value of key for where you have indicies
for key in self.template_event_dict:
self.template_event_dict[key] = \
self.template_event_dict[key][indices]
else:
indices = numpy.arange(len(tvec))
[docs]
def add_template_network_events(self, columns, vectors):
""" Add a vector indexed """
# initialize with zeros - since vectors can be None, look for the
# longest one that isn't
new_events = numpy.zeros(
max([len(v) for v in vectors if v is not None]),
dtype=self.network_event_dtype
)
# they shouldn't all be None
assert new_events is not None
new_events['template_id'] = self.template_index
for c, v in zip(columns, vectors):
if v is not None:
if isinstance(v, Array):
new_events[c] = v.numpy()
else:
new_events[c] = v
self.template_events = numpy.append(self.template_events, new_events)
[docs]
def add_template_events_to_network(self, columns, vectors):
""" Add a vector indexed """
# Just call through to the standard function
self.template_events = self.template_event_dict['network']
self.add_template_network_events(columns, vectors)
self.template_event_dict['network'] = self.template_events
self.template_events = None
[docs]
def write_to_hdf(self, outname):
self.events.sort(order='template_id')
th = numpy.array(
[p['tmplt'].template_hash for p in self.template_params])
f = H5FileSyntSugar(outname)
self.write_gating_info_to_hdf(f)
# Output network stuff
f.prefix = 'network'
network_events = numpy.array(
[e for e in self.network_events], dtype=self.network_event_dtype)
for col in network_events.dtype.names:
if col == 'time_index':
f['end_time_gc'] = (
network_events[col]
/ float(self.opt.sample_rate[self.ifos[0].lower()])
+ self.opt.gps_start_time[self.ifos[0].lower()]
)
else:
f[col] = network_events[col]
starts = []
ends = []
for seg in self.segments[self.ifos[0]]:
starts.append(seg.start_time.gpsSeconds)
ends.append(seg.end_time.gpsSeconds)
f['search/segments/start_times'] = starts
f['search/segments/end_times'] = ends
# Individual ifo stuff
for i, ifo in enumerate(self.ifos):
tid = self.events['template_id'][self.events['ifo'] == i]
f.prefix = ifo
ifo_events = numpy.array([e for e in self.events
if e['ifo'] == self.ifo_dict[ifo]], dtype=self.event_dtype)
if len(ifo_events):
f['snr'] = abs(ifo_events['snr'])
f['event_id'] = ifo_events['event_id']
try:
# Precessing
f['u_vals'] = ifo_events['u_vals']
f['coa_phase'] = ifo_events['coa_phase']
f['hplus_cross_corr'] = ifo_events['hplus_cross_corr']
except Exception:
f['coa_phase'] = numpy.angle(ifo_events['snr'])
f['chisq'] = ifo_events['chisq']
f['bank_chisq'] = ifo_events['bank_chisq']
f['bank_chisq_dof'] = ifo_events['bank_chisq_dof']
f['auto_chisq'] = ifo_events['auto_chisq']
f['auto_chisq_dof'] = ifo_events['auto_chisq_dof']
f['end_time'] = ifo_events['time_index'] / \
float(self.opt.sample_rate[ifo]) + \
self.opt.gps_start_time[ifo]
f['time_index'] = ifo_events['time_index']
f['slide_id'] = ifo_events['slide_id']
try:
# Precessing
template_sigmasq_plus = numpy.array(
[t['sigmasq_plus'] for t in self.template_params],
dtype=numpy.float32
)
f['sigmasq_plus'] = template_sigmasq_plus[tid]
template_sigmasq_cross = numpy.array(
[t['sigmasq_cross'] for t in self.template_params],
dtype=numpy.float32
)
f['sigmasq_cross'] = template_sigmasq_cross[tid]
# FIXME: I want to put something here, but I haven't yet
# figured out what it should be. I think we would also
# need information from the plus and cross correlation
# (both real and imaginary(?)) to get this.
f['sigmasq'] = template_sigmasq_plus[tid]
except Exception:
# Not precessing
template_sigmasq = numpy.array(
[t['sigmasq'][ifo] for t in self.template_params],
dtype=numpy.float32)
f['sigmasq'] = template_sigmasq[tid]
# FIXME: Can we get this value from the autochisq instance?
# cont_dof = self.opt.autochi_number_points
# if self.opt.autochi_onesided is None:
# cont_dof = cont_dof * 2
# if self.opt.autochi_two_phase:
# cont_dof = cont_dof * 2
# if self.opt.autochi_max_valued_dof:
# cont_dof = self.opt.autochi_max_valued_dof
# f['cont_chisq_dof'] = numpy.repeat(cont_dof, len(ifo_events))
if 'chisq_dof' in ifo_events.dtype.names:
f['chisq_dof'] = ifo_events['chisq_dof'] / 2 + 1
else:
f['chisq_dof'] = numpy.zeros(len(ifo_events))
f['template_hash'] = th[tid]
f['search/time_slides'] = numpy.array(self.time_slides[ifo])
if self.opt.trig_start_time:
f['search/start_time'] = numpy.array([
self.opt.trig_start_time[ifo]], dtype=numpy.int32)
search_start_time = float(self.opt.trig_start_time[ifo])
else:
f['search/start_time'] = numpy.array([
self.opt.gps_start_time[ifo] +
self.opt.segment_start_pad[ifo]], dtype=numpy.int32)
search_start_time = float(self.opt.gps_start_time[ifo] +
self.opt.segment_start_pad[ifo])
if self.opt.trig_end_time:
f['search/end_time'] = numpy.array([
self.opt.trig_end_time[ifo]], dtype=numpy.int32)
search_end_time = float(self.opt.trig_end_time[ifo])
else:
f['search/end_time'] = numpy.array(
[self.opt.gps_end_time[ifo] -
self.opt.segment_end_pad[ifo]], dtype=numpy.int32)
search_end_time = float(self.opt.gps_end_time[ifo] -
self.opt.segment_end_pad[ifo])
if self.write_performance:
self.analysis_time = search_end_time - search_start_time
time_ratio = numpy.array([float(self.analysis_time) /
float(self.run_time)])
temps_per_core = float(self.ntemplates) / float(self.ncores)
filters_per_core = float(self.nfilters) / float(self.ncores)
f['search/templates_per_core'] = \
numpy.array([float(temps_per_core) * float(time_ratio)])
f['search/filter_rate_per_core'] = \
numpy.array([filters_per_core / float(self.run_time)])
f['search/setup_time_fraction'] = \
numpy.array([float(self.setup_time) / float(self.run_time)])
[docs]
def finalize_template_events(self):
# Check that none of the template events have the same time index as an
# existing event in events. I.e. don't list the same ifo event multiple
# times when looping over sky points and time slides.
existing_times = {}
new_times = {}
existing_template_id = {}
new_template_id = {}
existing_events_mask = {}
new_template_event_mask = {}
existing_template_event_mask = {}
for i, ifo in enumerate(self.ifos):
ifo_events = numpy.where(self.events['ifo'] == i)
existing_times[ifo] = self.events['time_index'][ifo_events]
new_times[ifo] = self.template_event_dict[ifo]['time_index']
existing_template_id[ifo] = self.events['template_id'][ifo_events]
new_template_id[ifo] = self.template_event_dict[ifo]['template_id']
# This is true for each existing event that has the same time index
# and template id as a template trigger.
existing_events_mask[ifo] = numpy.argwhere(
numpy.logical_and(
numpy.isin(existing_times[ifo], new_times[ifo]),
numpy.isin(existing_template_id[ifo], new_template_id[ifo])
)).reshape(-1,)
# This is true for each template event that has either a new
# trigger time or a new template id.
new_template_event_mask[ifo] = numpy.argwhere(
numpy.logical_or(
~numpy.isin(new_times[ifo], existing_times[ifo]),
~numpy.isin(new_template_id[ifo], existing_template_id[ifo])
)).reshape(-1,)
# This is true for each template event that has the same time index
# and template id as an exisitng event trigger.
existing_template_event_mask[ifo] = numpy.argwhere(
numpy.logical_and(
numpy.isin(new_times[ifo], existing_times[ifo]),
numpy.isin(new_template_id[ifo], existing_template_id[ifo])
)).reshape(-1,)
# Set ids (These show how each trigger in the single ifo trigger
# list correspond to the network triggers)
num_events = len(new_template_event_mask[ifo])
new_event_ids = numpy.arange(self.event_index[ifo],
self.event_index[ifo] + num_events)
# Every template event that corresponds to a new trigger gets a new
# id. Triggers that have been found before are not saved.
self.template_event_dict[ifo]['event_id'][
new_template_event_mask[ifo]] = new_event_ids
self.template_event_dict['network'][ifo + '_event_id'][
new_template_event_mask[ifo]] = new_event_ids
# Template events that have been found before get the event id of
# the first time they were found.
self.template_event_dict['network'][ifo + '_event_id'][
existing_template_event_mask[ifo]] = \
self.events[self.events['ifo'] == i][
existing_events_mask[ifo]]['event_id']
self.event_index[ifo] = self.event_index[ifo] + num_events
# Add the network event ids for the events with this template.
num_events = len(self.template_event_dict['network'])
new_event_ids = numpy.arange(self.event_index['network'],
self.event_index['network'] + num_events)
self.event_index['network'] = self.event_index['network'] + num_events
self.template_event_dict['network']['event_id'] = new_event_ids
# Move template events for each ifo to the events list
for ifo in self.ifos:
self.events = numpy.append(
self.events,
self.template_event_dict[ifo][new_template_event_mask[ifo]]
)
self.template_event_dict[ifo] = \
numpy.array([], dtype=self.event_dtype)
# Move the template events for the network to the network events list
self.network_events = numpy.append(self.network_events,
self.template_event_dict['network'])
self.template_event_dict['network'] = \
numpy.array([], dtype=self.network_event_dtype)
[docs]
class EventManagerMultiDet(EventManagerMultiDetBase):
def __init__(self, opt, ifos, column, column_types, psd=None, **kwargs):
super(EventManagerMultiDet, self).__init__(
opt, ifos, column, column_types, psd=None, **kwargs)
self.event_index = 0
[docs]
def cluster_template_events_single_ifo(
self, tcolumn, column, window_size, ifo):
""" Cluster the internal events over the named column
"""
# Just call through to the standard function
self.template_events = self.template_event_dict[ifo]
self.cluster_template_events(tcolumn, column, window_size)
self.template_event_dict[ifo] = self.template_events
self.template_events = None
[docs]
def finalize_template_events(self, perform_coincidence=True,
coinc_window=0.0):
# Set ids
for ifo in self.ifos:
num_events = len(self.template_event_dict[ifo])
new_event_ids = numpy.arange(self.event_index,
self.event_index+num_events)
self.template_event_dict[ifo]['event_id'] = new_event_ids
self.event_index = self.event_index+num_events
if perform_coincidence:
if not len(self.ifos) == 2:
err_msg = "Coincidence currently only supported for 2 ifos."
raise ValueError(err_msg)
ifo1 = self.ifos[0]
ifo2 = self.ifos[1]
end_times1 = self.template_event_dict[ifo1]['time_index'] /\
float(self.opt.sample_rate[ifo1]) + self.opt.gps_start_time[ifo1]
end_times2 = self.template_event_dict[ifo2]['time_index'] /\
float(self.opt.sample_rate[ifo2]) + self.opt.gps_start_time[ifo2]
light_travel_time = Detector(ifo1).light_travel_time_to_detector(
Detector(ifo2))
coinc_window = coinc_window + light_travel_time
# FIXME: Remove!!!
coinc_window = 2.0
if len(end_times1) and len(end_times2):
idx_list1, idx_list2, _ = \
coinc.time_coincidence(end_times1, end_times2,
coinc_window)
if len(idx_list1):
for idx1, idx2 in zip(idx_list1, idx_list2):
event1 = self.template_event_dict[ifo1][idx1]
event2 = self.template_event_dict[ifo2][idx2]
self.coinc_list.append((event1, event2))
for ifo in self.ifos:
self.events = numpy.append(self.events,
self.template_event_dict[ifo])
self.template_event_dict[ifo] = numpy.array([],
dtype=self.event_dtype)
[docs]
def write_to_hdf(self, outname):
self.events.sort(order='template_id')
th = numpy.array([p['tmplt'].template_hash for p in
self.template_params])
tid = self.events['template_id']
f = H5FileSyntSugar(outname)
self.write_gating_info_to_hdf(f)
for ifo in self.ifos:
f.prefix = ifo
ifo_events = numpy.array([e for e in self.events if
e['ifo'] == self.ifo_dict[ifo]],
dtype=self.event_dtype)
if len(ifo_events):
f['snr'] = abs(ifo_events['snr'])
try:
# Precessing
f['u_vals'] = ifo_events['u_vals']
f['coa_phase'] = ifo_events['coa_phase']
f['hplus_cross_corr'] = ifo_events['hplus_cross_corr']
except Exception:
f['coa_phase'] = numpy.angle(ifo_events['snr'])
f['chisq'] = ifo_events['chisq']
f['bank_chisq'] = ifo_events['bank_chisq']
f['bank_chisq_dof'] = ifo_events['bank_chisq_dof']
f['cont_chisq'] = ifo_events['cont_chisq']
f['end_time'] = ifo_events['time_index'] / \
float(self.opt.sample_rate[ifo]) + \
self.opt.gps_start_time[ifo]
try:
# Precessing
template_sigmasq_plus = numpy.array([t['sigmasq_plus'] for
t in self.template_params], dtype=numpy.float32)
f['sigmasq_plus'] = template_sigmasq_plus[tid]
template_sigmasq_cross = numpy.array([t['sigmasq_cross']
for t in self.template_params], dtype=numpy.float32)
f['sigmasq_cross'] = template_sigmasq_cross[tid]
# FIXME: I want to put something here, but I haven't yet
# figured out what it should be. I think we would also
# need information from the plus and cross correlation
# (both real and imaginary(?)) to get this.
f['sigmasq'] = template_sigmasq_plus[tid]
except Exception:
# Not precessing
template_sigmasq = numpy.array([t['sigmasq'][ifo] for t in
self.template_params],
dtype=numpy.float32)
f['sigmasq'] = template_sigmasq[tid]
# FIXME: Can we get this value from the autochisq instance?
cont_dof = self.opt.autochi_number_points
if self.opt.autochi_onesided is None:
cont_dof = cont_dof * 2
# if self.opt.autochi_two_phase:
# cont_dof = cont_dof * 2
# if self.opt.autochi_max_valued_dof:
# cont_dof = self.opt.autochi_max_valued_dof
f['cont_chisq_dof'] = numpy.repeat(cont_dof, len(ifo_events))
if 'chisq_dof' in ifo_events.dtype.names:
f['chisq_dof'] = ifo_events['chisq_dof'] / 2 + 1
else:
f['chisq_dof'] = numpy.zeros(len(ifo_events))
f['template_hash'] = th[tid]
if self.opt.psdvar_segment is not None:
f['psd_var_val'] = ifo_events['psd_var_val']
if self.opt.trig_start_time:
f['search/start_time'] = numpy.array(
[self.opt.trig_start_time[ifo]], dtype=numpy.int32)
search_start_time = float(self.opt.trig_start_time[ifo])
else:
f['search/start_time'] = numpy.array(
[self.opt.gps_start_time[ifo] +
self.opt.segment_start_pad[ifo]], dtype=numpy.int32)
search_start_time = float(self.opt.gps_start_time[ifo] +
self.opt.segment_start_pad[ifo])
if self.opt.trig_end_time:
f['search/end_time'] = numpy.array(
[self.opt.trig_end_time[ifo]], dtype=numpy.int32)
search_end_time = float(self.opt.trig_end_time[ifo])
else:
f['search/end_time'] = numpy.array(
[self.opt.gps_end_time[ifo] -
self.opt.segment_end_pad[ifo]], dtype=numpy.int32)
search_end_time = float(self.opt.gps_end_time[ifo] -
self.opt.segment_end_pad[ifo])
if self.write_performance:
self.analysis_time = search_end_time - search_start_time
time_ratio = numpy.array(
[float(self.analysis_time) / float(self.run_time)])
temps_per_core = float(self.ntemplates) / float(self.ncores)
filters_per_core = float(self.nfilters) / float(self.ncores)
f['search/templates_per_core'] = \
numpy.array([float(temps_per_core) * float(time_ratio)])
f['search/filter_rate_per_core'] = \
numpy.array([filters_per_core / float(self.run_time)])
f['search/setup_time_fraction'] = \
numpy.array([float(self.setup_time) / float(self.run_time)])
__all__ = ['threshold_and_cluster', 'findchirp_cluster_over_window',
'threshold', 'cluster_reduce', 'ThresholdCluster',
'threshold_real_numpy', 'threshold_only',
'EventManager', 'EventManagerMultiDet', 'EventManagerCoherent']