diff --git a/bilby/core/sampler/__init__.py b/bilby/core/sampler/__init__.py
index 0ed5291ac4144b2dbe9fbdb53fc3929f32fd4450..0b529ca2acc995ac53f3b78de065079b6a690959 100644
--- a/bilby/core/sampler/__init__.py
+++ b/bilby/core/sampler/__init__.py
@@ -16,11 +16,13 @@ from .ptemcee import Ptemcee
from .ptmcmc import PTMCMCSampler
from .pymc3 import Pymc3
from .pymultinest import Pymultinest
+from .fake_sampler import FakeSampler
IMPLEMENTED_SAMPLERS = {
'cpnest': Cpnest, 'dynesty': Dynesty, 'emcee': Emcee, 'nestle': Nestle,
'ptemcee': Ptemcee,'ptmcmcsampler' : PTMCMCSampler,
- 'pymc3': Pymc3, 'pymultinest': Pymultinest, 'pypolychord': PyPolyChord }
+ 'pymc3': Pymc3, 'pymultinest': Pymultinest, 'pypolychord': PyPolyChord,
+ 'fake_sampler': FakeSampler }
if command_line_args.sampler_help:
sampler = command_line_args.sampler_help
diff --git a/bilby/core/sampler/fake_sampler.py b/bilby/core/sampler/fake_sampler.py
new file mode 100644
index 0000000000000000000000000000000000000000..4d4acf818c37ad2402454915a4541efeb991d809
--- /dev/null
+++ b/bilby/core/sampler/fake_sampler.py
@@ -0,0 +1,70 @@
+from __future__ import absolute_import
+
+import numpy as np
+from .base_sampler import Sampler
+from ..result import read_in_result
+
+
+class FakeSampler(Sampler):
+ """
+ A "fake" sampler that evaluates the likelihood at a list of
+ configurations read from a posterior data file.
+
+ See base class for parameters. Added parameters are described below.
+
+ Parameters
+ ----------
+ sample_file: str
+ A string pointing to the posterior data file to be loaded.
+ """
+ default_kwargs = dict(verbose=True, logl_args=None, logl_kwargs=None,
+ print_progress=True)
+
+ def __init__(self, likelihood, priors, sample_file, outdir='outdir',
+ label='label', use_ratio=False, plot=False,
+ injection_parameters=None, meta_data=None, result_class=None,
+ **kwargs):
+ Sampler.__init__(self, likelihood, priors, outdir=outdir, label=label,
+ use_ratio=False, plot=False, skip_import_verification=True,
+ injection_parameters=None, meta_data=None, result_class=None,
+ **kwargs)
+ self._read_parameter_list_from_file(sample_file)
+ self.result.outdir = outdir
+ self.result.label = label
+
+ def _read_parameter_list_from_file(self, sample_file):
+ """Read a list of sampling parameters from file.
+
+ The sample_file should be in bilby posterior HDF5 format.
+ """
+ self.result = read_in_result(filename=sample_file)
+
+ def run_sampler(self):
+ """Compute the likelihood for the list of parameter space points."""
+ self.sampler = 'fake_sampler'
+
+ # Flushes the output to force a line break
+ if self.kwargs["verbose"]:
+ print("")
+
+ likelihood_ratios = []
+ posterior = self.result.posterior
+
+ for i in np.arange(posterior.shape[0]):
+ sample = posterior.iloc[i]
+
+ self.likelihood.parameters = sample.to_dict()
+ logl = self.likelihood.log_likelihood_ratio()
+ sample.log_likelihood = logl
+ likelihood_ratios.append(logl)
+
+ if self.kwargs["verbose"]:
+ print(self.likelihood.parameters['log_likelihood'], likelihood_ratios[-1],
+ self.likelihood.parameters['log_likelihood'] - likelihood_ratios[-1])
+
+ self.result.log_likelihood_evaluations = np.array(likelihood_ratios)
+
+ self.result.log_evidence = np.nan
+ self.result.log_evidence_err = np.nan
+
+ return self.result
diff --git a/bilby/gw/detector.py b/bilby/gw/detector.py
index 1db3f6336d47e464959a5a5ab71fbddab94734d8..6ceb555c662844de4030c74d3d1a32ff81720b9f 100644
--- a/bilby/gw/detector.py
+++ b/bilby/gw/detector.py
@@ -65,7 +65,7 @@ class InterferometerList(list):
raise ValueError("The {} of all interferometers are not the same".format(attribute))
def set_strain_data_from_power_spectral_densities(self, sampling_frequency, duration, start_time=0):
- """ Set the `Interferometer.strain_data` from the power spectal densities of the detectors
+ """ Set the `Interferometer.strain_data` from the power spectral densities of the detectors
This uses the `interferometer.power_spectral_density` object to set
the `strain_data` to a noise realization. See
@@ -86,6 +86,28 @@ class InterferometerList(list):
duration=duration,
start_time=start_time)
+ def set_strain_data_from_zero_noise(self, sampling_frequency, duration, start_time=0):
+ """ Set the `Interferometer.strain_data` from the power spectral densities of the detectors
+
+ This uses the `interferometer.power_spectral_density` object to set
+ the `strain_data` to zero noise. See
+ `bilby.gw.detector.InterferometerStrainData` for further information.
+
+ Parameters
+ ----------
+ sampling_frequency: float
+ The sampling frequency (in Hz)
+ duration: float
+ The data duration (in s)
+ start_time: float
+ The GPS start-time of the data
+
+ """
+ for interferometer in self:
+ interferometer.set_strain_data_from_zero_noise(sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=start_time)
+
def inject_signal(self, parameters=None, injection_polarizations=None, waveform_generator=None):
""" Inject a signal into noise in each of the three detectors.
diff --git a/examples/injection_examples/fake_sampler_example.py b/examples/injection_examples/fake_sampler_example.py
new file mode 100755
index 0000000000000000000000000000000000000000..63c2713737b17586a9c2d296837beec21fd83b43
--- /dev/null
+++ b/examples/injection_examples/fake_sampler_example.py
@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+"""
+Read ROQ posterior and calculate full likelihood at same parameter space points.
+"""
+from __future__ import division, print_function
+
+import numpy as np
+import deepdish as dd
+import bilby
+import matplotlib.pyplot as plt
+
+
+def make_comparison_histograms(file_full, file_roq):
+ # Returns a dictionary
+ data_full = dd.io.load(file_full)
+ data_roq = dd.io.load(file_roq)
+
+ # These are pandas dataframes
+ pos_full = data_full['posterior']
+ pos_roq = data_roq['posterior']
+
+ plt.figure()
+ plt.hist(pos_full['log_likelihood_evaluations'], 50, label='full', histtype='step')
+ plt.hist(pos_roq['log_likelihood_evaluations'], 50, label='roq', histtype='step')
+ plt.xlabel(r'delta_logl')
+ plt.legend(loc=2)
+ plt.savefig('delta_logl.pdf')
+ plt.close()
+
+ plt.figure()
+ delta_dlogl = pos_full['log_likelihood_evaluations'] - pos_roq['log_likelihood_evaluations']
+ plt.hist(delta_dlogl, 50)
+ plt.xlabel(r'delta_logl_full - delta_logl_roq')
+ plt.savefig('delta_delta_logl.pdf')
+ plt.close()
+
+ plt.figure()
+ delta_dlogl = np.abs(pos_full['log_likelihood_evaluations'] - pos_roq['log_likelihood_evaluations'])
+ bins = np.logspace(np.log10(delta_dlogl.min()), np.log10(delta_dlogl.max()), 25)
+ plt.hist(delta_dlogl, bins=bins)
+ plt.xscale('log')
+ plt.xlabel(r'|delta_logl_full - delta_logl_roq|')
+ plt.savefig('log_abs_delta_delta_logl.pdf')
+ plt.close()
+
+
+def main():
+ outdir = 'outdir_full'
+ label = 'full'
+
+ np.random.seed(170808)
+
+ duration = 4
+ sampling_frequency = 2048
+ noise = 'zero'
+
+ sampler = 'fake_sampler'
+ # This example assumes that the following posterior file exists.
+ # It comes from a run using the full likelihood using the same
+ # injection and sampling parameters, but the ROQ likelihood.
+ # See roq_example.py for such an example.
+ sample_file = 'outdir_dynesty_zero_noise_SNR22/roq_result.h5'
+
+ injection_parameters = dict(
+ chirp_mass=36., mass_ratio=0.9, a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
+ phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=0.659,
+ phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+
+ waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
+ reference_frequency=20.0, minimum_frequency=20.0)
+
+ waveform_generator = bilby.gw.WaveformGenerator(
+ duration=duration, sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+ waveform_arguments=waveform_arguments,
+ parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
+
+ ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+
+ if noise == 'Gaussian':
+ ifos.set_strain_data_from_power_spectral_densities(
+ sampling_frequency=sampling_frequency, duration=duration,
+ start_time=injection_parameters['geocent_time'] - 3)
+ elif noise == 'zero':
+ ifos.set_strain_data_from_zero_noise(
+ 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.BBHPriorDict()
+ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'iota', 'psi', 'ra',
+ 'dec', 'phi_12', 'phi_jl', 'luminosity_distance']:
+ priors[key] = injection_parameters[key]
+ priors.pop('mass_1')
+ priors.pop('mass_2')
+ priors['chirp_mass'] = bilby.core.prior.Uniform(
+ 15, 40, latex_label='$\\mathcal{M}$')
+ priors['mass_ratio'] = bilby.core.prior.Uniform(0.5, 1, latex_label='$q$')
+ priors['geocent_time'] = bilby.core.prior.Uniform(
+ injection_parameters['geocent_time'] - 0.1,
+ injection_parameters['geocent_time'] + 0.1, latex_label='$t_c$', unit='s')
+
+ likelihood = bilby.gw.GravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=waveform_generator)
+
+ result = bilby.run_sampler(
+ likelihood=likelihood, priors=priors, sampler=sampler, sample_file=sample_file,
+ injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+ # Make a corner plot.
+ result.plot_corner()
+
+ # Compare full and ROQ likelihoods
+ make_comparison_histograms(outdir + '/%s_result.h5' % label, sample_file)
+
+
+if __name__ == '__main__':
+ main()