Commit 244443de authored by Leo Pound Singer's avatar Leo Pound Singer
Browse files

Add bayestar-mcmc tool

parent 9f32596d
......@@ -5,7 +5,9 @@ Changelog
0.0.13 (unreleased)
===================
- No changes yet.
- Add ``bayestar-mcmc`` tool for pure Markov Chain Monte Carlo parameter
estimation, without sky map postprocessing but with options for holding
parameters at fixed values.
0.0.12 (2018-07-18)
===================
......
......@@ -70,6 +70,7 @@ BAYESTAR Rapid Sky Localization
ligo/skymap/tool/bayestar_localize_coincs
ligo/skymap/tool/bayestar_localize_lvalert
ligo/skymap/tool/bayestar_mcmc
ligo/skymap/tool/bayestar_realize_coincs
ligo/skymap/tool/bayestar_sample_model_psd
......
Markov-Chain Monte Carlo Localization (`bayestar-mcmc`)
=======================================================
.. argparse::
:module: ligo.skymap.tool.bayestar_mcmc
:func: parser
......@@ -49,7 +49,7 @@ from .. import moc
from .. import healpix_tree
from .. import version
from .. import core
from ..core import log_posterior_toa_phoa_snr
from ..core import log_posterior_toa_phoa_snr as _log_posterior_toa_phoa_snr
from ..util.numpy import require_contiguous
from .ez_emcee import ez_emcee
......@@ -57,7 +57,18 @@ __all__ = ('derasterize', 'localize', 'rasterize')
log = logging.getLogger('BAYESTAR')
log_posterior_toa_phoa_snr = require_contiguous(log_posterior_toa_phoa_snr)
_log_posterior_toa_phoa_snr = require_contiguous(_log_posterior_toa_phoa_snr)
# Wrap so that ufunc parameter names are known
def log_posterior_toa_phoa_snr(
ra, sin_dec, distance, u, twopsi, t, min_distance, max_distance,
prior_distance_power, cosmology, gmst, sample_rate, epochs, snrs,
responses, locations, horizons):
return _log_posterior_toa_phoa_snr(
ra, sin_dec, distance, u, twopsi, t, min_distance, max_distance,
prior_distance_power, cosmology, gmst, sample_rate, epochs, snrs,
responses, locations, horizons)
def log_post(params, *args, **kwargs):
......
#
# Copyright (C) 2013-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.
#
"""
Markov-Chain Monte Carlo sky localization.
"""
from argparse import FileType
from . import (
ArgumentParser, waveform_parser, prior_parser, mkpath)
def parser():
parser = ArgumentParser(parents=[waveform_parser, prior_parser])
parser.add_argument(
'input', metavar='INPUT.{hdf,xml,xml.gz,sqlite}', default='-',
nargs='+', type=FileType('rb'),
help='Input LIGO-LW XML file, SQLite file, or PyCBC HDF5 files. '
'For PyCBC, you must supply the coincidence file '
'(e.g. "H1L1-HDFINJFIND.hdf" or "H1L1-STATMAP.hdf"), '
'the template bank file (e.g. H1L1-BANK2HDF.hdf), '
'the single-detector merged PSD files '
'(e.g. "H1-MERGE_PSDS.hdf" and "L1-MERGE_PSDS.hdf"), '
'and the single-detector merged trigger files '
'(e.g. "H1-HDF_TRIGGER_MERGE.hdf" and '
'"L1-HDF_TRIGGER_MERGE.hdf"), '
'in any order.')
parser.add_argument(
'--pycbc-sample', default='foreground',
help='(PyCBC only) sample population')
parser.add_argument(
'--coinc-event-id', type=int, nargs='*',
help='run on only these specified events')
parser.add_argument(
'--output', '-o', default='.',
help='output directory')
parser.add_argument(
'--condor-submit', action='store_true',
help='submit to Condor instead of running locally')
group = parser.add_argument_group(
'fixed parameter options',
'Options to hold certain parameters constant')
group.add_argument('--ra', metavar='DEG', help='Right ascension')
group.add_argument('--dec', metavar='DEG', help='Declination')
group.add_argument('--distance', metavar='Mpc', help='Luminosity distance')
return parser
def identity(x):
return x
def main(args=None):
opts = parser().parse_args(args)
import logging
log = logging.getLogger('BAYESTAR')
# BAYESTAR imports.
from ..io import events, hdf5
from ..bayestar import condition, condition_prior, ez_emcee, log_post
# Other imports.
from astropy.table import Table
import numpy as np
import os
from collections import OrderedDict
import sys
# Squelch annoying and uniformative LAL log messages.
import lal
lal.ClobberDebugLevel(lal.LALNDEBUG)
# Read coinc file.
log.info(
'%s:reading input files', ','.join(file.name for file in opts.input))
event_source = events.open(*opts.input, sample=opts.pycbc_sample)
mkpath(opts.output)
if opts.condor_submit:
if opts.coinc_event_id:
raise ValueError(
'must not set --coinc-event-id with --condor-submit')
cmd = ['condor_submit', 'accounting_group=ligo.dev.o3.cbc.pe.bayestar',
'on_exit_remove = (ExitBySignal == False) && (ExitCode == 0)',
'on_exit_hold = (ExitBySignal == True) || (ExitCode != 0)',
'on_exit_hold_reason = (ExitBySignal == True ? strcat('
'"The job exited with signal ", ExitSignal) : strcat('
'"The job exited with signal ", ExitCode))',
'request_memory = 1000 MB',
'universe=vanilla', 'getenv=true',
'executable=/usr/bin/env',
'JobBatchName=BAYESTAR', 'environment="OMP_NUM_THREADS=1"',
'error=' + os.path.join(opts.output, '$(CoincEventId).err'),
'log=' + os.path.join(opts.output, '$(CoincEventId).log'),
'arguments="' + ' '.join(arg for arg in sys.argv
if arg != '--condor-submit') +
' --coinc-event-id $(CoincEventId)"',
'-append', 'queue CoincEventId in ' + ' '.join(
str(coinc_event_id) for coinc_event_id in event_source),
'/dev/null']
os.execvp('condor_submit', cmd)
# Loop over all coinc_event <-> sim_inspiral coincs.
if opts.coinc_event_id:
event_source = OrderedDict(
(key, event_source[key]) for key in opts.coinc_event_id)
for int_coinc_event_id, event in event_source.items():
coinc_event_id = 'coinc_event:coinc_event_id:{}'.format(
int_coinc_event_id)
log.info('%s:preparing', coinc_event_id)
epoch, sample_rate, toas, snr_series, responses, locations, horizons \
= condition(event, waveform=opts.waveform, f_low=opts.f_low,
enable_snr_series=opts.enable_snr_series,
f_high_truncate=opts.f_high_truncate)
min_distance, max_distance, prior_distance_power, cosmology = \
condition_prior(horizons, opts.min_distance, opts.max_distance,
opts.prior_distance_power, opts.cosmology)
gmst = lal.GreenwichMeanSiderealTime(epoch)
args = (min_distance, max_distance, prior_distance_power, cosmology,
gmst, sample_rate, toas, snr_series, responses, locations,
horizons)
max_abs_t = 2 * snr_series.data.shape[1] / sample_rate
xmin = [0, -1, min_distance, -1, 0, 0]
xmax = [2 * np.pi, 1, max_distance, 1, 2 * np.pi, 2 * max_abs_t]
names = 'ra dec distance inclination twopsi time'.split()
transformed_names = 'ra sin_dec distance u twopsi time'.split()
forward_transforms = [identity, np.sin, identity,
np.cos, identity, identity]
reverse_transforms = [identity, np.arcsin, identity,
np.arccos, identity, identity]
kwargs = {}
# Fix parameters
for i, key in reversed(list(enumerate(['ra', 'dec', 'distance']))):
try:
value = getattr(opts, key)
except AttributeError:
pass
else:
if key in ['ra', 'dec']:
# FIXME: figure out a more elegant way to address different
# units in command line arguments and posterior samples
value = np.deg2rad(value)
kwargs[transformed_names[i]] = value
del (xmin[i], xmax[i], names[i], transformed_names[i],
forward_transforms[i], reverse_transforms[i])
log.info('%s:sampling', coinc_event_id)
# Run MCMC
chain = ez_emcee(log_post, xmin, xmax,
args=args, kwargs=kwargs, vectorize=True)
# Transform back from sin_dec to dec and cos_inclination to inclination
for i, func in enumerate(reverse_transforms):
chain[:, i] = func(chain[:, i])
# Create Astropy table
chain = Table(rows=chain, names=names)
log.info('%s:saving posterior samples', coinc_event_id)
hdf5.write_samples(
chain,
os.path.join(opts.output, '{}.hdf5'.format(int_coinc_event_id)),
path='/bayestar/posterior_samples')
......@@ -89,6 +89,7 @@ namespace_packages = ligo
[entry_points]
bayestar-localize-coincs = ligo.skymap.tool.bayestar_localize_coincs:main
bayestar-localize-lvalert = ligo.skymap.tool.bayestar_localize_lvalert:main
bayestar-mcmc = ligo.skymap.tool.bayestar_mcmc:main
bayestar-realize-coincs = ligo.skymap.tool.bayestar_realize_coincs:main
bayestar-sample-model-psd = ligo.skymap.tool.bayestar_sample_model_psd:main
ligo-skymap-combine = ligo.skymap.tool.ligo_skymap_combine:main
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment