Commit 6f142fdc authored by Vivien Raymond's avatar Vivien Raymond

Merge branch 'remove-bayestar-cruft' into 'master'

Remove old unmaintained versions of BAYESTAR-related files

See merge request !647
parents a3fce044 9a412584
Pipeline #49569 passed with stages
in 83 minutes and 47 seconds
......@@ -57,11 +57,8 @@ python/lalinference/_bayespputils.so
python/lalinference/_distance.so
python/lalinference/_lalinference*
python/lalinference/_moc.so
python/lalinference/bayestar/_sky_map.so
python/lalinference/bayestar/sky_map.so
python/lalinference/git_version.py
python/lalinference/lalinference.py
python/lalinference_average_skymaps
python/lalinference_burst_pp_pipe
python/lalinference_coherence_test
python/lalinference_compute_roq_weights
......
......@@ -20,7 +20,6 @@ AC_CONFIG_FILES([ \
python/lalinference/bayestar/Makefile \
python/lalinference/imrtgr/Makefile \
python/lalinference/io/Makefile \
python/lalinference/io/events/Makefile \
python/lalinference/plot/Makefile \
python/lalinference/popprior/Makefile \
python/lalinference/rapid_pe/Makefile \
......
usr/bin/cbcBayes*
usr/bin/imrtgr_*
usr/bin/lalinference_average_skymaps
usr/bin/lalinference_burst_pp_pipe
usr/bin/lalinference_coherence_test
usr/bin/lalinference_compute_roq_weights
......
......@@ -182,7 +182,6 @@ rm -Rf ${RPM_BUILD_DIR}/%{name}-%{version}%{?nightly:-%{nightly}}
%license COPYING
%{_bindir}/cbcBayes*
%{_bindir}/imrtgr_*
%{_bindir}/lalinference_average_skymaps
%{_bindir}/lalinference_burst_pp_pipe
%{_bindir}/lalinference_coherence_test
%{_bindir}/lalinference_compute_roq_weights
......
......@@ -15,7 +15,6 @@ pybin_scripts = \
rapidpe_compute_intrinsic_grid \
rapidpe_calculate_overlap \
imrtgr_imr_consistency_test \
lalinference_average_skymaps \
lalinference_burst_pp_pipe \
lalinference_coherence_test \
lalinference_compute_roq_weights \
......
......@@ -9,31 +9,8 @@ pymoduledir = $(pkgpythondir)/bayestar
pymodule_PYTHON = \
__init__.py \
command.py \
decorator.py \
deprecation.py \
filter.py \
ligolw.py \
sky_map.py \
postprocess.py \
timing.py \
$(END_OF_LIST)
if HAVE_CHEALPIX
if SWIG_BUILD_PYTHON
pymodule_LTLIBRARIES = _sky_map.la
_sky_map_la_CPPFLAGS = $(AM_CPPFLAGS) $(SWIG_PYTHON_CPPFLAGS) -I$(top_srcdir)/src
_sky_map_la_CFLAGS = $(AM_CFLAGS) $(SWIG_PYTHON_CFLAGS) -Wno-error
_sky_map_la_LDFLAGS = $(top_builddir)/src/liblalinference.la -shared -module -avoid-version
endif
endif
all-local: _sky_map.so
_sky_map.so:
rm -f $@ && $(LN_S) .libs/$@
CLEANFILES = _sky_map.so
endif
This diff is collapsed.
This diff is collapsed.
#
# Copyright (C) 2013-2017 Leo Singer
#
# 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.
#
"""
Collection of Python decorators.
"""
from functools import wraps
from collections import Hashable
from astropy.utils.misc import NumpyRNGContext
__all__ = ('memoized', 'with_numpy_random_seed', 'as_dict')
try:
from functools import lru_cache
except ImportError:
# FIXME: Remove this when we drop support for Python < 3.2.
def memoized(func):
"""Memoize a function or class by caching its return values for any
given arguments."""
cache = {}
@wraps(func)
def memo(*args, **kwargs):
# Create a key out of the arguments.
key = (args, frozenset(kwargs.items()))
if isinstance(args, Hashable): # The key is immutable.
try:
# Look up the return value for these arguments.
ret = cache[key]
except KeyError:
# Not found; invoke function and store return value.
ret = cache[key] = func(*args, **kwargs)
else: # The key is mutable. We can't cache it.
ret = func(*args, **kwargs)
# Done!
return ret
# Return wrapped function.
return memo
else:
memoized = lru_cache(maxsize=None)
def with_numpy_random_seed(func, seed=0):
"""Decorate a function so that it is called with a pre-defined random seed.
The random seed is restored when the function returns."""
@wraps(func)
def wrapped_func(*args, **kwargs):
with NumpyRNGContext(seed):
ret = func(*args, **kwargs)
return ret
return wrapped_func
def as_dict(func):
"""Decorate a generator function to turn its output into a dict."""
@wraps(func)
def wrapped_func(*args, **kwargs):
return dict(func(*args, **kwargs))
return wrapped_func
This diff is collapsed.
#
# Copyright (C) 2013-2017 Leo Singer
#
# 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.
#
"""
LIGO-LW convenience functions.
"""
# LIGO-LW XML imports.
from glue.ligolw.utils import process as ligolw_process
from glue.ligolw import ligolw
from glue.ligolw import array as ligolw_array
from glue.ligolw import param as ligolw_param
from glue.ligolw import table as ligolw_table
from glue.ligolw import lsctables
from lalinspiral.inspinjfind import InspiralSCExactCoincDef
def get_template_bank_f_low(xmldoc):
"""Determine the low frequency cutoff from a template bank file,
whether the template bank was produced by lalapps_tmpltbank or
lalapps_cbc_sbank. bayestar_sim_to_tmpltbank does not have a command
line for the low-frequency cutoff; instead, it is recorded in a row
of the search_summvars table."""
try:
template_bank_f_low, = ligolw_process.get_process_params(
xmldoc, 'tmpltbank', '--low-frequency-cutoff')
except ValueError:
try:
template_bank_f_low, = ligolw_process.get_process_params(
xmldoc, 'lalapps_cbc_sbank', '--flow')
except ValueError:
try:
search_summvars_table = ligolw_table.get_table(
xmldoc, lsctables.SearchSummVarsTable.tableName)
template_bank_f_low, = (
search_summvars.value
for search_summvars in search_summvars_table
if search_summvars.name == 'low-frequency cutoff')
except ValueError:
raise ValueError("Could not determine low-frequency cutoff")
return template_bank_f_low
def sim_coinc_and_sngl_inspirals_for_xmldoc(xmldoc):
"""Retrieve (as a generator) all of the
(sim_inspiral, coinc_event, (sngl_inspiral, sngl_inspiral,
... sngl_inspiral) tuples from found coincidences in a LIGO-LW XML
document."""
# Look up necessary tables.
coinc_table = ligolw_table.get_table(
xmldoc, lsctables.CoincTable.tableName)
coinc_def_table = ligolw_table.get_table(
xmldoc, lsctables.CoincDefTable.tableName)
coinc_map_table = ligolw_table.get_table(
xmldoc, lsctables.CoincMapTable.tableName)
# Look up coinc_def ids.
sim_coinc_def_id = coinc_def_table.get_coinc_def_id(
InspiralSCExactCoincDef.search,
InspiralSCExactCoincDef.search_coinc_type,
create_new=False)
def events_for_coinc_event_id(coinc_event_id):
for coinc_map in coinc_map_table:
if coinc_map.coinc_event_id == coinc_event_id:
for row in ligolw_table.get_table(
xmldoc, coinc_map.table_name):
column_name = coinc_map.event_id.column_name
if getattr(row, column_name) == coinc_map.event_id:
yield coinc_map.event_id, row
# Loop over all coinc_event <-> sim_inspiral coincs.
for sim_coinc in coinc_table:
# If this is not a coinc_event <-> sim_inspiral coinc, skip it.
if sim_coinc.coinc_def_id != sim_coinc_def_id:
continue
# Locate the sim_inspiral and coinc events.
sim_inspiral = None
coinc = None
for event_id, event in events_for_coinc_event_id(
sim_coinc.coinc_event_id):
if event_id.table_name == ligolw_table.Table.TableName(
lsctables.SimInspiralTable.tableName):
if sim_inspiral is not None:
raise RuntimeError(
'Found more than one matching sim_inspiral entry')
sim_inspiral = event
elif event_id.table_name == ligolw_table.Table.TableName(
lsctables.CoincTable.tableName):
if coinc is not None:
raise RuntimeError(
'Found more than one matching coinc entry')
coinc = event
else:
raise RuntimeError(
"Did not expect coincidence to contain an event of "
"type '{}'".format(event_id.table_name))
sngl_inspirals = tuple(
event for event_id, event in events_for_coinc_event_id(
coinc.coinc_event_id))
yield sim_inspiral, coinc, sngl_inspirals
@lsctables.use_in
class LSCTablesContentHandler(ligolw.LIGOLWContentHandler):
"""Content handler for reading LIGO-LW XML files with LSC table schema."""
@ligolw_array.use_in
@ligolw_param.use_in
@lsctables.use_in
class LSCTablesAndSeriesContentHandler(ligolw.LIGOLWContentHandler):
"""Content handler for reading LIGO-LW XML files with LSC table schema and
gridded arrays."""
# -*- coding: utf-8
#
# Copyright (C) 2013-2015 Leo Singer
#
# 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.
#
from __future__ import division
"""
Functions for predicting timing accuracy of matched filters.
"""
import logging
import numpy as np
import lal
import lalsimulation
from scipy import interpolate
from scipy import linalg
from lalinference.bayestar.deprecation import warn
warn('ligo.skymap.bayestar.filter')
log = logging.getLogger('BAYESTAR')
_noise_psd_funcs = {}
class vectorize_swig_psd_func(object):
"""Create a vectorized Numpy function from a SWIG-wrapped PSD function.
SWIG does not provide enough information for Numpy to determine the number
of input arguments, so we can't just use np.vectorize."""
def __init__(self, str):
self.__func = getattr(lalsimulation, str + 'Ptr')
self.__npyfunc = np.frompyfunc(getattr(lalsimulation, str), 1, 1)
def __call__(self, f):
fa = np.asarray(f)
df = np.diff(fa)
if fa.ndim == 1 and df.size > 1 and np.all(df[0] == df[1:]):
fa = np.concatenate((fa, [fa[-1] + df[0]]))
ret = lal.CreateREAL8FrequencySeries(
None, 0, fa[0], df[0], lal.DimensionlessUnit, fa.size)
lalsimulation.SimNoisePSD(ret, 0, self.__func)
ret = ret.data.data[:-1]
else:
ret = self.__npyfunc(f)
if not np.isscalar(ret):
ret = ret.astype(float)
return ret
for _ifos, _func in (
(("H1", "H2", "L1", "I1"), 'SimNoisePSDaLIGOZeroDetHighPower'),
(("V1",), 'SimNoisePSDAdvVirgo'),
(("K1"), 'SimNoisePSDKAGRA')
):
_func = vectorize_swig_psd_func(_func)
for _ifo in _ifos:
_noise_psd_funcs[_ifo] = _func
def get_noise_psd_func(ifo):
"""Find a function that describes the given interferometer's noise PSD."""
return _noise_psd_funcs[ifo]
class InterpolatedPSD(interpolate.interp1d):
"""Create a (linear in log-log) interpolating function for a discretely
sampled power spectrum S(f)."""
def __init__(self, f, S, f_high_truncate=1.0):
assert f_high_truncate <= 1.0
f = np.asarray(f)
S = np.asarray(S)
# Exclude DC if present
if f[0] == 0:
f = f[1:]
S = S[1:]
# FIXME: This is a hack to fix an issue with the detection pipeline's
# PSD conditioning. Remove this when the issue is fixed upstream.
if f_high_truncate < 1.0:
log.warn(
'Truncating PSD at %g of maximum frequency to suppress '
'rolloff artifacts. This option may be removed in the future.',
f_high_truncate)
keep = (f <= f_high_truncate * max(f))
f = f[keep]
S = S[keep]
super(InterpolatedPSD, self).__init__(
np.log(f), np.log(S),
kind='linear', bounds_error=False, fill_value=np.inf)
self._f_min = min(f)
self._f_max = max(f)
def __call__(self, f):
f_min = np.min(f)
f_max = np.max(f)
if f_min < self._f_min:
log.warn('Assuming PSD is infinite at %g Hz because PSD is only '
'sampled down to %g Hz', f_min, self._f_min)
if f_max > self._f_max:
log.warn('Assuming PSD is infinite at %g Hz because PSD is only '
'sampled up to %g Hz', f_max, self._f_max)
return np.where(
(f >= self._f_min) & (f <= self._f_max),
np.exp(super(InterpolatedPSD, self).__call__(np.log(f))), np.inf)
class SignalModel(object):
"""Class to speed up computation of signal/noise-weighted integrals and
Barankin and Cramér-Rao lower bounds on time and phase estimation.
Note that the autocorrelation series and the moments are related,
as shown below.
Create signal model:
>>> from . import filter
>>> sngl = lambda: None
>>> H = filter.sngl_inspiral_psd(
... 'TaylorF2threePointFivePN', mass1=1.4, mass2=1.4)
>>> S = get_noise_psd_func('H1')
>>> W = filter.signal_psd_series(H, S)
>>> sm = SignalModel(W)
Compute one-sided autocorrelation function:
>>> out_duration = 0.1
>>> a, sample_rate = filter.autocorrelation(W, out_duration)
Restore negative time lags using symmetry:
>>> a = np.concatenate((a[:0:-1].conj(), a))
Compute the first 2 frequency moments by taking derivatives of the
autocorrelation sequence using centered finite differences.
The nth frequency moment should be given by (-1j)^n a^(n)(t).
>>> acor_moments = []
>>> for i in range(2):
... acor_moments.append(a[len(a) // 2])
... a = -0.5j * sample_rate * (a[2:] - a[:-2])
>>> assert np.all(np.isreal(acor_moments))
>>> acor_moments = np.real(acor_moments)
Compute the first 2 frequency moments using this class.
>>> quad_moments = [sm.get_sn_moment(i) for i in range(2)]
Compare them.
>>> for i, (am, qm) in enumerate(zip(acor_moments, quad_moments)):
... assert np.allclose(am, qm, rtol=0.05)
"""
def __init__(self, h):
"""Create a TaylorF2 signal model with the given masses, PSD function
S(f), PN amplitude order, and low-frequency cutoff."""
# Find indices of first and last nonzero samples.
nonzero = np.flatnonzero(h.data.data)
first_nonzero = nonzero[0]
last_nonzero = nonzero[-1]
# Frequency sample points
self.dw = 2 * np.pi * h.deltaF
f = h.f0 + h.deltaF * np.arange(first_nonzero, last_nonzero + 1)
self.w = 2 * np.pi * f
# Throw away leading and trailing zeros.
h = h.data.data[first_nonzero:last_nonzero + 1]
self.denom_integrand = 4 / (2 * np.pi) * h
self.den = np.trapz(self.denom_integrand, dx=self.dw)
def get_horizon_distance(self, snr_thresh=1):
return np.sqrt(self.den) / snr_thresh
def get_sn_average(self, func):
"""Get the average of a function of angular frequency, weighted by the
signal to noise per unit angular frequency."""
num = np.trapz(func(self.w) * self.denom_integrand, dx=self.dw)
return num / self.den
def get_sn_moment(self, power):
"""Get the average of angular frequency to the given power, weighted by
the signal to noise per unit frequency."""
return self.get_sn_average(lambda w: w**power)
def get_crb(self, snr):
"""Get the Cramér-Rao bound, or inverse Fisher information matrix,
describing the phase and time estimation covariance."""
w1 = self.get_sn_moment(1)
w2 = self.get_sn_moment(2)
I = np.asarray(((1, -w1), (-w1, w2)))
return linalg.inv(I) / np.square(snr)
# FIXME: np.vectorize doesn't work on unbound instance methods. The
# excluded keyword, added in Numpy 1.7, could be used here to exclude the
# zeroth argument, self.
def __get_crb_toa_uncert(self, snr):
return np.sqrt(self.get_crb(snr)[1, 1])
def get_crb_toa_uncert(self, snr):
return np.frompyfunc(self.__get_crb_toa_uncert, 1, 1)(snr)
......@@ -3,8 +3,6 @@ MOSTLYCLEANFILES =
EXTRA_DIST =
include $(top_srcdir)/gnuscripts/lalsuite_python.am
SUBDIRS = events
if HAVE_PYTHON
pymoduledir = $(pkgpythondir)/io
......
BUILT_SOURCES =
MOSTLYCLEANFILES =
EXTRA_DIST =
include $(top_srcdir)/gnuscripts/lalsuite_python.am
if HAVE_PYTHON
pymoduledir = $(pkgpythondir)/io/events
pymodule_PYTHON = \
__init__.py \
base.py \
detector_disabled.py \
ligolw.py \
gracedb.py \
hdf.py \
magic.py \
sqlite.py
endif
# Copyright (C) 2017-2018 Leo Singer
#
# 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.
#
from __future__ import absolute_import
import os
import pkgutil
import six
__all__ = ()
# Import all symbols from all submodules of this module.
for _, module, _ in pkgutil.iter_modules([os.path.dirname(__file__)]):
six.exec_('from . import {0};'
'__all__ += getattr({0}, "__all__", ());'
'from .{0} import *'.format(module))
del module
# Clean up
del os, pkgutil, six
from lalinference.bayestar.deprecation import warn
warn('ligo.skymap.io.events')
# Copyright (C) 2017 Leo Singer
#
# 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.
#
"""
Base classes for reading events from search pipelines.
"""
from abc import ABCMeta, abstractproperty
from collections import Mapping
import six
__all__ = ('EventSource', 'Event', 'SingleEvent')
def _fmt(obj, keys):
kvs = ', '.join('{}={!r}'.format(key, getattr(obj, key)) for key in keys)
return '<{}({})>'.format(obj.__class__.__name__, kvs)
class EventSource(Mapping):
"""Abstraction of a source of coincident events."""
def __str__(self):
try:
length = len(self)
except NotImplementedError:
contents = '...'
else:
contents = '...{} items...'.format(length)
return '<{}({{{}}})>'.format(self.__class__.__name__, contents)
def __repr__(self):
try:
len(self)
except NotImplementedError:
contents = '...'
else:
contents = ', '.join('{}: {!r}'.format(key, value)
for key, value in self.items())
return '{}({{{}}})'.format(self.__class__.__name__, contents)
class Event(six.with_metaclass(ABCMeta)):
"""Abstraction of a coincident trigger."""
@abstractproperty
def singles(self):
"""Tuple of single-detector triggers."""
raise NotImplementedError
@abstractproperty
def template_args(self):
"""Dictionary of template parameters."""
raise NotImplementedError
__str_keys = ('singles',)
def __str__(self):
return _fmt(self, self.__str_keys)
__repr__ = __str__
class SingleEvent(six.with_metaclass(ABCMeta)):
"""Abstraction of a single-detector trigger."""
@abstractproperty
def detector(self):
"""Instrument name (e.g. 'H1')"""
raise NotImplementedError
@abstractproperty
def snr(self):
"""Signal to noise ratio (float)"""
raise NotImplementedError
@abstractproperty
def phase(self):
"""Phase on arrival (float)"""
raise NotImplementedError
@abstractproperty
def time(self):
"""Time on arrival (float, GPS)"""
raise NotImplementedError
@abstractproperty
def zerolag_time(self):
"""Time on arrival in zero-lag data, without time slides applied
(float, GPS)"""
raise NotImplementedError
@abstractproperty
def psd(self):
"""PSD (REAL8FrequencySeries)"""
raise NotImplementedError
@property
def snr_series(self):
"""SNR time series (COMPLEX8TimeSeries)"""
return None