Commit ffa02461 authored by Leo Pound Singer's avatar Leo Pound Singer

Add bayestar-localization to build

- Minimal changes to lalapps build
- cHEALPix now an optional dependency of lalinference
- OpenMP acceleration now optional in lalinference
- Building of Python C extension conditional on Python support
Original: b29f20ff11e62095dbd44e850b248ecc58b08a13
parent e5fa0a53
......@@ -122,6 +122,17 @@ src/inspiral/lalapps_trigbank
src/inspiral/lalapps_trigger_hipe
src/inspiral/lalapps_trigscan
src/inspiral/lalapps_write_ihope_page
src/inspiral/bayestar/bayestar_aggregate_found_injections
src/inspiral/bayestar/bayestar_lattice_tmpltbank
src/inspiral/bayestar/bayestar_littlehope
src/inspiral/bayestar/bayestar_localize_coincs
src/inspiral/bayestar/bayestar_localize_gracedb
src/inspiral/bayestar/bayestar_localize_lvalert
src/inspiral/bayestar/bayestar_plot_found_injections
src/inspiral/bayestar/bayestar_prune_neighborhood_tmpltbank
src/inspiral/bayestar/bayestar_realize_coincs
src/inspiral/bayestar/bayestar_sim_to_tmpltbank
src/inspiral/bayestar/python_config.sed
src/inspiral/posterior/ChiSquare_test
src/inspiral/posterior/SPINspiral/lalapps_spinspiral
src/inspiral/posterior/lalapps_coherence_test
......
......@@ -34,6 +34,7 @@ AC_CONFIG_FILES([\
src/pulsar/TimingTests/Makefile \
src/pulsar/TDS_isolated/Makefile \
src/inspiral/Makefile \
src/inspiral/bayestar/Makefile \
src/inspiral/posterior/Makefile \
src/inspiral/posterior/SPINspiral/Makefile \
src/inspiral/posterior/mpi/Makefile \
......
## Process this file with automake to produce Makefile.in
SUBDIRS = posterior
SUBDIRS = bayestar posterior
INSPIRALSRC = inspiral.c inspiralutils.c inspiral.h
TMPLTBANKSRC = tmpltbank.c inspiralutils.c inspiral.h
......
## Process this file with automake to produce Makefile.in
if HAVE_PYTHON
bin_SCRIPTS = \
bayestar_aggregate_found_injections \
bayestar_lattice_tmpltbank \
bayestar_localize_coincs \
bayestar_localize_gracedb \
bayestar_localize_lvalert \
bayestar_plot_found_injections \
bayestar_prune_neighborhood_tmpltbank \
bayestar_realize_coincs \
bayestar_sim_to_tmpltbank \
bayestar_littlehope
endif
python_config.sed: Makefile
@-rm -f python_config.sed
@echo "s+@PYTHONLIBDIR@+@pkgpythondir@+g" > python_config.sed
@echo "s+@PYTHONPROG@+@PYTHON@+g" >> python_config.sed
bayestar_aggregate_found_injections: bayestar_aggregate_found_injections.in python_config.sed
$(AM_V_GEN)sed -f python_config.sed $< > $@
bayestar_lattice_tmpltbank: bayestar_lattice_tmpltbank.in python_config.sed
$(AM_V_GEN)sed -f python_config.sed $< > $@
bayestar_localize_coincs: bayestar_localize_coincs.in python_config.sed
$(AM_V_GEN)sed -f python_config.sed $< > $@
bayestar_localize_gracedb: bayestar_localize_gracedb.in python_config.sed
$(AM_V_GEN)sed -f python_config.sed $< > $@
bayestar_localize_lvalert: bayestar_localize_lvalert.in python_config.sed
$(AM_V_GEN)sed -f python_config.sed $< > $@
bayestar_plot_found_injections: bayestar_plot_found_injections.in python_config.sed
$(AM_V_GEN)sed -f python_config.sed $< > $@
bayestar_prune_neighborhood_tmpltbank: bayestar_prune_neighborhood_tmpltbank.in python_config.sed
$(AM_V_GEN)sed -f python_config.sed $< > $@
bayestar_realize_coincs: bayestar_realize_coincs.in python_config.sed
$(AM_V_GEN)sed -f python_config.sed $< > $@
bayestar_sim_to_tmpltbank: bayestar_sim_to_tmpltbank.in python_config.sed
$(AM_V_GEN)sed -f python_config.sed $< > $@
bayestar_littlehope: bayestar_littlehope.in python_config.sed
$(AM_V_GEN)sed -f python_config.sed $< > $@
EXTRA_DIST = \
bayestar_aggregate_found_injections.in \
bayestar_lattice_tmpltbank.in \
bayestar_localize_coincs.in \
bayestar_localize_gracedb.in \
bayestar_localize_lvalert.in \
bayestar_plot_found_injections.in \
bayestar_prune_neighborhood_tmpltbank.in \
bayestar_realize_coincs.in \
bayestar_sim_to_tmpltbank.in \
bayestar_littlehope.in
CLEANFILES = \
$(bin_SCRIPTS) \
python_config.sed
#!@PYTHONPROG@
#
# Copyright (C) 2013 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
"""
Match sky maps with injections in an inspinjfind-style sqlite database and print
summary values for each sky map:
* event ID
* false alarm rate
* searched area
* searched posterior probability
* angle between true sky location and maximum a posteriori estimate
* runtime in seconds
"""
__author__ = "Leo Singer <leo.singer@ligo.org>"
# Command line interface.
from optparse import Option, OptionParser
from lalinference.bayestar import command
parser = OptionParser(
formatter=command.NewlinePreservingHelpFormatter(),
description=__doc__,
usage="%prog DATABASE.sqlite FILE1.fits[.gz] FILE2.fits[.gz] ...")
opts, args = parser.parse_args()
try:
dbfilename = args[0]
fitsfilenames = args[1:]
except IndexError:
parser.error("not enough command line arguments")
# Imports.
import os
import numpy as np
import healpy as hp
import sqlite3
from pylal.progress import ProgressBar
from lalinference.bayestar import fits
sql = """
SELECT DISTINCT sim.longitude AS ra, sim.latitude AS dec, ci.combined_far AS far
FROM coinc_event_map AS cem1 INNER JOIN coinc_event_map AS cem2
ON (cem1.coinc_event_id = cem2.coinc_event_id)
INNER JOIN sim_inspiral AS sim ON (cem1.event_id = sim.simulation_id)
INNER JOIN coinc_inspiral AS ci ON (cem2.event_id = ci.coinc_event_id)
WHERE cem1.table_name = 'sim_inspiral' AND cem2.table_name = 'coinc_event'
AND cem2.event_id = ?"""
def angle_distance(theta0, phi0, theta1, phi1):
cos_angle_distance = np.cos(phi1 - phi0) * np.sin(theta0) * np.sin(theta1) + np.cos(theta0) * np.cos(theta1)
if cos_angle_distance > 1:
return 0.
elif cos_angle_distance < -1:
return np.pi
else:
return np.arccos(cos_angle_distance)
def find_injection(sky_map, true_ra, true_dec):
"""
Given a sky map and the true right ascension and declination (in radians),
find the smallest area in deg^2 that would have to be searched to find the
source, the smallest posterior mass, and the angular offset in degrees from
the true location to the maximum (mode) of the posterior.
"""
# Compute the HEALPix lateral resolution parameter for this sky map.
npix = len(sky_map)
nside = hp.npix2nside(npix)
# Convert from ra, dec to conventional spherical polar coordinates.
true_theta = 0.5 * np.pi - true_dec
true_phi = true_ra
# Find the HEALPix pixel index of the mode of the posterior and of the
# true sky location.
mode_pix = np.argmax(sky_map)
true_pix = hp.ang2pix(nside, true_theta, true_phi)
# Compute spherical polar coordinates of true location.
mode_theta, mode_phi = hp.pix2ang(nside, mode_pix)
# Sort the pixels in the sky map by descending posterior probability and
# form the cumulative sum. Record the total value.
indices = np.argsort(sky_map)[::-1]
cum_sky_map = np.cumsum(sky_map[indices])
# Find the index of the true location in the cumulative distribution.
idx = (i for i, pix in enumerate(indices) if pix == true_pix).next()
# Find the smallest area that would have to be searched to find the true
# location.
searched_area = (idx + 1) * hp.nside2pixarea(nside, degrees=True)
# Find the smallest posterior mass that would have to be searched to find
# the true location.
searched_prob = cum_sky_map[idx]
# Permute the cumulative distribution so that it is indexed the same way
# as the original sky map.
cum_sky_map[indices] = cum_sky_map
# Find the angular offset between the mode and true locations.
offset = np.rad2deg(angle_distance(true_theta, true_phi, mode_theta, mode_phi))
# Done.
return searched_area, searched_prob, offset
progress = ProgressBar()
progress.update(-1, 'opening database')
db = sqlite3.connect(dbfilename)
print 'objid,far,searched_area,searched_prob,offset,runtime'
for fitsfilename in progress.iterate(fitsfilenames):
sky_map, metadata = fits.read_sky_map(fitsfilename)
coinc_event_id = metadata['objid']
try:
runtime = metadata['runtime']
except KeyError:
runtime = None
true_ra, true_dec, far = db.execute(sql, (coinc_event_id,)).fetchone()
searched_area, searched_prob, offset = find_injection(sky_map, true_ra, true_dec)
print ','.join(str(item) for item in (coinc_event_id, far, searched_area, searched_prob, offset, runtime))
#!@PYTHONPROG@
#
# Copyright (C) 2013 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
__doc__ = """
Create a template bank that samples a regular lattice in (tau0, tau3), starting
from an initial (mass1, mass2), with lattice points spaced according to the
metric at the initial point.
"""
__author__ = "Leo Singer <leo.singer@ligo.org>"
# Command line interface.
from optparse import Option, OptionParser
from lalinference.bayestar import command
parser = OptionParser(
formatter = command.NewlinePreservingHelpFormatter(),
description = __doc__,
usage="%prog [options] [INPUT.xml[.gz]] [-o OUTPUT.xml[.gz]]",
option_list = [
Option("--initial-mass1", metavar="Msun", type=float, default=1.4,
help="Mass of first component of an initial lattice point (default=1.4)"),
Option("--initial-mass2", metavar="Msun", type=float, default=1.4,
help="Mass of second component of an initial lattice point (default=1.4)"),
Option("--min-mass", metavar="Msun", type=float, default=1.,
help="Minimum component mass (default=1)"),
Option("--max-mass", metavar="Msun", type=float, default=3.,
help="Maximum component mass (default=3)"),
Option("-o", "--output", metavar="OUTPUT.xml[.gz]", default="/dev/stdout",
help="Name of output file (default=stdout)")
]
)
opts, args = parser.parse_args()
infilename = command.get_input_filename(parser, args)
# Python standard library imports.
import os
# LIGO-LW XML imports.
from glue.ligolw import ligolw
from glue.ligolw.utils import process as ligolw_process
from glue.ligolw import table as ligolw_table
from glue.ligolw import utils as ligolw_utils
from glue.ligolw import lsctables
# glue and LAL imports.
import lal
import lalinspiral.sbank.tau0tau3
import lalsimulation
# BAYESTAR imports.
from lalinference.bayestar import timing
# Other imports.
import numpy as np
import scipy.linalg
# Open output file.
xmldoc = ligolw.Document()
xmldoc.appendChild(ligolw.LIGO_LW())
# Write process metadata to output file.
process = ligolw_process.register_to_xmldoc(xmldoc, parser.get_prog_name(),
opts.__dict__)
# Create a SnglInspiral table and initialize its row ID counter.
sngl_inspiral_table = lsctables.New(lsctables.SnglInspiralTable)
xmldoc.childNodes[0].appendChild(sngl_inspiral_table)
sngl_inspiral_table.set_next_id(lsctables.SnglInspiralID(0))
f_low = 10.
f_high = 2048.
df = 0.1
initial_mchirp = lalinspiral.sbank.tau0tau3.m1m2_to_mchirp(opts.initial_mass1, opts.initial_mass2)
initial_eta = opts.initial_mass1 * opts.initial_mass2 / (opts.initial_mass1 + opts.initial_mass2)**2
initial_chi = 0.
initial_chirp_times = lalsimulation.SimInspiralTaylorF2RedSpinChirpTimesFromMchirpEtaChi(initial_mchirp, initial_eta, initial_chi, f_low)
initial_theta0_theta3 = initial_chirp_times[:2]
# Sampled PSD.
S = lal.CreateREAL8Vector(int(f_high // df))
S.data = [lalsimulation.SimNoisePSDaLIGOZeroDetHighPower(i * df)
for i in range(len(S.data))]
# Allocate noise moments.
moments = [lal.CreateREAL8Vector(int((f_high - f_low) // df)) for _ in range(29)]
# Compute noise moments.
lalsimulation.SimInspiralTaylorF2RedSpinComputeNoiseMoments(
*(moments + [S, f_low, df]))
# Compute Fisher information matrix. Note factor of 2:
# The function SimInspiralTaylorF2RedSpinFisherMatrixChirpTimes
# returns the Fisher matrix for an SNR of 1/sqrt(2). The factor
# of 2 makes this the Fisher matrix at unit SNR.
I = lalsimulation.SimInspiralTaylorF2RedSpinFisherMatrixChirpTimes(
*(initial_chirp_times + [f_low, df] + moments)).data * 2
# Blockwise separation of Fisher matrix. Parameters are in the following order:
# theta0, theta3, theta3S, t0, phi0
IA = I[0:2, 0:2] # (theta0, theta3) block
IB = I[0:2, 3:5] # cross block
ID = I[3:5, 3:5] # (time, phase) block
# Metric. We are dropping the theta3S terms completely and projecting out
# time and phase.
metric = IA - np.dot(IB, scipy.linalg.solve(ID, IB.T, sym_pos=True))
# Eigendecomposition of metric
metric_eigenvalues, metric_eigenvectors = np.linalg.eigh(metric)
# Shift between adjacent lattice points.
# FIXME: square root or no?
delta_theta0_theta3 = np.dot(metric_eigenvectors, np.diag(1 / np.sqrt(metric_eigenvalues)))
# FIXME: Determine appropriate boundaries to looping over lots of points that
# we are going to skip.
#
# T. Cokelaer (2007, http://dx.doi.org/10.1103/PhysRevD.76.102004) describes
# relationships between the component mass limits and the (tau0, tau3)
# boundaries.
n = 800
i0, i1 = np.mgrid[-n:n+1, -n:n+1]
i = np.column_stack((i0.flatten(), i1.flatten()))
# FIXME: Come up with a more natural way to specify the template spacing.
skip = 10
theta0_theta3 = np.dot(i, skip * delta_theta0_theta3.T) + initial_theta0_theta3
for th0, th3 in theta0_theta3:
th3S = 0
tau0 = th0 / (2 * np.pi * f_low)
tau3 = -th3 / (2 * np.pi * f_low)
mchirp, eta, chi = lalsimulation.SimInspiralTaylorF2RedSpinMchirpEtaChiFromChirpTimes(
th0, th3, th3S, f_low
)
# Skip if either mchirp, eta, or chi are unphysical, unless this is the
# initial point, which may be slightly unphysical just due to roundoff
if np.all([th0, th3] == initial_theta0_theta3):
mchirp = initial_mchirp
eta = initial_eta
mass1 = opts.initial_mass1
mass2 = opts.initial_mass2
mtotal = mass1 + mass2
elif not (mchirp >= 0 and 0 <= eta <= 0.25 and -1 <= chi <= 1):
continue
else:
mtotal = mchirp * eta ** -0.6
mass1, mass2 = sorted(np.roots([1, -mtotal, eta * mtotal**2]))
# Skip this one unless both component masses are in the appropriate range.
if not(opts.min_mass <= mass1 <= opts.max_mass and opts.min_mass <= mass2 <= opts.max_mass):
continue
# Create new sngl_inspiral row and initialize its columns to None,
# which produces an empty field in the XML output.
sngl_inspiral = lsctables.SnglInspiral()
for validcolumn in sngl_inspiral_table.validcolumns.iterkeys():
setattr(sngl_inspiral, validcolumn, None)
# Populate the row's fields.
sngl_inspiral.event_id = sngl_inspiral_table.get_next_id()
sngl_inspiral.mass1 = mass1
sngl_inspiral.mass2 = mass2
sngl_inspiral.tau0 = tau0
sngl_inspiral.tau3 = tau3
sngl_inspiral.mtotal = mtotal
sngl_inspiral.mchirp = mchirp
sngl_inspiral.eta = eta
sngl_inspiral.chi = chi
sngl_inspiral.f_final = timing.get_f_lso(mass1, mass2)
# Add the row to the table in the document.
sngl_inspiral_table.append(sngl_inspiral)
# Record process end time.
ligolw_process.set_process_end_time(process)
# Write output file.
ligolw_utils.write_filename(xmldoc, opts.output,
gz=(os.path.splitext(opts.output)[-1]==".gz"))
#!@PYTHONPROG@
#
# Copyright (C) 2013 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
"""
Synthesize triggers for simulated sources using a miniature matched-filter
detection pipeline. The input file (or stdin if the input file is omitted)
should be an optionally gzip-compressed LIGO-LW XML file of the form
produced by lalapps_inspinj. The output file (or stdout if omitted) will be an
optionally gzip-compressed LIGO-LW XML file containing single-detector triggers
and coincidences. A single template that has the same intrinsic parameters as
the injection is used.
"""
__author__ = "Leo Singer <leo.singer@ligo.org>"
# Determine list of known detectors for command line arguments.
import lal
available_ifos = sorted(det.frDetector.prefix
for det in lal.lalCachedDetectors)
# List of interpolation methods
available_interp_methods = [
"catmull-rom", "lanczos", "nearest-neighbor", "quadratic-fit"]
# Command line interface.
from optparse import Option, OptionParser
from lalinference.bayestar import command
parser = OptionParser(
formatter = command.NewlinePreservingHelpFormatter(),
description = __doc__,
usage="%prog [options] --template-bank TMPLTBANK.xml[.gz] [INPUT.xml[.gz]] [-o OUTPUT.xml[.gz]]",
option_list = [
Option("-o", "--output", metavar="OUTPUT.xml[.gz]", default="/dev/stdout",
help="Name of output file (default=stdout)"),
Option("--detector", metavar='|'.join(available_ifos), action="append",
help="Detectors to use. May be specified multiple times.",
choices=available_ifos),
Option("--trigger-window", type=float, default=0.1, metavar="SECONDS",
help="Search for a trigger across this many seconds before and after the time of the injection (default=0.1)"),
Option("--interp-method", metavar='|'.join(available_interp_methods),
default="lanczos", choices=available_interp_methods,
help="Trigger interpolation method (default=lanczos)"),
Option("--interp-window", metavar="SAMPLES", type=int, default=2,
help="Trigger interpolation window (default=2)"),
Option("--waveform",
help="Waveform to use for injections"),
Option("--snr-threshold", type=float, default=8.,
help="Single-detector SNR threshold (default=8)"),
Option("--min-triggers", type=int, default=3,
help="Emit coincidences only when at least this many triggers are found (default=3)"),
Option("-R", "--repeat-first-injection", type=int, default=None,
help="Instead of performing each injection once, just perform the first injection this many times."),
Option("--template-bank", metavar="TMPLTBANK.xml[.gz]",
help="Name of template bank file (required)")
]
)
opts, args = parser.parse_args()
infilename = command.get_input_filename(parser, args)
command.check_required_arguments(parser, opts, 'waveform', 'template_bank')
# Python standard library imports.
import os
# LIGO-LW XML imports.
from glue.ligolw import ligolw
from glue.ligolw.utils import process as ligolw_process
from glue.ligolw.utils import search_summary as ligolw_search_summary
from glue.ligolw import table as ligolw_table
from pylal import ligolw_thinca
from glue.ligolw import utils as ligolw_utils
from glue.ligolw import lsctables
# glue and LAL imports.
from glue import segments
import glue.lal
import lal, lalsimulation
# BAYESTAR imports.
from lalinference.bayestar import timing
from lalinference.bayestar import filter
from lalinference.bayestar import ligolw as ligolw_bayestar
# Other imports.
import numpy as np
template_approximant, template_amplitude_order, template_phase_order = timing.get_approximant_and_orders_from_string(opts.waveform)
# FIXME: sample rate could be a command line option; template duration and data
# duration should be determined from chirp time
sample_rate = 4096 # sample rate in Hz
template_duration = 128 # template duration in seconds
template_length = sample_rate * template_duration # template length in samples
data_duration = 512 # data duration in seconds
data_length = sample_rate * data_duration # data length in samples
# Open output file.
out_xmldoc = ligolw.Document()
out_xmldoc.appendChild(ligolw.LIGO_LW())
# Write process metadata to output file.
process = ligolw_process.register_to_xmldoc(out_xmldoc, parser.get_prog_name(),
opts.__dict__, ifos=opts.detector, comment="Little hope!")
# Add search summary to output file.
all_time = segments.segment([glue.lal.LIGOTimeGPS(0), glue.lal.LIGOTimeGPS(2e9)])
search_summary_table = lsctables.New(lsctables.SearchSummaryTable)
out_xmldoc.childNodes[0].appendChild(search_summary_table)
summary = ligolw_search_summary.append_search_summary(out_xmldoc, process,
inseg=all_time, outseg=all_time)
# Read template bank file.
xmldoc = ligolw_utils.load_filename(opts.template_bank)
# Determine the low frequency cutoff from the template bank file.
template_bank_f_low = ligolw_bayestar.get_temlate_bank_f_low(xmldoc)
template_bank = ligolw_table.get_table(xmldoc,
lsctables.SnglInspiralTable.tableName)
# Read injection file.
xmldoc = ligolw_utils.load_filename(infilename)
# Extract simulation table from injection file.
sim_inspiral_table = ligolw_table.get_table(xmldoc,
lsctables.SimInspiralTable.tableName)
# Create a SnglInspiral table and initialize its row ID counter.
sngl_inspiral_table = lsctables.New(lsctables.SnglInspiralTable)
out_xmldoc.childNodes[0].appendChild(sngl_inspiral_table)
sngl_inspiral_table.set_next_id(lsctables.SnglInspiralID(0))
# Create a time slide entry. Needed for coinc_event rows.
try: