diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 65a54eb66cef71ea4552596cfab13ed3040a5690..f6ceb56370c370db095756dbc954cbdde6b791f4 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -18,6 +18,8 @@ python-2:
stage: test
image: bilbydev/test-suite-py2
before_script:
+ # Source the .bashrc for MultiNest
+ - source /root/.bashrc
# Install the dependencies specified in the Pipfile
- pipenv install --two --python=/opt/conda/bin/python2 --system --deploy
script:
@@ -30,6 +32,8 @@ python-3:
stage: test
image: bilbydev/test-suite-py3
before_script:
+ # Source the .bashrc for MultiNest
+ - source /root/.bashrc
# Install the dependencies specified in the Pipfile
- pipenv install --three --python=/opt/conda/bin/python --system --deploy
script:
diff --git a/bilby/core/sampler/base_sampler.py b/bilby/core/sampler/base_sampler.py
index fdbc5a32553d2fc3c211a82df07066456dcf384b..a8ddc2fb0cf68effcfdcc5ea8bfc57ccdfe66e65 100644
--- a/bilby/core/sampler/base_sampler.py
+++ b/bilby/core/sampler/base_sampler.py
@@ -1,3 +1,4 @@
+from __future__ import absolute_import
import datetime
import numpy as np
from pandas import DataFrame
@@ -439,7 +440,7 @@ class MCMCSampler(Sampler):
def print_nburn_logging_info(self):
""" Prints logging info as to how nburn was calculated """
- if type(self.kwargs['nburn']) in [float, int]:
+ if type(self.nburn) in [float, int]:
logger.info("Discarding {} steps for burn-in".format(self.nburn))
elif self.result.max_autocorrelation_time is None:
logger.info("Autocorrelation time not calculated, discarding {} "
diff --git a/bilby/core/sampler/emcee.py b/bilby/core/sampler/emcee.py
index a641586e250fa7d320a023373006217c0ddb42fb..190e8eec7df6ee915748d3c372de701551142102 100644
--- a/bilby/core/sampler/emcee.py
+++ b/bilby/core/sampler/emcee.py
@@ -1,4 +1,4 @@
-from __future__ import absolute_import
+from __future__ import absolute_import, print_function
import numpy as np
from pandas import DataFrame
from ..utils import logger, get_progress_bar
diff --git a/bilby/core/sampler/ptemcee.py b/bilby/core/sampler/ptemcee.py
index 2efcab08c8ee7ae2808bfec4c6f7b798319358de..5ac414c0307d91f2e00963cfc0143f9e7c7639c6 100644
--- a/bilby/core/sampler/ptemcee.py
+++ b/bilby/core/sampler/ptemcee.py
@@ -1,4 +1,4 @@
-from __future__ import absolute_import
+from __future__ import absolute_import, division, print_function
import numpy as np
from ..utils import get_progress_bar, logger
@@ -69,16 +69,18 @@ class Ptemcee(Emcee):
total=self.nsteps):
pass
- self.result.nburn = self.nburn
+ self.calculate_autocorrelation(sampler.chain.reshape((-1, self.ndim)))
self.result.sampler_output = np.nan
+ self.print_nburn_logging_info()
+ self.result.nburn = self.nburn
+ if self.result.nburn > self.nsteps:
+ logger.warning('Chain not burned in, no samples generated.')
self.result.samples = sampler.chain[0, :, self.nburn:, :].reshape(
(-1, self.ndim))
+ self.result.betas = sampler.betas
+ self.result.log_evidence, self.result.log_evidence_err =\
+ sampler.log_evidence_estimate(
+ sampler.loglikelihood, self.nburn / self.nsteps)
self.result.walkers = sampler.chain[0, :, :, :]
- self.result.log_evidence = np.nan
- self.result.log_evidence_err = np.nan
- logger.info("Max autocorr time = {}"
- .format(np.max(sampler.get_autocorr_time())))
- logger.info("Tswap frac = {}"
- .format(sampler.tswap_acceptance_fraction))
return self.result
diff --git a/bilby/core/sampler/pymc3.py b/bilby/core/sampler/pymc3.py
index ae713532529d55edeb2380b2781d98dd19b4f752..a797836f26617118acffa67915388d9cbd90fe29 100644
--- a/bilby/core/sampler/pymc3.py
+++ b/bilby/core/sampler/pymc3.py
@@ -1,4 +1,4 @@
-from __future__ import absolute_import
+from __future__ import absolute_import, print_function
from collections import OrderedDict
import inspect
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index 6b8752fd3ef1b6cbf169760944b50dcd156cf991..4d2457ab24159df52b2697ea86eca5275b401ff0 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -405,3 +405,81 @@ def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=1
else:
logger.warning('No data loaded.')
return None
+
+
+def save_to_fits(posterior, outdir, label):
+ """ Generate a fits file from a posterior array """
+ from astropy.io import fits
+ from astropy.units import pixel
+ from astropy.table import Table
+ import healpy as hp
+ nside = hp.get_nside(posterior)
+ npix = hp.nside2npix(nside)
+ logger.debug('Generating table')
+ m = Table([posterior], names=['PROB'])
+ m['PROB'].unit = pixel ** -1
+
+ ordering = 'RING'
+ extra_header = [('PIXTYPE', 'HEALPIX',
+ 'HEALPIX pixelisation'),
+ ('ORDERING', ordering,
+ 'Pixel ordering scheme: RING, NESTED, or NUNIQ'),
+ ('COORDSYS', 'C',
+ 'Ecliptic, Galactic or Celestial (equatorial)'),
+ ('NSIDE', hp.npix2nside(npix),
+ 'Resolution parameter of HEALPIX'),
+ ('INDXSCHM', 'IMPLICIT',
+ 'Indexing: IMPLICIT or EXPLICIT')]
+
+ fname = '{}/{}_{}.fits'.format(outdir, label, nside)
+ hdu = fits.table_to_hdu(m)
+ hdu.header.extend(extra_header)
+ hdulist = fits.HDUList([fits.PrimaryHDU(), hdu])
+ logger.debug('Writing to a fits file')
+ hdulist.writeto(fname, overwrite=True)
+
+
+def plot_skymap(result, center='120d -40d', nside=512):
+ """ Generate a sky map from a result """
+ import scipy
+ from astropy.units import deg
+ import healpy as hp
+ import ligo.skymap.plot # noqa
+ import matplotlib.pyplot as plt
+ logger.debug('Generating skymap')
+
+ logger.debug('Reading in ra and dec, creating kde and converting')
+ ra_dec_radians = result.posterior[['ra', 'dec']].values
+ kde = scipy.stats.gaussian_kde(ra_dec_radians.T)
+ npix = hp.nside2npix(nside)
+ ipix = range(npix)
+ theta, phi = hp.pix2ang(nside, ipix)
+ ra = phi
+ dec = 0.5 * np.pi - theta
+
+ logger.debug('Generating posterior')
+ post = kde(np.row_stack([ra, dec]))
+ post /= np.sum(post * hp.nside2pixarea(nside))
+
+ fig = plt.figure(figsize=(5, 5))
+ ax = plt.axes([0.05, 0.05, 0.9, 0.9],
+ projection='astro globe',
+ center=center)
+ ax.coords.grid(True, linestyle='--')
+ lon = ax.coords[0]
+ lat = ax.coords[1]
+ lon.set_ticks(exclude_overlapping=True, spacing=45 * deg)
+ lat.set_ticks(spacing=30 * deg)
+
+ lon.set_major_formatter('dd')
+ lat.set_major_formatter('hh')
+ lon.set_ticklabel(color='k')
+ lat.set_ticklabel(color='k')
+
+ logger.debug('Plotting sky map')
+ ax.imshow_hpx(post)
+
+ lon.set_ticks_visible(False)
+ lat.set_ticks_visible(False)
+
+ fig.savefig('{}/{}_skymap.png'.format(result.outdir, result.label))
diff --git a/examples/injection_examples/plot_skymap.py b/examples/injection_examples/plot_skymap.py
new file mode 100644
index 0000000000000000000000000000000000000000..5f1db3fc1c293dede1c59bc20445b554be135cdc
--- /dev/null
+++ b/examples/injection_examples/plot_skymap.py
@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+"""
+Example script which produces posterior samples of ra and dec and generates a
+skymap
+"""
+from __future__ import division, print_function
+import bilby
+
+duration = 4.
+sampling_frequency = 2048.
+outdir = 'outdir'
+label = 'plot_skymap'
+injection_parameters = dict(
+ mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
+ phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659,
+ phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-0.2108)
+
+waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
+ reference_frequency=50.)
+
+waveform_generator = bilby.gw.WaveformGenerator(
+ duration=duration, sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+ parameters=injection_parameters, waveform_arguments=waveform_arguments)
+
+ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+ifos.set_strain_data_from_power_spectral_densities(
+ sampling_frequency=sampling_frequency, duration=duration,
+ start_time=injection_parameters['geocent_time'] - 3)
+ifos.inject_signal(waveform_generator=waveform_generator,
+ parameters=injection_parameters)
+
+priors = bilby.gw.prior.BBHPriorSet()
+for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi',
+ 'mass_1', 'mass_2', 'phase', 'geocent_time', 'luminosity_distance',
+ 'iota']:
+ priors[key] = injection_parameters[key]
+
+likelihood = bilby.gw.GravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=waveform_generator,
+ time_marginalization=True, phase_marginalization=True,
+ distance_marginalization=False, prior=priors)
+
+result = bilby.run_sampler(
+ likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
+ injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+# make some plots of the outputs
+result.plot_corner()
+
+bilby.gw.utils.plot_skymap(result, nside=2**6)
diff --git a/test/sampler_test.py b/test/sampler_test.py
index 963813f007d32fa9c66214c36ed9a5f19a1207b0..d4cac1ae3e35cceb27284dd29b42738ea521d9ec 100644
--- a/test/sampler_test.py
+++ b/test/sampler_test.py
@@ -5,7 +5,6 @@ from bilby.core.result import Result
import unittest
from mock import MagicMock
import numpy as np
-import inspect
import os
import copy
@@ -393,5 +392,63 @@ class TestPymultinest(unittest.TestCase):
self.assertDictEqual(expected, self.sampler.kwargs)
+class TestRunningSamplers(unittest.TestCase):
+
+ def setUp(self):
+ np.random.seed(42)
+ bilby.core.utils.command_line_args.test = False
+ self.x = np.linspace(0, 1, 11)
+ self.model = lambda x, m, c: m * x + c
+ self.injection_parameters = dict(m=0.5, c=0.2)
+ self.sigma = 0.1
+ self.y = self.model(self.x, **self.injection_parameters) +\
+ np.random.normal(0, self.sigma, len(self.x))
+ self.likelihood = bilby.likelihood.GaussianLikelihood(
+ self.x, self.y, self.model, self.sigma)
+
+ self.priors = dict(
+ m=bilby.core.prior.Uniform(0, 5), c=bilby.core.prior.Uniform(-2, 2))
+
+ def tearDown(self):
+ del self.likelihood
+ del self.priors
+ bilby.core.utils.command_line_args.test = False
+
+ def test_run_cpnest(self):
+ _ = bilby.run_sampler(
+ likelihood=self.likelihood, priors=self.priors, sampler='cpnest',
+ nlive=100, save=False)
+
+ def test_run_dynesty(self):
+ _ = bilby.run_sampler(
+ likelihood=self.likelihood, priors=self.priors, sampler='dynesty',
+ nlive=100, save=False)
+
+ def test_run_emcee(self):
+ _ = bilby.run_sampler(
+ likelihood=self.likelihood, priors=self.priors, sampler='emcee',
+ nsteps=1000, nwalkers=10, save=False)
+
+ def test_run_nestle(self):
+ _ = bilby.run_sampler(
+ likelihood=self.likelihood, priors=self.priors, sampler='nestle',
+ nlive=100, save=False)
+
+ def test_run_ptemcee(self):
+ _ = bilby.run_sampler(
+ likelihood=self.likelihood, priors=self.priors, sampler='ptemcee',
+ nsteps=1000, nwalkers=10, ntemps=10, save=False)
+
+ def test_run_pymc3(self):
+ _ = bilby.run_sampler(
+ likelihood=self.likelihood, priors=self.priors, sampler='pymc3',
+ draws=50, tune=50, n_init=1000, save=False)
+
+ def test_run_pymultinest(self):
+ _ = bilby.run_sampler(
+ likelihood=self.likelihood, priors=self.priors,
+ sampler='pymultinest', nlive=100, save=False)
+
+
if __name__ == '__main__':
unittest.main()